Multi-fidelity Surrogate Modelling of Wall Mounted Cubes

This paper focuses on the application of multi-fidelity surrogate modelling to characteristics of a flow as it changes with a parameter. This provides insight into the potential of combining multi-fidelity modelling approaches with varying fidelities of computational fluid dynamics methods to a parameter space exploration. A limited number of trusted high-fidelity large eddy simulation data points, in combination with an extended study using lower-fidelity Reynolds averaged Navier–Stokes modelling is used as the input for the surrogate model. Multi-fidelity surrogate models are implemented to bridge the low-fidelity and high-fidelity models providing an improved surrogate model over using a single fidelity alone. The flow around tandem wall mounted cubes at varying inlet yaw angle is used as an aerodynamic test case for this methodology. Results presented show that the multi-fidelity surrogate modelling provides a significant improvement over single fidelity modelling for the prediction of global flow properties. This methodology is then extended to combine multiple local flow features into the multi-fidelity model to build up fuller descriptions of the flow at angles not included in the training data for the model. The results of this are presented for both one-dimensional line plots at a range of locations along the center line of the flow and for two-dimensional slices of the velocity field. The multi-fidelity surrogate model produces results at locations in the parameter space away from the high fidelity training data that match closely to large eddy simulation results.


Introduction
In many instances of engineering for aerodynamic applications it is desirable to investigate how changing a certain parameter of the geometry or the flow will affect a quantity of interest (QoI). Examples of this include measurements of the lift and drag coefficients for a vehicle as the flow angle is varied such as the experimental studies of Bello-Millán et al. (2016) and Howell et al. (2002) on automotive applications. The geometry of a design may also be varied as a parameter such as in the numerical study on parameterised rib shaping by Pourfattah et al. (2017) or the geometry of a wingtip design as by Büscher et al. (2006), where reduced order modelling was used to build up a database over a range of designs. Designs that can be effectively parameterised such as aerofoils (Sobieczky 1999) are good candidates for such approaches.
To investigate this, it is typical to conduct an experimental or computational campaign (or a combination of the two (Deeb 2021)), varying one parameter between each individual run and recording the resulting QoI. From this data a surrogate model can be constructed to interpolate between the available results (Yondo et al. 2019;Sun and Wang 2019). This can then be used in an optimisation procedure to modify the design of a device as done by Yamazaki and Mavriplis (2011) and Tao and Sun (2019). To build a satisfactory representation of how the QoI varies with the chosen parameter, the number of runs conducted must be high enough to capture the changes in the QoI in parameter space. However, often the methods that provide the most accurate results, such as experimental runs in a wind tunnel or conducting turbulence-resolving simulations such as large eddy simulations (LES), are also the most costly to obtain. This constraint limits the number of runs that may be performed and therefore limits the resolution in parameter space that can reasonably be achieved. In contrast, other methods with a computational cost low enough to perform runs that cover the parameter space at a high resolution are limited in their accuracy by the modelling assumptions used.
Multi-fidelity (MF) methods aim to supplement the accurate high-fidelity (HF) results at sparse locations in parameter space, with potentially less accurate but less expensive low-fidelity (LF) results that cover the parameter space with high resolution. A MF model establishes a relationship between the LF and HF data to improve the prediction across the parameter space (Toal 2015). The review of Fernández et al. (2016) describes a commonly used correlation between the LF and HF data that was introduced by Kennedy and O'Hagan (2000) and is the basis of the commonly used co-kriging method (Le Gratiet and Garnier 2014). To model more complex relationships between the LF and HF (Perdikaris et al. 2017) proposed a generalized autoregressive scheme that introduces a non-linear relationship between the LF and HF components of the model. This work applies the MF method of Perdikaris et al. (2017) to generate surrogate models for an aerodynamic application, enabling the interrogation of the flow at unseen locations. LF data is collected from a series of Reynolds averaged Navier-Stokes (RANS) calculations and HF data is collected from LES calculations of the flow around a series of tandem wall mounted cubes. Creating a surrogate model for a single global parameter is useful for a design process, however understanding details of the flow can further aid in the development of aerodynamic designs. To this aim the MF method is extended to create a surrogate model that can predict more complex descriptions of the flow for both onedimensional lines and two-dimensional surfaces.
The remainder of this paper begins with Sect. 2 describing the MF framework and the regression models used within it, followed by a description of the tandem cube test case in Sect. 3. Section 4 presents an overview of the CFD results of the flow around the cubes followed by the application of the MF regression (MFR) models to the CFD results. Conclusions of the paper and possible avenues for further exploration are outlined in Sect. 5.

Multi-fidelity Regression
When constructing a surrogate model using MF modelling the goal is to not only establish a link between the input parameter and the QoI, but to also establish and use the relationship between the LF and HF data (Peherstorfer et al. 2018). A commonly used correlation between the LF and HF data can be expressed as (Fernández et al. 2016): where Y L is the LF and Y H is the HF data, whilst (X ) and (X ) are the multiplicative and additive correlation surrogate components. However, in many of the applications relevant to aerodynamic flows the relationship between the two fidelities is unlikely to be linear as is represented in Eq. 1. For such cases (Meng and Karniadakis 2020) suggested a generalized autoregressive scheme, which is expressed as: where F(⋅) is an undefined linear or nonlinear function, mapping the LF to the HF. This can then be combined with the additive correlation to form: In the current work two types of regression methods have been used as the base for the MF regression. Firstly Gaussian process regression (GPR) as in Cutajar et al. (2019), and secondly a multi-layer perceptron (MLP) neural network as in Meng and Karniadakis (2020). Figure 1 shows the architecture of an example of a MF-MLP.
The prediction from the LF regression at the known HF points in parameter space ( Y L (X H,i ) ) form additional input features to the HF part of the model. It is this inclusion to the model that allows the bridging of the two fidelities. A summary of the procedure used for the multi-fidelity regression is presented in Algorithm 1 In Fig. 1, L and H are the unknown weights and biases of the neural networks, which are learnt by minimising the respective loss functions defined as (Simon 1999): where Y * H and Y * L are the desired target values, Y H and Y L are the values produced by the perceptrons, and N H and N L are sample sizes, for the HF and LF respectively. L2 regularization has been used to help limit the occurrence of over-fitting with a value of = 0.0001 . The loss function is optimized using the L-BFGS solver from the family of quasi-Newton methods due to its fast convergence on smaller datasets (Liu and Nocedal 1989). To find an optimal architecture for the MLPs a random search is conducted over 60 possible hyperparameter configurations with leave one out cross validation performed for each. The configuration with the lowest mean squared error was chosen as the optimal solution and results are presented using those hyperparameters.
The same MF framework is used for the MF-GPR but with F L and F H determined using GPR rather than an MLP. For GPR the collection of possible functions mapping the features to the targets will have a Gaussian distribution, such that (Williams and Rasmussen 2006): where k is the kernel providing a measure of the statistical correlation between two points of the input space (X, X � ) . The covariance matrix is built so that K ij = k(X i , X j ; ) , with being the hyper-parameters.
For both MF-MLP and MF-GPR models embedding theory is used following the method described by Lee et al. (2019) and Meng and Karniadakis (2020) to better capture the link between the LF and HF data. This has been shown to provide a particular (4) benefit when there is a phase shift between the LF and HF data which is likely to occur in aerodynamic results. This modifies 3 to: Where is a shift in the parameter space and m is the number of shifts used. In this work values of = 0.02 and m = 2 are used. The MF framework presented is initially used for a single target, however to build up a fuller representation of the flow multiple targets should be considered. One approach would be to generate separate models for each target parameter of the flow, however in the case of a velocity field for example the multiple targets will be correlated with each other to varying degrees. This means that using separate models for each target may result in a reduction of information (Liu et al. 2018). To utilise the correlations between the outputs of the model multi-target regression is implemented (Borchani et al. 2015;Waegeman et al. 2019;Xu et al. 2019). When applied to the MF framework the relationship between the LF and HF represented in Eq. 3 components is also enhanced leading to: Where each of the multiple LF targets ( Y L,i ) are used as additional input features to the HF part of the MF model. This means that each of the HF target values ( Y H,i ) is correlated to the other HF targets as well as the LF targets.
In the current work the scikit-learn toolkit (Pedregosa et al. 2011) and its python API are used to build the MF framework described above for both MF-GPR and MF-MLP.

Array of Tandem Cubes Test Case
In this work we investigate the application of the MF surrogate modelling to computational fluid dynamics (CFD) results of the flow around an infinite array of tandem mounted cubes on a flat plate at a Reynolds number Re = 22, 000 . Previous experimental and numerical studies have been performed on this configuration. The experimental investigations of Havel (2000, 2004) showed that, with cubes separated by a distance of 4H in the streamwise direction, the flow separates at the leading edge of the upstream cube and reattaches at the wall before reaching the second cube, allowing a second horseshoe vortex system to be formed upstream of the second cube. The numerical study of Paik et al. (2009) compared the results of computations using RANS and various formulations of Detached Eddy Simulations (DES). They found that the RANS failed to capture the reattachment and subsequent second horseshoe vortex at the upstream junction of the second cube. Other key features of the flow, such as the recirculating flow on the top surface of the downstream cube and the turbulent statistics, were poorly predicted using RANS. The spacing between the cubes in the tested configuration is L = 4H , where H is the cube height, with the full set-up shown in Fig. 2.
This data set provides a useful basis for this work as it contains a number of features relevant to more complex geometries. The flow is fully turbulent and exhibits a complex semideterministic time-signal in the wake. This signal is a combination of coherent flow structures and stochastic motion, and as such represents a simple example of the issue at hand. The flow is highly three-dimensional in nature and is beyond the reach of LF models such as standard engineering turbulence models based on the RANS equations. The previous study of Paik et al. (2009) found a standard RANS approach significantly over-predicts the size of the flow recirculation region between the cubes and only with LES was the correct flow predicted, when compared to the reference experimental data. This makes the case a good candidate for MF modelling using RANS as a source for the LF data and LES providing the HF data. The coupling between the cubes is expected to change in a complex manner with variations in the parameter space, making the generation of a surrogate non-trivial.
Multiple configurations of the flow were run to build up a database of the flow as parameters were modified. Configurations of the case were run with the angle of the inlet flow direction, , varied as shown in Fig. 3.
The multi-fidelity method presented in Sect. 2 is applicable to any data coming from multiple sources at different fidelities and costs. This might include (but is not limited to) direct numerical simulations (DNS) or LES providing HF data, simple intrusive experimental techniques such a hot wire anemometer at the LF level with more advanced techniques such as particle image velocimetry providing HF data. LF levels might also include RANS modelling, lower dimensional approaches, or analytical methods with assumptions.
In this study LES was used to provide the HF data as it provides accurate data without the increased cost of DNS or long lead time of setting up a wind tunnel experiment. The LES data sets are acquired using the filtered incompressible Navier Stokes equations, given by: where the triangular brackets denote a filtering operation, ⟨u i ⟩ is the resolved velocity field, ⟨ p⟩ is the resolved pressure field, is the fluid density and is the kinematic viscosity  of the fluid. ⟨ u i u j ⟩ are the turbulent stress terms accounting for both the resolved and sub grid scales (SGS). In the current work the SGS turbulent stresses are modelled using the k-equation SGS model (Yoshizawa 1986). Equations are solved using the PISO algorithm with second order schemes employed in space and time. The mesh size of the LES case is 32 million cells. With the LES data spanning results from 0 • to 30 • in 5 • increments, there are a total of 7 cases using LES, 4 of which are used in constructing the MF models and the other 3 are withheld for testing.
The LF data for this study was collected using incompressible RANS calculations solving: where, in the final term, u i u j are the turbulent Reynolds stresses modelled here using the k − SST model (Menter et al. 2003). A total of 21 RANS cases were run from 0 • to 40 • at 2 • intervals, each on a mesh with a cell count of N c ≈ 2 × 10 6 .
Both RANS and LES calculations were conducted with the finite volume code Open-FOAM-v1812 (Weller et al. 1998). For all calculations the inlet was defined by a fixed constant velocity with a zero pressure gradient boundary. The outlet was defined with a zero gradient velocity and a fixed pressure of 0. All lower walls and surfaces of the cubes have a no-slip velocity boundary condition applied. The top surface has a freestream boundary condition and the left and right sides of the domain are periodic. The inlet velocity is defined based on the yaw angle such that U x = U 0 cos( ) , U y = 0 and U z = U 0 sin( ).
Structured block meshes were used for both the RANS and LES calculations with the RANS mesh shown in Fig. 4. This structure allowed for localised refinement at the wall and cube surfaces to achieve a y+ ≈ 1 . For the RANS cases this resulted in a total mesh count N c ≈ 3 × 10 6 . For the LES mesh a finer grid was used in the directions tangential to the walls resulting in a total mesh count of N c ≈ 20 × 10 6 . In all cases, the computational cost of training and inference of the ML models was negligible relative to the cost of running the CFD.

CFD Results
The result obtained from running a RANS and an LES around the base configuration of the cubes has been compared with experimental data. A visualisation of the velocity field and the recirculation regions along the centerline of the cubes in the z-normal direction for each method is shown in Fig. 5.
It can be clearly seen in Fig. 5 that the LES matches closely to the experimental result, with the separated flow between the two cubes reattaching at a similar location before the second cube, allowing for the formation of a second horseshoe vortex at the base of the second cube as well as the first. In the RANS however, there is no reattachment between the cubes and the flow remains separated for the entire span between the cubes. This failure of this RANS model to capture the reattachment in this case is due to the turbulence model under-predicting the level of turbulence as the flow separates around the first cube. This under-prediction means that there is less energy transferred from the high energy freestream flow into the separation region which enables the reattachment seen in the LES and experimental results. Whilst in principle a different RANS model may provide a better prediction of this flow an investigation in to this was not conducted as the difference in performance between the RANS and LES provides an opportunity for testing the MFM framework. As there is currently no RANS model that can be universally applied to match the performance of scale resolving methods in all scenarios the development of MF surrogates will be beneficial in many cases. Figure 6 shows contours of the velocity magnitude of the flow around the cubes at a height of y = 0.1H calculated from the LES cases at = 0 • , 5 • , 10 • , 15 • , 20 • , 25 • , 30 • . Both the instantaneous and the mean velocity contours are shown. The instantaneous contours show a combination of the coherent flow structures and stochastic motion from the resolved turbulence. The horseshoe vortices can be seen clearly in front of the first cube but in front of the second cube include a much higher degree of turbulent motion. Contour slices of the mean velocity magnitude at y = 0.1H are shown in Fig. 7 for each angle of attack tested using RANS. This provides an overview of the flow behaviour at varying angles of attack.
The coupling between the cubes here varies in a non-linear fashion for each of the different configurations in yaw. One example of this is the horse-shoe vortex ahead of the second cube as can be seen by the slices closer to the wall in Fig. 7. Due to the increased separation exhibited in the RANS calculations the horse-shoe vortex is not present ahead of the second cube for low angles of attack ( 0 • < ≤ 6 • ), however as the wake of the first cube moves around as the yaw increases ( > 6 • ), the second horseshoe vortex begins to form. This effect is very different from in the HF LES cases where the second horseshoe vortex is present at all angles of attack. This leads to a complex non-linear relationship between the LF RANS data and the HF LES data that will need to be predicted by the MF regression model.

Single Parameters
The MF surrogate modelling strategy has initially been tested on single parameters of the flow using the RANS results to provide the LF and LES results to provide the HF data. The coefficient of drag on the second cube, Cd 2 , and a probe at a single point in the flow of the Reynolds or time averaged velocity, U x (6, 0.6, 0) , are calculated from the RANS and LES. This probe location was chosen as it exhibits a reasonable amount of non-trivial variation as the inlet angle is varied. Both of these parameters have been  Fig. 8 use GPR whilst the results in Fig. 9 are generated using MLPs. Single-fidelity surrogate models based on the LF and HF data respectively are plotted alongside the MF surrogate model to provide a comparison. Selected HF data is withheld from all modelling to provide a test to see how well the surrogate models are performing at unseen yaw angles. The plots on the right of Figs. 8 and 9 show the relationship between the LF surrogate model and the HF output of the MF surrogate model. This shows the relationship between the LF and HF components predicted by the MF model and how this changes the surrogate model.
For the GPR the kernel used for the LF part of the model is the Matérn kernel, whilst the HF part uses a combination of the Matérn kernel, k M , and the linear kernel, k L , such that k hf = k M k L k M k L + 1 , following a similar construction to that of Cutajar et al (.2019). When using the MLP, the number of layers and number of neurons in each hidden layer were determined from the hyperparameter optimisation.
The results in Fig. 8 show that both Cd 2 and the probe of the velocity U x (6, 0.6, 0) are better predicted using the MF-GPR than either the LF or HF GPR model. This is particularly evident when comparing the surrogate models to the untrained HF data point at = 5 • . The MF-GPR model provides a much better prediction for this data point and the surrounding behaviour, than the HF-GPR model alone. Although it does not fully capture the Cd 2 at = 5 • , with it lying just inside the − bound, it does provide a downward trend towards this data point. The MF-GPR model for the probe of the velocity U x (6, 0.6, 0) matches closely to the unseen HF data point at = 5 • . This may be easier for the model to predict due to the more linear relationship between the LF and HF data. In both cases the standard deviation calculated from the MF-GPR model is largest either side of = 10 • suggesting that the uncertainty of the model is greatest here. Conversely, the standard deviation is reduced where the HF and LF data are similar suggesting an increased confidence in the MF model.
The MF-MLP model provides a significant improvement over the HF-MLP model alone for both Cd 2 and the velocity probe U x (6, 0.6, 0) shown in Fig. 9. The drop in Cd 2 at = 5 • is also better captured by the MF-MLP model than the MF-GPR model. This is captured even with the HF model only being trained on the four HF samples at 10 • increments in the parameter space which do not capture this drop. It is the inclusion of Y L to the HF part of the model that allows for this feature to be modelled. Further improvements would be possible through increasing the number of training HF samples, especially in the range 0 • < < 10 • . Additionally only including 4 HF data samples results in a 4-fold cross validation used for the hyperparameter optimisation. Finding a better set of hyperparameters through increasing the number of folds possible in the cross validation may also bring improvements.
The plot on the right of Fig. 9 shows that the relationship between the LF and HF data predicted by the MF-MLP model is complex and non-linear. It is this learnt relationship that is allowing for trends unseen in the HF training data can be captured by the MF models. The difference in the relationship between LF and HF is providing the improvement over the MF-GPR for this case. The local flow parameter U x (6, 0.6, 0) is also well predicted by the MF-MLP model, matching well to the test samples at = 5 • , 15 • , 25 • . The relationship between the LF and HF predicted for this parameter is much closer to a linear relationship which may indicate why both the MF-GPR and MF-MLP models perform slightly better in predicting this relationship with the limited data available.
at flow angles past the final HF data point at = 30 • any prediction is an extrapolation. In such instances the model should not be relied upon. Within the MF framework it may be the case that the LF data extends past the final HF data point. This may result in an improvement in the ability of the model to extrapolate. In Fig. 8 the GPR confidence  interval for the MFR is smaller than for the HF-GPR suggesting that there is an improvement to the model in this region, however they still increase considerably past the final HF data point at = 30 • .

Velocity Profiles
In order to build up a fuller description of the flow across the parameter space, a greater number of features and targets can be included in the MF regression. Initially this is done by using velocity profiles at locations along the centerline of the cubes as the regression targets rather than a single parameter. By inputting multiple features to the MF model, more complex relationships can be learnt by exploiting the correlation between the outputs, which provides a better prediction than modeling them individually (Liu et al. 2018).
Within the MF framework this should also allow for an improvement in the predicted relationship between the LF and HF data. Profiles of the mean streamwise velocity, at intervals of x∕H = 1 along the centerline of the domain ( z = 0 ), for the same locations in the parameter space as in Sect  The plots in Figs. 10 and 11 show profiles of velocity predicted by the MF models at = 5 • , 15 • , 25 • and are compared to the RANS and LES profiles at these angles, unseen by the MFR model. These locations in parameter space are the furthest from the samples provided for training, providing a good test of the performance of the MFR models. As in the single parameter case a leave one out cross validation is performed to optimise the hyperparameters for the MF-MLP models.
The results of the MF-GPR in Fig. 10 and MF-MLP in Fig. 11 trained on the velocity profiles both show a reconstruction of the velocity profiles that match more closely to the HF results of the LES than the LF RANS profiles at these locations at most of the locations. However, the predictions of the profiles at different locations along the x direction perform with varying levels of success. In particular, the profiles between the cubes at = 5 • show the largest difference between the MF modelling and the test LES data. However, the MF models here still provide a better prediction than the LF RANS calculation. Both MFR models also struggle to correctly predict the velocity deficit at the trailing edge of the second cube x∕H = 8.
Some profiles generated using the MF-MLP model in Fig. 11 appear noisy. In particular the profiles at locations x∕H = 4 , x∕H = 7 and x∕H = 9 for = 5 • , at x∕H = 5 for = 15 • and the profile at x∕H = 3 for = 25 • show significant noise. The universal approximation theorem suggest that the same relationship is possible with a MF-MLP as with the MF-GPR however the constraints of reasonable hyperparameters for the width and depth of the network may be limiting this for these profiles.

Velocity Slices
The method is further extended by including more features and targets into the training data. Rather than creating models for individual velocity profiles, a slice of the velocity at the y∕H = 0.1 plane is used to provide features and targets for each sample in parameter space. Results here are only presented using the MF-GPR method as this performed better in the previous cases with the current set-up without the need of the same level of hyperparameter optimisation. Velocity slices at y∕H = 0.1 are extracted from the results at the same HF and LF sample locations in parameter space as described in Sect. 4.2.1 and 4.2.2. As above the HF LES results are interpolated onto the LF RANS mesh to provide consistent targets. The mean streamwise velocity, at each of the 63140 cell center locations intersected by the slice, U(x∕H, y∕H = 0.1, z∕H) , provides the targets for the MF regression. Figures 12, 13 and 14 show contour plots of the inferred velocity in the x direction for the HF LES, LF RANS and the field generated from the MF-GPR model for each flow angle = 5 • , 15 • , 25 • respectively. To enable a comparison between the cases only the contour at U x ∕U 0 = 0.5 is shown. These show that the generated velocity slices constructed from the MF-GPR model match closely to the unseen slices of the velocity calculated using LES. This shows that the MFR model is identifying and reproducing the main features from the LES at other angles. One feature of the flow that is not as well captured by the MF-GPR model is the split in wake after the second cube. At = 15 • this feature is over-predicted, and at = 25 • it is under-predicted. However, the MFR model still provides a better representation of the flow in these regions than the LF result of the RANS.
For this case, with a large amount of data being included in the HF part of the model from the LF part, it may be beneficial to investigate the use of regression models such as convolutional neural networks. It is expected that this would have the most potential for improvement if used for the HF part of the model Y H and so a heterogeneous approach to the modelling methods may be explored.

Conclusions
The MF regression model approach developed by Perdikaris et al. (2017) has been applied to the turbulent flow around an array of tandem wall mounted cubes at varying inlet yaw angles. The method has been extended for multi-target regression to generate surrogate models for both profiles and slices of the flow field. LF and HF data is extracted from RANS and LES calculations respectively at various locations in the parameter space and the MFR methodology is used to create a surrogate model that exploits the non-linear link between the RANS and LES data to provide a better model than a single fidelity alone.
Tests of the MFR on single parameters of the flow show that the methodology provides a significant improvement over single fidelity modelling of these properties. The method has been extended beyond modelling a single parameter of the flow to both one-and twodimensional spatial fields. In both cases the method is shown to predict a velocity field at unseen locations with accuracy comparable to an LES calculation compared to an additional RANS calculation. Comparisons are drawn between the use of MF-GPR and MF-MLP for their performance on this case with MF-GPR found to more reliably produce smooth spatial fields through the multi-target provided to the model. We note that, although the profiles generated with the MF models match closely to the true velocity profiles of the LES, there is no guarantee that they will uphold physical constraints such as the no slip condition at the wall or continuity. In future work we will look to embed known physical constraints and boundary conditions into the models. Further work could also extend the parameter space to multiple dimensions with additional configurations of the flow, such as variation in the relative positioning of the cubes, as additional parameters. Different features of the flow may be included in the model for each of the fidelities. The extension to additional fidelities from lower order models or experimental results may be investigated, coupled with the inclusion of the uncertainty of each of the datasets.