Michel Visonneau

T-S Tollmien-Schlichting

LCTM Local Correlation Transition Models RANS Reynolds Averaged Navier-Stokes

This research work aims at assessing transition calculations feasibility in the context of the industrial codes. We discuss the strengths and weaknesses of the Local Correlation Based Transition Models (LCTM) γ and γ-Re θ for two-and three-dimensional configurations, focusing on both modeling and numerical aspects. The purpose is to analyze the potential of these RANSbased models as predictive tools, that can handle automatically and autonomously laminar-toturbulence transition. We evaluate the transition models performance on an extensive number of test cases, covering a wide range of transition mechanisms. In the case of 3D configurations, a crucial point is the modeling of transition due to stationary crossflow modes, which are the main three-dimensional transition mechanism in a low free-stream turbulence environment. For this purpose, we present in this dissertation an original re-calibration of the T c1 crossflow criterion proposed by Menter & Smirnov in 2014 and based on a local formulation of the renowned empirical C1 criterion of Daniel Arnal. The re-calibrated criterion notably extends the existing γ model for stationary crossflow transition prediction around three-dimensional non-wing-like geometries. v

Résumé

Ce travail de thèse vise à évaluer la faisabilité des calculs de transition dans le cadre des codes industriels. Les points forts ainsi que les faiblesses des modèles de transition basés sur des corrélations locales, γ et γ -Re θ , sont analysés sur des configurations bi et tridimensionnelles, en se concentrant sur la modélisation physique et les aspects numériques. L'objectif est d'analyser le potentiel de ces modèles RANS en tant qu'outils prédictifs, capables de gérer automatiquement et de manière autonome la transition d'un écoulement du régime laminaire au régime turbulent. Nous évaluons les performances des modèles de transition sur un grand nombre de cas de test, couvrant un large éventail de mécanismes de transition. Dans le cas des configurations 3D, un point crucial est la modélisation de la transition due aux modes transversaux et stationnaires, qui sont les principaux mécanismes de transition tridimensionnel dans un environnement à faible turbulence. A cet effet, nous présentons dans cette thèse un recalibrage original du critère transversal T c1, proposé par Menter & Smirnov en 2014 et basé sur une formulation locale du célèbre critère transversal C1 de Daniel Arnal. Ce critère recalibré étend notamment le modèle γ existant pour la prédiction de transition due aux modes transversaux et stationnaires autour des géométries tri-dimensionelles complexes. Prof. Christophe Airiau from Paul Sabatier University in Toulouse is sincerely acknowledged for sharing his Falkner-Skan-Cooke solver with us. Further, I would like to

Introduction

Nowadays computational fluid dynamics (CFD) simulations have a fundamental role in design decisions for marine vehicles and aircraft. In high fidelity simulations, such as Reynolds Averaged Navier-Stokes (RANS) approaches, it is common practice to assume the flow completely turbulent. This assumption is appropriate for full scale ships, submarines and aircraft. This is also the case at model scale, for which the flow is often tripped in the proximity of the leading edge, in order to force turbulence within the boundary layer. Nevertheless, there is a wide range of practical applications that deal with low/moderate Reynolds number, where transitional effects are important, if not dominant. For instance, aerial and marine unmanned vehicles, small submarines, but also wind energy systems, as wind turbines, is another growing area for which moderate Reynolds effects are important. Marine unmanned vehicles span a large range of sizes, from miniature devices to small ships, and a multitude of missions, for example ship tracking and oceans surveillance. Many of these will be operating in a Reynolds number regime where laminar-to-turbulence transition plays a crucial role. Its impact might be problematic from practical perspectives. Efficiency, maneuvering, speed, time on station, power requirements are fundamental requirements that must be established to determine if a vehicle can effectively perform the mission it has been designed for. Large differences in drag will be experienced by the vehicles if operating in the laminar, transitional or fully turbulent regime. Appendages forces can be very different depending on the flow state, impacting stability and control. Propulsion losses can be significant for propellers operating at low speeds, experiencing over 30% of efficiency loss. In order to evaluate accurately vehicles resistance and performance, the understanding, accurate prediction, and further controlling of the laminar-to-turbulence transition is required. Transitional effects need to be accounted also for model scale testing, if tripping is not employed, and for methodologies that aim to exploit laminar flows, such as laminar wing design. Not to mention, the improvement of aircraft performance, the reduction of their environmental footprint and lowering of the economical costs, which require an important decrease of the overall drag around the body. Transition modeling is a very hard task. RANS turbulence models are developed for high Reynolds number and fail to predict transition, notably they lead to a very short extended laminar region. Until the beginning of the 21st century, all the hopes were placed in Large Eddy Simulations (LES) and their capabilities to account for transitional effect. Even if practical flows at moderate Reynolds numbers are becoming approachable by LES codes, this method still requires high computational resources, Spalart, [START_REF] Spalart | Strategies for turbulence modelling and simulations[END_REF]. These means are not always available in an industrial milieu. In addition, a proper specification of the turbulence free-stream environment, external disturbances and structures, can be very challenging. A further difficulty LES-related is some sensitivity to the subgrid-scale model, as discussed in Lardeau et al., [START_REF] Lardeau | Large eddy simulation of transitional separated flow over a flat plate and a compressor blade[END_REF]. In order to avoid this dependency, Wall-Resolved LES (WRLES) should be preferred over Wall-Modeled LES (WMLES). WRLES, however, are not a practical approach in transportation flows, because of the stringent grid requirements, and will not be a real option until later than 2045, as mentioned by Spalart, [START_REF] Spalart | Hybrid RANS-LES methods[END_REF]. Similar issues are also shared by Direct Numerical Simulations (DNS) approach, despite their optimal performance for transitional flows as discussed in Durbin et al., [START_REF] Durbin | DNS of bypass transition[END_REF]. Because of DNS and LES computational limitations, from the 2000s, the CFD community has started devoting research efforts to the construction and formulation of transition models in a RANS environment. The aim is to build RANS-based CFD models able to handle automatically and autonomously transition: a relatively simple, but accurate, efficient and robust tools for its prediction. The main implication with RANS is that the averaging eliminates the linear amplification of the disturbances and the discarded linear effects might seem not compatible with transition physics. Nevertheless, one has to consider that, in many applications, transition occurs in small flow area and it is dictated by geometry features, pressure gradients, and flow separations. A RANS model is able to capture these effects with sufficient engineering accuracy, upon the inclusion of proper correlations in the models formulation. Local correlation-based transition models (LCTM) have been proposed based on this idea. The equations do not intend to model the real physics, the aim is to identify the region of the flow that is laminar and the transition location. The physics is entirely contained in empirical correlations that account for a specific transition process. This formulation gives some degrees of freedom, because different mechanisms can be included in the model, as long there is an empirical correlation that describes them, without demanding a considerable equations reformulation. Indeed, transition is not unique and it can occur in many several ways. This is one of the main difficulties related to its modeling. The approach proposed within the LCTM idea is very distant from transition models based on linear stability theory. Among them, the most successful is the e N method of Smith & Gamberoni, [5], and Van Ingen, [START_REF] Van Ingen | A suggested semi-empirical method for the calculation of the boundary layer transition region[END_REF]. This method has acquired a huge popularity over the years, and it is so far one of the transition models which performs better. Nevertheless, it is not free of empiricism. The critical N factor is not universal and it depends on the free-stream turbulence environment, wind tunnel and specimen characteristics used in the experiments. Moreover, it requires several non-local operations inconvenient for High Performance Computing (HPC) infrastructures, that massively make use of parallel domain decomposition. If we evaluate the various methods with respect to their suitability for industrial 3D design calculations in unstructured and parallelized CFD codes, the only models that traditionally meet these requirements are the low-Re models from Jones & Launders, [START_REF] Jones | The calculation of low-Reynolds-number phenomena with a two-equation model of turbulence[END_REF]. At first, these models have been tested by Rodi & Scheuerer, [START_REF] Rodi | Calculation of heat transfer to convection-cooled gas turbine blades[END_REF], for transition predictions, proving to be inaccurate. At best, they can predict bypass transition, which is dominated by diffusion effects from the free-stream into the boundary layer. A proper recalibration is needed for these models in order to account for transition, and this is not feasible without a substantial re-formulation of the equations. The alternative to low-Re models are the empirical correlations, that generally correlate the transition momentum thickness to turbulence intensity and pressure gradient. Notable examples are the one from Abu-Ghannam & Shaw, [START_REF] Abu-Ghannam | Natural transition of boundary layers-the effects of turbulence, pressure gradient, and flow history[END_REF], and Mayle, [START_REF] Mayle | The role of laminar-turbulent transition in gas turbine engines[END_REF]. These correlations, however, require the computation of integral boundary layer quantities, precluding them from being CFD-friendly. The wish for models formulated locally, built as a set of equations as generic as possible, has led to the LCTM philosophy. The widespread transition models γ -Re θ , by Langtry, [START_REF] Langtry | A correlation-based transition model using local variables for unstructured parallelized CFD codes[END_REF] and γ, by Menter et al., [START_REF] Menter | A one-equation local correlation-based transition model[END_REF], are based on it. This dissertation is devoted to the analysis and discussion of these two models performance on different canonical configurations. This research was conducted within the NATO collaborative group AVT-313, which has gathered researchers from NATO countries, both from experimental and numerical sides. The scope of this activity is to assess the existing capabilities, in regards to CFD codes, to predict transitional flows of interest, as well as an assessment of relevant experimental data for the validation of transition predictions. The interest of NATO lies in the further development of unmanned vehicles that are expected to dramatically change the military battle space of the future. The constitution of this collaborative group reflects the present deep interest towards laminar-to-turbulence transition modeling in a RANS framework. As matter of fact, two AIAA Workshops dedicated to transition have been organized in the last two years and another Workshop transition-devoted will be held within the 34th SNH Symposium in July 2022. Within the AVT-313, two other PhD dissertations were defended in the last six months, by Dongyoung Kim at University of Iowa1 and by Rui Lopes at IST Lisbon, [START_REF] Lopes | Simulation of Transition from Laminar to Turbulent Regime in Practical Applications of Incompressible Flow[END_REF]. Despite the research efforts dedicated to this topic, the current understanding is far from complete. RANS transition models are still in a validation stage. Already existing CFD tools need to be improved for capturing the main transitional phenomena and this can be achieved only if capable and accurate physical models are available in the codes. Robust transition models have to be established together with accurate turbulence models. As originally published in 2014, γ and γ -Re θ transition models accounted only for bypass and Tollmien-Schlichting dominated transition. Some stationary crossflow criteria have been proposed in the last years, mainly thought and built to be implemented within γ -Re θ model formulation. In this dissertation, we propose a re-calibration of the T c1 criterion proposed by Menter & Smirnov in 2014, [15], that extends the existing γ model to stationary-crossflow dominated transition. The T c1 is based on a local formulation of the C1 crossflow criterion by Arnal, [START_REF] Arnal | Théorie de l'instabilité laminaire et critères de transition en écoulement bi et tridimensionnel[END_REF], reconstructed using the solutions of the Falkner-Skan-Cooke boundary layer equations. It uses the directional change in the wall normal direction of the normalized vorticity vector as indicator of crossflow strength. As initially proposed by Menter & Smirnov, the criterion does not work accurately on the 6:1 prolate spheroid, which is the main configuration we have focused our attention on. For this reason, we propose a new recalibrated version of the T c1 crossflow criterion that now accounts for the sweep angle φ effect and enhances the criterion performance on non-wing-like geometries.

Manuscript organization

• Chapter 2 is devoted to a general introduction to transition physics. We will focus on the main transition mechanisms, that will be further discussed in this dissertation. We present a literature review of the existing transition modeling approaches. From the e N method, we will cover some of the different strategies, highlighting their strengths and limitations. We will conclude the chapter with the introduction of the LCTM concept: a set of transport equations built using local information, based on experimental correlations and to be coupled to already existing turbulence models.

• Chapter 3 is dedicated to the mathematical formulation of the RANS models under study.

The full description of γ -Re θ and γ models is presented with the intention of explaining the role of the different terms that enter the formulation. The description of k -ω SST turbulence is provided as well. The k -ω SST 2003 version from Menter et al., [START_REF] Menter | Ten years of industrial experience with the sst turbulence model[END_REF], is chosen for the turbulence coupling.

• Chapter 4 presents a description of ISIS-CFD incompressible solver, where both transition models have been implemented. Algorithms and methodologies employed for the spatial equations discretization are presented in order to introduce the numerical set-up of the simulations discussed in this work.

• Chapter 5 presents the assessment of the two transition models for two-dimensional configurations. The flow over the ERCOFTAC flat plate, for both natural and bypass transition, and Eppler387 and NACA0015 airfoils are considered. The latter two profiles were tested in order to assess the capability of the transition models to predict laminar separation bubbles. The flow around both airfoils presents an extended laminar region, and the laminar separation bubble strongly affects the pressure and skin friction distributions around the body. These three cases were chosen because of the availability of experimental data to compare with. All the results have been submitted as contributions to the 2D "AVT-313 Incompressible Laminar-to-Turbulent Flow Transition Study Comparison Workshop", which is an ongoing exercise within the NATO/AVT-313 research group. Numerical results computed by γ -Re θ and γ, as implemented in different solvers, have been compared in order to assess their consistency and robustness. An overall summary of these results will be presented at the 34th SNH Symposium in July 2022.

• Chapter 6 is devoted to the discussion of the T c1 crossflow criterion recalibration. We first introduce the original criterion as presented by Menter & Smirnov to explain the rationale behind our choice of pursuing a new calibration. The new criterion is presented along with its strengths and limitations.

• Chapter 7 is dedicated to the discussion of the numerical results around three-dimensional configurations: the 6:1 prolate spheroid at three different angles of incidence α = 5 , 15 • were analyzed. As a support to the discussion, we will also present numerical results computed by γ -Re θ using the widely spread helicity-based crossflow criterion, by Grabe et al., [START_REF] Grabe | Transport modeling for the prediction of crossflow transition[END_REF]. We will also consider the possible transition mechanism(s) that γ and γ -Re θ do not account for.

• In Chapter 8, we summarize the conclusions of this work. We will recapitulate the observed most crucial aspects, discussing the limitations of γ and γ -Re θ models as predictive tools. We outline possible paths to enhance their performance. Some suggestions are given on possible future works. Special emphasis is put on the ideal experimental data-set that we ought to have to further validate the models. Indeed, one of the purposes of this research work is to give experimentalists ideas about flow measurements that we are keen to have for further improving RANS correlation-based transition models.

Chapter 2 Background

General Introduction to Linear Stability

The purpose of this section is to give a brief and general introduction to hydrodynamic stability theory, which analyses how laminar flows respond to small disturbances. For a complete treatment of this subject the reader should refer to specialized texts as Schlichting, [START_REF] Schlichting | Boundary layer theory[END_REF], Chandrasekhar, [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF], Drazin & Reid, [START_REF] Drazin | Hydrodynamic stability[END_REF], Charru, [START_REF] Charru | Hydrodynamic instabilities[END_REF], or Schmid & Henningson, [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF]. The organization of this section follows the notes of the "Hydrodynamic Stability" course held by Dr. Carlo Cossu at Centrale Nantes in 2021, [START_REF] Cossu | Introduction to hydrodynamic instability[END_REF]. Transition studies are based on linear stability theory, the mathematical analysis of disturbances of small amplitude superimposed to a base laminar flow. Linear stability analysis allows to compute the unstable frequencies, their amplification, and their amplitude, until non-linear effects become important. In order to assess the stability of a flow, it is appropriate to think in terms of the kinetic energy E V of the perturbations in a volume V:

E V = 1 2 V - → u • - → u dV. (2.1)
and understand under which conditions it is amplified. Let us derive an evolution equation for E V from the Navier-Stokes equations. For an incompressible and viscous flow, the dimensionless Navier-Stokes equations read:

∇ • - → u = 0, (2.2)
∂ - → u ∂t + - → u • ∇ - → u = - 1 ρ ∇p + 1 Re ∆ - → u . (2.3) We introduce the definition - → u = - → u + - → U ,
where -→ U is the basic flow and -→ u is the perturbations vector. Removing the equations satisfied by the basic flow, the equations for the perturbations -→ u are found:

∇ • - → u = 0, (2.4) ∂ - → u ∂t + - → u • ∇ - → U + - → U • ∇ - → u = - 1 ρ ∇p + 1 Re ∆ - → u . (2.5)
The evolution equation of the disturbances kinetic energy E V is obtained multiplying Eq. (2.5) by -→ u and dropping out the non-linear terms. Rearranging the equation in order to obtain each term in a gradient form, and applying the divergence-free condition of the flow, the equation reads:

1 2 ∂ - → u • - → u ∂t = -( - → u ⊗ - → u ) : ∇ - → U P -∇• u • u 2 u + - → U + p - → u T - 1 Re - → u •∇ - → u - 1 Re ∇ - → u : ∇ - → u D , (2.6) where - → u ⊗ - → u = u i u j and ∇ - → u : ∇ - → u = ∂u i ∂x j ∂u i
∂x j , where the Einstein's convention is used for repeated indices. Integrating over a generic volume V Eq.(2.6), we find the Reynolds-Orr equation:

dE V dt = - V ( - → u ⊗ - → u ) : ∇ - → U dV - 1 Re V ∇ - → u : ∇ - → u dV. (2.7)
The integral of the term T in Eq.(2.6) is transformed into the boundary flux term using Gauss' theorem. It disappears in Eq.(2.7) when no-slip conditions are enforced at the wall and periodic boundary conditions are imposed as inflow and outflow. The terms resulting from the integration are: the viscous dissipation term D' -V DdV and the production term

P' = -V ( - → u ⊗ - → u ) : ∇ - → U dV . D'
either works as a sink for viscous flows, transforming the perturbation energy into heat, either it is zero for non-viscous flows. P' plays instead a fundamental role. A positive production is a necessary condition for E V to increase: P' > 0 corresponds to the situation where the tensor --→ u ⊗ -→ u extracts energy from the non-zero basic flow shear and transfer it to the perturbations. In particular, since a positive production term requires a basic flow shear (∇U = 0), it can be concluded that this is the core mechanism responsible of the instabilities. Averaging --→ u ⊗ -→ u , we obtain the Reynolds stress tensor.

Orr

∂u ∂x + ∂v ∂y + ∂w ∂z = 0, (2.8
)

∂u ∂t + U ∂u ∂x + v dU dy = - ∂p ∂x + 1 Re ∆u , (2.9
)

∂v ∂t + U ∂v ∂x = - ∂p ∂y + 1 Re ∆v , (2.10
)

∂w ∂t + U ∂w ∂x = - ∂p ∂z + 1 Re ∆w . (2.11)
Taking the divergence of Eq.(2.9)-(2.11) and using Eq.(2.8), an equation for the perturbation pressure can be built and used to obtain an equation for the normal velocity1 , v :

∂ ∂t + U ∂ ∂x ∆ - d 2 U dy 2 ∂ ∂x - 1 Re ∆ 2 v = 0.
(2.12)

In order to completely describe a three dimensional configuration, a second equation for the normal vorticity ,

η η = ∂u ∂z - ∂w ∂x , (2.13)
is needed. The equation reads as:

∂ ∂t + U ∂ ∂x - 1 Re ∆ η = - d 2 U dy 2 ∂v ∂z . (2.14)
We define D the differentiation with respect to the inhomogeneous coordinate direction y.

Eq.(2.12) and Eq.(2.14) paired with the boundary conditions v = Dv = η = 0 at the wall and in the far field, and with the initial condition v ( -→ x , 0) = v 0 and η ( -→ x , 0) = η 0 provide the description of the evolution of an arbitrary disturbance in space and time. Let us introduce in Eq.(2.12) and Eq.(2.14) wavelike solutions of the form: ṽ(y)e i(αx+βz-ωt) , η(y)e i(αx+βz-ωt) , (

where ṽ and η are the amplitude functions that depends on y. α, β ∈ R are the wave numbers along x and z directions, respectively, and ω ∈ C is the wave frequency. Non-real frequency and real wave numbers are the setting of the temporal problem: the spatial structure of the wavelike perturbation rests unchanged, while the amplitude of the wave grows or decays in time.

The final resulting equation is:

(-iω + iαU )(D 2 -k 2 ) -iα d 2 U dy 2 - 1 Re (D 2 -k 2 ) 2 ṽ = 0, (2.16)
(-iω + iαU ) - 1 Re (D 2 -k 2 ) η = -iβ dU dy ṽ, (2.17)
with boundary conditions ṽ = Dṽ = η = 0 at the wall and in the free-stream. Eq.(2.16) is the classical Orr-Sommerfeld equation, while Eq.(2.17) is known as the Squire's equation. The frequency ω is the eigenvalue of the Orr-Sommerfeld equation associated to the eigenvector ṽ.

Squire's Theorem. Given Eq.(2.16), a relation between two-dimensional and three-dimensional solutions can be stated. Rewriting the frequency as:

ω = αc, (2.18)
where c is the phase speed, Eq.(2.16) can be rewritten as:

(U -c)(D 2 -k 2 ) - d 2 U dy 2 - 1 iαRe (D 2 -k 2 ) 2 ṽ = 0, (2.19)
where k is the norm of the wave vector {α, 0, β}. We compare Eq.(2.19) to the Orr-Sommerfeld equation in for 2D wavelike solutions, i.e. β = 0, that can be written as:

(U -c)(D 2 -α 2 sq ) - d 2 U dy 2 - 1 iα sq Re sq (D 2 -α 2 sq ) 2 ṽ = 0. (2.20)
We observe that Eq.(2.19) and Eq.(2.20) have the same solutions if:

α 2 sq = k = (α 2 + β 2 ), (2.21
)

Re sq = Re α α sq . (2.22)
The relations in Eq.(2.21) and (2.22) are said Squire's transform. They relate a modal 3D solution at a given Reynolds number Re to a 2D modal solution at the lower Reynolds number Re sq = Re α αsq < Re. Following Squire's transformation, Squire's theorem asserts that parallel shear flows become first unstable to 2D wavelike perturbations at a Reynolds number smaller than any Reynolds number for which unstable three dimensional perturbations exist. Squire's transformation allows us to limit to the case of two-dimensional perturbations. In this case w = 0 and ∂ ∂z → 0, the wall-normal vorticity is zero and Squire's equation, Eq.(2.17), is identically satisfied. For imposed values of α and Re, given the base flow U (y), c = c r + ic i is the complex eigenvalue of linear system Lṽ = cṽ. In the (Re, α) diagram, the neutral curve is c i = 0, for c i > 0 the disturbances are amplified, growing exponentially. In the stable region c i < 0 disturbances are damped. An example of (Re, α) diagram for a Blasius boundary layer flow is shown in Fig. (2.1). The neutral curve c i = 0, on which the critical Reynolds number lies (red dot) divides the unstable (shaded) region from the stable region. The unstable region is bounded by an upper and a lower branch that define the unstable wave-band α for each Re.

Let us now distinguish between inviscid and viscous case. If the flow is inviscid, Eq.(2.19) reduces to:

(U -c)(D 2 -α 2 )ṽ - d 2 U dy ṽ = 0, (2.23)
which is known as the Rayleigh equation. In this case, Rayleigh proved that a necessary condition for instability, i.e. a perturbation with c i > 0, is that the velocity U (y) has at least an inflection point. This inflection point is a maximum of the spanwise mean vorticity (Fjotorft's criterion) 2 .

Examples of inflectional instabilities are free shear layers, jets and wakes.

In the viscous case, Rayleigh's criterion is not valid anymore. Nevertheless, solving the Orr-Sommerfeld equation for the plane Poiseuille flow, a slightly unstable eigenmode can be found even if the parabolic profile does not present any inflection point. This unstable eigenmode is called a Tollmien-Schlichting (T-S) wave. It takes the name from the first two researchers who showed its existence, Eckert, [START_REF] Eckert | A kind of boundary-layer 'flutter': The turbulent history of a fluid mechanical instability[END_REF]. The existence of an unstable mode comes from the destabilizing effect of the viscosity at large Reynolds number. Going back to the Reynolds-Orr equation, the viscous effects are felt through the dissipation term (-1 Re V ∇u : ∇u dV ), which reduces the growth rate of the instability. At high Re, viscosity induces positive production terms, because of the no-slip conditions at the wall, which produce non zero Reynolds stresses. Compared to inflectional instabilities, the typical growth of viscous instabilities is smaller, because they scale over viscous time scales which are larger at high Reynolds number, Tviscous Tconvective = Re, and they disappear for Re → ∞. For three-dimensional perturbations, the situations is more complex, because the solution η of Eq.(2.17) is not identically zero. Instead, it depends on the solution ṽ of Eq. (2.16). ṽ acts on the evolution of η through the non-symmetric coupling term iβ dU dy ṽ in Eq.(2.17). The existence of this term makes the overall Orr-Sommerfeld-Squire system non-normal, i.e. the operator does not commute with its adjoint LL + = L + L. The non-normality of the operator is a necessary condition for the existence of transient energy growth, that is the amplification of a perturbation at finite times. In the Orr-Sommerfeld-Squire equations, the strong non-normal coupling causes the energy amplification of wall-normal vorticity, visualized as streamwise streaks, induced by spanwise periodic distributions of wall-normal velocity, the streamwise vortices. The latter displace low-speed fluid away from the wall and push high-speed fluid towards it. This process, the lift-up mechanism, creates spanwise alternated regions of high-and low-speed fluid, the streaks. A smoke visualization of the streamwise streaks is presented in Fig. (2.2). The appearance of these streaks strongly affects the basic flow state at finite times, making it unstable to secondary instabilities that can cause breakdown of the laminar flow at Re lower than the Re c obtained by linear stability theory.

Transition Mechanisms

The previous section has served to introduce the Orr-Sommerfeld-Squire equations. They are needed to identify the type of instabilities at different flow conditions and their possible growth. In order to asses if these instabilities might trigger or not laminar-to-turbulent transition, it is more appropriate to identify the different transition mechanisms. Laminar-toturbulence transition process can be subdivided in different stages:

• the receptivity stage describes the means by which the ambient disturbances (free-stream turbulence intensity, T u = 2k 3U 2 , and noise, surface roughness, vibrations...) enter the laminar boundary layer as steady or unsteady fluctuations of the basic state. This process provides the initial conditions of amplitude and frequencies for the breakdown of the laminar flow.

• the linear stability stage describes the initial growth of these small disturbances until they reach sizes where non-linear effects become important. The amplification, whose form depends on the ambient disturbances environment (exponential growth, transient growth...), is tuned by local gradients, such as pressure and temperature. When the waves amplitude reaches a finite value, energy is redistributed among disturbances and the flow is transformed into a new base, possibly quasi-steady state (non-linear saturation stage). This new base flow is unstable to the so-called secondary instabilities (secondary instabilities stage), whose amplification is responsible of the breakdown to turbulence. Multiple types of disturbances can co-exist and possibly interact.

• 3D non-linear interactions is the last stage. Non-linearities and higher instabilities excite a wider number of scales and frequencies in the flow and rapidly lead to the breakdown to turbulence. The breakdown stage is characterized by the emergence of turbulent spots, that grow in size and eventually merge into a fully turbulent boundary layer. Their appearance is accompanied by a strong increase in skin friction.

Let us now address the different type of transition mechanisms that can occur depending on the external disturbance environments. A paragraph is dedicated to the different transition mechanisms, highlighting those mechanisms that will be treated more in detail in the rest of the thesis.

Natural Transition

In literature, it is referred to natural transition when transition emanates from exponential instabilities. It is commonly accepted that this kind of transition occurs for low free-stream turbulence level T u of less than 1%, Mayle, [START_REF] Mayle | The role of laminar-turbulent transition in gas turbine engines[END_REF]. Different types of instabilities can lead to natural transition, such as attachment line instabilities at the leading edge, T-S instabilities, crossflow and centrifugal (Görtler) instabilities on concave wall. The mentioned mechanisms can be identified on a swept wing in different region, see The waves normally occur in zone of zero or mild adverse pressure gradient, because favorable pressure gradients have a stabilizing effect. T-S waves can occur in both 2D and 3D boundary layers, in the latter case they are often referred to as streamwise instabilities, because identified by their propagation direction. As discussed in Loiseau et al., [START_REF] Loiseau | Investigation of the roughnessinduced transition: global stability analyses and direct numerical simulations[END_REF], the presence of localized or distributed surface roughness can promote or delay the transition process. Below a given threshold, based on Reynolds number, shape and spacing of the roughness elements, as well as their height with respect to the boundary layer thickness, 3-D roughness elements induce streamwise velocity streaks that stabilize T-S waves and delay natural transition process. Above this threshold, on the contrary, these large-amplitude velocity streaks can become unstable rapidly leading to the breakdown to turbulence downstream the roughness elements.

Görtler vortices can also occur in two-and three-dimensional boundary layers on concave surfaces. They manifest themselves as streamwise counter-rotating vortices inside the boundary layer. They are named after Henry Görtler, who first analytically predicted their occurrence.

If the boundary layer thickness δ is comparable to the radius of curvature R, concave surface boundary layer flows are subjected to centrifugal instabilities. These are caused by the imbalance between radial pressure gradient and centrifugal force. Instability of velocity profiles on a curved wall comes from Rayleigh's instability criterion for a basic swirling flow, with an arbitrary dependence on the angular velocity ω(r) on the distance r from the center of rotation. In this frame, the Rayleigh inflection criterion can be restated as if 1 r 3 dr 4 ω(r) 2 dr < 0, the flow is unstable, where r 2 ω(r) is the angular momentum. In three dimensional boundary layers, if the sweep angle is large compared to 1/ √ Re, Görtler mechanism becomes unimportant, as discussed in Hall, [30]. It is worthwhile to mention that Görtler vortices can also appear in the vicinity of the stagnation point on a convex surface, as discussed by Hirschel et al., [START_REF] Hirschel | Three-dimensional attached viscous flow[END_REF] and shown in Fig. (2.5), because of the local concavity of the streamlines at the leading edge. This behavior could eventually produce some streamwise vorticity that can destabilize the laminar boundary layer further downstream. A more detailed review of Görtler vortices can be find in Saric, [START_REF] Saric | Görtler vortices[END_REF]. Leading edge contamination phenomenon was observed for the first time in 1950s during some experiments on swept laminar-flow wings. It is common around swept wings with large leading edge radius and for large sweep angles. Indeed, around an unswept wing, the boundary layer starts at the leading edge and develops on the upper and lower surface. The favorable pressure gradients stabilize the flow within the boundary layer and it remains laminar until the change of sign of the pressure gradient. On the contrary, around a swept wing, the geometrical sweep of the wing induces a strong velocity component in the direction of the leading edge, causing the boundary layer to be turbulent right at the attachment line. The boundary layer is strongly influenced by the diverging streamlines at the leading edge, see Fig. (2.6). The turbulence, traveling in the spanwise direction, can also stem from the turbulent boundary layer around the fuselage or because of an high free-stream turbulence environment, hence the name "contamination". Once these disturbances reach a certain amplitude traveling in the spanwise direction along the attachment line, transition will occur, also affecting transition in the chordwise direction of the wing, unless relaminarization due to strong favorable pressure gradients. The local flow at the attachment line is commonly approximated as an Hiemenz flow, with velocity profiles similar to Blasius profile.

The key parameter to study the stability to small disturbances of the laminar flow is the sweep Reynolds number:

R = W e (νS) (2.24)
where W e is the spanwise component of the velocity at the boundary layer edge, S is the mean strain rate of the irrotational flow and ν the kinematic viscosity. The linear stability limit calculated by Hall, [START_REF] Hall | On the instability of a three-dimensional attachment-line boundary layer: weakly nonlinear theory and a numerical approach[END_REF], and confirmed by DNS simulation performed by Spalart, [START_REF] Spalart | Direct numerical study of leading-edge contamination[END_REF], is R=583, i.e. the flow is laminar and stable if R <583. This limit value is considerably lower, in case of boundary layer contamination, coming from the fuselage, which the wing is attached to, Juillen & Arnal , [START_REF] Juillen | Etude de la transition et de la contamination de bord d'attaque sur ailes en fleche[END_REF]. Two different reference systems can be defined: (x, y, z) system is related to the wing and (x 1 , y, z 1 ) to the external streamline. y represents the direction normal to the wall. U (y) and W (y) are the projections of the local velocity Q(y) along x, z, while U 1 (y) and W 1 (y) are the projections along x 1 , z 1 . For y → ∞, U → U e and W → W e , while W 1 → 0 and The continuity and momentum equations for the averaged flow written in the (x, y, z) reference system are:

U 1 → (U 2 e + W 2 e ).
∂U ∂x + ∂V ∂y = 0, ( 2
U ∂W ∂x + V ∂W ∂y = ν ∂ 2 W ∂y 2 , (2.27)
where it is considered that ∂ ∂z = 0 , coming from the hypothesis of infinite swept wing; The two profiles U (y) and W (y) are obtained through similarity solutions. The streamwise U 1 and the crosswise W 1 velocity profiles can be obtained by projection, as:

U 1 = U cos φ + W sin φ, W 1 = -U sin φ + W cos φ, (2.28)
where φ = arctan W e U e .

(2.29)

W 1 has its maximum close to the middle of the boundary layer and it goes towards zero at the wall, enforced by the no-slip conditions, and at the boundary layer edge. Thus, the profile has an inflection point and, according to Rayleigh criterion, this velocity profile is dynamically unstable. The external inviscid streamline will have an inflection point at the abscissa where the distribution of U e (x) has a maximum. At this location the streamwise pressure gradient changes sign and the transversal velocity profile W 1 (y) reverses at the wall assuming an S form.

If the pressure gradient is strong enough, W 1 (y) reverses completely. , [START_REF] Yiming | Boundary-layer transition of advanced fighter wings at high-speed cruise conditions[END_REF].

The eigenfunctions associated to crossflow instabilities are counterotating vortices, but when they reach a finite amplitude, they produce a system of corotating streamwise vortices aligned with the local velocity vector. This system of corotating vortices superimposed to the mean flow represents the new base flow unstable to secondary instabilities, that will lead to transition. These vortices are known as crossflow vortices and they are experimentally visualized as streaks on the surface. According to linear stability theory, the range of unstable frequencies is wider compared to streamwise instabilities. Indeed, crossflow instabilities can amplify zero frequency disturbances, the so-called stationary crossflow instabilities, opposed to non-zero frequency traveling crossflow. For traveling instabilities, it is meant any non-stationary instabilities of a convective nature (i.e., those that propagate downstream and amplify in amplitude in the propagation process), such as Tollmien-Schlichting and non-stationary crossflow instability. The dominant role of one of these two kinds of crossflow instabilities, either stationary or traveling, depends on the receptivity mechanism. It is commonly accepted that traveling crossflow instabilities trigger transition in high disturbance environments and are excited by free-stream turbulence. Stationary crossflow appears in low disturbance environments and are excited by small roughness elements Arnal & Délery, [START_REF] Délery | Laminar-turbulent transition and shock wave/boundary layer interaction[END_REF]. In literature a turbulence intensity of T u ∼ 0.15% is defined as the limiting value to distinguish between traveling and stationary crossflow, Deyhle & Bippes, [START_REF] Deyhle | Disturbance growth in an unstable three-dimensional boundary layer and its dependence on environmental conditions[END_REF]. Nevertheless, experimentally it was observed that in the range 0.1% < T u < 0.3% the type of CF transition cannot be assessed with exactitude. Indeed, the experiments, performed by White et al., [START_REF] White | Stages of swept-wing transition[END_REF], conducted on a swept wing model with variable-amplitude roughness elements and a turbulence intensity T u = 0.3%, have shown that turbulence intensity itself cannot be considered as the discriminating element of the crossflow type. Surface roughness, turbulence intensity and their interaction must be considered together.

In terms of flow visualization on the surface, stationary wave dominated transition is characterized by a sharp saw-tooth transition pattern, while a diffuse span-wise invariant transition front is indicative of traveling waves, Saric et al., [START_REF] Saric | Stability and transition of three-dimensional boundary layers[END_REF].

Separation induced transition

Another important transition mechanism is separation induced transition. The boundary layer detaches under the influence of an adverse pressure gradient and the separated shear layer undergoes transition because of the disturbances amplification in the unstable laminar layer. The momentum transfer in the normal direction, caused by the turbulent mixing, eventually eliminates the reversed velocities near the wall causing the boundary layer to reattach. This re-circulation zone is called Laminar Separation Bubble (LSB), whose schematic representation is shown in Fig .(2.9) on the left. The zone bounded by the surface and the dividing streamline ST R represents the re-circulatory flow forming the bubble. The zone between the divided streamline and the outer edge of the boundary layer is the separated shear layer, where the flow undergoes transition at the location T . R is the reattachment location. The effect of the laminar separation bubble on the overall pressure coefficient distribution is shown in Fig .(2.9) on the right. As the flow separates at location S, the edge of the bubble ahead of the transition position, T , is a zero pressure gradient streamline, as denoted by the pressure plateau. The primary instability in a separation bubble is inflectional. It originates upstream of the separation point, because of the strong adverse-pressure gradient, and it is convected downstream, as mentioned in Diwan & Ramesh, [START_REF] Diwan | On the origin of the inflectional instability of a laminar separation bubble[END_REF]. The laminar separation bubble length depends on Reynolds number and free-stream turbulence intensity. As T u and Re increase, the bubble length is reduced and the suction peak increases in magnitude, as described in Jahanmiri, [START_REF] Jahanmiri | Laminar separation bubble: its structure, dynamics and control[END_REF]. On the other hand, for decreasing Re the viscous damping effects become significant, suppressing transition and, eventually, delaying reattachment, Saxena, [START_REF] Saxena | The laminar separation bubble[END_REF]. The angle of attack α plays a major role on the separation point, increasing α the separation point moves upwards, because of the stronger pressure gradient at higher incidences. The distinction between short and long bubbles depends on their effect on the overall pressure distribution. A short separation bubble affects the pressure distribution locally, because it reattaches shortly after the separation point. Nevertheless, a short separation bubble might "burst" into a long one or into an unattached shear layer, causing a considerable loss of lift. Figure 2.9: Schematic representation of a laminar separation bubble: streamlines and velocity profile (left) from O'Meara & Muller, [START_REF] O'meara | Laminar separation bubble characteristics on an airfoil at low Reynolds numbers[END_REF] and the experimental pressure coefficient distribution (right) from Lee et al., [START_REF] Lee | Mechanisms of surface pressure distribution within a laminar separation bubble at different Reynolds numbers[END_REF].

Bypass Transition

The instabilities discussed above arise as exponential growing eigenmodes. Their growth, together with the subsequent secondary instabilities, give an understanding of the transitionto-turbulence process. Nevertheless, transition does not always emanate from exponential instabilities. We talk about "bypass" transition, when the appearance of T-S waves, spanwise vorticity (Görtler, crossflow vortices) and the three dimensional breakdown is bypassed. Turbulent spots directly appear within the boundary layer. In Schmid, [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF], it is proposed as definition for bypass transition transition emanating from nonmodal growth mechanisms. This is the case when the linear stability operator, built from the linearization of the Navier-Stokes equation, is non-normal, and a spectral analysis cannot describe the full dynamics of the equations. Transition mechanism is then characterized by the appearance within the boundary layer of elongated structures in the streamwise directions, the streamwise streaks described in Sec.(2.1) and shown in Fig. (2.2). This is the case either in high free-stream turbulence environment, T u > 1% is set as limit, either in the case of initial disturbances originating at the wall because of surface roughness. This kind of transition is typical around turbines, compressors or the rudder of a ship. In Mayle et al., [START_REF] Mayle | The path to predicting bypass transition[END_REF], it is shown that the first turbulent spot in the boundary layer appears when the skin friction profile deviates from the laminar profile, thus every other stage, through which the instabilities go before the appearance of turbulent spots, occur in a laminar boundary layer. Laminar fluctuations are the precursor of bypass transition and are generated from the unsteady pressure field caused by the "streaky structures" that appear in the laminar boundary layers. They are an effect of the free-stream turbulence intensity, as observed by Klebanoff, [START_REF] Klebanoff | Effect of free-stream turbulence on a laminar boundary layer[END_REF]. In Fig. (2.10), taken from Mayle & Schulz, [START_REF] Mayle | The path to predicting bypass transition[END_REF], it is shown the maximum growth rate of laminar fluctuations for different initial free-stream turbulence intensities in the boundary layer along a flat plate. As observed in Langtry, [START_REF] Langtry | A correlation-based transition model using local variables for unstructured parallelized CFD codes[END_REF], it can be noticed that only for the lowest turbulence intensity, T u = 0.9%, the fluctuations undergo a linear growth, while linear amplification is "bypassed" for higher T u. The transition mechanisms mentioned above will be treated in the following, in a more "engineering framework". For the sake of completeness, two additional transition mechanisms are mentioned, following Mayle, [START_REF] Mayle | The role of laminar-turbulent transition in gas turbine engines[END_REF], that will not be treated in this thesis:

• Reverse Transition or flow relaminarization. When the flow is strongly accelerated the vortex lines in the streamwise direction are stretched, the length scales of the turbulence are reduced and the vorticity is intensified, and energy is transferred to smaller scales. This mechanism involves a balance of production, dissipation and convection of turbulence kinetic energy in the boundary layer, as explained in Mayle, [START_REF] Mayle | The role of laminar-turbulent transition in gas turbine engines[END_REF]. Relaminarization can occur at low free-stream turbulence level and for the acceleration parameter K = ν/U 2 (dU/dx) greater than 3 • 10 -6 .

• Periodic unsteady transition is caused by the periodic passing of wakes from upstream airfoils and obstructions. This kind of transition appears to bypass the linear stability stage: turbulent spots are formed and immediately coalesce, grow and propagate downstream to form a fully turbulent boundary layer.

Transition Modeling

This section is dedicated to a review of some of the approaches used to model and predict laminar-to-turbulence flow transition. The list of the models is not exhaustive, but it is meant to highlight the advantages and disadvantages of the approaches used from the 1990's until today. The aim is to understand why the CFD community is moving towards automatic and autonomous RANS approach to model transition in an industrial framework.

DNS and LES

DNS (Direct Numerical Simulations) are performed by solving the full non-linear, time dependent Navier-Stokes equations. They provide a complete space-time history of the flow field. Since 1980s, as discussed in Kleiser, [START_REF] Kleiser | Numerical simulation of transition in wall-bounded shear flows[END_REF], DNS have become a valuable resource for transition modeling, complementing the traditional and experimental approaches. Indeed, they produce extensive set of full-instantaneous and ensemble data that can be used for a better understanding of transition mechanisms. In order to show the potentiality of DNS simulations, DNS results performed by C. Liu et al., [START_REF] Liu | DNS study on physics of late boundary layer transition[END_REF], are presented in Fig. (2.11). These are simulations for boundary layer transition on a flat plate. Through DNS, they showed that turbulence is not generated by vortex breakdown, which is the last stage of transition in classical theory, see Sec.(2.2), but rather by shear layer instabilities near the wall surface that generate small length scales. Their paper, besides challenging the classical transition theory, shows how DNS are precious tools to understand the flow physics, being less expensive and having more power in data acquisition than a classical experiment. The ability of DNS to produce accurate data that can be further used for the calibration and validation of transition models is undeniable. Nevertheless, DNS are constrained by computer resources, not only to produce the data, but also to extract them, and by algorithmic limitations. Highly resolved meshes, especially in the spanwise direction, are needed to assure an accurate prediction. If we consider as an estimate for the total number of grids point N tot ∼ Re 9 4 to resolve the dissipation length scale, as estimated in Rogallo 3 et al., [START_REF] Rogallo | Numerical simulation of turbulent flows[END_REF], these mesh requirements are not at stage where DNS can serve as a practical tool for engineering purposes. Because of the significant computational cost of DNS, LES (Large Eddie Simulations) have become an alternative approach for the resolution of turbulent and transitional flows. Using LES, only the important large scales are resolved , while the effect of small scales (subgrid) on the large scales are modeled. The application of LES for transitional flow dates from 1990s. Compared to DNS, this approach results in a considerable decrease of CPU hours required for the simulation, because they need less resolved grids. For the simulation of boundary layer transition, LES performed on a mesh of 24 × 48 × 24 points, assure results in reasonable agreement with DNS simulations performed on a grid of 144 × 144 × 224 points, Piomelli et al., [START_REF] Piomelli | On the large-eddy simulation of transitional wall-bounded flows[END_REF]. The problem with LES stands in the sensitivity to the Smagorinsky constant that is used to calibrate the subgrid eddy viscosity, Lardeau et al., [START_REF] Lardeau | Large eddy simulation of transitional separated flow over a flat plate and a compressor blade[END_REF]. The constant coefficient Smagorinsky models fail to predict the transition location, because of the turbulent viscosity ν t . This is active in the laminar and transition region, damping the disturbances in the boundary layer and preventing it from transitioning. In this regard, dynamic LES models perform better, qualitatively, because the Smagorinsky constant is computed locally. It is reduced to zero in the laminar boundary layer, allowing the disturbances growth in the boundary layer. Nevertheless, Wall-Modeled LES (WMLES) do not perform as well as DNS quantitatively. For instance, in Sayadi et al., [START_REF] Sayadi | Large eddy simulation of controlled transition to turbulence[END_REF], it is shown that SGS models fail to produce enough subgrid shear stresses to accurately predict skin friction and mean velocity profiles in 3 The Reynolds number in this estimate is based in terms of large eddy characteristic velocity and length scales.

the transition region, in spite of the use of a very fine grids. In order to avoid these SGSlimitations, one should preferred Wall-Resolved LES (WRLES) to WMLES, but, as mentioned in Spalart, [START_REF] Spalart | Hybrid RANS-LES methods[END_REF], this is not a practical approach. Choi et al., [START_REF] Choi | Grid-point requirements for large eddy simulation: Chapman's estimates revisited[END_REF], gave an estimation of the grid points requirements for LES. In order to simulate the turbulent flow over a flat plate of length L x , the number of grid points necessary to resolve the Kolmogorov length scale are proportional to Re Lx : WMLES requires N tot ∼ Re Lx , while WRLES needs N tot ∼ Re 13/7

Lx . With respect to the same Reynolds number definition, the estimation for DNS becomes N tot ∼ Re 37/14

Lx . Spalart, [START_REF] Spalart | Hybrid RANS-LES methods[END_REF], also mentions a possible future development of DNS-LES hybrid approaches. This is the case of the cylinder simulations for Re ∈ [2.5 × 10 5 , 8.5 × 10 5 ] performed by Rodriguez et al., [START_REF] Rodríguez | On the flow past a circular cylinder from critical to super-critical Reynolds numbers: Wake topology and vortex shedding[END_REF], that, even if labeled as LES, the resolution in the separation-induced transition region is at a DNS level. Their very good results, in terms of transition prediction, confirm the need of a boundary layer well-resolved in order to capture the complex physics phenomena that occur within this region. Nevertheless, we are still outside the range of practical applications, being the simulations from Rodriguez run on unstructured grids from 38M up to 105M, depending on the Reynolds number. Other than the stringent grid requirements, another limitation of both of these approaches is the specification of boundary conditions and initial solution, i.e. base flow, regular perturbations and background noise. Transition is an initial value problem and different value of free-stream noise can strongly affect the transition predictions. The inlet conditions are not trivial, and waves that respect the profile of the disturbances related to the transition process should be imposed.

e N Method One of the most popular approach to predict transition is the e N method, developed by Smith & Gamberoni, [5], and Van Ingen, [START_REF] Van Ingen | A suggested semi-empirical method for the calculation of the boundary layer transition region[END_REF]. An historical review can be found in Van Ingen, [START_REF] Van Ingen | The e N method for transition prediction. historical review of work at TU Delft[END_REF]. The e N method makes use of the spatial problem, where the frequency ω ∈ R and the wave number α ∈ C. This spatial setting does not change the general form of the equations as presented in Sec.(2.1). Nevertheless, the spatial stability problem is given by an eigenvalue problem were now the eigenvalue α appears nonlinearly up to the 4th power. The N factor is the total growth rate of the most unstable disturbance of frequency f at the transition location x. For 2D, incompressible flows, N is defined as N = log( A A 0 ) max = max f -x x 0 α i dx , where x 0 is the x-coordinate where the perturbation enters the unstable zone, α i is the local spatial amplification rate. A and A 0 are the amplitudes at x and x 0 , respectively. N is the envelope of curves formed by the change of the total amplification rate for a variety of perturbations with different frequencies. Transition occurs when N reaches its critical value. The critical N c factor at which transition occurs is computed solving the incompressible linearized stability equations and then calibrating with the experimental conditions, free-stream environment, wall surface roughness. For T-S waves in low free-stream turbulence environment, Mack in 1984 proposed the empirical formula for the critical N factor:

N c = -8.43 -2.4 ln(T u).
(2.30)

The experimental definition of the critical N factor makes the model semi-empirical. A valid and clear set of experimental data is needed for every new flow configurations, because the stability limit is specific to the case. It is built on the basis of experimental data, with respect to the free-stream turbulence intensity and surface roughness. e N method are used as a part of transition prediction module along with a boundary layer code. In practice, RANS are used to compute the fully laminar flow around the configurations of interest until convergence of the lift force. The obtained pressure distribution is used as input data for the boundary layer code, which provides the viscous data that are needed by the stability code: determination of the boundary layer edge velocities, calculation of streamlines, which are mapped to the surface, and extraction of BL profile along the streamlines. The calculation of these data requires very dense grids. Based on the stability analysis results, which provide the estimated frequencies and wave lengths, the local N T S and N CF are computed for each streamlines. The transition location is determined along each streamline and communicated back to the RANS solver. The turbulence model is then activated downward the transition front. e N method is a useful physics-based tool for predicting the position of boundary layer transition, nevertheless the need to determine the N factor using semi-empirical method represents its biggest limitation. In Krumbein, [START_REF] Krumbein | e N transition prediction for 3D wing configurations using database methods and a local, linear stability code[END_REF], this problem is clearly discussed. It is explained that the application of the e N method for free flight configurations is justified by the extensive flight test campaigns runs that have allowed the determination of the ranges of values of the N factors for T-S and CF instabilities. However, when simulating wind tunnel flows, the N factors strongly depend on the wind tunnel specific characteristics and must be determined for each wind tunnel configuration. The determination of the N factor, thus, requires sufficient experimental data. We often do not dispose of these data, or they cannot be clearly interpreted or they are even contradictory.

PSE

A more general approach than the application of linear stability equations, are the parabolized stability equations (PSE). This approach allows the direct description of the disturbances growth in the boundary layer and their impact on the base flow. This set of equations is obtained through the parabolic approximation of Navier-Stokes equations in the streamwise direction. The disturbances are introduced by amplitude functions which vary in the x -zplanes, as well as the base flow. Supposing that the variations in the streamwise directions are weaker than those in the wall normal direction, non parallel effect are taken into account. A non-linear version of PSE can also be used for transition prediction, allowing the calculation of the disturbances growth in the non-linear stage. For their practical application, the main challenge stands in the definition of the initial amplitude and phases of the various modes which have to be imposed at the boundary. They have to be estimated from the initial disturbances environment under appropriate assumptions on the receptivity process. Non-linear PSE can reproduce DNS results at a lower cost, as discussed in Bagheri et al., [START_REF] Bagheri | The stabilizing effect of streaks on Tollmien-Schlichting and oblique waves: A parametric study[END_REF], however non linear stability approaches are time consuming and hard to use in practical applications, Zang et al., [START_REF] Zang | The transition prediction toolkit: LST, SIT, PSE, DNS and LES[END_REF].

Low Reynolds turbulent closure approach

We underlined various disadvantages of the transition prediction tools described above, the stringent mesh requirements, the elevated computational costs and their problematic use in practical applications. Moving to a RANS frame seems a natural choice, because of the simplicity of models implementation. The interest does not lie in models that are capable of describing transition dynamics, but in approaches able to predict the statistical characteristics of transitional flows. In a first attempt, in Rodi & Scheuerer, [START_REF] Rodi | Calculation of heat transfer to convection-cooled gas turbine blades[END_REF], the capability of the RANS low-Re versions of the k -turbulence model, by Jones and Launder, [START_REF] Jones | The calculation of low-Reynolds-number phenomena with a two-equation model of turbulence[END_REF], to account for transition was tested. These models were found to perform fairly well for bypass transition, despite the fact that they were not calibrated for transitional flows. Low-Re k -models rely on the diffusion of the fluctuations turbulence kinetic energy from the free-stream into the boundary and its interaction with the model source terms to predict transition, Langtry, [START_REF] Langtry | A correlation-based transition model using local variables for unstructured parallelized CFD codes[END_REF]. It was soon observed that their ability to predict transition was coincidental, as discussed in Schmidt & Patankar, [START_REF] Schmidt | Simulating boundary layer transition with low Reynolds number k -e turbulence models. I-An evaluation of prediction characteristics. II-An approach to improving the predictions[END_REF]. Indeed, transition predictions by these RANS models was due to the similarities between the viscous sublayer and the developing laminar boundary layer. In addition, the production of turbulence kinetic energy is bounded through limiters based on stability criteria and correlated to the free-stream turbulence intensity, and it is extremely sensitive to the type and implementation of the initial and boundary conditions. Their behavior mimics bypass transition, without any other modification of these models implementation specific to transition, however low-Re turbulence models cannot be sensitive enough to other transition mechanisms, as discussed in Zheng, [START_REF] Zheng | Turbulent transition simulation using the kω model[END_REF]. The need of a proper calibration of these turbulence models for transition was further reproposed by [START_REF] Rumsey | Arbitrary steady-state solutions with the k-epsilon model[END_REF] and Patterson-Reif (2006), [START_REF] Rumsey | Arbitrary steady-state solutions with the k-epsilon model[END_REF]. They showed that different numerical strategies and mesh densities and other numerical parameters could lead to different converged transition solutions. They concluded that without a proper calibration these models cannot predict transitioning flow behavior with its specific flow dynamics.

k -k l Model

An alternative to k -and k -ω RANS turbulence models is the k -k l model. It is based on the idea that transition dynamics can be modeled accounting for the development of the pre-transitional laminar fluctuations. These are considered through the quantity k l , the laminar kinetic energy, that is the magnitude of these streamwise fluctuations. The idea was proposed at first by Mayle & Schulz, [START_REF] Mayle | The path to predicting bypass transition[END_REF], as already discussed in Sec.(2.2). Later, it was improved by Walters & Leylek, [START_REF] Walters | A New Model for Boundary Layer Transition Using a Single-Point RANS Approach[END_REF], [START_REF] Walters | Computational Fluid Dynamics Study of Wake-Induced Transition on a Compressor-Like Flat Plate[END_REF], who developed a transport equation for the laminar kinetic energy to overcome the issue of using non-local variables. In the k -k l model the equation for the laminar turbulence kinetic energy is added to the equations for turbulent kinetic energy k and dissipation rate ω. k and ω equations are not identical to the standard k -ω models. k -k l does not include specific transition features, besides the definition of ν t,l , that is the large scale turbulent viscosity used in the production term of the laminar kinetic energy. In the Reynolds stresses terms, the total eddy viscosity is given by the sum of the small-scale and large scale eddy viscosity. This model performs well for bypass transition, Choudhry et al., [START_REF] Choudhry | A study of long separation bubble on thick airfoils and its consequent effects[END_REF], as expected. Nevertheless, it requires larger computational resources and time, compared to models calibrated to predict transition, and does not account for natural transition mechanisms, as discussed in Aftab et al., [START_REF] Aftab | Turbulence model selection for low Reynolds number flows[END_REF].

Intermittency models

Intermittency models have known a huge success in the last twenty years. This approach is based on the definition of an empirical correlation for the intermittency, γ. A peculiar characteristic of a transitional flow is its intermittent behavior. Turbulent spots appear intermittently in the transition process and are convected downstream in the boundary layer, see Fig. (2.4). γ is the fraction of time during which the flow is turbulent: if γ=0 the flow stays laminar, for γ = 1 the flow is fully turbulent. Based on this definition in 1958, Dhawan & Narasimha, [START_REF] Dhawan | Some properties of boundary layer flow during the transition from laminar to turbulent motion[END_REF], developed the empirical correlation of Eq. (2.31). This describes the streamwise evolution of the intermittency, based on the transition point onset x t , the free stream velocity U , the formation rate of the turbulent spots n and their propagation rate σ. This empirical correlation is often used as a basis for the development of intermittency models.

γ =    1 -exp --(x-xt) 2 nσ U , (x x t ) 0, (x x t ) (2.31)
The intermittency γ factor can be used in two different frames:

• Algebraic intermittency models, the algebraic definition of the intermittency factor is incorporated in the turbulence model through its multiplication by the eddy viscosity. The intermittency factor modulates the turbulence. Upstream the transition location γ = 0 and ν t = 0, and a laminar boundary layer develops. As the transition criterion is met, γ starts increasing as the eddy viscosity. Defined the pressure gradient parameter λ θ as:

λ θ = θ 2 ν dU ds , (2.32)
where θ is the momentum thickness, and dU ds is the velocity derivative along the streamlines direction.

The transition onset location is based on an additional empirical correlation, which correlates the T u and λ θ to the transition momentum thickness Reynolds number, like the correlation of Abu-Ghannam & Shaw, [START_REF] Abu-Ghannam | Natural transition of boundary layers-the effects of turbulence, pressure gradient, and flow history[END_REF], Mayle, [START_REF] Mayle | The role of laminar-turbulent transition in gas turbine engines[END_REF] and Schiele, [START_REF] Schiele | Calculating turbulent and transitional boundary-layers with two-layer models of turbulence[END_REF]. The use of intermittency along with an empirical correlation based on the momentum thickness provides information about the flow physics in the transition region, that would be otherwise condensed in one point (x c x t ).

• Intermittency transition models, where the intermittency is the solution of an additional transport equation, where the source terms mimic the behavior of algebraic intermittency models. The advantage of this approach is that it provides the intermittency distribution across the boundary layer and not only in the flow direction, leading to more accurate transition prediction. In this frame, the local correlation based transition models were conceived.

LCTM

The LCTM framework was first proposed by Langtry in his doctoral thesis, [START_REF] Langtry | A correlation-based transition model using local variables for unstructured parallelized CFD codes[END_REF]. In his dissertation, Langtry clearly summarized the main requirements for a transition model to be CFD-compatible. They are hereafter quoted:

• Be formulated locally (no search or line-integration operations);

• Allow the calibrated prediction of the onset and the length of transition;

• Allow the inclusion of different transition mechanisms;

• Avoid multiple solutions (same solution for initially laminar or turbulent boundary layer);

• Do not affect the underlying turbulence model in fully turbulent regimes;

• Allow a robust integration down to the wall with similar convergence as the underlying turbulence model;

• Applicable to three-dimensional boundary layers;

• Be formulated independent of the coordinate system.

The necessity to avoid non-local information in correlation-based models was first addressed by Menter et al., [START_REF] Menter | Transition modelling based on local variables[END_REF], who proposes the use of quantities constructed locally to trigger the production of intermittency in the boundary layer. Within this new formulation, the link between the correlations and the intermittency equation is achieved through the introduction of the vorticity Reynolds number, Re V . This quantity is used to trigger the transition onset, as first proposed by Van Driest & Blumer, [START_REF] Van Driest | Boundary layer transition-freestream turbulence and pressure gradient effects[END_REF].

Re V measures the strength of the local inertial stresses to the local viscous stresses:

Re V = τ turb τ visc = -ρu v -µ du dy = -ρl 2 ( du dy ) 2 -µ du dy = ρl 2 µ du dy ∼ ρy 2 S µ . (2.33)
l is the scale length of the turbulence eddies, which is assumed to be proportional to the wall distance y. The derivative of the velocity in the y direction is expressed through the mean strain rate S in the boundary layer. Langtry & Sjölander , [START_REF] Langtry | Prediction of transition for attached and separated shear layers in turbomachinery[END_REF], showed that the peak growth of disturbances in the boundary layer coincides with the maximum of Re V . Physically, as stated in Menter & Langtry, [START_REF] Menter | Transition modelling for turbomachinery flows[END_REF], y 2 S is responsible for the growth of disturbances inside the laminar boundary layer, while µ of their damping. For constant µ, as the laminar boundary layer thickness grows with y 2 S, transition will take place when a critical value of Re V is reached. The scaled profile of the vorticity Reynolds number for a Blasius boundary layer is shown in Fig. (2.13). The scaling is chosen in order to have a maximum of one in the boundary layer. It is achieved dividing Re V by the corresponding momentum thickness Reynolds number and a calibration constant. In the two LCTM models, γ -Re θ , Menter & Langtry, [START_REF] Menter | Transition modelling for turbomachinery flows[END_REF] and γ by Menter et al., [START_REF] Menter | A one-equation local correlation-based transition model[END_REF], the vorticity Reynolds number enters the intermittency transport equation through its ratio to the critical momentum thickness Reynolds number Re θc . This critical ratio is used to trigger transition onset in the streamwise direction. The vorticity Reynolds number is a local property which can be easily computed at each grid point, also in an unstructured solver. The critical momentum thickness Reynolds number Re θc is computed as an empirical correlation. γ -Re θ is based on the solution of two transport equations in addition to the turbulence equations: one for the intermittency γ and one for the momentum Reynolds number Re θt . The intermittency is used to turn on the production of turbulence kinetic energy downstream the transition onset point within the boundary layer. γ is set to be equal to one in the free-stream, in order not to interfere with the turbulence model in the free-stream. A different value of intermittency in the free-stream would affect the turbulence decay ahead of the body, at the edge of the boundary layer, and near the stagnation point. The equation for Re θt is needed to transport non-local information, as turbulence intensity T u and non-dimensional pressure gradient effects, from the free-stream inside the boundary layer in order to trigger the transitional correlation. As the free-stream T u enters the empirical correlations, the model is not Galilean invariant. The one equation local correlation transition model, γ model, was proposed in 2015. The transport equations are reduced from two to one, simplifying drastically the formulation and Galilean invariance is achieved. In this new formulation, the critical momentum thickness Reynolds number results from an empirical correlation based on local quantities: the local turbulence intensity and the local pressure gradient parameter. Both transition models were calibrated on the ERCOFTAC T3 series of flat plate experiments, [START_REF] Coupland | Ercoftac classic database[END_REF], and the Schubauer and Klebanoff flat plate experiments, [START_REF] Schubauer | Turbulent boundary layer[END_REF]. The T3 series has high freestream turbulence intensities (bypass transition), while Schubauer and Klebanoff are for natural transition. It is important to underline that in both models the transport equations for the transition quantities are not meant to describe transition physics. Physical information are instead contained in the empirical correlations incorporated in the models' formulation, used to trigger the transition process. This approach allows the inclusion of the different transition mechanisms through the definition of additional specific correlations. The original models, as published in the two mentioned articles by Menter, cover for both bypass transition and streamwise transition at low free-stream turbulence levels.

Re V /(2.2Re θ ) η 0 0.2 0.4 0.6 0.

AFT

To conclude, the AFT model from Coder & Maughmer, [START_REF] Coder | Application of the amplification factor transport transition model to the shear stress transport model[END_REF], [START_REF] Coder | Further development of the amplification factor transport transition model for aerodynamic flows[END_REF], solves an additional equation for the amplification factor ñ, that characterizes the collection of the instabilities linearly amplified within the boundary layer. The basic idea is to use locally computed surrogates of integral values, such as the local shape factor H 12 , which is the ratio of the streamwise displacement thickness to the streamwise momentum thickness. The purpose is to define a physical model based on fundamental boundary-layer theory. The local shape factor is based on the wall normal gradient of the wall normal momentum. As for e N stability methods, transition onset is controlled through the critical parameter N c evaluated as a function of the free-stream turbulence intensity. It is expressed based on Mack's relation with an additional modification and it reads:

N c = -8.43 -2.4 ln τ 100 , (2.34)
τ = 2.5 tanh T u(%) 2.5 (2.35)
The equation for ñ, which does account for any transition criteria for transition mechanisms other than natural transition, is coupled to a turbulence model, either of k -ω family or Spalart-Allmaras. The production and/or the destruction terms in the turbulence equations are modified, in order to balance the production of turbulence variables in the laminar boundary layer. The strength of the models is its reduced sensitivity to inlet conditions, namely the eddy viscosity ratio, as discussed in Lopes et al., [START_REF] Lopes | Assessment of rans transition models[END_REF]. Nevertheless, efforts are still ongoing to make the model Galilean invariant and to include 3D transition mechanisms, as crossflow.

To our best knowledge, a crossflow criterion for AFT has not been published yet. The main difficulty relies in the definition of the critical N c factor for this transition mechanism, for which no universal empirical definition exists. This critical factor enters the model within the function that controls the suppression of turbulence production in the laminar region in the transport equation for k.

Main Remarks and Objectives of the Research

A brief general introduction to the basic concept of linear stability theory served the purpose to present the main laminar-to-turbulence transition mechanisms. Transition to turbulence is not a unique process, but it depends on the initial disturbances and the type of flow under analysis. Depending on the free-stream turbulence intensity, we can distinguish two main classes of transition mechanisms: natural transition, as Tollmien-Schlichting (T-S) waves or crossflow (CF) vortices, and bypass transition. A free-stream turbulence of 1% is often taken as the boundary between these two mechanisms. The transition scenario depends on whether or not the instabilities exponential growth (natural transition) is able to compete with the growth of streamwise streaks within the boundary layer (bypass transition). Natural or bypass transition, vortex breakdown and the generation of turbulence spots are generic processes that do not strongly depend on the shear flow under consideration. The first spot appearance coincides with the increase of skin friction, which deviates from its laminar profile, until a fully developed turbulence state is attained. The position at which the skin friction starts to increase is the transition position, at which the transition momentum thickness Reynolds number exceeds its critical value. In practical and industrial applications, a model able to predict transition position and characteristics of the flow in the transition region is fundamental. For instance, separation induced transition can lead to the appearance of laminar separation bubbles (LSB). LSB have degrading effects on aerodynamic performance, and models able to predict them accurately (prediction of separation, transition and reattachment point) are needed in order to control, and eventually eliminate them. Laminar flow is less stable and more prone to separation than turbulent flow. Thus, the most effective method to avoid the flow separation is to force the flow to become turbulent, for instance with turbulators on laminar wing. Several transition modeling strategies have been described above. DNS and LES are approaches that allow the characterization of transition physics. Especially, DNS is a valuable resource to understand transition mechanisms in lack of experimental results. However, with the intention of applying transition models to industrial application, for Reynolds number ranging from 10 5 up to 10 7 , both DNS and LES are too computationally expensive. For this reason, CFD model developers have started to think about transition prediction/modeling in the frame of RANS methods. RANS formulation is obtained by averaging the Navier-Stokes equations and the effect of linear amplification of the disturbances is suppressed. Nevertheless, the purpose of modeling transition in a RANS framework is to account for transitional effects on the quantities that describe the flow physics, such as pressure and skin friction distribution. From this perspective models conceived in a RANS environment can be used for transition, upon proper calibration. A RANS transition model should require only local quantities calculations, it should be easily implemented in unstructured and parallelized CFD solvers and should allow the integration of the different transition mechanisms. This is the philosophy of the LCTM concept proposed by Langtry in his doctoral thesis, [START_REF] Langtry | A correlation-based transition model using local variables for unstructured parallelized CFD codes[END_REF], that has led to the development of γ -Re θ and γ transition models. This thesis focuses on the analysis of the performance of these two RANS transition models under modeling and numerical points of view. 2D and 3D flows transitional flows are treated in the following discussion. The mathematical formulation of γ -Re θ and γ is given in the next chapter.

RANS: Reynolds-averaged Navier Stokes equations

RANS (Reynolds Averaged Navier Stokes) equations are the approach commonly used in industrial codes. The set of equations, resulting from averaging the Navier-Stokes equations, Eq.(2.2)-(2.3), aims to represent the action of turbulence on the average flow through the least expensive statistical approach. Given a generic scalar quantity φ, an averaging operator has the following property:

φ = φ. (3.1)
Reynolds averaging can involve a variety of forms. For the turbulence modeling, the most pertinent is the time average, because of the in-homogeneity of the flows commonly involved in engineering problems. The time average of a quantity φ(x, t) is defined by:

φ = lim T -→∞ 1 T t+T t φ(x, t)dt. (3.2)
With such a time average, the quantity φ does not depend on time. Once an averaging operator is defined, φ(x, t) can be decomposed as:

φ(x, t) = φ + φ (x, t), (3.3)
where φ (x, t) represents the fluctuating part. The decomposition in Eq.(3.3) can be used for the instantaneous velocity components, u i , and the pressure, such that:

u i = u i + u i and p = p + p . (3.4)
The resulting averaged Navier-Stokes equations for each component u i are:

∂u i ∂x i = 0, (3.5
)

∂u i ∂t + u j ∂u i ∂x j = - 1 ρ ∂p ∂x i + ∂ ∂x j ν ∂u i ∂x j + ∂u j ∂x i + ∂ ∂x j (-u i u j ). (3.6)
An unclosed term appears from the averaging, the Reynolds stress tensor, whose components are τ ij = -u i u j . This tensor quantifies the influence of the velocity fluctuations over the mean flow. The primary problem with RANS is that there are more unknowns (pressure, three velocity components, six components of the symmetric tensor τ ij ) than equations and additional relations are needed to close the system. The mathematical models providing an implicit or explicit link between τ ij and the flow quantities are said turbulence models. This link can be achieved either solving an equation for each component of the Reynolds stress tensor, either seeking to construct an expression for τ ij . This is the case of the Boussinesq's hypothesis, which assumes the existence of a linear relationship between the Reynolds stress tensor and the strain rate tensor

S ij = 1 2 ∂u i ∂x j + ∂u j ∂x i .
For isotropic linear turbulence closures, we have the following definition of the Reynolds stress tensor components:

-u i u j = ν t ∂u i ∂x j + ∂u j ∂x i - 2 3 kδ ij , (3.7)
where k = 1 2 u i u j and ν t is the turbulent viscosity. In the following we are going to present the formulation of the so-called k-ω SST model from Menter, based on Boussinesq's hypothesis. The k -ω framework is based on the resolution of the turbulence kinetic energy, k, and the turbulence frequency, ω, in order to evaluate the turbulence eddy viscosity ν t , that appears in Eq.(3.7).

Turbulence Model: k -ω SST

The k -ω Shear Stress Transport (SST) two equations eddy-viscosity model is a blend between k -and k -ω formulation: k -formulation is used in the free-stream, in order to avoid the dependence of k -ω on the free-stream turbulence intensity. The model switches to the k -ω formulation near the wall, because of its superior behavior within the viscous sublayer. A first version of the model was proposed in 1994 by Menter, [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF]. Hereafter, we are going to refer to the version of 2003 by Menter et al., [START_REF] Menter | Ten years of industrial experience with the sst turbulence model[END_REF]. The turbulence model solves two different equations, one for the turbulence kinetic energy k:

∂(ρk) ∂t + ∂ ∂x j ρU j k -(µ + σ k µ t ) ∂k ∂x j = τ ij S ij -β * ρωk, (3.8)
and one for the turbulence frequency, ω:

∂(ρω) ∂t + ∂ ∂x j ρU j ω -(µ + σ ω µ t ) ∂ω ∂x j = γρΩ 2 -βρω 2 + 2(1 -F 1 ) ρσ ω2 ω ∂ω ∂x j ∂k ∂x j . (3.9)
The blending between the two formulations is achieved through the function F 1 , which reads as:

F 1 = tanh min max 2 √ k 0.09dω , 500 µ ρd 2 ω , 4ρσ ω2 k CD kω d 2 4 , (3.10)
where

CD kω = max 2ρσ ω2 ω ∂k ∂x j ∂ω ∂x j .10 -20 .
(3.11)

LCTM

The function F 1 operates blending the model constants φ = {β, γ, σ k , σ ω }. Given two constants φ 1 and φ 2 , where the subscript indicates which formulation is used, the blend is achieved through the definition:

φ = F 1 φ 1 + (1 -F 1 )φ 2 .
(3.12)

F 1 controls also the activation of the cross-diffusion term in ω equation:

2(1 -F 1 ) ρσ ω2 ω ∂ω ∂x j ∂k ∂x j . (3.13)
The latter term results from the derivation of the equation for ω from the equation, given the relation ω = /kβ * . The eddy viscosity is defined as:

µ t = ρk/ω max{1, ΩF 2 /(a 1 ω)} , (3.14)
where the function F 2 is defined as follows:

F 2 = tanh max 2 √ k 0.09dω , 500 µ ρd 2 ω , (3.15)
with d is the distance to the wall. The constants are:

• for the inner model:

σ k1 = 0.85, σ ω1 = 0.5, β 1 = 0.075, γ 1 = 0.555,
• for the outer model:

σ k2 = 1, σ ω2 = 0.5, β 2 = 0.440.
Finally, a 1 = 0.31 and β * = 0.09. Finally, the boundary conditions are:

• at the wall: k = 0 and ω = 10 • 6µ/β 1 ρ(∆y) 2 ;

• as outer limit :

ω ∞ = λ U ref L ref , µ t∞ = 10 -3 µ and k ∞ = µ t∞ ω∞ ρ
, where λ assumes values between 1 and 10. ∆y is the distance of the first point away from the wall and such that y + < 1.

LCTM

γ -Re θ Transition Model γ Equation

The intermittency transport equation is:

∂(ργ) ∂t + ∂(ρu j γ) ∂x j = P γ -E γ + ∂ ∂x j µ + µ t σ γ ∂γ ∂x j . (3.16)
The transition source is defined as:

P γ = c a1 F length ρS(γF onset ) 0.5 (1 -c e1 γ), (3.17)
where S is the mean strain magnitude. F length is the empirical correlation which controls the transition length, i.e. the magnitude of the source term, and depends on Re θt , solution of the second transport equation. P γ ∼ 0 until the local vorticity Reynolds number exceeds the critical momentum thickness Reynolds number, as determined through F onset .

The function F onset which triggers the transition is defined as follows:

Re V = ρy 2 S µ , (3.18)
F onset,1 = Re V 2.193Re θc , (3.19)
F onset,2 = min(max(F onset,1 , F 4 onset,1 ), 2), (3.20)
F onset,3 = max 1 - R T 2.5 3 , 0 , (3.21)
F onset = max(F onset,2 -F onset,3 , 0), (3.22)
where

R T = ρk ωµ (3.23)
is the eddy viscosity ratio. Dividing the vorticity Reynolds number by 2.193Re θc ensures Re V to have a maximum of one within the boundary layer. F onset,2 assures a rapid change from zero to one as the onset transition criterion F onset,1 > 1 is met. As mentioned before, the production of γ promotes the building up of turbulence through the activation of the production term of turbulence kinetic energy, k. F onset,3 in Eq.(3.21) is meant to keep active the production term of γ throughout the transition process, as the eddy viscosity ratio R T increases. Indeed, during the transition process Re V can diminish, because of the change in the velocity profiles, causing the transition process to stall.

Re θc = ρθU µ c
is the critical Reynolds number at which the intermittency starts to increase and turbulence is built up. It is found upstream than Re θt and their link is obtained through an empirical correlation, that will be discussed in the following. The last term in Eq.(3.17) is needed to limit the maximum value of γ, in order for the intermittency not to exceed one. This limiting value depends on the constant c e1 . Finally, the term ρS is used to achieve the right units for the source term. The destruction/relaminarization term is defined as:

E γ = c a2 ρΩF turb γ(c e2 γ -1), (3.24)
where Ω is the vorticity rate magnitude. As explained in the thesis of Langtry, Ω is preferred to S in the destruction of γ in order to avoid the destruction of intermittency due to free-stream strain rates. E γ works as a sink, forcing the intermittency to be equal to its minimum in the laminar boundary layer. The constant c e2 enforces the lower bound for γ. The function F turb is meant to deactivate the destruction/relaminarization term outside of the laminar boundary layer. It is defined as:

F turb = e -R T 4 4
.

(3.25)

The constants of the equations obtained by numerical calibration are:

c a1 = 2, c e1 = 1, c a2 = 0.06, c e2 = 50, σ γ = 1. (3.26)
The boundary conditions for γ are zero normal flux at the wall and 1 at the inlet. γ is imposed to one in the free-stream, in order to avoid any contamination of the free-decay of turbulence variables in the free-stream. The constant c e2 guarantees a lower bound on γ of 0.02. This enforcement has a small impact on the k -ω equations in the viscous sub-layer, indeed such a small value of γ leads to a negligible production of turbulence.

LCTM

Re θt Equation

The transport equation for the transition momentum thickness Reynolds number Re θt is:

∂(ρRe θt ) ∂t + ∂(ρu j Re θt ) ∂x j = P θt + ∂ ∂x j σ θt µ + µ t ∂Re θt ∂x j . (3.27)
The source term is designed to enforce Re θt to be equal to the local value of Re θt obtained from an empirical correlation. It is defined as follows:

P θt = c θt ρ t (Re θt -Re θt )(1 -F θt ), (3.28)
where

t = 500µ ρU 2 , (3.29)
F θt = min max F wake e -y δ 4 , 1 - γ -c e2 1 -1/c e2 2 , 1 , (3.30)
θ BL = Re θt µ ρU , δ BL = 15 2 θ BL , δ = 50Ωy U δ BL , (3.31
)

Re ω = ρωy 2 µ , F wake = e -Reω 1e+5 2 . (3.32)
where t is a time scale, needed for dimensional reasons, with U velocity norm. It scales with the convection and diffusion term of the transport equation. The function F θt is equal to one, i.e. P θt ∼ 0, inside the boundary layer in order to allow the scalar Re θt to diffuse in the boundary layer from the free-stream. It is zero in the free-stream. The function F wake is introduced in order to deactivate the blending function F θt in the wake regions downstream the body. The constants obtained by numerical optimization are:

c θt = 0.03, σ θt = 2. (3.33)
c θt controls the magnitude of the source term, σ θt of the diffusion term. For a value of σ θt higher than two, the transition model becomes less sensitive to the history effect, the past history of pressure gradient and turbulence intensity. The value of two was obtained based on the T3 flat plate series, where history effects are supposed to be significant. The boundary conditions for Re θt are zero flux normal at wall. Its value at the inlet depends on the empirical correlation reported in the coming paragraph, and it depends on the inlet turbulence intensity.

Empirical Correlations

The model presents three different empirical correlations:

• one correlation for the transition momentum thickness Reynolds number Re θt , based on the turbulence intensity and the pressure gradient;

• one correlation for F length as a function of Re θt , numerically calibrated on simulations of the flow on the flat plate and constructed by curve fitting with respect to the experimental data;

• a final correlation for the critical Reynolds number Re θc based on Re θt , obtained by numerical experimentation on a flat plate.

The empirical correlation for Re θt is based on the Thwaithes pressure gradient coefficient:

λ θ = ρθ 2 µ dU dS , (3.34)
where dU dS is the velocity derivative along the streamwise direction:

dU dS = u i u j U 2 ∂u i ∂x j , (3.35)
and the local turbulence intensity T u:

T u = 100 2k/3 U . (3.36)
The local transition momentum thickness Reynolds number Re θt is defined as:

Re θt = 1173.51 -589.428T u + 0.2196 T u 2 F (λ θ ), if T u 1.3, Re θt = 331.50(T u -0.5658) -0.671 F (λ θ ), if T u > 1.3, (3.37)
where

F (λ θ ) = 1 -(-12.986λ θ -123.66λ 2 θ -405.689λ 3 θ )e -T u 1.5 1.5 , if λ θ 0, F (λ θ ) = 1 + 0.275(1 -e -35λ θ )e -T u 0.5 , if λ θ > 0.
(3.38)

For numerical robustness the following limitations should be imposed:

-0.1 λ θ 0.1, T u 0.0027,
Re θt 20.

(3.39)

For T u = 0.0027, the transition momentum thickness Reynolds number is equal to 1450, which is the largest experimental value observed in experiments around a flat plate. Being Re θt = U θt ν , its calculation is done iteratively, because θ t appears in both right and left hand sides.

The function F length controls the transition length and it is defined as: In order to avoid any nonphysical skin friction increase shortly after transition F length is imposed to be equal to its maximum value in the viscous sub layer, F length = 40, through the following functions:

F length =                  [39.8189 + (-0.011927Re θt ) + (-0.000132567Re θt 2 )], if
F sublayer = e -Rω 0.4 2 , Re ω = ρωy 2 500µ , (3.41)
F length = F length (1 -F sublayer ) + 40F sublayer . (3.42)
This redefinition is needed because of the decrease of Re θt in the boundary layer, for instance, at high Reynolds numbers. Indeed, for low values of Re θt , F length increases locally, causing the increase of the intermittency source term. This behavior results in a sharp increase of the skin friction. As mentioned in Menter & Langtry, [START_REF] Menter | Transition modelling for turbomachinery flows[END_REF], this is caused by a sharp change in the y + in the viscous sublayer where the intermittency is forced to be equal to its minimum by the destruction term.

The correlation between the critical transition momentum Reynolds number Re θc and Re θt is defined as follows:

Re θc =        Re θt - 3
γ sep = min s 1 max 0, Re V 3.235Re θc -1 F reattach , 2 F θt , (3.44)
where

F reattach = e -R T 20 4
, s 1 = 2.

(3.45)

F reattach deactivates the term when the flow reattaches. γ sep accounts for the fact that Re V significantly exceeds Re θc when laminar separation occurs. The flow reattachment location is controlled by the constant s1, which is set to be equal to two. In order to account for the flow separation, in the coupling with the turbulence model, the intermittency appears in the turbulence kinetic energy equation as:

γ ef f = max(γ, γ sep ). (3.46)
Through the effective intermittency we allow the intermittency to be bigger than one, when separation-induced transition occurs. It boosts the production of turbulence kinetic energy downward the laminar separation point.

Coupling with k -ω SST turbulence model

The transition model interacts with the turbulence model through an adaptation of the production and destruction terms of the turbulence kinetic energy transport equation, as follows:

∂(ρk) ∂t + ∂(ρu j k) ∂x j = Pk -Dk + ∂ ∂x j µ + µ t σ k ∂k ∂x j , (3.47)
Pk = γ ef f • P k , (3.48) Dk = min(max(γ ef f , 0.1), 1)D k .
(3.49)

The production term P k is computed using the Kato-Launder modification, [START_REF] Kato | The modelling of turbulent flow around stationary and vibrating square cylinders[END_REF], in order to reduce the excessive production of turbulence kinetic energy in regions with large normal strain:

P k = µ t SΩ. (3.50)
D k is the original destruction term of the SST model, Eq.(3.8). The equation for the turbulence frequency ω, Eq.(3.9), is unchanged. Last modifications are needed to redefine the blending function F 1 , which is responsible of the switching from k -ω to k -. The new blending function is forced to be equal to one in the laminar region and it is redefined as:

F 1 = max(F 1orig , F 3 ), (3.51)
F 3 = e ( Ry 120 ) 8 , R y = ρy √ k µ , (3.52)
where F 1orig is the definition of F 1 in the SST model, Eq.(3.10).

γ Model γ Equation

The one transport equation for γ reads as:

∂(ργ) ∂t + ∂(ρu j γ) ∂x j = P γ -E γ + ∂ ∂x j µ + µ t σ γ ∂γ ∂x j . (3.53)
The transition source term is defined as:

P γ = F onset [F length (ρS(1 -γ)γ)], (3.54)
The constant F length , which is not an empirical correlation, as it was for γ -Re θ , is set equal to 100. F onset is defined as follows:

F onset1 = Re V 2.2Re θc , (3.55)
F onset,2 = min(F onset1 , 2.0), (3.56)
F onset,3 = max 1 - R T 3.5 3 , 0 , (3.57)
F onset = max(F onset,2 -F onset,3 , 0). (3.58)
The critical Reynolds number Re θc is computed as an empirical correlation , that depends on local values of turbulence intensity T u L and pressure gradient λ θ,L . As in γ -Re θ the local vorticity Reynolds number is formulated as a function of the strain rate S, Eq.(3.18), and the various F onset functions have the same role, even if formulated differently. The relaminarization/destruction term depends on the magnitude of the absolute vorticity rate and it is defined as:

E γ = c a2 ρΩγF turb (c e2 γ -1), (3.59)
where the definition of F turb has slightly changed and it reads as:

F turb = e -( R T 2 ) 4 .
(3.60)

LCTM

The model constants obtained by numerical calibration are:

c e2 = 50, c a2 = 0.06, σ γ = 1. (3.61)
The boundary conditions are zero normal flux at the wall and γ equal to one at the inlet. As in γ -Re θ , the equation enforces γ = 0.2 in the viscous sub-layer.

Local formulation of T u L and λ θ,L .

The major difference between γ -Re θ and γ is the definition of the critical Reynolds number. Re θc is now computed as an empirical correlation Re θc = f (T u L , λ θL ). T u L and λ θL are calculated inside the boundary layer using local quantities. Their estimation through local quantities avoids the transport equation for the Re θt and makes the model Galilean invariant.

The turbulence intensity within the boundary layer is expressed as:

T u L = min 100 2k/3 ωy , 100 . (3.62)
The expression is similar to the definition of free-stream turbulence intensity. y is the distance from the wall and ω the specific turbulence dissipation rate. ωy provides a velocity scale in the boundary layer, indeed, U ∼ Sy ∼ ωy.

λ θ represents a much more critical parameter to estimate. The pressure gradient parameter is defined as:

λ θ = - θ 2 µ 1 U dP ds = ρθ 2 µ dU dS , (3.64)
where dU dS is the velocity derivative in the streamwise direction. For a flat plate, Eq.(3.64) reduces to:

λ θ = ρθ 2 µ du dx = - ρθ 2 µ dv dy , (3.65)
exploiting the incompressibility for two-dimensional flows. dv dy is the wall-normal derivative of the wall normal velocity component v. For a general geometry, it can be computed as:

dv dy = ∇( - → n • - → u ) • - → n . (3.66)
Eq.(3.65) is valid only in the free-stream. However, we can assume that this definition is still valid at the edge of the boundary layer and it can be used as an indicator of the pressure gradient imposed by the free-stream flow everywhere inside the boundary layer. In order to provide a local approximation of λ θ , the momentum thickness θ is replaced by the normal distance to the wall y. This is doable because the transition model activates at the middle of the boundary layer, where y ensures a good scaling of θ, since at that position y = δ/2 ∼ θ. The final formula for λ θ,L reads as:

λ θ,L = -7.57 • 10 -3 dv dy y 2 ν + 0.0128. (3.67)
0.0128 is added in order to account for the fact that dv dy is never zero in the middle of the boundary layer, neither for zero pressure gradient flow. This is due to the growth of the boundary layer thickness. The scaling is meant to ensure that in the middle of the boundary layer λ θ,L matches the solution λ θ of the Falkner-Skan equations. For numerical reasons, λ θ is bounded as follows: λ θ,L = min(max(λ θ,L , -1), 1).

(3.68)

The pressure is not directly used and the pressure parameter is estimated through the velocity derivative inside the boundary layer, based on 2D considerations. These assumptions are valid for all the streamwise transition processes, at least for Tollmien-Schlichting waves (T-S). A different scaling of λ θ,L is needed in order to account for three-dimensional crossflow instabilities.

Empirical correlation

Once the local quantities are estimated, the correlation is defined as follows:

Re θc (T u L , λ θ,L ) = C T U 1 + C T U 2 e -C T U 3 T u L F P G (λ θ,L ) , (3.69)
where

C T U 1 = 100, C T U 2 = 1000, C T U 3 = 1. (3.70)
The constants of the correlation were calibrated on a series of Falkner-Skan profiles and account for different behaviors:

• C T U 1 defines the minimal value of Re θc , the exponential goes to zero for high value of turbulence intensity;

• C T U 1 + C T U 2
defines the maximal value of Re θc ;

• C T U 13 establishes the exponential rate, i.e. the speed at which Re θc decreases for increasing T u L ;

• F P G (λ θ,L ) accounts for the streamwise pressure gradient.

F P G is empirical and calibrated on a series of Falkner-Skan profiles, as well. It is defined as:

F P G (λ θ,L ) = min(1 + C P G1 λ θ,L , C lim P G1
), for λ θ,L 0, min(1 + C P G2 λ θ,L + C P G3 min[λ θ,L + 0.0681, 0], C lim P G2 ), for λ θ,L < 0.

(3.71) With respect to its definition, the role of each constant in Eq.(3.71) is clear: C P G1 controls the critical momentum thickness Reynolds number for favorable pressure gradient, C P G2 for adverse pressure gradients. C P G3 is activated when separation occurs, at λ θ = -0.0681. The constants are set to:

C P G1 = 14.68, C P G2 = -7.34, C P G3 = 0.0, (3.72)
C lim P G1 = 1.5, C lim P G2 = 3.0. (3.73)
In order to avoid negative values, F P G is bounded as:

F P G = max(F P G , 0). (3.74)
This new definition of Re θc accounts for the significant change of the relative difference between Re V and Re θc as a function of the pressure gradient. This behavior is not considered in γ -Re θ formulation.

Re V to Re θ scaled ratio as a function of the stream-wise shape factor is shown in Fig. (3.1).

The black squares are from Langtry assumption in γ -Re θ , while the red line represents the Falkner-Skan solutions. Resolving the Falkner-Skan equations, we observe that this ratio changes considerably depending on the shape factor, therefore the pressure gradient, by as much as a factor of 1.5. This pressure gradient effect on the scaled ratio can be accounted for either through an additional F onset correction either through a correlation, as it is the case of Eq.(3.69) in γ. Coupling with k -ω SST turbulence model As for γ -Re θ , γ transition model interacts only with the turbulence kinetic energy transport equation of the SST turbulence model, while the equation for ω rests unchanged. The production and destruction terms of k are modified as follows:

∂(ρk) ∂t + ∂(ρu j k) ∂x j = P k + P lim k -D k + ∂ ∂x j µ + σ k µ t ∂k ∂x j , (3.75)
The primary production term P k is defined as:

P k = γP k , (3.76)
while, the destruction term as:

D k = max(γ, 0.1) • D k . (3.77)
The P k in Eq.(3.76) is computed using Kato-Launder modification. The additional production term P lim k is defined as follows:

P lim k = 5C k max(γ -0.2, 0)(1 -γ)F lim on max(3C SEP µ -µ t , 0)SΩ, (3.78)
where

F lim on = min max Re V 2.2 • Re lim θc -1, 0 , 3 , (3.79)
Re lim θc = 1100, (3.80)

C k = 1.0, C SEP = 1.0. (3.81)
In Eq.(3.78), the max(γ -0.2, 0) ensures that the additional production term is activated only when the transition is triggered, i.e. when the intermittency is bigger than 0.2. max(3C SEP µ-µ t , 0) switches the additional term off in fully turbulent regions where 3C SEP µ < µ t . The limit Re lim θc = 1100 ensures the activation of the production term for high Reynolds numbers and/or separating flows, when Re V becomes larger than this limit. This additional source term makes the transition prediction more reliable when it develops under low turbulence intensities, as it is the case of separation-induced transition. The blending function F 1 is modified as in γ -Re θ , see Eq.(3.51) and Eq.(3.52).

k γ modification

In Sandhu et al., [START_REF] Sandhu | A local correlation-based zero-equation transition model[END_REF], a modification of γ model was presented, where it is reduced to zeroequation transition model, the k γ model. The transport equation for the intermittency is incorporated into the transport equation for k through the definition of an additional production term, to predict transition directly. In this version, the intermittency is not directly calculated, but it is estimated using the eddy viscosity ratio R T as:

γ ∼ γ = (1 -e R T ) n , (3.82)
where the constant n is numerically calibrated in order to assure γ = 0 in a laminar region and to match with γ from Menter in the fully turbulent region. This formulation for γ makes the model sensitive to the initial conditions, especially for low free-stream turbulence environment, as discussed in the model's publication. Indeed, a laminar or fully turbulent flow initialization can affect the transition location. k γ performs very similarly to γ model. It can be easier to implement and cheaper in terms of computational cost. It shows faster convergence in case where both models perform similarly, predicting similar transition locations. Efforts are still ongoing for the validation of this model, and its dependence on the inlet conditions has not been clarified yet. For these reasons, it has not been considered further.

Considerations on γ and γ -Re θ models

It is worth it to underline that despite the wide use of the Local Correlation Transition Models, γ -Re θ and γ, in the CFD community, there are still open questions on their modeling approach. First of all, for both models, the relation between the vorticity Reynolds number and transition momentum thickness Reynolds number is only valid for laminar flows, see Eq.(3.19) and Eq.(3.55). The capability of these models to predict the transition process does not lie on transition physics, instead it is related to the design of the transition onset functions. As the transition onset criteria are met, the flow is forced to be turbulent by the function F onset,3 , see Eq.(3.21) and Eq.(3.57).

In the γ -Re θ formulation, the Re V to Re θ ratio defined in Eq. (3.19) is only valid for zero pressure gradient. Indeed, resolving the FSK equations, it was shown that this ratio changes considerably as a function of the shape factor, see Fig. (3.1). This might also affect the transition prediction when laminar separation occurs. For instance, the constant 3.235, used in Eq.(3.44), is set with respect to the value of the ratio Re V 2.193Re θ at H 12 ∼ 3.5. Pressure gradient effects are accounted for in γ formulation through the definition of an empirical correlation calibrated on Falkner-Skan profiles for the critical Reynolds number, Eq.(3.69).

Both models were calibrated on a very large number of flat plate test cases, in the transitional Reynolds number range of Re < 10 7 , the T3 series from the ERCOFTAC database and Schubauer and Klebanoff experiments. This calibration resulted in a large set of models constants. For instance, the constants s 1 in Eq. (3.44) or C SEP in Eq.(3.78) allow to control the length of the separation bubble, and higher values than those imposed by default lead to a shorter reattachment length. Each constant has its proper role and it might be possible that one set of constants cannot be sufficient to cover all possible transition dominated flows.

Another big limitation is the models sensitivity to the initial conditions, namely the free-stream turbulence level T u and the eddy viscosity ratio R T . In γ -Re θ , the momentum thickness θ t , computed from the Re θt transport equation, is sensitive to the free stream turbulence intensity. This dependency is nonphysical and it might cause significant deviation from the experimental results for some flows. For instance, if the boundary layer thickness estimated by the transport equation is thicker, this could cause an early transition. The turbulence viscosity ratio R T controls the decay of turbulence kinetic energy in the far-field, but also enters both the transition models via different functions, as F onset,3 , see Eq.(3.21) and Eq.(3.57), or F turb , see Eq.(3.25) and Eq.(3.60), controlling the transition process. In light of this sensitivity, the initial settings, inlet boundary conditions and initial conditions, are crucial for the transition predictions.

Both transition models were implemented in ISIS-CFD as proposed in their original publications. It is important to keep these considerations in mind for the following discussion of the numerical results.

Chapter 4

ISIS-CFD

This chapter aims to describe ISIS-CFD incompressible flow solver, the algorithms, and methodologies employed for the spatial equations discretization. The solver has been developed by the research group METHRIC, acronym for "Modélisation des Ecoulements Turbulents à Haut Reynolds Incompressibles et Couplages", which is part of the research laboratory LHEEA at the École Centrale de Nantes. The solver is commercialized by the company Cadence Design System, Inc. as part of the software FINE T M /Marine computing suite since 2007. The flow solver is combined with the mesh generator HEXPRESS T M , which generates full hexaedral unstructured meshes. ISIS-CFD solves the unsteady Reynolds averaged Navier-Stokes equations (URANS) and it is based on a fully-unstructured (face based) finite volume discretization. Specific functionalities account for multiphase flows and industrial applications, see Queutey & Visonneau, [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF], and Leroyer & Visonneau, [START_REF] Leroyer | Numerical methods for ranse simulations of a self-propelled fish-like body[END_REF]. Other than URANS turbulence models, Hybrid RANS/LES (HRLES) models based on Detached Eddy Simulation (DES-SST,DDES-SST, IDDES) are also available within the code and further details can be found in Guilmineau et al., [START_REF] Guilmineau | Cross wind effects on a simplified car model by a DES approach[END_REF]. The velocity and pressure coupling rests on the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm with the flux reconstruction based on Rhie &Chow method: at each time step, the velocity field is updated by the momentum balance. The pressure is obtained using the mass conservation equation transformed into a pressure equation. For the integration of the free surface and the simulation of a two-phase fluid medium, the code adopts a volume of fluid (VOF) type interface capture method. For URANS, temporal discretization is done through the Backward Difference of 2 nd order. Different methods of adaptive refinement and sliding meshes are available in the commercialized version. The adaptive grid refinement (AGR) method is based on the creation of optimal meshes through the recursive local division of the existing cells into smaller ones, in regions where more precision is demanded. For instance, this is an ideal strategy for flows around ships, that present several phenomena that are highly localized in space and require locally fine grids to be correctly resolved. The refinement procedure is called repeatedly during the flow computation, setting a given number of time steps between each call. First, the procedure checks whether the current grid satisfies the refinement criterion, and, if not, the grid is adapted by the refinement or coarsening of the earlier refinement. For steady flows, the procedure eventually converges and the grid is no longer changed when the procedure is called. For unsteady flow, the grid changes permanently. The refinement can be either isotropic or directional, i.e. based on a vector field, according to the flow features. Detailed description of AGR strategy can be found in Wackers et al., [START_REF] Wackers | Adaptive grid refinement for hydrodynamic flows[END_REF], [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF] and [START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF].

RANS equations formulation has already been introduced in Chapter 3. Hereafter, we aim to define the numerical framework within ISIS-CFD solver operates. In the following we will detail the spatial schemes used for the equations discretization. We will conclude with the description of the velocity-pressure coupling algorithm.

Numerical Framework

Given a generic scalar quantity Q, a general expression for its transport equation can be written in the following form:

∂ ∂t V ρQ dV + V ρQ ( - → U - - → U d ) • - → n dS = V P Q dV - V D Q dV + S Γ Q ∇(Q) • - → n dS. (4.1)
P Q and D Q represent the production and dissipation terms, respectively.Γ Q characterizes an isotropic or anisotropic diffusion coefficient. Within finite volume approach, the computational domain is discretized with an unstructured grid. Each individual cell volume V is considered as a control volume where the integral formulation of the conservation equations has to be satisfied. All variables are positioned at the center of the control volumes (collocated cell-centered arrangement). No hypothesis is made on the control volumes topology. In particular, a volume is made of an arbitrary number of constitutive faces noted f , composed by an arbitrary number of nodes. After discretizing the Navier-Stokes equations, in order to determine the evolution of variables, such as pressure and velocity, the evaluation of the volume and surface integrals presented above is required. To this purpose, volume centers, face centers and face vectors of each of the different volumes have to be determined. Given the spatial discretization of the finite volume method, it turns out that these integrals can be calculated by adding the contributions of each of the faces that constitute the volume V . The integral volume on V of the quantity Q can be approximated by:

V Q dV = Q C V dV = Q C V, (4.2)
where Q C is the value of Q at the center of the cell. The surface integral can be determined using the approximation:

S Q - → n dS = f S f Q f - → n f dS f = f -→ S f Q f , (4.3)
where -→ S f = -→ n S f is the surface-oriented vector. Q f is the value of Q at the center of the face f and it can be reconstructed from the values at the center of the cells on each side of the considered face f . We are going to refer to these values as Q L and Q R , where the subscript L stands for left, R for right.

Evaluating the volume and surface integrals using Eq.(4.2) and Eq.(4.3), the semi-discrete formulation of Eq.(4.1) yields:

∂ ∂τ (ρV Q)| C + ∂ ∂t (ρV Q)| C + f (F C f -F D f ) = (S V Q ) + f (S f Q ), (4.4)
with

F C f = ρ( - → U - - → U d ) f • -→ S f = ṁf Q f , F D f = (Γ Q ) f ( --→ ∇Q f • - → i k )(S f ) k . (4.5)
ṁf is the mass flux across the face f and -→ i k the identity vector with the k-th component different from zero. The terms F C f and F D f are respectively the convection and diffusion fluxes across the face. S V Q and S f Q are volume and surface source terms. τ is a local fictitious time variable, introduced to enhance the diagonal dominance of the system of linearized equations. Temporal derivatives in Eq.(4.5) are evaluated using upwind second-order scheme, also known as implicit Euler scheme or BDF2. The time derivative is approximated as:

∂Q ∂t ∼ e c Q c + e p Q p + e q Q q , ( 4.6)
where the subscript c stands for the current time step, p for the previous and q the time step anterior to p. The coefficients {e c , e p , e q } can be calculated as

e c = 2∆t c + ∆t p ∆t c (∆t p + ∆t c ) , e p = ∆t c + ∆t p ∆t c ∆t p , e q = ∆t c ∆t p (∆t c + ∆t p ) . (4.7)
The time step, ∆t, can be non-constant and not necessarily uniform.

The fictitious time derivative is evaluated by:

∂Q ∂τ = (Q c -Q c0 ) ∆τ , (4.8)
where Q c0 is the previous estimation of the quantity Q c , within the non-linear loop.

Spatial discretization schemes

Centered face reconstruction

With the exception of convection terms and mass fluxes, the variables at the center of a face are reconstructed using centered approximations, based on linear interpolations between cell and face centers. For the discretization to be effectively of order two, it is necessary that the reconstruction at the faces has the same precision. We assume the presence of a known discontinuity for a variable Q and its gradient across the interface Γ of center f . This interface divides the computational domain into the disjoint subdomains Ω + and Ω -. Notations are presented in Fig. (4.2). We will suppose that the discrete location of discontinuity coincides with the mesh faces. Then, Q discontinuity across the interface can be approximated as:

[Q] Q + f -Q - f . (4.9)
Jump conditions along Γ are defined by the following relations:

[Q] = a, (4.10) [c --→ ∇Q] = - → b , (4.11)
where a, -→ b , and c are known functions. c is assumed to be continuous on each subdomain, but it may be discontinuous across the interface Γ. Knowing a-priori the gradient of Q on both sides of the face f from the cell-centered approximation:

-

-→ ∇Q = 1 V V --→ ∇QdV f Q f -→ S f , (4.12)
it is possible to reconstruct the cell-centered quantities Q L and Q R , by retaining a second order Taylor expansions:

Q L Q f -- --→ C L f f • - → ∇Q f - Q f --h -( - → ∇Q f -• - → n - - → ∇Q f -• - → e -), (4.13)
and

Q R Q f + - --→ f C R f • - → ∇Q f + Q f + -h -( - → ∇Q f + • - → n - - → ∇Q f + • - → e + ), (4.14)
with

h -= --→ C L f f • - → n , h + = --→ f C R f • - → n , - → e -= - → n -β + - --→ C L f f , - → e + = - → n -β + --→ f C R f . (4.15)
The subscripts ± refers to Ω ± domains. The coefficients β ± are chosen to fulfill conditions

- → e + • - → n = 0 and - → e -• - → n = 0.
Deduction of centered differencing reconstructions for Q f + and Q f -and the normal gradients to the face

- → ∇Q f + • - → n and - → ∇Q f -• - → n ,

Upwinded face reconstruction

To avoid the appearance of nonphysical oscillations in the solution and to reinforce the stability of the numerical scheme, off-center reconstructions are used to treat the convection terms. The order of precision of these reconstructions is between 1 and 2, and depends on both the mesh and the physical nature of the problem treated. The 1st order Upwind Differencing Scheme (UDS) accounts for the direction of the flow. It retains as quantity at the center of the face the variable located upstream with respect to the speed flow through the face. The scheme is unconditionally stable and very diffusive, therefore it requires very fine meshes to obtain precise solutions. In ISIS-CFD is a combination of upwind and centered schemes is available and it is known as the hybrid scheme. It results from a local blending between linearly interpolated values of neighboring cells. The blending factor is not uniform, but it depends on the number of Peclet evaluated at the face as:

P e f = ṁf ----→ C L f C R f 2SΓ Q . (4.16)
The accuracy of this reconstruction is not necessarily uniform but locally adapted. Using UDS, the relative orientation between the face and the velocity is taken into account.

A blending factor d is computed using an exponential scheme, which ensures a smooth transition. Q f is approximated as:

Q f = d L Q L + d R Q R + d L --→ C L f f • - → ∇Q | L + d R --→ f C R f • - → ∇Q | R, (4.17)
with

d = exp(P e f ) (1 + exp(P e f ))
and

d L = 1 -d R . (4.18)
This method however exhibits a lack of stability and does not guarantee the boundedness of the solution, which is necessary, for example, for multi-phase flows in which the concentration must remain between 0 and 1. To overcome these limitations new schemes, based on the Normalized Variable Diagram analysis (NVD) introduced by Leonard, [START_REF] Leonard | Simple high-accuracy resolution program for convective modelling of discontinuities[END_REF], have been proposed. One of these schemes is the Gamma Differencing Scheme (GDS), from Jasak, [START_REF] Jasak | High resolution NVD differencing scheme for arbitrarily unstructured meshes[END_REF], which proposes a transition between a second order Central Differential Scheme (CDS) and the Upwind Differencing Scheme (UDS). It is supposed that the value of Q is available in three points, U, C and D, located along the convection direction. To reconstruct the quantity at the face f located between C and D, the normalized variables system is defined as: The GDS scheme proposes the reconstruction based on these three points, where Qf = f ( QC ).

Q = Q -Q U Q D -Q U = f (Q) (4.19) U C D UD Flow direction CD
To avoid any nonphysical oscillation, it is necessary to bound Q C between min{Q U , Q D } and max{Q U , Q D }. In order to preserve the boundedness at the face, the following constraints should be verified:

• QC < 0 , Qf = QC ,
• 0 ≤ QC ≤ 1, Qf is bounded by Qf ≥ QC and by unity,

• for QC > 1, Qf = QC .

These conditions are met in the gray area of the plots shown in ). The implementation of AVLSMART for non-uniform grids is detailed in Table (4.1): the base scheme is the third-order QUICK scheme, by Leonard, [START_REF] Leonard | A stable and accurate convective modelling procedure based on quadratic upstream interpolation[END_REF]. The interpolation factor f x , that accounts for non-uniform grids, is defined, referring to notations in Fig. (4.2), as the ratio:

f x = || -→ f D|| || --→ CD|| . (4.20) QC Qf Q f Scheme ] -∞, 0] QC Q C UDS ]0, 1/4[ 9 2 QC C D Q D + (1 -C D )Q C , C D = 5 2 (1 -f x ) QC 1-QC UDS to QUICK ]1/4, 3/4[ 3 8 + 3 4 QC C D Q D + (1 -C D )Q C , C D = 1 4 (1 -f x ) 3-2 QC 1-QC QUICK [3/4, 1[ 3 4 + 1 4 QC C D Q D + (1 -C D )Q C , C D = 3 2 (1 -f x ) QUICK to UDS [1, ∞[ QC Q C UDS Table 4
.1: AVLSMART scheme implementation for non-uniform grids.

Velocity-Pressure coupling algorithm

In order to resolve the pressure, a pressure-equation is obtained by transforming the continuity equation. The approach adopted reconstructs the velocity flows at the interfaces using a pseudo-physical approach, which in turn uses the semi-discretized formulation of momentum balance. The method is based on a Semi-Implicit Method for Pressure Linked Equations (SIM-PLE) type algorithm, see Caretto et al. [START_REF] Caretto | Two calculation procedures for steady, three-dimensional flows with recirculation[END_REF], Patankar et al. [START_REF] Patankar | A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[END_REF], and Issa, [START_REF] Issa | Solution of the implicitly discretised fluid flow equations by operator-splitting[END_REF], formulated in the form of a pressure equation with a flow reconstruction like the one described in the work of Rhie and Chow, [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF]. The peculiarities of this formulation are the management of pressure gradient discontinuities in presence of a density discontinuities, as well as the treatment of pseudo-unsteady terms. The latter are interpolated in order to have a solution independent of the time step ∆t and of the local fictitious time step ∆τ , when an overall steady solution is expected. The presented method uses hydrostatic discontinuity hypothesis for the pressure and continuity hypothesis for the velocity field, i.e. viscous effects and surface tension jumps across the interface Γ are neglected. These hypothesis are taken into account by writing the two following conditions at the faces:

[p] = 0, -→ ∇p • - → n ρ = 0. (4.21)
Isolated the pressure gradient along with the gravity acceleration in Eq.(4.4) and considering that the time derivatives are treated as in Eq.(4.6) and Eq.(4.8), the semi-discretized momentum equation takes the form:

(e c + 1/∆τ C )(V ρ - → U ) c C + eV ρ - → U ) p C + (eV ρ - → U ) q C + a C - → U c C + nb a nb - → U c nb + - → S C + (V -→ ∇p) c C = (ρV ) c C - → g + (ρV ) c C - → U c0 C /∆τ C . (4.22)
a C and a nb are the diagonal and extra-diagonal matrix coefficients from the implicit part of convection and diffusion terms. V C represents the volume of the cell C, while -→ U nb stands for the velocity variables of its neighbors.

-→ S is a source term containing all explicit contributions, and external force fields other than the gravity and pressure. The velocity variables at the cell centers can be expressed in a compact way as:

- → U c C = -C p C ( - → Û C + (∇p/ρ) c C -- → g ) + C p C - → U c0 C /∆τ C -C p C [(eV ρ - → U ) p C + (eV ρ - → U ) q C ]/(ρV ) c C , (4.23)
with

C p C = (e c + 1/∆τ C + a C /(ρV ) c C ) -1 . (4.24)
The vector -→ Û C called pseudo-velocity is homogeneous to gravity acceleration and contains part of convection, diffusion and the source terms. It is defined as:

- → Û C = nb a nb - → U c nb + - → S C /(ρV ) c C . (4.25)
According to the method presented by Rhie and Chow, it is assumed that the velocities at the interfaces, referred to as -→ U c f can be expressed in a similar way as in Eq.(4.23). The quantities at the face center f , that substitute the variables at the cell center in Eq.(4.23) and (4.24), are interpolated by the available cell quantities (C L f , C R f ), except from the pressure gradient. The latter is discretized in the center of the face to avoid the so-called odd-even decoupling problem. The flux F ( -→ U ) through a face f reads as:

F ( - → U ) = -C p f (F ( - → Û ) + F (∇p/ρ) c f -- → g ) + C p f F ( - → U c0 f )/∆τ f -C p f (eV ρ) p F ( - → Û p ) + (eV ρ - → U ) q f F ( - → Û q ) /(ρV ) c f . (4.26)
When the velocity reconstruction at the cell face is substituted in the continuity equation, the pressure equation can be written as:

-D(C p -→ ∇p/ρ) =D(C p - → Û ) -D(C p - → U c0 /∆τ ) -D(C p - → g ) + D C p (eV ρ) p - → Û p + (eV ρ) q - → Û q /(ρV ) c f . (4.27)
Provided that the pressure equation is satisfied, and guaranteed that the volumetric flux in Eq.(4.26) is conservative, the discretized equation over the control volume reads as:

- f C p f S f p R -p L hρ = f C p f F ( - → Û ) + S f E p hρ -F ( - → g ) + F i . (4.28)
h is the distance to the face, as in Eq.(4.15) and E p represents the explicit part of the pressure gradient flux. Unsteady and pseudo-unsteady fluxes are gathered into the flux F i . ρ is the reconstructed density taking into account the hydrostatic discontinuity. The matrix assembled from all the control volume is sparse, symmetric and positive definite.

Chapter 5

2D Configurations: Results

Introduction

This chapter is devoted to the analysis of the performance of both LCTM models γ and γ -Re θ for two-dimensional configurations. We consider the ERCOFTAC flat plate, NACA 0015 and Eppler 387 profiles. 2D test cases were used both to validate the models implementation in ISIS-CFD solver, and as contributions to the "AVT-313 Incompressible Laminar-to-Turbulent Flow Transition Study Comparison Workshop", [START_REF] Eça | AVT-313 incompressible laminar-to-turbulent flow transition study comparison workshop[END_REF], an ongoing exercise within the NATO/AVT-313 collaborative group. The aim of the workshop is to perform a verification and validation exercise for two-dimensional transition mechanisms, in order to assess the reliability of numerical predictions. Solution verification substantiates the correctness of the input data, the numerical accuracy of the solution obtained, and the correctness of the output data for a particular simulation, as stated in Oberkampf and Roy, [START_REF] Oberkampf | Verification and validation in scientific computing[END_REF]. On the other hand, models validation quantifies the accuracy of the numerical predictions by comparing numerical results with experiments. This chapter is dedicated to the discussion of some of the results of this verification and validation exercise. For each test case, we will give details on the experimental results used for the comparison and we will discuss the numerical predictions with their numerical uncertainties. The procedure to estimate the numerical uncertainties is detailed in Sec.(5.2.1). Test cases conditions are given in Table (5.1), that includes the angle of attack, if any, the transition model used for the computations, the turbulence intensity T u in , at the inlet and the type of transition mechanism we are validating the model for. Details on the procedure for setting the inlet conditions are given in Sec. (5.3). For the flat plate T3A test case we will also introduce the multivariate metric for validation tool proposed by the ASME V& V 20 Standard Committee [START_REF]Supplement to ASME V& V 20-2009-Multivariate metric for Validation[END_REF] as a supplement of the V& V 20-2009 Standard Committee [START_REF]Standard for Verification and Validation in Compuatational Fluid Dynamics and Heat Transfer[END_REF]. These multivariate results were presented at the V& V ASME Validation and Verification symposium in Las Vegas (2019).

The quantities of interest used to asses the performance of the transition models are:

• The skin friction coefficient C f defined as:

C f = τ w 1 2 ρU 2 ∞ , (5.1)
where τ w is the shear stress at the wall and it is defined as:

τ w = µ ∂u ∂y y=0 , (5.2)
where u is the flow velocity along the boundary layer edge. As mentioned in Chapter 2, the skin friction coefficient characterizes the transition region. The minimum skin friction C f represents the start of the transition process. The transition region is identified by its sharp increase. When separation-induced transition occurs, the C f assumes negative values in the separated region. The reattachment position corresponds to the point where the C f changes sign and becomes positive again. This quantity well characterizes the transition process, however, it is very hard to measure it experimentally.

• The second quantity is the pressure coefficient C p defined as

C p = (p -p ∞ ) 1 2 ρU 2 ∞ . (5.3)
C p distribution is particularly significant when separation-induced transition occurs, because it allows to identify the laminar separation bubble, as shown in Fig.

(2.9). In the numerical results, the maximum pressure at the inlet is taken as reference value p ∞ .

Whenever possible, we are going to refer also to the drag C D , its friction and pressure components, and the lift C l coefficients. They are defined as:

C D = D 1 2 ρU 2 ∞ L , (5.4)
C l = L 1 2 ρU 2 ∞ L , (5.5) (5.6)
where D and L are the drag and lift forces, respectively. L is the reference length, which is substituted by the reference surface S for 3D configurations. Mesh Considerations All the meshes for the 2D test cases were provided by Rui Lopes and Luis Eça, IST Lisbon. They are generated using in-house generation tools, whose description can be found in Eça, [START_REF] Eça | Practical grid generation tools with applications to ship hydrodynamics[END_REF]. The meshes are multi-block structured grids, all characterized by an O-topology at the leading edge and a y + max lower than one. Not having taken care personally of the mesh generation, we did not perform an analysis of the influence of the first near wall cell on the transition predictions. Nevertheless, this analysis was performed by Rui Lopes and discussed in his PhD dissertation, [START_REF] Lopes | Simulation of Transition from Laminar to Turbulent Regime in Practical Applications of Incompressible Flow[END_REF]. Hereafter, we report the main remarks of this study, conducted for the flat plate on nine different sets of five grids.

Geometry

Each set is generated such that:

• the coarsest grid has an averaged y + of one for set 9;

• the coarsest grid has y + max ∼ 1 for set 5;

• the coarsest grid has y + max ∼ 1 for set 1.

For each grid set, the level of grid refinement and total cell count is kept constant. The number of cells on the surface of the coarsest grid is 512 up to 1024 on the finest. Only the size of the cells in the wall-normal direction is varied, as indicated above.

This study showed a clear influence of the y + max value on the results, independently of the used transition model. It affects the C f distribution on the body, and, consequently, the friction drag. For γ -Re θ coupled with k -ω SST (2003), the y + max effect is clearly visible on the transition and turbulent region. Notably, the C f value increases in the turbulent region for decreasing y + values. On the contrary, the transition point position is independent of the wall-normal grid refinement. The same conclusions are valid for γ model. Nevertheless, the latter model shows a stronger dependence on the grid refinement in the streamwise direction compared to γ -Re θ . In conclusion, y + max influences the solution the most in the turbulent region, indicating that the mesh requirements for transition models are very similar to those for turbulence models. y + max ∼ 1 is the upper bound limit in order to obtain accurate numerical results. Nevertheless, for the majority of the test cases discussed hereafter, y + max < 0.5 has been used. The use of such dense meshes are related to the will of reducing to minimum values the discretization error, which is the main component of the numerical uncertainty. The definition of numerical uncertainties is given in the following section, devoted to the basic concepts of the validation and verification procedure.

Validation and Verification

Given the quantity φ, we define the error E φ as:

E φ = φ obtained -φ true .
(5.7)

φ true is the true value of the quantity φ and φ obtained is the obtained value of the quantity φ.

By "obtained value" we mean that it can be derived from any source, e.g. numerical solution, computational simulation, or experimental measurements. The error E φ can be decomposed in its different contributors as:

E φ = E model + E num + E input . (5.8)
E model is the modeling error and it is given by the difference of the true value and the exact solution of the mathematical equations. E input is input/parameter error relates to possible differences between numerical and experimental setting, e.g. in terms of boundary conditions, computational domain and geometry features. E num is the numerical error, whose contributions are:

• Round-off error, consequence of the finite precision of computers and algorithms.

• Iterative error, resulting from the iterative resolutions of the non-linear system of equations. This error can be assessed using either the L ∞ or L 2 norms.

Given the solution φ and being N cell the total number of cells we have

L 2 (φ) = N cell i res(φ i ) 2 N cell , (5.9)
and

L ∞ (φ) = max i=1,...,N cell |res(φ i )|.
(5.10)

• Discretization error, due to the spatial and temporal discretization. It depends on grid refinement and temporal steps.

• Statistical error, that occurs in scientific computing approaches that are stochastic in nature and requires time or ensemble averaging to determine mean system response. For steady and deterministic simulation, this kind of error cannot occur.

Once defined the error E φ , we can distinguish between validation and verification. On one hand, the output of a validation exercise are modeling error and uncertainties. It requires an attentive comparison with experimental measurements. On the other hand, solution verification only needs the computational simulation. It assesses the accuracy of the numerical solution with respect to the exact solution of the mathematical model, which is generally unknown. For this reason, we are going to refer to uncertainties rather than errors, because in lack of the knowledge of the true value, it is more appropriate to characterize the accuracy of the results as uncertainty.

We now introduce the definition of the numerical uncertainties, whose main contributor is considered the discretization error. Round off and iterative errors are assumed to be negligible, if they are two times smaller in respect of the discretization one. All the calculations presented hereafter are performed in double precision, which allows to neglect the round-off error. Concerning the iterative convergence, the simulations are run until the L 2 norm of the normalized residuals gain 4 orders minimum. This is a fairly acceptable criterion for turbulent/transitional flows. With respect to this gain, the iterative error is not always negligible with respect to the discretization one. Nevertheless, it is not considered in our simulations.

Discretization Error and Numerical Uncertainty

The behavior of the error of the discrete solution relates to the convergence order p of the discretization method. Relying on Richardson extrapolation (RE), the exact solution φ 0 can be expanded in respect of the discrete solution φ i as

φ 0 = φ i (x) + h p i α + o(x, h q i ).
(5.11) o(x, h q i ) includes all the terms of order q > p. The index i relates to the different grids on which the solution φ i is computed and α is the error constant. The parameter h i is the typical cell size and in its most general formulation it is defined as:

h i = N cell i=1 ∆V i N cell , (5.12)
with ∆V i being the volume of the i-th cell.

Based on the expansion in Eq.(5.11), the discretization error δ RE can be written as

δ RE = φ i -φ 0 = αh p i . (5.13)
This expansion can be manipulated to obtain an expression for the order of convergence p. The expected p for a finite volume method is p = 2. Theoretically, p can be assessed using the results φ i as computed on three topologically similar grids. Often, such a low number of grids is not enough, because of data scatter caused by the noise of the solutions. In addition, Eq.(5.11) is only valid for meshes in the asymptotic range. The grids need to be fine enough to meet this assumption. For these reasons, it is more appropriate to approximate the order of convergence using a minimum number of five grids and a least squares fit. p is obtained by minimizing the function:

S RE (φ 0 , α, p) = Ng i=1 [φ i -(φ 0 + αh p i )] 2 , (5.14)
where N g is the number of grids used for the approximation. Once estimated φ 0 , α and the order of convergence p, the reliability of the fit in Eq.(5.14) is evaluated with respect to its standard deviation σ RE , defined as:

σ RE = Ng i=1 (φ i -(φ 0 + αh p i )) 2 N g -3 , (5.15)
and the data range parameter ∆(φ i ), defined as:

∆(φ i ) = (φ i ) max -(φ i ) min N g -1 . (5.16)
The fit can be considered reliable if:

1 ≤ p < 2.1, (5.17)
σ RE < ∆(φ i ). (5.18)
Because of the sensitivity of the convergence order p to any noise in the data used for the fit, p might be negative or might assume very small or high values, outside the range defined in Eq.(5.17). The estimation of p using Eq.(5.14) might not be sufficient and additional functions can be used for the least square fit, in order to meet the conditions of Eq.(5.17) and (5.18). Details on the estimation procedure of p are given in Appendix A, following Eça et al., [START_REF] Eça | Evaluation of numerical error estimation based on grid refinement studies with the method of the manufactured solutions[END_REF], [START_REF] Eça | A verification and validation exercise for the flow over a backward facing step[END_REF], [START_REF] Eça | A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies[END_REF]. Once the discretization error δ RE has been estimated, we can define the numerical uncertainty u num . This quantity permits to define the interval which contains the exact solution φ 0 with 95% of confidence, such that:

φ i -u num (φ i ) ≤ φ 0 ≤ φ i + u num (φ i ). (5.19)
The numerical uncertainty is calculated in respect of the solution on the finest grid following a modified Grid Convergence Index (GCI) procedure, which takes into account the standard deviation of the fit, other than the discretization error. If the conditions of Eq.(5.17) and Eq.(5.18) are met, the estimated numerical uncertainty reads as:

u num (φ i ) = F s δ RE + σ RE .
(5.20)

The error estimator δ RE is converted into a numerical uncertainty through the safety factor F s , following the approach of Roache, [START_REF] Roache | Verification and validation in computational science and engineering[END_REF]. Other expression of u num (φ i ) for the different fits are given in Appendix A.

Numerical uncertainty is an important tool that allows to define the credibility of the solution of the mathematical model. For all the test cases, the flow solution is computed on five topologically similar grids. The solutions on the five meshes are used to compute the numerical uncertainties of all the quantities hereafter presented. Before discussing the results, let us discuss the inlet boundary conditions, which play a crucial role in transition simulations.

Inlet Boundary Conditions

Transition is an initial value problem: the transition process depends on the initial disturbances, on their frequency and their amplitude. These initial conditions are translated within the RANS framework in the values of turbulence intensity T u, and eddy viscosity ν t , imposed at the inlet. T u quantifies the velocity fluctuations, while ν t is linked to the turbulence reference length L turb by the relation:

ν t = T uU ref L turb = k/ω. (5.21)
They depend on k, ω and their free-decay ahead of the body. Indeed, in the free-stream, the destruction terms of the transport equations are active. The turbulence quantities undergo a decay that can have a strong impact on the numerical solution within the boundary layer, as discussed by Spalart & Rumsey, [START_REF] Spalart | Effective inflow conditions for turbulence models in aerodynamic calculations[END_REF]. In the following, we present the free-decay equations and discuss the possible strategies to set to the inlet conditions.

Turbulence Quantities Decay Free decay rates of turbulence quantities can be studied from the solution of the k -ω SST (2003) equations in the approaching flow field. The transport equations (3.8) and (3.9) for an incompressible and steady flow reduce to:

U j ∂k ∂x j = ν t S 2 -β * ωk + ∂ ∂x i (ν + σ k ν t ) ∂k ∂x i , (5.22)
U j ∂ω ∂x j = γΩ 2 -βω 2 + ∂ ∂x i (ν + σ ω ν t ) ∂ω ∂x i + CD kω , (5.23)
where the constants are β = 0.0828 and β * = 0.09. For an uniform flow and aligned with x, neglecting the diffusion and cross-diffusion term, the equations simplify to:

dk * dx * = -β * k * ω * , (5.24
)

dω * dx * = -β(ω * ) 2 .
(5.25)

The superscript * indicates the dimensionless variables,

k * = k/U 2 ∞ , ω * = ωL/U ∞
, where U = U ∞ , and x * = x/L. The analytical solutions of Eq.(5.24) and (5.25) are given by:

k * = k * in (1 + β(x * -x * in )ω * in ) -β * β , (5.26
)

ω * = ω * in (1 + β(x * -x * in )ω * in ) -1 .
(5.27)

The solution for ν * t reads as:

ν * t ν = ν * t in ν 1 + β((x * -x * in )ω * in ) β * β -1
.

(5.28)

The subscript in indicates the quantity value at the inlet of the computational domain, positioned at x in . The decay of k * and ω * , as well as ν * t /ν is exponential.

INLET BOUNDARY CONDITIONS

ν * t /ν is the quantity which undergoes the decay at the slowest rate. Rewriting the specific turbulence dissipation rate at the inlet as:

ω * in = k * in ν ν * t in Re, (5.29)
and substituting Eq.(5.29) in Eq.(5.26), we observe that the decay depends on the Reynolds number and can be controlled through the eddy viscosity ratio. Increasing the eddy viscosity ratio at the inlet, the decay of k * and ω * can be contained.

In order to supervise the decay of turbulence quantities ahead of the body, high values of R T are often imposed at the inlet. As discussed by Spalart & Rumsey, these high values of eddy viscosity ratio might pollute the flow field in non-turbulent region. In their publication, they propose as indicative value for ν t /ν ∼ 2 × 10 -7 × Re, for Re in the range 10 6 -10 7 . Nevertheless, these values cannot be normally used for transitional flows, because it would cause such a rapid decay ahead of the body to make the turbulence kinetic energy in the vicinity of the body insignificant for the transition prediction. Keeping in mind that R T values have to set with care, different strategies can be adopted to impose the right boundary inlet conditions for transition models:

1. Setting to specific values of eddy viscosity ratio at the inlet, in order to obtain the required value of turbulence intensity close to body. From Eq.(5.26), (5.28), and (5.29), the value at the inlet in can be computed theoretically, once the targeted values of k and νt ν next to the body are chosen. Normally, k is recovered from the experimental set-up, if any is at disposal.

2.

Defining ambient values k amb and ω amb , as proposed in Spalart & Rumsey, [START_REF] Spalart | Effective inflow conditions for turbulence models in aerodynamic calculations[END_REF]. These values are recovered from the turbulence intensity and eddy viscosity ratio in the freestream. A production term based on these ambient values is added to k and ω equations, in order to overcome the respective destruction terms. These terms are zeroed in the fully turbulent region, in order not to affect the turbulent solution. However, they exist in the laminar boundary layer. For this reason, they cannot be set to significant values, otherwise they could corrupt the solution of the transition model. In addition, this approach requires a re-definition of the no-slip boundary conditions, in order to recover the ones from the original turbulence model.

3.

Freezing the turbulence equations until a certain position, x A , ahead of the body. The values imposed at the inlet are convected until the position x A , where, the dissipation term is reactivated until the leading edge. From x A on, the turbulence quantities undergo their free-decay, as described in Eq.(5.26) and Eq.(5.28), but along a considerably reduced region. This strategy, which is the most commonly used, is exploitable only when the incoming flow is aligned with the axis x.

The 1st and 3rd approaches are the ones generally employed in transition simulations. In our case the 1st method has been opted for and it is used in all our computations. It is worthwhile to mention that both strategies have a common drawback: they are not practical for real applications. Indeed, when dealing with bodies with appendages, different and specific values of turbulence intensities are expected in the vicinity of the different surfaces. The choice of T u in and

ν * t in
ν is crucial and strongly affects transition predictions. The dependence on the inlet turbulence intensity is desirable and respects the physics of an initial value problem. For instance, transition due to T-S waves is fed by the free-stream turbulence.

Different T u can accelerate or delay the transition process. On the other hand, the dependence on the eddy viscosity ratio values is not physical and it is related to the presence of R T within the production and destruction term of the intermittency in both model formulation. In γ -Re θ , this dependence is even stronger, because of the calculation of the empirical momentum thickness Reynolds number, as mentioned in Chapter 3. It is worthwhile to underline that, whatever strategy is chosen, a common issue is that experiments reports only provide the turbulence intensity T u value in one point of the domain. If any free-stream condition is reported. This is not a sufficient information, because turbulence needs at least two variables to be fully described. This uncertainty might cause some inconsistency within a validation exercise. In fact, a rigorous validation exercise cannot be performed if a description of the turbulence quantities experimental decay is not provided.

Once given the setting of the inlet boundary conditions and defined the numerical uncertainty, we can analyze the numerical predictions for two-dimensional flows. The first configuration we are going to focus on is the flat plate, for bypass and natural transition. Then, we will consider the flow around the two airfoils, Eppler387 and NACA0015, at low Reynolds number. We will discuss computations using both γ and γ -Re θ transition models, coupled to k -ω SST (2003) with Kato-Launder production limiter.

Flat Plate

The flat plate was the first test case discussed within the AVT-313 collaborative group. We present hereafter the simulations performed for the Reynolds number Re = 10 7 , based on the length of the flat plate, L, and for different inlet boundary conditions. For the comparison, we are going to refer to the experimental data available at the ERCOFTAC database, [START_REF] Coupland | Ercoftac classic database[END_REF], specifically, the T3A and T3A-test cases. The inlet boundary conditions are chosen in order to meet the experimental turbulence decay ahead of the body. We remind that the ERCOFTAC database was used by Langtry and Menter for the calibration of both γ and γ -Re θ transition models. A comparison with the experimental results is therefore meaningless. Nevertheless, this is a valid test case to start investigating modeling features and identifying some differences between the two transition models γ and γ -Re θ approaches. To these purposes, for one of the set of boundary condition we have used the multivariate metric validation tool to estimate the modeling error. The procedure is described hereafter.

Multivariate Metric

The aim is to determine the interval that contains the modeling error δ model with 95% of confidence. Given the validation comparison error E = S -D, where S are numerical solutions and D experimental measurements, δ model is contained within the interval:

E -u val δ model E + u val .
(5.30)

The quantity u val contains all the possible source of uncertainties. If the error sources are independent, it can be defined as:

u val = u 2 input + u 2 d + u 2 num .
(5.31)

u num and u d take into account the contributions of the numerical and experimental uncertainties, respectively. u input is the uncertainty related to the parameters/input. The input uncertainty is defined as1 :

u 2 input = m i=1 ∂S ∂X i u X i 2 .
(5.32)

m is the cardinality of the input/parameter space and u X i is the standard uncertainty of the input parameter X i , that should be provided by prior experimental tests. ∂S ∂X i is defined as the sensitivity coefficient and it measures the sensitivity of the simulation results S to changes in the input. E can be considered as an estimate of δ model and u val is the standard uncertainty of that estimate. Analyzing E with respect to u val , it can be stated that:

• If |E| u val
, then E ≈ δ model and the main source of error is the modeling one, that has to be reduced to improve the numerical results;

• if |E| < u val , the modeling error is within the noise level imposed by the numerical, input and experimental uncertainties. If E val is small enough, it can be asserted that the model has been validated with u val precision. If it is not the case, this result suggests that the model and/or the experiments should be improved.

We can now define the multivariate metric tool E mv . If we assume that all the estimates of E and u val are distributed as a gaussian, without any lack of generality, given a set of n validation points, a reference value E ref for the error can be defined as:

E ref = n + √ 2n.
(5.33)

Built U val the n × n matrix, whose entrances are given by U num + U d + U input , the validation metrics E mv is defined as follows:

E mv = E T U -1 val E (5.34)
If E mv /E ref is larger than one, there is an indication that the model is not able to reproduce the experimental data within the range of the validation uncertainty at each set point. The point wise information is still essential to assess the level of the validation uncertainty. For instance, it could happen that E mv /E ref is smaller than one for a single point, because u val shows a big value at that specific point. This multivariate metric tool was used for the T3A test case. In order to evaluate the input uncertainty, two additional simulations varying the inlet eddy viscosity ratio R T were run for both transition models. The eddy viscosity is one of the parameters that controls the decay of the turbulence variables ahead of a body. Unfortunately, experimental uncertainties are not reported in the ERCOFTAC database. Uncertainties estimated for this exercise are based on the reported variability between the values of C f measured with different techniques, which is about 2%.

Computational Domain

The computational domain is a rectangle of length 1.5L and height 0.25L. The inlet and outlet are placed 0.25L away from the leading edge of the plate, at x = 0, and from the trailing edge of the plate, at x = L, respectively. The plate is modeled as a zero-thickness surface and the effects of the leading edge curvature on the transition process are not taken into account. No-slip boundary conditions are imposed at the surface of the plate. The surfaces upstream and downstream are treated as slip surfaces, i.e. zero normal derivatives for all variables, and zero normal velocity at the boundary. Within the AVT-313, two different sets of five grids with different topology at the leading edge were initially proposed: the H-Grid and O-Grid topology sets. The H-Grid set consists of Cartesian grids made of a single block, with cells clusters at the leading and trailing edge of the plate in the flow direction. This set was soon discarded, because of the highly stretched and high aspect ratio cells downstream of the plate that resulted in a very slow iterative convergence. Instead, we opted for the O-Grid set, that, contrary to H topology, leads to a much faster convergence, in particular on such fine grids as the ones used for this exercise. Mesh details are reported in the Table (5.2), where the y + max refers to the computations using γ. Similar values of y + max are obtained for γ -Re θ . The grid refinement ratio r i in Table (5.2) is defined as:

r i = h i h 1 = N 1 N i 1/2 , (5.35)
where N i is the number of cells on the surface of the i-th grid. The entire mesh domain and a close up of the leading edge are shown in Fig. (5.1).

N cells N cells-plate r i = h i /h

T3A

Let us discuss the results for the T3A test case. The inlet conditions for this test case are summarized in Table (5.3). We report the free-stream turbulence intensity T u in and the turbulence intensity in the vicinity of the leading edge T u, as obtained by varying the inlet eddy viscosity ratio. ), the input uncertainty is set to zero (u input = 0). This is the so-called strong formulation, because the input parameters are considered as part of the model, Eça et al., [START_REF] Eça | On the interpretation and scope of the v&v 20 standard for verification and validation in computational fluid dynamics and heat transfer[END_REF]. The numerical uncertainty is plotted for each point of the validation set. As expected, given the refinement level in the wall-normal direction of the grids used for the exercise, the numerical uncertainty is everywhere negligible, except in the transition region. These higher uncertainties are related to the steep gradients of C f . For both γ and γ -Re θ , the C f solution in these points, as computed on the five grids, does not converge monotonically. In Lower is the eddy viscosity ratio at the inlet, lower values of C f are predicted in the transition region. This behavior is related to the fact the transition process is kept active by the function

F onset,3 = max 1 -R T a 3
, which is tuned by the eddy viscosity ratio. Because of this latter function, the intermittency production term and, consequently, the values of turbulence kinetic energy are affected.

(a) E(C f ) with u val = (u 2 num + u 2 d ). (b) E(C f ) with u val = (u 2 num + u 2 input + u 2 d ).

T3A-

The second set of boundary conditions is for of the so-called T3A-test case. Compared to the T3A test case, the free-stream turbulence is considerably reduced. The type of transition mechanism is at the limit between bypass and natural transition. Detailed inlet boundary conditions are reported in Table (5.8).

Re T u in (%) ν t in /ν T u(%) T3A-10 7 1.00135 25 0.6465 Table 5.8: Flat Plate: T3A-. Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%). As the flow transitions, the differences between the two numerical predictions become evident. γ predicts a noticeable smaller transition region extent. The increase of C f is considerably steeper compared to γ -Re θ results and experimental measurements. This very steep gradient results in a significantly higher numerical uncertainty at this position, due to the non-monotonic convergence of the solutions on the five tested grids. The turbulence intensity contours in the proximity of the transition point are shown in The discussion on the flat plate test cases has highlighted how the two models perform differently, because of the different correlations that account for transition. The different models performances on the flat plate for these boundary conditions were already known at the time of The following sections are devoted to flow around airfoils at incidence and for low Reynolds numbers. The multivariate metric for validation is an important tool, that allows to point out where improvements are required, but it is useful only when experimental and input uncertainties are clearly defined in the experimental database used for the validation. For this reason, E mv will not be used for the next test cases.

Eppler387

The second configuration chosen for the validation of the transition models performance for 2D configurations is the Eppler387 foil. This geometry is frequently used for sailplanes. This validation and verification exercise was performed in order to assess the liability of the transition models to predict separation-induced transition. As discussed in Chapter 2, this kind of transition process occurs when the laminar boundary layer separates under the effect of the pressure gradient. The flow transitions within the separated shear layer as a result of an inviscid instability mechanism. The flow, eventually, reattaches and the level of turbulence intensity in the boundary layer determines the length of the separation bubble. For this test case, numerical results from different participants of the collaborative group have been presented by Luis Eça at the V&V Verification and Validation Symposium held virtually in Baltimore in May 2020. The initial conditions used for this exercise are reported in Table (5.9). Within the group, two different angles of attack were tested α = 1 • and α = 7 • , for a Reynolds number Re = 3 × 10 5 , based on the chord. However, in the following, we are going to discuss only the numerical results for α = 1 • . We only mention results for α = 7 • . The experimental pressure distribution on the surface of the airfoil is taken from McGhee et al., [START_REF] Mcghee | Experimental results for the Eppler 387 airfoil at low Reynolds numbers in the Langley low-turbulence pressure tunnel[END_REF], and Cole et al., [START_REF] Cole | Experimental measurements of the laminar separation bubble on an Eppler 387 airfoil at low Reynolds numbers[END_REF]. Both experimental campaigns were conducted to study the laminar separation bubble that this particular geometry exhibits at low Reynolds numbers. For a Reynolds number of 3 × 10 5 and for the angles of attack ranging from α = -1 • up to α = 3 • , measurements show a laminar separation bubble, which moves upward to the leading edge as the angle of attack is increased. For these conditions, the flow on the lower surface remains laminar. At higher angles of attack than α = 7 • , no laminar separation is experimentally observed, and the flow naturally transitions on the upper surface. The airfoil shape definition follows the geometry details given in McGhee et al.. The experimental uncertainty is set equal to 0.25%, based on the available experimental information.

Computational Domain

The computational domain extends 12c upstream of the leading edge, which is placed at x/c = 0, and 23c downstream of the trailing edge, x/c = 1. Slip conditions are imposed on the top and bottom boundaries, placed at y = ±12c. The mesh originally consists of a multiblock structured grid. The blocks next to the leading edge form a C-topology that extends slightly downstream of the trailing edge. A second 'layer' of blocks that surround the C-shape gives rise to an O-block topology. The angle of attack is imposed by rotating the geometry with respect to its center. The combination of C and O topology has the advantage of preventing the propagation of very thin cells from the boundary layer through the wake, that would be observed using C-grid topology. The presence of cells with high aspect ratio in the wake would strongly penalize the iterative convergence, in particular, when using a segregated solver for the pressure-velocity coupling, as in ISIS-CFD. On the other side, combining the two approaches, we avoid the possible misalignment between the faces and centers of the cells in the wake. The computational domain is shown Fig. (5.8), as well as a close-up around the airfoil, at the leading edge and in the wake, in order to show the mesh topology. Details on the meshes are given in Table (5.10). For each grid, we report the y + max value, as computed with γ transition model. Similar results are obtained for γ -Re θ .

Re

T u in ν t in /ν T u 3 × 10 5 1% 0.003 0.0078% Table 5.9: Eppler387: α = 1 • , 7 • . Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%). Numerical results computed by the two transition models are compared to measurements. We also show the prediction by k -ω SST (2003) turbulence model, that totally fails to predict the laminar separation that occurs on the upper surface of the foil. Taking the experimental pressure distribution as a reference for the identification of the "critical" locations, laminar separation occurs ahead of the pressure plateau, at about 45% of the chord. Transition occurs as the pressure gradient starts increasing behind the separation point at about 70% of the chord. Eventually, the flow reattaches at the end of the transition region, around 75% of the chord. Compared to γ, γ -Re θ model predicts a slightly longer bubble. Even though the laminar separation occurs around the same position for both models (x/c ∼ 0.45) the transition and reattachment points are pushed towards the trailing edge using γ -Re θ . This difference is related to the handling of separation-induced transition. In γ implementation, an additional turbulence production term P lim k is added in the transport equation for the turbulence kinetic energy. In γ -Re θ , separation-induced transition is accounted for by the definition of γ eff , which is multiplied to the turbulence kinetic energy production term in the k equation. Being γ eff = max(γ, γ sep ), this effective intermittency is bigger than 1 when separation occurs, i.e when Rev 3.235Re θc > 1. As this limit is met, the production of turbulence kinetic energy is accelerated in order to compensate numerically the low level of turbulence kinetic energy. Comparing the two transition model results, we observe that the turbulence kinetic energy production term as constructed in γ is stronger than the one in γ -Re θ , resulting in a smaller bubble extent. γ predicts the same transition and reattachment position as in the experiments, but it overpredicts the C p value at this location. This behavior is not justified by none of the reported uncertainties. Overall, the validation comparison error E for the pressure distribution is relatively small. Outside the separated region, it can be observed that the numerical results deviate quantitatively from the measurements in the laminar region. Even if the input uncertainties have not been calculated, this deviation is unlikely to be related to the eddy viscosity ratio, which is considerably small. It is most probably associated to the actual angle of attack used in the experimental set-up. A small variation of the angle of attack of ±0.1 • would affect the pressure distribution around the foil in the first half of the chord, where the flow is laminar. Nevertheless, this overestimation does not affect the transition process, which is determined by the pressure gradients. C f distribution, as predicted numerically, is shown in Fig. (5.9b). For this quantity we do not dispose of measurements. From this perspective the differences between the two transition models are even clearer, namely the shifting in the transition locations, min(C f ), and in the reattachment point, i.e. the position at which the computed C f crosses zero and becomes positive again. For both C f and C p the numerical uncertainty is negligible, given the extremely fine grids used, except in the transition region where it is affected by the steep gradients. The overall dimensionless forces coefficient predicted by γ and γ -Re θ , and k -ω SST (2003) are reported in Table (5.11). Clearly, the different features of the laminar separation bubble predicted by the two transition models significantly affects the drag coefficient C d . The larger extent of the bubble predicted by γ -Re θ contributes to an higher pressure component of the drag coefficient C d,p . On the other side, γ predicts an higher viscous component of the drag, C d,f , because of the shorter predicted laminar region. Significant improvements with respect to the turbulent solution are observed. Indeed, a turbulence model overestimates significantly the drag coefficient, because it does not account for the laminar flow on the lower surface and the laminar separation bubble. Contrary to the drag, the lift coefficient is not affected by the used transition model. Transition models predict an higher lift, because of the different bubble extension with respect to measurements. We report hereafter some of the results of the AVT-313 Workshop devoted to 2D configurations. The figures are taken from the presentation of Luis Eça at 45th AVT Panel Business Meeting held in Trondheim in October 2019. Fig. (5.10a) and (5.10b) show the pressure distribution around the foil as computed using γ and γ -Re θ models and different CFD solvers. Fig. (5.11) and (5.12) show the lift and pressure and friction drag coefficients plotted against the grid step size r i = h i /h 1 . All the results were computed on the same grids using the boundary conditions reported in Table (5.9). All the solvers predict very similar C p distributions around the foil, independently of the used transition model. Nevertheless, from the performed grid convergence study on the forces, it can be noticed that different solvers converge towards different values. In particular, this difference is important for γ -Re θ results. γ implementation appears to be more robust among the different solvers.

Further analysis for this configuration have revealed the dependency of the numerical predictions on the turbulence kinetic energy production limiter used within the TKE transport equation. The production limiter controls the growth of turbulence at the stagnation point, further affecting its evolution at the edge of the boundary layer. Different production limiters can be used within the k -ω SST turbulence model, as found in literature. Other than the Kato-Launder formulation, which is used in all the computations presented in this dissertation, and defined as:

P k = µ t SΩ, (5.36)
the production term can also be limited using the following formulation High eddy-viscosity ratio imposed in the free-stream can amplify the effect of the production limiter, affecting also the transition location up to a 10% of differences among results with different limiters, as discussed in the doctoral thesis of Lopes. The production limiter happens to have more impact on γ -Re θ results than γ, because it affects Re θt in the free-stream and its diffused value within the boundary layer. A similar behavior of γ -Re θ was also observed in Diakakis et al., [START_REF] Diakakis | Assessment of transition modeling for high Reynolds flows[END_REF].

P k = max(µ t S 2 , aβkω). ( 5
Despite the important impact of the production limiter on the transition process, it cannot be asserted which limiter is the best to use and the choice is up to the developer. Nevertheless, both transition models have been originally calibrated using Kato-Lauder formulation, and that is where our choice lies.

α = 7 •

To conclude our analysis on the flow around the Eppler387 foil, we show the computed pressure coefficient computed by γ for α = 7 • in Fig. (5.13). As mentioned at the beginning of this section, for high angles of attack none of the experiments showed any laminar separation.

The flow transitions naturally on the upper surface, at about 40% of the chord, as mentioned in Cole et al.. This flow topology is not predicted numerically. For the boundary conditions of Table(5.9), γ and γ -Re θ (not shown) predict a laminar separation. The flow separates at about the same position where it is experimentally supposed to transition. This behavior of the numerical solution is not unexpected, indeed, the turbulence kinetic energy undergoes such a strong decay ahead of the foil, that the value of TKE at the leading edge is very small. Such a small value is "forgotten" by the mathematical model and it does not play any role in the transition process. For such low values of free-stream turbulence intensity, the transition models cannot predict any transition mechanisms other than separation induced, because γ ef f and P k l im will activate to compensate the low value of turbulence intensity in the free-stream. We need to consider that we do not know the experimental value of the free-stream turbulence intensity used in the wind tunnel, other than the formulation at the lowest turbulence intensity as possible. For this reason the numerical results might deviate from the measurements. Other solvers perform as ISIS-CFD for these boundary conditions.

NACA0015

The third configuration that was selected to pursue the validation of the transition model for 2D flows is the airfoil NACA0015. This geometry was chosen following the exchange of new experimental data among the participants of the AVT-313 collaborative group. The results for the angles of attack α = 5 • and α = 10 • were chosen as test case for the 2D "AVT-313 Incompressible Laminar-to-Turbulent Flow Transition Study Comparison Workshop". In this dissertation, we are only going to discuss the results for α = 5 • . Hereafter, we give additional detail on the experimental technique and performed measurements.

Experimental Set Up

The experiments, presented in Miozzi et al., [START_REF] Miozzi | Skin friction and coherent structures within a laminar separation bubble[END_REF], were carried out in the CEIMM cavitation tunnel, featuring a square test section of side 600mm, at CNR-IMM in Rome. The model has a chord length of 120mm and the span equal to the test section side.

The model was mounted in the middle of the tunnel test section. It is allowed to rotate around its geometrical center in order to simulate different angles of attack, α = 3 • , 5 • , 7 • , 10 • . The Reynolds number is Re = 1.8 × 10 5 , based on the chord c of the foil and on a free-stream velocity U ∞ = 1.5m/s. This set of experiments was run in order to analyze and characterize the laminar separation bubble which is observed on the upper surface of the foil. To this purpose, the skin friction vector field is needed to reconstruct the flow topology, in order to understand the boundary layer laminar-to-turbulent transition. The skin friction streamlines are used to identify the separation and reattachment points, following the theory developed by Surana, [START_REF] Surana | Exact theory of three-dimensional flow separation. Part 1. Steady separation[END_REF]. The model is coated with Temperature-Sensitive Paint(TPS), a thin luminescent polymer used to measure temperature fields on the surface. A detailed description of the technique is presented in Liu et al., [START_REF] Liu | Pressure-and temperature-sensitive paints, encyclopedia of aerospace engineering[END_REF]. Skin friction data are derived from the temporal evolution of the temperature maps on the surface. These maps contain footprints of the coherent structures that appear and act next to the wall, under the hypothesis that the wall heat flux is mainly due to the convective action of the flow, as discussed in Miozzi et al., [START_REF] Miozzi | Skin friction and coherent structures within a laminar separation bubble[END_REF]. Transition location at the wall within the laminar separation bubble is identified by the maximum gradient of the streamwise temperature distribution ∂T /∂x| max , for ∂T /∂x < 0. This location corresponds with the negative kink of the time averaged skin friction distribution. Decreasing the angle of attack α, different flow regimes are identified. For high angles of attack, the flow is quasi-steady and the size of the LSB is rather small. For this case, two dimensional instabilities appear to be predominant in the flow separation, transition and reattachment. As the angle of attack is decreased, i.e. for weaker adverse pressure gradient, the flow regime is dominated by intermittent and wedge-shaped structures. Coherent structures with a more emphasized three-dimensional shape are observed within the bubble.

Computational Domain

The computational domain was selected in order to match the experimental tunnel characteristics. The mesh domain is a rectangle of length 21.7c and width 5c, the center of the foil (half-chord) is located at the origin of the reference system and the axis x is aligned with the incoming flow. The inlet is located 7.2c upstream of the center of the airfoil and the outlet is located 14.5c away from the trailing edge. Top and lower boundary are treated as slip walls and they are placed at y = ∓2.5c. The inlet conditions are reported in Table (5.12), as well as the value of turbulence intensity in the vicinity of the leading edge.

Re

T u in ν t in /ν T u 1.8 × 10 5 0.5% 0.0018 0.01% Table 5.12: NACA0015: α = 5 • . Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%).

The meshes topology is the same as the one described for Eppler387, but for a lower grid refinement level. Details of the meshes are given in Table (5.13). For each grid, we report the y + max value as computed with γ transition model. In Fig. (5.14) we show the computational domain size, and a close-up of the mesh topology around the airfoil and in the wake. For the coming analysis, we focus only on the upper surface of the foil. The experimental time-averaged skin friction distribution on the upper surface of the foil is shown in Fig. (5.15). In the separated region, the C f distribution evolves from weaker to more intense negative values. This change represents the passage from the dead water region, see Fig. (2.9), influenced by transition happening in the separated shear layer, to the reattachment zone. The simulations were run on a 2D configuration, and, with the purpose of comparing experimental to numerical results, the measured C f profile was averaged in the spanwise direction.

N cells N cells-foil r i = h i /

NACA0015

The reported experimental uncertainty corresponds to the standard deviation of this average. Let us notice first that, even if the grid refinement level is smaller, the numerical uncertainties are negligible, except in the transition region. The low numerical uncertainty justifies the use of coarser meshes, always keeping a y + max < 1. As a general trend, it can be observed that both transition models under predict the skin friction all along the upper surface of the foil and this deviation is justified neither by the numerical nor the significant experimental uncertainties. The experimental laminar separation bubble is smaller compared to simulations. The numerical transition and reattachment points are pushed towards the center of the foil. Nevertheless, the extent of the transition region as predicted by computations is in good agreement with experiments. The negative peak, which is predicted by both transition models, might not be visible in the experimental distribution as an effect of the average in the spanwise direction. This discrepancies are related to the fact that a threedimensional experiment has been numerically treated as two-dimensional. In addition, at the time of the measurements, there was no technique that allowed to measure negative values of skin friction on the surface, which are instead recovered in the post processing, introducing additional experimental uncertainties. The major difference between experiments and simulations is observed in the turbulent region, from x/C ∼ 0.2, where the increase of C f towards the trailing edge is not reproduced by the simulations. The behavior of the measured C f is most probably related to some uncertainty of the experimental results in this second-half portion of the foil. Such an increase of C f should be justified by a favorable pressure gradient, that seems unlikely to be found in this region of the foil. The estimation of C f might be affected by the impossibility of keeping the heat flux uniform for x/C higher than 0.15, because of the thinning foil profile, as mentioned in Miozzi et al., [START_REF] Miozzi | Taylor hypothesis applied to direct measurement of skin friction using data from temperature sensitive paint[END_REF]. In order to compare the transition models prediction among them, we present as well the pressure coefficient C p distribution on the upper surface, computed by γ -Re θ and γ and k -ω SST (2003) in Fig. (5.17). As for the C f profile, differences are observed in the separated region. γ -Re θ predicts a smaller transition region extent, compared to γ, whose predicted transition point is slightly pushed upwards compared. The reattachment points are nearly coincident. The shifted transition positions leads to an more extended laminar region, and the overall viscous drag coefficient, as predicted by γ -Re θ model, is expected to be smaller than γ predictions.

NACA0015

Nevertheless, this shifting is overcompensated by the higher skin friction predicted by γ -Re θ behind the reattachment point. Finally, the viscous drag components predicted by the two models are very similar, despite the differences in the C f profile, as it can be observed by comparing the forces, shown in Table (5.14). This higher value of C f predicted by γ -Re θ downward the reattachment zone comes from the stronger turbulence intensity T u predicted by γ -Re θ . Contours of T u in the vicinity of the reversed flow region, as indicated by the velocity streamlines, are plotted in Fig. (5.18). The production term of turbulence kinetic energy in γ -Re θ is the double compared to the one of γ. This results in higher values of turbulence intensities as the flow reattaches and higher values of local skin friction. This behavior is different from what was observed for the Eppler387, where the low values of turbulence intensity delayed the transition process in γ -Re θ compared to γ results. Nevertheless, if we refer to the numerical results, the adverse pressure gradient around the NACA 0015 is stronger than around Eppler387, and it compensate the low free-stream turbulence intensity. Despite the quantitatively different values of skin friction and pressure coefficient predicted by the two transition models, the simulated flow physics is very similar and the resulting overall forces differ of less than 0.15%.

C d × 10 3 C d,p × 10 3 C d,f × 10 3 C l γ -

Computational costs

In the presented discussion, we have analyzed the performance of the two correlation based transition models for 2D transitional flows. Transition models are clearly needed to predict features that a turbulence models fails to represent, for instance the laminar-separation bubble around airfoil at low Reynolds numbers. Nevertheless, these modeling capabilities come at a price. The number of non-linear iterations up to convergence for each test case and and each transition model is reported in the histogram in Fig. (5.19). We consider the simulations run on the coarsest mesh of the set used for each test case. Indeed, on each finer grid, we used the interpolated solution from its coarser grid as initial solution. This technique allows to save computational time. In order to consider a simulation converged, we require the normalized residuals in L 2 norm to gain at least four orders, and the flow solution and the forces to converge, i.e. the solution does not change anymore at each nonlinear iteration. Compared to turbulence models, both transition models demand an elevated number of non-linear iterations. γ -Re θ model needs more non-linear iteration compared to γ because of the additional transport equation for the momentum thickness Reynolds number. This behavior is expected, because the coupling between the transition and turbulence equations introduce additional difficulties in converging the flow solution. In particular, on grids as dense as those presented in this chapter, the solution requires time to converge. Very dense grids in the wall-normal direction are needed in order to capture accurately transition in the boundary layer. This constraint results in highly stretched grids with high aspect ratio cells in the vicinity of the walls, which affects negatively the convergence of the turbulence quantities in the boundary layer. Clearly, reducing the mesh size, we can gain computational time. For instance, the flow around Eppler387 and NACA0015 airfoils presents a very similar, unsteady topology. Using meshes twice coarser for NACA0015, compared to Eppler387, the number of non-linear iterations are considerably reduced in half. In addition, the use of such fine grids for laminar flows, which do not require such spatial resolution, affects negatively the robustness of the solver and the overall convergence rate. For instance, for the T3A and T3A-tests, the meshes are the same and the flow under consideration does not present a particular topology. Nevertheless, the T3A-requires a significant high number of non-linear iterations, because of the larger extent of the laminar region, that affects the simulation convergence. The high number of non-linear iterations is also associated to a noisy iterative convergence: the normalized residuals decrease rapidly at the beginning of the simulation, start oscillating and stall. The evolution of the residuals in L The iterative convergence is disturbed on one hand by the presence of the several min and max limiters within the transition models formulation. These functions, as for example F onset,3 , of Eq.(3.21) and shown in Fig. (5.21), are continuous, but not smooth, and do not promote stability. F onset,3 , as part of the F onset functions, activates or totally switches off, i.e. the min is over zero, the production of intermittency P γ within the boundary layer. Thus, P γ might start oscillating between low values and zero, affecting not only the residuals of γ, but also of the other transported variables. On the other hand, the periodic behavior of the residuals is also related to the flow physics under consideration. For instance, for NACA0015 the steady approximation is inadequate for an unsteady phenomenon, such as flow separation, and clearly affects the robustness of the solver.

Conclusion 2D Test Cases

Within the AVT-313 group, the two transition models were tested for three different geometries: ERCOFTAC flat plate, Eppler387 and NACA0015. These simulations were run to asses the models performance for 2D transition mechanisms, from bypass to separation-induced transition. Over all, the two correlation-based transition models perform very well. They are capable of predicting flow features that would be totally neglected if turbulence models were used. For instance, for both Eppler387 and NACA0015 profiles, a transition model is able to simulate the laminar-separation occurring on the upper surfaces of the foils at incidence and at low Reynolds number. The presence of a laminar separation bubble strongly affects the pressure distribution around the foils, and leads to a significant increase of the pressure drag. A RANS simulations able to characterize such kind of phenomenon is an important tools. For instance, it can be used as a reference for the placements of flow control devices, such as turbulators, to improve the airfoil performance at low Reynolds numbers. A thorough investigation of the two transition models performance has revealed that the correlations that account for transition, as differently defined within the two models formulation, lead to different predictions. In particular, the main difference lies in how the model reacts to the imposed free-stream/inlet conditions. The value of eddy-viscosity ratio in the free-stream appears to be a crucial parameter, that can strongly affects the transition process. As a matter of fact, the multivariate metric exercise performed on the flat plate test case T3A has shown how different values of eddy viscosity ratio at the inlet can affect the transition solution. γ -Re θ appears to be more influenced by the free-stream eddy viscosity ratio R T values more than γ. Indeed, R T controls the decay of turbulence intensity ahead of the body, which affects Re θt computed in the free-stream, possibly accelerating or delaying the transition process. γ performs differently from γ -Re θ , because in its formulation are used only local quantities. Despite their good modeling performance, the transition models are not entirely robust and require a significant number of non-linear iterations to convergence. This behavior is attributable to the use of very dense grid that affects the convergence of the turbulence quantities within the boundary layer. These meshes are also deleterious in term of convergence rate, when the solution presents an extended laminar flow region.

The robustness of the flow is also affected by the several non-smooth min, max functions within the transition equation formulations. These functions cause the non-linear residuals to stagnate and oscillate. To this matter, Lopes, in his doctoral thesis, has proposed smoothing approximations for these functions by exponential penalty functions, as proposed in Piotrowski et al., [START_REF] Piotrowski | Investigation of a local correlation-based transition model in a Newton-Krylov algorithm[END_REF]. This approach prevents the non-linear residuals from stalling, nevertheless, its use requires a re-calibration of the models formulation. Being γ and γ -Re θ models based on correlations built empirically, every change in their formulation can affect the transition process. For this reason, we opted not to implement this strategy in ISIS-CFD. Unfortunately, such fine grids are needed for accurate predictions, especially in the transition region, where, despite the mesh density, the numerical uncertainties are still not negligible because of the steep gradients that the flow variables are subjected to. At this stage, it is not possible to define which approach is the most reliable, because it is arduous to perform a rigorous validation exercise in lack of detailed experimental measurements. We are missing experimental uncertainties and, above all, of well-defined inlet/free-stream conditions, in order to characterize the turbulence decay ahead of the body. Nevertheless, for the grid refinement level used for the presented computations, our preference lies with the one-equation γ model. Its formulation is simpler and it presents a faster non-linear convergence. In addition, based on the results for two-dimensional configurations, it has been proven to be the most robust formulation between the two correlation-based models. For this reason, the next chapters will be mainly devoted to γ transition model.

Chapter 6

Crossflow Criterion Calibration

We have discussed and analyzed the performance and accuracy of the local-correlation based transition models γ and γ -Re θ against two-dimensional test cases. The two transition models were implemented as proposed in their original publications. In none of γ -Re θ publications, the doctoral thesis of Langtry, and the article from Menter & Langtry of 2012, nor in the one for γ by Menter et al. in 2015, any criteria other than streamwise transition and separation-induced transition was presented. Nevertheless, moving to three dimensional tests cases, further modifications are necessary to account for the different behavior of 3D boundary layers. Examples of 3D flows of interest for transition are the swept wing, rotating disks and rotating cones. They all present a characteristic stability behavior: the presence of streamwise vorticity within the boundary layer that produces a strong spanwise modulation of the basic state and gives rise to secondary instabilities, Saric, [START_REF] Saric | Stability and transition of three-dimensional boundary layers[END_REF]. In the following, we are going to focus our attention on stationary crossflow transition, that is the main mechanism which a 3D boundary layer is subjected to in a low free-stream turbulence environment. We aim to define a local crossflow criterion to be included in the transition models equations. A correlation that takes into account these instabilities should be a function of the multiple parameters that play a role in the transition process. First of all the pressure gradient in the stream-wise direction, the Reynolds number, based on the boundary layer thickness, and the crossflow strength are essential. Turbulence intensity and surface roughness are also crucial, but they are not considered in the following discussion. Within γ -Re θ formulation, transition due to crossflow waves can be accounted for either by modifying the production term of the transported Re θt , Eq.(3.28), or by modifying F onset in the transport equation for the intermittency γ, Eq.(3.17). In the first case, the aim is to lower Re θt value within the boundary layer, while using the second approach the production of intermittency is enhanced. Within γ model, only the second strategy is possible. Independently of the strategy, a crossflow criterion is needed. In the last six years, different efforts have been made by research groups to build a crossflow criterion in the LCTM frame. The majority rests on the well-known C1 crossflow criterion, proposed by Arnal et al. in 1984, [16], and based on the crossflow Reynolds number Re δ2 :

Re δ 2 = U 1e δ 2 ν where δ 2 = - ∞ 0 W 1 U 1e dy. (6.1)
Referring to the notations in Sec.(2.2), U 1 and W 1 are the streamwise and crosswise velocity components. U 1e = (U 2 e + W 2 e ) is the velocity at the edge of the boundary layer. Menter & Smirnov, [15], proposed in 2014 the T c1, a C1-based criterion, built to extend γ transition model to CF predictions. This model variant was presented in occasion of the 19th STAN/DGLR Symposium in 2014, but never officially published. A non-official version of this publication was shared by Menter with the author of this dissertation. The T c1 is a local approximation of the C1 criterion, evaluated through auxiliary functions expressed by the solution of the Falkner-Skan-Cooke (FSC) equations. These are an extension of the two-dimensional Falkner-Skan equations by considering the spanwise component of the velocity for a flow over an infinitely yawed wedge at zero angle attack, Cooke, [START_REF] Cooke | The boundary layer of a class of infinite yawed cylinders[END_REF]. These equations, that will be introduced in Sec.(6.2), are solved under the assumption of zero-pressure gradient in the spanwise direction ∂ ∂z = 0 . In the approach proposed by Menter & Smirnov, the wall-normal change of the normalized vorticity vector is used as an indicator of the crossflow strength. This quantity preserves the Galilean invariance of γ model. Another C1-based criterion is the local C1 approach proposed by Grabe, [START_REF] Grabe | Transport modeling for the prediction of crossflow transition[END_REF], and suggested as CF extensions of γ -Re θ model. This crossflow criterion is based on the ratio of the maximum crossflow vorticity Reynolds number to the displacement thickness Reynolds number in the crosswise direction. To our best knowledge, this criterion has been implemented only in DLR-TAU-Code solver and it has not been successfully used for non-swept wing geometries. For these reasons, it was not further studied by the authors of this dissertation.

In this chapter, we are going to introduce first the C1 criterion by Arnal and the Falkner-Skan-Cooke equations. Then, the T c1 criterion, as originally formulated and calibrated on FSC solutions by Menter & Smirnov, will be discussed. Finally, the T c1 new re-calibration, proposed by the author of this dissertation, is presented.

Crossflow Transition: C1 criterion (1984)

The first researcher who proposed a crossflow criterion based on the crossflow Reynolds number Re δ 2 was Beasley in 1973, [START_REF] Beasley | Calculation of the laminar boundary layer and prediction of transition on a sheared wing[END_REF]. Beasley empirically found as Reynolds number value at the transition onset Re δ 2 t = 150. Later on, in 1984, Arnal found that crossflow transition does not occur for a unique value of Re δ 2 . Instead, the Reynolds number at transition onset varies as a function of the longitudinal shape factor H 12 , which accounts for the streamwise pressure gradient. In a three-dimensional boundary layer, H 12 is defined as:

H 12 = δ * θ , (6.2)
δ * = δ 0 1 - U 1 U 1e dy, (6.3) θ = δ 0 1 - U 1 U 1e U 1 U 1e dy, (6.4)
with δ * and θ, displacement and momentum thickness, respectively. In [START_REF] Arnal | Théorie de l'instabilité laminaire et critères de transition en écoulement bi et tridimensionnel[END_REF] The upper bound of H 1 = 2.7 is imposed, because no crossflow transition is expected for H 1 > 2.7, where transition is dominated by longitudinal instabilities. Re δ 2 t expression is based on experimental results at low T u from Poll, [START_REF] Poll | Some aspects of the flow near a swept attachment line with particular reference to boundary layer transition[END_REF], around a cylinder of large diameter, Boltz et al., [START_REF] Boltz | Effects of sweep angle on the boundary-layer stability characteristics of an untapered wing at low speeds[END_REF], around NACA64 2 A 015 symmetric profile, and on the measurements around the ONERA D profile by Schmitt et al., [START_REF] Schmitt | Ecoulements subsoniques et transsoniques sur une aile en flèche variable[END_REF], and at ONERA/CERT/DERAT, [126]. As mentioned at the beginning of this chapter, the C1 criterion is widely used for the prediction of crossflow transition, nevertheless it is necessary to mention that this criterion does not take into account the receptivity stage. No information about the turbulence intensity in the free-stream or the surface roughness are included in its formulation. In this regard, in the same publication, Arnal proposed a second criterion, the so-called C2. This one is based on the momentum thickness Reynolds number relative to the velocity profiles projected onto the direction of the wave vector. This new criterion results in an envelope of curves Re δ 1 t = f (H 12 , T u). Its application, however, requires the knowledge of the longitudinal and transversal velocity profiles in the boundary layer and it cannot be easily implemented in a CFD solver. For this reason, it will not be discussed further.

To be implemented within a RANS transition model, the C1 criterion is reconstructed via the FSC equations. These are the key to define the link between non-local quantities such as the crossflow Reynolds number and the local quantities needed to build a local formulation. The FSC equations are described in the following section.

Falkner-Skan-Cooke

Let us consider the three-dimensional boundary layer equations system for an infinite swept wing, i.e. the derivative along the span is zero, reported in Eq.(2.25)-(2.27). We suppose that the inviscid chordwise velocity at the boundary layer edge U e follows a potential law over the coordinate normal to the leading edge x and that the spanwise velocity W e is constant, such that:

U e ∼ U ∞ x L m , (6.6
)

W e = const. (6.7)
L is the characteristic length and U ∞ the reference velocity. m is the streamwise pressure gradient and it is defined as:

m = x U e dU e dx . (6.8)
Once defined the Blasius similarity variable η as:

η = y U ∞ (m + 1) 2νL x L (m-1)/2 , (6.9)
and introduced the stream function Ψ:

Ψ = 2U ∞ νL m + 1 x L (m+1)/2 f (η), (6.10)
with U = ∂Ψ ∂y , V = -∂Ψ ∂x , we observe that the continuity equation is automatically satisfied. Eq.(2.26) becomes:

f + f f + β(1 -f 2 ) = 0, (6.11)
where β is the Hartree parameter associated to m by the relation:

β = 2m m + 1 . (6.12)
The dash in Eq.(6.11) denotes the differentiation with respect to η.

On the other hand, given W = W e g(η), Eq.(2.27) becomes:

g + f g = 0. ( 6.13)
The system of equations given in Eq.(6.11)-(6.13) are the Falkner-Skan-Cooke equations. The corresponding boundary conditions are:

f, f , g → 0 for η → 0, (6.14)
and

f , g → 1 for η → ∞. (6.15)
The solutions f and g can be combined into the dimensionless streamwise and crosswise velocity components, non-dimensionalized with respect to the velocity magnitude at the edge of the boundary layer

U 1e = U 2 e + W 2 e : U 1 /U 1e = f cos(φ) 2 + g sin(φ) 2 , (6.16
)

W 1 /U 1e = (g -f ) cos(φ) sin(φ). (6.17)
φ is the sweep angle, i.e. the angle of the inviscid flow direction with respect to the chordwise direction at the edge of the boundary layer. It is constant along the wall normal height of the boundary layer, and it is defined as:

tan(φ) = W e U e . (6.18)
According to the definition in Eq.(6.17), W 1 = 0 for φ = 0 • and φ = 90 • , but also for zero pressure gradient flows β = 0, because g -f = 0. The last condition does not occur in real physical flows, because crossflow velocity also exists for zero-pressure gradients 3D flows.

The FSC equations will be used in the following for the calibration of the T c1 crossflow criterion.

A solver for their resolution was kindly provided by Prof. Christophe Airiau, who is gratefully acknowledged. The source code is included into the FSC project and it is derived from the FundAeroSuite < 2018 >. The code, as provided in April 2021, is written in fortran90. 4th order Runge-Kutta integral scheme is used for the resolution of the main equations.

Menter T c1 Crossflow Criterion

In this section, we introduce the T c1 criterion as first published by Menter & Smirnov, in [15]. Exploiting the provided FSC solver, we have been able to retrace the steps of the calibration performed by Menter & Smirnov. By analyzing their work, we will discuss the rationale behind our re-calibration of the T c1 criterion. The C1 criterion can be rewritten as:

Re δ 2 f (H 12 ) 150 = 1, (6.19)
where f (H 12 )150 is Re δ 2 t of Eq.(6.5). Transition occurs when the ratio of the crossflow Reynolds number Re δ 2 , Eq.(6.1), to the crossflow Reynolds number at transition onset Re δ 2 t , Eq.(6.5), reaches a value of 1 or beyond. Menter & Smirnov propose to split the C1 criterion in three different terms, in order to identify and quantify each parameter that affects the transition process. The C1 criterion can then be rewritten in the form:

Re δ2 f (H 12 ) ∼ F (H 12 )XRe stream . (6.20)
The function F (H 12 ) = 1 f (H 12 ) takes into account the pressure gradient in the streamwise direction. X is a measure of the crossflow strength and Re stream is the Reynolds number relative to the streamwise velocity component. The latter is taken to be the maximum value of the vorticity Reynolds number in the boundary layer, Re Vmax , as in the original γ model's formulation for T-S and bypass transition. Based on Eq.(6.20), the local approximation of the C1 criterion, referred to as T c1, can be written as:

T c1 = 1 150 G • Ψ • Re Vmax , (6.21)
where

Ψ ∼ X = Re δ2 Re Vmax , (6.22)
G ∼ F (H 12 ) = 1 f (H 12 ) . (6.23)
The T c1 criterion, as formulated above, is mathematically equivalent to Arnal's C1 criterion. Let us discuss in detail each of the quantity that appears in Eq.(6.21). The function Ψ approximates the ratio of the crossflow to the streamwise strength. The idea proposed by Menter & Smirnov is to use the wall-normal change of the normalized vorticity.

This physical quantity describes the three dimensionality of the boundary layer. Indeed, Ψ = 0 for 2D flows, because ∂ ∂z = 0 and w = 0. Ψ is defined as:

Ψ =| - → ψ | •y, (6.24)
where y is the wall normal distance. The components of the vector -→ ψ = {ψ i } are defined as:

ψ i = ∂ω i ∂x j n j , where ω i = ω i | - → ω | . (6.25)
The scalar quantity Ψ can be interpreted as an indicator of the crossflow strength because it is proportional to the local change of the flow angle. Once defined Ψ, it has to be assessed how well it represents X. To this purpose, the FSC equations are solved in the parameter range:

0 <β ≤ 1, (6.26)
0 • <φ < 90 • . (6.27)
β is restricted to positive values, because crossflow instabilities occur for accelerated flow, for favorable pressure gradient. β = 1 is the case of 90 • wedge, the 2D stagnation flow. It is the highest possible acceleration parameter. φ is the sweep angle of Eq.(6.18), see Fig.

(2.7) as reference. In Menter & Smirnov, the upper bound for φ is set to 60 • . This is a reasonable limit for most of aerodynamic flows. However, if we aim to extend this criterion to other geometries that can present higher sweep, such as the 6:1 prolate spheroid geometry, treated in the next chapter, this limit has to be extended to its maximum value. All the quantities are evaluated at the particular point η = η max , such that the T c1(η max ) reaches its maximum value in the boundary layer. The quantity Ψ is approximated considering only the derivatives with respect to the normal direction η expressed as a function of y. The only vorticity components that "survive" in the FSC framework are:

ω x ∼ ∂W 1 /U 1e ∂η ∂η ∂y = (g -f ) sin(φ) cos(φ) ∂η ∂y /U 1e , (6.28
)

ω z ∼ - ∂U 1 /U 1e ∂η ∂η ∂y = -f sin(φ) 2 + g cos(φ) 2 ∂η ∂y /U 1e , (6.29)
where we have used the notations given in Sec.(6.2).

The ratio R(β, φ) = aΨ X is shown in Fig. (6.2), for a = 0.4 chosen to match the two indicators for β → 0. The maximum deviation of the two indicators ratio with respect to the targeted value of 1 is about 35%. Nevertheless, it occurs at the corners of the domain and it lies in the limits of the experimental C1 correlation from Arnal, i.e. deviation of the experimental results from the correlation Re δ 2 t . The principal complication arises from the introduction of the function f (H 12 ). R(β, φ) = a Ψ/(XF (H 12 )) is shown in Fig. (6.3), where the new constant a = 0.684 is chosen in order to have a value of ∼ 1 for β → 0. We notice that, at the upper corners of the domain, the ratio R significantly departs from the value of 1. This discrepancy becomes unacceptable and it needs a correction. The approximation is hereafter repeated:

λ θ = - θ 2 ν dv dy ∼ - y 2 ν dv dy . (6.30)
We remind that this approximation is valid at the middle of the boundary layer where θ ∼ δ/2 ∼ y. Using Eq.(6.30), where the derivative dv dy is expressed as in Eq.(3.66), the final definition of the local pressure gradient λ θ,CF is:

λ θ,CF = 0.1111 • λ θ + 2.3, (6.31)
λ θ,CF = min[max(λ θ,CF , 0), 0.7]. (6.32)
The correction in Eq.(6.31) is needed in order to take into account that dv/dy is not zero for zero-pressure gradient flows. The strategy proposed by Menter & Smirnov is to construct a polynomial function in λ θ,CF , to reduce the spread of the ratio R, projected on the β-plane, around the target value of 1. This results in the correction function G(λ θ,CF ): An alternative form for the correction polynomial is also presented in Nichols, [START_REF] Nichols | Addition of a local correlation-based boundary layer transition model to the CREATE T M -AV Kestrel unstructured flow solver[END_REF], as taken from ANSYS-CFX Solver Guide (Release 17.0, March 2018), and implemented in CREATE TM -AV Kestrel unstructured finite volume solver KCFD. In this last publication, the function G is formulated as:

g(λ θ,CF ) = 8.8λ 3 θ,CF -9.1λ 2 θ,CF + 3.7λ θ,CF + 1 g(λ θ,CF ) = min[max(g(λ θ,CF ), 1), 2.3] G(λ θ,CF ) = 0.684 g(λ θ,CF ) (6.
g(λ θ,CF ) = 27864λ 3 θ,CF -1962λ 2 θ,CF + 54.3λ θ,CF + 1, g(λ θ,CF ) = min[max(g(λ θ,CF ), 1), 2.3)], G(λ θ,CF ) = 0.684 g(λ θ,CF ) , (6.34)
where λ θ,CF is defined as:

λ θ,CF = -7.5710 -3 dv dy y 2 ν + 0.0174, (6.35)
λ θ,CF = min[max(λ θ,CF , 0), 0.0477]. (6.36)
The new defined G leads to the same distribution of the corrected ratio R as the one shown in Fig. (6.4b). The upper bound 0.0477 in Eq.(6.36) is the maximum value that the new λ θ,CF , as formulated in Eq.(6.35), can reach if computed from FSC solutions.

Considerations on the crossflow transition criterion as proposed by Menter & Smirnov

Simulations run with this γ CF extension were presented in [15], the same paper where the T c1 transition criterion formulation was presented for the first time. The two tested configurations in the publication were the infinite swept wing NLF(2)-0415 with a geometrical sweep angle Φ = 45 • and angle of attack α = -4 • , and the 6:1 prolate spheroid at α = 15 • angle of attack and Re = 6.5 × 10 6 . Hereafter we discuss their numerical results computed by γ + CF , as presented in their paper. Nevertheless, we are not aware about the strategy used by Menter & Smirnov to include the T c1 within γ model formulation. Results for NLF(2)-0415 are shown in Fig. (6.5). Comparing γ, γ + CF correlation (T c1 criterion) and experimental results from Dagenhart et al., [START_REF] Dagenhart | Crossflow stability and transition experiments in swept-wing flow[END_REF], we observe that the γ transition model with this crossflow inclusion predicts quite good results. The model is able to predict the sensitivity of the transition location (x/c) tr to the Reynolds number, feature that the γ without crossflow criterion is not able to capture. The results obtained around the 6:1 prolate spheroid are presented in Fig. (6.6). γ predictions with and without crossflow extension are compared to the experimental results from Kreplin, [START_REF] Kreplin | Wall shear stress measurements on an inclined prolate spheroid in the DFVLR 3m × 3m low speed wind tunnel[END_REF]. We do not go into detail about the flow physics around the prolate spheroid, because it will be discussed in the following chapter dedicated to 3D numerical results. Nevertheless, we observe that accounting for crossflow effects results in a better prediction compared to results without crossflow inclusion. Menter & Smirnov T c1 criterion promotes transition in a narrower region compared to measurements and expectations. This behavior was not entirely understood at the time of the publication. It was supposed that the behavior was related to the use of meshes not sufficiently refined. In fact, the mesh used for the validation was fairly coarse: a total number of 2.2M elements and a y + < 0.9 at the surface. The need of finer meshes for predicting transition and for computing the 2nd derivative of the velocity, that enters the model through the function Ψ (Eq.(6.24)-(6.25)), is undeniable. However the behavior of the crossflow criterion for the prolate spheroid is also related to its calibration against the Falkner-Skan-Cooke solutions. Compared to a swept wing, the 6:1 prolate spheroid, is a geometry significantly more swept. At the symmetry plane on the windward side, if we exclude the stagnation point, the geometrical sweep angle is Φ = π 2 -α, where α is the angle of attack, as explained by Arnal, [START_REF]Three-dimensional boundary layer: laminar-turbulent transition[END_REF]. For the specific case of α = 15 • , we would have a geometrical sweep angle Φ ∼ 75 • and we expect a local sweep angle φ of this order on the windward side, where crossflow transition can occur. The T c1 criterion from Menter & Smirnov is not calibrated to cover cases for sweep angles higher than φ = 60 • . In order to extend the range of applicability of the T c1 criterion, a new calibration that accounts for higher sweep angle is needed and it is proposed in the following section.

T c1 Crossflow Criterion New Calibration

Recapitulating what has been discussed above, the first limitation of the T c1 calibration performed by Menter & Smirnov is the neglect of sweep angles higher than φ = 60 • . This limit, which is reasonable for swept wing geometries, that hardly experience higher sweep angles, is dictated by the will of correcting the indicator functions ratio R = 0.684Ψ/(XF (H 12 ), shown in Fig. (6.3), by a one variable function. This is the 3rd order polynomial g(λ θ,CF ) defined in Eq. (6.33). Indeed, the two independent FSC parameters, β and φ, are not known in the local formulation. With these requirements, a one parameter function, that reduces the spread of the ratio R = 0.684Ψ/(XF (H 12 )) around the value of one, can hardly be determined because of the different behavior of the surface R at the upper corners of the domain, shown in Fig. (6.3). A possible strategy to achieve a better calibration is, hence, to include in the formulation of the T c1 criterion a local version of the sweep angle φ and a local approximation of the Hartree parameter β. The local sweep angle is introduce to better evaluate the pressure gradient parameter λ θ , which is defined as:

λ θ = θ 2 ν dU e dx = β ∞ 0 U 1 U 1e 1 - U 1 U 1e dη 2 , (6.37)
using FSC notations. We can construct the pressure gradient parameter λ θ as a two parameter function, depending on the sweep angle φ, and dv dy , as defined in Eq.(3.66). The use of dv dy comes as a natural choice to be coherent with the γ formulation for 2D transition mechanisms. The new λ θ is constructed approximately at the middle of the boundary layer, as a 3rd order polynomial in the two variables dv dy and φ. It is obtained using a least squares method fit. The new pressure gradient parameter, here referred to as λ CF , is formulated as: For numerical reasons the value of λ CF is further bounded as

λ CF =0.0473 -0.
λ CF = min max λ CF , 0 , 0.16 . (6.39)
Once constructed the pressure gradient parameter, we can build the non-linear function G(λ CF , φ) that approximates the shape parameter function in the T c1 criterion, Eq.(6.21). The new function G reads as

G(λ CF , φ) =1.992 -0.7328 φ -0.00573 λ CF +0.02344 φ 2 -0.1868 φλ CF -0.08126 λ 2 CF +0.05222 φ 3 + 0.02332 φ 2 λ CF + 0.0490310 -5 φλ 2 CF +0.03326 λ 3 CF , (6.40)
where both λ CF and φ are normalized with respect to their mean value and standard deviation. We need now to build a local approximation of the sweep angle φ. The definition of φ as the local angle between the external potential flow direction, U 1e , aligned with the reference coordinate system (x 1 , y, z 1 ), and the local wing-attached reference system (x, y, z), is not CFDcompatible. Indeed, it would require the definition of the local wing-attached reference system and the identification of the boundary layer edge. Hence, following Högberg & Henningson, [START_REF] Högberg | Secondary instability of cross-flow vortices in Falkner-Skan-Cooke boundary layers[END_REF], the local sweep angle is defined with respect to the reference system (x p , y, z p ), identified by the direction of the pressure gradient vector at each point, ---→ (∇p), and the reference coordinate system (x, y, z), identified through the velocity vector -→ u = (u, v, w). In this new coordinate system,

x p is aligned with the pressure gradient, y is normal to the surface, and z p is perpendicular to the plane (x p , y). Indeed, by FSC assumption, the pressure gradient is zero in the spanwise direction. On the other hand the use of the local velocity vector is an acceptable approximation, because the maximum value of the T c1 criterion is reached close to the boundary layer edge, recovering the original definition of φ.

Based on these assumptions, the local sweep angle φ L definition does not depend on the geometry and, following Choi et al., [START_REF] Choi | Recent improvement of a correlation-based transition model for simulating three-dimensional boundary layers[END_REF], it can be computed as follows:

φ L = arccos - → u wt • ---→ (∇p) wt || - → u wt |||| ---→ (∇p) wt || , φ L = min[φ L , π -φ L ]. (6.41)
-→ u wt and ---→ (∇p) wt are the tangential projection at the wall of the local velocity vector and the pressure gradient. They are computed as:

- → u wt = - → u -- → u • - → d || - → d || - → d || - → d || (6.42)
---→ (∇p) wt = ---→ (∇p) - ---→ (∇p) • - → d || - → d || - → d || - → d || , (6.43)
where -→ d is the wall distance vector. Finally, being ∂p ∂z = 0, because of the infinite swept wing assumption, and ∂p ∂y = 0 within the boundary layer, we obtain:

cos(φ L ) = ∂p ∂x , 0, 0 • u, 0, w ∂p ∂x , 0, 0 u, 0, w = u (u 2 + w 2 ) , (6.44)
that at the boundary layer edge becomes

cos(φ L ) e = u e (u 2 e + w 2 e ) . (6.45)
Eq.(6.45) recovers the FSC definition of the sweep angle φ, as defined in Eq.(6.18).

Inclusion within γ formulation

The recalibrated T c1 criterion is further included in the γ formulation by modifying the F onset function in the transport equation for the intermittency γ, Eq. (3.53). The new F onset,CF that triggers the production of γ is based on the T c1 criterion and summed up to the F onset function of the original formulation. F onset,CF is defined as:

F onset1,CF = (G • Ψ • Re V ) c 150 (6.46) F onset2,CF = min[max(F onset1,CF , 0), 2] (6.47)
F onset3,CF = max(1 -(R T /a) 3 , 0) (6.48)
F onset,CF = max(F onset2,CF -F onset3,CF , 0) (6.49)
a = 1.5, c = 0.6 (6.50)

The functions F onset2,CF and F onset3,CF were defined following the original model implementation, the function F onset3,CF is needed to active the transition criterion beyond the transition point. The constant a by which the eddy viscosity ratio R T is divided in Eq.(6.48) affects the transition location, pushing slightly upward and forward the transition front on the surface. Higher is the constant a, stronger is the limitation on the function F onset3,CF , resulting in a slower transition. Through the proportionality constant c, we account for the difference between the critical crossflow Reynolds number, at which the intermittency starts to increase, and the crossflow Reynolds number at transition. The T c1 criterion gives the value of Re δ 2 at the transition onset, when the velocity profile deviates from a laminar profile. Nevertheless, the transition model under study is constructed to activate the production of intermittency when the critical value is met, which is lower than the one encountered at transition onset.

The new F onset function in the transport equation for the intermittency γ, Eq.(3.53), becomes:

F onset = F length F onset + F length,CF F onset,CF , (6.51)
where F length = 100, F onset from Eq.(3.55), and F length,CF = 5. The latter parameter has been set to considering that crossflow instabilities develop on a longer length compared to T-S waves (F length = 100) as referring to the numerical calibration performed by DLR for the their crossflow criterion. The constants, a and F length,CF , and c were obtained from the computations performed around the 6:1 prolate spheroid at α = 15 • , 30 • angles of attack, on a mesh fine enough to consider negligible the discretization error.

Limitations of the T c1 recalibration

It is important to keep in mind that the presented criterion was calibrated under the stringent assumptions of the FSC equations, i.e. zero-pressure gradient in the spanwise direction and zero crossflow for zero pressure gradient flows. No further modification has been included to account for the fact that these hypothesis are no more valid for real three-dimensional boundary layers. The choice of the normal derivative of the normalized vorticity as indicator of the crossflow strength, the function Ψ, has been proven to give a fairly accurate approximation of the ratio X = Re δ2 /Re Vmax within the FSC framework. On the other hand, the pressure-gradient parameter λ θ approximation is a more critical function. First of all, the hypothesis that λ θ can be written as a function of the wall-normal derivative of the wall normal velocity comes from 2D boundary layers considerations and it does not represent the real physics. Secondly, the momentum thickness θ is replaced by the wall distance y, which provides at the middle of the boundary layer a proper scaling of the momentum thickness θ. For this latter reason, the function λ θ has been calibrated approximately at the middle of the boundary layer and then evaluated at the location where the local T c1 reaches its maximum. The position in the boundary layer where λ θ is calibrated is a crucial question. However, at this stage, we are not able to give an alternative approach to the one used for the calibration presented in this dissertation.

Concerning the numerical calibration for the definition of the constants appearing in the model formulation, it is important to underline that there is no universal rationale behind their setting. Different values may be chosen for crossflow transition modeling on other configurations, than those used for the present numerical calibration. For instance, through F length,CF , we are able to control the strength of the intermittency production term, pushing the transition front upward or downward. Nevertheless, the experimental data used for the calibration covered a swept wing flows and complex flows around three-dimensional configurations, as the prolate spheroid, and they might be sufficient for the constants calibration.

Galilean invariance

The idea behind this re-calibration of the T c1 proposed by Menter & Smirnov is to include more physics within the transition criterion. This is achieved by including a local approximation of the sweep angle in the model formulation. The approximation proposed for φ L is widely used, but the explicit use of the local velocity vector makes the model not Galilean invariant. In order to achieve a weak Galilean invariance of the new formulation, the local sweep angle is numerically implemented using the relative velocity vector -→ u rel defined as

- → u rel = - → u -- → u wall , (6.52)
instead of the local velocity vector. This is a fair modification, when dealing with boundary-layer transition. In ISIS-CFD, the velocity at the wall is known throughout the simulation. From the beginning of the simulation, for each cell center, close to a no-slip wall, the face index of its correspondent point at the wall is stocked in a table. By "its correspondent point at the wall" is meant the point which minimizes the distance from the cell center to the wall. As in ISIS-CFD, the search of the point at the wall is a feature common to most of industrial codes. This implementation has its own limitations: the research of the point at the wall might be troubling at the junction between multiple bodies with possibly different velocities. Indeed, there could be several points on the different surfaces that minimize the distance to the wall. This is the case of the rotor blades of an helicopter, which can undergo crossflow transition. However, a strong Galilean invariant formulation would require the discard of the sweep angle. Indeed, even if the local velocity vector is substituted by another variable, the dependence on the axes aligned with the velocity would be intrinsic and its use outside boundary layers can hardly be defended as Galilean invariant.

T-S/CF interaction In real flow physics there is a considerable interaction between streamwise and crossflow instabilities, that accelerates the transition process. For instance, Reed, [START_REF] Reed | An analysis of wave interactions in swept-wing flows[END_REF], showed that the interaction between CF and T-S waves produces a double exponential growth of T-S. In the presented γ model variant, the final F onset function results in the sum between the contribution coming from the streamwise transition criterion and the new crossflow criterion. However, the interaction between CF and T-S is not really accounted for in the model formulation. Indeed, F onset,1,CF , Eq.(6.46), does not activate at the same position as F onset,1,T S , Eq.(3.55): the T c1 criterion reaches its maximum close to the leading edge of the boundary layer, while the vorticity Reynolds number ratio for streamwise transition is constructed to activate at the middle of the boundary layer. A not weighted sum of the functions that trigger the production of intermittency is not enough to account for the physical interaction of the two mechanisms.

Choi et al., [START_REF] Choi | Recent improvement of a correlation-based transition model for simulating three-dimensional boundary layers[END_REF], proposed to use the local sweep angle as a measure of the T-S/CF interaction: higher is the local sweep angle, stronger is the interaction between the instabilities. In their model, presented as CF extension of γ -Re θ formulation, they propose the inclusion of a function F inter as a linear combination of their F onset,1,CF , based on the crossflow Reynolds number, and F onset,1,T S from the original γ -Re θ , Eq. (3.19): The constants 36 43 and 77 86 are obtained by the stability diagram for the 6:1 prolate spheroid shown in Fig. (2.12). In their formulation, the function F length , which is not a constant, but it depends on the momentum thickness Reynolds number at transition onset, Re θt . This approach has not yet been explored for the present implementation, but further studies will be conducted in the next future.

F onset = max(F onset
Influence of T u The receptivity stage is not considered in the crossflow criterion: neither the surface roughness, nor the turbulence intensity are accounted for in the criterion formulation, despite their importance in the determination of the strength of stationary crossflow waves. Nevertheless, the chosen value of inlet turbulence intensity will affect the transition front position, pushing it downstream for lower level of turbulence. This behavior represents the model reaction to the variation of turbulence intensity, which can enhance or lessen the turbulence production term, but it does not represent the true physics.

Implementation in a CFD solver Within ISIS-CFD the derivative of the vorticity that enters the T c1 criterion through the indicator Ψ is calculated as the Hessian of the velocity. Each entrance of the Hessian matrix is built through a least squares 3rd order accurate interpolation. The evaluation of the hessian of a scalar quantity in a the center of the cell C 0 makes use of n points that provide the centers of the neighboring cells C i=1,...,n . As a first set of neighbors cells (C 1 ) i the volumes which share a vertex with the cell C 0 are taken. As a second set of neighbors the volumes (C 2 ) i which share a face with (C 1 ) i are chosen. This approximation is fairly good for 2nd derivatives calculated at the edge of the boundary layer, which is our region of interest. The calculation of a 2nd derivative within the boundary layer might be troubling and requires a mesh enough refined in both streamwise and spanwise directions. Nevertheless, these stringent requirements on the mesh refinement are the same demanded by the transition models for accurate predictions with negligible discretization error.

The re-calibrated T c1 criterion for CF transition presented above has been included in the oneequation transition model γ. In the following, we are going to refer to this CF variant as the "γ + CF " model. The numerical performance of the latter will be discussed in Chapter 7, devoted to simulations for 3D boundary layers, around the 6:1 prolate spheroid and the sickle wing. Some of the results computed by γ + CF model will also be compared to γ -Re θ simulations. Within γ -Re θ formulation, CF transition is accounted for through an helicity-based crosswise criterion. The concept of helicity as alternative to the C1 criterion is presented in the next section.

The length scale l represents the momentum thickness at the point where the helicity Reynolds number reaches its maximum within the boundary layer. The constant C He,max = 0.6944 is calibrated for infinite swept wing flows. H 12 in Eq.( 6

F onset3,CF = max 1 - R T 2 3
, 0 (6.64)

F onset,CF = max(F onset2,CF -F onset3,CF , 0) (6.65)
where the constants C = 0.7 and F length,CF = 5 result from numerical calibration. The explicit use of the helicity makes the model not Galilean invariant, nevertheless it can be opted for the use of the relative velocity, as proposed for the recalibrated T c1 criterion.

Conclusions on the T c1 New Calibration

In this chapter, it has been proposed a new calibration of the crossflow criterion originally presented by Menter & Smirnov . The criterion is meant to predict crossflow transition in a low turbulence environment and it is based on Arnal's C1 criterion. Following Menter, the formulation of Arnal correlation is split into functions which account for the different parameters that affect the transition process. For each of this individual term, a local function related to the global information in Arnal correlation is built. The local version of the C1 criterion is referred to as the T c1. The directional change of the vorticity across the boundary layer, the quantity Ψ in Eq.(6.21), is used as indicator of the ratio of the crossflow to streamwise strength, Re δ 2 /Re V . The indicator Ψ is proportional to the local change of the flow angle and it can be interpreted as a measure of the crossflow strength. Other than the crossflow strength, we need to account for the pressure gradient in the streamwise direction effects. This is the role of the function G in Eq.(6.21). Menter & Smirnov propose as G a third-order polynomial in the one variable λ θ,CF , a local approximation of the pressure gradient parameter. The latter one is calculated as a linear combination of the normal derivative of the normal velocity components, dv dy . However, their calibration does not take into account the effect of the sweep angle on the pressure gradient parameter. For this reason, we propose a new calibration that envisages the use of the sweep angle within the formulation of the T c1 criterion. The new function G is constructed as a nonlinear polynomial in two variables. The first one is a local reconstruction of the sweep angle, φ L . This is obtained using the pressure gradient and the local velocity vector. The explicit use of the velocity makes the formulation non Galilean invariant. Nevertheless, a strategy has been proposed to restore a weak Galilean invariance in the boundary layer. The second variable is the pressure gradient parameter λ CF , approximated using the wall-normal derivative of the wall-normal velocity dv dy and, again, φ L .

6:1 Prolate Spheroid

The 6:1 prolate spheroid is the configuration on which we are going to focus our attention the most. This geometry is commonly used both in hydrodynamic and aerodynamic, because it serves as simple surrogate for axisymmetric bodies as airplane fuselage and submarine hulls, as well as engine cowling on helicopters. The 6:1 prolate spheroid at incidence is one of the most investigated test case, both experimentally and numerically, because it exhibits all the complex physics associated with cross-flow separation. This separation, due to the streamlines curvature, generates pairs of counter rotating vortices that are further convected downstream. The complexity of crossflow separation stands in its not unique character. It does not originate from an unique singularity point nor line. In addition, being the flow nearly attached to the surface ahead separation, it does not strongly interact with the local flow field. Crossflow separation can be identified as a boundary layer rolling around itself. We refer to crossflow separation as open separation, characterized by an envelope of converging streamlines, following the definition by Surana, [START_REF] Surana | Exact theory of three-dimensional flow separation. Part 1. Steady separation[END_REF], and used also by Wetzel et al., [START_REF] Wetzel | Measurement of three-dimensional crossflow separation[END_REF]. An extensive review on the 6:1 prolate spheroid can be found in Andersson et al., [START_REF] Andersson | Instabilities in the wake of an inclined prolate spheroid[END_REF], and Fu, [START_REF] Fu | Challenges in naval ship and submarine hydrodynamics[END_REF]. Experimental measurements can be found in Fu et al., [START_REF] Fu | The flow structure in the lee of an inclined 6: 1 prolate spheroid[END_REF], and Wetzel et al., [START_REF] Wetzel | Unsteady crossflow separation location measurements on a maneuvering 6: 1 prolate spheroid[END_REF]. These experiments, performed in the 1990's were used by different groups for the assessment of (U)RANS, DES and LES, Kim et al., [START_REF] Kim | Application of modern turbulence models to vortical flow around a prolate spheroid[END_REF], Wikström et al., [START_REF] Wikström | Large eddy simulation of the flow around an inclined prolate spheroid[END_REF], and Fureby et al., [START_REF] Fureby | LES of the flow past a 6: 1 prolate spheroid[END_REF], [START_REF] Fureby | RANS, DES and LES of the flow past the 6: 1 prolate spheroid at 10 and 20 angle of incidence[END_REF]. In the present context, however, we are interested in the characterization of the flow at the surface and the prediction of laminar-to-turbulence transition within the boundary layer. This is often a flow feature that is not discussed. In fact, being the focus concentrated more on separation than flow transition, the boundary layer is often tripped. Turbulence is forced close to the leading edge, in order to avoid transitional effects. When the geometry is not tripped, the majority of the experimental database is concentrated in the vicinity of the separation location and in the wake. So far, the only experimental database that fully describes transition at the surface, through the measurements of wall shear stress magnitude and direction, was provided by Kreplin in 1985. These are the measurements we are going to refer to in the following discussion.

Experimental Set Up

The measurements performed by Kreplin, [START_REF] Kreplin | Wall shear stress measurements on an inclined prolate spheroid in the DFVLR 3m × 3m low speed wind tunnel[END_REF], in 1985, were conducted in the 3m × 3m low speed wind tunnel at DLR Gottingen. The length of the specimen was 2.4m and different Reynolds numbers, ranging from Re = 1.5 × 10 6 up to Re = 7.2 × 10 6 , and angles of attack, from α = 5 • to α = 30 • , were tested. The equation of the 6:1 prolate spheroid surface is given by: The magnitude of the wall shear stress τ w is derived from the sum of the heat transfer rates of the films of each probes. A rough estimation of the error bound for the wall shear stress magnitude is given to be ∆τ w = ±20%. This high error bound is related to the fact that the hot-film probes were calibrated on flat tunnel wall for a 2D turbulent boundary layer, and, around this geometry at incidence, the laminar region is quite extended. Although the probes calibration raises some doubts on the reliability of the values of C f predicted around the 6:1 prolate spheroid, we are confident in the accuracy of the measured transition location.

x 2 + 36y 2 + 36z 2 = 1 4 . ( 7
For high Reynolds number, such as the one under consideration Re = 6.5 × 10 6 , transition occurs under the interaction of Tollmien-Schlichting and crossflow instabilities. For the angles of attack α = 15 • , 30 • , zones of pure crossflow in the middle of the inclined prolate spheroid are observed and they become wider as the inclination is increased. At low Reynolds numbers, such as Re = 1.5 × 10 6 , the main transition mechanism are T-S instabilities, as reported in Krimmelbein, [START_REF] Krimmelbein | Transition prediction for three-dimensional flows using parallel computation[END_REF]. For this reason, they were not considered in the present computations.

Unfortunately, no specific indication on the free-stream conditions is given in the experimental report, but it is mentioned that the turbulence intensity T u varies from 0.1% to 0.3% depending on the Reynolds number. As already mentioned for 2D test cases, the unknown exact initial conditions can be very problematic for transition predictions, being the transition models strongly dependent on the inlet turbulence values.

Computational Domain

The grids were provided by Rui Lopes, from IST Lisbon. A set of 5 multiblock structured grids was generated with the GridPro software, with an O-topology encircling the spheroid. The computational domain is a box of total length 200L and width 100L, where L is the length of the 6:1 prolate spheroid. The geometrical center of the body is positioned at x/L = 0 and its distance from the boundaries is approximately 100L. The incidence angles are imposed by rotating the spheroid with respect to its center, as well as an inner O-block around it, and the flow is aligned with x-axis. The grids are for half of the geometry, making use of a symmetry plane. The finest grid has 760 cells in the longitudinal direction, N x , measured along the upper side of the surface, and 176 cells in the transversal direction, N φ , measured along the plane located at half of the longitudinal length of the surface. The size of the first near wall cell in the direction normal to the surface is ∆y ∼ 2.3 × 10 -6 . The finest grid has a total of 126016 cells on the surface of the spheroid and 42.5M volume cells. The remaining four grids are obtained from the finest one using coarsening factors of 0.875, 0.75, 0.625 and 0.5. This means that the number of cells along each direction and for each grid can be found by multiplying the values for the finest grid by the coarsening factor. Mesh details are given in Table(7 The inlet conditions for all the angles of attack under study are reported in Table (7.3): the turbulence intensity, the eddy viscosity ratio at the inlet, and the value of turbulence intensity in the vicinity of the body. These inlet conditions were chosen within the AVT-313 group.

Re

T u in (%) ν t in /ν T u(%) 6.5 × 10 6 0.5 250 0.15 γ mesh dependence

In a first attempt, a slightly different set of meshes was put at disposal for running the transition simulation. The details of the coarsest set of this mesh, named Grid5bis, are given in Table (7.4). The mesh Grid5bis is part of a mesh set coarser than the one used for the presented simulations, whose details were reported in Table (7.2). In particular, if we compared Grid5bis and Grid5, we will notice that even if the first one has more cells in the transversal direction, N φ = 88 on Grid5, the cells in the longitudinal direction are considerably less, N x = 380 on Grid5. Both grids have similar size of the first near wall cell in the direction normal to the surface, which is about ∆y ∼ 4.6 × 10 -6 . The C f contours obtained on both grids using γ without any crossflow inclusion are shown Fig. (7.5) for α = 15 • . x/L represents the axial position along the surface, while φ is the azimuthal angle. φ = 0 • corresponds to the symmetry plane on the windward side, while φ = 180 • is the symmetry plane on the leeward side. On Grid5bis, we expect the discretization error to pollute the solution, nevertheless the results are quite surprising. γ predicts two different flow physics: the flow on Grid5bis is turbulent, while on Grid5 transition due to TS is predicted where expected. The intermittency variable γ is imposed to be equal to 1 everywhere except at the wall, in order not to interfere with the turbulence model in the free-stream. The transition models work first by installing turbulence everywhere, and then destroying it throughout the simulations resolving the transport equation for γ. On a mesh too coarse in the longitudinal direction, the intermittency γ is over-resolved, because of the significant numerical error, causing the flow to be practically turbulent. These results show the dependence of the transition model γ on the mesh density. Notably, the number of cells in the streamwise direction plays an important role. A sufficient number of cells from the start of the laminar region up to the point where transition is completed needs to be provided, in order to assure an accurate solution in the laminar region. Best practices for meshing clearly differ from case to case and depend on the physics. After this parenthesis of the mesh dependence of γ on the numerical results, let us discuss the models performance for the different angles of incidence. Computational Costs All the computations around the 6:1 prolate spheroid were run on national HPC resources. For each grid, the interpolated solution from its correspondent "one level coarser" grid is used as initial solution. The computations by both γ +CF and γ -Re θ +CF on the coarsest grid, Grid5, were run on 280 processors. The simulation for γ-Re θ +CF required ∼ 13000 CPU hours, while for γ + CF , ∼ 9000 CPU hours. Simulations on finer meshes were run on an higher number of processors, up to a maximum of 784.

α = 15 •

Let us start by the angle of incidence α = 15 • . Results for γ + CF and γ -Re θ + CF were computed on all the five grids. The results of the grid convergence of the friction and pressure drag coefficients are shown in Fig. (7.6) and (7.7). We observe that we obtain an estimated order p = 2 on the friction component, while the pressure drag presents a quasi-2nd order of convergence with respect to the grid refinement. The behavior of the pressure is related to the noisy convergence of the non-linear residuals on such fine meshes, as discussed for 2D cases. We also notice that C d,f increases with the grid refinement for γ -Re θ + CF results, in contrast to γ + CF results, for which C d,f decreases. This is related to the behavior discussed above, namely the overestimation of the intermittency on coarser grids. These results show the stronger dependency of γ on the spatial discretization compared to γ -Re θ model. The same observation is valid in the two-dimensional case. Using γ +CF the predictions are very similar in terms of transition front predictions, other than forces, on the three finest grids, Grid1, Grid2, Grid3 from Table (7 ). This is a mesh effect, rather than an erroneous prediction. Indeed, some crossflow occurs close to the vertical symmetry plane, but, because of the symmetry plane, the streamlines are forced to follow the gridlines, delaying the transition process. This laminar tongue is grid dependent, i.e. it gets shorter on coarser grids. Over all, the predicted C f is noticeably lower than what predicted in the experiments. Only from the central section of the spheroid, the quantitative deviation in the skin friction magnitude between experiments and predictions is justified by the experimental uncertainty. This quantitative difference is most probably related to the input uncertainty. As already mentioned the free-stream T u is not exactly known. A different turbulence intensity in the vicinity of the prolate spheroid nose is expected to increase or decrease the C f values at the surface. Most probably, the T u used in the simulations is lower than the experimental. We want to dedicate this section to a summing-up of the conclusions on the numerical results obtained around the 6:1 prolate spheroid at the angles of attack α = 5 • and α = 15 • . Indeed, as mentioned in the beginning, these results were presented in occasion of the AVT-313 workshop, held virtually in March 2021. All the participants to the workshop have presented their own simulations, computed using different solvers and different crossflow variants.

Without entering the merits of the used computational solvers, the other models used for this exercise by the different research groups were:

• both γ-Re θ and γ extended to CF inclusion using the helicity-based crossflow criterion, as presented in Sec.(6.5);

• γ-Re θ , with the crossflow variant proposed by Langtry in 2015. As mentioned in the introduction to Chapter 6, this crossflow inclusion plays on the reduction of the destruction term of the transported momentum thickness Re θt . This strategy relies on a helicity-based crossflow criterion. This variant is the only in literature which accounts for surface roughness h rms , that plays a major role in the destabilization of stationary CF. For instance, as shown in Crouch and Ng, [START_REF] Crouch | Variable n-factor method for transition prediction in threedimensional boundary layers[END_REF], the N-factor for stationary CF decreases as the roughness level increases. Around the NLF(2)-0415 airfoil, Radeztsky et al., [START_REF] Radeztsky | Effect of isolated micron-sized roughness on transition in swept-wing flows[END_REF], showed that this dependence is quasi-linear. In Kreplin report, the surface roughness is not mentioned. For the simulations, h rms is assumed to be 3.3 µm, which is typical for painted surfaces.

• e N two N-factor approach, based on the stability diagram in Fig. (2.12). This N T S -N CF diagram was explicitly built for the experiences run in the wind-tunnel in Gottingen around the 6:1 prolate spheroid.

Solutions were computed on different grids, both structured and unstructured. Different inlet boundary conditions combinations have been used in order to obtain T u ∼ 0.15% at the nose of the prolate spheroid.

The results show a general trend. All the transition models, combined with their respective crossflow criteria, fail to predict transition on the windward side, for approximately φ < 30 • . The e N method is the only strategy able to account for transition in this region. This behavior on the windward side shed light on the deficiencies of the current correlation-based transition model. The models do not account for any transition mechanisms other than crossflow and T-S waves, nor for the interaction between different transition mechanisms. Around the 6:1 prolate spheroid, the possible transition mechanisms, missed in the correlationbased transition models, are:

• the non-linear interaction effects between T-S waves and crossflow instabilities;

• attachment line instabilities;

• traveling crossflow;

• Görtler vortices.

Close to the windward symmetry plane, no-crossflow transition can happen, because the crossflow component of the velocity is zero. Let us discuss separately each mechanism and analyze if their occurrence around the 6:1 prolate spheroid is a real possibility.

The non-linear interaction between T-S waves and crossflow instabilities is surely the main feature which is not accounted for in the models. On the 6:1 prolate spheroid, for α = 5 • , 30 • , important interaction between T-S and CF is expected, especially for the lowest angle of incidence. Among the transition modeling approaches used in the workshop, the only one able to capture the two mechanisms interactions is the e N method. Indeed, it accounts for the lower critical N T S and N CF amplification factors, when there is some sort of interaction between the two instabilities. As mentioned at the end of Chapter 6, an effort to include the interaction between the CF and T-S within γ -Re θ has been made by Choi et al., [START_REF] Choi | Recent improvement of a correlation-based transition model for simulating three-dimensional boundary layers[END_REF]. In their simulations, the crossflow criterion is based as well on the C1 from Arnal transition, using different quantities than those proposed in this dissertation. In their variant, the waves interaction is accounted for through the introduction of the weighted sum of the transition onset functions, the one for T-S and the one for CF. The F length functions is reformulated as well, depending on the transition mechanism. As shown in their publication, with this formulation a better agreement with experiments is achieved on the windward side of the spheroid with a transition front pushed closet to the windward symmetry line. The angles of incidence tested in their publication are α = 10 • and α = 20 • , but we can suppose that the same results can be achieved for α = 5 • and α = 15 • .

Leading edge instabilities are related to 2-D boundary layer developing at the attachmentline. They can cause the flow to become turbulent along the streamwise direction, as mentioned in Chapter 2. Attachment line instabilities represent a very possible mechanism that can occur at the windward symmetry plane. Indeed, crossflow can move the transition line forward, but not all the way down to the attachment line.

As mentioned in Chapter 2, the key parameter for the attachment line instabilities is R, computed with respect to the crossflow velocity component at the edge of the boundary layer. Unfortunately, the crossflow velocity is not a quantity commonly available in a CFD solver, and computations cannot be helpful to check if the stability limit R = 583 is exceeded or not. Nevertheless, some theoretical hypothesis can be made. If we approximate the spheroid at the location of its maximum width by a swept cylinder, R can be written2 with respect to the Reynolds number based on the radius r and the leading edge sweep angle Φ, defined as Φ = π/2 -α, where α is the angle of attack. The final formulation reads as

R = Q ∞ r sin(Φ) tan(Φ) 2ν . (7.2)
The radius r is L/12, where L is the reference length of the spheroid, and Q ∞ is the modulus of the velocity in the free-stream. Making the calculations, we will have R bigger than the stability limit for both angles of attack: R ∼ 1800 for α = 5 • and R ∼ 1000 for α = 15 • . Clearly, this is just an approximation, and R around the 6:1 prolate spheroid is lower, considering that it has more streamline divergence than a swept cylinder in its first half. Nevertheless, this calculation shows that attachment line instabilities might be a possibility.

Undoubtedly, stability-based method as e N are able to predict physical mechanisms that the correlation-based transition models under study miss. The reason is that only using linear stability we account for the physical direction of propagation and amplification of the instabilities. Let us consider the streamwise transition criterion based on Re θ . As discussed by Arnal, in [START_REF]Three-dimensional boundary layer: laminar-turbulent transition[END_REF], transition along the attachment line is the typical case where transition is entirely governed by streamwise instabilities. Nevertheless, the streamwise criterion based on the characteristic boundary layer thickness cannot predict it. Quoting Arnal, in two dimensional flows, this criterion gives good results because an implicit relation exists between Re θ and the streamwise direction. However, this relation fails on complex three dimensional configurations where the streamlines are not parallel. For instance, laminar boundary layer calculations performed by DLR3 showed that the characteristic boundary thickness remains constant along the plane φ = 0 • , which is the windward symmetry plane of the 6:1 prolate spheroid. The constant boundary layer thickness prevents the streamwise criterion based on Re θ from predicting transition in this region.

Traveling crossflow Although e N predictions assumes stationary crossflow instability to be the dominant player in the transition process in the experiment by Kreplin, transition can also occur due to traveling crossflow instability. For the range of free-stream turbulence intensity, T u = 0.15% -0.35%, as indicated in the experimental report, there is reason to believe that the traveling crossflow instabilities can play an important role in initiating the breakdown process. For instance, in the classic experiments by Deyhle & Bippes, [START_REF] Deyhle | Disturbance growth in an unstable three-dimensional boundary layer and its dependence on environmental conditions[END_REF], at DLR, the traveling crossflow instabilities controlled the transition behavior for T u > 0.15%, which is the level of freestream turbulence used for the presented computations. Even if free-stream turbulence influences the transition process indirectly, compared to the two-dimensional case, for a turbulence level T u > 0.2% and smooth surfaces the traveling instability waves dominate. Nevertheless, all the crossflow transition criteria discussed in this chapter ignore the role of traveling crossflow waves, because calibrated on experiences run at very low free-stream turbulence levels.

Görtler vortices are the least possible transition mechanism that can occur on the windward side of the spheroid, being a prerogative of concave surface. Nevertheless, they can originate at the leading edge of a blunt body, because of the local concavity of the external streamlines at the leading edge. These vortices appear in the laminar steady flow and their spanwise modulation can further destabilize T-S waves. Footprints of Görtler vortices were also observed experimentally close to the leading edge of the NACA0015 profile, as discussed in Miozzi et al., [START_REF] Miozzi | Skin friction and coherent structures within a laminar separation bubble[END_REF]. Nevertheless, this scenario is likely in a high level free-stream turbulence environment, for which the presence of streamwise vorticity within the boundary layer can be justified. These conclusions on the numerical results for α = 5 • and α = 15 • are more transition-physics related, rather than directed to a transition modeling discussion. Nevertheless, establishing the causative factors for transition is a fundamental pre-requisite to further move the discussion on the modeling capabilities of the different approaches.

Clearly, the presented analysis shows one of the main drawback of the correlation based transition models: they miss three-dimensional physical processes, that can be only predicted using linear stability theory.

α = 30 •

The 6:1 prolate spheroid series is concluded with the angle of incidence α = 30 • . At this angle of attack, crossflow instabilities dominate the flow transition on a large portion of the prolate spheroid. As mentioned in the introduction to this chapter, for this test case, we are going to discuss only the numerical predictions by γ, using the re-calibrated T c1 criterion presented in this dissertation, as inclusion of crossflow transition, referred to as γ + CF . The experimental C f contours are shown in Fig. (7.22). Taking as a reference the middle section of the 6:1 prolate spheroid, the one of maximum width, the flow starts transitioning on the windward side at φ ∼ 50 • and transition completes at φ ∼ 70 • . On the leeward side, the situation appears slightly more complicated. A laminar separation bubble occurs near the nose of the spheroid, the flow separates and then rapidly reattaches. This separation-induced transition promotes the appearance of T-S waves on this side of the spheroid. The flow undergoes a second separation further downstream. The large skin friction values on the leeward side result from the induced velocities due to the separated vortex flow from the windward side, as mentioned in Kreplin The transition model fails, because Re V in the first 20% of the spheroid length does not exceed the limit value of 2200, Eq.(3.79), to activate the additional production term of turbulence kinetic energy. The function that accounts for separation (Eq.(3.78)) never activates. This behavior might be related to the diverging streamlines predicted at the leeward symmetry plane, close to the leading edge, Fig. (7.25). The streamwise criterion fails because of the differences between the physical distance along which the waves propagate and the distance along which the characteristic boundary layer thickness grows, as mentioned in Sec.(7.2.3). The unpredicted laminar separation bubble pollutes further downstream the numerical solution. The flow at the leeward side numerically transitions around X/L = 0.139, but the skin friction is systematically underestimated with respect to measurements, as shown in Fig. ). The separation line for the flow coming from the leeward side, which is the attracting portion of the skin friction line, can be roughly identified with respect to the local minimum of the magnitude of the wall shear stress. This is now a fair approximation because the skin friction lines are converging from different sides of the prolate spheroid. This second separation location occurs at the leeward side approximately around X/L ∼ 0.55, but predicted skin friction lines are converging further downstream. From the C f contours the predicted minimum is located at X/L = 0.9. The separation might be pushed downward because of the underpredicted laminar-to-turbulence transition within the boundary layer on the leeward side. The local maximum peak of C f at X/L = 0.565 at φ ∼ 130 • , Fig. (7.24d), due to the velocity fluctuations induced by the separated vortex flow, is consider-ably pushed downward close to the trailing edge. Due to the large separation the flow undergoes at such an high angle of attack, we cannot exclude that the RANS turbulence model does not affect the overall results at separation. Turbulence models typically fails in massively separated region. To this matter, efforts have been going on to combine transitional RANS model into an hybrid RANS-LES method, Hodara, [START_REF] Hodara | Hybrid RANS-LES closure for separated flows in the transitional regime[END_REF], and Kim and Carrica (private conversation). As for turbulence models, the key issues for hybrid RANS/LES models is the location of the RANS to LES switch. A switch to LES within the boundary layer or close to its edge can strongly affect the numerical results leading to nonphysical results at the surface. The boundary layer should be shielded in order to avoid the contamination of the transition process because of the change of turbulence quantities in LES region, i.e. decay of the modeled turbulence kinetic energy within the grey region.

Sickle Wing

The last validation case is the sickle wing configuration. This wing design is conceived in order to generate a 3D boundary layer with increasing crossflow in the spanwise direction. Crossflow modes are highly amplified by the large spanwise gradients created in correspondence of the sweep kinks of the sickle shaped planform and move the transition front upward. Our interest toward this geometry relies on the will to assess how much these strong spanwise gradients challenge the FSC assumptions, therefore the T c1 crossflow criterion. The flow around the sickle wing is computed only using γ model with the recalibrated T c1 criterion for crossflow inclusion, γ + CF variant. Experimental data for different Reynolds numbers and angles of attack are available in Petzold et al., [START_REF] Petzold | Transition on a wing with spanwise varying crossflow and linear stability analysis[END_REF], and Kruse et al., [START_REF] Kruse | Transition prediction results for sickle wing and NLF (1)-0416 test cases[END_REF]. Both experimental campaigns were run in an atmospheric, closed circuit type, low-speed wind tunnel: the DNW-NWB for the measurements of Kruse

Computation Domain

Results for this configuration were computed on an hexahedral unstructured mesh, generated using the Hexpress grid generation package. In the attempt to reproduce the wind tunnel dimension from Kruse et al., the inlet is located at x/C = -4.6 and the outlet at x/C = 8.455. The top and bottom walls are located at z/C = ±1.633 and the side wall at y/C = 2.8. In order to obtain the best agreement with the experimental results, the tunnel walls are treated as slip boundary conditions. No slip conditions are imposed on the body. The body is turned with respect to the angle of incidence under study and the flow is aligned with x-axis. The grid was generated performing a systematic grid refinement, a simplified variant of the anisotropic adaptive gird refinement implemented in ISIS-CFD solver. Starting from an initial mesh of 10M cells, with approximately 50 cells in the boundary layer, the mesh has been further refined in the longitudinal. In order to avoid an excessive refinement of the cells at the edges of the sickle wing, a maximum cell size equal to 1/4th of the cells size in the x-direction has been imposed. The refinement procedure has been limited to the region next to the body, in order to avoid unnecessary cells in the free-stream region. Different grid refinement levels have been tested, nevertheless, the one described above is relative to the grid used for the computations presented in this chapter. The computational domain is presented in Fig. (7.27), as well as the mesh at the surface. Details on the mesh are reported in Table (7.5).

Computational Costs

The simulations for the sickle wing were run as well on HPC national resources. The simulation on the grid, whose details were given in Table (7.5), was performed on 784 processors for approximately ∼ 27440 CPU hours. The numerical solution by γ without crossflow inclusion has been used as initial solution for the computation by the crossflow variant of γ, γ + CF .

Results

The skin friction contours on the upper and lower surface of the sickle wing are shown in Fig. (7.28). Predictions by γ without and with crossflow inclusion, γ and γ + CF variants, respectively, are compared to the extracted experimental transition lines. As a general remark, all the γ variants predict a jagged transition front. This is not due to the unsteadiness of the flow, but rather a consequence of the unstructured grid topology and the hanging nodes in the grid close to the boundary layer. The differences between predictions by γ without and with crossflow inclusion are noticeable, especially on the upper surface, where transition is dominated by crossflow instabilities, Fig. (7.28a) vs Fig. (7.28c). γ predictions are in good agreement with measurements on Sec. A. The transition model without crossflow inclusion is able to reproduce the laminar separation bubble that occurs towards the trailing edge. Moving upward, towards Sec. B, the measured transition line considerably changes in the spanwise direction. The transition process is here CF-dominated. As mentioned in Kruse et al., [START_REF] Kruse | Transition prediction results for sickle wing and NLF (1)-0416 test cases[END_REF], the sudden shift upward of the transition front in the middle of Section B is related to disturbances caused by the discrete roughness of the pressure tap row. This shift is not observed in the experiments from Petzold et al., where the transition line is uniform and upwards compared to measurements from Kruse et al.. The numerical simulation does not account for the discrete roughness effects, and predictions by γ + CF model on the upper surface are closer to the measurements from Petzold et al.. We believe that differences between Petzold et al. and Kruse et al. results stand in the experimental turbulence intensity in the free-stream and the different strategy used in the identification of the transition location. On Section C, crossflow instabilities become weaker, as indicated by the downward shift on the transition front. This weakening is fairly represented in the numerical results by γ + CF . The features we are interested to are the turbulent wedges in the regions of sweep changeover, from one section of the sickle wing to the other. As mentioned before, these are related to the strong spanwise gradients. The recalibrated T c1 criterion within γ + CF model predicts fairly well these turbulent wedges, in spite of its calibration on FSC solutions, for which zero-spanwise gradients are assumed. The agreement between the experiments and simulations with crossflow inclusion is quite satisfactory, nevertheless the effect of the discretization error on the solution is quite evident. In particular at the junction between Section A and Section B, where the spanwise gradients are considerably strong. The turbulent wedges predicted in the computations are underestimated and diffused on their sides. Despite the overall high number of cells of the mesh used for the computations, the grid at the surface is too coarse in correspondence of the kinks to predict accurately such strong gradients, see Fig. (7.27b). As an example of a proper mesh that should be used for the current computation, we show the grid used by Kim et al., [START_REF] Kim | Boundary layer transition models for naval applications: capabilities and limitations[END_REF], for their transition simulations presented in the paper, see Fig. (7.29). In order to predict accurately the gradients in the spanwise direction, a cluster of cells in correspondence of the junctions is needed. Such grid cannot be generated in Hexpress. Further refinement is needed using, for instance, either systematic or adaptive grid refinement techniques. The mesh should be refined in several steps in order to obtain a mesh of the same kind of the one presented in Fig. (7.29). Limiting boxes should be defined in order to avoid additional refinement at the leading and trailing edges. The main limitation of the proposed strategy is the random access memory (RAM) needed for running the computations, post-process and visualize the results. In fact, in a first attempt, this strategy has been tested on an initial mesh of approximately 3M cells. In order to obtain a sufficient mesh density in correspondence of the junctions, the mesh has been refined in three different steps. The final grid resulted in 1000M cells. The results on this mesh would require more than 1TB of RAM to be post-processed and this is unfeasible for the computing resources that we currently dispose of.

Conclusion 3D Test Cases

This chapter was devoted to the analysis of the performance of the correlation-based transition model γ in 3D boundary layers. The recalibrated T c1 crossflow criterion was proven to perform fairly well on the two geometries tested, the 6:1 prolate spheroid and the sickle wing. Compared to the original T c1 version from Menter & Smirnov, γ + T c1-MS, the γ + CF with the recalibrated T c1 criterion performs better. Their behavior was compared on the 6:1 prolate spheroid at α = 15 • . γ + T c1-MS, as implemented in ISIS-CFD, fails to predict transition in the middle/windward side of the spheroid, where the transition process is dominated by stationary crossflow instabilities. Within the recalibrated T c1, we account for the sweep angle impact on the pressure gradient in the streamwise direction and on the overall crossflow transition process. Including additional physical parameters within the correlations, we can improve considerably their reliability. The calibration of the T c1 in respect of the solutions of the Falkner-Skan-Cooke equations, which assume zero-spanwise gradients, does not strongly affect the results. The results for the flow around the sickle wing, that experiences strong spanwise gradients at the junctions between the segments, are quite encouraging. With the recalibrated T c1, as crossflow variant, we are able to assure satisfying predictions of the transition front on the upper surface of the wing, despite the discretization error pollution. The mesh used for this simulation is too coarse at the junctions and the intermittency is considerably diffused. At the end, predictions on the 6:1 prolate spheroid at different incidences and on the sickle wing are in good agreement with experiments and these results pave the way to the future employment of T c1 as crossflow criterion in transition modeling. The test cases presented in this chapter covered two very different geometries. Nevertheless, additional numerical simulations are needed to fully validate the presented T c1 crossflow criterion. We cannot exclude that the set of numerically calibrated constant, i.e. c, F length , a, see Eq.(6.50), might need some adjustments on different geometries. Nonetheless, the presented calibration is a very promising starting point and the number of constants is limited. γ + CF was further compared to γ -Re θ + CF model. As crossflow criterion for γ -Re θ it was opted for the helicity criterion from Grabe et al.. The performance of the two crossflow combinations is very similar, in spite of few differences in the predictions of the transition process. Neither γ + CF nor γ -Re θ + CF predict transition down to the windward symmetry plane of the 6:1 prolate spheroid for α = 5 • , 15 • . The reason is that none of the two models variant accounts for other transition mechanisms than crossflow and T-S/streamwise instabilities, nor for their non-linear interactions. The models fail because of missing physics in their formulation. The possible missing transition mechanisms around this geometry are, other than the non-linear interactions between CF and T-S waves, attachment line instabilities and traveling crossflow. On one hand, the correlation-based transition models fail to predict the streamwise instabilities at the attachment line, because there is no relation between Re θ and the direction of propagation of the disturbances, as mentioned by Arnal. On complex three-dimensional geometries, the streamlines are far from being parallel. In this case, the only option is to solve the stability equations and compute the N factor by integrating the local amplification rates along the external streamlines, quoting Arnal, [START_REF]Three-dimensional boundary layer: laminar-turbulent transition[END_REF]. On the other hand, the two crossflow correlations, the T c1 and the helicity-based, neglect traveling crossflow waves, because empirically calibrated on experiences run in low/ultra-low free-stream turbulence environments. At the end, regardless of the actual cause of transition in the experiments, it is our belief that including additional correlations to account for different transition mechanisms is crucial to pursue the validation of these RANS transition models. To this purpose, it is fundamental for the measurements to inform us on the instability mechanism(s) causing transition in order to predict the transition behavior and build empirical correlation to be further included in the transition model's formulations.

Chapter 8

Conclusions

This thesis presented a thorough assessment of the RANS local correlation-based transition models (LCTM) γ and γ -Re θ . These models are complementary to the turbulence ones and considerably enhance their performance for low/moderate Reynolds numbers. In fact, turbulence models fail to predict transition effects and lead to very inaccurate flow predictions. In the last twenty years, different research groups have devoted their attention to laminar-to-turbulence transition modeling. They proposed different models and strategies, but, in our opinion, the LCTM concepts is the most appealing one. The strength of these models stands in their local formulation and their generic character, that allows the inclusion of different transition mechanisms without considerably changing the structure of the solved transport equations. γ and γ -Re θ models coupled to k -ω SST (2003), for turbulence resolution, have been tested on a certain number of canonical test cases, which represent a perfect blend of simplicity and complexity. The choice of these configurations was also dictated by the availability of experimental measurements. Indeed, we are still in a validation stage of these models, and, in order to assess their performance and reliability, experiments are needed to be compared with. As two-dimensional configurations, we have analyzed the results around the ERCOFTAC flat plate for both natural and bypass transition. The flow around two different airfoils at incidence, Eppler387 and NACA0015, has also been considered. Approaches that account for transition are essential and this is a major capability for models, especially if we are moving towards a future in which simulations are supposed to substitute experiments. Clearly, the two transition models perform differently with respect to their different calibration. Nevertheless, the impact on the overall results is small. In fact, the purpose of a RANS transition model is not to describe the correct physical process, but rather to identify the critical points. For instance, in the case of a laminar separation bubble, we are interested in separation, transition and, eventually, reattachment locations, and both models were proven to perform fairly well. γ and γ -Re θ models are very promising, but their modeling capabilities comes at a cost. They cannot be considered entirely robust, because they require a considerably high number of nonlinear iterations up to convergence. This behavior is due to the use of very fine meshes, in both streamwise and wall-normal directions, and the presence of min and max limiters in the equations that are not smooth functions. In addition, steady assumptions of clearly unsteady phenomena can negatively affect the iterative convergence, which is very noisy. Numerical strategies can be opted for to reduce the number of overall non-linear iterations and improve convergence. Other possible solutions are the smoothing of the functions that cause the residuals to stagnate, being aware that this might require a recalibration of the correlations, and move from a segregated solver to a coupled one for the equations resolutions. Another option would be the use of an adaptive grid refinement approach to locally refine the meshes where spatial resolution is required. The other big limitation of the transition models is their dependency on the turbulence boundary conditions. The dependence of the transition predictions on the boundary condition is desirable, because it reflects the response of the flow to the environment. Nevertheless, in order to reproduce the experimental results, at least two quantities are needed to describe the turbulence in the free-stream. This is never the case, and the experimental information are often reduced to the average free-stream turbulence intensity, if any data is available. This information is not sufficient to describe the free decay of turbulence intensity ahead of the body. Numerically, this is predicted by the underlying turbulence model and depends on the values of inlet eddy viscosity R T . Higher is the imposed value of R T , more moderate is the free-decay of turbulence. However, such high values of eddy viscosity ratio can pollute the solution in the laminar region, the velocity profiles and the skin friction distribution. On the other hand, for very low R T , the decay of T u can be so considerable that its value at the leading edge is so small to not affect anymore the transition process. So far, there is no systematic solution to this problem. Few strategies to set the inlet conditions were discussed, but they are not practical for real multi-bodies geometries. We are not at the stage where we can assert with certitude which transition model between γ and γ -Re θ performs better. Nevertheless, within the AVT-313 group, we have remarked that γ model is more consistent among the different solvers. The numerical results obtained with this model are not strongly affected by its implementation in different CFD environments. In addition, its simpler mathematical formulation positively affects the iterative convergence and the overall number of non-linear iterations. Nevertheless, γ has to be used with caution, especially when it comes to mesh density. Contrary to γ -Re θ , indeed, it is highly sensitive to the spatial discretization. Because of the higher robustness of γ model, we have focused the analysis for three-dimensional flows on the performance of its crossflow variant.

For three-dimensional configurations, we have discussed the numerical results around the 6:1 prolate spheroid at three different angles of attack, and the sickle wing. As originally published, γ and γ -Re θ did not account for stationary crossflow transition, the most common transition mechanism in three-dimensional boundary layers in a low free-stream turbulence environment. While for γ -Re θ different crossflow transition criteria have been proposed by different groups, only one was found in literature for γ, the so-called T c1 criterion, proposed by Menter & Smirnov in 2014. The T c1 is a local version of the C1 crossflow criterion proposed by Arnal, reconstructed using the Falkner-Skan-Cooke solutions. It uses the wall-normal directional change of the vorticity vector as indicator of crossflow strength. As published by Menter & Smirnov, the criterion does not work accurately on the 6:1 prolate spheroid, the main configuration we focused our attention on. In this dissertation, we have proposed a new recalibrated version of the T c1, which uses a local approximation of the sweep angle φ, computed with respect to the pressure gradient and the local velocity vector. This modification accounts for the sweep angle effect on the pressure gradient parameter and on the overall transition process. With this new recalibaration the performance of the original T c1 criterion is considerably enhanced for the prediction of crossflow transition around non-wing-like geometries. The explicit use of the local velocity vector for the integration of the sweep angle within the T c1 new reformulation makes this formulation not Galilean invariant. Nevertheless, this is a minor problem. The velocity vector can be substituted by the relative velocity. This might be not an elegant strategy, nevertheless it is very efficient, and weakly restores the Galilean invariance. This is an acceptable strategy because we are considering laminar-to-turbulence transition within the boundary layer. γ + CF , where the new recalibrated Tc1 is used for crossflow inclusion, was proven to perform very well for all the tested flows, where transition is stationary CF-dominated. Despite the calibration on the solutions of the Falkner-Skan-Cooke equations, that assume zero-spanwise gradients, we have obtained quite satisfactory results also on the sickle wing. This geometry is designed to experience strong pressure gradients in the spanwise direction. Even if the criterion has not been tested on a exhaustive set of configurations, it has been shown to be robust enough. The criterion needs to continue to be validated, possibly on test cases for which we dispose of measurements. Nevertheless, it can already be used as starting basis for future developments. Indeed, there is big room for the improvement of the local correlation-based transition models γ and γ -Re θ . We should start from the inclusion of additional transition mechanisms, that are still not accounted for, such as traveling crossflow, non-linear interactions between Tollmien-Schlichting and crossflow waves, and attachment line instabilities. To this purpose, Falkner-Skan-Cooke solutions are not sufficient to calibrate these transition criteria. New experiments are needed to characterize the transition process and to define further empirical correlations to be included within the models formulation.

Possible useful experiments for CFD developer Let us consider, for instance, traveling crossflow, for which there is no reason why the crossflow strength indicators, Ψ and He, discussed in this dissertation, should not work. Since the free-stream turbulence levels in flight are very low, transition is "believed" to be dominated by the stationary modes of crossflow instability, which further explains why the N T S -N CF procedure based on stationary crossflow predictions is so common. The role of traveling crossflow in transition, vis-à-vis stationary crossflow vortices, is determined by the receptivity mechanisms, i.e. the combination of free-stream T u levels and model roughness, as a first level approximation. The most tell-tale sign for knowing if transition is dominated by stationary crossflow (primary) instability or some traveling kind (either TS/traveling crossflow/both) is whether the transition front includes stationary wedges/sawtooth like pattern or whether it varies smoothly along the spanwise direction. Sublimating chemicals, such as naphthalene based flow visualization, can be typically employed to visualize this pattern. We do not require new experiments around new configurations, but it will be very useful to repeat transition measurements on the 6:1 prolate spheroid at different levels of free-stream turbulence intensity. This might be accomplished by selective removal/addition of turbulence screens. In addition, if a model front half with different levels of surface finish is used, surface roughness effect can be evaluated as well. Indeed, experiments under systematically varied environments showed that surface roughness represents the key parameter responsible for the initiation of stationary crossflow vortices. Repeated measurements around the prolate spheroid will help characterizing the sensitivity of the transition behavior, as well as providing data that Kreplin's experiment did not. Additional surface-based instrumentation, as microphones, can be also used to provide unsteady measurements that would give additional information about the frequency content. Also knowing the orientation of the waves, via suitable clusters of sensors, would help determine whether unsteady Tollmien-Schlichting waves or traveling crossflow modes are the important ones ahead of the transition front. More detailed information on the free-stream conditions ought to be given. The turbulence intensity at a given location of the domain is not a sufficient information and the turbulence dissipation might be the additional quantity to be measured in order to characterize the turbulence environment and free-decay. Skin friction and pressure distribution are the needed data to validate transition models. Comparison between measured and predicted forces and moments would also allow a more rigorous validation, especially for those configurations where transition is accompanied by important flow separation, such as the 6:1 prolate spheroid at high incidences. That is the minimum physics that needs to be included in a useful data set to be a benchmark of quality for further numerical simulations.

Future Work

The correlation-based transition models γ and γ -Re θ are very promising tools and they can bring CFD simulations to an higher level. Kim et al.,[START_REF] Kim | Boundary layer transition models for naval applications: capabilities and limitations[END_REF], have presented the numerical results computed by γ -Re θ for the fully appended MARIN submersible geometry BB2. They showed that the transition model performs qualitatively well in the identification of laminar and turbulent flow regions with respect to the predicted vortical structures. These results pave the way for the use of these models for real complex geometries. Nevertheless, repeated or new experiments are needed to improve the models. The aim is to add more physics within the correlations used to account for the different transition mechanisms. The ONR, Office of Naval Research, has recently launched the HIPRO, High Reynolds number Prolate spheroid, program, whose objective is to measure the flow around the 6:1 prolate spheroid at different angles of incidence and Reynolds numbers. The aim is to provide the CFD community with new and high quality data set. Even if the majority of the experimental campaigns will be run for tripped models, an experimental campaign for the untripped model is expected to be run at the David Taylor Basin at Navsea Corderock by Philippe Bardet and his team in Spring 2022. Experiments will be run around a 3m specimen for a wide range of Reynolds number, i.e. the maximum possible velocity in the basin is about 20m/s. Nonintrusive optical diagnostic techniques will be used to measure instantaneous, time resolved, flow quantities in the boundary layer and the shear stress at the wall. These techniques are based on Molecular Tagging Velocimetry (MTV) schemes, that overcome the PIV limitations in reaching the smallest scale of turbulence within the boundary layer. Details on MVT can be found in Fort et al., [START_REF] Fort | Development of long distance 2D micro-molecular tagging velocimetry (µMTV) to measure wall shear stress[END_REF]. In the meantime, a joint effort between RANS and DNS specialists would be very fruitful to further improve laminar-to-turbulence transition modeling. Indeed, using DNS we can have a deep insight on the transition physics, that RANS cannot provide us with. Also data from e N computations can be exploited for a better understanding of the physical processes. Even if these results are not numerical error free, because based on RANS solutions, they are very accurate and capable of accounting for three-dimensional effects, that the correlation based transition model do not see. For instance, velocity profiles obtained by highly resolved boundary layer codes can be used to compute the R limit for attachment line instabilities. The final purpose is to incorporate within the models formulation as much physics as possible, through enhanced empirical correlations. Among the missing physical mechanisms, the most urgent to include are, undoubtedly, traveling crossflow and attachment line instabilities. Their inclusion within the models will allow to make predictions that are useful to some degree, regardless of the specific experiments that are used for the models validation. Our wish for this document is to be further used by experimentalists as a reference of the key data that CFD model developers are keen to have to improve transition modeling. Validation of the T c1 criterion need to be pursued. To this matter, we are planning to run new simulations for the sickle wing, at different angles of incidence and Reynolds numbers. These computations are planned to be conducted on meshes exchanged within the AVT-313 group. Further developments are also needed to make these model numerically robust. It is clear that the improvement of their iterative convergence is strictly related to the coarsening of the meshes used for transition simulations. We have discussed that the models work first building turbulence and then destroying it. Thus, an ideal approach would be to construct an adaptive grid refinement criterion that follows the transition front. Using this strategy, the mesh will be unrefined in region where such spatial resolution is not needed. The highest cells density will be concentrated in the transition region. We are currently working on an adaptive grid refinement criteria based on the second derivatives of the intermittency, quantity which characterizes the transition process. Preliminary simulations are currently going on. With the purpose of improving the iterative convergence, we are also currently working on an a coupled numerical approach for the resolution of the transport equations, momentum, transition and turbulence. Within a coupled formulation, we expect a noticeable decrease of the overall non-linear iterations needed up to convergence. Abstract : This research work aims at assessing transition calculations feasibility in the context of the industrial codes. We discuss the strengths and weaknesses of the Local Correlation Based Transition Models (LCTM) γ and γ-Re_θ for two-and three-dimensional configurations, focusing on both modeling and numerical aspects. The purpose is to analyze the potential of these RANS-based models as predictive tools, that can handle automatically and autonomously laminar-to-turbulence transition. We evaluate the transition models performance on an extensive number of test cases, covering a wide range of transition mechanisms.

In the case of 3D configurations, a crucial point is the modeling of transition due to stationary crossflow modes, which are the main threedimensional transition mechanism in a low free-stream turbulence environment. For this purpose, we present in this dissertation an original re-calibration of the Tc1 crossflow criterion, proposed by Menter and Smirnov in 2014 and based on a local formulation of the renowned empirical C1 criterion of Daniel Arnal. This re-calibrated criterion notably extends the existing γ model for stationary crossflow transition prediction around threedimensional non-wing-like geometries.

39 4. 1 45 4. 2 49 5. 1 63 5. 2 64 5. 4 65 5. 5 89 6. 2 93 6. 4 96 6. 6 96 6. 7 97 6. 8 R 110 7.5 6 : 1

 39145249163264465589293496696797811061 Example of the spatial discretization. (Left) Typical unstructured control volume. (Center) Cell C 0 and its neighborhood. (Right) Centered face reconstruction notations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Evolution of Q in the inflow direction. . . . . . . . . . . . . . . . . . . . . . . . . 47 4.3 NVD diagram: UDS, CDS and GDS schemes. . . . . . . . . . . . . . . . . . . . . 48 4.4 NVD diagram: AVLSMART scheme for non-uniform grids. . . . . . . . . . . . . Flat Plate. O-grid topology mesh. Entire computational domain and close-up of the leading edge of the coarsest mesh. . . . . . . . . . . . . . . . . . . . . . . . . Flat Plate: T3A. Skin friction profiles predicted by γ and γ -Re θ , and k -ω SST (2003) turbulence model along the flat plate. Numerical results are compared to experimental measurements (black filled squares) and the Blasius solution (dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.3 Flat Plate: T3A. Skin friction profiles predicted by γ and γ -Re θ along the flat plate with numerical and input uncertainties. Numerical results are compared to experimental measurements (black filled squares). . . . . . . . . . . . . . . . . . . Flat Plate: T3A. Validation comparison error E(C f ) in the transition region with numerical and input uncertainties. Numerical results are predicted by γ and γ -Re θ models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Flat Plate: T3A. Streamwise velocity profiles at three different sections x/L plotted on the experimental results. Experimental, numerical, and input uncertainties are reported at each point of the validation set. Predictions are computed using γ and γ -Re θ models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.6 Flat Plate: T3A-. Skin friction profiles predicted by γ and γ -Re θ , and k -ω SST (2003) turbulence model along the flat plate. Numerical results are compared to experimental measurements (black filled squares) and the Blasius solution (dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.7 Flat Plate: T3A-. Turbulence intensity contours in the proximity of the transition location. The dotted line represents the boundary layer edge as computed using the formula from Griffin et al., [112]. The boundary layer edge definition is given in Appendix B. Predictions by both γ and γ -Re θ models. . . . . . . . . . . . . 5.8 Eppler387: α = 1 • . Illustration of the entire computational domain, close-up of the mesh topology around the foil and of the wake of the coarsest mesh. . . . . . 5.9 Eppler387: α = 1 • . C p and C f distributions on the airfoil as computed using γ and γ -Re θ transition models. Numerical prediction by k -ω SST (2003) turbulence model are given as a reference. Predictions are compared to measurements from Mcghee et al.(black squares) and Cole et al.(white squares). . . . . . . . . . . . . 5.10 Eppler387: α = 1 • . Pressure coefficient C p distributions with numerical uncertainties around the foil. Numerical results are computed by different flow solvers, indicated by different symbols, using γ and γ -Re θ transition models. . . . . . . 5.11 Eppler387: α = 1 • . Numerical lift coefficient plotted against the grid refinement h i /h 1 . Results are computed using γ and γ -Re θ and different flow solvers. . . 5.12 Eppler387: α = 1 • . Numerical pressure C d,p and friction C d,f drag coefficient plotted against the grid refinement h i /h 1 . Results are computed using γ and γ -Re θ and different flow solvers. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Eppler387: α = 7 • . C p distribution on the airfoil as computed using γ transition model. Predictions are compared to measurements from Mcghee et al. (black squares) and Cole et al. (white squares). . . . . . . . . . . . . . . . . . . . . . . 5.14 NACA0015: α = 5 • . Illustration of the entire computational domain, close-up of the mesh topology around the foil and of the wake of the coarsest mesh. . . . . . 5.15 NACA 0015: α = 5 • . Measured time-averaged C f distribution on the surface. . . 5.16 NACA 0015: α = 5 • . Measured ans predicted C f distributions on the upper surface of the foil computed by γ -Re θ and γ along with the experimental measurements. Results are presented with and without uncertainties. Predictions by k -ω SST (2003) are given as a reference. . . . . . . . . . . . . . . . . . . . . . . 5.17 NACA 0015: α = 5 • . Predicted C p profiles on the upper surface of the foil computed by γ -Re θ and γ models. . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 NACA0015: α = 5 • . Contours of turbulence intensity T u within the laminar separation bubble. The velocity streamlines indicate the re-circulation zone, while the dashed line on the top is the boundary layer edge as approximated using the formula in Appendix B. Prediction by γ and γ -Re θ transition models. . . . . . 5.19 Number of non-linear iterations on the coarsest meshes up to convergence for the test cases discussed in the present chapter: Flat Plate T3A and T3A-(Grid5), Eppler387 (E387) at α = 1 • (Grid5) and NACA0015 (N-0015) at α = 5 • (Grid9). 5.20 NACA0015: α = 5 • . Evolution of the normalized residuals in L 2 norm of all the transported quantities. Computation run on the coarsest grid of the set (Grid9) using γ model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx LIST OF FIGURES 5.21 F onset,3 , as in the γ -R θ transition model, vs the eddy viscosity ratio R T . The position of the kink is pushed backward or forward with respect to the constant by which R T is divided. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 6.1 C1 crossflow transition criterion proposed by Arnal, Re δ 2 t = f (H 12 ). Experimental results are denoted by the symbols: Poll (white circles), Schmitt (white squares), ONERA/CERT/DERAT (black circles) and Boltz (black squares). Figure reproduced from Arnal, [16]. . . . . . . . . . . . . . . . . . . . . . . . . . . . Crossflow indicators ratio R(β, φ) = 0.4Ψ/X, as a function of the Hartree parameter, β, and the sweep angle, φ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.3 Ratio R(β, φ) = 0.684Ψ/(XF (H 12 )), as a function of the Hartree parameter, β, and the sweep angle, φ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Crossflow indicators ratio R distribution, projected on the β-plane and plotted as a function of the new λ θ,CF . Planes are colored by the contours of the sweep angle φ. Left: uncorrected ratio R = (0.684Ψ)/(XF (H 12 )). Right: corrected ratio R = (G(λ θ,CF )Ψ)/(XF (H 12 )) . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.5 Transition location on the infinite swept wing NLF(2)-0415 as a function of Re. Figure taken from the paper of Menter & Smirnov, [15]. . . . . . . . . . . . . . . Skin friction distribution C f around the 6:1 prolate spheroid. Figure taken from the paper of Menter & Smirnov, [15]. . . . . . . . . . . . . . . . . . . . . . . . . . Pressure gradient parameter λ θ plotted as a function of the Hartree parameter β and the sweep angle φ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . = (G(λ θ,CF )Ψ)/(XF (H 12 )) vs the fitted pressure gradient parameter λ CF dv dy , φ and the sweep angle φ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 7.1 6:1 prolate spheroid. Sketch of the specimen with probes position. Figure taken from Kreplin, [129]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.2 6:1 prolate spheroid: α = 15 • . Full domain and close up on the O-block surrounding the body for the coarsest mesh. . . . . . . . . . . . . . . . . . . . . . . 109 7.3 6:1 prolate spheroid: α = 15 • . Closed up on the first half of the spheroid: grid at the surface and at the symmetry plane for the coarsest mesh. . . . . . . . . . . . 109 7.4 6:1 prolate spheroid: α = 15 • . Coarsest mesh in the wall normal direction at the section X/L = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Prolate Spheroid: α = 15

  Figure taken from Kruse et al., [153]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 7.27 Sickle Wing: α = -2.6 • . Grid used for the flow simulation around the sickle wing. Total computational domain and mesh at the surface. In Fig.(7.27b), the solid red lines represent the kinks, where the geometrical sweep angle changes. . 130 7.28 Sickle wing: α = -2.6• . Predicted skin friction contours on the upper and lower surface of the sickle wing by γ and γ + CF models. The black dotted line is the extracted experimental transition location from Petzold et al., the black solid line is from Kruse et al.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 7.29 Sickle wing. Computational mesh used from Kim et al., for their transitional simulation around the sickle wing. Figure reproduced by Kim et al., [13]. . . . . 132 7.30 Sickle wing: α = -2.6 • . Comparison of measured and calculated pressure coefficient distributions at each sweep section of the model, from Section A to Section C. Pressure is extracted along the dotted lines (pressure taps) represented in Fig.(7.26). Numerical results are computed by γ and γ + CF transition variants. 134 List of Tables 4.1 AVLSMART scheme implementation for non-uniform grids. . . . . . . . . . . . 5.1 2D Test cases details: geometry, angle of attack α and Reynolds number. The turbulence intensity at the inlet, T u in , and the transition mechanisms are reported as well. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Flat Plate. Details of the five structured O-grids: total number of cells, N cells , number of cells on the flat plate surface, N cells-plate , grid refinement ratio, r i = h i /h, and y + max values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Flat Plate: T3A. Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5.4 Flat Plate: T3A. Multivariate metric for the skin friction in different regions and all along the flat plate for γ model results. . . . . . . . . . . . . . . . . . . . . . . 5.5 Flat Plate: T3A. Multivariate metric for the skin friction in different regions and all along the flat plate for γ -Re θ model results. . . . . . . . . . . . . . . . . . . 5.6 Flat Plate: T3A. Multivariate metric for the streamwise velocity at different sections of the flat plate for γ model results. . . . . . . . . . . . . . . . . . . . . . 5.7 Flat Plate: T3A. Multivariate metric for the streamwise velocity at different sections of the flat plate for γ -Re θ model results. . . . . . . . . . . . . . . . . . 5.8 Flat Plate: T3A-. Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Eppler387: α = 1 • , 7 • . Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Eppler387. Details of the five structured C-O-grids: total number of cells, N cells , number of cells on the flat plate surface, N cells-plate , grid refinement ratio, r i = h i /h, and y + max values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Eppler387: α = 1 • . Experimental and numerical efforts predicted using γ and γ -Re θ and k -ω SST (2003). . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii xxiv LIST OF TABLES 5.12 NACA0015: α = 5

-

  Sommerfeld and Squire Equations. Let us consider the case of a parallel base flow -→ U = U (y) -→ e x . The parallel flow assumption implies that the stability of the flow at a certain position (x, z) is independent of the flow history. It depends on the local condition at the specific location. If this mean flow is introduced into Eq.(2.5) and non-linear terms are neglected, the resulting equations written by components are:

Figure 2 . 1 :

 21 Figure 2.1: (Re, α), Reynolds number and streamwise wavenumber, diagram for Blasius boundary layer flow: contours of constant growth rate c i . The red dot corresponds to the critical Reynolds number, Re c = 519.4, and streamwise wavenumber α c = 0.303. They lie on the neutral curve c i = 0 that divides the unstable (shaded) and stable region.Figure reproduced from Schmid, [23].

  Figure 2.1: (Re, α), Reynolds number and streamwise wavenumber, diagram for Blasius boundary layer flow: contours of constant growth rate c i . The red dot corresponds to the critical Reynolds number, Re c = 519.4, and streamwise wavenumber α c = 0.303. They lie on the neutral curve c i = 0 that divides the unstable (shaded) and stable region.Figure reproduced from Schmid, [23].

Figure 2 . 2 :

 22 Figure 2.2: Smoke visualization of streamwise streaks appearing in a boundary layer subjected to moderate level of freestream turbulence intensity T u = 2.2%. The flow, coming from the left, is captured at a certain distance from the leading edge and it transitions (breakdown and turbulence spot formations) in the right edge of the image. Figure reproduced from Matsubara et al., [26].

Fig.( 2

 2 .3) as a reference.

Figure 2 . 3 :

 23 Figure 2.3: Different instability mechanisms on the swept wing. Figure reproduced from Shahriari, [27].

Figure 2 . 4 :

 24 Figure 2.4: Schematic representation of the transition process due to T-S waves. Figure is reproduced from White & Corfield, [28].

Figure 2 . 5 :

 25 Figure 2.5: Görtler vortices at near a stagnation point (S) lying on a convex surface. R is the curvature radius of the external streamline. Figure reproduced from Hirschel et al., [31].

Figure 2 . 6 :

 26 Figure 2.6: Schematic representation of the flow near the leading edge of a swept wing, reproduced from Poll, [36].

Figure 2 . 7 :

 27 Figure 2.7: Infinite swept wing notations. Figure re-adapted from Arnal et al., [16].

Figure 2 . 8 :

 28 Figure 2.8: Development of velocity profiles on a swept wing when suffering from an inflection of the pressure gradient. W 1 (y) and U 1 (y) are the crossflow and main-flow profile projected onto a coordinate system relative to the external inviscid streamline. β is the angle between the wall shear stress and the external streamlines. Figure re-adapted from Yiming et al., [38].

Figure 2 . 10 :

 210 Figure 2.10: Maximum growth of laminar fluctuations in a laminar boundary before transition for a flow on a flat plate. Figure reproduced from Mayle & Schulz , [47].

  (a) Small length scale vortices generation at different time steps visualized as λ2 isosurfaces. (b) λ2 isosurfaces and velocity perturbation contours.

Figure 2 . 11 :

 211 Figure 2.11: Visualization of small length vortices at different time steps through λ 2 isosurfaces: view from the top in Fig.(2.11a), view in the (y, z) plane for fixed x in Fig.(2.11b) with the correspondent velocity perturbations contours. Simulations are performed at Re = 1000 on a grid of 1920 × 241 × 128 points. Figures are reproduced from Liu et al., [51].

For

  3D configurations the strategy often used is the N T S -N CF stability boundary, where N T S is the critical N factor for T-S waves, N CF for crossflow instabilities. An example of stability boundary is shown in Fig.(2.12), reproduced from Stock,[START_REF] Stock | e N transition prediction in three-dimensional boundary layers on inclined prolate spheroids[END_REF]. Through this diagram, we are able to account for the interaction between T-S and CF. In this specific case, it can be observed that N CF factor decreases linearly as N T S increases. In terms of transition process, this translates into a reduction of the stable region because of the simultaneous excitation of the two different disturbances. The intersection of the boundary describes which kind of instability is triggering transition. For a simplified two-N strategy, N CF is computed taking into account CF modes with frequency f = 0 (stationary crossflow) and different wave length λ, and it depends mainly on the surface roughness, while the N T S is computed considering free-stream T-S modes with wave number β = 0 for different frequencies f. No curvature effects are taken into account, nor traveling crossflow modes at higher turbulence level.

Figure 2 . 12 :

 212 Figure 2.12: N T S and N CF stability diagram for the 6:1 prolate spheroid. Figure is reproduced from Stock, [58].

Figure 2 . 13 :

 213 Figure 2.13: Scaled vorticity Reynolds number Re V for a Blasius boundary layer, as computed using Falkner-Skan equations. Re θ is the corresponding momentum thickness Reynolds number and 2.2 is the calibration constant set to have a maximum of one in the boundary layer.

Chapter 3 Local

 3 Correlation Transition Models: γ -Re θ & γ This chapter presents the mathematical formulation of the RANS Local Correlation Transition Models γ -Re θ and γ. The equations of the turbulence model k -ω SST 2003, which the transition model are coupled to, are given as well.

Figure 4 . 1 :

 41 Figure 4.1: Example of the spatial discretization. (Left) Typical unstructured control volume. (Center) Cell C 0 and its neighborhood. (Right) Centered face reconstruction notations.

  rests on Eq.(4.14) and Eq.(4.13) plus the discontinuity constraints of Eq.(4.10) and Eq.(4.11).

Figure 4 . 2 :

 42 Figure 4.2: Evolution of Q in the inflow direction..

  Fig.(4.3). In Fig.(4.3a), it can be observed that first wind order scheme (UDS) fulfills all these conditions, whereas the Centered Differencing Scheme (CDS) does not respect these criteria in the interval QC ∈]∞, 0[ ∪ ]1, +∞[. This implies that CDS is only useful in the range 0 ≤ QC ≤ 1.

  UDS and CDS schemes in NVD diagram. GDS scheme in NVD diagram.

Figure 4 . 3 :

 43 Figure 4.3: NVD diagram: UDS, CDS and GDS schemes.

Figure 4 . 4 :

 44 Figure 4.4: NVD diagram: AVLSMART scheme for non-uniform grids.

  (a) Entire computational domain. (b) Close-up of the leading edge.

Figure 5 . 1 :

 51 Figure 5.1: Flat Plate. O-grid topology mesh. Entire computational domain and close-up of the leading edge of the coarsest mesh.

Re 3 :

 3 T u in (%) ν t in /ν T u(%) Flat Plate: T3A. Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%).

Figure 5 . 2 :

 52 Figure 5.2: Flat Plate: T3A. Skin friction profiles predicted by γ and γ -Re θ , and k -ω SST (2003) turbulence model along the flat plate. Numerical results are compared to experimental measurements (black filled squares) and the Blasius solution (dashed line).

  Fig.(5.3b), the input uncertainties are included.

  (a) C f profile with numerical uncertainties. (b) C f profile with numerical and input uncertainties.

Figure 5 . 3 :

 53 Figure 5.3: Flat Plate: T3A. Skin friction profiles predicted by γ and γ -Re θ along the flat plate with numerical and input uncertainties. Numerical results are compared to experimental measurements (black filled squares).

Figure 5 . 4 :

 54 Figure 5.4: Flat Plate: T3A. Validation comparison error E(C f ) in the transition region with numerical and input uncertainties. Numerical results are predicted by γ and γ -Re θ models.

  (a) x/L = 0.01006 position along the skin friction profile. (b) Streamwise velocity profile. (c) x/L = 0.02035 position along the skin friction profile. (d) Streamwise velocity profile. (e) x/L = 0.05273 position along the skin friction profile. (f) Streamwise velocity profile.

Figure 5 . 5 :

 55 Figure 5.5: Flat Plate: T3A. Streamwise velocity profiles at three different sections x/L plotted on the experimental results. Experimental, numerical, and input uncertainties are reported at each point of the validation set. Predictions are computed using γ and γ -Re θ models.

  The skin friction distributions in the first half of the flat plate are shown in Fig.(5.6). Numerical results computed by γ and γ -Re θ are compared to measurements. Turbulent and Blasius solutions are given as a reference, as well.

Figure 5 . 6 :

 56 Figure 5.6: Flat Plate: T3A-. Skin friction profiles predicted by γ and γ -Re θ , and k -ω SST (2003) turbulence model along the flat plate. Numerical results are compared to experimental measurements (black filled squares) and the Blasius solution (dashed line).

  Fig.(5.7). Compared to γ -Re θ , the γ derivative of turbulence kinetic energy in the boundary layer are much stronger. Despite the later activation of the F onset function for γ results, as the intermittency starts to be produced within the boundary layer, the flow becomes rapidly turbulent. The differences between γ and γ -Re θ predictions stand in the F length definition. F length , as a function of Re θt , is very small at high Re θt , due to the low free-stream turbulence intensity. On the other hand, in γ formulation, F length is set to a constant and it is equal to 100. This behavior was already noticed in the original publication of γ model formulation by Menter et al.. A comparison of the velocity profiles at different sections is worthless, because the two models predictions significantly differ in the transition and turbulent region, as a result of the different production of turbulence kinetic energy within the boundary layer.

γ

  model calibration stage. Nevertheless, in spite of the impact on the local numerical results, qualitatively the two models performs very similarly and the predicted transition positions by γ and γ -Re θ coincide. (a) γ. (b) γ -Re θ .

Figure 5 . 7 :

 57 Figure 5.7: Flat Plate: T3A-. Turbulence intensity contours in the proximity of the transition location. The dotted line represents the boundary layer edge as computed using the formula from Griffin et al., [112]. The boundary layer edge definition is given in Appendix B. Predictions by both γ and γ -Re θ models.

  (a) Entire computational domain. (b) Close-up of the grid around the airfoil. (c) Close-up of wake. (d) Entire computational domain.

Figure 5 . 8 :

 58 Figure 5.8: Eppler387: α = 1 • . Illustration of the entire computational domain, close-up of the mesh topology around the foil and of the wake of the coarsest mesh.

5. 5 . 1 α = 1 •

 511 Let us discuss the results for α = 1 • . γ -Re θ and γ predictions are shown in Fig.(5.9). The pressure coefficient distribution C p is shown in Fig.(5.9a) along with the numerical uncertainties.

  (a) Cp distribution on the airfoil with numerical uncertainties. (b) C f distribution on the airfoil with numerical uncertainties.

Figure 5 . 9 :

 59 Figure 5.9: Eppler387: α = 1 • . C p and C f distributions on the airfoil as computed using γ and γ -Re θ transition models. Numerical prediction by k -ω SST (2003) turbulence model are given as a reference. Predictions are compared to measurements from Mcghee et al.(black squares) and Cole et al.(white squares).

C

. 37 )

 37 a is a constant and it varies with respect to k -ω SST version: a = 10 for k -ω SST version 2003 and a = 20 for k -ω SST version 1994.

( a )

 a Pressure coefficient -Cp profiles using γ model. (b) Pressure coefficient -Cp profiles using γ -Re θ model.

Figure 5 .

 5 Figure 5.10: Eppler387: α = 1 • . Pressure coefficient C p distributions with numerical uncertainties around the foil. Numerical results are computed by different flow solvers, indicated by different symbols, using γ and γ -Re θ transition models.

  C l computed using γ. (b) C l computed using γ -Re θ .

Figure 5 . 11 :

 511 Figure 5.11: Eppler387: α = 1 • . Numerical lift coefficient plotted against the grid refinement h i /h 1 . Results are computed using γ and γ -Re θ and different flow solvers.

  (a) C d,p computed using γ. (b) C d,p computed using γ -Re θ . (c) C d,f computed using γ. (d) C d,f computed using γ -Re θ .

Figure 5 . 12 :

 512 Figure 5.12: Eppler387: α = 1 • . Numerical pressure C d,p and friction C d,f drag coefficient plotted against the grid refinement h i /h 1 . Results are computed using γ and γ -Re θ and different flow solvers.

( a )

 a Cp computed using γ.

Figure 5 .

 5 Figure 5.13: Eppler387: α = 7 • . C p distribution on the airfoil as computed using γ transition model. Predictions are compared to measurements from Mcghee et al. (black squares) and Cole et al. (white squares).

( a )

 a Entire computational domain. (b) Close-up of the grid around the airfoil.(c) Close-up of wake.

Figure 5 .

 5 Figure 5.14: NACA0015: α = 5 • . Illustration of the entire computational domain, close-up of the mesh topology around the foil and of the wake of the coarsest mesh.

5. 6 . 1 α

 61 = 5 •

Figure 5 .

 5 Figure 5.15: NACA 0015: α = 5 • . Measured time-averaged C f distribution on the surface.

  (a) C f profiles without uncertainties. (b) C f profiles with respective uncertainties.

Figure 5 .

 5 Figure 5.16: NACA 0015: α = 5 • . Measured ans predicted C f distributions on the upper surface of the foil computed by γ -Re θ and γ along with the experimental measurements. Results are presented with and without uncertainties. Predictions by k -ω SST (2003) are given as a reference.

Figure 5 .

 5 Figure 5.17: NACA 0015: α = 5 • . Predicted C p profiles on the upper surface of the foil computed by γ -Re θ and γ models.

  (a) γ -Re θ . (b) γ.

Figure 5 .

 5 Figure 5.18: NACA0015: α = 5 • . Contours of turbulence intensity T u within the laminar separation bubble. The velocity streamlines indicate the re-circulation zone, while the dashed line on the top is the boundary layer edge as approximated using the formula in Appendix B. Prediction by γ and γ -Re θ transition models.

  2 of each transport equation for NACA0015 at α = 5 • using γ model are shown in Fig.(5.20) as example.

Figure 5 . 19 :

 519 Figure 5.19: Number of non-linear iterations on the coarsest meshes up to convergence for the test cases discussed in the present chapter: Flat Plate T3A and T3A-(Grid5), Eppler387 (E387) at α = 1 • (Grid5) and NACA0015 (N-0015) at α = 5 • (Grid9).

Figure 5 .

 5 Figure 5.20: NACA0015: α = 5 • . Evolution of the normalized residuals in L 2 norm of all the transported quantities. Computation run on the coarsest grid of the set (Grid9) using γ model.

Figure 5 . 21 :

 521 Figure 5.21: F onset,3 , as in the γ -R θ transition model, vs the eddy viscosity ratio R T . The position of the kink is pushed backward or forward with respect to the constant by which R T is divided.

(6. 5 )

 5 Re δ 2 t expression is shown in Fig.(6.1), where it can be observed how for low values of H 12 , the curve approaches the value of 150 proposed byBeasley.

Figure 6 . 1 :

 61 Figure 6.1: C1 crossflow transition criterion proposed by Arnal, Re δ 2 t = f (H 12 ). Experimental results are denoted by the symbols: Poll (white circles), Schmitt (white squares), ON-ERA/CERT/DERAT (black circles) and Boltz (black squares). Figure reproduced from Arnal, [16].

Figure 6 . 2 :

 62 Figure 6.2: Crossflow indicators ratio R(β, φ) = 0.4Ψ/X, as a function of the Hartree parameter, β, and the sweep angle, φ.

Figure 6 . 3 :

 63 Figure 6.3: Ratio R(β, φ) = 0.684Ψ/(XF (H 12 )), as a function of the Hartree parameter, β, and the sweep angle, φ.

  Fig.(6.4) shows the uncorrected ratio R = (0.684Ψ)/(XF (H 12 )) with the cubic polynomial g(λ θ,CF ) on the right, Fig.(6.4a, and the corrected R = (G(λ θ,CF )Ψ)/(XF (H 12 )) on the left, Fig.(6.4b). The spread around the target value of 1 is consistently reduced to around 25% for sweep angles φ lower than 60• , which is the upper limit chosen by Menter & Smirnov for the calibration of their criterion. Nevertheless, the error committed on higher angles of attack is still very high.

  (a) Distribution of the uncorrected indicator ratio R. Black line is the correction polynomial g(λ θ,CF ). (b) Distribution of the corrected indicator ratio R.

Figure 6 . 4 :

 64 Figure 6.4: Crossflow indicators ratio R distribution, projected on the β-plane and plotted as a function of the new λ θ,CF . Planes are colored by the contours of the sweep angle φ. Left: uncorrected ratio R = (0.684Ψ)/(XF (H 12 )). Right: corrected ratio R = (G(λ θ,CF )Ψ)/(XF (H 12 ))

Figure 6 . 5 :

 65 Figure 6.5: Transition location on the infinite swept wing NLF(2)-0415 as a function of Re.Figure taken from the paper of Menter & Smirnov, [15].

  Figure 6.5: Transition location on the infinite swept wing NLF(2)-0415 as a function of Re.Figure taken from the paper of Menter & Smirnov, [15].

Figure 6 . 6 :

 66 Figure 6.6: Skin friction distribution C f around the 6:1 prolate spheroid. Figure taken from the paper of Menter & Smirnov, [15].

Figure 6 . 7 :

 67 Figure 6.7: Pressure gradient parameter λ θ plotted as a function of the Hartree parameter β and the sweep angle φ.

  λ CF used in the polynomial fitting is evaluated at η = η max . The benefit of introducing the sweep angle φ is a considerable reduction of the spread of the ratio R = (G(λ θ,CF )Ψ)/(XF (H 12 )) around the value of one, as it can be observed in Fig(6.8), where R is plotted against the two dependent variables λ CF dv dy , φ an φ. The maximum deviation is less than 10%.

Figure 6 . 8 :

 68 Figure 6.8: R = (G(λ θ,CF )Ψ)/(XF (H 12 )) vs the fitted pressure gradient parameter λ CF dv dy , φ and the sweep angle φ.

. 1 )

 1 Measurements at the surface were performed using surface hot film probes: wall shear stress magnitude is derived from the heat transfer rates of the films of each probe. The probes are positioned at 12 different stations along the major axes, as shown in Fig.(7.1).

Figure 7 .

 7 Figure 7.1: 6:1 prolate spheroid. Sketch of the specimen with probes position. Figure taken from Kreplin, [129].

  .2). The coarsest grid from different points of view at α = 15 • is shown in Fig.(7.2) and (7.3). The mesh in the wall normal direction at the plane X/L = 0.5 is shown in Fig.(7.4).

  (a) Full Computational Domain. (b) Zoom on the O-block surrounding the spheroid.

Figure 7 .

 7 Figure 7.2: 6:1 prolate spheroid: α = 15 • . Full domain and close up on the O-block surrounding the body for the coarsest mesh.

Figure 7 .

 7 Figure 7.3: 6:1 prolate spheroid: α = 15 • . Closed up on the first half of the spheroid: grid at the surface and at the symmetry plane for the coarsest mesh.

Figure 7 . 4 : 6 :

 746 Figure 7.4: 6:1 prolate spheroid: α = 15 • . Coarsest mesh in the wall normal direction at the section X/L = 0.

N 1

 1 cells N surface y + max Prolate Spheroid: α = 15 • . Mesh details of Grid5bis and Grid5. N x is measured along the upper side of the surface and N φ is measured along the plane located at half of the longitudinal length of the surface.

  (a) X/L -φ plane, Grid5bis: C f distribution. (b) X/L -φ plane, Grid5: C f distribution.

Figure 7 . 5 : 6 : 1

 7561 Figure 7.5: 6:1 Prolate Spheroid: α = 15 • . Skin friction C f distribution on the 6:1 prolate spheroid computed using γ model without crossflow inclusion on Grid5 and Grid5bis of Table(7.2).

  C d,p .

Figure 7 . 6 : 6 : 1

 7661 Figure 7.6: 6:1 Prolate Spheroid: α = 15 • . Convergence of pressure and viscous drag for γ + CF formulation with grid refinement. p is the estimated convergence order of the discretization error, as described in Sec.(5.2.1).

  C d,f .

78 (

 78 b) C d,p .

Figure 7 . 7 : 6 : 1

 7761 Figure 7.7: 6:1 Prolate Spheroid: α = 15 • . Convergence of pressure and viscous drag for γ -Re θ + CF formulation with grid refinement. p is the estimated convergence order of the discretization error, as described in Sec.(5.2.1).

  .2). The results for α = 15• shown in the following are computed on Grid1.The experimental C f contours are shown in Fig.(7.8) in the X/L -φ plane. For this angle of attack and Reynolds number, a zone of CF-dominated transition is observed in correspondence of the middle of prolate spheroid. The kink at about 20% of the length of the spheroid, at φ = 130 • , marks the change of the transition process. Indeed, transition on the leeward side is driven by T-S waves. The numerical skin friction contours computed by γ +CF and γ -Re θ +CF are presented in Fig.(7.9) and(7.11), along with the pressure coefficient C p distribution. We also report, in Fig.(7.10), the numerical predictions obtained with the original crossflow T c1 criterion proposed by Menter & Smirnov,[15] and described in Sec.(6.3). This γ model variant is referred to as γ + T c1-MS. In the original publication from Menter & Smirnov, neither the definition of the F onset,CF , based on their version of the T c1 criterion, nor the modification of the F onset function in the transport equation for the intermittency γ, Eq.(3.53), are mentioned. For this reason, we have included the T c1-MS criterion within γ formulation as described in Eq.(6.46)-(6.50). Finally, the differences between the two γ versions, γ + CF and γ + T c1-MS, as implemented in ISIS-CFD, stand in the construction of the function G, Eq.(6.46), which accounts for the pressure gradient in the streamwise direction and the definition of the pressure gradient parameter. γ + T c1-MS results are only included in the present discussion for α = 15 • . The purpose is to show the different performance between the two γ versions and understand the benefits of the inclusion of the local sweep angle within the new T c1 formulation. If we refer to the skin friction contours shown in Fig.(7.5b) 1 , as flow solution without any CF inclusion, we understand that a correlation that accounts for crossflow instabilities is fundamental. A zone of pure crossflow transition is observed at the middle of the spheroid down to the windward side. CF waves contribute to transition also close to the leeward symmetry plane. The new T c1 criterion in γ + CF , Fig.(7.9), performs significantly better than the original T c1 in γ + T c1-MS, Fig.(7.10), in terms of skin friction predictions. Including the local sweep angle, γ + CF predicts transition in a bigger region on the middle and windward side of the spheroid.

Figure 7 . 8 : 6 : 1

 7861 Figure 7.8: 6:1 Prolate Spheroid: α = 15 • . Measured skin friction C f distribution in the X/L-φ plane.

Figure 7 . 9 : 6 : 1

 7961 Figure 7.9: 6:1 Prolate Spheroid: α = 15 • . Predicted skin friction C f and pressure C p coefficients distributions computed using γ + CF model.

Figure 7 .

 7 Figure 7.10: 6:1 Prolate Spheroid: α = 15 • . Predicted skin friction C f and pressure C p coefficients distributions computed using γ + T c1-MS model.

  (a) X/L -φ plane: C f distribution.. (b) X/L -φ plane: Cp distribution.

Figure 7 .

 7 Figure 7.11: 6:1 Prolate Spheroid: α = 15 • . Predicted skin friction C f and pressure C p coefficients distribution computed using γ -Re θ + CF model.

  (a) X/L = 0.139. (b) X/L = 0.223.

Figure 7 .

 7 Figure 7.12: 6:1 Prolate Spheroid: α = 15 • . Experimental and numerical girthwise distribution of the C f coefficient at the sections X/L = 0.139 and X/L = 0.223. Numerical results are computed with γ +CF and γ -Re θ +CF . At each validation point, the experimental uncertainty is reported.

  (a) X/L = 0.309. (b) X/L = 0.480.

Figure 7 .

 7 Figure 7.13: 6:1 Prolate Spheroid: α = 15 • . Experimental and numerical girthwise distribution of the C f coefficient at the sections X/L = 0.309 and X/L = 0.480. Numerical results are computed with γ +CF and γ -Re θ +CF . At each validation point, the experimental uncertainty is reported.

Figure 7 .

 7 Figure 7.15: 6:1 Prolate Spheroid: α = 15 • . Skin friction lines as computed by γ + CF seen from different points of view, leeward (LW), top and windward (WW) sides.

7. 2 . 2 α = 5 •

 225 The 6:1 prolate spheroid at α = 5 • undergoes a different transition process than α = 15• . Nevertheless, the conclusions on the transition models γ +CF and γ -Re θ +CF performance are very similar. For the lowest angle of attack, the zone of pure crossflow transition is considerably reduced, and concentrated in the middle of the spheroid. Transition on the windward side is due to the non-linear interaction between TS and CF instabilities, and, as mentioned at the end of Chapter 6, this is a physical feature that the new γ + CF variant does not account for, nor γ -Re θ + CF . The experimental C f contours in the x/L -φ plane are shown in Fig.(7.16). Predictions by γ + CF are shown in Fig.(7.17) and γ -Re θ + CF in Fig.(7.18). The pressure distribution computed by the two transition models is very similar, Fig.(7.17b) and Fig.(7.18b). The region of the favorable pressure gradient is considerably reduced with respect to α = 15 • , as expected because of the lower angle of attack. The inclusion of a crossflow criterion contributes to move upward the transition front in the region included in the range 70 • < φ < 160 • , nevertheless transition on the windward side is still not predicted. The flow remains laminar until the section X/L = 0.395. At this position transition is predicted by both transition models, promoted by the crossflow transition criteria. At this position the flow should be laminar, according to the experimental results. The same behavior was observed at α = 15 • . The transition front is pushed outward, towards the leeward symmetry plane, by γ + CF when compared γ -Re θ + CF , Fig.(7.19a). The latter is in better agreement with the experimental transition front as predicted at the next section X/L = 0.480, as plotted in Fig.(7.19b). Transition close to the symmetry plane at the leeward side is triggered by the T-S onset criteria of the two models, respectively. As similarly observed for α = 15 • , from the section X/L = 0.652 on, Fig.(7.20), the experimental measurements deviate considerably from the numerical predictions. Towards the trailing edge, Fig.(7.20b), the transition front predicted by γ + CF is slightly fuller than γ -Re θ + CF , inverting the trend observed at the previous sections, Fig.(7.20a). Both transition models, predict higher values of skin friction in the turbulent region. In Fig.(7.21), the skin friction lines as computed by the γ + CF are plotted. Results for γ -Re θ + CF are very similar and therefore omitted. As supposed by the pressure distribution contours, the wall streamlines curvature in the first-half of the spheroid confirm that the crossflow is considerably weaker, compared to α = 5 • . The two predictions are very similar, and no characteristic topological sign of separation is observed. Indeed, the flow detaches at the end of the body, as indicated by the minimum C f .

Figure 7 .

 7 Figure 7.16: 6:1 Prolate Spheroid: α = 5 • . Measured skin friction coefficient C f distribution in the X/L -φ plane.

  (a) X/L -φ plane: C f distribution. (b) X/L -φ plane: Cp distribution.

Figure 7 . 17

 717 Figure 7.17: 6:1 Prolate Spheroid: α = 5 • . Predicted skin friction C f and pressure C p coefficients distributions computed using γ + CF model.

: 6 : 1

 61 Figure 7.17: 6:1 Prolate Spheroid: α = 5 • . Predicted skin friction C f and pressure C p coefficients distributions computed using γ + CF model.

  (a) X/L -φ plane: C f distribution. (b) X/L -φ plane: Cp distribution.

Figure 7 .

 7 Figure 7.18: 6:1 Prolate Spheroid: α = 5 • . Predicted skin friction C f and pressure C p coefficients distributions computed using γ -Re θ + CF model.

  (a) X/L = 0.395. (b) X/L = 0.480.

Figure 7 .

 7 Figure 7.19: 6:1 Prolate Spheroid: α = 5 • . Experimental and numerical girthwise distribution of the C f coefficient at the sections X/L = 0.395 and X/L = 0.480. Numerical results are computed with γ + CF and γ -Re θ + CF models. At each validation point, the experimental uncertainty is reported.

  (a) X/L = 0.652. (b) X/L = 0.936.

Figure 7 .

 7 Figure 7.20: 6:1 Prolate Spheroid: α = 5 • . Experimental and numerical girthwise distribution of the C f coefficient at the sections X/L = 0.652 and X/L = 0.936. Numerical results are computed with γ + CF and γ -Re θ + CF models. At each validation point, the experimental uncertainty is reported.

Figure 7 .

 7 Figure 7.21: 6:1 Prolate Spheroid: α = 5 • . Skin friction lines as computed by γ + CF model seen from different points of view, leeward (LW), top and windward (WW) sides.

7. 2 . 3

 23 Conclusions of the NATO-AVT/313 Workshop for α = 5 • and α = 15 •

  . The numerical predictions by γ + CF are shown in Fig.(7.23), the skin friction C f and the pressure distribution C p . The new crossflow criterion performs very well on the windward side. Transition predictions are in good agreement with experiments. Despite the lower skin friction contours predicted by the transition model, the quantitative deviation between measured and predicted C f is justified by the experimental uncertainties, as shown in Fig.(7.24). Despite the good performance at the windward side, γ + CF model fails to predict the correct transition features on the leeward side. The measured laminar separation within the flow undergoes transition right after the nose of the spheroid is not reproduced in the simulations. The laminar separation is clearly visible by the girthwise distribution of C f . The positive peak of C f at φ ∼ 150 • in Fig.(7.24a) and Fig.(7.24b) represent the transition location within the bubble.

  (

  7.24c) to Fig.(7.24f). Two different converging lines are observed on the surface, one close to the middle of the spheroid, the other on the leeward side, as shown in Fig.(7.25

Figure 7 .

 7 Figure 7.22: 6:1 Prolate Spheroid: α = 30 • . Measured skin friction C f distribution in the X/L -φ plane.

  (a) X/L -φ plane: C f distribution. (b) X/L -φ plane: Cp distribution.

Figure 7 .

 7 Figure 7.23: 6:1 Prolate Spheroid: α = 30 • . Predicted skin friction C f and pressure C p coefficients distributions on the 6:1 prolate spheroid computed using γ + CF model.

  (a) X/L = 0.053. (b) X/L = 0.139. (c) X/L = 0.395. (d) X/L = 0.565. (e) X/L = 0.738. (f) X/L = 0.936.

Figure 7 .

 7 Figure 7.24: 6:1 Prolate Spheroid: α = 30 • . Experimental and numerical girthwise distribution of the C f coefficient at different sections. Numerical results are computed with γ + CF model. At each validation point, the experimental uncertainty is reported.

Figure 7 .

 7 Figure 7.25: 6:1 Prolate Spheroid: α = 30 • . Skin friction lines as computed by γ + CF model seen from different points of view, leeward (LW), top and windward (WW) sides.

  et al.and the MUB for Petzold et al.experiments. The free-stream conditions are taken from Petzold et al., who reported an average turbulence intensity T u = 0.17%. In Kruse et al. experiments, it is mentioned that the average free-stream disturbance level is lower than in Petzold et al.. In both cases, we observe a noticeable extension of the laminar region.Experimental Set Up A sketch of the geometry under study is shown in Fig.(7.26). The wing consists of five sections. The unswept peniche raises the model above the tunnel wall, in order to avoid the contamination of the laminar boundary layer by the turbulent boundary layer developed at the tunnel wall. The other three swept sections, A, B, C are dedicated to transition measurements. Each segment have the same thickness and span length, but with increasing sweep 30 • , 45 • , 55 • , respectively. The wing in completed by the tip, positioned at its outer extremity. Because of this particular design, the sickle wing experiences crossflow from the root towards its tip. The predominance of stationary crossflow instabilities over T-S waves in the transition process depends on the Reynolds number and angle of incidence. For the current simulations, we are going to refer to the flow around the wing operating near its design point, Re = 2.75 × 10 6 and α = -2.6 • . For these conditions, T-S and CF waves occur at the same time. Experimental data include pressure tap measurements approximately along the midsection of each segment. Transition location is detected by means of infrared thermography images and it is based on the increase of the convective heat flux due to turbulent mixing. Its position is given by the local minimum of temperature gradients in the measurements from Kruse et al.. In the paper presented by Petzold et al., the transition location is defined as the middle point between the x-coordinate value of the position of temperature maximum, before the temperature drop, and the point of minimum temperature right after the temperature drop. Each section presents a peculiar transition pattern. On the upper surface, transition on section A is dominated by T-S waves, as suggested by the measured straight transition line. Moving upstream in the spanwise direction, the transition process becomes to be CF-dominated. The zig-zag pattern, visible in the infrared images of the transition location (not shown here), on section B and the first half of section C, indicates that transition is dominated by stationary crossflow. Transition on the lower surface occurs under a moderate adverse pressure gradient and it is T-S dominated. In the experimental report by Petzold et al., the reported surface roughness is about 1.47 µm, however it is neglected in our the numerical simulations.

Figure 7 . 26 :

 726 Figure 7.26: Sickle Wing. Sketch of the specimen, all the dimensions are in mm.Figure taken from Kruse et al., [153].

  Figure 7.26: Sickle Wing. Sketch of the specimen, all the dimensions are in mm.Figure taken from Kruse et al., [153].

5 :

 5 Sickle Wing: α = -2.6 • . Mesh details.(a) Total computational domain. (b) Mesh at the surface.

Figure 7 . 27 :

 727 Figure 7.27: Sickle Wing: α = -2.6 • . Grid used for the flow simulation around the sickle wing. Total computational domain and mesh at the surface. In Fig.(7.27b), the solid red lines represent the kinks, where the geometrical sweep angle changes.

  (a) γ: C f on the upper surface. (b) γ: C f on the lower surface. (c) γ + CF : C f on the upper surface. (d) γ + CF : C f on the lower surface.

Figure 7 . 28 :

 728 Figure 7.28: Sickle wing: α = -2.6 • . Predicted skin friction contours on the upper and lower surface of the sickle wing by γ and γ + CF models. The black dotted line is the extracted experimental transition location from Petzold et al., the black solid line is from Kruse et al..

Figure 7 . 29 :

 729 Figure 7.29: Sickle wing. Computational mesh used from Kim et al., for their transitional simulation around the sickle wing. Figure reproduced by Kim et al., [13].

Figure 7 . 30 :

 730 Figure 7.30: Sickle wing: α = -2.6 • . Comparison of measured and calculated pressure coefficient distributions at each sweep section of the model, from Section A to Section C. Pressure is extracted along the dotted lines (pressure taps) represented in Fig.(7.26). Numerical results are computed by γ and γ + CF transition variants.

Titre:

  Modélisation de la Transition du Régime Laminaire au Régime Turbulent d'un Ecoulement Incompressible dans un Environnement RANS pour des Configurations 2D et 3D Mots clés : RANS, LCTM, Modes Transversaux, Transition Résumé : Ce travail de thèse vise à évaluer la faisabilité des calculs de transition dans le cadre des codes industriels. Les points forts ainsi que les faiblesses des modèles de transition basés sur des corrélations locales, γ et γ-Re_θ, sont analysés sur des configurations bi et tridimensionnelles, en se concentrant sur la modélisation physique et les aspects numériques. L'objectif est d'analyser le potentiel de ces modèles RANS en tant qu'outils prédictifs, capables de gérer automatiquement et de manière autonome la transition d'un écoulement du régime laminaire au régime turbulent. Nous évaluons les performances des modèles de transition sur un grand nombre de cas de test, couvrant un large éventail de mécanismes de transition. Dans le cas des configurations 3D, un point crucial est la modélisation de la transition due aux modes transversaux et stationnaires, qui sont les principaux mécanismes de transition tridimensionnel dans un environnement à faible turbulence. A cet effet, nous présentons dans cette thèse un recalibrage original du critère transversal Tc1, proposé à l'origine par Menter et Smirnov en 2014 et basé sur une formulation locale du célèbre critère transversal C1 de Daniel Arnal. Ce critère recalibré étend notamment le modèle γ existant pour la prédiction de la transition transversal autour des géométries complexes. Title : Laminar-to-Turbulence Transition Modeling of Incompressible Flows in a RANS Framework for 2D and 3D Configurations Keywords : RANS, LCTM, Crossflow, Transition

  Reynolds number and streamwise wavenumber, diagram for Blasius boundary layer flow: contours of constant growth rate c i . The red dot corresponds to the critical Reynolds number, Re c = 519.4, and streamwise wavenumber α c = 0.303. They lie on the neutral curve c i = 0 that divides the unstable (shaded) and stable region. Figure reproduced from Schmid, [23]. . . . . . . . . . . . . . . . . . . . . Development of velocity profiles on a swept wing when suffering from an inflection of the pressure gradient. W 1 (y) and U 1 (y) are the crossflow and main-flow profile projected onto a coordinate system relative to the external inviscid streamline. β is the angle between the wall shear stress and the external streamlines. Figure re-adapted from Yiming et al., [38]. . . . . . . . . . . . . . . . . . . . . . . . . . . Visualization of small length vortices at different time steps through λ 2 isosurfaces: view from the top in Fig.(2.11a), view in the (y, z) plane for fixed x in Fig.(2.11b) with the correspondent velocity perturbations contours. Simulations are performed at Re = 1000 on a grid of 1920 × 241 × 128 points. Figures are reproduced from Liu et al., [51]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.12 N T S and N CF stability diagram for the 6:1 prolate spheroid. Figure is reproduced from Stock, [58]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.13 Scaled vorticity Reynolds number Re V for a Blasius boundary layer, as computed using Falkner-Skan equations. Re θ is the corresponding momentum thickness Reynolds number and 2.2 is the calibration constant set to have a maximum of one in the boundary layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Evolution of the vorticity Reynolds number to the momentum thickness Reynolds number vs the shape factor H 12 as presented in Menter & Langtry , [74], (black squares), for γ -Re θ , and as computed resolving Falkner-Skan equations (red line).

	xiv xviii		NOMENCLATURE LIST OF FIGURES
	Re δ 2 Crossflow Reynolds Number 2.11
	Re θ	Momentum Thickness Reynolds Number
	ρ	Fluid Density	[kg/m 3 ]
	Nomenclature θ Momentum Thickness C d Drag Coefficient List of Figures	[m]
	C f	Skin Friction Coefficient, τw 1 2 ρU 2 ∞
	C l	Lift Coefficient
	α β δ δ * γ λ θ µ µ t ν ν t Ω ω Ω ij Re θt Transported Momentum Thickness Reynolds Number Angle of Attack Hartree Parameter Boundary Layer Thickness Displacement Thickness Intermittency Pressure Gradient Parameter, (ρθ 2 /µ)(dU/ds) Fluid Dynamic Viscosity Dynamic Eddy Viscosity Fluid Kinematic Viscosity Kinematic Eddy Viscosity Magnitude of the Mean Vorticity Tensor, 2Ω ij Ω ij Specif Turbulence Dissipation Rate Mean Vorticity Tensor, R Sweep Reynolds Number, W e / (νS) -→ n Surface Normal Vector Φ Geometrical Sweep Angle φ Azimuthal Angle of the 6:1 Prolate Spheroid φ Sweep Angle φ L Local Sweep Angle Approximation Ψ Crossflow Strength Indicator in T c1 Formulation Re Reynolds Number xiii C p Pressure Coefficient, (p-p∞) 1 2 ρU 2 S ij Mean Strain Rate Tensor, ∂U i /∂x j +∂U j /∂x i 2 T u Turbulence Intensity, 2k 3U 2 u, v, w Velocity Components, x, y, z Reference Coordinate System y Wall Normal Distance c Chord 2.1 (Re, α), xvii	[ • ] [m] [m] [1/s] [kg/(m s)] [kg/(m s)] [m 2 /s] [m 2 /s] [1/s] [1/s] [ • ] [ • ] [ • ] [ • ] [1/s] [m/s] [m]

∂U i /∂x j -∂U j /∂x i 2 [1/s] ∞ C d,

f Friction Drag Coefficient C d,p Pressure Drag Coefficient dv/dy Wall-normal Derivative of the Wall-normal Velocity Component, ∇( -→ u • -→ n ) • -→ n h i Typical Cell Size of Grid i H 12 Shape Factor in Streamwise Direction He Helicity, -→ u • (∇ × -→ u ) [m/s 2 ] k Turbulence Kinetic Energy [m 2 /s 2 ] R T Eddy Viscosity Ratio, ν t /ν Re V Vorticity Reynolds Number, ρy 2 Ω/µ S Magnitude of the Mean Strain Rate Tensor, 2S ij S ij [1/s] 2.2 Smoke visualization of streamwise streaks appearing in a boundary layer subjected to moderate level of freestream turbulence intensity T u = 2.2%. The flow, coming from the left, is captured at a certain distance from the leading edge and it transitions (breakdown and turbulence spot formations) in the right edge of the image. Figure reproduced from Matsubara et al., [26]. . . . . . . . . . . . . . . . 2.3 Different instability mechanisms on the swept wing. Figure reproduced from Shahriari, [27]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Schematic representation of the transition process due to T-S waves. Figure is reproduced from White & Corfield, [28]. . . . . . . . . . . . . . . . . . . . . . . 2.5 Görtler vortices at near a stagnation point (S) lying on a convex surface. R is the curvature radius of the external streamline. Figure reproduced from Hirschel et al., [31]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Schematic representation of the flow near the leading edge of a swept wing, reproduced from Poll, [36]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Infinite swept wing notations. Figure re-adapted from Arnal et al., [16]. . . . . . 2.8 2.9 Schematic representation of a laminar separation bubble: streamlines and velocity profile (left) from O'Meara & Muller, [45] and the experimental pressure coefficient distribution (right) from Lee et al., [46]. . . . . . . . . . . . . . . . . . 2.10 Maximum growth of laminar fluctuations in a laminar boundary before transition for a flow on a flat plate. Figure reproduced from Mayle & Schulz , [47]. . . . . .

  • . Skin friction C f distribution on the 6:1 prolate spheroid computed using γ model without crossflow inclusion on Grid5 and Grid5bis of Table(7.2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 7.6 6:1 Prolate Spheroid: α = 15 • . Convergence of pressure and viscous drag for γ + CF formulation with grid refinement. p is the estimated convergence order Prolate Spheroid: α = 30 • . Predicted skin friction C f and pressure C p coefficients distributions on the 6:1 prolate spheroid computed using γ + CF model. . 126 7.24 6:1 Prolate Spheroid: α = 30 • . Experimental and numerical girthwise distribution of the C f coefficient at different sections. Numerical results are computed with γ + CF model. At each validation point, the experimental uncertainty is reported.127 7.25 6:1 Prolate Spheroid: α = 30 • . Skin friction lines as computed by γ + CF model seen from different points of view, leeward (LW), top and windward (WW) sides. 128 7.26 Sickle Wing. Sketch of the specimen, all the dimensions are in mm.

of the discretization error, as described in Sec.(

5

.2.1). . . . . . . . . . . . . . . . . 112 7.7 6:1 Prolate Spheroid: α = 15 • . Convergence of pressure and viscous drag for γ -Re θ + CF formulation with grid refinement. p is the estimated convergence order of the discretization error, as described in Sec.(5.2.1). . . . . . . . . . . . . 112 7.23 6:1

  • . Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.13 NACA0015. Details of the five structured C-O Grids: total number of cells, N cells , number of cells on the flat plate surface, N cells-plate , grid refinement ratio, r i = h i /h, and y + max values. . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.14 NACA0015: α = 5 • . Numerical efforts predicted using γ, γ -Re θ and k -ω SST (2003) models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7.1 3D Test cases details: geometry, angle of attack α and Reynolds number. The transition model used for the computations is provided in the last two columns. 106 7.2 6:1 Prolate Spheroid: Mesh details. . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.3 6:1 Prolate Spheroid: α = 5 • , 15 • , 30 • . Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

	Chapter 1

  Evolution of the vorticity Reynolds number to the momentum thickness Reynolds number vs the shape factor H 12 as presented in Menter & Langtry ,[START_REF] Menter | Transition modelling for turbomachinery flows[END_REF], (black squares), for γ -Re θ , and as computed resolving Falkner-Skan equations (red line).

								3.3. LCTM
	1.8
					Langtry&Menter
	1.6				FSK
	1.4
	θ )
	Re V /(2.2Re 1.2 1
	0.8
	0.6
	2.2 0.4	2.4	2.6	2.8	3	3.2	3.4	3.6
				H 12
	Figure 3.1:

Table 5 .

 5 1: 2D Test cases details: geometry, angle of attack α and Reynolds number. The turbulence intensity at the inlet, T u in , and the transition mechanisms are reported as well.

		α	Re	γ γ -Re θ	T u in	Transition Mechanism
	Flat Plate-T3A-		10 7		0.0100135	Natural
	Flat Plate-T3A		10 7		0.0536609	Bypass
	NACA0015	5 • 1.8 × 10 5		0.005	Separation-Induced
	Eppler 387	1 • 7 •	3 × 10 5		0.01	Separation-Induced Natural

Table 5 . 2

 52

	y + max

: Flat Plate. Details of the five structured O-grids: total number of cells, N cells , number of cells on the flat plate surface, N cells-plate , grid refinement ratio, r i = h i /h, and y + max values.

Table 5 .

 5

		Laminar Region Transition Region Turbulent Region Total
	E mv /E ref (U input = 0)	2.36	2.23	0.76	2.06
	E mv /E ref (U input = 0)	2.35	1.97	0.76	1.95

4: Flat Plate: T3A. Multivariate metric for the skin friction in different regions and all along the flat plate for γ model results.

Table 5 .

 5 5: Flat Plate: T3A. Multivariate metric for the skin friction in different regions and all along the flat plate for γ -Re θ model results.Velocity profiles at three different sections are presented in Fig.(5.5). The multivariate metric results are given inTable(5.6) and (5.7), for γ and γ -Re θ , respectively. As the skin friction, the multivariate metric results in the transition region is affected by the inclusion of the input uncertainties. In particular, E mv /E ref at x/L=0.02035 for γ -Re θ considerably reduces if input errors are taken into account. This behavior is related to the other source of dependence on the eddy viscosity ratio within γ -Re θ formulation, namely, the value of the empirical Re θt diffused within the boundary layer. Being Re θt a function of the free-stream turbulence intensity, any difference of the latter in the free-stream can accelerate or decelerate the transition process.

		X/L = 0.01006 X/L = 0.02035 X/L = 0.05237
	E/E ref (U input = 0)	1.61	1.36	1.76
	E/E ref (U input = 0)	1.61	1.32	1.76

Table 5 .

 5

		1.57	1.89	1.58
	E/E ref (U input = 0)	1.57	1.62	1.58

6: Flat Plate: T3A. Multivariate metric for the streamwise velocity at different sections of the flat plate for γ model results. X/L = 0.01006 X/L = 0.02035 X/L = 0.05237 E/E ref (U input = 0)

Table 5 . 7

 57

: Flat Plate: T3A. Multivariate metric for the streamwise velocity at different sections of the flat plate for γ -Re θ model results.

Table 5 .

 5

	5.5. EPPLER387

10: Eppler387. Details of the five structured C-O-grids: total number of cells, N cells , number of cells on the flat plate surface, N cells-plate , grid refinement ratio, r i = h i /h, and y + max values.

Table 5 .

 5

		d C d,p C d,f	C l
	EXP	9.3			0.465
	γ	8.66 3.2 5.46 0.484
	γ -Re θ	9.3	4	5.3 0.486
	k -ω SST 12.6	10	2.6 0.474

11: Eppler387: α = 1 • . Experimental and numerical efforts predicted using γ and γ -Re θ and k -ω SST (2003).

Some results from the AVT-313 Workshop on the Eppler387 profile at α = 1 •

Table 5 .

 5

	h 1	y + max

13: NACA0015. Details of the five structured C-O Grids: total number of cells, N cells , number of cells on the flat plate surface, N cells-plate , grid refinement ratio, r i = h i /h, and y + max values.

Table 5 .

 5 14: NACA0015: α = 5 • . Numerical efforts predicted using γ, γ -Re θ and k -ω SST (2003) models.

	Re θ	15.825228 9.344608	6.48062	0.658284
	γ	15.84759	9.558022	6.289568 0.6534532
	k -ω SST	16.6643	6.399988	10.264318	0.4987

  , Arnal et al. proposed a new analytical expression for Re δ 2 t , formulated as:

	Re δ 2 t =	 	300 π arctan	0.106
		

(H 12 -2.3)

2.052 , for 2.3 < H 12 2.7, 150, for H 12 2.3.

Table 7 .

 7 • . Coarsest mesh in the wall normal direction at the section X/L = 0. N cells N surface h i /h 1 y + 2: 6:1 Prolate Spheroid: Mesh details.

	max

Table 7 .

 7 3: 6:1 Prolate Spheroid: α = 5 • , 15 • , 30 • . Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%).

Table 7 .

 7 Table(7.6). 6: Sickle Wing: α = -2.6 • . Inflow turbulence intensity T u in (%) and eddy viscosity ratio ν t in /ν, and value of turbulence intensity in the proximity of the leading edge T u(%).

	Re	T u in (%) ν t in /ν T u(%)
	2.75 × 10 6	0.20	2.24	0.17

The PhD dissertation has not been published yet, some of the results can be found in Kim et al.,[START_REF] Kim | Boundary layer transition models for naval applications: capabilities and limitations[END_REF]

Details of the manipulation to obtain the equation for the normal velocity can be found in Schmid,[START_REF] Schmid | Stability and Transition in Shear Flows[END_REF].

Derivation of both criteria is omitted and can be found in Schmid,[START_REF] Schmid | Stability and Transition in Shear Flows[END_REF].

Following the ASME V&V

20-2009 standard, [104].

Results without CF using γ -Re θ are not shown in the dissertation, but they are very similar to γ.

As reported in Perraud et al.,[START_REF] Perraud | Overview of transition prediction tools in the elsA software[END_REF].

This information was shared within an AVT-313 meeting.

Acknowledgements

First of all, I would like to thank my supervisor Michel Visonneau and all the METHRIC team for their help, guidance and support during these years. I will carry this experience with me throughout my future career.

Their work is the result of the cooperative effort of CNR-INM and DLR. The majority of the computations were performed using HPC resources from GENCI (Grand Equipement National de Calcul Intensif) (Grant-A0102A00129), which is gratefully acknowledged.

Algorithm

The discretization of mass and momentum conservation yields to a set of algebraic equations: one for each volume and for each transport/conservation equations. These non-linear and coupled equations are solved by a segregated algorithm, following the steps reported below:

1 Flow field initialization Q 0 at t = t 0 ; 2 New time step t = t + ∆t;

Transition equations

All the simulations presented in the next chapters are for single fluid steady flow. The time derivative is zeroed ∂ ∂t = 0 , nevertheless the fictitious time derivative ∂ ∂τ is kept to increase the diagonal dominance of the the linear solver. This is used for the resolution of transition, turbulence and velocity components. About 30 Gauss-Seidel iterations are needed to converge, with a diagonal dominance increased by 50%. The pressure system is resolved using an algebraic multigrid technique, BoomerAMG. This approach reduces the solving time compared to the Krylov subspace solvers-like PGMRES (Preconditioned Generalized Minimal Residual) algorithm. Convective fluxes of transition, turbulence and momentum equations are discretized using AVLSMART scheme. On parallel machines, the computational domain is split into multiple connected domains with approximately the same number of unknowns, using the M TiS partitioning algorithm, by Karipis et al., [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF], [START_REF] Tis | Family of multilevel partitioning algoritms[END_REF]. Communication of data at the faces between domains is performed according to the Message Passing Interface standard (MPI), Clarke et al., [START_REF] Clarke | The MPI Message Passing Interface Standard[END_REF].

Alternative Approach to Arnal's C1 criterion

In the last six years, the research of a criterion for stationary crossflow prediction easily applicable to arbitrarily shaped 3D configurations and not based on the FSC profiles has led to the proposition of the Helicity Reynolds number Re He CF extension. Re He is defined as:

where the helicity He,

is used as crossflow strength indicator. Indeed, He is the streamwise vorticity within the boundary layer. Crossflow instabilities result in a system of co-rotating vortices in the boundary layer, physically and He can be seen as a measure of the flow's 'twist'. Its use was first proposed by Muëller and Herbst, in [START_REF] Müller | Modelling of crossflow-induced transition based on local variables[END_REF], and again by Langtry et al., [START_REF] Langtry | Extending the γ-re θ correlation based transition model for crossflow effects[END_REF]. Both Muëller and Langtry CF extensions involved the modification of the production term of Re θt : the production P θt in Eq.(3.28) is lowered by a sink term, that activates when the crossflow criterion based on the helicity Reynolds number is met. In Langtry model, the crossflow criterion also accounts for surface roughness effects. Also Grabe et al., [START_REF] Grabe | Transport modeling for the prediction of crossflow transition[END_REF], proposed an helicity-based crossflow criterion, that, differently from the variants mentioned above, envisages the modification of the function F onset in the transport equation for the intermittency γ. These helicity-based crossflow criteria were all calibrated on numerical database constructed considering the experimental results on the ONERA D airfoil from Schmitt et al., [START_REF] Schmitt | Ecoulements subsoniques et transsoniques sur une aile en flèche variable[END_REF], and in ONERA/CERT/DERAT, [126], results on the infinite-swept NLF (2)-0415 wing from Dagenhart et al., [START_REF] Dagenhart | Crossflow stability and transition experiments in swept-wing flow[END_REF], but also on the 3D 6:1 Prolate Spheroid from Kreplin, [START_REF] Kreplin | Wall shear stress measurements on an inclined prolate spheroid in the DFVLR 3m × 3m low speed wind tunnel[END_REF].

We will use the helicity criterion formulated by Grabe as CF extension for γ -Re θ computations.

The detailed mathematical formulation of the criterion is given below.

Helicity criterion: Crossflow Inclusion in γ -R θ model

The critical crosswise Reynolds number ratio reads as

where the correlation Re + He,t is numerically determined. It is curve fitted through the numerical data obtained for different configurations at several operating conditions, ONERA D profile, NLF(2)-0415, and the 6:1 prolate spheroid. It is expressed as function of the shape factor H 12 , as follows:

Re + He,t = max(-456.83H 12 + 1332.7, 150). (6.58)

The shape factor H 12 is approximated through the pressure gradient parameter λ θ , that is defined as:

y. (6.60)

CONCLUSIONS ON THE T C1 NEW CALIBRATION

The T c1 criterion has been calibrated using the FSC solutions, f , g, the shape factor H 12 and the Hartree parameter β. We are aware of the stringent assumptions of FSC equations, i.e. zero-pressure gradient in the spanwise directions and zero crosswise velocity component for zeropressure gradient (β = 0) boundary layers. Nevertheless, in lack of detailed experimental results, we are convinced that the calibration on Falkner-Skan-Cooke solutions is a robust approach. It does not depend on numerical discretization errors nor input uncertainties, especially given the dependence of the model on the inlet eddy viscosity ratio and meshes, see Chapter 5.

Chapter 7

3D Configurations: Results

Introduction

The present chapter is dedicated to the analysis of the performance of the local correlationbased transition models for three-dimensional boundary layers. We will present a comparison between simulations and experimental data for two different geometrical configurations: the 6:1 prolate spheroid at three different angles of attack, α = 5 • , 15 • , 30 • , and Re = 6.5 × 10 6 , and the sickle wing at α = -2.6 • , and Re = 2.75 × 10 6 . The test cases are schematically summarized in Table (7.1). The last two columns refer to the transition model chosen for the numerical simulations. Clearly, our attention will mainly focus on the validation of the T c1 crossflow criterion, as re-calibrated in Chapter 6, and coupled to γ transition model (γ + CF ). Nevertheless, we will also analyze the numerical results around the 6:1 prolate spheroid at the two lowest angles of attack computed by γ -Re θ . As crossflow variant for γ -Re θ we have opted for the helicity crossflow criterion (γ -Re θ + CF ) presented in Sec.(6.5). We are going to start from the results around the 6:1 prolate spheroid. We will first consider the two lowest angles of attack. We will compare both predictions by γ + CF and γ -Re θ + CF models with the experimental measurements from Kreplin, [START_REF] Kreplin | Wall shear stress measurements on an inclined prolate spheroid in the DFVLR 3m × 3m low speed wind tunnel[END_REF]. These results have been also presented at the AVT-313 workshop dedicated to transition modeling for 3D boundary layers, held virtually in March 2021. A section will be devoted to the workshop discussion. The rest of the chapter is dedicated to the numerical prediction by γ + CF . We will discuss the flow around the 6:1 prolate spheroid at α = 30 • and we will conclude with the flow around the sickle wing at incidence.

Geometry

Approaching the trailing edge of the body, none of the models is able to predict transition for φ < 30 • , and the flow close to the windward symmetry plane remains laminar until the trailing edge. Measurements and predictions start deviating at the section X/L = 0.652, Fig. 7.14a, up to the last section X/L = 0.936, Fig. 7.14b. None of the crossflow criteria, the T c1 nor the He is active in this region of the prolate spheroid. Nevertheless, at least at the windward symmetry plane, transition is hardly due to pure crossflow. This incorrect flow feature at the windward side is related to the fact that both of the formulations, γ + CF and γ -Re θ + CF , do not account for any transition mechanisms other than T-S and CF instabilities. The possible missing physical mechanisms are discussed later in the chapter. Contrary to the C f , the pressure distribution C p is not affected by the chosen transition model, nor by the crossflow inclusion (results not shown). The skin friction lines, as computed by the γ + CF model are plotted in Fig. (7.15). Results for γ -Re θ + CF are very similar and therefore omitted. Both models, indeed, predict an envelope of converging wall streamlines, which are representative of an open-separation. In this case, the separation line cannot be uniquely identified by skin friction minima, according to Simpson, [START_REF] Simpson | Aspects of turbulent boundary-layer separation[END_REF]. Indeed, skin friction lines are converging from both the windward (WW) and leeward sides (LW). The angle of incidence is mild, and the boundary layer remains attached until nearly the trailing edge.

Given the better results obtained by the recalibrated T c1 criterion as included in γ, in the following discussion the results from the original criterion by Menter & Smirnov will not be considered for the next comparisons at other angles of attack. The computations for α = 5 • , 30 • are computed on Grid3 of Table (7.2). Indeed, considering the predicted transition onset, the computations on this grid are sufficiently accurate, as mentioned before.

Appendix A

Evaluation of the order of convergence p and the numerical uncertainty

When dealing with a set of noisy data, the determination of the order of convergence p is extremely sensitive to any small perturbation. The scatter in the data is due to different factors, namely the lack of geometry similarity of the tested grids and the use of switches and damping functions in the turbulence/transition models. We cannot rely only on the discretization error δ RE , as defined in Eq.(5.13), which gives an error estimate when the convergence is monotonic, but alternative error estimators need to be introduced. In terms of least squares fit, the complete set of functions used for the estimations of φ 0 , p and α are:

and

The weights w i are based on the typical cell size and, given the number of grids N g , they are defined as follows:

where

The non-weighted functions are obtained imposing w i = 1 and they are referred to as S RE , S 1 , S 2 and S 1,2 . The introduction of weighted functions is due to the need to give more emphasis to the finest of the grids in practical applications.

The procedure for the error estimation is given below, step by step:

1. First, the observed order of convergence p is determined by minimizing the functions S RE and S w RE . If both solutions lead to p > 0, the best fit is chosen as the one with the smallest standard deviation σ RE , that is hereafter repeated:

If any of the solutions yields to a negative p, it is discarded. If both of them give a negative p, the data are classified as anomalous. An anomalous behavior represents either monotonic divergence, or oscillatory convergence, i.e. the solution is converging above and below the exact solution, or non-monotonic convergence.

2. Once the best fit has been chosen, the order of convergence p, the constant α and the estimate of the exact solution φ 0 are derived. If p is in an acceptable range, i.e. 0.95 ≤ p < 2.1, the numerical uncertainty is estimated from δ RE .

3. If the estimated p is bigger than 2.1, then four fits are determined by solving the equations S w 1 , S 1 , S w 2 and S 2 , because such an high p would lead to an extremely small δ RE . The best fit is chosen with respect to the smallest standard deviation.

4. If p is smaller than 0.95, than δ RE would produce an error estimate too conservative, i.e.

p → 0 leads to δ RE → ∞. To the four equations S w 1 , S 1 , S w 2 and S 2 , the equations S w 1,2 , S 1,2 are added in order to avoid to be under-conservative. 5. For anomalous data, the best fit is obtained among the fits S w 1 , S 1 , S w 2 , S 2 , S w 1,2 and S 1,2 with the smallest standard deviation. The estimation of the numerical uncertainty depends on the error estimators φ used for the fit. Reminding the definition of the data range parameter

the uncertainty is defined with respect to the apparent convergence condition determined by the least square fit:

• If the fit is reliable, i.e. 0.5 ≤ p < 2.1 and σ ≤ ∆ φ , the safety factor F s , following the Grid Convergence Index (GCI) procedure, is equal to 1.25 and

• If σ > ∆ φ , then F s = 3 and

where φ f it is the value obtained from the fit for the same grid density.

Appendix B Boundary Layer Thickness

The boundary layer thickness is identified using the method proposed by Griffin et al., [START_REF] Griffin | General method for determining the boundary layer thickness in nonequilibrium flows[END_REF]. The definition for the boundary layer thickness reads as:

where U I (y) is the local reconstruction of the inviscid mean streamwise velocity profile. δ n indicates the wall-normal distance y at which n% on U I (y) is attained, for the boundary layer edge it is taken as y = δ 99 . The reconstruction of the inviscid solution U I (y) makes use of the stagnation pressure P o , defined for an incompressible steady flow as:

P is the static pressure and U m = U 2 I +V 2 , with V mean wall normal velocity. The reconstructed streamwise velocity U I (y) is obtained using Bernoulli's equation. Indeed, Bernoulli's equation can be applied locally, considering an hypothetical irrotational flow, whose V (y) and P (y) have an equivalent profile as the one they would in a viscous flow. U I (y) reads as:

where P o,ref is the total pressure defined as P o,ref = max(P o ).