G Rubino

M Visonneau

Laminar-to-Turbulence Transition Modeling around the 6:1 Prolate Spheroid at Different Angles of Attack

, is a local approximation of the well-known C1 criterion by Daniel Arnal (1984). The C1 criterion is evaluated locally through auxiliary functions expressed by the solution of the Falkner-Skan-Cooke (FSC) equations. In the following, an original recalibration of the T c1 is presented. A local approximation of the sweep angle is included in the criterion formulation in order to achieve better results on non-wing-like geometries. The performance of both γ and γ -Re θ transition models, with their respective crossflow transition criterion, is therefore discussed for the 6:1 prolate spheroid highlighting strengths and weaknesses of the transition models current formulations.

INTRODUCTION

Understanding, predicting and controlling laminar to turbulence transition is nowadays a main challenge in Computational Fluid Dynamics (CFD). Indeed, there are several practical applications that deal with low/moderate Reynolds numbers, such as aerial and marine unmanned vehicles, small submarines, but also wind turbines. For these applications, transitional effects are important, if not dominant, and Reynolds averaged Navier-Stokes (RANS) turbulence models fail to predict transitional flow features. In the last decade, the local-correlation transition model (LCTM) concept by Menter & Langtry has known a huge success. The RANS transition models γ -Re θ , Menter and Langtry (2012), and γ, [START_REF] Menter | A one-equation local correlation-based transition model[END_REF], are today widely used. The latter, γ, was proposed by Menter et al. as a drastic simplification of γ -Re θ : the transport equations for the transition variables are reduced from two to one, because only local quantities are used in the empirical correlations, and the model formulation is Galilean invariant.

RANS transition models are widely preferred in practical applications to Direct Numerical Simulations (DNS) and Large Eddy Simulations (LES) approaches, because of the reduced computational costs. Nevertheless, RANS approaches are not able to provide the real transition physics. As a matter of fact, the aim of these transition models is not to describe transition physics, but rather to identify the different regions of the flow and to predict accurately the transition location. All the physics is contained in empirical correlations, that account for a specific transition mechanism. As initially presented, γ and γ -Re θ models only accounted for bypass and streamwise natural transition due to Tollmien-Schlichting (T-S) waves. In the past years, several approaches have been proposed by different research groups to account for transition due to stationary crossflow (CF) waves. These instabilities occur in three-dimensional boundary layers and in low free-stream turbulence environments. So far, the helicity He based criterion, proposed by [START_REF] Grabe | Transport Modeling for the Prediction of Crossflow Transition[END_REF] as an extension of γ -Re θ , is the one that has known the biggest success. It uses the helicity as indicator of the crossflow strength. Indeed, crossflow instabilities produce a system of co-rotating vortices within the boundary layer, aligned with the local velocity vector. The calculation of the helicity uses the local velocity vector, thus the criterion is not Galilean invariant. For this reason, in 2014, Menter & Smirnov proposed the T c1 criterion, a crossflow criterion that uses the derivative of the normalised vorticity in the wall-normal direction as indicator of the crossflow strength. The T c1 is the local reconstruction of the renowned C1 criterion proposed by [START_REF] Arnal | Théorie de l'instabilité laminaire et critères de transition en écoulement bi et tridimensionnel[END_REF] and based on the crossflow Reynolds number Re δ 2 :

Re δ 2 = U 1 e δ 2 ν where δ 2 = - ∞ 0 W 1 U 1 e dy. (1)
U 1 and W 1 are the streamwise and crosswise velocity components. U 1 e = (U 2 e +W 2 e ) is the velocity at the edge of the boundary layer. The C1 criterion reads as:

Re δ 2 Re δ 2t = 1, (2)
where the Reynolds number value at the transition location, Re δ 2t , is not unique, but it depends on the longitudinal shape factor H 12 . In a three-dimensional boundary layer, H 12 is defined as:

H 12 = δ * θ , (3)
δ * = δ 0 1 - U 1 U 1 e dy, (4)
θ = δ 0 1 - U 1 U 1 e U 1 U 1 e dy, (5)
where δ * and θ are the displacement and momentum thickness, respectively. Finally, the crossflow Reynolds number at transition onset is given by the expression:

Re δ 2t = 300 π arctan 0.106 (H 12 -2.3) 2.052 , for 2.3 < H 12 2.7, 150, for H 12 2.3.

(6) The T c1 criterion, as formulated by Menter & Smirnov, does not use explicitly the velocity vector. Therefore it preserves the Galilean-invariant formulation of γ transition model.

In the following discussion, the transition models γ -Re θ and γ mathematical equations are presented, as well as the respective crossflow transition criteria, the one based on He and the T c1. Attention is focused especially on the latter. First, the T c1 formulation proposed by Menter & Smirnov is discussed, highlighting its limitation. Thus, the new original re-calibration of the T c1 is proposed, that extends the criterion application to more complex 3D geometries.

The transition models performance is then tested on the 6:1 prolate spheroid. This geometry serves as simple surrogate for axisymmetric bodies, such 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, both experimentally and numerically, test case, because it exhibits all the complex physics associated with crossflow separation. The complexity of crossflow separation stands in its character: it does not originate from a unique singularity point nor line, and it does not strongly interact with the local flow field. Indeed, the flow is nearly attached to the surface ahead separation. Crossflow separation can then be identified as a boundary layer rolling around itself and it is usually referred to as open separation. It is characterized by an envelope of converging streamlines, following the definition in Surana et al. (2006), used also in Wetzel et al. (1998). An extensive review on the 6:1 prolate spheroid can be found in [START_REF] Andersson | Instabilities in the wake of an inclined prolate spheroid[END_REF] and [START_REF] Fu | Challenges in Naval Ship and Submarine Hydrodynamics[END_REF]. Experimental measurements can be found in Fu et al. (1994) and Wetzel and Simpson (1998). These experiments, performed in the 1990s, were used by different groups for the assessment of (U)RANS, DES and LES, [START_REF] Kim | Application of modern turbulence models to vortical flow around a prolate spheroid[END_REF], [START_REF] Wikström | Large eddy simulation of the flow around an inclined prolate spheroid[END_REF], [START_REF] Fureby | LES of the Flow Past a 6: 1 Prolate Spheroid[END_REF], [START_REF] Fureby | 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, the interest lies on 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, forcing turbulence close to the leading edge and avoiding transitional effects. Even when the geometry is not tripped, the majority of the experimental data in literature are the vicinity of the separation location and/or in the wake. So far, the only experimental database that is exploitable1 and fully describes transition at the surface was provided in Kreplin et al. (1985). The transition path is reconstructed through the measurements of wall shear stress magnitude and direction. The experimental data from Kreplin are here used for the validation of the two transition models with their crossflow criteria. γ -Re θ and γ predictions are compared for the 6:1 prolate spheroid at 5 • and 15 • angles of attack for Re = 6.5 × 10 6 . These two test cases were chosen within the NATO/AVT-313 collaboration group in occasion of the 2021 Workshop devoted to the 6:1 prolate spheroid geometry. γ transition model with the recalibrated T c1 criterion is also tested for 30 • angle of attack, for which transition on the windward side is entirely dominated by crossflow waves.

ISIS-CFD AT A GLANCE

The solver ISIS-CFD, available as a part of the FINE TM /Marine computing suite distributed by Cadence Design Systems, Inc., is an incompressible multiphase unsteady Reynolds-averaged Navier-Stokes (URANS) solver, mainly devoted to marine hydrodynamics. It is based on a fully-unstructured (face-based) finite volume discretization with specific functionalities needed for multiphase flows and industrial applications, see [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF], [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF], Wackers et al. (2012) and Wackers et al. (2017). The method features several sophisticated turbulence models: apart from the classical two-equation k-ε and k-ω models, the anisotropic two-equation Explicit Algebraic Reynolds Stress Model (EARSM), as well as Reynolds Stress Transport Models (SSG/LRR). All models are available with wall-function or low-Reynolds near wall formulations. Hybrid RANS/LES turbulence models based on Detached Eddy Simulation (DES-SST, DDES-SST, IDDES) are also implemented and have been validated on automotive flows characterized by large separations, see Guilmineau et al. (2013). The two γ -Re θ , γ transition models equipped with crossflow criteria have been recently included and validated in [START_REF] Rubino | Laminar-to-Turbulence Transition Modeling of Incompressible Flows in a RANS Framework for 2D and 3D Configurations[END_REF]. Laminar-to-turbulence transition simulations presented in the discussion are for single fluid steady flows. Convective fluxes of transition, turbulence and momentum equations are discretized using AVLSMART scheme, a NVD diagram discretization scheme, based on the third-order QUICK scheme, [START_REF] Leonard | A stable and accurate convective modelling procedure based on quadratic upstream interpolation[END_REF].

LOCAL CORRELATION TRANSITION MODELS γ AND γ -RE θ MATHEMATICAL FORMULATIONS

This section is devoted to the mathematical formulation of the local correlation transition models (LCTM) γ -Re θ and γ. The aim is to recall the main terms, important for the upcoming discussion.

γ -Re θ Formulation

The intermittency, γ, transport equation is given by:

∂ (ργ) ∂t + ∂ (ρu j γ) ∂ x j = P γ -E γ + ∂ ∂ x j µ + µ t σ γ ∂ γ ∂ x j , (7)
where the intermittency production term P γ is defined as:

P γ = c a1 F length ρS(γF onset ) 0.5 (1 -c e1 γ). (8)
S is the mean strain magnitude, c a1 and c e1 are model constants. P γ ∼ 0 until the local vorticity Reynolds number Re V exceeds the critical momentum thickness Reynolds number Re θ c , as determined through the function F onset . The latter is built as:

Re V = ρy 2 S µ , (9)
F onset,1 = Re V 2.193Re θ c , (10)
F onset,2 = min(max(F onset,1 , F 4 onset,1 ), 2), ( 11)

F onset,3 = max 1 - R T 2.5 3 , 0 , (12)
F onset = max(F onset,2 -F onset,3 , 0), (13)
where

R T = ρk ω µ (14)
is the eddy viscosity ratio. 11), assures a rapid change from zero to one, as the onset transition criterion F onset,1 > 1 is met. F onset,3 , Eq.( 12), keeps active the production term of γ throughout the transition process, as the eddy viscosity ratio R T increases. F length is an empirical correlation function of Re θ t , and it controls the transition length. The destruction/relaminarization term E γ is defined as:

E γ = c a2 ρΩF turb γ(c e2 γ -1), ( 15
)
where Ω is the vorticity rate magnitude, c a2 and c e2 are constants. The function F turb is meant to deactivate the destruction/relaminarization term outside of the laminar boundary layer and it depends on the eddy viscosity ratio R T . The boundary conditions for γ are zero normal flux at the wall and 1 at the inlet to avoid any contamination of the free-decay of turbulence variables in the free-stream. γ is enforced to be equal to 0.02 at the wall. This lower bound does not impact the production of turbulence kinetic energy in the viscous sub-layer.

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 .
(16) The source term P θ t is defined as:

P θ t = c θ t ρ t (Re θ t -Re θ t )(1 -F θ t ). (17)
Through P θ t , Re θ t is forced to be equal to the local value outside the boundary layer, Re θ t , obtained from an empirical correlation. In the boundary layer the function

F θ t is 1, such
λ θ = ρθ 2 µ dU dS , (18)
where dU dS is the velocity derivative along the streamline direction:

dU dS = u i u j U 2 ∂ u i ∂ x j , (19)
and the turbulence intensity Tu:

Tu = 100 2k/3 U . ( 20
)
The boundary conditions for Re θ t are zero flux normal at wall. Its value at the inlet depends on the empirical correlation for Re θ t , supposing zero-pressure gradient.

The intermittency that enters in the coupling with the chosen turbulence model is the effective intermittency defined as :

γ e f f = max(γ, γ sep ), (21)
where

γ sep = min s 1 max 0, Re V 3.235Re θc -1 F reattach , 2 F θ t , (22)
with s 1 constant. This modification is needed to account for separation-induced transition. Through Eq.( 21), the intermittency is allowed to be bigger than one when laminar separation occurs, boosting the production of turbulence kinetic energy downward the separation point.

The turbulence model used for the coupling is the kω SST (2003), [START_REF] Menter | Ten years of industrial experience with the SST turbulence model[END_REF]. Hereafter, only the modification due to the transition model coupling are reported. The detailed turbulence model formulation can be found in the original reference. γ e f f is multiplied by the production and destruction terms of turbulence kinetic energy transport equation, as follows:

∂ (ρk) ∂t + ∂ (ρu j k) ∂ x j = Pk -Dk + ∂ ∂ x j µ + µ t σ k ∂ k ∂ x j , (23) Pk = γ e f f • P k , (24) Dk = min(max(γ e f f , 0.1), 1)D k . ( 25
)
The production term P k is computed using the Kato-Launder modification (P k = µ t SΩ), [START_REF] Kato | The modelling of turbulent flow around stationary and Vibrating square cylinders[END_REF]. D k = β * ωk is the original destruction term of the turbulence model. The equation for the turbulence frequency ω rests unchanged:

∂ (ρω) ∂t + ∂ ∂ x j ρU j ω -(µ + σ ω µ t ) ∂ ω ∂ x j = γρΩ 2 -β ρω 2 + 2(1 -F 1 ) ρσ ω2 ω ∂ ω ∂ x j ∂ k ∂ x j . (26)
Nevertheless, the blending function F 1 in Eq.( 26), responsible of the switching from kω to kε, is reformulated. It is forced to be equal to one in the laminar region and it is redefined as:

F 1 = max(F 1orig , F 3 ), (27)
F 3 = e ( Ry 120 ) 8 , R y = ρy √ k µ , (28)
where F 1orig is the definition of F 1 in the k -ω SST ( 2003) model.

γ model

In γ transition model formulation, the transport equations are reduced from two to one. All the empirical correlations that enter the model are now computed using only local quantities. The transport equation for the intermittency reads as:

∂ (ργ) ∂t + ∂ (ρu j γ) ∂ x j = P γ -E γ + ∂ ∂ x j µ + µ t σ γ ∂ γ ∂ x j .
(29) The intermittency source term P γ is constructed as:

P γ = F onset [F length (ρS(1 -γ)γ)], (30)
where F length is set to the constant value 100. The function F onset is built as:

F onset,1 = Re V 2.2Re θ c , (31)
F onset,2 = min(F onset,1 , 2.0), (32)
F onset,3 = max 1 - R T 3.5 3 , 0 , (33)
F onset = max(F onset,2 -F onset,3 , 0). ( 34
)
In this formulation, Re θ c in Eq.( 31) is still given by an empirical correlation, but it is function of local quantities: the local turbulence intensity Tu L and pressure gradient parameter, λ θ ,L . The turbulence intensity within the boundary layer is expressed as:

Tu L = min 100 2k/3 ωy , 100 , (35)
while the local pressure gradient parameter is expressed by: λ θ ,L = -7.57 • 10 -3 dv dy

y 2 ν + 0.0128, ( 36
)
where dv dy is the wall-normal derivative of the wall-normal velocity component v of the velocity vector -→ u . It is computed as:

dv dy ≡ ∇( - → n • - → u ) • - → n . (37)
The constant 0.0128 in Eq.( 36) accounts for the fact that dv dy is not zero in the middle of the boundary layer, where transition occurs, for zero-pressure gradient flows. The use of dv dy comes from two-dimensional flow considerations. The relaminarization/destruction term E γ depends on the magnitude of the absolute vorticity rate and it is defined as:

E γ = c a2 ρΩγF turb (c e2 γ -1). ( 38
)
Finally, in the coupling with kω SST (2003) turbulence model, as for γ -Re θ , only the production and destruction terms of k, Eq.( 23), are modified. The turbulence kinetic energy production term is given by the sum of two different terms. A primary production term P k is defined as:

P k = γP k , ( 39
)
where P k is also computed using Kato-Launder modification. Then, an additional production term P lim k is included to account for separation-induced transition and it is defined as:

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

F lim on = min max Re V 2.2 • 1100 -1, 0 , 3 , (41)
and C SEP and C k are constants. The modified destruction term is given by:

D k = max(γ, 0.1) • D k , (42)
where D k comes from the original kω SST (2003) turbulence model. The blending function between kω and kε formulations, F 1 , is reformulated as in γ -Re θ .

The boundary conditions for γ are the same as for the transport equation for the intermittency in γ -Re θ .

The transition models equations, as presented above, account for bypass and streamwise transition, either due to T-S waves or separation-induced. Nevertheless, the present study aims to analyze the flow around the 6:1 prolate spheroid at incidence, configuration for which crossflow instabilities play a dominant role. The two crossflow criteria chosen in this study are discussed in the following section. In order to describe the philosophy behind T c1 crossflow criterion and its new recalibration, the Falkner-Skan-Cooke (FSC) equations are introduced. This set of equations is the key link between local and non-local quantities (Re δ 2 , H 12 , ...).

FALKNER-SKAN-COOKE (FSC) EQUATIONS

The Falkner-Skan-Cooke (FSC) are the extension of the two-dimensional laminar boundary layer Falkner-Skan equations by considering the spanwise velocity component for a flow over an infinitely yawed wedge at zero angle attack, [START_REF] Cooke | The boundary layer of a class of infinite yawed cylinders[END_REF]. Fig.

(1) presents an infinite swept wing geometry. Two different coordinate systems are identified: the wing-attached (x, y, z), with x aligned with the chordwise direction, and the reference system (x 1 , y, z 1 ), where x 1 is aligned with the external inviscid streamline. The three-dimensional boundary layer equations system for the infinite swept wing, i.e. the derivatives along the span is zero ∂ ∂ z = 0 , are:

∂U ∂ x + ∂V ∂ y = 0, (43)
U ∂U ∂ x +V ∂U ∂ y = U e dU e dx + ν ∂ 2 U ∂ y 2 , ( 44
)
U ∂W ∂ x +V ∂W ∂ y = ν ∂ 2 W ∂ y 2 . ( 45
)
with boundary conditions:

U = V = W = 0 at y = 0, (46)
U →U e and W → W e as y → ∞.

(47)

Figure 1: Infinite swept wing reference systems used for the formulation of the FSC equations.

It is supposed that the inviscid chordwise velocity at the boundary layer edge U e follows a potential law in the coordinate x, normal to the leading edge, and the spanwise velocity (parallel to the leading edge) W e is constant. These two components can be written as:

U e ∼ U ∞ x L m , (48)
W e = const, (49)
where L is the characteristic length and U ∞ is the free-stream longitudinal velocity. m is the streamwise pressure gradient and it is expressed as:

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

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

Ψ = 2U ∞ νL m + 1 x L (m+1)/2 f (η), (52)
such that U = ∂ Ψ ∂ y , V = -∂ Ψ ∂ x
, the continuity equation is automatically satisfied and Eq.( 44) becomes:

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

β = 2m m + 1 . ( 54
)
The dash in Eq.( 53) denotes the differentiation with respect to η. On the other hand, given W = W e g(η), Eq.( 45) becomes:

g + f g = 0. (55)
Eq.( 53)-( 55) are the Falkner-Skan-Cooke equations, with boundary conditions:

f , f , g → 0 for η → 0, (56) f , g → 1 for η → ∞. ( 57
)
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 1 e . Their expression is given by:

U 1 /U 1 e = f cos(φ ) 2 + g sin(φ ) 2 , ( 58
) W 1 /U 1 e = (g -f ) cos(φ ) sin(φ ), ( 59
)
where φ 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. This angle is constant along the wall normal height of the boundary layer and it is defined such that:

tan(φ ) = W e U e . (60)
According to the definition in Eq.( 59), W 1 = 0 for φ = 0 • and φ = 90 • , but also for zero-pressure gradient flows, because gf = 0. The last condition does not occur in real physical flows, because crossflow velocity also exists for zero-pressure gradients three-dimensional flows.

TC1 CRITERION

Arnal's C1 criterion in Eq.( 2) can be rewritten in the form:

Re δ 2 f (H 12 )150 = 1, (61)
where f (H 12 )150 is the value of the crossflow Reynolds number Re δ 2 t at the transition onset. [START_REF] Menter | Development of a RANS-based Model for Predicting Crossflow Transition[END_REF] propose to split the C1 criterion, as expressed in Eq.( 61), in three different terms, each one accounting separately for the parameters that affect the crossflow transition process. The criterion can, then, be expressed as:

Re δ 2 f (H 12 ) ∼ F(H 12 )XRe stream . ( 62
)
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 accounts for the Reynolds number effect. Based on Eq.( 62), a local approximation of the C1 criterion, here referred to as T c1, is given by:

T c1 = 1 150 G • Ψ • Re V max , (63)
where

Ψ ∼ X = Re δ 2 Re V max , (64)
G ∼ F(H 12 ) = 1 f (H 12 ) . ( 65
)
The T c1 criterion, as formulated in Eq.( 63), is mathematically equivalent to Arnal's C1 criterion. The function Ψ measures the ratio of the crossflow to the streamwise strength. As proposed by Menter & Smirnov, it is constructed using the wall-normal change of the normalized vorticity. Ψ is defined as:

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

ψ i = ∂ ω i ∂ x j n j , where ω i = ω i | - → ω | . ( 67
)
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. In addition, this quantity describes the three-dimensionality of the boundary layer, being Ψ = 0 for 2D flows, given that ∂ ∂ z = w = 0. The T c1 local approximation in Eq.( 63) is evaluated using the solutions of FSC equations, which are the key to define the link between non-local and local quantities. The aim is to understand how well the local functions that enter the T c1 approximate the boundary layer quantities appearing in the C1 criterion.

The FSC equations are solved in the parameters range:

0 <β ≤ 1, ( 68
) 0 • <φ < 90 • . ( 69
)
β is restricted to positive values, because crossflow instabilities occur for accelerated flow, in a favorable pressure gradient region. β = 1 is the case of 90 • wedge, i.e. the 2D stagnation flow, and it is the highest possible acceleration parameter. φ is the sweep angle, as defined in Eq.( 60). The functions appearing in Eq.( 63) are evaluated at the particular position η = η max in the wall-normal direction, such that T c1(η max ) reaches its maximum value. Within the FSC framework, the quantity Ψ is approximated only considering the derivatives with respect to the normal direction η. Thus, the vorticity components that exist are ω x and ω z , computed using the velocities expressed in Eq.( 58) and Eq.( 59). At first, the comparative ratio R(β , φ ) = 0.684Ψ/(XF(H 12 )) is considered and it is shown in Fig.

(2). The constant 0.684 is chosen in order to achieve the value R ∼ 1 for β → 0. It is observed that at the upper corners of the domain the ratio R significantly departs from the targeted value of 1. This behavior suggests that the crossflow indicator Ψ does not fairly reproduce the quantity XF(H 12 ). A function G, that accounts for the pressure gradient in the streamwise direction, i.e. the non-linear term F(H 12 ), is needed to correct the ratio R and reduce the deviation from the targeted value R = 1. The main obstacle to construct G is that the two independent FSC parameters, β and φ , are not known in the local formulation. In their publication, Menter & Smirnov propose to construct G as a function of the dimensionless quantity dv dy y2 ν , that approximates the pressure gradient parameter λ θ . G is expressed 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 MS (λ θ ,CF ) = 0.684 g(λ θ ,CF ) , (70)
where the cubic polynomial g(λ θ ,CF ) is fitted through the clouds of points. The pressure gradient parameter λ θ ,CF is expressed as:

λ θ = - dv dy y 2 ν , (71)
λ θ ,CF = 0.1111 • λ θ + 2.3, ( 72
)
λ θ ,CF = min[max(λ θ ,CF , 0), 0.7], (73)
where the wall-normal derivative of v in Eq.( 71) is defined as in Eq.( 37). 3b), projected onto the plane (R, β ), where β is expressed through λ θ ,CF . The spread around the target value of 1 is reduced to around 25% for sweep angles φ lower than 60 • . This is the upper limit for φ chosen by Menter & Smirnov for the resolution of FSC and the calibration of their criterion. Nevertheless, the error committed on higher angles of attack is still important. The choice of limiting the sweep angle φ to less than 60 • is acceptable if one is interested in the prediction of crossflow around wing-like geometries. Nevertheless, this is a considerable low limit for complex 3D configurations, such as the 6:1 prolate spheroid, which is a geometry significantly more swept compared to any swept wing 2 . In order to widen the range of applications of the T c1 criterion to more complex 3D geometries, the strategy proposed by the authors is to include in the local formulation of the T c1 criterion a local approximation of the sweep angle φ . This is an important parameter to account for, because the pressure gradient parameter λ θ strongly depends on the sweep angle φ , as shown in ν and φ fitted using a least squares method. Its formulation is referred to as λ CF and it is given by: (74) For numerical reasons, the value of λ CF is further bounded as:

λ CF = min max λ CF , 0 , 0.16 . (75)
Then, the function G, in Eq.( 63), that accounts for the streamwise pressure gradient, is constructed as the surface G new (λ CF , φ ). It is expressed as:

G new (λ 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.04903 10 -5 φ λ 2 CF + 0.03326 λ 3 CF , (76)
where both λ CF and φ are normalized with respect to their mean value and standard deviation. For the construction of the function G new , λ CF is evaluated at the position η max where the T c1 reaches its maximum.

The new ratio R = G new Ψ/(XF(H 12 )) is shown in Fig.( 5) with respect to the dependent variables λ CF and φ . It is observed that the benefit of introducing the sweep angle φ is a considerable reduction of the spread of the ratio R around the value of one, with a maximum deviation of less than 10%. Now, the definition of φ as the angle between the external potential flow direction, U 1 e , aligned with the reference coordinate system (x 1 , y, z 1 ), and the wing-attached reference system (x, y, z), is not CFD-compatible. For this reason, following Högberg and Hennigson (1998), the 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), since, by FSC assumption, the pressure gradient is zero in the spanwise direction. A local sweep angle φ L definition can be computed as in [START_REF] Choi | Recent improvement of a correlation-based transition model for simulating three-dimensional boundary layers[END_REF], as follows:

φ L = arccos - → u wt • --→ (∇p) wt || - → u wt |||| --→ (∇p) wt || , φ L = min[φ L , π -φ L ].
(77) -→ u wt and --→ (∇p) wt are the tangential projection at the wall of the local velocity vector and the pressure gradient. 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, where φ L recovers the original definition of φ from Eq.( 60).

This approximation of the local sweep angle uses the local velocity vector, therefore it makes the recalibrated T c1 criterion not Galilean invariant. In order to restore invariance of the new formulation with respect to Galilean transformations, it is proposed to use the relative velocity vector -→ u rel instead of the local velocity vector in Eq.( 77). It is defined as:

- → u rel = - → u -- → u wall . ( 78
)
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 3 is stored in a table. This implementation has its own limitations, because the research of the point at the wall might be troublesome at the junctions between multiple bodies with possibly different velocities, as the rotor blades of an helicopter.

Once established a "Galilean-invariant" like formulation, the T c1 criterion is finally included in γ formulation, modifying the F onset function in the transport equation for the intermittency γ, Eq.( 30). A new F onset,CF that triggers the production of γ and based on the T c1 criterion is summed up to the F onset function of the original formulation. F onset,CF is defined as:

T c1 = G new Ψ Re V c 150 , ( 79
)
F onset1,CF = T c1 (80) F onset2,CF = min[max(F onset1,CF , 0), 2],
(81)

F onset3,CF = max(1 -(R T /a) 3 , 0), (82) F onset,CF = max(F onset2,CF -F onset3,CF , 0), (83) a = 1.5, c = 0.6. ( 84
)
Through the proportionality constant c in Eq.( 79), it is accounted for the difference between the critical crossflow Reynolds number, at which the intermittency starts to increase, and the crossflow Reynolds number at transition location. The new F onset function in the transport equation for the intermittency γ, Eq.( 30), becomes:

F onset = F length F onset + F length,CF F onset,CF , (85)
where F length,CF = 5. The latter parameter has been set considering that crossflow instabilities develop on a longer length compared to T-S waves for which F length = 100 . The constants, a and F length,CF , c were obtained from numerical calibration on the 6:1 prolate spheroid at 15 • .

HELICITY-BASED CRITERION

For the extension to crossflow transition in γ -Re θ model, it has been considered and implemented in ISIS-CFD the helicity-based criterion by DLR. The onset of crossflow transition depends on the critical crosswise Reynolds number ratio: Re He,max Re +

He,t = 1, (86)
where Re + He,t at transition is given by a correlation numerically determined. It is based on the numerical data obtained for different configurations the ONERA D 3 The point that minimizes the distance from the cell center to the wall. profile, NLF(2)-0415, and the 6:1 prolate spheroid. Re + He,t is expressed as a function of the shape factor H 12 , as follows:

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

(87)

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

λ + = ρl 2 µ d| - → u e | d - → s , (88)
l = 1 C He,max 2 15 y. ( 89
)
The length scale l represents the momentum thickness at the point where the helicity Reynolds number reaches its maximum within the boundary layer. C He,max is a constant. H 12 in Eq.( 87) is then substituted by the correlated H + 12 based on Cliquet's correlation derived for zero sweep angle flows, [START_REF] Cliquet | Application of laminar-turbulent transition criteria in Navier-Stokes computations[END_REF]. The inclusion of the helicity-based crossflow criterion within γ -Re θ formulation is achieved through the definition of an additional F onset,CF function that is summed up to F onset of Eq.( 10). F onset,CF is given by:

F onset1,CF = Re He CRe + He,t (90)
F onset2,CF = min[max(F onset1,CF , F 4 onset1,CF ), 2] (91) F onset3,CF = max 1 - R T 2 3 , 0 (92)
F onset,CF = max(F onset2,CF -F onset3,CF , 0) (93)
where the constants C = 0.7 and F length,CF = 5 result from numerical calibration. The explicit use of the velocity within 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.

6:1 PROLATE SPHEROID SIMULATIONS

This section is devoted to the discussion of the computations performed with ISIS-CFD around the 6:1 prolate spheroid at incidence. As mentioned in the introduction, it is referred to the experiments performed by Kreplin in 1985 at DLR for the models validation. These experiments were performed in the 3m × 3m low speed wind tunnel at DLR Gottingen, around a 6:1 prolate spheroid of length of 2.4m. Measurements at the surface were obtained using surface hot film probes: the 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 longitudinal axis of the specimen. 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 uncertainty is related to the fact that the hot-film probes were calibrated on flat tunnel walls for a 2D turbulent boundary layer. Around the 6:1 prolate spheroid at incidence, the laminar region is quite extended.

Presented computations are run for Re = 6.5 × 10 6 and three different angles of attack. Zones of pure crossflow are observed in the middle of the inclined prolate spheroid and they become wider as the inclination increases. For α = 5 • , 15 • , a comparison is proposed between γ with the recalibrated T c1, referred to as γ + CF model, and γ -Re θ with the helicity-based criterion, referred to as γ -Re θ +CF. Only for α = 15 • , simulations by γ using the T c1 criterion, as published by Menter & Smirnov, here referred to as γ + T c1 MS model, are presented. A summary of the notations used for the different models with their respective crossflow criterion is given in Table(1). Simulations for the angle of attack α = 30 • were run using only γ +CF model.

Unfortunately, no specific indication on the free-stream conditions was given in the experimental report. It is mentioned that Tu varies from 0.1% to 0.3% close to the nose of the body. Not knowing the free-stream Tu value is a big limitation, that does not allow to perform a rigorous validation exercise. Indeed, experimental conditions cannot be repeated with exactitude. For the presented exercise, the free-stream conditions were chosen in order to have a value of Tu ∼ 0.15% in the vicinity of the leading edge. This value is obtained from:

Tu = Tu in 1 + 3(x * -x * in )β Tu 2 in Re 2(ν t in /ν) -β * β 0.5 . (94)
Table 1: Notations for the transition models and crossflow criteria presented in the discussion.

Computational Domain

The grids used for the simulations 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. The flow is aligned with x-axis.

The grids are for half of the geometry, making use of a symmetry plane. All the mesh details are given in Table (3.

h i = ∑ N cell i ∆V i N cell
is the typical cell size, with ∆V i the volume of the i-th cell and N cell is the total number of cells. 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 counts 126016 cells on the surface of the spheroid, and a total of 42.5M volume cells. The remaining four grids are obtained from the finest one using the coarsening factors of 0.875, 0.75, 0.625 and 0.5. The coarsest grid from different points of view is shown in Fig.( 6)-( 8). All the computations around the 6:1 prolate spheroid were run on national HPC resources. Only for α = 15 • it has been performed a mesh convergence study. For each grid, the interpolated solution from its correspondent "one level coarser" grid is used as initial solution. Convergence is controlled by a gain of minimum four orders of the normalized residuals for all the turbulence and transition variables and also by forces convergence. Computations by γ + CF on the coarsest grid, Grid5, were run on 280 processors, ∼ 9000 CPU hours of simulation were needed up to convergence. On the same grid and for the same number of processors, simulations using γ -Re θ + CF required ∼ 13000 CPU hours to converge. Such long simulations time are due to the fact that the convergence of the non-linear residuals in L 2 norm is very noisy, because of the presence of several min, max limiters in the original γ and γ -Re θ formulations. The faster convergence of γ respect to γ -Re θ is due to the discard of the transport equation for Re θ in γ formulation. Nevertheless, it is worthwhile to mention that γ model is more sensitive to the mesh density than γ -Re θ , see [START_REF] Rubino | Laminar-to-Turbulence Transition Modeling of Incompressible Flows in a RANS Framework for 2D and 3D Configurations[END_REF].

Simulations on finer meshes were run on higher numbers of processors, up to a maximum of 784.

Results: α = 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.( 9) and (10). It is observed that an estimated order p = 2 on the friction component is obtained, while the pressure drag presents a quasi-2nd order of convergence with respect to the grid refinement. The behavior of the pressure is most probably related to the noisy convergence of the non-linear residuals on such fine meshes. For what concerns the γ + T c1 -MS, no information were given in the publication by Menter & Smirnov about the criterion inclusion in γ formulation, except from their definition of the T c1, here referred to as T c1 MS . The latter is given by:

T c1 MS = G MS ΨRe V 150 . (95)
In Eq.( 95), G MS formulationis given by Eq.( 70), using the λ θ approximation of Eq.( 71)-( 73). In T c1 MS no constant c was foreseen to account for the difference between the critical Reynolds number and the one at transition onset, as in Eq.( 79). Finally, the F onset,CF in γ + T c1 -MS model formulation used in the presented computations is given by:

F onset1,CF = T c1 MS (96) F onset2,CF = min[max(F onset1,CF , 0), 2],
(97) 99)

F onset3,CF = max(1 -(R T /a) 3 , 0), (98) F onset,CF = max(F onset2,CF -F onset3,CF , 0), (
a = 1.5. ( 100
)
The final F onset that enters the intermittency production term is given by Eq.( 34 the sweep angle, it is accounted for its impact on the overall transition process and that makes a huge difference when dealing with non-wing-like geometries. Focusing on the comparison between γ +CF and γ -Re θ +CF, the two transition models with their respective CF criterion perform very similarly, with few exceptions. γ -Re θ + CF is able to predict the upper side kink, that marks the change of the transition process, which is on the contrary absent in C f contours by γ +CF. This difference is clearly visible from the predicted C f girthwise distributions at the sections X/L = 0. 12d). 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 numerically crossflow transition. This laminar tongue is grid dependent and it becomes shorter on coarser grids. Over all, the predicted C f by γ + CF and γ -Re θ + CF 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 uncertainties. The underestimation of the C f is mainly due to the fact that, at this angle of attack, interaction between T-S and CF waves is expected and none of the two formulations really accounts for the modes interaction. The sum of the F onset and F onset,CF functions, is not enough to account for the two modes interactions.

Approaching the trailing edge of the body, none of the two models is able to predict flow transition for φ < 40 At least at the windward symmetry plane, transition is hardly due to pure crossflow. Actually, at the symmetry plane, no crossflow transition can happen, because the crossflow component of the velocity is zero. The incorrect flow prediction at the windward side is due to the fact that γ +CF or γ -Re θ +CF, as here presented, do not account for any transition mechanisms other than T-S waves and CF instabilities. The possible missing mechanisms are attachment line instabilities, traveling crossflow modes, and the non-linear interaction between T-S waves and crossflow instabilities. Indeed, on the 6:1 prolate spheroid, for the treated angle of attack, some interaction between T-S and CF is expected. Neglecting the interaction of these two modes, not only affects the transition front closer to the windward symmetry plane, but also the skin friction values are considerably underestimated. A possible strategy to account for T-S and CF interaction would be to use the local sweep angle φ L as a measure of the two modes interaction, as proposed in [START_REF] Choi | Recent improvement of a correlation-based transition model for simulating three-dimensional boundary layers[END_REF]. The idea is to weigh up the sum of the F onset and F onset,CF functions with respect to φ L , knowing that the interaction is stronger for φ L = 90 • and it weakens as the sweep angle decreases. This has not been tested yet in ISIS-CFD, but further studies are going on in this direction.

The other transition mechanism that could play an important role are traveling crossflow waves. Indeed, the level of turbulence intensity reported in the experimental report by Kreplin is in a range for which traveling modes might be important. As shown in the classic experiments in Deyhle and Bippes (1996), traveling crossflow waves prevail over stationary modes in the transition process for a turbulence level Tu > 0.2% and on smooth surface, as it is the case of the 6:1 prolate spheroid studied by DLR. Nevertheless, nor the recalibrated T c1 criterion nor the helicity-based can predict transition due to traveling crossflow, because the two criteria are built with respect to experiments run in a very low free-stream turbulence environment, where stationary crossflow waves are expected to be dominant in the transition process. In addition, both criteria do not explicitly depend on the turbulence level intensity, which is the parameter that plays the major role in transition due to traveling crossflow, making them unsuitable, as currently formulated, for this transition mechanism.

Finally, along the windward symmetry plane, leading edge instabilities may cause transition, as well. These instabilities are related to 2D boundary layer developing at the attachment-line and they can cause the flow to become turbulent along the longitudinal direction. The windward symmetry plane is a typical case where transition is dominated by streamwise instabilities, but the streamwise criterion Re θ -based cannot predict it. As explained in [START_REF] Arnal | Three-dimensional boundary layer: laminar-turbulent transition[END_REF], for complex 3D flows, where the streamlines are far from parallel, there is no implicit relation between Re θ and the physical distance along which the instabilities propagate. For this reason, criteria which involve boundary layer parameters fail to predict transition at the attachment lines, explaining both models failure along this symmetry plane.

Skin friction lines predicted by γ +CF are shown in Fig. (13). A very similar topology is obtained using γ -Re θ +CF. The model predicts an envelope of converging wall streamlines on the top of the spheroid coming from the windward and leeward sides, which is representative of an open-separation, according to the definition in Surana et al. (2006).

The following two sections are devoted to the analysis of the flow around the 6:1 prolate spheroid for α = 5 • , 30 • . Attention is focused on the results obtained by γ +CF, using the recalibrated T c1 crossflow criterion introduced above. The presented computations were computed on Grid3 of Table (3 18). The high value of skin friction close to the transition front, as observed experimentally, might be related to the experimental value of Tu in the freestream.

In the experiments, it is mentioned that Tu could be up to 0.3%, while the value imposed in the simulations is lower, i.e. Tu = 0.15%. 18a) and (18b) roughly indicates the transition location within the bubble. This is not visible in γ + CF results. This discrepancy is related to γ model original formulation. The criterion that accounts for separation-induced transition is the additional production term in Eq.( 40). Re V in the first 20% of the spheroid length does not exceed the limiting value of 2200 to activate the additional production term of turbulence kinetic energy, Eq.( 41). This behavior might be related to the diverging streamlines predicted at the leeward symmetry plane, close to the leading edge, 19). 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 minimum of the magnitude of the wall shear stress, according to [START_REF] Simpson | Aspects of turbulent boundary-layer separation[END_REF]. This is a fair approximation because the skin friction lines are converging from different sides of the prolate spheroid.

The second separation location occurs at the leeward side approximately around X/L ∼ 0.55, but predicted skin friction lines are converging further downstream. 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. (18d), due to the velocity fluctuations induced by the separated vortex flow, is considerably pushed downward at the rear side of the spheroid.

Due to the large separation, which the flow undergoes at such an high angle of attack, it cannot be excluded that the RANS turbulence model does not affect the overall results at separation. Indeed, RANS turbulence models typically fail in massively separated region. A possible strategy would be to use Hybrid RANS/LES, that uses RANS transition models within the boundary layer and LES in the separated region. To this matter, efforts have been going on to combine transitional RANS model into an hybrid RANS-LES method, for instance in [START_REF] Hodara | Hybrid RANS-LES closure for separated flows in the transitional regime[END_REF]. Nevertheless, as for turbulence models, the key issues for hybrid RANS/LES transition models is the location of the RANS to LES switch. Indeed, 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, such as the relaminarization of the boundary layer. 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.

CONCLUSION AND FUTURE WORKS

This work was devoted to the analysis of the local correlation transition models γ and γ -Re θ performances on the 6:1 prolate spheroid at incidence. For the chosen test cases, crossflow instabilities play a fundamental role in the flow transition. In order to account for it, an original model variant of γ model is presented, the so-called γ + CF model, that includes the T c1 crossflow criterion. This criterion, originally published by Menter & Smirnov, is based on the C1 criterion by Daniel Arnal. In the present discussion, a recalibration of the T c1 is proposed, that uses a CFD-compatible formulation of the local sweep angle. The inclusion of the latter extends the criterion to crossflow transition predictions on complex three-dimensional geometries. γ + CF variant is proven to perform very well on the 6:1 prolate spheroid at the different angles of attack. Indeed, within the newly re-calibrated T c1, it is accounted for the variation of the pressure gradient parameter in the streamwise direction with respect to the local sweep angle and its influence on the overall crossflow transition process. This modification strongly enhances the original T c1 (T c1 MS ) by Menter & Smirnov, as shown for α = 15 • .

For the angles of attack α = 5 • , 15 • , γ + CF is also compared to γ -Re θ model coupled to the DLR helicity-based criterion for crossflow transition, the so-called γ -Re θ + CF model.

Over all, the two transition models crossflow variants, γ + CF and γ -Re θ + CF, perform very similarly. Their predictions are quite satisfactory for α = 15 • , at which transition is dominated by crossflow modes. Nevertheless, deviation from experiments becomes important at lower incidence, as at α = 5 • . In fact, none of the transition models variant predicts the flow transition at the windward side, close to the symmetry plane, as expected in the experimental results. This unpredicted flow feature is due to missing physics within the transition models formulations. Transition close to the windward symmetry plane is hardly due to solely crossflow waves. Physical mechanisms that are not accounted for and might play an important role are non-linear interactions between crossflow and T-S waves, attachment line instabilities, and traveling crossflow instabilities. In spite of their similar performance for the treated test cases, the philosophy behind the two crossflow criteria used for these simulations is considerably different. The T c1 coupled to γ is constructed using Falkner-Skan-Cooke solutions, while the helicity-based He crossflow criterion, coupled to γ -Re θ , is constructed numerically.

A FSC-based calibration, despite the stringent hypothesis of this set of equations, i.e. zero-spanwise gradients, allows the inclusion of additional physical parameters within the criteria in lack of detailed experimental results and this is the reason why the authors have devoted their efforts in the recalibration of this criterion. Indeed, it is important to highlight that the presented transition models are still in a validation phase. The strategy to improve their performance is to include more transition physics within the models correlations, needed to trigger the different transition processes. For example, the inclusion of the local sweep angles not only enhances the T c1 performance, but it also paves the way to the further inclusion of a correlation that accounts for T-S/CF modes interaction. As proposed in Choi et al., φ can also be an indicator of the intensity of the modes interaction and it might be used to weigh up the sum of the different F onset functions. This is one of the possible developments in the coming future.

In general, there is big room for the improvement of the current transition models. However, the main roadblock is the lack of new detailed experimental data. Experiments are needed not only to define new criteria for the different transition mechanisms, but also to improve the correlations that the current models already have, i.e. a streamwise criterion that does not depend on θ and crossflow criteria that also account for surface roughness and turbulence intensity. The two latter parameters strongly affect the transition process. For instance, 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 Tu levels and roughness, as a first level approximation.

It is the authors belief that new experimental results around the 6:1 prolate spheroid are fundamental and could provide data that Kreplin's experiment did not. Repeated transition measurements on the 6:1 prolate spheroid at different levels of free-stream turbulence intensity would help characterizing the sensitivity of the transition behavior to Tu. This might be accomplished by selective removal/addition of turbulence screens ahead of the body. In addition, if a model front half with different levels of surface finish is used, surface roughness effects can be evaluated as well. The visualization of the flow pattern would also be very helpful in determining the transition mechanism. Sublimating chemicals, such as naphthalene based flow visualization, can be typically employed to visualize it.

Additional surface-based instrumentation, as microphones, can be also used to provide unsteady measurements that would give additional information about the frequency contents. 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. Skin friction and pressure distributions are the needed data to validate transition models. However, 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. This is the minimum physics that needs to be included in a useful data set to be a benchmark of quality for further numerical simulations.

Discusser: Rickard Bensow, Dept. Mechanics and Maritime Sciences, Chalmers University of Technology, 412 96 Göteborg, Sweden, rickard.bensow@chalmers.se

The paper gives detailed and rather thorough description of the transition models used and developed; although then by necessity quite technical and brief it is still appreciated to allow for reimplementation. However, the origin of some empirical correlations/fits are unclear or not referenced (e.g. Eq 74); it would be appreciated if the authors could review and complement the paper for a final version. Similarly, the expression ( 94) is just briefly explained. I assume it is a way to compute appropriate inlet turbulence levels to achieve a certain level at the body. Is there a reference for this expression? Question 1. If I understand it correctly, the finest grid uses 176 cells in the circumferential direction for half the body, i.e. approximately 1 • in resolution. Do you conclude, based on your experience, that this is sufficient to capture transition and separation, or would a higher surface resolution be desirable? Question 2.

You only briefly touch upon the difficulties in converging the simulations with the transition models. Did you note any difference in convergence between the models that is worthwhile to take into consideration in choosing between them? Question 3. In many aspects, the recalibrated model behaves very well, while in others it is clear the flow physics is too complex. You did a good job in describing the behavior and give explanations for the shortcomings. You mention a few items for future work, primarily related to trying to include the interaction between the transition mechanisms. Do you thus feel that the current model is good enough for each mechanisms on its own, or are there still improvements that are needed? Question 4. Finally, what kind of experimental data would be needed to help the model development to move forward?

AUTHOR'S REPLY

As observed by the reviewer, some of the correlations mentioned in the paper were not referenced. This choice was dictated by the context and the limited number of pages at the authors' disposal, nevertheless, a detailed description of all the empirical correlations, used in the original models' formulation, as well as for the recalibrated crossflow criterion, can be found in [START_REF] Rubino | Laminar-to-Turbulence Transition Modeling of Incompressible Flows in a RANS Framework for 2D and 3D Configurations[END_REF]. Regarding both λ CF (Eq.( 74)), and G new (Eq.( 76)), these two correlations were obtained using a linear least-squares based fitting. The results of the fitting process are the estimate of the model coefficients, i.e. the coefficients of the polynomial function that fits the data. By steps, the pressure gradient parameter λ θ and the quantity dv dy y 2 ν , used as an approximation of β , are computed in the FSC solver, evaluated at the middle of the boundary layer and stored in a table. Then, the dv dy y 2 ν , solution of the FSC equations, and the sweep angle φ , which varies in the range [0 • , 90 • ], are the variables used in the construction of the correlation λ CF , that best fits the surface λ θ (Fig.( 4)). The same approach is used for the definition of G new , the function that accounts for the streamwise pressure gradient effect. G new is constructed as a function of the variables φ and λ CF , where the latter is evaluated at the wall-normal position, η max , where the T c1 criterion, as reconstructed in the FSC frame, reaches its maximum value. The value of Tu close to the body, as expressed in Eq.( 94), is obtained by the resolution of the transport equation of k and ω in a uniform flow, i.e. only the destruction terms of k and omega exist, and manipulation of the solution for k. k and ω transport equations of the kω SST ( 2003) model for an incompressible and steady flow reduce to:

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

dk * dx * = -β * k * ω * , (103)
dω * dx * = -β (ω * ) 2 . ( 104
)
The superscript * indicates the dimensionless variables,

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

k * = k * in (1 + β (x * -x * in )ω * in ) -β * β , (105)
ω * = ω * in (1 + β (x * -x * in )ω * in ) -1 , (106)
where the subscript in indicates the quantity value at the inlet of the computational domain, positioned at x in .

Rewriting the specific turbulence dissipation rate at the inlet as

ω * in = k * in ν ν t in Re, (107)
substituting Eq.( 107) in Eq.( 105), given that Tu = ( 2k * 3 , Eq.( 94) is obtained by manipulation.

1. It is commonly accepted that transition models are extremely sensitive to mesh density. Unfortunately, no general mesh requirements can be outlined. For the present exercise, a set of five topologically similar grids was used for the simulations. The results presented for α = 15 • are computed on the finest structured mesh Grid1, which has a maximum of 176 cells in the circumferential direction, i.e. at the middle of the spheroid. Nevertheless, results were also computed on coarser meshes, those presented in Table (3), the coarsest of which had a maximum of 140 cells in the circumferential direction. During the grid refinement study, the authors noticed that transition results were strongly affected by the mesh density in the streamwise direction rather than in the circumferential one. For instance, it was noted that at coarser grid resolution, predictions were spuriously exhibiting transition at or near the windward symmetry plane depending on the mesh resolution along the streamwise direction and keeping unchanged the mesh resolution in the circumferential direction.

As mesh refinement was increased and the solution approached grid independence, laminar flow emerged on the windward side of the model. The explanation for premature transition at coarser mesh refinement levels is accepted to be the fortuitous cancellation of numerical and modeling errors due to mesh density. While a clear dependence of the transition prediction on the circumferential mesh density was not noted on this specific set of meshes, the resolution in the circumferential direction might affect separation. Unfortunately, separation was not considered in this exercise, and no further indications can be given.

2.

Compared to turbulence models, both transition models demand an elevated number of non-linear iterations. This behavior is expected because of the most complicated physics that has to be resolved. In addition, the coupling between the transition and turbulence equations introduces additional difficulties in converging the flow solution. In particular, on grids as dense as those used for this exercise, the solution requires a long time to converge. These very dense grids in the wall-normal direction are needed to accurately capture transition in the boundary layer, i.e. the large gradients of the intermittency γ. 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. 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. The high number of non-linear iterations is also associated with a noisy iterative convergence: the normalized residuals decrease rapidly at the beginning of the simulation, then start oscillating and stall. The iterative convergence is disturbed by the presence of the several min and max limiters within the transition models' formulation. These functions are continuous, but not smooth, and do not promote stability. Although the comparison between the two transition models' behavior, as well as iterative convergence, strongly depends on the selected flow solver, the one-equation γ tends to be easier to converge than the two-equation γ -Re θ model. However, the numerical robustness of the γ model degrades significantly on coarse grids. This behavior is related to the fact that in γ formulation all the correlations are locally computed in the middle of the boundary layer. On a mesh coarse both in the wall-normal and/or streamwise directions, it is observed what can be defined as an "hysteresis effect": once turbulence is installed, numerically the model is not able to destroy it in favor of laminar/transitional flow, because of the discretization errors committed on under-resolved meshes. Such kind of behavior is not observed when using γ -Re θ . In conclusion, the authors believe that γ is very promising, but at the current stage of transition models validation, if there are no experiments to compare with, γ -Re θ might be more reliable.

3-4. The authors believe that none of the discussed criteria is good enough and there is still a huge room for improvement. Unfortunately, this is not possible without any further experiments, that would help to have additional insight into the different transition mechanisms. For instance, many of the correlations used in the LCTM methods originate from correlations derived on a flat plate, with adverse/favorable pressure gradients generated by the shape of the outer wall of the test section. This is a notably different physical setup from the geometries we are interested in. It was also raised the problem that streamwise criteria based on characteristic boundary layer thickness might fail on complex three-dimensional configurations where the streamlines are not parallel, as is the case for the flow at the symmetry plane at the windward side of the 6:1 prolate spheroid. Additional data may also provide more insight to better characterize the mechanism for transition. Measurements of velocity profile, pressure distribution, actual turbulence intensity, and surface roughness would be useful quantities to both identify the transition mechanism, to further improve the existing transition criteria, and provide validation data.

Furthermore, better quantification of the experimental uncertainty may recontextualize the perceived accuracy of the CFD results.

To this matter, an additional paragraph was added to the conclusions to outline possible useful experiments around the 6:1 prolate spheroid that could help developers improve transition models. At the end, the authors stress the possibility of including the interaction between the transition mechanisms, as it is one of the few problems that might be resolved soon, without any specific additional experiments. This is a well-written and thorough paper that presents an improved model for predicting crossflow transition on complex three-dimensional bodies. The authors' enhancement of Menter and Smirnov's Tc1 crossflow criterion has significantly improved the crossflow transition prediction for the 6:1 prolate spheroid at an angle of attack. The geometry and flow conditions are relevant to unmanned underwater vehicles so the improved model is of practical value to the naval hydrodynamics community. The authors should be congratulated for this valuable contribution. As the authors mention, not many past experiments with the 6:1 prolate spheroid have focused on the laminar-to-turbulent transition. Only the database of Kreplin et al. from 1985 was found to be suitable for the model validation. Unfortunately, there is a large uncertainty in the experimental free-stream turbulence intensity, Tu, which is an important parameter for transition. It was mentioned that Tu varies from 0.1% to 0.3% close to the nose of the body in the experiment and the authors selected a leading-edge Tu of 0.15% for the computations. Question 1. How sensitive is the transition location and the skin friction to free-stream turbulence for this case, in the experimental uncertainty range of 0.1% Tu 0.3%? Have the authors considered varying the free-stream turbulence intensity in the computations to the see the sensitivity of the model prediction? As the free-stream turbulence is increased, will the model eventually predict the transition on the windward side that was seen in the experiments but missing in the computations, and could this occur within the bound of experimental uncertainty (i.e., Tu 0.3%)? Question 2.

DISCUSSION

The authors identify several additional physical transition mechanisms that are currently not included in the transition models that could explain the failure to predict transition on the windward side. Do the authors have suggestions for future experiments that could help further the development of transition models? Question 3. Finally, it is mentioned that iterative convergence is challenging for the simulations. Compared to RANS computations without transition, it looks like a large number of CPU hours are needed. Did the authors' modifications to the Tc1 criterion affect computational cost significantly? Do the authors have any ideas for ways to improve iterative convergence and reduce computational cost? AUTHOR'S REPLY 1. The value Tu = 0.15% was chosen within the AVT-313 group in occasion of the workshop dedicated to the 6:1 prolate spheroid at incidence. It corresponds to the turbulence intensity at which the critical N factor for T-S waves, N T S , is encountered, as computed in the stability diagram used by DLR for their e N computations around the 6:1 prolate spheroid at 15 • incidence. While DLR e N method gives very good results for this initial condition this is not the case for any LCTM models which are not able to predict transition near the windward symmetry plane. This non-predicted flow feature is hardly due to the the turbulence level used in the computations and it is generally accepted that an higher value of turbulence intensity would not lead to the transition prediction on the windward side of the prolate spheroid. To this matter, only one simulation for an higher turbulence level, i.e. Tu = 0.3%, was run using γ -Re θ model in ISIS (results are not shown hereafter). No significant change was observed in the transition front position, other than a local higher value of the skin friction coefficient. These higher values, however, do not reflect the real physics of the problem. Indeed, the higher C f is due to the LCTM T-S criterion, rather than the crossflow ones. In fact, none of the transition criteria used in this research depends on the turbulence intensity level in the boundary layer, as they are meant to predict stationary crossflow waves, that are not really sensitive to this parameter. In the present context, the observation about the inconsistency between numerical and experimental set-ups goes beyond the presented quantitatively comparisons. Indeed, the problem is that for higher values of incoming turbulence intensities the transition mechanisms might differ from those considered in the research, see traveling crossflow, for instance. The unknown experimental conditions do not allow performing of a rigorous validation exercise, complicating the identification of the prediction errors sources.

2. The authors added some suggestions for future experiments for the improvement of the transition models at the end of "Conclusions".

3.

In general, as for two-dimensional simulations, γ converges faster than γ -Re θ . Nevertheless, the computation of the hessian of the velocity degrades to some extent γ model convergence. This 2nd derivative of the velocity is computed as it follows: 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 region were the criterion activates. The calculation of a 2nd derivative within the boundary layer is troublesome and requires a mesh enough refined in both streamwise and spanwise directions. Nevertheless, these stringent requirements on the mesh resolution are the same demanded by the transition models for accurate predictions with negligible discretization error. At the end, in terms of computational time, the CPU needed by γ + T c1, despite the computation of the velocity hessian, are still less than γ -Re θ + CF, i.e. 9000 CPU vs. 13000 CPU, as added in the paragraph "Computational Domain". Undoubtedly, considering the cost of these simulations in terms of both wall-clock time and computational resources, some effort should be directed to determining what are the actual requirements on mesh resolution for obtaining accurate results. It is possible that grids of 40-50 million cells are not really needed to accurately capture the transition front, but similar results can be obtained with less dense meshes by better targeting certain regions for refinement. A possible strategy would be to define an adaptive grid refinement strategy, where the refinement dynamically follows the transition front, where an higher mesh density is needed to capture the large gradients of the intermittency, γ. Unfortunately, current trends in gridding for industrial applications may fail to capture large gradients of γ in the core of the boundary layer, as well as computing accurately the hessian of the velocity, but the possible consequences are still not clear at this stage. 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 Piotrowsky has proposed a smoothing approximations for these functions by exponential penalty functions, as proposed in [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.

Indeed, being γ and γ -Re θ models based on correlations built empirically, every change in their formulation can affect transition predictions.

Figure 2 :

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

  Fig.(3) shows the uncorrected ratio R = 0.684Ψ/(XF(H 12 )) with the cubic polynomial g(λ θ ,CF ), Fig.(3a), and the corrected R = (G MS (λ θ ,CF )Ψ)/(XF(H 12 )), Fig.(

  Fig.(4).

( a )

 a Uncorrected indicator ratio R with g(λ θ ,CF ). (b) Corrected indicator ratio R.

Figure 3 :

 3 Figure 3: Crossflow indicators ratio R distribution, projected on the β -plane and plotted as a function of the new λ θ ,CF , for family of FSC profiles. Planes are colored by the contours of the sweep angle φ . Top: uncorrected ratio R = (0.684Ψ)/(XF(H 12 )), black line g(λ θ ,CF ). Bottom: corrected ratio R = (G MS (λ θ ,CF )Ψ)/(XF(H 12 )).

Figure 4 :

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

λ

  CF =0.0473 -0.0001338 φ -

Figure 5 :

 5 Figure 5: R = G new Ψ/(XF(H 12 )) vs the fitted pressure gradient parameter λ CF dv dy y 2 ν , φ and the sweep angle φ .

  γ +CF kω SST (2003) γ Recalibrated T c1 with φ L and G new γ + T c1 -MS kω SST (2003) γ T c1 MS with G MS γ -Re θ +CF kω SST (2003)γ -Re θ He β = 0.0828 and β * = 0.09 are constants. Eq.(94) is obtained by the resolution of the kω (2003) SST transport equations for a uniform, steady flow, aligned with x, and neglecting the diffusion and cross-diffusion terms. For the presented simulations, the values imposed at the inlet, * in , are reported in Table(2).

  (a) Full Computational Domain.(b) O-block.

Figure 6 :

 6 Figure 6: 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: 6:1 prolate spheroid: α = 15 • . Coarsest mesh in the wall normal direction at the section X/L = 0.

Figure 8 :

 8 Figure 8: 6:1 prolate spheroid. Coarsest mesh at α = 15 • , closed up on the first half of the spheroid: surface mesh (black), mesh at the symmetry plane (red).

Figure 9 :

 9 Figure 9: 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, relying on Richardson extrapolation.

Figure 10 :

 10 Figure 10: 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, relying on Richardson extrapolation.

  ), as for the presented variant γ +CF. The C f distribution as computed by γ +T c1-MS model is shown in Fig.(11c).

( a )

 a Measured C f . (b) Predicted C f by γ -Re θ +CF. (c) Predicted C f by γ + c1 -MS. (d) Predicted C f by γ +CF.

Figure 11 :

 11 Figure 11: 6:1 Prolate Spheroid: α = 15 • . Measured and predicted skin friction C f distributions in the X/Lφ plane. Numerical results are computed using γ + CF and γ -Re θ + CF models. The results obtained with γ +T c1-MS, with the T c1 version of Menter & Smirnov, the T c1 MS , are also shown.

  139 and X/L = 0.223, Fig.(12a) and Fig.(12b). Actually, the kink is observed in the experiments at the section X/L = 0.223, Fig.(12b), which is downward the kink location predicted by γ -Re θ + CF. γ + CF predicts an uniform transition front. This behavior is due to the fact that C1-based criteria have the tendency to predict transition upstream when associated with T-S criteria, Bégou (2018). In general, both γ -Re θ +CF and γ +CF predict transition upstream with respect to the experiments. Close to the leeward symmetry plane, independently of the transition model, it can be observed a tongue of delayed transition. The laminar tongue appears at section X/L = 0.223 and it is visible until half of the prolate spheroid length, Fig.(12c) and Fig.(

  • . The flow close to the windward symmetry plane remains laminar until the rear side of the spheroid. Measurements and predictions start deviating at the section X/L = 0.652, Fig.(12e), up to the last section X/L = 0.936, Fig.(12f). None of the crossflow criteria, the recalibrated T c1 nor the He-based, is active in this region of the prolate spheroid.

  (a) X/L = 0.139. (b) X/L = 0.223. (c) X/L = 0.309. (d) X/L = 0.480. (e) X/L = 0.652. (f) X/L = 0.936.

Figure 12 :

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

  ). Indeed, considering the predicted transition front, computations on this grid are sufficiently accurate. For α = 5 • predictions by γ + CF model are also compared to the results obtained with γ -Re θ + CF. Hereafter, the results computed with the original criterion by Menter & Smirnov, γ + T c1 -MS, are no more considered.

Figure 13 :

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

  (a) Measured C f . (b) Predicted C f by γ -Re θ +CF. (c) Predicted C f distribution by γ + CF.

Figure 14 :

 14 Figure 14: 6:1 Prolate Spheroid: α = 5 • . Measured and predicted skin friction C f distributions in the X/Lφ plane. Numerical results are computed using γ +CF and γ -Re θ +CF models.

Figure 15 :

 15 Figure 15: 6:1 Prolate Spheroid: α = 5 • . Experimental and numerical girthwise distribution of the C f coefficient at different section. Numerical results are computed with γ + CF and γ -Re θ + CF models. At each validation point, the experimental uncertainty is reported.

Figure 16 :

 16 Figure 16: 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.

  (a) Measured C f . (b) Predicted C f by γ +CF.

Figure 17 :

 17 Figure 17: 6:1 Prolate Spheroid: α = 30 • . Measured and predicted skin friction C f distribution on the 6:1 prolate spheroid. Numerical results are computed using γ + CF model.

  Fig.(19), that cause the failure of the streamwise criterion, i.e. the differences between the physical distance along which the waves propagate and the distance along which the characteristic boundary layer thickness grows. The unpredicted laminar separation bubble pollutes further downstream the numerical solution. The flow at the leeward side numerically starts transitioning around X/L = 0.139, but the skin friction is systematically underestimated with respect to measurements, Fig.(18c) to Fig.(18f).

  (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 18 :

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

Figure 19 :

 19 Figure 19: 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.

  Discusser: Mark Bettle, Defence Scientist, Defence Research and Development Canada (DRDC)

  Dividing the vorticity Reynolds number by 2.193Re θ c ensures Re V to have a maximum of one within the boundary layer. Re θ c is an empirical correlation, determined with respect to numerical flat plate experiments. The correlation is function of the solution of the second transport equation Re θ t . F onset,2 , Eq.(

  that the non-local empirical quantity Re θ t is diffused within the boundary layer from the free-stream. The solution of this transport equation is then the local Reynolds number Re θ t . The empirical correlation for Re θ t is based on the pressure gradient parameter λ θ , defined as:

Table 2 :

 2 6:1 Prolate Spheroid: Computations details.

	5 •
	15 • 6.5 × 10 6	0.5	250	0.15
	30 •

α

Re Tu in (%) (ν t /ν) in Tu(%)

Table 3 :

 3 6:1 Prolate Spheroid: Mesh details.

	N cells	N surface h i /h 1 y + max
	Grid1 42.6M 126016	1	0.4
	Grid2 28.3M 95816	1.14 0.46
	Grid3 17.9M 70884	1.33 0.54
	Grid4 10.3M 48750	1.61 0.65
	Grid5 5.3M	31504	2	0.8

Other experimental campaigns for crossflow transition were presented in[START_REF] Ahn | An experimental study of flow over a 6 to 1 prolate spheroid at incidence[END_REF], nevertheless for the untripped geometry the effects of gravity acting on the oil mixture did not allow to show distinct flow patterns.

If we exclude the stagnation point, the windward symmetry plane looks like the attachment line of a swept cylinder with geometrical sweep angle Φ = π 2α, where α is the angle of attack, as explained in[START_REF] Arnal | Three-dimensional boundary layer: laminar-turbulent transition[END_REF].

ACKNOWLEDGMENTS

The computations were performed using HPC resources from GENCI (Grand Equipement National de Calcul Intensif) (Grant-A0102A00129), which is gratefully acknowledged. The help of Dr. Rui Lopes from IST Lisbon who provided us with the structured grids for the 6:1 prolate spheroid is warmly acknowledged. Finally, Prof. Christophe Airiau from Paul Sabatier University in Toulouse is sincerely acknowledged for sharing his Falkner-Skan-Cooke solver with the authors.