Progressive augmentation of Reynolds stress tensor models for secondary flow prediction by computational fluid dynamics driven surrogate optimisation

Generalisability and the consistency of the a posteriori results are the most critical points of view regarding data-driven turbulence models. This study presents a progressive improvement of turbulence models using simulation-driven Bayesian optimisation with Kriging surrogates where the optimisation of the models is achieved by a multi-objective approach based on duct flow quantities. We aim for the augmentation of secondary-flow prediction capability in the linear eddy-viscosity model k − ω SST without violating its original performance on canonical cases e.g. channel flow. Progressively data-augmented explicit algebraic Reynolds stress models (PDA-EARSMs) for k − ω SST are obtained enabling the prediction of secondary flows that the standard model fails to predict. The new models are tested on channel flow cases guaranteeing that they preserve the successful performance of the original k − ω SST model. Subsequently, numerical verification is performed for various test cases. Regarding the generalisability of the new models, results of unseen test cases demonstrate a significant improvement in the prediction of secondary flows and streamwise velocity. These results highlight the potential of the progressive approach to enhance the performance of data-driven turbulence models for fluid flow simulation while preserving the robustness and stability of the solver.


Introduction
Reynolds-averaged Navier-Stokes (RANS) equations are widely preferred over high-fidelity methods like direct numerical simulation (DNS) and large-eddy simulation (LES) for the industrial applications of computational fluid dynamics (CFD) due to their robustness and computational speed.In RANS, the physics of turbulence is predicted by a Reynolds stress tensor (RST) model; hence, the results obtained are dependent on the performance of these model predictions.Despite the popularity of RANS simulations, the common empirical models have been found to have shortcomings [1], particularly in capturing Prandtl's second kind of secondary flow [2].This limitation is accentuated in the most commonly used RANS turbulence models (e.g.two-equation models based on the Boussinesq assumption like k − ε and k − ω) due to their inability to predict complex turbulence anisotropy.
Although the development of reliable RANS turbulence models has remained stagnant for decades [3], the recent advances in data-driven techniques have motivated a new wave of studies aiming at improving the performance of RANS turbulence models [4].The majority of these studies have used available high-fidelity data of the RST values to train a model and improve the predictions obtained from empirical models [5][6][7][8][9][10].These studies have mainly focused on obtaining ways to correct the RST prediction [11][12][13][14][15], or to modify the available empirical models [16][17][18].Taking into consideration a different approach, CFD-driven models [19,20] have shown promising results in finding reliable variants of RANS turbulence modelling.Hence, CFD-driven models can guarantee consistency and robustness in a posteriori results (i.e. by a model-consistent training method [16]), as opposed to other data-driven RANS models (i.e. by a a priori-training method [16]) since they are optimised while performing simulations to guarantee the improvement of the results obtained by new models.
In order to improve these models by a CFD-driven approach, it is necessary to solve a complex optimisation problem.In this regard, the use of optimisation algorithms that aim to minimise or maximise one or more functions in a multi-dimensional parameter space is essential.Some of the commonly used optimisation algorithms include slope followers [21], simplex methods [22], multi-objective evolutionary [23], and particle swarm algorithms [24], among others.However, these algorithms require the evaluation of an objective function, which can be computationally expensive if CFD simulations are involved, especially in cases where a large number of test configurations are required.
To address this issue, development has been made towards the use of a relatively smaller set of simulations to create computationally efficient surrogate models, also known as response surfaces, that can then be fastly optimised [25].Response surfaces are mathematical models that approximate the behaviour of the objective function in the parameter space and can be used to predict the function values at untested configurations.This approach has been applied to various engineering applications, such as optimising complex internal-flow systems based on ultrasonic flow metering [26,27], improving the efficiency of gas cyclones [28], optimising the aero-structural design of plane wings [29], and improving the performance of ground vehicles [30], among others.
There are only a few studies that have investigated the use of optimisation algorithms and data-driven models for the improvement of the RANS turbulence models.Reference [19] combined CFD-driven optimisation with gene expression programming to obtain a correction for the RST modelled by the k − ω SST model.They showed that the CFD-driven model had an improved performance compared to the data-driven model trained on the same case.They concluded that CFD-driven models have a great potential for developing reliable improved RANS models even though their new model is limitedly optimised for wake mixing flow.In another study, Ref. [20] used response surfaces to find the best linear combination of candidate functions to correct the RST values modelled by the k − ω SST model for the cases with flow separation.They also showed that a CFD-driven approach yields models that perform better than the data-driven models trained on the same set of data.In Ref. [31], the authors resort to Kriging and obtained wall models for boundary layer flows subjected to system rotation in an arbitrary direction.The model was shown to predict deviations in the mean flow from the equilibrium law of the wall.In Ref. [32], the authors employed an evolutionary neural network and arrived at numerical strategies for the pressure Poisson equation with density discontinuities.Furthermore, the CFD-driven methodology is extended to a multi-objective optimisation for coupled turbulence closure models by [33].However, the application of surrogate-based optimisation (SBO) and Bayesian optimisation methods [34] to improve complex mathematical tools such as RANS turbulence modelling has not been widely explored.In this regard, SBO and Bayesian optimisation methods show potential in evaluating problems based on design and analysis of computer experiments (DACE) [35] with a low-to-medium number of parameters to optimise.
One of the most important topics in data-driven RANS modelling is the generalisability of the new models for unseen cases [36].It has been suggested that using a multi-case CFD-driven approach to consider different turbulent phenomena during the optimisation process can help with the generalisability of the new model [37].However, even these new models are still specific to either wall-free or wall-bounded flows; therefore, the generalisability problem requires further investigation.
The generalisability is the ability of data-driven models to perform well both with unseen cases (out of the training cases range) and with canonical cases like channel flow.The combination of a conventional progressive approach with data-driven approaches has been proposed to address this issue [38].In the progressive approach, the starting point is a baseline model that is already performing adequately for simple flows and adds more complexity step by step without violating the model's performance for flow cases that the model was calibrated against.In this study, we use the progressive approach to add the capability of secondary flow reconstruction to a linear eddy-viscosity model without violating its successful performance in a channel flow simulation.
Since the CFD-driven technique ensures the consistency of the a posteriori results and the progressive approach aids with the generalisability of the new models, in this study we combine the CFD-driven surrogate and Bayesian optimisation technique with the progressive approach to obtain an explicit algebraic Reynolds stress correction model (EARSCM) to the k − ω SST [39] model for exclusively predicting secondary flows.We likewise investigate the potential of a progressive approach in the development of generalisable CFD-driven RANS models.Since standard linear eddy-viscosity models have difficulty predicting secondary flows [2], we use the Pope's decomposition of RST [40] to add a non-linear term of the RST to the model, and we optimise the new models for the prediction of secondary flows induced by a duct flow with an aspect ratio (AR) of 1 (squared) at bulk Reynolds number of 3500 (depicted in Fig. 1).The new models' boundary-layer prediction is tested on the channel flow to ensure that they do not affect the successful performance of the original model.Considering the generalisability, the new models are verified on duct flow cases with different Reynolds numbers and aspect ratios.Finally, to test the new models against an extremely disadvantageous case regarding accurate flow prediction by RANS models, we verify the models in a case where the secondary flow is induced by roughness patches in a channel flow [41,42] with a nominally infinite Reynolds number.

Methodology
This section presents the structure of the progressive correction model for RANS modelling and the optimisation technique for training the correction model.By using the Reynolds decomposition of velocity and pressure, the Navier-Stokes equation for an incompressible steady flow can be written as [43] where i, j = 1, 2, 3 are indicating the streamwise (x), spanwise (y), and vertical (z) directions, respectively.u i and p are velocity and pressure decomposed to a temporal mean term, indicated by ⟨•⟩, and fluctuations, indicated by • ′ .The kinematic viscosity is denoted by ν, and ρ is the fluid density.
is the anisotropic part of RST, and pressure is modified with the isotropic part of the RST as ⟨P⟩ = ⟨p⟩ + 1 3 ρ⟨u ′ i u ′ i ⟩.In this study, the k − ω SST [39] model is used as a baseline model to be progressively corrected to predict secondary flows.In the standard k − ω SST model, A i j is modelled as where is the strain rate tensor, ν t is the turbulent viscosity which is calculated by values of turbulent kinetic energy (TKE) (i.e.k), and the specific dissipation rate (i.e.ω), which are modelled by the original two-equation model [44].

Progressive EARSCM
We extend the structure of the linear eddy-viscosity model (Eq.3) by following Pope's decomposition [40] of the RST and considering the two first terms of the decomposition, where is the rotation rate tensor, and α (n) are unknown functions of 5 invariants defined as, A comparison between Eq. 3 and Eq. 4 shows that the k − ω SST model is already providing the first term in Pope's decomposition of the RST (i.e., the turbulent viscosity); therefore following the progressive approach, only the second term is trained in this study.The reason behind this decision is to preserve the original performance of the k − ω SST model in the prediction of the turbulent viscosity (i.e., ν t ) while the second basis tensor's coefficient is determined with the exclusive purpose of secondary flow prediction.Hence, the new RST model can be written as where ν t and ω are modelled by the standard k − ω SST model, where Eq. 6 is used for the production term by the Reynolds stress tensor instead of Eq. 3, and the unknown function of α (2) is determined by the CFD-driven optimisation technique.It should be mentioned that the linear part of the new model is treated implicitly as turbulent viscosity, and the non-linear part of the RST is added explicitly to the RANS equations.The assumption behind the progressive approach is that using only the second basis tensor does not affect the incompressible parallel shear flow, either in the momentum equation or via production in the k-equation; therefore, the performance of k − ω SST is preserved in cases where secondary flow is not present.Inspired by a sparse regression of candidate functions (SpaRTA) [14] used for α (n) , we use a set of candidate functions to describe α (2) as where θ i are constant coefficients to be determined by the CFD-driven optimisation process.To achieve a more efficient sparse optimisation, the normalised candidate functions are normalised and defined as where µ i , σ i are the mean and the standard deviation of each candidate function D i , respectively.These statistics are calculated based on high-fidelity data from the optimisation case.Therefore, Eq. 7 is rewritten as, where an optimisation technique determines the coefficients C i based on the performance of the correction model for the reconstruction of the high-fidelity velocity field of an optimisation case.In this study, only 2D canonical flow cases are considered, since they are computationally cost-effective.However, reducing the dimensionality of the optimisation problem and using computational parallelisation is key if a complex 3D case is used in the training process.
Since considering 21 optimisation variables is impractical, making the model unnecessarily complex and exposed to solution instabilities, two approaches are chosen: 1. Selecting only the first m leading candidate functions.2. Reducing the dimensionality of the problem using a statistical technique.
Both of these approaches are compared, and for the purpose of dimensionality reduction, principal component analysis (PCA) is applied on the B i functions to obtain the first m principal components as, where the coefficients a ( j) i are obtained by performing PCA on the high-fidelity data from the optimisation case.It should be mentioned that PCA also determines which features among all the 20 features have a higher importance in the variability of α (2) by determining the a ( j)  i coefficients for each of them.By considering this transformation equation, Eq. 10 is rewritten as, Based on the high-fidelity data of the optimisation case, in Sec.3.1 it is shown that the first two principal components are enough to represent a high percentage of the variability of the dataset (i.e.m = 2).Therefore, two different structures for α (2) are considered: which are described as model I and model II, respectively.Where a unique set of three coefficients C 0 , C 1 , and C 2 are determined by a multi-objective optimisation technique for each of the models.

Optimisation methods
To approach the optimisation problem, a surrogate based on Kriging and Gaussian processes (GPs) is built.The surrogate is based on DACE, using observations obtained from a sequence of CFD solutions and post-processing of results.SBO has shown its advantages for these types of problems, minimising the number of required observations, and allowing the use of goal-seeking algorithms such as Bayesian optimisation [45].To assess the model's performance, two main objective functions are established that involve streamwise velocity and streamwise vorticity, defining a multi-objective optimisation problem.An effective sampling plan is also outlined, and an infill criterion constrained by quality metrics is specified.Finally, the methodology is implemented using the general-purpose software OpenFOAM [46].A flow diagram of the complete optimisation methodology is depicted in Fig. 2.

Objective functions
The principal issue that is assessed in this study is the complete lack of prediction of secondary flows in RANS turbulence models based on the Boussinesq assumption.In the chosen optimisation case, the addition of these secondary flow motions is strong enough to also incur changes in the streamwise direction of the flow.In addition, once the corrections are made, the mean streamwise flow is simultaneously modified.Since the aim of this study is to improve the overall flow prediction, it is important to not disregard possible predictions that could improve secondary flow but worsen streamwise velocity.Hence, two objective functions that physically describe the streamwise flow and the secondary flow prediction are defined and evaluated.
In this regard and to holistically evaluate the flow prediction, the normalised error of both the streamwise velocity and vorticity with respect to high-fidelity data, are evaluated and taken as objectives to be minimised.To yield a fair and accurate single value representing these quantities, the volumetric average of each field is computed as the objective functions, following where is the vorticity of the mean flow, where ε ilk is the alternating unit tensor, HF stands for high fidelity, and k − ω refers to the solution of the standard k − ω SST model.The definition of these functions yields a value of 1 when their output matches the prediction of k −ω SST and 0 when the predictions match the high-fidelity data.Subsequently, to evaluate the overall results of the method, a global fitness function is defined as

Sampling plan
Two main algorithms are used in order to sample observations and data in this study: Monte Carlo and Latin hypercube.Monte Carlo sampling is used when extracting a test subset to compute the surrogate quality.Since the chosen sampling plan for the initial observations is oriented to be space-filling, Latin Hypercube Sampling (LHS) [47] with optimised design by the enhanced stochastic evolutionary (ESE) algorithm [48,49] is used to obtain an initial sample of the surrogate with a minimum number of observations, and simultaneously, representing the real variability of the parameters.This deterministic sampling technique is a type of stratified Monte Carlo that divides each dimension space representing a variable into n 0 sections, and only one point is placed in each partition.The initial sample by LHS with ESE follows n 0 = 30K, where n 0 is the initial number of observations, and K is the number of design variables.

Surrogate construction
In the context of a sparse optimisation problem, surrogates are commonly used to approximate a function y = f (x) based on known observations.In this study, the Kriging method is chosen for generating the surrogate based on CFD-driven observations.Kriging is one of the most common surrogate construction methods due to its well-known implementation, ability to compute uncertainties and speed.Furthermore, it has been proven useful in physical, uncertainty quantification, and engineering applications [27,50].Kriging interpolates the observations as a linear combination of a deterministic term and a stochastic process, which is represented by where f (x) is the surrogate prediction, β is a linear deterministic model, f (x) is the known function, and Z(x) is the realisation of a stochastic process with zero mean and spatial covariance function given by Here, the spatial correlation function R determines how smooth the Kriging model is, how easily the response surface can be differentiated, and how much influence the nearby sampled points have on the model.In this study, the spatial correlation is defined following the squared exponential (Gaussian) function as where the correlation scalar θ l is used to define the variance of a Gaussian process at each point, with higher values indicating a stronger correlation between points.By maximising the maximum likelihood estimation, optimal values for hyper-parameters as θ l , mean, and standard deviation can be found [51].

Quality metrics
The quality of the surrogate refers to how accurately it approximates the true function being modelled.In order to evaluate this quality, a number of random test observations n t = n 0 x − → y, are taken.where n t is a function of the number of initial observations.Once the test observations have been generated, they are used to compare the surrogate's predictions to the true values of the function being modelled.The root-mean-squared error (RMSE) and Pearson's correlation coefficient squared (r 2 ), defined as are two common metrics used to evaluate the quality of surrogates.While the RMSE measures the difference between the surrogate's predictions and the true values, r 2 measures the strength of its linear relationship.A high RMSE and a low r 2 indicate that the surrogate is performing poorly and needs to be refined, while a low RMSE and a high r 2 indicate that the surrogate is performing accurately and can be used with confidence.Hence, these values serve as reference points to determine whether the surrogate model is accurate enough to be used in the optimisation process.In this study, thresholds of 0.2 for RMSE and 0.8 for r 2 are established as indicators of good surrogate quality, based on previous studies [52,53].

Infill space exploration
Efficient global optimisation (EGO) based on Bayesian optimisation strategies is used to improve the accuracy of surrogate models and decrease overall uncertainty [45,54].This is achieved by exploring the surrogate beyond the initial sampling.EGO is a well-known algorithm that employs both local and global searches to find the optimal solution by means of the expected improvement (E[I(x)]) function as a key metric to direct its search.The function calculates the potential improvement that can be obtained by evaluating a new observation point based on the current best solution and the overall uncertainty of the surrogate model [55].The implementation of this function provides the necessary sparsity-promoting behaviour to accordingly explore the design space in regions where their initial sampling was not sufficient.The expected improvement function is defined as where f min = min Y, and Φ(•) and ϕ(•) are respectively the cumulative and probability density functions of N(0, 1), following the distribution N(µ(x), σ 2 (x)).Using EGO to explore the surrogate beyond the initial sampling aids in decreasing the overall uncertainty and improving the precision of the surrogate model, especially in the unexplored and extreme regions of the design space.Consequently, the algorithm can effectively balance local and global searches and locate the optimal solution more proficiently.Hence, the following sampling point is determined by In order to fully explore the design space, new data is collected for each objective function separately.This involves taking one sample of n EGO = n 0 x − → y for each objective, resulting in a total of 180 new sampled points (n 0 for each of the objectives).With the addition of the initial LHS, the total number of observations sums to 270 per model, yielding a cost-efficient computational problem to be analysed by Kriging.However, since the design space explored is large and the use of Bayesian optimisation might sample observations in very close proximity, the gradients between points may be high, yielding possible noise.Therefore, the surrogate is regularised following the methodology by [56], where noise is considered to be Gaussian distributed and its homoscedastic variance is evaluated based on the performed observations to obtain a smooth response surface.

Multi-objective solution
The Kriging method offers a significant benefit since it allows the application of algorithms that can thoroughly explore the objective space.Consequently, the objective of this study is to find the optimal solution for the surrogate model.To achieve this, the multi-objective evolutionary algorithm (MOEA) NSGA-II has been chosen.This algorithm is known for its versatility, fast and efficient convergence, and ability to handle non-penalty constraints.It also has a wide-range search for solutions, which implies that it can explore the objective space more thoroughly and is less likely to get grounded in local optima [57].By using the MOEA NSGA-II algorithm, the improved solutions are ensured to be non-dominated, meaning that they are not worse than any other feasible solution in all objective functions.This also ensures that the solutions avoid local minima and are highly likely to discover the global optimum for the sampled response surface.This makes the algorithm able to search a large portion of the objective space and find the best possible solution for the surrogate model while avoiding suboptimal solutions.

High-fidelity data
Since the main purpose of this study is to obtain an EARSCM for capturing the secondary flow, the DNS data of a canonical case with this flow characteristic is chosen.A duct flow case of AR = 1 and bulk Reynolds number of Re b = 3500 is used for the training process (depicted in Fig. 1).The DNS data is obtained from Ref. [58] curated by Ref. [59].Following the progressive approach, the trained EARSCMs are tested on two cases of channel flow with friction Reynolds number of Re τ = 395 and Re τ = 5200 for which data is obtained from Refs.[60] and [61], respectively.Considering the generalisability of the new models, we likewise test them on three unseen cases containing secondary flows, including a duct secondary flow case with AR = 1 and higher bulk Reynolds number of Re b = 5700 [62], a duct secondary flow case with AR = 3, and a lower bulk Reynolds number of Re b = 2600 [62], and a roughness-induced secondary flow with nominally infinite Reynolds number (more information about the geometry and properties of this case is available at Refs. [15,42].

Results and discussion
In this section, the results are presented in two subsections.In the first one, the results of the optimisation process and training of the two best EARSCMs are presented, whereas, in the second, the performance of the trained models is evaluated on the test cases with different geometries, Reynolds numbers, and boundary conditions.Associated contours and velocity profiles accompany all results which are similarly described and discussed, followed by their corresponding figures.It is essential to qualitatively indicate the original performance of k − ω SST against high-fidelity data prior to any modifications.For both fidelity levels, the results for the square duct case show the progressive deceleration of the streamwise flow as it approaches the walls of the duct, exhibiting higher gradients in the near-wall regions as a consequence of the boundary layer and the fully developed turbulent flow.Whereas this behaviour can be predicted by k − ω SST, the presence of a secondary flow predicted by DNS is completely neglected by the RANS model (Fig. 3), as expected [2].Two antisymmetric rotational regions are generated diagonally at the duct's vertex, preserving their symmetry at the 4 quadrants of the duct.This secondary flow is strong enough to drive the streamwise flow towards the duct vertices, skewing the streamwise velocity field in the domain through an impinging-like flow behaviour.Since k − ω SST does not predict these rotational motions, the streamwise velocity distribution is not skewed and the spanwise and vertical components of the velocity are zero.Figure 4 compares the barycentric plots representing the shape of RSTs for the optimisation case.The location of each point is calculated based on the eigenvalues of the normalised anisotropic part of the RST (more information about barycentric plots available at Refs. [63,64]).While k − ω SST yields all results within the plain-strain limit as expected from a linear eddy-viscosity RANS model, DNS shows a more complex stress anisotropy distribution as a combination of states toward one-component turbulence fluctuations x1 c (also known as rod-like or cigar-like turbulence) and an offset from the plain-strain limit with a tendency towards isotropic (or spherical) turbulence x3 c .

Optimisation case: EARSCMs
Since k − ω SST is not able to predict the secondary flow of the optimisation case, CFD-driven optimisation is applied following the models described in Eqs.13a and 13b focusing on predicting secondary flows and correct streamwise velocity.
Regarding model II, PCA is applied to the set of candidate functions (listed in Eq. 8) using the high-fidelity data of the optimisation case.Figure 5 presents the PCA results and shows that only the first 2 principal components (φ 1 and φ 2 ) yield an explained variance ratio of 0.90 (Fig. 5a), where φ 1 is the main principal component explaining the variability.Figure 5b presents the coefficients corresponding to the two first principal components used for model II.As mentioned in Sec. 2, the two first principal components for model II are chosen.To compare model I in an analogous manner, the two first leading candidate functions are likewise chosen where the coefficients of these two models are determined by surrogate and Bayesian optimisation.
The surrogate model is generated using the RANS simulations results to find the optimal values of the optimisation variables.Regarding the optimisation process and since cases with high values of the objective functions are not of interest, the MOEA optimisation is performed by imposing the constraints j 1 ≤ 0.5 and j 2 ≤ 0.4, therefore, the analysed non-dominated solutions only take place in the regions of highest interest and potential global minimisation.To yield a smooth non-dominated solution front, 5000 samples are required by the MOEA in this study.This showcases the cost advantages of surrogates, where 270 CFD observations are required, yielding a cost reduction of 18.4 times.
Regarding both models' results, an optimal design space of [−2, −1, −1] ≤ C i ≤ [0, 1, 1] is chosen for model I, and an optimal design space of [−1.75, −0.25, −0.25] ≤ C i ≤ [−1.25, 0.25, 0.25] is chosen for model II after various progressive surrogate iterations.The iterations re-define the design space limits by analysing the non-dominated values of the coefficients after optimisation.If these values reach the imposed limits, the objective space is reset, and the design space is adapted.
Model I. Contour plots in Fig. 6b illustrate how the optimisation variables affect the objective functions.On the one hand, the surrogate visualisation for j 1 (Fig. 6a) shows a quasi-linear tendency for all variables.The global minimum for j 1 is shown towards negative values for all C i without clear visualisation of extrema.On the other hand, there is a clear global minimum seen for j 2 values, shown by a non-linear surrogate for C 0 and C 1 , where quasi-linear behaviour is predicted by the relationship between C 1 and C 2 (Fig. 6b).It is important to highlight that, although the lowermost values for the coefficients predict an improvement for j 1 , the tendency is not completely followed by j 2 prediction.This implies that a more accurate prediction of the streamwise flow does not necessarily correlate with the secondary flow prediction in this case.Regarding the surrogate quality (Fig. 7a), and Pareto front prediction and numerical validation (Fig. 7b), the surrogate predicts the test set with great accuracy for the whole design space tested, which is likewise reflected by a highly accurate Pareto front prediction, where all numerical tests simulated are able to predict the objective space within 1σ uncertainty bounds for both objective functions.Model II.In contrast with model I, the surrogate visualisation for both objective functions (Fig. 8) in model II shows a non-linear prediction for all variables.This behaviour is somewhat expected due to the added complexity of the model with the PCA.The global minimum for j 1 does not have a straightforward location.Instead, there are diverse regions where the global extrema can be found.Specifically, negative values of C 0 , and near-zero values of C 1 and C 2 predict generally lower values of j 1 and j 2 (Figs.8a and 8b).Similarly to model I, although the objective spaces for both functions follow similar tendencies, there are visible differences in the predicted locations of the minima, reflecting once more that a more accurate prediction of the streamwise flow does not necessarily correlate with the secondary flow prediction.Regarding the surrogate quality (Fig. 9a), and Pareto front prediction and numerical validation (8b), the surrogate predicts the test set with high accuracy for the whole design space tested.These lower-confidence regions are reflected in the non-dominated solutions nor in their numerical validation.Since the expected improvement function prioritises good confidence in the near-extrema region and the complexity of the model is higher, a certain level of uncertainty is expected in far-away regions of the global minimum.Nevertheless, the surrogate shows a good level of confidence at the Pareto front (within 2σ in Fig. 9b). Figure 9: Surrogate quality against test set (Fig. 9a), and Pareto front validation with 2σ uncertainty bounds (Fig. 9b) for model II.

Comparison of models
The final results for each model are shown in Table 1, where the best cases per model are shown, followed by their C i values and surrogate quality parameters.Both models significantly improved the standard k − ω SST by 61.3% and 60.2% for model I and II respectively.Differences in performance for the objective values are considered negligible.As seen in Fig 7a and 9a, the surrogate quality is slightly worse for model II.This is due to the higher complexity added by PCA and the inclusion of higher order terms in α 2 , yielding more non-linear behaviour in the design space and, thus, a higher RMSE and lower r 2 for model II.Regarding qualitative results for the optimisation case, both models are able to accurately predict the direction and symmetry of the secondary motions while improving the streamwise flow prediction.However, the streamwise vorticity error is generally higher in the near-wall regions (Fig. 10).As a consequence of the prediction of a correct secondary flow direction, the streamwise flow prediction likewise improves.Regions of slight overprediction of ⟨u 1 ⟩ are seen in the near-wall vicinity as well as the middle section of the computational volume, whereas regions of ⟨u 1 ⟩ underprediction are seen in the bulk flow close to the channel vertex.Regarding the prediction of ⟨ω 1 ⟩, antisymmetric predictions are seen with respect to the diagonal symmetry line with a general over and underprediction in the near-wall region.
As expected from the objective function results, qualitative differences between models (Fig. 11) are not clearly seen.For the streamwise velocity prediction, no significant differences can be seen in the contours of the error function, whereas for the vorticity prediction, slight variations in the error can be seen between models without a significant impact on the overall performance of the models.In summary, both models are able to predict the secondary flow with high accuracy.Concerning the quantitative analysis of the velocity profiles at = [0.25,0.5, 0.75], an overall and significant improvement in the prediction of both models against standard k − ω SST, can be seen.Although the prediction of both models displays minor differences between each other, they both predict with high accuracy the DNS data.Some light discrepancy can be seen in the magnitudes at near-wall regions, where the gradient of ⟨u 2 ⟩ is high.Nevertheless, gradients of velocity are accurately predicted, only displaying a slight underprediction in the velocity magnitude in near-wall regions.The turbulence shape of the developed models is represented in Figs. 12 and 13, where the reliasability for both models is conserved.A certain level of anisotropy is added to the models, yielding results in the vicinity of the plainstrain limit.Both models predict a very similar turbulence shape while model II slightly yields results farther away from the plain-strain limit.Nonetheless, results are distant from x1 c and DNS anisotropy.Although these results are expected, the integration of the whole Pope's decomposition with its 10 tensors and a more thorough definition of B could be able to further improve the RANS anisotropy prediction.Figure 14 presents the profiles of RST components obtained by new models and shows that new correction models improve the prediction of RST components.It should be mentioned that the Reynolds stress tensor is calculated as The results of R i j show that the introduced correction does not change the original R 11 prediction by k − ω SST while the 2 other principal components of the tensor (R 22 and R 33 ) are predicted in better agreement with highfidelity data.Predictions of the off-diagonal components of R i j are overall displaying an improvement in prediction, highlighting the non-zero prediction of R 23 , which shows a certain level of mismatch with high-fidelity data in the near-wall regions.
Qualitatively and quantitatively, both models improve the prediction and agreement with high-fidelity Reynolds stresses while not showing significant differences in the prediction of Reynolds stresses between models.

Verification and generalisability on test cases
In this subsection, the best-found models are tested against diverse canonical cases in which secondary flow is an important component as well as cases where secondary flow is not present.This verification is performed in order to provide a performance overview of the models' generalisability, stability, and robustness to preserve canonical flow features.Firstly, validation of the channel case and the boundary layer is performed at different Reynolds numbers.Secondly, the models are tested against diverse cases of ducts with different aspect ratios and Reynolds numbers.Lastly, the models are tested against a wall-modelled nominally infinite Reynolds number case in which secondary flow is roughness-induced.Similarly as in section 3, all results are shown and discussed with their respective velocity contour plots with stream functions and qualitative flow analysis of the velocity profiles.

Channel flow and law of the wall
The addition of T (2)  i j and further modifications in this study must not destabilise k − ω SST and yield unphysical results.Therefore, both models are tested in a channel flow case at friction Reynolds numbers 395 and 5200 to verify that the law of the wall is preserved.The canonical channel flow cases are driven by a constant pressure gradient in order to match the friction Reynolds number from DNS.
Results shown in Fig. 15 depict the streamwise velocity as u + = u/u τ and k + = k/u 2 τ in function of y + = u τ y/ν, where u τ is the friction velocity.The results show a consistent prediction of the law of the wall with k − ω SST without introducing any noticeable changes.Both models yield the same results for both the streamwise velocity and the turbulence kinetic energy, and no improvements or diminishments in the prediction of these variables are seen compared to standard k − ω SST. Figure 15: streamwise velocity and turbulent kinetic energy profiles for channel flow: Re τ = 395 (Fig. 15a), and Re τ = 5200 (Fig. 15b).High-fidelity data obtained from Ref. [60].

Duct cases
To showcase the generalisability of the models, diverse duct cases are tested.The first case is the flow through a duct of AR = 1 and a considerably higher Re b = 5700 compared to the optimisation case.
A qualitative depiction of the results is shown in Fig. 16, where both models are able to predict the gradients of the secondary flow and improve the prediction of the streamwise velocity compared to k − ω SST.In agreement with the baseline results of the models, the magnitude of the secondary motion is weaker than the high-fidelity data although the motion's antisymmetry and location of the vortices centres are predicted accurately.
Regarding the velocity profiles of the case (Fig. 17), a similar trend is observed where the gradients and magnitude of the profiles are predicted with high accuracy in both models, only displaying some discrepancies in the wall-near regions, where the high-fidelity data yields slightly higher velocity gradients.In order to verify the performance of the models in diverse cases with similar features, the models are likewise tested against a duct of AR = 3 and a slightly lower Re b = 2600 compared to the optimisation case.
The qualitative results shown in Fig. 18 indicate a similar performance for both models: gradients, symmetry, and direction of secondary flow are predicted accurately with slight inaccuracies regarding the secondary flow magnitude in the near-wall regions.The location of the vortices centres is likewise predicted with high accuracy, agreeing with high-fidelity data (Fig. 18).In terms of the quantitative analysis and the evaluation of the velocity profiles, the results shown in Fig. 19 likewise follow a similar trend compared to previous results in this study.The region in the vicinity of the duct vertex shows some discrepancy in the velocity prediction.Nevertheless, the bulk flow away from the near wall is accurately predicted.

Roughness-induced secondary flow case
Finally, the verification of the models is tested on a roughness-induced secondary flow case.The case is based on the studies by [41,42], and is characterised by displaying a nominally infinite Reynolds number.It should be noted that due to the very high Reynolds number, the use of wall models is inevitable with current hardware, therefore, the atmospheric wall models [65] have been used.The secondary flows generated in this case are recognised as Prandtl's second kind of secondary flow, similar to the square duct flow [2].The roughness-induced case is chosen to showcase the models' performance in a more complex and challenging to predict case, mostly by two-equation RANS turbulence models.
For a better analysis of roughness-induced secondary flow, dispersive velocity components are defined as, where ⟨u ′′ i ⟩ is the dispersive velocity components, and ⟨ũ i ⟩ is the spatial spanwise-averaged mean velocity.Following previous verification, the qualitative results of the streamwise dispersive velocity (⟨u ′′ 1 ⟩) are shown in Fig. 20.On the one hand, results show that standard k − ω SST does not predict any secondary motion that drives the flow to display a high-momentum path on top of the high-roughness patch and a low-momentum path at y/H = 0, as the high-fidelity data predicts.On the other hand, the enhanced models are both capable of predicting the secondary motion, although at a lower intensity compared to high-fidelity results.The direction of the rotation (clockwise) is correctly predicted although the vortex centre does not visibly match the high-fidelity counterpart.Even though the predicted secondary flow is not strong enough to move the low-momentum path to the top of the low-roughness patch (y/H = 0), the existence of a weak high-momentum path is visible on top of the high-roughness patch.Since the source of vorticity production lies at roughness-heterogeneity existing at the bottom wall (more information at Ref. [42]), the authors believe that further investigations about wall models regarding the secondary flow can help a better prediction of the secondary flow.Quantitive results reflect the previous qualitative analysis.Figure 21 shows the velocity profiles at y/H = [π/8, π/4, 3π/8], where the improvement of the models can be seen.In contrast to other verification cases in this study, both models' predictions improve the performance of standard k − ω SST, however, the secondary flow intensity is weaker compared to the high-fidelity data.The predicted tendencies are correct, although there is a discrepancy in the gradients and magnitudes predictions relative to previous verification cases.
linear eddy-viscosity model with the ability to predict secondary flows while maintaining its effective performance in canonical flow cases.The enhancement is based on the introduction of explicit algebraic Reynolds stress correction models with two levels of complexity: introducing a linear combination of the first two candidate functions and introducing a linear combination of the first two principal components obtained from PCA on all candidate functions.These two models are optimised to reach the best prediction of both secondary flow and streamwise flow in the optimisation case of a duct flow with an aspect ratio of 1 and Re b = 3500.
Considering the progressive approach, the enhanced models are tested against channel flow cases at different friction Reynolds numbers, where the enhanced models preserve the original performance of the k − ω SST model in a successful reproduction of the velocity and TKE profiles.
For the purpose of generalisability investigation, new models are tested against verification duct-flow cases with different Reynolds numbers and different aspect ratios.Both enhanced models show a successful prediction of secondary flow and improvement of streamwise velocity prediction in the unseen cases.The third and final verification case is a roughness-induced secondary flow case at a nominally infinite Reynolds number, where the enhanced models are able to predict the secondary motion, although at a lower intensity compared to high-fidelity data due to the limitations of wall modelling for this type of cases.
Overall, the results, summarised in Table 2, show that the developed models perform substantially better than the standard k − ω SST model in all verification cases.The enhanced models are able to predict secondary flow features, showing global improvements in the velocity field prediction.Both models show similar performance in successfully predicting secondary flows and streamwise velocity.Therefore, more invariants combined with PCA do not necessarily improve the performance of the models in the 2-dimensional cases considered in this study.Nevertheless, considering scaling up the complexity of the model in 3-dimensional cases, the combination with PCA dimensionality reduction could further improve the predictions.It should be noted that the optimisation process and the nature of the EARSCMs yield a certain level of turbulence anisotropy while always maintaining solution robustness and stability.This fact results in some discrepancies between the EARSCMs and high-fidelity data in the predicted turbulence shape.

Case
Model These findings suggest that the use of CFD-driven optimisation with surrogate modelling and Bayesian optimisation is a robust approach to enhance linear-eddy-viscosity turbulence models to predict more complex physics.The enhanced EARSCMs developed greatly improve the prediction of turbulence anisotropy in wall-bounded flows, especially in cases where the standard k − ω SST presents difficulties in accurately predicting the secondary flow.The progressive nature of this development as a first step focuses on the prediction of secondary flows, allowing further improvement of the models by adding more predictive physics and verifying them on other test cases.

Figure 1 :
Figure 1: Flow characteristics of a duct flow with AR = 1 and hight and width of 2H.Note the formation of secondary motions and their symmetry.

Figure 2 :
Figure2: Design optimisation strategy employed involving an initial sampling plan based on Latin hypercube sampling (LHS) which is solved by CFD.Subsequently, an initial surrogate is constructed using Kriging.Since it is likely that the initial sampling may not evaluate the surrogate in extrema, Bayesian strategies based on efficient global optimisation (EGO) and the evaluation of the expected improvement (E [I(x)]) function are then applied to further explore the surrogate and improve its quality.Finally, once the quality requirements have been met, a genetic algorithm is employed to search for non-dominated optima on the surrogate model.

Figure 3 :
Figure 3: Contours of streamwise velocity and secondary flow streamlines for duct flow optimisation case at AR = 1 and Re b = 3500 for RANS k − ω SST (left) and DNS (right) data.DNS data obtained from Ref. [58].The dotted white lines denote the location of the velocity profiles evaluated throughout this study to perform a quantitative comparison.It should be noted the varying thickness of the secondary flow streamlines denotes their magnitude.

Figure 4 :
Figure 4: Barycentric representation of RSTs' shape for k − ω SST (left) and DNS (centre) for training duct flow case with AR = 1 at Re b = 3500.DNS data obtained from Ref.[58].The plain-strain limit is represented by the dotted line, and the colour representation (right) follows the work by[64].
Explained variance ratio of each principal component.Coefficients associated with the first two principal components.

Figure 5 :
Figure 5: PCA results of the candidate functions calculated by DNS data of the optimisation case.

Figure 6 :
Figure 6: Surrogate contour visualisation for model I.It should be noted that each pair of variables is shown by holding the other variable at the mid-value of their range.
Surrogate for j 2 .

Figure 8 :
Figure 8: Surrogate contour visualisation for model II.It should be noted that each pair of variables is shown by holding the other variable at the mid-value of their range.

Figure 10 :
Figure 10: Qualitative streamwise and secondary flow prediction of the best cases for each of the models.Left: streamwise flow contours with stream function lines depicting secondary flow prediction and direction.Centre: streamwise flow error compared to DNS data.Right: streamwise vorticity error compared to DNS data.It should be noted that the stream function line thickness denotes the secondary flow intensity.In the colourbar, ϕ indicates either ⟨u 1 ⟩ or ⟨ω 1 ⟩.

Figure 11 :
Figure 11: Profiles of velocity components for duct flow case with AR = 1 and Re b = 3500.High-fidelity data obtained from Ref. [58].

Figure 12 :Figure 13 :
Figure 12: Barycentric map showing the physical reliasability and turbulence anisotropy of the developed models following the colour representation by [64].

Figure 14 :
Figure 14: Profiles of Reynolds stress components for duct case with AR = 1 and Re b = 3500.High-fidelity data obtained from Ref. [58].

Figure 16 :Figure 17 :
Figure 16: Streamwise flow contours with stream function iso-lines depicting secondary flow prediction and direction: Duct flow case with AR = 1 and Re b = 5700.High-fidelity data obtained from Ref. [62].

Figure 18 :
Figure 18: Streamwise flow contours with stream function iso-lines depicting secondary flow prediction and direction: Duct flow case with AR = 3 and Re b = 2600.High-fidelity data obtained from Ref. [62].

Figure 19 :
Figure 19: Profiles of velocity components for duct case with AR = 3 and Re b = 2600.High-fidelity data obtained from Ref. [62].

Table 2 :
Objective function results for optimisation case and all verification cases tested.Results are presented by the case abbreviation and characteristics: Duct flow cases as DF AR, Re and roughness-induced secondary flow as RI Re .It should be noted that DF 1, 3500 is the optimisation case.J = 1 indicates the k − ω SST original prediction and J = 0 indicates a match of the high-fidelity data.

Table 1 :
Objective function results, coefficients, and quality for best models and surrogates of each approach