Customized data-driven RANS closures for bi-fidelity LES-RANS optimization

Multi-fidelity optimization methods promise a high-fidelity optimum at a cost only slightly greater than a low-fidelity optimization. This promise is seldom achieved in practice, due to the requirement that low- and high-fidelity models correlate well. In this article, we propose an efficient bi-fidelity shape optimization method for turbulent fluid-flow applications with Large-Eddy Simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) as the high- and low-fidelity models within a hierarchical-Kriging surrogate modelling framework. Since the LES-RANS correlation is often poor, we use the full LES flow-field at a single point in the design space to derive a custom-tailored RANS closure model that reproduces the LES at that point. This is achieved with machine-learning techniques, specifically sparse regression to obtain high corrections of the turbulence anisotropy tensor and the production of turbulence kinetic energy as functions of the RANS mean-flow. The LES-RANS correlation is dramatically improved throughout the design-space. We demonstrate the effectiveness and efficiency of our method in a proof-of-concept shape optimization of the well-known periodic-hill case. Standard RANS models perform poorly in this case, whereas our method converges to the LES-optimum with only two LES samples.


Introduction
Numerical fluid-dynamic shape-optimization is an increasingly central tool in engineering practice. Typically Reynolds-Averaged Navier-Stokes (RANS) is used as the fluid model, due to its acceptable computational cost and accuracy. However in many important situations, the physics demands scale-resolving simulations of turbulence such as large-eddy simulation (LES). This includes, for instance, designs where separation occurs; junctions flows with complex turbulence behaviour (Simpson, 2001;Belligoli et al., 2019); and novel boundary layer flow control applications. LES is often applied for a final analysis and validation, but its high computational cost precludes its use within the design loop.
Multi-fidelity optimization (MFO) methods are supposed to address exactly this dilemma. In MFO, low-fidelity but cheap-to-evaluate physics models are used to explore the design-space rapidly; their predictions are then corrected by a few expensive high-fidelity simulations. One such class of methods are the multi-fidelity surrogate modelling methods , combined with a suitable sampling criteria (Jones et al., 1998).
There has been significant work on multi-fidelity surrogate modelling methods, almost all based on Gaussianprocess models for the surrogate (i.e. Kriging), differing mainly in the manner in which the high-fidelity correct the lower. Haftka (1991) and Chang et al. (1993) developed a variable-fidelity Kriging model using a multiplicative "bridge function" to correct the low-fidelity model. Gano et al. (2005) developed a hybrid bridge function method, which uses a second Kriging model to scale the low-fidelity model. Han et al. (2013) combined gradient information and a generalized hybrid bridge function. Cokriging was originally proposed in the geostatistics community by Journel and Huijbregts (1978) and then extended to deterministic computer experiments by Kennedy and O'Hagan (2000).  proposed an improved version of cokriging, which can be built in one step, and a hierarchical Kriging (HK) model , appears to be as simple and robust as the correction-based method and as accurate as cokriging method.
However, the success of all these approaches relies heavily on the correlation between the high-and low-fidelity models (Han et al., 2020;Zhang et al., 2018). Since models with higher fidelity typically have also significantly higher numerical costs, in practice the high-fidelity model is chosen as the cheapest model that predicts the relevant physics. Necessarily, the low-fidelity model will lack some of the relevant physics, leading to a poor correlation between the models, which limits the practical application of MFO methods. To this end, Gómez et al. (2020) proposed a novel multi-fidelity optimization method, injecting LES Reynolds stress statistics into RANS model to improve them, and then using both in a bi-fidelity optimization.
In separate developments of the past few years, data-driven methods for turbulence modelling based on supervised machine learning have been introduced to improve RANS predictions when reference data is available (Xiao and Cinnella, 2019;Duraisamy et al., 2019;Durbin, 2018;Kaandorp and Dwight, 2020). Notably, Parish and Duraisamy (2016) introduced a local multiplicative term to correct the turbulence production in the -equation of the − model. This term was chosen by solving an inverse problem to match LES reference data, and was then used to train a Gaussian process to correct the baseline turbulence model. Ling et al. (2016) trained a deep neural network to predict turbulence anisotropy and replace the baseline turbulence model. An alternative is to use random-forests for the same task (Wang et al., 2016;Kaandorp and Dwight, 2020). Even though these approaches generalize the linear eddyviscosity concept, they generate complex black-box closure models. More promising are approaches that generate compact explicit expressions for models. Notable are gene-expression programming (GEP) Sandberg, 2016, 2017), and deterministic sparse regression of turbulence anisotropy (SpaRTA) (Schmelzer et al., 2020). These methods generate models that can be rapidly implemented in existing CFD codes and evaluated at every iteration of a RANS solution, and potentially inspected to identify the physical mechanisms influencing the flow. The SpaRTA approach has shown the ability to consistently reproduce specified LES mean-flows in a RANS solver, using only three to five additional non-linear closure terms. It is important to understand that -at their present level of developmentthese methods do not generate general-purpose turbulence models, but rather deliver corrections for flows similar to the flow they were trained on. Thus, for flows with a high degree of similarity, e.g. modifications in the geometry, the methods generate models effectively customized to the flow at hand.
The key original contribution of this paper is a highly-efficient bi-fidelity optimization procedure, using LES as the high-fidelity model, and data-enhanced RANS as the low-fidelity model. In outline, we: 1. Perform a single LES simulation at the baseline geometry 0 , = 0; 2. Use the full-field LES data to generate a customized RANS model matching the LES at ; 3. Sample this custom RANS model throughout the design-space (cheap); 4. Combine custom-RANS and LES results in a bi-fidelity surrogate; 5. Perform a new LES at the point of maximum expected improvement +1 (expensive); 6. Generate new custom RANS model(s), based on all LES data available so far; 7. ← + 1, goto 3.
The key observation is that the custom RANS model correlates well with LES in a region of the design-space for which the LES training data is informative. Provided the physical processes occurring do not significantly change, the data-driven RANS remains a viable model. Also, whereas typical multi-fidelity methods use only the objective function from the high-fidelity simulation, our approach uses much more information. There are many variations on the basic pattern described above, for example: Which LES samples are used for training in Step 6?; What criteria are used for sampling RANS in Step 3, and LES in Step 5?.
In the following we demonstrate the method on a proof-of-concept optimization problem: a generalized periodichill (P-H) geometry, in which the steepness of the hills is varied (Gómez, 2018;Gómez et al., 2020). The baseline periodic-hill flow (as proposed by Mellen et al. (2000) based on an experimental study by Almeida et al. (1993)) is notoriously for poor performance of essentially all standard RANS models, due to its sensitivity on the flow separation location and large-scale low-frequency dynamics. We use a hierarchical Kriging bi-fidelity surrogate model, with an EGO-like sampling procedure for the LES (Jones et al., 1998;Zhang et al., 2018). The custom RANS model is sampled on a uniform grid, and updated at each step with the latest LES data. To achieve the LES optimal solution we require only two LES samples (and one further sample to verify the result). Our method therefore offers a promising path towards LES-quality optimization with only a handful of LES samples.
The paper is structured as follows: in Section 2 we describe the computational setting of the baseline RANS and LES simulations. Section 3 briefly describes SpaRTA, the machine-learning method used in this work to generate datadriven RANS models. SpaRTA is demonstrated on the baseline periodic-hill case (Mellen et al., 2000). Our proposed bi-fidelity optimization procedure is described in detail, and applied to the generalized P-H optimization case with the hill width as a design variable in Section 4.

Baseline incompressible LES and RANS simulations
The solvers used here for LES and RANS are both finite-volume methods, but otherwise numerically distinct -we describe both briefly. For the standard periodic-hill test-case we compare our LES results with the experimental data of Rapp (2009) and well-resolved LES reference data from Breuer et al. (2009), and demonstrate excellent agreement. We show that RANS results obtained with a standard -SST model are poor, which is in agreement with literature.

Baseline Reynolds-Averaged Navier-Stokes with -SST
Assuming incompressible flow with constant fluid density equal to unity, the steady RANS equations are where with , ∈ {1, 2, 3} are the components of the mean-flow velocity, is the mean pressure, and is the kinematic viscosity. The volume forcing serves to drive the flow through the doubly periodic domain. The effect of turbulence on the momentum equation is confined to the Reynolds stress tensor , which must be modelled. In this paper, we use the popular -SST model as a baseline model. The model defines transport equations for the turbulence kinetic energy ∶= 1 2 , and the specific turbulence dissipation rate : The eddy viscosity is modelled as , and the anisotropic part of Reynolds-stress tensor is modelled by the Boussinesq assumption and the isotropic part is absorbed into the pressure. All remaining terms and coefficients are omitted for brevity, see (Menter, 1994) for details. The computational mesh used for RANS simulation of the periodic-hill is two-dimensional, with 120 × 130 cells in stream-wise and wall-normal directions respectively. Periodic boundary conditions are set at inlet and outlet, and no-slip conditions on the walls. Volume forcing is used to drive the flow, with a PID controller used to achieve the target velocity. The simulation is performed at Re h = 10595 based on the hill height ℎ using second-order accurate SIMPLE solver. The pressure is solved using geometric-algebraic multigrid with a Gauss-Seidel smoother, while , and are solved using a Diagonal-Incomplete-LU-preconditioned BiCG method.

Large-eddy simulation for the periodic-hills
The high-fidelity and "ground-truth" model in this work is wall-modelled LES with our in-house finite-volume solver INCA. We use the incompressible staggered-grid version of the solver with a block-Cartesian background mesh and the conservative immersed boundary method of Meyer et al. (2010) for representing the geometry. Our LES follow a holistic modelling approach, where the numerical discretization and the subgrid-scale (SGS) turbulence closure are fully merged: The adaptive local deconvolution method (ALDM) of Hickel et al. (2006); Hickel and Adams (2007) is a nonlinear finite-volume discretization scheme tailored for LES of turbulent flows. Optimum model form and discretization parameters were learnt by a physics-informed genetic algorithm using turbulence data, spectral analysis and constraints from turbulence theory (Hickel et al., 2006). The near-wall turbulence is modelled rather than resolved (Hickel et al., 2012;Chen et al., 2014), dictated by the available computational resources and the Reynolds numbers in our case. A third-order explicit Runge-Kutta scheme is used for time integration and the pressure-correction Poisson equation is solved using BiCGstab with incomplete LU and algebraic-multigrid as preconditioners. Overall, the our LES solver is very similar to that used by Hickel et al. (2008) for wall-resolved LES of the periodic-hill flow, and largely identical to that used by Meyer et al. (2013) for wall-modelled LES of a low-speed aerospace application.
A parameter study was performed on the effects of the span-wise extent of the domain, the mesh resolution, and the solution averaging time (Gómez, 2018). A good compromise of accuracy and cost was achieved for a spanwise extent of = 4.5ℎ (with ℎ = hill height), about 1.2 million cells (giving + < 20 everywhere), and averaging over ≃ 55 flow-through times started after an initial transient of ≃ 35 flow-through times. With these parameters, converged statistics are obtained in less than 2 days on a dual Intel ( ) Xeon ( ) CPU E5-2670v2. Figure 1 shows our LES results for the mean velocity profiles at several streamwise locations for the standard periodic-hill geometry and the well-resolved LES data of Breuer et al. (2009) (13.1 million cells) and the reference experiments of Rapp (2009). We observe that our INCA LES results are in very good agreement with the reference data.
Also shown in Figure 1 are the RANS results for the baseline -SST model described in Section 2.1. The RANS simulation significantly over-predicts the length of the separated flow region, which distorts the flow and leads to large errors everywhere in the domain. Having flow-features very dissimilar to the high-fidelity LES, the baseline RANS is therefore likely inappropriate as a low-fidelity model in a bi-fidelity optimization or surrogate modelling procedure.

Enhancing RANS with data-driven turbulence modelling
In this section we briefly describe the SpaRTA procedure to generate improved RANS models from LES datasee Schmelzer et al. (2020) for a detailed discussion. The procedure has two main stages: 1. We use the -corrective-frozen-RANS approach (Schmelzer et al., 2020) to solve for corrective fields for the anisotropy tensor, and turbulence production. These fields, when directly injected into the RANS solver, lead the solver to reproduce the LES mean-flow.

Deterministic symbolic regression for corrective field modelling
Having identified corrective fields, to make predictions it is necessary to model these fields in terms of known (mean-flow) quantities. Following Pope (1975), we assume that the non-dimensional strain-rate and rotation-rate tensors̃ ∶= 1 ,Ω ∶= 1 2 − are sufficient to describe the corrective fields, leading to Pope's well-known basis-tensor series. The Cayley-Hamilton theorem dictates that the most general form of the function Δ (̃ ,Ω ) is (under the assumptions of analyticity and Galilean invariance): where ( ) are ten basis tensors, are the five invariants of̃ andΩ, and (⋅) are arbitrary scalar functions of these invariants. In this paper, we only consider the first three base tensors, and first two invariants, which are given by: To modelR we assume that it takes the form for some tensor-field , modelled similarly to Δ by (7). It remains to fit a function of the form (7) to the LES data, for which we use symbolic regression. For details refer to Schmelzer et al. (2020). In short, we form a large library of candidate basis functions, consisting of powers, and products of the invariants . We then regress the data -consisting of ( , ( ) , Δ ) triples at each spatial mesh point -against this library with sparsity-promoting 1 regularization. This selects a small set of basis functions that represent the data well. Concretely, let the regressing function be: where are coefficients of the basis functions in the library (⋅) for each basis tensor, and (⋅) is the implied operator, which is linear in . Then we solve the elastic-net regularized regression problem: where the first norm is over all mesh points of all training cases, and once again ⋆ denotes quantities evaluated from the LES data. We solve for a number of different values of and , to obtain a number of candidate sparsity patterns, and we discard the values of the nonzero coefficients . For each of the sparsity patterns, we perform ridge regression, using only the basis functions identified in the sparse regression. Let be a vector of non-zero coefficients identified by (8), then we solve for each sparsity pattern. The final model is selected based on a compromise between sparsity, goodness-of-fit, and ideally performance in cross-validation if data is available (Schmelzer et al., 2020). This two step procedure allows us to independently control the level of sparsity (with , ), and the magnitude of the resulting coefficients (with ), which can otherwise become very large. Applying this procedure to our LES data for the standard P-H case ( = 1), the resulting algebraic models for Δ andR are: Figure 3 shows the optimal = 2 Δ andR identified for the P-H case using the -frozen approach of the previous section. The predictions of the model (10) are also shown. The agreement is generally good, with the notable exception of the largeR near the wall, which is not reproduced by (10). This failure can be attributed to the nondimensionalization of , Ω by the specific dissipation rate , which becomes large close to the wall, pulling all basistensors ( ) down to zero. This effect can be overcome by using coefficients that go to infinity near the wall, as in (Kaandorp and Dwight, 2020). However, we consider this unnecessary, as corrections close to the wall, especially in Δ , have very little effect on the mean-flow. Implementing Δ andR in the RANS solver, we verify that the mean-flow of the training case is reproduced, and we make predictions for = 0.25 and = 4.0 and compare with our LES reference, see Figure 2. In all cases the mean-flow is significantly improved over the baseline -SST, but do not always reach the accuracy of the directly injected corrective fields shown ("frozen RANS" in Figure 2) which represent a best case scenario.

Bi-fidelity optimization loop with custom RANS as low-fidelity
We aim to efficiently solve the discretized PDE-constrained minimization problem min ∈Ψ ( ; ) subject to LES ( ; ) = 0, where are geometric design-variables in the design space Ψ, is the full flow-state, is the cost-function and LES is the (highest-fidelity) discretized PDE operator (which depends on via boundary-conditions). We assume that no significant modelling is required to evaluate given , which is the case for most cost-functions of interest. We also have available a low-fidelity discretized operator RANS [ ]( ; ), which depends on some closure model derived using the procedure of the previous section.
2. Identify the model-form error terms Δ andR using the -corrective-frozen-RANS approach of Section 3.1 and the LES data from Step 1. Subsequently find the algebraic models 0 Δ , 0 using Section 3.2, to obtain Note that in this work the low-fidelity model depends on the result of the high-fidelity (and the locations at which this is evaluated), in contrast to standard multi-fidelity models. This presents a challenge from the standpoint of statistical modelling of the bi-fidelity system -for example in cokriging, where a prior correlation relationship must be specified between the fidelities. This will be the subject of future work. We bypass the issue here by using hierarchical Kriging. The optimization algorithm is conducted using an in-house tool "SurroOpt" (Han, 2016).

Variable-fidelity surrogate modelling: Hierarchical Kriging
The hierarchical Kriging model ) (HK) is constructed in a recursive way. First, an ordinary Kriging model̂ Low is built for the low-fidelity objective function. Then a surrogate for the high-fidelity cost function is built using the scaled low-fidelity Kriging model̂ Low as the model trend. In this way the high-fidelity responses are treated as realizations of a random process: where is a stationary Gaussian process with zero mean and a covariance Cov[ ( ), ( * )] = 2 ( , * ), where 2 is the process variance, and (⋅, ⋅) and the correlation function. In this paper, we use a cubic spline function , which introduces hyper-parameters .
Similar to the standard single-fidelity Kriging model, the corresponding predictor and mean-squared error (MSE) of Kriging for the high-fidelity cost function arê Here, ∈ ℝ is the vector of high-fidelity observations at sample sites is the correlation matrix between the observed high-fidelity samples and is the correlation vector between the observed samples and the predicting point. Therefore, the Kriging model for high fidelity cost function can be regarded as a sum of scaled lower fidelity Kriging predictor and the discrepancy between the scaled lower-fidelity function and higher-fidelity function.
Model fitting of the HK model uses the maximum-likelihood estimate: given which an explicit expression for is: It remains to specify in what manner new high-fidelity samples are chosen. We use the maximum expected improvement (EI) principle (Jones et al., 1998;Zhang et al., 2018). Define the improvement with respect to a current best value min as (formally): Since the prediction of the HK model at any untried site is the normal distribution  (̂ ( ), ( ) 2 ), we choose the value of that maximizes the expected improvement ( ). This has a closed form expression: where (⋅) and Φ(⋅) are respectively the unit normal density and distribution functions.

Periodic hill optimization: Cost function and design-space exploration
The baseline periodic-hill geometry (hill width = 1) is the widely known ERCOFTAC Testcase 9.2 with the geometry definition of Mellen et al. (2000). We choose an objective function depending on the total drag on the hills and averaged turbulence kinetic energy over the whole flow fields at a specified constant mass-flow. To demonstrate optimization in this proof-of-concept study, we choose a cost function has two local optima within the design-space, ∈ Ψ = [0. 25,4]. The chosen function is: wherē , is the averaged turbulent kinetic energy over the flowfield, and , the volume forcing term in (1), acts as a proxy for drag. The cost function is normalized by the LES result of the baseline geometry, 0 ∶=̄ ⋆ − 2.5 ⋆ . The geometries of periodic-hills at the lower-and upper-boundary of the design-space are shown in Figure 5. Before the optimization, to evaluate the accuracy of the customized RANS model, we use the model (10) trained at = 1 from Section 3.2 to predict the flow at other values of , and compare with our LES ground-truth. In Figure 6, we see that the mean-velocity prediction is significantly improved compared with the baseline RANS simulation, across the majority of the design-space. Figure 7 shows that the prediction of the averaged turbulence kinetic energȳ and proxy-drag are both significantly improved -a consequence of more accurate Reynolds stresses and mean flow fields for all test geometries. The custom RANS performs best when ≥ 1, for which the flow shows a medium-sized separation bubble. For the steeper geometries < 1, the larger separation bubbles are slightly under-estimated by the custom closure model. In Figure 8, we observe that the custom RANS matches the main trend of LES very well. In terms of the objective function , the optimization problem has two local minimum points near = 0.75 and = 2.0, whereas the objective estimated by the custom RANS trained using the baseline data does only show the first one. With additional hi-fidelity LES samples, the custom RANS model will become more accurate during the optimization. The baseline RANS model, however, misses the trend and significantly under-predicts drag and Reynolds stress.

Optimization convergence
The optimization starts, as shown in Figure 9(a), with one LES as the high-fidelity sample (red solid circle) and 26 custom RANS simulations as the low-fidelity samples (red hollow circle). The initial HK response surface passes through the LES, and since there is only a single high-fidelity sample, the shape is dictated entirely by the custom RANS, and the scaling parameter in HK model is found to be 0.78. The expected improvement criterion selects a second LES sample at ≃ 0.93, Figure 9(a). LES is performed at this point, and a custom closure is built, which is found by Δ ≃ Δ = 1.2 (2) ; R ≃ = 0.45 (1) . Using this new closure model in RANS simulation, the low-fidelity samples near newly added sample are reevaluated. With the updated low-fidelity samples and two LES high-fidelity samples, the HK model is rebuilt, shown in Figure 9(b). Again the closure passes through the LES samples exactly, and the hierarchical kriging is getting more close to the true cost function in the region near the optimum. The scaling parameter in the HK model is found to be 0.93, closer to 1 because of the better match between low-and high-fidelity objective functions. A third iteration is shown in Figure 9(c). The convergence of the procedure in terms of is shown in Figure 10 (red solid line). As a comparison, we also performed a single-fidelity optimization using only the high-fidelity LES model, and a bi-fidelity optimization using the baseline RANS as the low-fidelity model. The convergence histories in terms of are shown in Figure 10 (black dash-dotted line and blue dashed line, respectively), where in each case we show the objective-function value of the most recently obtained sample. The convergence history of bi-fidelity optimization using the baseline RANS model as the lo-fi model, shows no improvement over the single-fidelity optimization -and is in fact significantly more expensive due to the large number of RANS solves needed. This is a typical example bi-fidelity methods with poor lo-fi to hi-fi correction resulting in no speed-up compared to hi-fi only. Consequently, more than 10 LES samples are required to locate the global optimum.
Finally, the optimum solution is found at = 0.75944 with the cost function value around −1.68. In summary, the optimization has converged to within the accuracy of the high-fidelity model with only 2 LES samples, with one further sample used to verify the result.

Conclusions
In this work, a novel bi-fidelity fluid-dynamic shape optimization methodology is proposed, in which LES is used as the high-fidelity simulation tool and an automatically customized RANS turbulence closure as the low-fidelity model. The full-field LES data obtained from the high-fidelity samples are used to train a RANS model for the flow being optimized. The SpaRTA method is used as a robust and effective way of tailoring RANS models with simple algebraic augmentation. The customized low-fidelity RANS model is successively improved as more LES data becomes available within the design space. The two fidelities are combined in a hierarchical Kriging surrogate, and the EGO procedure is used to add new high-fidelity samples.
Given the technical complexity of this procedure, we have demonstrated it on a proof-of-concept shape optimization of a periodic-hill geometry with varying hill-width. In this case, the method converges to the true optimum (within the tolerance of the LES simulations) in two iterations. Shape optimization for flow-problems where the details of turbulence are important, including separated flows, junction flows, and transitional flows, may become computationally feasible with this scheme. Future work will focus initially on application to a three-dimensional optimization problem of practical engineering interest.