Charlotte Her

| Hamid-Reza Rezaei

Sandrine Hughes

Saeid Naderi

Marilyne Duffraisse

| Marjan Mashkour

Hamid-Reza Naghash

Adrian Bălășescu

Gordon Luikart

Steve Jordan

Deniz Özüt

Aykut Kence

Michael W Bruford

| Anne Tresset

Jean-Denis Vigne

Pierre Taberlet

Catherine Hänni

François Pompanon
email: francois.pompanon@univ-grenoble-alpes.fr

Leca

Broad maternal geographic origin of domestic sheep in Anatolia and the Zagros

We investigated the controversial origin of domestic sheep (Ovis aries) using large samples of contemporary and ancient domestic individuals and their closest wild relatives: the Asiatic mouflon (Ovis gmelini), the urial (Ovis vignei) and the argali (Ovis ammon). A phylogeny based on mitochondrial DNA, including 213 new cytochrome-b sequences of wild Ovism confirmed that O. gmelini

The earliest archaeological evidence for the presence of domesticated sheep was found in the Taurus mountains, Southeastern Anatolia and dates to 10,500 years BP (International Council for Archaeozoology et al., 2005;Zeder, 2008), making this region a putative area of sheep domestication. Sheep later spread over the Old World during the Neolithic, following human migration routes. Toward Europe, they took two major trajectories: the Danubian route across continental Europe and the Mediterranean route, crossing the sea to large islands such as Cyprus and Sardinia from the ninth to the sixth millennia BP [START_REF] Tresset | Going over: the mesolithic-neolithic transition in North-West Europe[END_REF], 2011;[START_REF] Vigne | First wave of cultivators spread to Cyprus at least 10,600 y ago[END_REF]. Toward Asia, the most likely hypothesis is colonization through Turkmenistan and Mongolia [START_REF] Flad | Zooarcheological evidence for animal domestication in northwest China[END_REF][START_REF] Lv | Mitogenomic meta-analysis identifies two phases of migration in the history of eastern Eurasian sheep[END_REF][START_REF] Vigne | Early domestication and farming: what should we know or do for a better understanding[END_REF]. In the eastern part of the spread, the earliest evidence of sheep is to be found at around 9500-9000 years BP in central and southern Zagros [START_REF] Daujat | Faunal remains from middle Neolithic site of Qaleh Rostam[END_REF][START_REF] Kamjan | The archaeozoology of the Pre-Pottery and Pottery Neolithic site of Qasr-e Ahmad (Iran)[END_REF][START_REF] Roustaei | The Neolithic of the Iranian plateau. Recent research and prospects[END_REF]) and t about the same dates in the north-east of Iran [START_REF] Roustaei | Tappeh Sang-e Chakhmaq and the beginning of the Neolithic in north-east Iran[END_REF] and southern central Asia [START_REF] Dobney | The vertebrate assemblage from excavations at Jeitun, 1993 and 1994[END_REF]. The presence of domestic sheep is documented much later in China, around 4500-4000 years BP [START_REF] Flad | Zooarcheological evidence for animal domestication in northwest China[END_REF].

Based on genetic studies, several wild Asiatic species have been proposed as the ancestor of domestic sheep [START_REF] Nadler | G-band patterns as chromosomal markers, and the interpretation of chromosomal evolution in wild sheep (Ovis)[END_REF][START_REF] Pedrosa | Evidence of three maternal lineages in near eastern sheep supporting multiple domestication events[END_REF]. These are the argali, Ovis ammon, the Asiatic mouflon, Ovis gmelini (often previously named orientalis following the nomenclature of [START_REF] Nadler | G-band patterns as chromosomal markers, and the interpretation of chromosomal evolution in wild sheep (Ovis)[END_REF] and the urial, Ovis vignei. Cytogenetic [START_REF] Nadler | G-band patterns as chromosomal markers, and the interpretation of chromosomal evolution in wild sheep (Ovis)[END_REF] and archaeozoological [START_REF] Clutton-Brock | Domesticated animals from early times[END_REF][START_REF] Uerpmann | The ancient distribution of Ungulate Mammals in the Middle East[END_REF] data have suggested the elimination of O. ammon and O. vignei as potential ancestors, leaving O. gmelini as the best candidate wild species from which sheep were domesticated. However, fluctuations in the nomenclature [START_REF] Rezaei | Evolution and taxonomy of the wild species of the genus Ovis (Mammalia, Artiodactyla, Bovidae)[END_REF], the occurrence of a hybrid zone between O. gmelini and O. vignei and putative domestication events in the lower Indus valley [START_REF] Meadows | Early herders and their flocks[END_REF] raise the possibility of genetic contributions from other wild species. Thus, we aimed to confirm the hypothesis of a unique maternal O. gmelini origin, and to localize the domestication area by finding the present wild populations representing the most direct descendants of domesticated matrilineages. We also sought to identify the first maternal lineages of domestic sheep entering Europe during the Neolithic spread.

Mitochondrial DNA (mtDNA) has been extensively and successfully used to describe the genetic diversity of domestic animals and to study their origins and histories. Indeed, it allows the identification of the maternal evolutionary relationships with wild relatives because it is evolutionarily well conserved, yet sufficiently variable to capture genetic structure [START_REF] Bruford | DNA markers reveal the complexity of livestock domestication[END_REF][START_REF] Zeder | Documenting domestication: new genetic and archaeological paradigms[END_REF]. Fragments of the non-coding Control Region (CR) or Cytochrome b (Cytb) gene have been used to show the existence of five different haplogroups in domestic sheep, designated as A-E [START_REF] F I G U R E 2 ; Arora | Phylogenetic relationships between the domestic sheep (Ovis aries) and its potential wild ancestors: the argali (Ovis ammon), the urial (Ovis vignei) and the Asiatic mouflon (Ovis gmelini). The phylogenetic tree was built using maximum likelihood methods implemented in iqtree with cytb sequences (1217 pb). The numbers on the principal nodes correspond to support values (Ultrafast Bootstrap value/SH-aLRT value). The black circle arcs represent the different mitochondrial haplogroups[END_REF][START_REF] Chen | Origin, genetic diversity, and population structure of Chinese domestic sheep[END_REF][START_REF] Demirci | Mitochondrial DNA diversity of modern, ancient and wild sheep (Ovis gmelinii anatolica) from Turkey: new insights on the evolutionary history of sheep[END_REF][START_REF] Guo | A novel maternal lineage revealed in sheep (Ovis aries)[END_REF][START_REF] Hermes | Mitochondrial DNA of domesticated sheep confirms pastoralist component of Afanasievo subsistence economy in the Altai Mountains (3300-2900 cal BC)[END_REF][START_REF] Hiendleder | Molecular analysis of wild and domestic sheep questions current nomenclature and provides evidence for domestication from two different subspecies[END_REF][START_REF] Horsburgh | Genetic characterization of an archaeological sheep assemblage from South Africa's Western Cape[END_REF][START_REF] Ibrahim | The genetic profiles and maternal origin of local sheep breeds on Java Island (Indonesia) based on complete mitochondrial DNA D-loop sequences[END_REF][START_REF] Kolosov | Investigation of the D-loop sequence of mitochondrial DNA of the Volgograd sheep breed[END_REF][START_REF] Meadows | Five ovine mitochondrial lineages identified from sheep breeds of the near East[END_REF][START_REF] Pedrosa | Evidence of three maternal lineages in near eastern sheep supporting multiple domestication events[END_REF][START_REF] Resende | Mitochondrial DNA variation of domestic sheep (Ovis aries) in Kenya[END_REF][START_REF] Rezaei | Evolution and taxonomy of the wild species of the genus Ovis (Mammalia, Artiodactyla, Bovidae)[END_REF][START_REF] Singh | Extensive variation and sub-structuring in lineage A mtDNA in Indian sheep: genetic evidence for domestication of sheep in India[END_REF][START_REF] Tapio | Sheep mitochondrial DNA variation in European, Caucasian, and Central Asian areas[END_REF][START_REF] Zhao | Mitochondrial DNA diversity, origin, and phylogenic relationships of three Chinese large-fat-tailed sheep breeds[END_REF]. Analysis of complete mitochondrial genomes [START_REF] Meadows | Haplogroup relationships between domestic and wild sheep resolved using a mitogenome panel[END_REF] confirmed the validity of such mtDNA fragments for discriminating among haplogroups and tracing the history of maternal lineages.

We built two datasets, one to compare modern domestic and wild sheep and another to compare modern and ancient domestic sheep (see Table S1). For the first one, we carried out an extensive non-invasive sampling of wild Ovis corresponding to 132 wild sheep, including three O. ammon, 85 O. gmelini, 15 O. vignei, one Ovis musimon and 28 hybrids gmelini × vignei. Samples originated from 55 geographic localities representing most of their distribution area and encompassing the putative region of domestication according to archaeological data (see Table S2, Appendix S1). Moreover, we collected samples from 81 domestic sheep from 28 different countries as representative of the mtDNA worldwide domestic diversity (Table S2, Appendix S1). We sequenced 1217 bp of the Cytb locus for these samples and retrieved 421 supplementary sequences from GenBank (Table S2, Appendix S1). We built the second dataset on 120 bp of the more variable CR sequence to analyze intra-specific diversity (corresponding to the most recent divergences). Data for 1283 present-day sheep were collected from GenBank, and we sequenced 54 modern O. gmelini (Table S3) and 57 ancient domestic samples mostly distributed in Europe and dating back around 8000 and 5000 years BP (Figure 1b, Tables S4 andS5, Appendix S1). We built maximum likelihood and neighbor-joining [START_REF] Saitou | The neighbor-joining method: a new method for reconstructing phylogenetic trees[END_REF] trees using iqtree [START_REF] Nguyen | IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF] and the R package ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF], respectively (Appendix S1), to infer the phylogenetic relationships between the Middle East wild species and 362 Ovis aries representative of worldwide domestic diversity. Both phylogenies showed the monophyly of each of the three wild species (i.e. O. gmelini, O. vignei and O. ammon). All O. aries fall within the O. gmelini clade, confirming that domestic sheep evolved from ancestral matrilineages of the Asiatic mouflon [START_REF] Bunch | Phylogenetic analysis of snow sheep (Ovis nivicola) and closely related taxa[END_REF][START_REF] Demirci | Mitochondrial DNA diversity of modern, ancient and wild sheep (Ovis gmelinii anatolica) from Turkey: new insights on the evolutionary history of sheep[END_REF][START_REF] Hiendleder | The complete mitochondrial DNA sequence of the domestic sheep (Ovis aries) and comparison with the other major ovine haplotype[END_REF][START_REF] Hiendleder | Molecular analysis of wild and domestic sheep questions current nomenclature and provides evidence for domestication from two different subspecies[END_REF][START_REF] Meadows | Five ovine mitochondrial lineages identified from sheep breeds of the near East[END_REF][START_REF] Meadows | Haplogroup relationships between domestic and wild sheep resolved using a mitogenome panel[END_REF][START_REF] Tapio | Sheep mitochondrial DNA variation in European, Caucasian, and Central Asian areas[END_REF]. This latter clade was divided into two clusters, both containing haplotypes from wild and domestic sheep (Figure 2), showing the early divergence between A/B haplogroups and C/E haplogroups [START_REF] Meadows | Five ovine mitochondrial lineages identified from sheep breeds of the near East[END_REF]. Interestingly, individuals from the O. gmelini × vignei hybrid zone in central Iran had either an O. gmelini or an O. vignei maternal origin, as already shown by [START_REF] Demirci | Mitochondrial DNA diversity of modern, ancient and wild sheep (Ovis gmelinii anatolica) from Turkey: new insights on the evolutionary history of sheep[END_REF]. None of them fell into a clade containing domestic sheep (Figure 2), precluding the contribution of hybrid populations to domestics. However, urial Y-chromosome haplotypes, found in several Asiatic mouflons [START_REF] Deng | Paternal origins and migratory episodes of domestic sheep[END_REF], support a partial paternal contribution of the urial.

Assuming that the boundaries of the current distribution of wild species are similar to those at the time of domestication [START_REF] Uerpmann | The ancient distribution of Ungulate Mammals in the Middle East[END_REF], the exclusion of O. ammon, O. vignei and O. vignei × O. gmelini hybrids as possible ancestors precludes hypotheses that the Central Iranian Plateau, the Lower Indus Valley and even more eastern locations contributed to domestication. Based on archaeozoological data, the cradle for sheep domestication would extend from central Anatolia to northern and central Zagros (International Council for Archaeozoology et al., 2005;Zeder, 2008Zeder, , 2011)). This is confirmed by mtDNA data, as the highest nucleotide diversity was found in this region (Table S6) where, moreover, the mtDNA haplogroups carried by modern Asiatic mouflon were the closest to those of the domestics (Figure 1b′). This region overlaps the current distribution of the subspecies O. gmelini gmelini, whereas O. g. isphahanica and O. g. laristanica are spread further east and south, respectively (Figure 1a; Michel & Ghoddousi, 2020). The strong bottleneck experienced (even recently) by the Cyprus mouflon (Ovis ophion), which is close to O. g. gmelini [START_REF] Sanna | The first mitogenome of the Cyprus mouflon (Ovis gmelini ophion): new insights into the phylogeny of the genus Ovis[END_REF] and results from the early transportation of individuals from the continent [START_REF] Vigne | The transportation of mammals to Cyprus sheds light on early voyaging and boats in the Mediterranean Sea[END_REF], makes it very difficult to assess its potential contribution to the domestic diversity. Besides that, the presence of some haplotypes similar to those of domestic sheep in modern O. g. gmelini populations offers evidence that the latter gave rise to the domestic haplogroups without excluding later introgression of domestic stock into wild populations.

We used CR sequences from 1283 modern and 57 ancient domestic sheep to investigate global sheep mitochondrial diversity. Most modern sheep (56%) were from the B haplogroup, followed by the A (34%), C (9%) and D and E (<0.5%) (Figure S2). These proportions are similar to those of other studies [START_REF] Demirci | Mitochondrial DNA diversity of modern, ancient and wild sheep (Ovis gmelinii anatolica) from Turkey: new insights on the evolutionary history of sheep[END_REF][START_REF] Kandoussi | The origin of sheep settlement in Western Mediterranean[END_REF][START_REF] Lv | Mitogenomic meta-analysis identifies two phases of migration in the history of eastern Eurasian sheep[END_REF][START_REF] Meadows | Five ovine mitochondrial lineages identified from sheep breeds of the near East[END_REF][START_REF] Sanna | The first mitogenome of the Cyprus mouflon (Ovis gmelini ophion): new insights into the phylogeny of the genus Ovis[END_REF]. As already shown [START_REF] Handley | Genetic structure of European sheep breeds[END_REF][START_REF] Kijas | A genome wide survey of SNP variation reveals the genetic structure of sheep breeds[END_REF][START_REF] Meadows | Mitochondrial sequence reveals high levels of gene flow between breeds of domestic sheep from Asia and Europe[END_REF], sheep showed very weak genetic differentiation at the continental scale (no between-continent F ST >0.02). The A and B haplogroups have a worldwide distribution that is not homogeneous, with more than 78% of the A haplogroup located in Asia. In northern and southern Europe, A and B haplogroups were detected with a high preponderance of the B haplotypes (above 90%) and some rare C haplotypes (Table S6; [START_REF] Demirci | Mitochondrial DNA diversity of modern, ancient and wild sheep (Ovis gmelinii anatolica) from Turkey: new insights on the evolutionary history of sheep[END_REF][START_REF] Deng | Paternal origins and migratory episodes of domestic sheep[END_REF][START_REF] Lv | Mitogenomic meta-analysis identifies two phases of migration in the history of eastern Eurasian sheep[END_REF].

The genotyping of 57 ancient samples showed that there were at least three haplogroups present in Neolithic European domestic sheep, with the B haplogroup already predominant (i.e. 55 B, 1 A and 1 D haplotypes; Appendix S1, Tables S5 andS6). We probably did not detect other low-frequency haplogroups from that area. However, this result indicates that the A and D haplogroups have been present since the beginning of western sheep diffusion during the Neolithic. This predates by several millennia the first evidence of the A haplogroup in Europe, as it was already detected in medieval sheep from Finland but not in the oldest samples from the Iron period [START_REF] Demirci | Mitochondrial DNA diversity of modern, ancient and wild sheep (Ovis gmelinii anatolica) from Turkey: new insights on the evolutionary history of sheep[END_REF][START_REF] Niemi | Mitochondrial DNA and Ychromosomal diversity in ancient populations of domestic sheep (Ovis aries) in Finland: comparison with contemporary sheep breeds[END_REF][START_REF] Rannamäe | Maternal and paternal genetic diversity of ancient sheep in Estonia from the Late Bronze Age to the post-medieval period and comparison with other regions in Eurasia[END_REF]. Thus we confirm the hypothesis of a very early presence of this haplogroup in Europe, and support the idea that a later introduction of A haplotypes from Spanish populations during the sixteenth century is unlikely as these importations concerned mainly rams [START_REF] Niemi | Mitochondrial DNA and Ychromosomal diversity in ancient populations of domestic sheep (Ovis aries) in Finland: comparison with contemporary sheep breeds[END_REF].

Combining genetic data from O. aries and its closest wild relatives with archaeozoological data makes it possible to propose a subspecies of the Asiatic mouflon, O. gmelini gmelini, as the maternal origin of domestic sheep. The current distribution of the Asiatic mouflon and the geographic pattern of its genetic diversity are consistent with domestication first taking place in eastern Anatolia and in northern and central Zagros, as supported by archaeozoological data [START_REF] Clutton-Brock | Domesticated animals from early times[END_REF]International Council for Archaeozoology et al., 2005;[START_REF] Uerpmann | The ancient distribution of Ungulate Mammals in the Middle East[END_REF]. They preclude more eastern independent domestication centers, especially in the Indus Valley, as well as the contribution of O. gmelini × O. vignei hybrids. Studies on the diversity of chromosome Y have given valuable insights into the demographic history of sheep [START_REF] Deng | Paternal origins and migratory episodes of domestic sheep[END_REF], but have not provided information on the geographic origin. However, studies on chromosome Y [START_REF] Li | A Hu sheep genome with the first ovine Y chromosome reveal introgression history after sheep domestication[END_REF] and whole genome data [START_REF] Barbato | Genomic signatures of adaptive introgression from European mouflon into domestic sheep[END_REF] F I G U R E 1 Geographic distribution of wild Ovis species (a) and ancient domestic (b) and modern wild (b′) haplogroups. We obtained species distributions from the IUCN database (see references below). For Ovis gmelini subspecies, we combined information from sites described in the detailed report and general distribution of the species to draw approximate ranges of continental subspecies and underwent early diffusion toward Europe at the beginning of the Neolithic. Thus, our results globally indicate that sheep domestication was a process that occurred on a large geographical scale. It involved the capture of a wide genetic diversity over a large area, as already shown for other livestock species such as goats [START_REF] Daly | Ancient goat genomes reveal mosaic domestication in the Fertile Crescent[END_REF][START_REF] Naderi | The goat domestication process inferred from large-scale mitochondrial DNA analysis of wild and domestic individuals[END_REF], pigs [START_REF] Bosse | A genomics perspective on pig domestication[END_REF][START_REF] Larson | Ancient DNA, pig domestication, and the spread of the Neolithic into Europe[END_REF][START_REF] Larson | Patterns of East Asian pig domestication, migration, and turnover revealed by modern and ancient DNA[END_REF] and cattle [START_REF] Götherström | Cattle domestication in the Near East was followed by hybridization with aurochs bulls in Europe[END_REF][START_REF] Mctavish | New world cattle show ancestry from multiple independent domestication events[END_REF].

AC K NOW L E DGE M E N T S

We would like to thank Benjamin Gillet and Bénédicte Bertrand for help with the ancient DNA experiments. We also are grateful to all of the archaeozoologists who helped to obtain the ancient sheep samples, in particular to Daniel Helmer, Lenka Kovacikova, Rose-Marie Arbogast and Joris Peters. We would also like to thank our fellow archaeologists who allowed us access to the archaeozoological materials from Armenia (Ruben Badalyan and Christine Chataigner) and Romania (Radian Andreescu, Pavel Mirea, Dragomir Popovici, Valentin Parnic, Cătălin Bem and Florin Vlad). This study has been supported by a grant from the French Agence Nationale de la Recherche (Project Chronobos, ANR 05-GANI-004-02).

CON F L IC T OF I N T E R E ST S

The authors declare no conflict of interests.

DATA AVA I L A BI L I T Y STAT E M E N T

All Cytb and CR sequences produced in this study have been published on GenBank under the accession numbers MW684398-MW684451, MW684452-MW684508 and MW677233-MW677445.

ORC I D

Charlotte Her https://orcid.org/0000-0001-6418-6201

R E F E R E NC E S

  is the maternal ancestor of sheep and precluded mtDNA contributions from O. vignei (and O. gmelini × O. vignei hybrids) to domestic lineages. We also produced 54 new control region sequences showing shared haplogroups (A, B, C and E) between domestic sheep and wild O. gmelini which localized the domestication center in eastern Anatolia and central Zagros, excluding regions further east where exclusively wild haplogroups were found. This overlaps with the geographic distribution of O. gmelini gmelini, further suggesting that the maternal origin of domestic sheep derives from this subspecies. Additionally, we produced 57 new CR sequences of Neolithic sheep remains from a large

  area covering Anatolia to Europe, showing the early presence of at least three mitochondrial haplogroups (A, B and D) in Western colonization routes. This confirmed that sheep domestication was a large-scale process that captured diverse maternal lineages (haplogroups).K E Y W O R D Sancient DNA, animal domestication, area of early domestication, maternal genetic origins, mitochondrial DNA, phylogeny, sheep ancestor, wild sheep

  brought evidence for the occurrence of post-domestication introgressions from the wild toward domestic sheep. At least three different hap-(B, A and D) diverged before domestication

Funding information

Project Chronobos, Grant/Award Number: ANR 05-GANI-004-


Ginevra Rubino
email: ginevra.rubino@ec-nantes.fr

Michel Visonneau
email: michel.visonneau@ec-nantes.fr

Improved Crossflow Transition Predictions for the one-equation γ transition model

Keywords: RANS, Transition Modeling, Stationary Crossflow, Falkner-Skan-Cooke, γ model

A variant of the local correlation-based transition model (LCTM) γ is presented to account for the prediction of transition due to stationary crossflow instabilities. The presented T c1 crossflow criterion is a local approximation of the well-known C1 criterion by Arnal (1984), evaluated locally through auxiliary functions expressed by the solution of the Falkner-Skan-Cooke (FSC) equations. The criterion, originally proposed by Menter & Smirnov (2014), uses the wall-normal directional change of the normalized vorticity as indicator of the crossflow strength. Hereafter, an original calibration FSC-based of Menter & Smirnov criterion is proposed. A local approximation of the sweep angle is included in order to achieve better results on non-wing-like geometries. The capability of this model variant to predict stationary crossflow transition is therefore discussed. A validation study is presented using experimental data on the 6:1 prolate spheroid and the sickle wing.

Introduction

Nowadays understanding, predicting and controlling laminar-to-turbulence transition is one of the biggest challenges in Computational Fluid Dynamics (CFD). There are several practical applications that deal with low/moderate Reynolds numbers: 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 them. In the last decade, the local-correlation transition model (LCTM) concept has known a huge success. The RANS transition models γ -Re θ , proposed by 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], are today widely spread. They are preferred in practical applications to Direct Numerical Simulations (DNS) and Large Eddy Simulations (LES) approaches, because of the reduced computational costs. The aim of these RANS transition models is not to describe the 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. γ -Re θ model, officially published in 2012, is based on the solution of two transport equations: one for the intermittency γ and one for the momentum thickness Reynolds number Re θ t . The quantity γ is used to turn on the production of turbulence kinetic energy in the boundary layer downstream the transition onset location. The equation for Re θ t is needed to transport non-local information, namely the empirical momentum thickness Reynolds number Re θ t , from the free-stream inside the boundary layer. The solution of this second transport equation, Re θ t , appears in the definition of the critical momentum thickness Reynolds number, Re θ c . The latter is required in the empirical correlation that triggers the production of intermittency in the boundary layer. Within the NATO AVT-313 "Incompressible Laminar-to-Turbulent Flow Transition Study" collaboration group, γ -Re θ and γ models performances have been tested on several configurations, both two-and threedimensional. During these studies, it was observed that the two models perform very similarly in terms of transition predictions. Nevertheless, the one-equation γ model was found easier and faster to converge compared to γ -Re θ (see Lopes, [START_REF] Lopes | Simulation of transition from laminar to turbulent regime in practical applications of incompressible flow[END_REF], Rubino, [START_REF] Rubino | Laminar-to-turbulence transition modeling of incompressible flows in a rans framework for 2d and 3d configurations[END_REF]). For this reason, in the present discussion attention is focused on the one-equation γ model and its extension to crossflow transition predictions. Indeed, as originally presented, both γ-Re θ and γ models only accounted for bypass and streamwise natural transition, as due to Tollmien-Schlichting (T-S) waves. In the last six years, different research groups have devoted their effort to the definition of a criterion to account for transition due to stationary crossflow waves (CF). These are the dominant instabilities modes on three-dimensional configurations operating in very low free-stream turbulence environments.

An empirical correlation that accounts for crossflow transition should be a function of the multiple pa-rameters that play a role in the process: the pressure gradient in the streamwise direction, the Reynolds number (based on the boundary layer thickness), and the crossflow strength. Turbulence intensity and surface roughness also play a decisive role, but are often not accounted for. The majority of the crossflow criteria in literature uses the helicity He = -→ u • (∇ × -→ u ) as a measure of the crossflow strength. Its use was first proposed by Muëller & Herbst, [START_REF] Müller | Modelling of crossflow-induced transition based on local variables[END_REF], then by Langtry et al., [START_REF] Langtry | Extending the γ-re θ correlation based transition model for crossflow effects[END_REF]. Both Muëller et al. and Langtry et al. CF extensions were thought to be coupled to γ -Re θ model, involving the modification of the second transport equation for Re θ t . Therefore, their approaches cannot be directly used within the one-equation γ model formulation.

Grabe et al., [START_REF] Grabe | Transport modeling for the prediction of crossflow transition[END_REF], also proposed an helicity-based crossflow criterion for γ -Re θ . The latter crossflow model, differently from the two variants mentioned above, envisages the modification of the onset function of the production term in the transport equation for the intermittency γ. This helicity-based crossflow criterion is calibrated on a 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, [START_REF] Onera | Research done at DERAT[END_REF], the 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].

A different criterion to account for stationary crossflow transition is the one proposed by Menter & Smirnov, [START_REF] Menter | Development of a RANS-based model for predicting crossflow transition[END_REF]. This is the so-called T c1 crossflow criterion and it is a local reconstruction of the C1 criterion of Daniel Arnal, [START_REF] Arnal | Théorie de l'instabilité laminaire et critères de transition en écoulement bi et tridimensionnel[END_REF]. Arnal's C1 criterion is based on the crossflow Reynolds number Re δ2 :

Re δ 2 = U 1 e δ 2 ν where δ 2 = - ∞ 0 W 1 U 1 e dy, (1)
where 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 crossflow 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 . It is defined as: (

Re δ 2 t =              
) 3
The C1 criterion is empirically-based on the 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 mea- In the specific case of the T c1 crossflow criterion, it was observed that this criterion, as originally formulated by Menter & Smirnov, did not perform well only on non-wing-like geometries. Therefore, in order to improve the T c1 performance, an original re-calibration based on FSC velocity profiles is proposed to extend its range of application to more complex 3D geometries. The details of the re-calibration are given in this paper, that includes the following sections: description of the mathematical formulation of the one-equation γ model, introduction of the FSC equations, the detailed steps of the proposed original recalibration. Finally, the numerical results obtained with the re-calibrated T c1 criterion are discussed for two different three-dimensional configurations, exhibiting strengths and weaknesses of the proposed approach.

surements
γ model with the T c1 crossflow extension was implemented in the in-house ISIS-CFD solver, whose main features are described below. [START_REF] Guilmineau | Cross wind effects on a simplified car model by a DES approach[END_REF]. The flow solver is combined with the mesh generator HEXPRESS T M , which generates full hexaedral unstructured meshes.

ISIS-CFD AT A GLANCE

The 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, Leonard, [START_REF] Leonard | A stable and accurate convective modelling procedure based on quadratic upstream interpolation[END_REF].

The following section is devoted to γ transition model formulation, as originally published by Menter et al.

One-equation γ transition model

The one transport equation for γ reads as:

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

P γ = F onset,2D [F length,2D (ρS (1 -γ)γ)], (5)
where F length,2D is constant and set to the value of 100. It represents the transition length and controls the strength of the intermittency production term. The activation of the latter is controlled by the function F onset,2D , built as:

F onset,2D-1 = Re V 2.2Re θ c , (6)
F onset,2D-2 = min(F onset,2D-1 , 2.0), (7)
F onset,2D-3 = max 1 - R T 3.5 3 , 0 , (8)
F onset,2D = max(F onset,2D-2 -F onset,2D-3 , 0), (9)
where Re V = ρS y 2 µ is the vorticity Reynolds number. Re θ c is the critical momentum thickness Reynolds number, at which the instabilities are amplified in the boundary layer. R T = ν t /ν is the eddy viscosity ratio. 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). ( 10
)
F turb deactivates the destruction term outside the boundary layer and it is given by:

F turb = e -R T 2 4 . (11)
The model constants obtained by numerical calibration are :

c e2 = 50, c a2 = 0.06, σ γ = 1. ( 12
)
The boundary conditions for γ are zero normal flux at the wall and equal to one at the inlet. The constant c e2 , Eq.( 12), controls the lower limit of γ, i.e. the value at which the destruction term E γ in Eq.( 10

The pressure gradient parameter λ θ commonly used in empirical correlations is given by:

λ θ = - θ 2 µ 1 U dP ds = θ 2 ν dU ds , (14)
where dU ds is the derivative of the velocity in the streamwise direction at the edge of the boundary layer. For the flow on a flat plate, using the incompressibility constraint, Eq.( 14) can be expressed as:

λ θ = θ 2 ν du dx = - θ 2 ν dv dy , (15)
where v and y are the wall-normal velocity component and coordinate in the free-stream, respectively. For a generic geometry, given the wall-normal velocity component v of the velocity vector -→ u , dv dy can be expressed as:

dv dy ≡ ∇( - → n • - → u ) • - → n . ( 16
)
The quantity dv dy is used in this formulation as an indicator within the boundary layer of the pressure gradient imposed from the free-stream. In order to express λ θ locally, the momentum thickness θ is replaced by the wall distance y. The transition model is built to activate close to the center of the boundary layer and at this location, given the boundary layer thickness δ, y = δ 2 ∼ θ. Thus, y gives a proper scaling of θ inside the boundary layer. Finally, the local expression of pressure gradient parameter is expressed by:

λ * θ,L = -7.57 • 10 -3 dv dy y 2 ν + 0.0128. (17)
The scaling coefficients for λ * θ,L in Eq.( 17) are obtained such that in the middle of the boundary layer λ * θ,L well approximates λ θ obtained using the Falkner-Skan profiles. The constant 0.0128 is added in Eq.( 17) in order to account for the fact that dv dy is not zero in the middle of the boundary layer for zero pressure gradient flows (β = λ θ = 0), due to the growth of the boundary layer. The local pressure gradient parameter that enters the empirical correlation has to bounded for numerical reasons, resulting in the final formulation:

λ θ,L = min(max(λ * θ,L , -1), 1). ( 18
)
The empirical correlation is finally expressed as follows:

Re θ c (T u L , λ θ,L ) = C T U1 + C T U2 * e -C T U3 T u L F PG (λ θ,L ) , (19)
where the constants are:

C T U1 = 100, C T U2 = 1000, C T U3 = 1. (20)
F PG in Eq.( 19) is built from the empirical function F * PG , calibrated on a series of Falkner-Skan profiles. The latter is defined as:

F * PG (λ θ,L ) =            min(1 + C PG1 λ θ,L , C lim PG1 ) for λ θ,L 0 min(1 + C PG2 λ θ,L + C PG3 min[λ θ,L + 0.0681, 0], C lim PG2 ) for λ θ,L < 0 (21)
The constants are:

C PG1 = 14.68, C PG2 = -7.34, C PG3 = 0.0, (22)
C lim PG1 = 1.5, C lim PG2 = 3.0. ( 23
)
Finally, F PG is given by:

F PG = max(F * PG , 0), (24)
in order to avoid negative values. γ model interacts only with the turbulence kinetic energy k transport equation of the k -ω SST turbulence model 2003 version by Menter et al., [START_REF] Menter | Ten years of industrial experience with the sst turbulence model[END_REF]. 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 , (25)
The primary production term P k is defined as:

P k = γP k (26)
while, the destruction term is given by:

D k = max(γ, 0.1) • D k . ( 27
)
D k in Eq.( 27) comes from the original k -ω SST (2003) turbulence model formulation. P k in Eq.( 26)

is computed using 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. The additional production term P lim k accounts for separation-induced transition and it is defined as follows:

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

F lim on = min max Re V 2.2 • 1100 -1, 0 , 3 , (29)
C k = 1.0, C S EP = 1.0. ( 30
)
The equation for ω remains unchanged:

∂(ρω) ∂t + ∂ ∂x j ρU j ω -(µ + σ ω µ t ) ∂ω ∂x j = γρΩ 2 -βρω 2 + 2(1 -F 1 ) ρσ ω2 ω ∂ω ∂x j ∂k ∂x j . ( 31
)
The blending function between k -ω and kformulations, F 1 , is reformulated as:

F 1 = max(F 1orig , F 3 ), (32)
F 3 = e ( Ry 120 ) 8 , R y = ρy √ k µ , (33)
where F 1orig is the original function of k -ω SST model:

F 1 = tanh min max 2 √ k 0.09dω , 500 µ ρd 2 ω , 4ρσ ω2 k CD kω d 2 4 . (34)
The transition models equations, as presented above, account for bypass and streamwise transition, either due to T-S waves or separation-induced. The present study aims to incorporate in the model formulation an empirical correlation that accounts for instabilities amplifying in the crosswise direction. The crossflow criterion is discussed upon the introduction of Falkner-Skan-Cooke equations in the following section.

Falkner-Skan-Cooke Equations

In Fig. [START_REF] Menter | Transition modelling for turbomachinery flows[END_REF], an infinite swept wing geometry is presented. Two different coordinate systems can be defined:

the wing attached one, where x is aligned with the chordwise direction, and the second reference system, with x 1 aligned with the external inviscid streamlines. The three-dimensional boundary layer equations system for the infinite swept wing, where the derivative along the span is zero ∂ ∂z = 0 , reduces to:

∂U ∂x + ∂V ∂y = 0, (35)
U ∂U ∂x + V ∂U ∂y = U e dU e dx + ν ∂ 2 U ∂y 2 , (36)
U ∂W ∂x +V ∂W ∂y = ν ∂ 2 W ∂y 2 . ( 37
)
with boundary conditions:

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

Figure 1: Infinite swept wing notations used to formulate the FSC equations.

It is supposed 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 (parallel to the leading edge)

W e is constant. The two components can be written as:

U e ∼ U ∞ x L m , (40)
W e = const, (41)
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 . (42)
Once defined the Blasius similarity variable η as:

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

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

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

β = 2m m + 1 . ( 46
)
The dash in Eq.(45) denotes the differentiation with respect to η. Finally, given W = W e g(η), Eq. [START_REF] Petzold | Transition on a wing with spanwise varying crossflow and linear stability analysis[END_REF] becomes:

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

f, f , g → 0 for η → 0, (48)
f , g → 1 for η → ∞. ( 49
)
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 , (50)
W 1 /U 1 e = (g -f ) cos(φ) sin(φ). (51)
In Eq.(50) and Eq.(51) φ 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. It is defined such that:

tan(φ) = W e U e . (52)
According to the definition in Eq.( 51), W 1 = 0 for φ = 0°and φ = 90°, but also for zero pressure gradient flows, β = 0, 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 Calibration

The C1 criterion from Arnal of Eq.( 2) can be rewritten in the form:

Re δ 2 f (H 12 )150 = 1, (53)
where f (H 12 )150 is the value of the crossflow Reynolds number Re δ 2 t at transition onset. Menter & Smirnov in [START_REF] Menter | Development of a RANS-based model for predicting crossflow transition[END_REF] propose to split the C1 criterion in three different terms. The purpose is to identify and quantify each parameter that affects the crossflow transition onset. The C1 criterion can then be rewritten as:

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

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

Ψ ∼ X = Re δ2 Re V max , (56)
G ∼ F(H 12 ) = 1 f (H 12 ) . ( 57
)
The T c1 criterion, as formulated in Eq.( 55), is mathematically equivalent to Arnal's C1 criterion.

The function Ψ approximates 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. This quantity describes the three dimensionality of the boundary layer. Indeed, Ψ = 0 for 2D flows because w = 0, as well as all the derivatives in the spanwise direction z. Ψ is defined as:

Ψ =| - → ψ | •y, ( 58
)
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 | - → ω | . (59)
The scalar quantity Ψ can be interpreted as an indicator of the crossflow strength being proportional to the local change of the flow angle.

The local approximation in Eq.( 55) is evaluated using FSC solutions and it is used to develop a CFD crossflow criterion. An original re-calibration of T c1 FSC-based is presented in the following discussion.

The FSC equations are solved in the parameter range:

0 <β ≤ 1, (60) 0°<φ < 90°. ( 61
)
β is restricted to positive values, because crossflow instabilities occur for accelerated flow, in a favorable pressure gradient. β = 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.(52). For the calibration, Ψ and all the local variables are evaluated at the location η = η max in the wall-normal direction where T c1(η max ) reaches its maximum value. Within the FSC framework the quantity Ψ is approximated considering only the derivatives with respect to the normal direction η. The vorticity components that exist in the FSC framework are:

ω x ∼ ∂W 1 /U 1 e ∂η ∂η ∂y = (g -f ) sin(φ) cos(φ) ∂η ∂y U 1 e , (62)
ω z ∼ - ∂U 1 /U 1 e ∂η ∂η ∂y = -f sin(φ) 2 + g cos(φ) 2 ∂η ∂y U 1 e . ( 63
)
The ratio R(β, φ) = aΨ X is shown in Fig. [START_REF] Menter | A one-equation local correlation-based transition model[END_REF], where the constant a is set to 0.4 in order to match the two indicators Ψ and X for β → 0. The maximum deviation of the two indicators ratio with respect to the targeted value of 1 is about 35%. Nevertheless, this deviation occurs at the corners of the domain and lies in the limits of the experimental correlation C1 from Arnal, i.e. deviation of the experimental results from the correlation Re δ 2 t . Therefore, Ψ fairly approximates the ratio X.

The principal complication arises from the introduction of the function f (H 12 ).

The ratio R(β, φ) = a Ψ/(XF(H 12 )) is shown in Fig. [START_REF] Lopes | Simulation of transition from laminar to turbulent regime in practical applications of incompressible flow[END_REF], where the new constant a = 0.684 is chosen in order to have a value of R ∼ 1 for β → 0. It can be noticed that, at the upper corners of the domain, the ratio R now significantly departs from the targeted value of 1. This discrepancy is unacceptable and needs a correction. Menter & Smirnov propose a one-parameter function G MS as a correction of the ratio R = 0.684Ψ/(XF(H 12 )). Nonetheless, their calibration was performed for sweep angles smaller than 60°, which is a very low limit for complex non-wing-like 3D geometries. ν . The wall-normal derivative of v is defined as in Eq.( 16). The new λ θ for the crossflow inclusion is defined approximately at the middle of the boundary layer, where y = δ/2 ∼ θ, as a 3rd order polynomial in the two variables dv dy y 2 ν 310 and φ. The fitted surface is obtained using a least squares method. The surface λ * CF is given by:

λ θ = θ 2 ν dU 1 e dx = β ∞ 0 U 1 U 1,e 1 - U 1 U 1,e dη 2 . (64)
λ * CF =0.0473 -0.0001338 φ -0.02524 dv dy y 2 ν +5.493e -6 φ 2 -2.148e -5 φ dv dy y 2 ν + 0.001067 dv dy y 2 ν 2 -4.031e -8 φ 3 -2.81210 -7 φ 2 dv dy y 2 ν + 1.053e -5 φ dv dy y 2 ν 2 +0.0002366 dv dy y 2 ν 3 . (65)
For numerical reasons, the λ CF used in criterion has to be bounded. It is defined as:

λ CF = min max λ * CF , 0 , 0.16 . (66)
Then, the function G, in Eq.( 55), that accounts for the streamwise pressure gradient, is constructed as the 315 surface G = G(λ CF , φ). Its expression is given by:

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.04903 10 -5 φ λ 2 CF +0.03326 λ 3 CF , (67)
where both λ CF and φ are normalized with respect to their mean value and standard deviation. 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.

Indeed, it would require the definition of the 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 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 Choi et al., [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 ]. (68)
-→ 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 φ of Eq.( 52).

Galilean invariance. The approximation proposed for φ L uses the local velocity vector and it makes γ model not Galilean invariant. In order to achieve a "weak Galilean invariant" formulation, such that the invariance with respect to Galilean transformations is preserved despite the use of the velocity vector, the local sweep angle is numerically implemented using the relative velocity vector -→ u rel defined as:

- → u rel = - → u -- → u wall , (69)
instead of the local velocity vector. This is a fair modification, when dealing with boundary-layer transition.

In ISIS-CFD solver, 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, or a surface treated with wall functions, 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, and it is performed independently of the use of a transition model. This implementation has its own limitations. Indeed, 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.

A "strong Galilean invariant" formulation, that does not use explicitly the velocity or the streamlines direction, 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, as in its "strong" sense.

Implementation in ISIS-CFD solver. The derivative of the vorticity, which enters the T c1 criterion through the indicator Ψ, is calculated as the Hessian of the local velocity vector. 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 at 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 the region of interest.

The calculation of a 2nd derivative within the boundary layer might be troubling and requires a mesh refined enough 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.

Inclusion of Tc1 Crossflow Criterion within γ formulation

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

F onset1,CF = G Ψ Re V c 150 , (70)
F onset2,CF = min[max(F onset1,CF , 0), 2], (71)
F onset3,CF = max(1 -(R T /a) 3 , 0), (72)
F onset,CF = max(F onset2,CF -F onset3,CF , 0), (73)
a = 1.5, c = 0.6. ( 74
)
Through the proportionality constant c in Eq.(70), 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. The new F onset function that substitutes F onset,2D in the transport equation for the intermittency γ, Eq.( 5), is given by:

F onset = F length,2D F onset,2D + F length,CF F onset,CF , (75)
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,2D = 100 . The constants, a and F length,CF , c were obtained from numerical calibration.

3D Simulations

The present section is devoted to the analysis of the performance of T c1 crossflow criterion for three-

Inlet Boundary Conditions

Initial disturbances, their frequency and their amplitude are translated within the RANS framework in the value of turbulence intensity T u, and eddy viscosity ν t , that are imposed at the inlet. T u quantifies the velocity fluctuations and it is defined as:

T u = 2k 3 /U, (76)
where U is the velocity norm. ν t is linked to the turbulence reference length L turb by the relation:

ν t = T uU re f L turb = k/ω. ( 77
)
Both quantities depend on the turbulence kinetic energy k and the turbulence frequency rate ω at the inlet, and on their free-decay ahead of the body. Indeed, in the free-stream, the destruction terms of the turbulence 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]. Free decay rates of turbulence quantities can be studied from the solution of the k -ω SST equations in the approaching flow field. The transport equations for k and ω for an incompressible and steady flow reduce to:

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

dk * dx * = -β * k * ω * , ( 80
)
dω * dx * = -β(ω * ) 2 . ( 81
)
The superscript * indicates the dimensionless variables, defined as

k * = k/U 2 ∞ , ω * = ωL/U ∞ , ν * t = k * /ω * ,
with U = U ∞ , and x * = x/L. The analytical solutions of Eq.( 80) and (81) are given by:

k * = k * in (1 + β(x * -x * in )ω * in ) -β * β , (82)
ω * = ω * in (1 + β(x * -x * in )ω * in ) -1 . ( 83
)
The solution for ν * t reads as:

ν * t ν = ν * t in ν 1 + β((x * -x * in )ω * in ) β * β -1 . (84)
The subscript in indicates the variable at the inlet of the computational domain. From Eq.( 82), ( 83) and

(84), it can be observed that the rate of the decay is exponential. ν * t /ν is the quantity which undergoes the decay at the slowest rate. Rewriting the specific turbulence dissipation rate ω * in as:

ω * in = k * in ν ν * t in Re, (85)
and substituting it in Eq.( 82) and (84), it is observed that the decay depends on the Reynolds number and can be controlled through the eddy viscosity ratio. By increasing the eddy viscosity ratio at the inlet, the decay of k * and ω * can be contained. In order to limit the decay of turbulence quantities ahead of the body, a possible strategy is, therefore, to impose high values of eddy viscosity ratio R T at the inlet. This method has to be used with care because, as discussed in Spalart & Rumsey, high values of eddy viscosity ratio can pollute the flow field in non-turbulent region.

Finally, to impose the initial conditions, specific values of eddy viscosity ratio, chosen by the CFD users, are specified at the inlet, in order to obtain the desired value of turbulence intensity close to the body. The latter should be recovered from the experiments, but it is not always an available information. The value of turbulence kinetic energy at the inlet k in is then obtained from Eq.(82).

6:1 Prolate Spheroid

The 6:1 prolate spheroid is a geometry 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 cases, both experimentally and numerically, because it exhibits all the complex physics associated with crossflow transition and crossflow separation. For the validation of the numerical results here presented, it is referred to the experiments performed by Kreplin in 1985 at DLR, [START_REF] Kreplin | Wall shear stress measurements on an inclined prolate spheroid in the DFVLR 3m × 3m low speed wind tunnel[END_REF], that fully describe transition at the surface, through the measurements of wall shear stress magnitude and direction.

Experimental Set Up

Experiments were performed in the 3m×3m low speed wind tunnel at DLR Gottingen, around the 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 wall for a 2D turbulent boundary layer, and, around the 6:1 prolate spheroid at incidence, the laminar region is quite extended.

Presented computations are run for the Reynolds number Re = 6.5 × 10 6 , for which transition occurs under the interaction of T-S and CF instabilities. For angles of attack higher than α = 5°, zones of pure crossflow are observed in the middle of the inclined prolate spheroid and they become wider as the inclination is increased. Unfortunately, no specific indication on the free-stream conditions was given in the experimental report, but it is mentioned that T u varies from 0.1% to 0.3% depending on the Reynolds number. Not knowing the free-stream T u 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 are chosen in order to obtain a value of T u ∼ 0.15% in the vicinity of the leading edge. This value is obtained resolving the turbulence equations as described in Sec.(6.1).

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. Thus, the flow is aligned with the 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 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. Mesh details are given in Table [START_REF] Menter | Transition modelling for turbomachinery flows[END_REF], where h i is the typical cell size. It is defined as:

h i = N cell i ∆V i N cell , (86)
where ∆V i is the volume of the i-th cell, and N cell is the total number of cells. The coarsest grid from different points of view is shown in Fig.( 6) and [START_REF] Grabe | Transport modeling for the prediction of crossflow transition[END_REF].

N cells N surface h i /h

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 γ + CF on the coarsest grid, Grid5, were run on 280 processors, for ∼ 9000 CPU hours (a) Surface mesh (black), mesh on the symmetry plane (red).

(b) Surface mesh (black), mesh on the symmetry plane (red), mesh on the plane X/L = 0 (blue). Convergence results of friction and pressure drag coefficient with respect to the grid refinement are shown in Fig. [START_REF] Schmitt | Ecoulements subsoniques et transsoniques sur une aile en flèche variable[END_REF]. The estimation of the convergence order p of the discretization method relies on Richardson Extrapolation (RE), following the procedure from Ec ¸a et al., [START_REF] Ec ¸a | Evaluation of numerical error estimation based on grid refinement studies with the method of the manufactured solutions[END_REF], [START_REF] Ec ¸a | A verification and validation exercise for the flow over a backward facing step[END_REF], [START_REF] Ec ¸a | A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies[END_REF]. Relying on 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 ). (87)
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, as defined in Eq.( 86).

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

δ RE = φ i -φ 0 = αh p i . (88)
This expansion can be manipulated to obtain an expression for the order of convergence p. It is observed that an estimated order of convergence p = 2 of the discretization error is obtained on the friction component, as theoretically expected for a finite volume method. The pressure drag presents a quasi-2nd order convergence with respect to the grid refinement. The behavior of the pressure is not worrying, but related to the noisy convergence of the non-linear residuals. It is also noticed that C d, f decreases with the grid refinement. γ transition model works first installing turbulence in the boundary layer and further destroying it throughout the simulation. The intermittency is overestimated on coarser grids, resulting in a slightly wider turbulent region. In terms of drag coefficients, the results on the three finest grids are very similar.

The same observation is valid for the skin friction C f contours and transition onset predictions. Results for α = 15°are shown on Grid1 of Table [START_REF] Menter | Transition modelling for turbomachinery flows[END_REF].

h i /h 1 0 0.5 1
T c1 MS = G MS ΨRe V 150 . ( 89
)
In Eq.( 89), G MS formulation is given by

g(λ CF-MS ) = 8.8λ 3 CF-MS -9.1λ 2 CF-MS + 3.7λ CF-MS + 1 g(λ CF-MS ) = min[max(g(λ CF-MS ), 1), 2.3] G MS (λ CF-MS ) = 0.684 g(λ CF-MS ) , (90)
where λ CF-MS is:

λ * CF-MS = -0.1111 • dv dy y 2 ν + 2.3, ( 91
)
λ CF-MS = min[max(λ * CF-MS , 0), 0.7]. ( 92
)
In the presented simulations, it was decided for the following formulation for the F onset,CF-MS in γ + T c1 -MS model:

F onset1,CF-MS = T c1 MS /c (93) F onset2,CF-MS = min[max(F onset1,CF-MS , 0), 2], (94)
F onset3,CF-MS = max(1 -(R T /a) 3 , 0), ( 95
)
F onset,CF-MS = max(F onset2,CF-MS -F onset3,CF-MS , 0), (96)
a = 1.5, c = 0.6. ( 97
)
The final F onset that enters the intermittency production term is the same as in Eq.( 9), where F onset,CF-MS substitutes F onset,CF . In conclusion,the differences between γ + CF and γ + T c1-MS variants lie in the definition of the function pressure function, G vs G MS , as well as the pressure gradient parameter approximation, λ CF vs λ CF-MS , because of the introduction of the sweep angle φ.

In general, it is clearly visible that a correlation that accounts for crossflow instabilities is fundamental.

The skin friction contours as computed with and without crossflow criterion considerably differ. A zone of pure crossflow transition is observed at the middle of the spheroid down to the windward side. sweep angle Φ = π 2 -α, with α angle of attack, as discussed by Arnal, in [START_REF]Three-dimensional boundary layer: laminar-turbulent transition, in: Computation of Three-Dimensional Boundary Layers Including Separation[END_REF]. For the specific case of α = 15°, the geometrical sweep angle at the symmetry plane would be Φ ∼ 75°and it could be expected a local sweep angle φ of this order on the windward side. 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, [START_REF] Bégou | Prévision de la transition laminaire-turbulent dans le code elsa par la méthode des paraboles[END_REF]. In general, using γ + CF the intermittency production term activates more upstream than using γ, as well as the destruction of intermittency, E γ . Thus, at each section chosen for the comparison, it results that E γ , which is proportional to γ(1 -γ), is stronger in the computations by γ + CF than the ones from γ, explaining the lower C f . Close to the 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. This incorrect flow prediction at the windward side is due to the fact that γ + CF does not account for any transition mechanisms other than T-S waves and CF instabilities separately. Possible missing mechanisms are the non-linear interaction between T-S waves and crossflow instabilities, attachment line instabilities, and traveling crossflow modes. On the 6:1 prolate spheroid, for this angle of attack, some interaction between T-S and CF is expected. Nevertheless, the sum of the two F onset,* functions, as given in Eq.( 75), is probably not enough to account for the two modes interactions.

Leading edge instabilities are related to 2D boundary layer developing at the attachment-line. They can cause the flow to become turbulent along the longitudinal direction. This is a typical case where transition is dominated by streamwise instabilities, but the streamwise criterion Re θ -based cannot predict it. As explained by Arnal, [START_REF]Three-dimensional boundary layer: laminar-turbulent transition, in: Computation of Three-Dimensional Boundary Layers Including Separation[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 propagates. For this reason, criteria which involve boundary layer parameters fail to predict transition at the attachment lines.

Traveling crossflow waves prevail over stationary modes for a turbulence level T u > 0.2% and on smooth surfaces, as shown 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]. The level of turbulence intensity reported in the experimental report by Kreplin is in a range for which traveling modes might be important and can occur around the prolate spheroid. The T c1 criterion is based on Arnal's C1 criterion, which was empirically established based on experiences at very low free-stream turbulence levels, and it is not expected to predict traveling modes.

Skin friction lines predicted by γ + CF are shown in Fig. [START_REF] Arnal | Théorie de l'instabilité laminaire et critères de transition en écoulement bi et tridimensionnel[END_REF]. The model predicts an envelope of converging wall streamlines on the top of the spheroid coming from the windward and leeward sides, which are

representative of an open-separation, according to the definition of Surana, [START_REF] Surana | Exact theory of three-dimensional flow separation. Part 1. Steady separation[END_REF]. The following sections are devoted to the analysis of the flow around the 6:1 prolate spheroid for α = 5°, 30°.

The presented computations are computed on Grid3 of Table [START_REF] Menter | Transition modelling for turbomachinery flows[END_REF]. Indeed, considering the predicted transition front, computations on this grid are sufficiently accurate. Hereafter, the results computed with the original criterion by Menter & Smirnov, γ + T c1 -MS , are not considered in the comparison. Crossflow modes contribute to the transition predictions on the windward side. In general, γ + CF predicts a transition profile fuller towards the windward side, down to φ ∼ 60°, compared to γ results, see Fig. [START_REF] Boltz | Effects of sweep angle on the boundary-layer stability characteristics of an untapered wing at low speeds[END_REF].

Nevertheless, the deviation between measurements and predictions by γ + CF on the windward side is still considerable at all the analyzed sections. This is due to the fact that the interaction between T-S and CF

Sickle Wing: CASE A

The last validation case is the sickle wing configuration. The 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. The interest of the authors toward this geometry is to assess how much these strong spanwise gradients challenge the FSC assumptions, therefore the T c1 crossflow criterion. Experimental data for different Reynolds numbers and angles of attack are available in Kruse et al., [START_REF] Kruse | Transition prediction results for sickle wing and NLF (1)-0416 test cases[END_REF]. The experimental campaign was run in an atmospheric, closed circuit type low-speed wind tunnel, the DNW-NWB in Braunschweig.

The free-stream conditions are not given in Kruse's report. Nevertheless, it is mentioned that the average free-stream disturbance level is lower than in the experiments run from Petzold et al., [START_REF] Petzold | Transition on a wing with spanwise varying crossflow and linear stability analysis[END_REF]. The latter were run around the same geometry, for the same Reynolds and incidence conditions, but in a different wind tunnel. Petzold et al. reported an average turbulence intensity of T u = 0.17% that is used for the present computations.

Experimental Set Up

The wing consists of five sections: an unswept peniche raises the model above the tunnel wall, in order α Re T u in (%) (ν t /ν) in T u(%) -2.6°2.75 × 10 6 0.20 2.24 0.17

Computational Costs

The simulations for the sickle wing were run on HPC national resources as well. The simulation on the grid, whose details were given in Table (3), was performed on 784 processors for approximately ∼ 27440 730 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. [START_REF] Leonard | A stable and accurate convective modelling procedure based on quadratic upstream interpolation[END_REF]. Predictions by γ without and with crossflow inclusion, γ and γ + CF variants, respectively, are compared to the measurements. The latter consist of the extracted experimental transition lines on both surfaces. 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. approach, its main drawback is that it does not allow to control the surface aspect ratio of the grid at specific spanwise locations. Thus, the risk is that the overall number of cells might explode, reaching values that 765 are inconceivable for a RANS simulation. As example of a proper mesh that should be used for transition simulations around the sickle wing the reader should refer to Kim et al., [START_REF] Kim | Boundary layer transition models for naval applications: capabilities and limitations[END_REF].

On the lower surface, the transition line is straight over the whole span.

Conclusion

This work was devoted to the analysis of γ +CF transition model performance for 3D boundary layers. The recalibrated T c1 criterion here presented is used for the inclusion of transition due to stationary crossflow vortices. γ + CF variant was proven to perform well on the two tested geometries, the 6:1 prolate spheroid

  The one-equation local correlation γ transition model was proposed in 2015 as a drastic simplification of γ -Re θ model. The transport equations are reduced from two to one and the critical momentum thickness Reynolds number is computed through a fully local empirical correlation.

  around the ONERA D profile by Schmitt et al. and at ONERA/CERT/DERAT. The T c1 local approximation proposed by Menter & Smirnov uses the wall-normal directional change of the normalized vorticity vector, hereafter defined Ψ, as crossflow strength indicator. This quantity represents a measure of the three-dimensionality of the boundary layer. The quantity Ψ enters the critical crosswise Reynolds-number ratio of Eq.(2), which is evaluated locally through auxiliary functions expressed by the solution of the Falkner-Skan-Cooke (FSC) equations. The FSC are an extension of two-dimensional laminar boundary layer Falkner-Skan equations by considering the spanwise velocity component for a flow over an infinitely yawed wedge, i.e. the gradients in the spanwise directions are zero, at zero angle attack, Cooke, [16]. In Menter & Smirnov paper, the Tc1 criterion is used in conjunction with γ model, modifying the production term in the transport equation for the intermittency. In this context, it is opted to pursue this research using the FSC-based T c1 criterion, rather than the helicitybased criterion by Grabe et al. , as possible crossflow variant of the one-equation γ model. The rationale behind this choice is that a FSC-based calibration widens the "room of maneuvers" for model/code developers with respect to a calibration based on experimental databases reconstructed numerically. The use of an FSC framework allows to introduce more physics in the correlations used to trigger the transition process and to better understand a model behavior in default of detailed experimental data for transitional flows.

  ) changes sign. The value c e2 = 50 enforces γ = 1/c e2 = 0.02 in the viscous sub-layer down to the wall. Re θ c is an empirical correlation and it is function of the local turbulence intensity T u L and pressure gradient parameter, λ θ,L . The turbulence intensity within the boundary layer is expressed as:

Figure 2 :

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

Figure 3 :

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

Figure 4 :

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

  For the construction of the function G, λ CF is evaluated at the position η max , where the T c1 reaches its maximum. The new ratio R = GΨ/(XF(H 12 )) is plotted in Fig.(5) versus the two dependent variables λ CF and φ. It can 320 be observed that the benefit of introducing the definition of φ is a considerable reduction of the spread of the ratio R around the value of one, with a maximum deviation of less than 10%.

Figure 5 :

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

  dimensional configurations. Predictions by the crossflow T c1 model variant, hereafter referred to as γ + CF, are compared to results from the transition model without crossflow, indicated by the title γ, and experimental results. Two different geometrical configurations are considered: 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 inlet boundary conditions are imposed with respect to a targeted value of turbulence kinetic energy T u recovered from the experimental data. The strategy to impose the boundary conditions is explained in the next section.

485

  of simulation to converge. Convergence is controlled by a gain of minimum four orders of the normalized residuals for all the turbulence and transition variables and by forces convergence. It is worthwhile to mention 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 γ formulation, which are not smooth functions. This problem has already been observed also for 2D simulations, Lopes,[START_REF] Lopes | Simulation of transition from laminar to turbulent regime in practical applications of incompressible flow[END_REF]. Simulations on finer meshes were run on an 490 higher number of processors, up to a maximum of 784. (a) Full Computational Domain. (b) Zoom on the O-block surrounding the spheroid.

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 : 6 :

 76 Figure 7: 6:1 prolate spheroid. Coarsest mesh at α = 15°, closed up on the first half of the spheroid: surface mesh, mesh on the symmetry plane, mesh in the wall normal direction on the plane X/L = 0.

  C d,p .

Figure 8 : 6 : 1

 861 Figure 8: 6:1 Prolate Spheroid: α = 15°. Convergence of pressure and viscous drag for γ + CF formulation with respect to the grid refinement. p is the estimated convergence order of the discretization error relying on Richardson Extrapolation.

  CF instabilities also contribute to transition close to the leeward symmetry plane. The predictions by γ + CF are in very good agreement with the measurements within the range 30°< φ < 120°. Discrepancies are observed close to the two symmetry planes. The proposed calibration considerably enhances the performance of the T c1 criterion, see Fig.(9c) vs Fig.(9d), predicting a considerable fuller transition front. By including the sweep angle, it is accounted for its impact on the overall transition process. Compared to a swept wing, the 6:1 prolate spheroid, is a geometry significantly more swept. If the region next to the stagnation point is excluded, the windward symmetry plane looks like the attachment line of a cylinder with a geometrical

  (a) X/L -φ plane. Measured C f distribution. (b) X/L -φ plane. Predicted C f distribution by γ. (c) X/L -φ plane. Predicted C f distribution by γ + T c1-MS. (d) X/L -φ plane. Predicted C f distribution by γ + CF.

Figure 9 : 6 : 1

 961 Figure 9: 6:1 Prolate Spheroid: α = 15°. Measured and predicted skin friction C f distributions. Numerical results are computed using γ and γ + CF models . The results obtained with γ + T c1-MS, with the T c1 version of Menter & Smirnov, are also shown.

Fig.( 10 )

 10 Fig.(10) to Fig.(12) present the girthwise skin friction distribution at different sections. Results from Menter & Smirnov Tc1 version, γ + T c1-MS, are not considered in these plots. At X/L = 0.223, it is visible that γ + CF predicts an uniform and fuller transition front with respect to the experiments, without any visible kink. For φ < 120°, γ + CF is able to predict the transition that γ completely neglects. The different estimation of C f close to the leeward side (φ > 130°) by γ + CF and γ, from X/L = 0.223 on, is related 560

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

Figure 10 : 6 : 1

 1061 Figure 10: 6:1 Prolate Spheroid: α = 15°. Experimental and numerical girthwise distribution of C f at the sections X/L = 0.223 and X/L = 0.309. Numerical results are computed with γ + CF and γ models. At each validation point, the experimental uncertainty is reported.

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

Figure 11 : 6 : 1

 1161 Figure 11: 6:1 Prolate Spheroid: α = 15°. Experimental and numerical girthwise distribution of C f at the sections X/L = 0.395 and X/L = 0.565. Numerical results are computed with γ + CF and γ models. At each validation point, the experimental uncertainty is reported.

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

Figure 12 : 6 : 1

 1261 Figure 12: 6:1 Prolate Spheroid: α = 15°. Experimental and numerical girthwise distribution of C f at the sections X/L = 0.652 and X/L = 0.936. Numerical results are computed with γ + CF and γ models. At each validation point, the experimental uncertainty is reported.

Figure 13 : 6 : 1

 1361 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.

6. 2 . 5 .

 25 Results: α = 5°T he 6:1 prolate spheroid at α = 5°undergoes a different transition process than α = 15°, nevertheless, the conclusions on the γ + CF transition model 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 most probably occurs because of the non-linear interaction between T-S and CF instabilities. The experimental C f contours in the x/L-φ plane are shown in Fig.(14a). Predictions by γ + CF are shown in Fig.(14c). The inclusion of a crossflow criterion contributes to move upward the transition front in the region included in the range 70°< φ < 160°, if compared to γ results, shown Fig.(14b). Transition on the windward side is not predicted and this feature is related to the missing transition mechanism(s) that the model does not account for.

  (a) X/L -φ plane. Measured C f distribution. (b) X/L -φ plane. C f distribution by γ . (c) X/L -φ plane: predicted C f distribution by γ + CF.

Figure 14 : 6 : 1

 1461 Figure 14: 6:1 Prolate Spheroid: α = 5°. Measured and predicted skin friction C f distributions. Numerical results are computed using γ + CF and γ models.

Fig.( 15 )

 15 Fig.(15) to Fig.(16) show the girthwise distribution of the skin friction coefficient at different X/L planes, for X/L 0.480. Only the second-half of the prolate spheroid is considered, because the flow is laminar in 620

  625

  mode is not accounted for by γ + CF variant. Predictions by γ + CF and γ become very similar towards the trailing edge, Fig.(16), close to the location where the flow is forced by the geometry to detach. Indeed, the skin friction lines as computed by γ + CF, shown in Fig.(17), do not show any characteristic topological sign of crossflow separation prior to the end of the body.(a) X/L = 0.480. (b) X/L = 0.565.

Figure 15 : 6 : 1

 1561 Figure 15: 6:1 Prolate Spheroid: α = 5°. Experimental and numerical girthwise distribution of C f at the sections X/L = 0.480 and X/L = 0.565. Numerical results are computed with γ + CF and γ models. At each validation point, the experimental uncertainty is reported.

Figure 16 : 6 : 1

 1661 Figure 16: 6:1 Prolate Spheroid: α = 5°. Experimental and numerical girthwise distribution of C f at the sections X/L = 0.738 and X/L = 0.936. Numerical results are computed with γ + CF and γ models. At each validation point, the experimental uncertainty is reported.

Figure 17 : 6 : 1

 1761 Figure 17: 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) X/L -φ plane: Measured C f distribution. (b) X/L -φ plane. Predicted C f distribution.

Figure 18 : 6 : 1

 1861 Figure 18: 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.

  (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 19 : 6 : 1

 1961 Figure 19: 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 20 : 6 : 1

 2061 Figure 20: 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.

  to avoid the contamination of the laminar boundary layer by the turbulent boundary layer at the tunnel wall. The central three swept sections, A, B, C are dedicated to transition measurements. Each segment has the same thickness and span length, but with increasing sweep 30°, 45°, 55°, respectively. The wing is completed by a 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, the flow around the wing operating near its design point, at Re = 2.75 × 10 6 and α = -2.6°, is considered. For these conditions, T-S and CF waves occur at the same time. Each section presents a peculiar transition pattern. The 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. Transition lines are extracted by infrared images using a temperature gradient detection scheme. The local transition position is recovered from the minimum on the temperature gradients, Kruse et al.,[START_REF] Kruse | Determination of the critical cross flow n-factor for the low-speed wind tunnel braunschweig (dnw-nwb)[END_REF]. On the upper surface, transition on section A is dominated by T-S waves. Nevertheless, moving upstream in the spanwise direction, the transition process becomes to be CF-dominated. The zig-zag pattern, observed in the experiments, 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. Experimental data included pressure tap measurements approximately along the midsection of each segment. The average model surface roughness of about 1.47 µm is neglected in the numerical simulations.6.3.2. Computational DomainThe computational domain presented in Fig.(21) reproduces the wind tunnel dimension from the experiments of 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 agreementwith the experimental results, the tunnel walls are treated as slip boundary conditions. No slip conditions are imposed on the body. The wing is turned with respect to the angle of incidence and the flow is aligned with x-axis. The mesh was generated using the hex-based unstructured grid generator HEXPRESS, and further refined using the systematic grid refinement (SGR) functionality implemented in the in-house flow solver ISIS-CFD. This SGR feature is a simplified variant of the adaptive grid refinement functionality: the grid is systematically refined without being controlled by any specific flow-feature. This approach allows to significantly improve the quality of the grid in the boundary layer and it is very useful when low-Re near wall grids need to be generated. Starting from an initial mesh of 10M of cells, with approximately 50 cells in the boundary layer, the mesh has been further refined with SGR in the longitudinal x-direction. 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 those 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. Different grid refinement levels have been tested, nevertheless the one described above is relative to the grid used for the presented computations hereafter. The overall mesh is shown in Fig.(21a), as well as the mesh at the surface, Fig.(21b). In the latter figure, the solid red lines represent the boundaries between the three central swept sections, at which 725 the geometrical sweep angle changes. Details on the mesh are shown in

  (a) Total computational domain.(b) Mesh at the surface.

Figure 21 :

 21 Figure 21: Sickle Wing. Grid used for the flow simulation around the sickle wing. Total computational domain and representation of the mesh at the surface. In Fig.(21b), the solid red lines represent the kink region, where the geometrical sweep angle changes.

  The differences between predictions by γ and γ + CF models are noticeable, especially on the upper surface, where transition is dominated by crossflow instabilities, Fig.(22a) vs Fig.(22c).γ predictions are in good agreement with measurements on Sec. A. Indeed, 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, because of the amplification of stationary crossflow instabilities. 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. The numerical simulation does not account for these roughness effects and γ + CF predicts a more uniform transition front. 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 model. The turbulent wedges in the regions of sweep changeover are reasonably well predicted by γ + CF, despite the T c1 calibration on FSC solutions, for which zerospanwise gradients are assumed. They are nonetheless polluted by the discretization error committed at the junctions between the segments, where the spanwise gradients are considerably strong. In fact, the turbulent wedges predicted in the computations by γ + CF model are diffused on their sides, and they are not as sharp as in the measurements. This discretization error is due to the fact that the grid at the surface is generally too coarse in the spanwise direction. Notably in correspondence of the junctions, i.e. the red lines at the surface in Fig.(21b), the mesh is not refined enough to predict accurately the strong spanwise gradients, that cause the turbulent wedges.

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

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

  Transition is dominated by T-S waves. The crossflow criterion should contribute to the small kinks in the wing planform. Nevertheless, γ + CF reproduces these features to a lesser extent, and the transition line is predicted upstream of the experimental front, Fig.(22d). This behavior is related to the discretization error that pollutes the overall results and an excessive diffusion of the intermittency γ. The measured and calculated pressure coefficient distributions are presented in Fig.(23), as extracted along each midsection. The numerical results are computed with γ and γ + CF transition model. Over all, a very good agreement is achieved between numerical and experimental results at each analyzed section. The numerical C p is underestimated by both γ variant at Sec.B, Fig.(23b), on the lower surface, i.e. the minimum computed C p is lower than the measured one.Nevertheless, such quantitative differences are not relevant to transition predictions, that are affected by the pressure gradients. At all sections on the upper surface, the favorable pressure gradient region is considerably extended. On Sec.A, Fig.(23a), the pressure distribution indicates the presence of the small laminar separation bubble mentioned above. The laminar separation point is approximately located at x/C ∼ 0.75. This flow feature is only predicted by γ model, without crossflow inclusion. The latter also predicts laminar separation bubbles approximately at the same position for all the other sections that are not observed in the experimental results. On the lower side, where transition is dominated by adverse pressure gradients, only γ, without CF, is able to reproduce the characteristic wiggles between x/C ∼ 0.35 and x/C ∼ 0.40. These flow features are not predicted by γ + CF, nevertheless their absence does not affect the overall transition process, if evaluated with respect to the C f contours.

Figure 23 :

 23 Figure 23: 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. Numerical results are computed by γ and γ + CF models.

790

  and the sickle wing, at different angles of attack. Within the new 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, as published by Menter & Smirnov. γ + CF performs very well around the 6:1 prolate spheroid at different angles of incidence. The unpredicted flow transition on the windward side is not crossflow related, but is due to missing physics within the mathematical model formulation. Possible physical mechanisms that γ + CF does not account for are the non-linear interactions between crossflow and T-S waves, attachment line instabilities, and traveling crossflow. In spite of the T c1 calibration with respect to the solutions of the Falkner-Skan-Cooke equations, which assume zero-spanwise gradients, the results for the flow around the sickle wing are very encouraging. γ + CF model is able to predict the turbulent wedges at the junctions of the sickle shaped planform, due to the local amplification of stationary crossflow modes. However, the discretization error resulting from the coarseness of the mesh at the junctions between the three central swept section pollutes the result. Strategies to create a proper mesh within Hexpress are currently under study. Future works will be devoted to the inclusion of the non-linear interaction effects of crossflow and Tollmien-Schlichting modes within γ + CF formulation. It is the authors belief that including additional correlations accounting for other transition mechanisms is crucial to further validate the modeling performance of RANS γ model.

  

	300 π arctan	0.106 (H 12 -2.3) 2.052 , for 2.3 < H 12 2.7,
	150, for H 12 2.3.

Table 1 :

 1

				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
	Geometry	α	Re	T u in (%) (ν t /ν) in T u(%)
		5°6
	:1 Prolate Spheroid 15°6.5 × 10 6	0.5	250	0.15
		30°T
		able 2: 6:1 Prolate Spheroid: Computations details.

6:1 Prolate Spheroid: Mesh details.

The inlet conditions for the angles of attack under study are reported in Table

(

2): the turbulence intensity T u in and the eddy viscosity ratio (ν t /ν) in at the inlet, and the value of turbulence intensity T u predicted in the vicinity of the body .

  Table(3). The inlet conditions are reported in Table(4).

	N cells N surface	y + max
	54M 462243 0.436709

Table 3 :

 3 Sickle Wing. Mesh details.

Table 4 :

 4 Sickle Wing. Inlet conditions for the computations.

This is a fair approximation because the skin friction lines are converging from different sides of the prolate spheroid.

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.

The help of Prof.Dr.Ing. Rolf Radespiel from TU Braunschweig, who provided us with these experimental data from the Sickle Wing, is kindly acknowledged. Finally, Prof. Christophe Airiau from Paul Sabatier University in Toulouse is sincerely acknowledged for sharing his Falkner-Skan-Cooke solver with the authors.