Transition Modeling for Low Pressure Turbines Using Computational Fluid Dynamics Driven Machine Learning

: Existing Reynolds Averaged Navier–Stokes-based transition models do not accurately predict separation induced transition for low pressure turbines. Therefore, in this paper, a novel framework based on computational ﬂuids dynamics (CFD) driven machine learning coupled with multi-expression and multi-objective optimization is explored to develop models which can improve the transition prediction for the T106A low pressure turbine at an isentropic exit Reynolds number of Re 2 is = 100,000. Model formulations are proposed for the transfer and laminar eddy viscosity terms of the laminar kinetic energy transition model using seven non-dimensional pi groups. The multi-objective optimization approach makes use of cost functions based on the suction-side wall-shear stress and the pressure coefﬁcient. A family of solutions is thus developed, whose performance is assessed using Pareto analysis and in terms of physical characteristics of separated-ﬂow transition. Two models are found which bring the wall-shear stress proﬁle in the separated region at least two times closer to the reference high-ﬁdelity data than the baseline transition model. As these models are able to accurately predict the ﬂow coming off the blade trailing edge, they are also able to signiﬁcantly enhance the wake-mixing prediction over the baseline model. This is the ﬁrst known study which makes use of ‘CFD-driven’ machine learning to enhance the transition prediction for a non-canonical ﬂow.


Introduction
As air travel is predicted to see continued growth, the aviation industry is focused on further improving the efficiency of gas turbines.Any technological advances which incrementally increase gas turbine efficiency can save millions of dollars, reduce carbon emissions, and create a sustainable future.One place to bring about such a change is at the grassroots level of engine component design by understanding and enhancing the ability to model the flow physics of these components in an increasingly holistic, yet cheaper manner.
Low-pressure turbines (LPT) are one of the most important components which drive the efficiency of gas turbines.If the efficiency of the LPT increases by 1%, the specific fuel consumption will decrease by 0.6-0.8%[1], translating into considerable fuel savings and emissions reductions, considering the volume of operations.LPT design and analysis is challenging as there is a large variation in flow physics between the aircraft taking off and entering cruise conditions.The Reynolds number of LPTs ranges from 0.5 × 10 5 at cruise to 5 × 10 5 at takeoff [2].This causes a change in the boundary layers from turbulent at takeoff to transitional at cruise, which is the major portion of the flight.At lower Reynolds numbers (cruise), boundary-layer transition and separation play an important role in determining engine performance.There are three modes of transition: natural, bypass and separated-flow [3].In LPTs, bypass and separated-flow transitions are the most common transition modes.Natural transition begins with weak instabilities in the laminar boundary layer, which lead to three-dimensional instabilities resulting in turbulent spots [4].However, when free-stream turbulence intensities are greater than 1%, boundary layers typically proceed from laminar to fully turbulent states without the occurrence of linear instabilities, thus bypassing various stages of natural transition.This is called bypass transition.Transition may also occur in the free shear-layer flow that results when a laminar boundary layer separates.The flow could reattach and form a laminar-separation/turbulent-reattachment bubble or could remain detached resulting in an open separation.Separated-flow transition occurs in regions of adverse pressure gradients near the leading edge or after peak suction in LPTs.Separated-flow transition also plays an important role in periodic unsteady transition wherein unsteady incoming wakes from previous blade rows interact with the separation bubble on the LPT suction-side [5].Separated boundary layers with longer separation bubbles produce greater loss and cause larger deviation in exit flow angles [6], and they are thus to be avoided.At low Reynolds numbers, transition controls the extent of the separation bubble and a reliable prediction of these aspects is crucial for the design of blades with high-lift.
With the increase of computational power, high-fidelity simulations [1,7] that are able to accurately predict the transition phenomena have become feasible.However, they are still not viable as an iterative industrial design tool, due to the associated high cost.Thus, industry remains quite heavily dependent on RANS (Reynolds Averaged Navier-Stokes) calculations as they are a cost-effective way of modeling turbulent flows for engineering solutions.Several transition models are at the disposal of turbomachinery designers [8].The γ − Re θ model [9] uses two transport equations for the intermittency (γ) and transition momentum thickness Reynolds number (Re θ ), in addition to the turbulent kinetic energy (TKE) and specific dissipation.A class of models uses the laminar kinetic energy (LKE) concept [10,11], wherein an additional transport equation is solved for the LKE.The kv 2 ω model solves an additional transport equation for v 2 , which represents wall-normal velocity fluctuations responsible for transition initiation.However, all of these existing RANS-based transition models struggle to accurately predict the onset and extent of separation on the suction-side of the aft portion of the LPT blade.Therefore, novel methods and models have to be explored in order to enhance the accuracy of transition prediction in the context of RANS.
Data-driven turbulence modeling is a field that has gained tremendous momentum in recent years [12].In fact, many leading researchers [13][14][15] view data-driven turbulence modeling as a promising avenue for bridging the gap between the accuracy of high and low fidelity simulations, while maintaining the same order of cost as low-fidelity calculations.With the help of machine learning processes, high-fidelity datasets can be harvested to extract meaningful physics-based insights to improve the accuracy of RANS calculations.Many different data-driven approaches have been applied to a large variety of industrially relevant cases.Neural networks embedding Galilean invariance to modify the anisotropy tensor were proposed by Ling et al. [16].With these neural networks, they improved the Reynolds stress and velocity field predictions over linear eddy viscosity models for duct flows and flow over wavy walls.Milani et al. [17] extended the work of Ling et al. [16] and developed a deep neural network to model the scalar-flux of turbulent diffusivity for two different jets in crossflow.Milani et al. [18,19] also used random forests to improve heatflux modeling for film cooling applications.A more comprehensive overview of existing data-driven methods for turbulence model development can be found in [12,20,21].
It is evident that neural networks are the most popular way to improve the Reynolds stress and mean flow prediction.However, one major drawback lies in their 'black box' nature-they are difficult to interpret as physics-based information remains hidden in the network.An alternative that addresses this issue is symbolic regression, which produces interpretable and explicit model forms.Ferreria [22] introduced gene expression programming (GEP) which was a new evolutionary algorithm for solving regression based problems.Weatheritt and Sandberg [23] developed GEP for turbulence modeling, which creates Galilean invariant turbulence closures, also known as explicit algebraic Reynolds stress models (EARSM).The GEP algorithm for turbulence model development produces tangible symbolic expressions which can be correlated to physical phenomena and can be directly inserted into RANS codes.For example, Akolekar et al. [24] correlated different models from GEP to physical mechanisms such as laminar separation and wake-wake interactions in LPTs with unsteady inflow conditions using a zonal model development approach.
GEP also offers the added advantage of less computational cost over neural networks, in both the development and implementation phases.GEP has been successfully applied to a number of canonical and industrial (complex flows) turbulence modeling problems by enhancing the flow prediction capability for RANS and LES.Akolekar et al. [25] developed EARSMs for enhancing the prediction of LPT wake flow originating from open and closed separated boundary layers using 'frozen' [26] high-fidelity databases (the 'frozen' methodology only solves the RANS equations for those quantities which cannot be directly obtained from a high-fidelity database-for example, the specific dissipation rate (ω) is not readily available from a high-fidelity database and is thus derived from a RANS calculation while keeping all other values constant-i.e., as derived from high-fidelity calculations).They also tested these EARSMs successfully on high-pressure turbine (HPT) flows [27], thus demonstrating the robustness of the models developed.Weatheritt et al. [28] developed GEP-based EARSMs for HPTs and assessed their performance in an a priori sense.Recently, Zhao et al. [29] introduced a novel 'CFD-driven' (computational fluid dynamics) framework for turbulence model development based on GEP.This algorithm uses RANS feedback as part of the model training process.It was found that this method has high potential for offering robust turbulence closures to enhance the wake mixing prediction across different high and low pressure turbine operating conditions and blade geometries [30].Waschkowski et al. [31] recently proposed a multi-objective and multi-expression extension to the 'CFD-driven' framework.The multi-objective approach makes use of different cost functions simultaneously, whereas the multi-expression capability develops coupled turbulence closures across different regions of the flow field.
Efforts to enhance the accuracy of transition modeling with machine learning have only recently begun.Zafar et al. [32] developed a convolutional neural network based transition model capable of accurately predicting the transition onset location for incompressible, two-dimensional attached flows in a physically informed manner, without requiring direct computations, from linear instability theory.Yang and Xiao [33] developed a framework to improve transition models using field inversion and machine learning for airfoils.There are no studies, however, which focus on enhancing transition prediction of RANS models using machine learning for LPTs.There is thus immense scope to improve the performance of existing RANS models for transition prediction with machine learning techniques.GEP is an ideal algorithm to formulate new transition models as expressions can be formulated which yield physical insights into the modeling process.
This study, therefore, proposes a novel machine learning framework based on GEP to generate models which enhance the prediction of separated-flow transition in LPTs.First, a brief overview of the LKE transition model [11], which is the baseline model used in this study, is presented.The functional forms of the LKE model's production and transfer terms are then modified using non-dimensional pi groups [34] based on local flow quantities.An extensive discussion ensues on the scientific ideology of making use of these particular pi groups.Using the recently developed multi-objective cost function based GEP algorithm [31] and the 'CFD-driven' [29] framework, models which enhance the separated-flow transition are proposed.The performance of these models is evaluated in terms of Pareto analysis and physical characteristics such as wall-shear stress, blade loading and turbulent wake mixing.This is the first known study which makes use of 'CFD-driven' machine learning to enhance the separation-induced transition prediction based on multi-objective and multi-expression optimization using non-dimensional pi groups.

Transition Modeling
The LKE transition model [11] is used as the baseline model.This section discusses the transport equations used for the LKE transition model.It also discusses the intricacies of making possible improvements in the LKE model using non-dimensional pi groups.

Laminar Kinetic Energy Model
This study makes use of a transition model developed by Pacciani et al. [11] which is based on the LKE concept, as introduced by Mayle and Schulz [35], and takes into account the pre-transitional rise of the fluctuating kinetic energy.The modeled transport equation for the LKE (k l ), with ν as the kinematic viscosity, is The gradual rise of pre-transitional energy in the separated region is modeled using a production term (P l ) based on the mean shear rate (S) and a laminar eddy viscosity concept [36].It is expressed as where ν l is the laminar eddy viscosity, C 1 = 0.01, f 1 is defined as a function of the inlet free-stream turbulence intensity (Tu), δ Ω is the maximum in the normal-to-the-wall direction shear-layer vorticity thickness which is assumed as a length scale representative of separated laminar flow and l is a proper length scale, characteristic of the transition mode to be represented (natural, bypass or separated-flow).Ω is the vorticity magnitude and U is the magnitude of velocity.The choice of a vorticity-based length scale for the LKE production term was motivated by the fact that the most unstable wavelength in the separated shear-layer is related to the shear-layer thickness.The fluctuating energy is then gradually transferred to the turbulence field via a transfer term (R) which is defined as with C 2 = 0.3, C 3 = 8, C 4 = 10 and β * = 0.09.f 2 is a damping function and R y is the wall distance Reynolds number with y as the wall-normal wall distance.The quantity ψ can be considered as the transition parameter as it controls the transfer of energy from the laminar to the turbulent state, which occurs when the wall distance Reynolds number R y reaches the threshold value of C 4 .The term R appears in both laminar and TKE (k) transport equations which indicates no net change in the total fluctuating kinetic energy.This term implies that there is a transfer of energy from k l to k.The TKE transport equation is written as where σ k = 0.5 and the TKE production (P k ) is modeled as where τ ij is the Reynolds stress tensor from the Boussinesq approximation [37], S is the magnitude of the strain-rate tensor and ν T is the turbulence eddy viscosity.The closure damping functions are where Re T is the turbulence Reynolds number and is equal to k/(νω).The transport equation of the specific dissipation rate is given as with α = 3/40, β = 5/9 and σ ω = 0.5.

Proposed Model Formulation
In the LKE transition model, the LKE production (P l ) and transfer (R) terms play a crucial role in determining the transition onset and extent.It is proposed in this study to modify certain parts of these terms.The scientific ideology as well as the method of achieving this is discussed below.

Scientific Ideology
The LKE model, with some success, was applied to the prediction of separated-flow transition in LPT cascades by Pacciani et al. [11,38], where the vorticity thickness of the separated shear-layer was assumed as a length scale.However, the eddy viscosity equation (Equation ( 2)) has two significant drawbacks.Firstly, it involves non-local quantities (e.g., the shear-layer vorticity thickness).Secondly, it does not account for the free-stream turbulence level.This latter deficiency was addressed by introducing the function f 1 , which is expressed in terms of the inlet free-stream turbulence intensity Tu.Such a modification yielded satisfactory results on isolated cascades, as shown by Pacciani et al. [39].However, the ν l term has to consequently be modified in order to include only local quantities.
The transfer term R is formulated following Walters and Cokljat [10].Equation (3) has proven to be reliable in RANS calculations of separated-flow transition as it results in separation bubble lengths, transition onset locations and pressure recovery prior to reattachment that are in good agreement with experimental data [11,38].On many accounts, it is observed that the transition onset is dependent on the Reynolds number, pressure gradient and free-stream turbulence intensity.It is therefore argued that the term R y − C 4 in the ψ function is too crude for triggering the transition onset.A better formulation, which is representative of the flow physics, has to be found for ψ.
In addition, the LKE model [11] was conceived for separated-flow transition prediction as the adopted length scale is a useful measure of the most unstable Kelvin-Helmholtz (KH) modes.The length scale also bears a loose relation with fluctuations developing in attached boundary-layers (Tollmien-Schlichting waves and Klebanoff streaks).A step forward towards a true LKE-based multi-mode model requires a more general framework that is able to address the growth of pre-transitional fluctuations in attached and separated boundary layers under the effects of the turbulent scales coming from the free-stream.Therefore, we propose the use of non-dimensional pi groups in order to modify certain terms of the LKE model.The use of machine learning (see Section 3) will make the process more efficient and holistic.

Non-Dimensional Pi Groups
The transition problem is based on a number of local independent dimensional variables.In order to best capture the effect of these independent dimensional variables on the dependent variables (for example, the ψ term), the Buckingham pi theorem [34] is used to create dimensionless groups, called pi groups.Seven non-dimensional pi groups are considered as possible input parameters to modify the LKE production term (Equation ( 2)) and the transfer term (Equation ( 3)) as where ν is the kinematic viscosity and the turbulence length scale is defined as The particular choice for the pi groups is motivated by their widespread use as transition sensor functions in a variety of models (Π 4 , Π 5 and Π 6 [8]), by their nature of scale-determining terms for laminar fluctuations (Π 1 , Π 2 and Π 6 ) or for turbulent fluctuations coming from the free-stream (Π 3 and Π 7 ).Π 1 can be regarded as the ratio of two time scales: the time scale for the convection of disturbances in the laminar boundary layer (1/omega) and the time scale for the diffusion of large-scale fluctuations (e.g., KH instabilities) from the edge of the shear-layer (k l /νω 2 ).Π 2 and Π 6 are relevant scales for the separated shear-layer thickness with the advantage of being local parameters.In the construction of such groups, the use of the vorticity magnitude in favor of the strain rate appeared somewhat arbitrary.It was therefore decided to include both formulations, letting the machine learning framework determine the respective roles.Π 3 is a scale-determining term for the turbulent fluctuations coming from the edge of the boundary layer.Π 4 is a turbulent Reynolds number based on wall distance.It accounts for the rise of turbulent kinetic energy inside the boundary layer.Its use as a transition onset parameter is extremely common in transition frameworks coupled with kω turbulence models [8][9][10].In this role, it is a fundamental component of the original model formulation and it has proven effective in switching on the transfer term when transition is initiated.Π 5 is a well-established shear-sheltering parameter.Shear-sheltering is the phenomenon by which small-scale, high-frequency, freestream perturbations are damped out by the laminar boundary layer.Π 7 is again the ratio of two time scales, where the laminar fluctuation time scale is now compared with the turbulent time scale (1/ω).
Based on the shortcomings of the LKE model (Section 2.2.1), we propose to modify the laminar eddy viscosity (ν l ) term from the LKE production term (P l ) and the transition parameter (ψ) from the transfer term (R) using the non-dimensional pi groups.If a local formulation for ν l has to be found in terms of non-dimensional pi groups and local flow quantities, it can be represented as where f 1a = f (Π i ).The machine learning processes in the model development process are tasked to find an expression for f 1a based on symbolic regression.There is no restriction on the order or the number of terms in the expression of f 1a .In order to ensure realizability of the production term, only the terms which ensure f 1 > 0 are selected for model evaluation, which is enforced by a limiter.It was also found that, to keep an upper limit on the value of f 1a , an additional limiter min( f 1 , f 1a ) is also used.This limiter ensures stability and ease of convergence in the calculations.It also ensures that f 1a is activated only in the regions of interest, and not everywhere, as is evident in Section 5.2.3.
In its current form, ψ can be represented in terms of the non-dimensional pi groups as ψ = max(Π 4 − C 4 , 0).Initially, an attempt was made to find an expression for ψ without the use of a limiter.However, it was found that it was creating non-physical and unstable solutions in the laminar boundary layer as it was activating the transfer term which should be zero in the laminar boundary layer.Another consideration was to model the entire transfer term (R); however, the machine learning algorithm was not proposing a transfer function which would tend to zero when moving closer towards the wall.A consideration was also made on whether only the LKE production term or transfer term can be modified, without modifying the other one.However, since there is strong coupling between the two terms, both have to be modified together.Therefore, using any linear or non-linear combination of the non-dimensional pi groups, an expression is to be found for

Machine Learning Framework
This section discusses the various aspects of the machine learning transition modeling framework.The advantages of the 'CFD-driven' algorithm and the extensions by Waschkowski et al. [31] are explained.The methodology of this framework is also elucidated.

Advantages
GEP-based EARSM development had been centered around the 'frozen' approach [24,25] until the integration of CFD with GEP was made by Zhao et al. [29].The 'CFD-driven' approach is feasible with GEP, because the GEP algorithm returns explicit model equations that can be instantly and automatically be implemented into RANS solvers during the model development process.Since RANS calculations are then performed as part of the training loop, the resulting 'CFD-driven' models are tested in their target environment and are thus ready to be directly implemented into industrial design tools.Another advantage of the 'CFD-driven' training is that the definition of the cost function is more flexible as compared to the 'frozen' approach.Rather than being restricted to the quantities that can be directly calculated from the model equations, such as the anisotropy tensor or TKE production, the cost functions used in the 'CFD-driven' approach can be tailored to capture any important flow feature of interest.One can specify any mean or secondary statistic profiles in order to derive a suitable model.In fact, if large databases of high-fidelity data are unavailable and only few experimental or high-fidelity profiles are available, this approach has the potential to offer suitable models based on this limited data.This approach also ensures that the models generated are numerically stable as they have already been executed in CFD calculations.This 'CFD-driven' framework has so far been successfully applied to turbulence model improvement for wake mixing problems [29,30].
All current machine learning frameworks in the context of data-driven turbulence modeling, including GEP-based model development (both 'frozen' and 'CFD-driven'), are limited to developing one model per training run and only use one cost function (or a single training objective).However, the recent study of Waschkowski et al. [31] introduced multi-expression and multi-objective capabilities as extensions to the 'CFD-driven' approach to overcome these limitations.By developing expressions for multiple models, e.g., an anisotropy and a turbulence heat flux model, simultaneously and assigning a shared cost function value to the trained models, multi-expression training is able to take the coupling effects of these models on the RANS predictions into account.This extension derived significantly improved turbulence models for a vertical natural convection flow with strongly coupling momentum and thermal fields.The second extension by Waschkowski et al. [31] made the use of multiple cost functions feasible for the 'CFDdriven' framework.A multi-objective optimization algorithm based on the concept of Pareto dominance was implemented.Previously, with only one cost function available, several quantities of interest were combined into a composite function, which required weighting these quantities before training.The multi-objective extensions allows a trade-off between the cost functions after analyzing the training results.Additionally, this extension yields a diverse set of candidate solutions, which was demonstrated to be an important asset in order to select physically plausible models.Initially, a population of candidate models is created.In contrast to previously developing EARSMs, the GEP algorithm formulates two expressions per candidate solution-one each for f 1a and f 2a .The 'CFD-driven' approach reads the candidates' functional forms into RANS solvers and performs a series of RANS-based CFD calculations in parallel.The converged CFD solutions are then evaluated against one or more cost functions.This work makes use of two cost functions-the normalized mean squared errors between RANS and time-averaged high-fidelity data of the wall-shear stress and the pressure coefficient.Section 5.2.1 discusses these cost functions in greater detail.Finally, based on the cost function values or fitness of the candidate models, a new generation of candidates is formed through natural selection and genetic modifications (refer to [23,40] for additional details).The updated models are evaluated again through integrated CFD calculations.This training process iterates and the population of candidate models evolves over numerous generations.

Numerical Setup and Configuration
The numerical setup and the flow configuration are briefly discussed in this section.The T106A LPT linear cascade was used to develop transition models with the 'CFD-driven' multi-objective framework.Steady, compressible and 2D RANS calculations with steady inflow conditions were conducted using the TRAF solver [41].The focus of this study was the midspan section of the T106A LPT blade, which is the highest source of loss [6] in an LPT cascade.The midspan section can be approximated as a 2D (quasi-3D) case.
The flow configuration was chosen so that the inflow conditions, Reynolds number and Mach number match the experimentally-validated direct numerical simulation (DNS) data of Michelassi et al. [1].The isentropic exit Reynolds number and Mach numbers were Re 2is = 100,000 and Ma 2is = 0.4, respectively.The inlet free-stream turbulence intensity (Tu) was set to 4%.
Figure 2 shows the LPT domain used with the real (C) and axial chords (C ax ) shown.At the inlet, the stagnation pressure and temperature along with the direction of the velocity vector (46.1 • anti-clockwise from the horizontal at the blade leading edge) are specified.The upper and bottom edges of the domain have pitchwise periodic boundary conditions.The grid, which has been used for previous studies [30,40], has approximately 83,000 points in the domain, with an O block around the blade and an H-block stitched on for the wake region.There are more than 30 points in the boundary layer with a y + < 1 for the first-offthe-wall points on both sides of the blade.An extensive grid convergence study for this grid was conducted by Akolekar [40].All the viscous terms for the RANS calculations are discretized using second-order central difference schemes, whereas the inviscid fluxes use a second-order cell-centered scheme for discretization.The axial distance between the inlet boundary and the blade LE was established to be consistent with the position of the turbulence generating grid in the test rig where the measurements used for validation in [38,39] were taken.Such a configuration ensured a good mesh quality and was not altered in the present work.The location of the outlet boundary of the computational domain was chosen to allow wake traversing for a range of axial distances where a possible downstream blade row is expected to sit.The presence of the H-type stitched block ensures mesh density of sufficient quality to prevent the wake from being smeared.The wake profile comparisons, which give indications on the predicted level of wake mixing, are thus expected to be significantly reliable.

Results and Discussions
Using single and multi-objective optimization, models are developed using the machine learning framework outlined in Section 3 for the T106A LPT linear cascade at Re 2is = 100,000.This section discusses the performance of the models developed in terms of mean quantities such as the wall-shear stress, pressure coefficient and stagnation pressure loss profiles.Pareto analysis is also conducted to discuss the performance of the models developed via multi-objective optimization.

Single-Objective Optimization
Initially, a single-objective cost function was used to develop expressions for f 1a and f 2a .The baseline RANS case, cost function used and the performance of the developed model are discussed in this section.

Baseline RANS Calculation
Figure 3 shows the suction-side wall-shear stress (τ w ) from the baseline RANS case (i.e., using the LKE model with the relations outlined in Section 2.1) and DNS data [1].At Re 2is = 100,000, the boundary layer separates and then reattaches, thus depicting a closed separation.In Figure 3, it is also very clear that the current baseline model does not predict the wall-shear stress accurately.What matters the most from the perspective of the boundary layer loss and subsequent downstream wake mixing loss is the accuracy of the prediction of the wall-shear stress close to and downstream of the separation onset (x/C ax ≈ 0.84).As shown in Figure 3, the baseline LKE model offers a suitable prediction of the separation onset and extent, although the predicted reattachment occurs slightly earlier than for the high-fidelity data.The transition location, which occurs at the minimum wall-shear stress (x/C ax = 0.972), is predicted quite closely with respect to the DNS data.However, there is still plenty of scope to improve the prediction of the wall-shear stress in the separated flow regions, especially in the region leading to and immediately after the transition location (i.e., overshooting in the wall-shear stress profile).

Cost Function
Initially, a single-objective cost function was used in the 'CFD-driven' framework.This cost function revolved around minimizing the difference between the wall-shear stress obtained from RANS with respect to the one from DNS.Since the main objective of this study was to improve the separation prediction, the wall-shear stress data from 0.75 ≤ x/C ax ≤ 1 was selected.The cost function (J SO ) using a single-objective function (SO-Model) is

Analysis of the SO-Model
The 'CFD-driven' framework was able to provide suitable functional forms for f 1a and f 2a that significantly reduced the value of J SO .For the baseline case, the value of the cost function is 4.147 × 10 −6 , and, for the SO-model, the value of J SO is 0.192 × 10 −6 .In order to obtain this solution, 90 generations, similar to the 'CFD-driven' model development for LPT wake mixing [30], were required with a population of 125 candidates.The new formulations of f 1a and f 2a were applied across the entire blade, even though they were developed using reference data from 0.75 ≤ x/C ax ≤ 1.The wall-shear stress plot is shown in Figure 4a and the pressure coefficient plot (C p ) is shown in Figure 4b.The framework produces a model that results in a close match with the DNS wall-shear stress in the separated-flow region.It even gets the reattachment location predicted correctly.Since the region immediately upstream of x/C ax = 0.75 was not included in the cost function, there is not much change in the wall-shear stress behavior there.The SO-model, however, also improves the prediction of the leading-edge separation, even though that region was not included as part of the training data.Even though the wall-shear stress prediction has improved with the SO-model, the pressure coefficient (C p ) distribution has worsened over the baseline case, especially in the separated flow region on the suction-side.As the baseline case and DNS both show, there should be a plateau in the pressure coefficient plot when the flow separates.The SO-model, unlike the baseline case, does not predict a strong peak in the separated flow region in the wall-shear stress plot.The strong peak in the wall-shear stress is a peculiarity of the LKE framework, as it is the result of the turbulence 'ignition' by the transfer term.In the authors' experience, the position of the minimum wall-shear stress and the subsequent ramp in wall-shear stress is crucial to get the correct pressure coefficient distribution in the separated region.As the SO-model does not show this behavior, it tends to worsen the pressure coefficient predictions over the baseline case.

Multi-Objective Optimization
Next, the multi-objective optimization approach is used to develop transition models.This section discusses the cost functions and the performance of the solutions obtained with respect to the cost functions and physical characteristics.

Cost Functions
The main motivation to explore the avenue of multiple cost functions was due to the worsening in the pressure coefficient prediction in the separation region on the suctionside over the baseline case when a single objective function was used.In order to at least maintain the pressure coefficient distribution as obtained from the baseline case, an additional cost function was created based on the pressure coefficient from close to peak suction right to the blade trailing edge (TE)(0.6 ≤ x/C ax ≤ 1-this region includes the separation plateau).The two cost functions used are The multi-objective optimization was run for 90 generations with a population of 125 individuals.In order to obtain the final solution (at n = 90), the multi-objective approach used about 1500 core hours, which was about three times that of the singleobjective approach.This is because the multi-objective approach evaluates many more possible candidate solutions as opposed to a single best solution in the single-objective approach.If the data from baseline RANS calculations are fed into the cost functions, then J MO1 = 5.24 × 10 −6 and J MO2 = 3.816 × 10 −4 .

Pareto Analysis
The best way to analyze the evolution of the models generated with the multi-objective optimization approach is with Pareto analysis.In order to compare two candidate solutions, this optimization approach uses the concept of Pareto domination.One candidate is said to dominate the other candidate when its cost function values are all lower or equal (≤) to the dominated candidate's values with one value being strictly lower (<).When Pareto domination is applied to the entire population via pairwise comparison, each candidate solution is assigned a rank r.Within one rank, no solution dominates another solution.
However, all solutions of rank r = k 1 dominate all solutions of rank r = k 2 , where k 1 and k 2 are integer values and k 1 < k 2 .The solutions of rank r = 1 are not dominated by any other solution in the population and define what is commonly referred to as the Pareto front.Figure 5 shows the Pareto front evolution up to 90 generations (n).The value of J MO1 and the natural logarithm of J MO2 are shown on the x and y axes, respectively.As the generations go on increasing, the number of candidate solutions on the Pareto front go on increasing, which indicates that a better set of solutions (in terms of the cost functions) are being found with each generation.The Pareto front shifts closer to the x and y axes, indicating that the objective functions are being minimized with each progressive generation.At n = 90, all the solutions on the Pareto fronts have a value of J MO1 , which is more than two times smaller than the baseline case.However, not all solutions offer reductions in the value of J MO2 over the baseline case.The shaded green area in Figure 5 denotes the region in which the value of J MO2 is less than the baseline value of J MO2 = 3.816 × 10 −4 .

Analysis of Solutions
At the last generation (n = 90), every model on the Pareto front was analyzed for its performance in terms of the wall-shear stress and the pressure coefficient.It was found that only those models which reduced the value of J MO2 over the baseline case were able to produce numerically stable solutions.Of these, solutions A and B best represented the flow physics on the blade and in the wake.Although solutions C-F reduced the value of J MO2 over the baseline case and had lower values of J MO1 than solutions A and B, they either predicted an early transition onset location or were unable to provide a strong peak in the separated region in the wall-shear stress profile (i.e., similar to the SO-model in Figure 4a).Absence of a strong peak in the separated flow region of the wall-shear stress profile or predicting an early onset of transition as compared to the DNS data (even with a strong peak in the separated region) resulted in excessive turbulent diffusion and some numerical instabilities in the wake.The model formulation of solution A is The formulation of the f 1a function highlights the leading role of the Π 2 term which is a vorticity scale closely related to the shear-layer vorticity thickness.This is consistent with the original modeling assumption.Π 2 and Π 4 introduce the effect of the freestream turbulence level in LKE production.The Π 7 term is large in the wall proximity and rapidly decays with the distance from the wall.It appears with a negative sign, and thus it tends to suppress LKE production near the wall.This is consistent with the fact that in separated laminar flows the amplification of pre-transitional fluctuations occurs in the separated shear-layer and thus far from the wall.The original formulation of the ψ function employs the wall distance turbulent Reynolds number (i.e., Π 4 ) to control the transition onset.In Equation ( 13), it is replaced by Π 5 .Some recently developed transition prediction frameworks [8,10] share this particular feature with the proposed one.The transition onset parameter, which is originally a constant (i.e., C 4 ), is replaced in Equation ( 13) by a function of Π 1 and Π 6 .The effect is to delay transition as the shear-layer thickness increases (e.g., long bubbles).Additionally, the limiter min( f 1 , f 1a ) was also used as part of the multiobjective calculations.The use of this limiter improved the stability of the solution in the wake region.Due to the use of this limiter, the leading-edge separation prediction did not change with respect to the baseline case.
Figures 6 and 7 show the pressure coefficient and suction-side wall-shear stress, respectively, for the baseline case, DNS data and the solutions A (MO-A) and B (MO-B).Both solutions not only predict the separation plateau feature in the pressure coefficient plot, as opposed to the SO-model, but they also offer an improved prediction in the separation plateau over the baseline case.This is because both the MO-A and MO-B models predict a peak (of smaller magnitude than the baseline case) in the separated flow region of the wallshear stress profile, a necessary characteristic of the LKE modeling framework to predict a separation plateau in the pressure coefficient.These models also offer significantly better predictions of the wall-shear stress over the baseline case.They reduce the magnitude of both peaks (one at transition onset and the other one after turbulent reattachment) in the wall-shear stress by more than two times over the baseline case.The MO-A and MO-B models predict reattachment at x/C ax = 0.987 and x/C ax = 0.988, respectively, which is quite close to the DNS value of x/C ax = 0.986.This is an improvement over the baseline case which predicts reattachment at x/C ax = 0.982.Consequently, these models have a positive impact on the wake mixing prediction.Figure 8 shows the stagnation pressure loss at 25%C downstream of the blade TE (see Figure 2).The stagnation pressure loss (Ω * ) is defined as where P t,1 is the inlet stagnation pressure, P 2 is the exit static pressure, P t (y * ) is the stagnation pressure at the location of interest and y * is a pitchwise normalized coordinate defined as y * = (y − y max )/(y max − y min ).Since the MO-A and MO-B models are able to predict the flow coming off the blade TE quite accurately with respect to the high-fidelity data, the resulting stagnation pressure loss profiles (and thus the wake mixing) are also better predicted than the baseline case.These models improve the prediction of the peak wake loss error and increase the width of the profile, indicating additional diffusion in the wake over the baseline case.It is to be noted that the Boussinesq approximation was used in the wake region for the baseline and MO models.The wake mixing prediction can be further improved by using machine-learned EARSMs.However, further enhancing the wake mixing prediction was not addressed in this work, as this is considered separately in [25,29,30].The 'CFD-driven' canonical-wake EARSM of Akolekar et al. [30] is apt for further improving the wake mixing prediction for this particular case.The canonical-wake model has been formulated such that it can improve the wake evolution when the flow coming off the blade TE is accurately predicted.
In summary, including an additional cost function based on the suction-side pressure coefficient was an immense benefit to formulate models which preserve the integrity of the LKE modeling framework and improve the prediction of the physical characteristics of the flow, not only in the separation region but also downstream of the blade TE.

Conclusions
This study used a multi-objective 'CFD-driven' algorithm to enhance the prediction of separated flow transition of the T106A LPT blade at Re 2is = 100, 000.The existing LKE transition model struggles to accurately predict the wall-shear stress in the separated flow region, and thus the wake mixing downstream of the blade TE is also poorly predicted.Using the 'CFD-driven' multi-objective framework and seven non-dimensional pi groups, model formulations for the transition parameter ψ and the f 1a term in the laminar eddy viscosity are proposed.Initially, a single-objective approach was explored, using the wallshear stress from DNS data in the separation region as a reference.The resulting SO-model was able to predict the wall-shear stress surprisingly close to the DNS.However, as there was no strong peak in the separation region of the wall-shear stress profile (which is a result of turbulence ignition by the transfer term), the prediction of the separation plateau in the pressure coefficient plot worsened over the baseline case.
In order to account for this deficiency of the SO-model, multi-objective optimization was used and an additional cost function which included the suction-side pressure coefficient was introduced.Pareto analysis was used to assess the performance of the models developed.It was found that only those models which reduced the value of the pressure coefficient cost function over the baseline case produced numerically stable solutions.Two of these models were able to preserve the integrity of the LKE framework (i.e., strong peak in the separated region) as well as improved the prediction of physical characteristics of the transition process such as the reattachment location and magnitudes of the peaks in the wall-shear stress profile.Thus, the overall structure of the separation bubble was improved.Since these models were able to closely match the flow coming off the blade TE, the resulting turbulent diffusion in the wake was increased.This resulted in better prediction of wake loss profiles with respect to high-fidelity data.
Overall, a promising framework for multi-objective and multi-expression 'CFD-driven' model development using non-dimensional pi groups was explored.To the authors' knowledge, this is the first known study which makes use of data-driven machine learning tools to enhance the prediction of separation-induced transition in LPTs.Additional analyses and development is required to test this framework and resulting models across different operating conditions and geometries.

Figure 1 Figure 1 .
Figure 1 depicts the methodology of the transition modeling framework used in this study.

Figure 2 .
Figure 2. The T106A blade geometry and computational domain with boundary conditions.

Figure 3 .
Figure 3. (a) Suction-side wall-shear stress from the baseline LKE model and DNS.(b) Zoomed view of the wall-shear stress in the separated flow region, near the blade trailing edge.

Figure 4 .
Figure 4. (a) Suction-side wall-shear stress and (b) pressure and suction-side pressure coefficient for the SO-Model, baseline case and DNS.

Figure 5 .
Figure 5. Evolution of the Pareto front across 90 generations.

Figure 6 .
Figure 6.(a) Pressure coefficient with the various multi-objective optimized models.(b) Zoomed view of the suction-side pressure coefficient in the separated flow region, near the blade trailing edge.

Figure 7 .
Figure 7. (a) Suction-side wall-shear stress with the various multi-objective optimized models.(b) Zoomed view of the suction-side wall-shear stress in the separated flow region, near the blade TE.

Figure 8 .
Figure 8. Stagnation pressure loss profiles 25%C downstream of the blade TE.