Applying Bayesian Optimization with Gaussian Process Regression to Computational Fluid Dynamics Problems
Y. Morita, S. Rezaeiravesh, N. Tabatabaei, R. Vinuesa, K. Fukagata, P. Schlatter
AApplying Bayesian Optimization with Gaussian Process Regression toComputational Fluid Dynamics Problems
Y. Morita a,c, ∗ , S. Rezaeiravesh a,b, ∗ , N. Tabatabaei a,b , R. Vinuesa a,b , K. Fukagata c , P. Schlatter a,b a SimEx/FLOW, Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden b Swedish e-Science Research Centre (SeRC), Stockholm, Sweden c Department of Mechanical Engineering, Keio University, Yokohama, Japan
Abstract
Bayesian optimization (BO) based on Gaussian process regression (GPR) is applied to different CFD (com-putational fluid dynamics) problems which can be of practical relevance. The problems are i) shape opti-mization in a lid-driven cavity to minimize or maximize the energy dissipation, ii) shape optimization ofthe wall of a channel flow in order to obtain a desired pressure-gradient distribution along the edge of theturbulent boundary layer formed on the other wall, and finally, iii) optimization of the controlling parame-ters of a spoiler-ice model to attain the aerodynamic characteristics of the airfoil with an actual surface ice.The diversity of the optimization problems, independence of the optimization approach from any adjointinformation, the ease of employing different CFD solvers in the optimization loop, and more importantly,the relatively small number of the required flow simulations reveal the flexibility, efficiency, and versatilityof the BO-GPR approach in CFD applications. It is shown that to ensure finding the global optimum of thedesign parameters of the size up to 8, less than 90 executions of the CFD solvers are needed. Furthermore,it is observed that the number of flow simulations does not significantly increase with the number of designparameters. The associated computational cost of these simulations can be affordable for many optimizationcases with practical relevance.
Keywords:
Bayesian optimization, Gaussian process regression, Computational fluid dynamics, Turbulentboundary layers, Spoiler-ice model
1. Introduction
In computational fluid dynamics (CFD), the Navier–Stokes equations or variations of them are numeri-cally solved in order to simulate fluid flows. At most of the Reynolds numbers associated to environmentalflows and relevant to engineering applications, flows are turbulent, see e.g.
Ref. [14]. Various approacheshave been developed to numerically simulate turbulent flows, ranging from low-fidelity Reynolds-averagedNavier–Stokes (RANS) simulation to high-fidelity approaches including large eddy simulation (LES) and di-rect numerical simulation (DNS), see Ref. [51]. Moving from RANS toward DNS, the influence of modeling isreduced, the role of numerics becomes more dominant, and above all, the computational cost increases [56].In many applications where the fluid flows are involved, we need to optimize different parameters in order tomeet (within a margin) a certain set of objectives while satisfying a set of constraints, see Refs. [64, 17] andthe references therein. The aim of the present paper is to illustrate how Bayesian optimization (BO) [55, 20],which has been mostly utilized in the field of machine learning and data sciences, can be applied to differenttypes of optimization problems arising in CFD. Previous application of Bayesian optimization to CFD prob-lems is limited to a few studies, see Refs. [61, 39, 34]; therefore, the diversity of the applications consideredin the present study can reveal the flexibility of the method as well as shed light on its pros and cons. ∗ Principal Corresponding Author
Email addresses: [email protected] (Y. Morita), [email protected] (S. Rezaeiravesh), [email protected] (N. Tabatabaei), [email protected] (R. Vinuesa), [email protected] (K. Fukagata), [email protected] (P. Schlatter)
Preprint submitted to Elsevier January 26, 2021 a r X i v : . [ phy s i c s . f l u - dyn ] J a n o put the work in context, we shortly review the main approaches which have been employed forthe purpose of optimization in CFD. Such approaches are iterative and can, in general, be divided intogradient-based, gradient-free, and gradient-enhanced, see Refs. [41, 64]. The gradient-based methods arethe appropriate choice when the evaluation of the derivatives of the cost function with respect to the designparameters is computationally efficient. At each iteration, the first-order gradient is needed to determinethe searching direction and the second-order gradients ( i.e. Hessian matrix) or the approximations of themare required to obtain the optimal step size, see e.g.
Ref. [41]. Therefore, depending on the approach,a minimum order of smoothness for the objective function is necessary. The gradients can be computedby automatic-differentiation (finite differences), or more efficiently by solving forward or adjoint sensitivityequations. The use of the latter is recommended when the number of design parameters is more than thenumber of objectives, see [41, 26, 64]. As its main advantage, the gradient-based optimizations exhibit fastconvergence and can handle large numbers of design parameters. However, it is always likely to be trapped bylocal optima. Besides this, in the cases of unsteady Navier–Stokes equations the need for saving the forwardsolution fields which are required to solve the adjoint equations may lead to memory issues. To rectify theseissues, extra treatments are needed, see [3, 27]. Various types of application of gradient-based methodsfor optimization in CFD can be found in the literature, including shape optimization [37], optimization ofinitial conditions for fast transition to turbulence [29], and topology optimization in turbulent flows usingthe RANS approach [17].In gradient-free approaches, the simulator (here, the CFD solver) is treated as a blackbox and run toevaluate realizations of the response at the samples of the design parameters. This approach which goesback to Box and Draper [7] can be implemented using different methods, see the review by Larson et al.[30]. The simplest method is the grid search, where a set of manually selected samples are tested to find theoptimum. Evolutionary algorithms [2], which mimic the nature’s survival by iteratively simulating “selectionand mutation” of the design parameters, are also considered as an effective approach; see e.g.
Ref. [69] wherethe algorithm is applied to modify the RANS stress–strain relationship. However, the most frequently-usedmethod of this type in CFD is called the response-surface method (RSM), in which a metamodel or surrogateis constructed using a finite set of training data comprised of parameter samples and corresponding responses.The size of the training data set is determined as a balance between the required accuracy of the surrogateand the cost involved in running the simulator. In many CFD applications, Gaussian process regression(GPR), or the Kriging method [48, 25], has been employed to construct the surrogate. As detailed inSection 2, GPR is among the non-parametric surrogates, naturally allows for noisy (uncertain) data, andmore importantly, can estimate the uncertainties involved in their predictions. Examples of RSM-basedoptimization can be found, for instance, in Viana et al. [65]. Contrary to the gradient-based methods, RSMis suitable for finding the global optima, but it suffers from the curse of dimensionality. To reduce therequired number of samples for high-dimensional parameters yet constructing an accurate response surface,gradient information can be added to the training data, see the review by Laurent et al. [32]. In particular,the gradient-enhanced Kriging (GEK) has been extensively used in the last two decades for optimizationin CFD, see for instance, [11, 31]. The study by Laurenceau et al. [31] showed that at large number ofparameters (=45), GEK outperforms RSM-GPR, while at small number of parameters (=6), inclusion ofthe gradients makes no significant improvement in the computational performance of the optimization.In the present study, the Bayesian optimization based on Gaussian process regression, hereafter referredto as BO-GPR, is employed, see Refs. [20, 55]. This approach can be classified among the gradient-freeapproaches, although a gradient-enhanced version of it has also been proposed, see Ref. [71]. In BO, asequence of samples for the design parameters is taken which converges to the global optimum. Therefore,as a main difference with the RSM approach, the surrogate in BO is sequentially updated. That is, infact, the clear connection of BO to the Bayesian formalism, see e.g.
Ref. [57]. As detailed in Section 2,at each iteration, the decision about the next sample is taken based on combining the exploitation andexploration, where the latter directly takes into account the predictive uncertainty of the surrogate. Thisactive involvement of the uncertainties in the algorithm is another distinguishing characteristic of the BOapproach over RSM methods. In general, the predictive uncertainties are a combination of the uncertainCFD data and the constructed surrogate. The use of the BO is a timely choice considering the application ofthe uncertainty quantification (UQ) [57] in different aspects of CFD and numerical simulation of turbulence2ver the last decade, see e.g.
Refs. [4, 43, 72, 49, 50]. In fact, most of those studies are classified as UQforward problems, the outcomes of which can be nicely employed in BO, which is, in fact, an inverse UQproblem. Having the design parameters defined over an admissible space, the BO can be non-intrusivelylinked to any CFD solver. This provides a great flexibility noting that a CFD solver is not necessarilyequipped with adjoint solvers to compute gradients, and in any case, computing the adjoint sensitivitiescould be expensive.The application of BO in CFD has been very limited. Talnikar et al. [61] developed a parallel versionof BO to minimize the drag in a turbulent channel flow simulated by LES. Nabae et al. [39] used BO tooptimize the phase-speed of streamwise traveling wave-like wall deformation for drag reduction in a turbulentchannel flow simulated by DNS. More recently, BO was employed by Mahfoze et al. [34] to optimize theblowing amplitude and coverage in order to reduce the friction drag to achieve a certain net-power savingin a zero-pressure-gradient turbulent boundary layer (TBL) simulated by DNS. To study different aspectsof BO, a diverse set of optimization problems are considered in the present paper. The objectives, numberof design parameters, and CFD solver are different in these problems.This paper is structured as follows. In Section 2, a detailed review is given to the the theoretical aspectsof the adopted Bayesian optimization methodology. In Section 3, BO is applied to optimize the initialconditions of a set of coupled non-linear ordinary differential equations so that an energy norm at some latertime is maximized. As the second example, shape optimization of a laminar lid-driven cavity flow is studiedin Section 4, where the objective is to either minimize or maximize the energy dissipation. Another shapeoptimization which has practical applications, is considered in Section 5. The objective is to optimize theshape of a wall boundary of a channel, so that the streamwise pressure gradient of the turbulent boundarylayer (TBL) over the other wall matches a given profile. Section 6 is devoted to the optimization of aspoiler-ice airfoil profile. The parameters controlling the configuration of two spoilers which are located onthe pressure and suction sides of the airfoil are optimized so that the distribution of the pressure on theairfoil surface becomes the same as if a real ice was formed on the front part of the airfoil. In the end,concluding remarks are provided in Section 7.
2. Bayesian Optimization
Consider the set of design parameters q = { q , q , · · · , q d } which are allowed to vary over the admissiblespace Q ⊂ R d (note that Q is not the set of rational numbers in this paper). Running the simulator, i.e. theCFD code, realizations for the quantities of interest (QoIs) or responses are obtained, which in general canbe noisy. The aim of the optimization is to minimize (or maximize) an objective function r ( q ) which dependson the simulator outputs over space Q , and find optimal parameter q opt , q opt = arg min q ∈ Q r ( q ) . (1)Ideally, a functional dependency between the observed values r and parameters q could be establishedas [48, 25], r ( q ) = f ( q ) + ε , (2)where, ε denotes the observation noise, and without loss of generality it is assumed to be identically dis-tributed for all samples as N (0 , σ ). However, the complexities in the flow solver make acquiring a closed formfor f ( q ) infeasible. Alternatively, a surrogate ˜ f ( q ) can be constructed adopting a non-intrusive approach inwhich the simulator is treated as a blackbox. To this end, a finite set of training data D = { ( Q ( i ) , R ( i ) ) } ni =1 is employed, where Q and R are respectively samples of q and corresponding realizations of r . In general,the number of the training data, n , is limited due to the cost of running the simulator.The two main components of the Bayesian optimization, i.e. the choice of method for constructing ˜ f ( q ),and drawing the sequence of training samples, are now shortly reviewed. The Gaussian process regres-sion (GPR) [48, 25] is employed to build ˜ f ( q ). In this approach, ˜ f ( q ) is assumed to be random for which aprior distribution in the form of a Gaussian process is assumed. A Gaussian process GP is a multi-variateGaussian distribution defined as, ˜ f ( q ) ∼ GP ( m ( q ) , c ( q , q (cid:48) ; Θ ε , β )) , (3)3n which the mean and covariance are respectively given by, m ( q ) = E [ ˜ f ( q )] , (4) c ( q , q (cid:48) ; Θ ε , β ) = E [( ˜ f ( q ) − m ( q ))( ˜ f ( q (cid:48) ) − m ( q (cid:48) ))] . (5)In these definitions, E [ · ] denotes the expected value, Θ ε specify the parameters building the noise structure,and β are the hyperparameters in the kernel function which represents the covariance of ˜ f ( q ). Given thetraining data D along with the associated observation uncertainty, the posterior and posterior predictiveof ˜ f ( q ) (conditioned on D and Θ ε ) can be inferred. Simultaneously, the hyperparameters β are optimized.The posterior-predictive of ˜ f ( q ) at a set of test samples Q ∗ ∈ Q , where Q ∗ = { q ∗ ( i ) } n ∗ i =1 , has a multivariateGaussian distribution N ( m ( R ∗ |D , Θ ε , Q ∗ ) , v ( R ∗ |D , Θ ε , Q ∗ )), where, m ( R ∗ |D , Θ ε , Q ∗ ) = C ( Q ∗ , Q )( C ( Q , Q ) + C N ) − R T , (6) v ( R ∗ |D , Θ ε , Q ∗ ) = C ( Q ∗ , Q ∗ ) − C ( Q ∗ , Q )( C ( Q , Q ) + C N ) − C ( Q , Q ∗ ) . (7)Here, without loss of generality, the mean of ˜ f ( q ) in (3) is assumed to be zero, R = [ r (1) , r (2) , · · · , r ( n ) ], C ( Q , Q (cid:48) )is a n × n (cid:48) matrix where [ C ( Q , Q (cid:48) )] ij = c ( q ( i ) , q (cid:48) ( j ) ) and c ( · , · ) is a kernel function dependent on hyperpa-rameters β . Moreover, C N is the covariance matrix of the noise in the training data. In this derivation,the noise samples are assumed to be independent and identically distributed (iid) as ε ∼ N (0 , σ ). A moregeneral case in which the noise is allowed to be observation-dependent has been developed by Goldberg et al.[23], and recently applied for the purpose of uncertainty quantification in CFD by Rezaeiravesh et al. [50].The use of such a model for BO will be considered in future studies.Now, a short overview is given on how the sequence of samples for q can be drawn from Q so that theyconverge to q opt . At the k -th iteration of optimization the training data set is D k := { ( q ( i ) , r ( i ) ) } ki =1 . Thedecision about the next sample q ( k +1) is made through maximizing an acquisition function (ACF) α ( q ; D k ),see e.g. Ref. [20]. The most popular ACF is the expected improvement (EI) which was first suggested byMoˇckus [36] and then utilized in BO by Jones et al. [28]. The EI for the minimization problem is defined as,EI( q ) := E [ I ( q )] = E [max (cid:0) r ( q † ) − r ( q ) , (cid:1) ] , (8)where q † = arg min q ∈D k r ( q ), denotes the best estimate for optimum among the k observations (iterations).Note that the improvement I ( q ) is random since r ( q ) is so, see Eq. (3). In the cases of using a standardGPR for ˜ f ( q ), Jones et al. [28] derived a closed-form expression for the EI which reads as,EI( q ) = (cid:40)(cid:0) r ( q † ) − m ( q ) + ξ (cid:1) Φ( ζ ) + s ( q ) ϕ ( ζ ) , s ( q ) > , s ( q ) = 0 , (9)where s ( q ) = (cid:112) v ( q ), and m ( q ) and v ( q ) are the mean and standard deviation of the posterior predictive ofthe GPR at any q , as given by Eqs. (6) and (7), respectively. Further, Φ( ζ ) and ϕ ( ζ ) respectively representthe CDF (cumulative distribution function) and PDF (probability density function) of the standard Gaussiandistribution for ζ defined by, ζ = (cid:40)(cid:0) r ( q † ) − m ( q ) + ξ (cid:1) /s ( q ) , s ( q ) > , s ( q ) = 0 . (10)In Eq. (9), the first term in the summation specifies exploitation, the use of the best value so far, and thesecond term identifies exploration parts of the Q where the surrogate has highest uncertainty. The ad-hocparameter ξ can be used in practice as the controller of the exploitation in comparison with exploration.In the present study, the Bayesian optimization algorithm based on GPR is implemented using the Python libraries
GPy [24] and
GPyOpt [63]. The developed
Python codes are non-intrusively linked to theCFD solvers, Nek5000 [19] and OpenFOAM [70] through appropriate bash drivers. As a result, the wholeoptimization loop is fully automated. All optimizations started from a random sample for q taken from Q .4or constructing the kernel in the GPR surrogates, the Matern52 function is used, see e.g. [48, 25]. Tooptimize the GPR hyperparameters, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is utilized,see e.g. [41].Recognizing the convergence of the BO and hence imposing a stopping criterion for that can be, ingeneral, not so straightforward, see e.g. [9, 54]. In fact, this is a common characteristic of all gradient-freeoptimizations. In some cases, monitoring (cid:107) q ( k ) − q ( k − (cid:107) at the k -th iteration can help. But noting thatthe exploration guides the samples to different parts of the admissible space Q , obtaining a small differencebetween the successive samples, especially for large-dimensional parameters, is difficult to achieve in practice.We have observed that a more reliable measure is tracking the best value of the objective, i.e. r ( q † ) where q † = arg min q ∈D k r ( q ). If this measure remains unchanged after a sufficiently large n , then we can considerthe optimization to be converged. However, the value of n is not known a priori and it is basically limitedby the computational cost of the simulations and available computational budget. Despite this, examiningthe response surface constructed by the surrogate can help identifying parts of Q with insufficient samplesat which the predictive uncertainty of the surrogate is high. Guiding the sequence of samples toward thoseparts can be, for instance, done by adjusting ξ in Eq. (9). Despite this, the cases studied in Sections 5and 6 are not directly affected by the above discussion, since the optimum value for the objective is thedeviation of a computed QoI from a target value, and is therefore a priori known (albeit potentially notexactly achievable given the design space). In Section 4, where this is not the case, monitoring the bestvalue found when sufficiently many samples are drawn will help.
3. An introductory example
In many applications in fluid dynamics, in particular hydrodynamic stability, we are interested in findingthe optimum initial condition which leads to fastest growth of the unstable modes and hence fastest transitionfrom laminar state to turbulence. Kerswell et al. [29] have considered this problem through using adjoint-based optimization. As a “toy” problem to examine their approach, Kerswell et al. [29] introduced thefollowing dynamical system comprised of two ordinary differential equations (ODE),d x d t = (cid:20) − x + 10 x x (10 exp( − x / − x )( x − (cid:21) , (11)with the initial condition x (0) = x . The objective is to find the initial data x parameterized as x = a (cos πθ, sin πθ ) (with fixed a ) which maximizes (cid:107) x ( T ) (cid:107) /a , the normalized energy gain at time T . There-fore, the design parameter is q = θ which is allowed to vary over the admissible range Q = [0 , x ( t ) start from the same distance (cid:107) x (cid:107) = a . Depending on the valueof a , the trajectories for x ( t ) will end up at different attractors of the dynamical system at time T . In thegradient-based optimization adopted in [29], the adjoint equation of Eq. (11) is derived using the Lagrangemultiplier approach and is then solved backward from t = T to t = 0 in each iteration of the optimization.Here, we apply the BO algorithm explained in Section 2 to find the optimum θ for different values of a adopted in [29]. In each case the BO-GPR approach is capable of finding the global optimum with the sameaccuracy as in Ref. [29] without the need for any adjoint solution, see Fig. 1. This is a significant advantageof the BO-GPR method, noting that solving the same problem with the adjoint-based method could bemore involved and for some values of a the convergence of the algorithm would be a bit difficult to reachdue to numerical instabilities. On the other hand, the convergence in all cases shown in Fig. 1 is achievedwith less than 15 samples.
4. Shape optimization of a lid-driven cavity
We now move to problems involving the Navier–Stokes equations. Consider an incompressible two-dimensional cavity flow, where all walls are at rest except the upper wall (lid) which moves with the constantvelocity U in the positive horizontal direction, x . The aim is to optimize the shape of the left and rightwalls so that the energy dissipation is either minimized or maximized with the constraint of keeping the5 .0 0.5 1.0 1.5 2.00.000.010.020.030.04 || x ( T ) || / a SamplesTrue FunctionPredicted MeanOptimum || x ( T ) || / a SamplesTrue FunctionPredicted MeanOptimum (a) (b) || x ( T ) || / a SamplesTrue FunctionPredicted MeanOptimum || x ( T ) || / a SamplesTrue FunctionPredicted MeanOptimum (c) (d)
Figure 1: Gain (cid:107) x ( T ) (cid:107) /a plotted versus θ at T = 2 . a = 10 − , (b) a = 0 .
9, (c) a = 0 . a = 1 . a are the same as in Kerswell et al. [29]. The optimum θ maximizing the gain is shown by the solid marker,while other samples drawn during the Bayesian optimization are represented by hollow markers. The predicted mean by theconstructed GPR is close to the true gain obtained by solving (11) by fourth-order Runge-Kutta scheme given the optimum θ .In constructing the GPR, the σ in the observational noise in (2) was considered to be zero. The shaded areas specify the 95%confidence intervals built around the posterior predictive mean. total volume of the fluid confined in the cavity fixed and also retain the length and position of the upperand lower walls. Therefore, the objective is to find the wall shape Ω which either minimizes or maximizes, D = 2Re (cid:90) Ω S ij S ij d x , (12)keeping (cid:82) Ω d x = m fixed. In Eq. (12), S is the strain rate tensor, where S ij = ( ∂u i /∂x j + ∂u j /∂x i ) / i, j =1 , u i denotes the i -th component of the velocity vector. In order to apply the Bayesian optimizationof Section 2, we parameterize the geometry using a third-order polynomial, ¯ x = a + a ¯ y + a ¯ y + a ¯ y ,where ¯ x = x/l , ¯ y = y/h , l is the length of the upper and lower walls, h denotes the height of the sidewalls, and, x and y specify horizontal and vertical coordinates, respectively. Considering the origin of thecoordinates at ( x = 0 , y = 0), the constraint of keeping l fixed leads to, a ,L = 0 and a ,L = − ( a ,L + a ,L )for the left wall, and, a ,R = 1 and a ,R = − ( a ,R + a ,R ) for the right wall. As a result, four designparameters are left to optimize, q = ( a ,L , a ,L , a ,R , a ,R ). The admissible space of q is the tensor-productof Q a ,L = Q a ,R = [ − . , .
7] and Q a ,L = Q a ,R = [ − . , .
5] which are arbitrarily chosen to allow fordifferent wall shapes to occur. For any sample taken from Q , first the shape of side walls is determined.At this stage, the height of the walls, h , is adjusted to satisfy the constant-volume constraint. Then, themesh is generated using gmsh [22]. The simulations are performed using Nek5000, an open-source spectral-element-based flow solver developed by Fischer et al. [19]. In Nek5000, the weak form of the incompressibleNavier–Stokes equations are integrated in time on Gauss-Lobatto-Legendre (GLL) points through expressingthe velocity and pressure fields in terms of Lagrange interpolants of the Legendre polynomials, see the detailsin Deville et al. [16]. To avoid the discontinuity in the velocity at the top corners where the moving lid andthe side walls meet, the smoothing suggested in Ref. [42] is applied. In all simulations, 30 elements in both x and y directions are considered with 10 GLL points in each spatial directions per element. Simulations start6 igure 2: Contours of the velocity magnitude (normalized with the lid velocity U ) and streamlines of the cavity flow with(top) minimum and (bottom) maximum energy dissipation. Third-order polynomials are used to parameterize the side wallsof the cavity. The flow is steady with Re = 2000. at t = 0 from a uniform velocity field and continue until t = 50 l/U , which is chosen to ensure the steady-state value for dissipation (12) has been reached. Then, decision about next sample of the design parametersis taken. This algorithm is repeated until optimal design parameters are obtained. For this example we onlyconsider steady laminar flow, however an extension to turbulent flow is possible in a straight-forward way.For a cavity flow at Reynolds number Re = U l/ν = 2000 ( ν is the kinematic viscosity), Fig. 2 showsthe obtained shapes of the cavity with minimum and maximum energy dissipation. The obtained shapesare consistent with the flow physics, for instance, in case of the maximum dissipation two vortices adjacentto the lid are generated which prevent the fluid underneath to receive a high velocity. The obtained optimalshapes can be compared to results reported by Nakazawa [40] using an adjoint method for optimization,although the setups in the two studies are slightly different. In Ref. [40] the lid velocity was assumed tovary with x (it might be hard to interpret it physically) and the Reynolds number Re was higher than thecritical value of about 8000 [6] (although no treatment for turbulence modeling was mentioned). For thecase of minimizing the dissipation, the shape in Fig. 2(top) agrees well with [40]. But we found the shapein Fig. 2 (bottom) has a slightly higher dissipation compared to the maximum-dissipation case reported byNakazawa [40] in which both side walls were deformed towards the inside of the cavity. This can be anindication for the ability of the Bayesian optimization in exploring different possible geometries and findingthe global optimum without quickly getting stuck in a local optimum.Following the discussion in Section 2 regarding the convergence of the BO, the reported optimal shapesin Fig. 2 are based on the parameters which were found to be the best values within a sample size of n = 50.In fact, the parameters for minimizing and maximizing the dissipation were already found at the 20-thand 15-th iterations, respectively. However, the extra samples were taken to ensure the parameter space issufficiently explored. It is recalled that the definition of the convergence criterion in Bayesian optimizationis not straightforward, see [9, 54]. Therefore, in order to impose an appropriate stopping criterion, the wholehistory of the samples has to be considered. 7 . Shape optimization of a diverging channel Studying in a controlled fashion pressure-gradient (PG) turbulent boundary layers (TBLs) is importantsince they are present in a wide range of industrial applications, for instance, the flow around the curvedboundaries such as car bodies, airplanes, airfoils, etc. The so-called Clauser pressure-gradient parame-ter β [12, 13] is widely accepted as one of the most relevant non-dimensional parameters for PG TBLs [52],which is defined as, β = δ ∗ τ w d P e d x , (13)where δ ∗ , τ w and d P e / d x are, respectively, the displacement thickness, the magnitude of the wall-shear stressand the pressure gradient along the boundary-layer edge in the streamwise direction. Since the history ofthe PG β ( x ) has a significant impact on the characteristic of the TBLs, flows with a constant β in thestreamwise direction are very relevant for the study of PG TBLs [5, 18]. Despite their importance, thereare few studies in the literature on flows with constant pressure gradient, mainly due to the difficulty ofsetting up this configuration with sufficient accuracy. In experiments, a desired β distribution at the edgeof a TBL is usually achieved by changing the wall shape of wind tunnels through a trial-and-error process,see e.g. Sanmiguel Vila et al. [52]. This is, in fact, very time consuming, since β is needed to be measuredin the streamwise direction at each iteration. Consequently, it is very difficult to achieve a long constant- β region in the streamwise direction. If a proper shape of the wall could be known prior to an experiment, notonly time and cost could be saved, but also a more accurate constant- β distribution would be eventuallyobtained.The objective of this section is to efficiently obtain the optimal shape of the upper wall of a channelusing the Bayesian optimization described in Section 2 so that a target β along the edge of the TBL overthe lower wall is achieved. We consider both zero-pressure-gradient (ZPG) and adverse-pressure-gradient(APG) TBLs with constant values for the target β . Besides these, an optimization is considered where thetarget β distribution varies in the streamwise direction. Such distribution is taken from the flow arounda NACA0012 airfoil. This demonstrates the applicability of the BO method to a wide variety of targetdistribution for the pressure gradient β . The considered flows have high Reynolds number and are hence turbulent. Therefore, we performincompressible Reynolds-averaged Navier–Stokes (RANS) simulations using the open-source software Open-FOAM [70]. In particular, the simpleFoam solver is used which is based on the SIMPLE scheme [47] forvelocity-pressure coupling. The k − ω shear-stress transport (SST) turbulence model introduced by Menteret al. [35] is used, since for TBLs this model has shown better agreement with experimental data comparedto other models, see e.g. Vinuesa et al. [66]. Note that the default values are used for the RANS modelcoefficients.A two-dimensional domain is considered to mimic the conditions that may exist in a wide wind tunnel, seethe schematic in Fig. 3. The domain has a fixed-shape flat-plate lower wall and an upper wall whose shapeis subject to optimization. The domain length in the streamwise direction and the domain height at theinlet are denoted by L x and L in y , respectively. In particular, we consider L x /δ in99 = 1000 and L in y /δ in99 = 40,where δ in99 is the 99% boundary layer thickness at the inlet. In general, δ is defined as the distance fromthe wall at which the mean streamwise velocity becomes 99% of the local edge velocity, U e . In the presentstudy, U e represents the local maximum velocity of the lower-wall’s TBL. Note that the diagnostic-plotscaling method proposed by Vinuesa et al. [67] is not applicable to RANS due to the lack of informationabout the instantaneous velocity. Note also that δ in99 pertaining to the upper and lower TBLs are the samesince the inflow condition for the ZPG-TBLs over the upper and lower walls are taken to be the same. Theinflow profiles are assigned based on the DNS data of ZPG-TBL from Schlatter and ¨Orl¨u [53], see the detailsin Ref. [38]. The upper wall is constructed by using a spline function to smoothly connect the intersection8 DNSinflow
Figure 3: Schematic of the domain used in the shape optimization of the diverging channel. The black dashed lines illustratethe TBLs over the upper and lower walls. Red circles denote the control points, whose heights q i are to be optimized. Theoptimization range is set to avoid any inlet or outlet effects on the TBLs. of the channel inlet and the upper wall, i.e. (cid:0) x = 0 , y = L in y (cid:1) , and d control points. The control pointsare equally spaced in the streamwise direction. The heights y of these points are subject to optimization,therefore the ( x, y ) coordinates of the i -th control point are defined as ( iL x /d, q i ) where i = 1 , , · · · , d ,see Fig. 3. When applying the BO of Section 2, for each sample of q a new geometry for the domain isobtained for which a structured hexahedral mesh is generated using OpenFOAM standard meshing function, blockMesh . In order to resolve the near-wall region of the TBL, a fine mesh resolution is needed near the wall.The mesh is constructed in such a way that y +1 = u τ y /ν remains below unity at the walls. Here, y is thedistance of the first computational cell center from the lower wall, ν is the kinematic viscosity, u τ = (cid:112) τ w /ρ is the wall-friction velocity with τ w and ρ respectively denoting the wall-shear stress and fluid density.The aim of the optimization is to obtain a pressure-gradient distribution sufficiently close to a giventarget β t . Therefore, we can formulate a minimization problem whose objective function is defined as anerror between the computed β distribution and the target distribution, r ( q ) = (cid:107) β ( q , x ) − β t ( x ) (cid:107) L , (14)where, (cid:107) · (cid:107) L denotes a standard L -Lebesgue norm evaluated over the domain of x . Note that regionsclose to the inlet and outlet of the domain have to be excluded when evaluating r ( q ) to avoid any inlet oroutlet effects on the β distributions, see Fig. 3. The remaining range of x over which the norm in Eq. (14)is computed is referred to as optimization range which may be chosen adequately for each optimization caseas explained in Section 5.3. For additional details about the problem setup, refer to Ref. [38]. As explained in the Section 5.1, optimization is conducted for three different distributions of β : twoconstant values of zero (ZPG) and one (APG), and a non-constant distribution corresponding to the flowaround a NACA0012 airfoil. These optimization cases are named as Constant-0, Constant-1 and Wing,respectively, and are discussed in the following subsections. The conditions of these cases are summarizedin Table 1. Prior to the optimization, numerical experiments are conducted to decide about the number ofcontrol points, admissible range of parameters, optimization range, and the total number of iterations. β distribution For the Constant-0 case, essentially a correction for the displacement effect of the growing boundarylayers is expected, which yields a comparably simple shape of the upper wall. Thus we chose to only use9 able 1: Summary of the three optimization cases. The admissible space Q is non-dimensionalized by δ in99 , i.e. ˜ Q = Q /δ in99 . Thetotal number of iterations is denoted by n . Case name Constant-0 Constant-1 Wing β t Q ˜ Q = [40 , Q = [40 ,
46] ˜ Q = [55 , Q = [75 , Q = [85 , Q = [95 , Q = [40 ,
42] ˜ Q = [40 , Q = [40 ,
45] ˜ Q = [40 , Q = [40 ,
50] ˜ Q = [40 , Q = [45 ,
60] ˜ Q = [50 , < x/δ in99 <
900 300 < x/δ in99 <
900 150 < x/δ in99 < n
25 50 802 points along the length of the channel. The optimization of the Constant-0 case is conducted and theoptimum is found at the 20-th iteration. Note that the optimum is taken as the best value over the first n = 25 iterations. Due to the definition of the objective function, see Eq. (14), the validity of the computedoptimal parameters can be examined. The results of the optimization are shown in Fig. 4. It is observedthat the global minimum of r is found in the considered admissible space and the associated β distributionis highly accurate (within ± .
005 of the target value of zero). The 95% confidence interval of r ( q ), which isshown in Fig. 4(d), is smaller close to the optimum because of the higher number of samples.Similar to the Constant-0 case, optimization of the Constant-1 case, i.e. for β t = 1 is conducted. Sincethe target distribution of β is constant but larger than zero, the optimum shape of the upper wall isexpected to be only moderately complex; thus, 4 control points are used. The optimum is found at the 40-thiteration, whose results are shown in Fig. 5. The optimal shape of the upper wall is steeper than that of theConstant-0 case, since the target β t is higher. In the Constant-1 case, the computed optimal β distributionremains almost constant in the optimization range and deviates from the target β t by less than ± .
02. It isnoteworthy that since the inflow is taken from a ZPG-TBL, a development length is needed to reach β = 1.This development length is chosen to be 300 δ in99 . We observed that forcing a faster development of the β through imposing a shorter development length, could lead to the separation of the TBL over the upper walland hence make it more difficult to achieve a constant β for the lower-wall’s TBL. As separation is typicallyunwanted in a wind-tunnel setting, we prefer to have a longer development length avoiding this issue.To validate the results of the constant- β TBLs obtained from the optimization, the streamwise devel-opment of the skin-friction coefficient c f and the shape factor H are compared to the reference data, seeFig. 6. The skin-friction coefficient and the shape factor are defined as c f = 2( u τ /U e ) and H = δ ∗ /θ ,where δ ∗ and θ respectively denote the displacement and momentum thickness. As the reference data forthe constant- β TBLs, the correlations based on the data compiled by Vinuesa et al. [68] and the experi-mental data from Sanmiguel Vila et al. [52] are adopted. The Reynolds number based on the momentumthickness, Re θ = U e θ/ν , is used as the horizontal axis of the plots in Fig. 6, so that the direct comparisonis possible. Clearly, the optimized cases show excellent agreement with the reference data in the regionwhere β is constant (Re θ (cid:38) β = 0 to β = 1 profilebecause of the development of β in the streamwise direction, see Fig. 5(b). Note that the experimental β distribution of Sanmiguel Vila et al. [52] has a variation in the streamwise direction (about ± β distributions in the experiments, as mentioned in Section 5.1. Asa conclusion, not only the applicability of the optimization method to obtain constant- β distributions isconfirmed but also the validity of the resulting constant- β TBLs is approved. β distribution The same optimization procedure can be applied when β t is not constant in the streamwise direction. Thetarget β distribution is taken from the flow around a NACA0012 airfoil. Note that due to the geometricalsymmetry with respect to the airfoil chord, β distributions of the suction and pressure sides are identical forNACA0012 at zero angle of attack. The distribution of the target β is taken from the numerical simulation byTanarro et al. [62]. Since the target β is originally defined based on the normalized coordinate x/c , where c
200 400 600 800 1000 x/δ in99 y / δ i n . . . . . . U/U in e
40 41 42 43 44 q /δ in99 q / δ i n . . . . . . . . . . P r e d i c t e d M e a n o f r ( q ) (a) (b) x/δ in99 − − β × − Result of the optimizationTarget distribution
40 41 42 43 44 q /δ in99 q / δ i n . . . . . . . . . . % C I (c) (d) Figure 4: Results of the optimization of the Constant-0 case using 2 control points: (a) contours of the streamwise velocity,(b) contours of the mean prediction of the objective function (Eq. (14)), (c) comparison of the computed β distribution withthe target value, (d) contours of the 95% confidence interval of the objective function. The red markers show the optimalparameter q opt /δ in99 = (42 . , . δ of the TBL over the lower wall. The black hollow markersand dashed lines in (b) and (d) show the sampled parameters and their trace, respectively. Note that the streamwise velocityis non-dimensionalized by the edge velocity at the inlet, U in e . denotes the chord length, it has to be mapped to the channel normalized coordinate in the streamwisedirection, i.e. x/δ in99 . Here, the original data are mapped to the first 95% of the channel length, i.e. 850 is less than 0 . 54. This promising result proves that the Bayesian optimization methodcan also be used for more complex pressure-gradient distributions. This is, in fact, very useful for futurewind-tunnel experiments to study wall-bounded turbulent flows with arbitrary β distributions, noting thatthe common approach in practice is based on trial-and-error, which can be inaccurate and time consuming.The results of the present section have been obtained by RANS, which may have its own deficiencyin fidelity compared to e.g. large eddy simulations or even direct numerical simulations. However, the11 200 400 600 800 1000 x/δ in99 y / δ i n . . . . . . U/U in e x/δ in99 . . . . β Result of the optimizationTarget distribution (a) (b) Figure 5: Results of the optimized flow in Constant-1 case using 4 control points: (a) contours of the streamwise velocity, (b)comparison of the computed β distribution with the target value. The optimal parameter is q opt /δ in99 = (62 . , . , . , . Re θ c f × − Constant-0Constant-1Vila et al. ( β ≈ β = 0)Vinuesa et al. ( β = 1) Re θ . . . . . H Constant-0Constant-1Vila et al. ( β ≈ β = 0)Vinuesa et al. ( β = 1) (a) (b) Figure 6: Validation of the optimized Constant-0 and Constant-1 cases: (a) streamwise development of the skin-frictioncoefficient c f , (b) streamwise development of the shape factor H . Reference correlations and experimental data are based onthe data compiled by Vinuesa et al. [68] and Sanmiguel Vila et al. [52], respectively. underlying flow simulation method does not affect the performance of the optimization procedure. Therefore,the current results could be even achieved using higher-fidelity methods, or even – if properly implemented– in an experimental setup. 6. Optimization of a spoiler-ice airfoil model Finally, in this section we present a more applied flow case which shows the potential of the BO frameworkalso for industrial flow cases where optimal designs are sought. When operating in cold climate, the performance of the wind turbines can be reduced by icing [8, 60].In an extreme condition, heavy icing can even lead to a complete stop of the turbine. Icing is also a riskycondition for airplane wings since it may induce flow separation at a small angles of attack, which thenmight lead to stall and consequently loss of lift. Ice accretion is a complex process which depends on bothaerodynamics and thermodynamics. The process is affected by many parameters, for instance, ambienttemperature, surface properties, relative speed of the airfoil, and the time span over which the icing event12 200 400 600 800 1000 x/δ in99 y / δ i n . . . . . . U/U in e x/δ in99 β Result of the optimizationTarget distribution (a) (b) Figure 7: Results of the optimized flow in Wing case using 8 control points: (a) contours of the streamwise velocity, (b)comparison of the computed β distribution with the target β distribution. The target distribution is taken from the simulationdata by Tanarro et al. [62]. The optimal parameter values are q opt /δ in99 = (40 . , . , . , . , . , . , . , . (a) (b) Figure 8: The schematics of the horn ice (a) and streamwise ice (b). The black, red and blue lines respectively illustratethe surface of a clean airfoil, the actual ice shape, and an example of the spoiler-ice configuration employed to model theaerodynamic effects of the actual ice. takes place. As a result, numerous possibilities for an ice profile exist. Bragg et al. [8] categorized theleading-edge ice on the airfoils into four types: roughness, horn, streamwise, and spanwise-ridge ice. In thisstudy, we focus on the horn and streamwise ice, see Fig. 8.The main characteristic of the horn ice is the presence of the large protuberances which induce a large flowseparation downstream of the ice, and hence dramatically affect the aerodynamic performance, includingthe lift, drag and moment coefficients. On the other hand, streamwise ice is smoothly formed along thestreamlines and is less dangerous in the sense of having adverse effects on the aerodynamic performance.Modeling the iced-airfoil is a challenging task in both CFD and experiments since the ice profile essentiallyhas infinite degrees of freedom. To rectify this, the simplified “simulated-ice” concept has been consideredas a representative of the actual ice profile, see Refs. [45, 44, 46, 59]. In particular, the thin plates which arecalled the “spoilers” and have the same height as the original ice can be used to reproduce the aerodynamicperformance of the ice, see Fig. 8. The resulting model is called the “spoiler-ice” model. Although theheight, location, and angle of the spoilers have been considered as important parameters which affect theairfoil aerodynamic performance, see e.g. Refs. [45, 44], Tabatabaei et al. [58] suggested that the thicknessof the spoilers should be also taken into account as an effective parameter.13 . . . . . . x/c − . . . . y / c Original iceClean airfoil Figure 9: Original iced airfoil used in the present study. The ice is formed on a clean NACA64618 airfoil, which is shownwith a black solid line. The original ice shape is illustrated with a red solid line and a filled gray region. Coordinates arenon-dimensionalized by the chord length of the clean airfoil, c . Note that the origin of the x coordinate is set as the leadingedge of the clean airfoil, and this is the convention throughout Section 6. Although parametric studies of the spoiler-ice model have been conducted to investigate associated effectson the aerodynamic performance, see e.g. Refs. [45, 44], it is not yet clearly known which configuration ofthe spoiler is appropriate to accurately represent the aerodynamic performance of an arbitrary given realice profile. The BO method introduced in Section 2 can be utilized to find the optimal parameters of thespoiler configuration. Although the spoiler-ice concept has been previously used mainly for the horn ice,see e.g. Ref. [46], here we also consider to apply the idea to the streamwise ice. Therefore, an originalarbitrary ice shape is chosen so that it has both horn and streamwise ice on the suction and pressure sideof the airfoil, respectively. Two spoilers on the suction and pressure sides (the upper and lower spoilers,respectively, shown in Fig. 10(a)) are used to obtain a similar effect of the original iced airfoil. Since aspoiler-ice with the same height as the original ice shape has been used in the previous studies [45, 44, 46],here the optimization is first conducted with the constraint of keeping the heights of the spoilers fixed.However, it is shown in Fig. 14 that the results of the optimization without such a constraint are in betteragreement with those of the original ice. It will thus be shown that a more complete parameter optimizationmay indeed be relevant for accurate simulations. In this section, we describe the setup and implementation of the icing flow case, including the basicsimulation, the meshing and optimization loop. As clean airfoil, a NACA64618 profile is considered, which has been used in a 5MW NREL (NationalRenewable Energy Laboratory) wind-turbine blade. For the purpose of the optimization, an arbitraryoriginal ice shape is required as a reference which is taken from the ice formed on the leading edge of theairfoil studied in Ref. [21], see Fig. 9.The ice on the upper (suction) and lower (pressure) sides of the airfoil is of the horn and streamwisetypes, respectively. As mentioned in Section 6.1, the flow around the original iced airfoil can be modeledby two spoilers each located on one side of the airfoil. For the purpose of optimization, the configuration ofthe spoilers is controlled by seven parameters: heights, angles and widths of the upper and lower spoilers, h u , h l , θ u , θ l , w u and w l , respectively, and the displacement of the lower spoiler, ( d l ), see Fig. 10(a).The displacement of the upper spoiler is not taken into account since the original ice shape on the suctionside is horn ice, which has a clear connection between the position of the horn of the original ice and theposition of the spoiler-ice. On the other hand, the pressure side of the original ice is of streamwise type,hence it is not, in general, clear where the spoiler-ice should be located. To rectify this, the displacementparameter d l is needed for the lower spoiler to optimize its position as well as the other parameters.Prior to the optimization, different numerical experiments, detailed in Section 6.2.2, with different combi-nations of the parameters are conducted to decide about the admissible range of the parameters. Fig. 10(b)shows the configuration of the spoilers when the design parameters are at either end of their admissiblerange. Note that the upper limit of the admissible range of the heights, h u and h l , is the same as the highestpoint of the original ice profile. 14 a)Design parameters θ u , θ l h u , h l w u , w l d l Lower limit of the rangeUpper limit of the range (b) Figure 10: Schematic representation of the spoiler ice on the airfoil to model the original ice profile: (a) design parameters,(b) configuration when the design parameters are at the either end of their admissible range. The red circles in (a) specify therotation center of θ l and θ u . The direction of d l is defined to be parallel to the tangential direction at the root of the lowerspoiler. The red dashed line in (b) shows the original ice profile. Note that in the particular configuration at (b) all parametersare fixed at the middle of their admissible range except the specified design parameters. Numerical simulations are conducted considering a cross section of an iced-airfoil inside a (virtual) windtunnel. Two-dimensional steady and incompressible RANS are performed using OpenFOAM with the samesolver settings and turbulence model as mentioned in Section 5.2. The computational domain, see Fig. 11,is made based on the test section of the Minimum-Turbulence-Level (MTL) wind tunnel, which is a closed-loop circuit tunnel housed at KTH. Additional information about the MTL wind tunnel can be found in e.g. Ref. [33].A uniform inflow condition for the velocity is applied at the inlet of the domain. The resulting Reynoldsnumber based on the airfoil chord is Re c = U ∞ c/ν = 4 . × , where U ∞ and c , respectively, denotethe freestream velocity and the chord length of the airfoil. No-slip/no-penetration boundary conditions areimposed at the top and bottom walls of the domain, as well as on the airfoil, ice, and the spoiler ice.The commercial software Ansys ICEM CFD y +1 is kept below 1 all around the airfoil. The resulting total number of computational cells isabout 1 million. Note that the number of computational cells slightly differs for the simulations associatedwith a particular configuration of the spoiler ice. To put the computational cost in perspective, each flowsimulation takes about one hour running in parallel on 24 cores of Intel Xeon E5-2690v3 Haswell processors.To measure the similarity of the aerodynamic performance of a spoiler-ice model and the actual iceprofile, the pressure coefficient ( c p ) distribution around the airfoil is considered, since it is directly related15 igure 11: Schematic of the domain used in the iced-airfoil simulations. The chord length of the airfoil and the freestreamvelocity are denoted by c and U ∞ , respectively. to the aerodynamic characteristics such as the lift and drag coefficients. The objective function in theoptimization is thus defined as the L -norm of the error between the computed c p distributions on thesuction and pressure sides of the airfoil, and the corresponding distribution of the original iced airfoil. Thepressure coefficient is defined as c p = 2( p − p ∞ ) / ( ρU ∞ ), where p and p ∞ denote the local pressure on thesurface of the airfoil and the freestream pressure at the inlet, respectively. The objective function reads, r ( q ) = Σ f =1 (cid:115)(cid:90) c . c | c p f ( q , x ) − c p f ,t ( x ) | d x , (15)where Σ f =1 specifies the summation over the suction and pressure sides of the airfoil and c p,t denotes the c p distribution of the original iced airfoil. Note that the objective function is evaluated only downstream ofthe original ice, meaning that the optimization range is set to 0 . < x/c < c p,t . The resulting velocity contours and pressure-coefficient distributionare shown in Fig. 12. Because of the ice formed on the leading edge, there are separation bubbles rightafter the ice on both suction and pressure sides. Moreover, the c p near the leading edge exhibits suddenchanges in the streamwise direction. Note that x/c (cid:46) . x/c = 0 . 086 and 0 . . < x/c < 1, as represented in Fig. 12. Optimization of the design parameters is conducted with two configurations: spoilers with fixed heightand flexible height. As mentioned in Section 6.1, glaze ice has been simulated by spoiler ice whose heightis the same as the highest point of the original ice [45, 44, 46]. As a result, one optimization is conductedwith the constraints of having the heights of the spoilers fixed, and the case is referred to as “fixed-heightoptimization”. Imposed by these constraints, the number of design parameters is reduced to d = 5. Sincethe height of the spoiler ice is a relevant parameter on the aerodynamic performance [60, 15], the secondoptimization is conducted without the above constraints, and the case is called “flexible-height optimization”.Allowing the height of the spoilers to change, the number of design parameters becomes d = 7. The maximumnumber of iterations in the BO-GPR algorithm of Section 2 is chosen as n = 50 and 90 for fixed and flexible-height optimizations, respectively. These limits are chosen based on the numerical experiments conductedprior to the optimization. The optimal parameters are found for the 13-th and the 43-rd iterations forthe fixed and flexible-height optimizations, respectively. The flow fields and c p distributions of the optimalconfigurations are compared to those of the original iced airfoil in Figs. 13 and 14, in order to validate theresults of the optimizations.In Fig. 13, the separation bubbles and vortex structures are illustrated as they appear behind the iceprofile and the spoiler ice on both pressure and suction sides. The general structure is in agreement withthe previous studies [44, 59], where the primary and secondary vortices were observed around the spoiler-icethrough numerical simulations. In fact, it is not even necessary that the flows downstream of the original16 . 25 0 . 00 0 . 25 0 . 50 0 . 75 1 . 00 1 . x/c − . − . . . . y / c . . . . . . . U/U ∞ (a) − . 25 0 . 00 0 . 25 0 . 50 0 . 75 1 . 00 1 . x/c − − c p Downstreamof the ice Original ice (suction side)Original ice (pressure side) (b) Figure 12: Results of the original iced-airfoil case: (a) contours of the streamwise velocity, (b) the pressure-coefficient distri-bution. The velocity is non-dimensionalized by the freestream velocity at the inlet. The black dash-dot line in (a) shows thebeginning of the optimization range. The gray area in (b) is the range excluded from the optimization. Note that the originof the x coordinate is set at the leading edge of the clean airfoil, so the iced airfoil starts from negative x/c , see Fig. 9. Therange of the horizontal axes are chosen to be the same so that comparison between two figures becomes easier. ice and the spoiler ice have similar vortex structures since the primary purpose of using the spoiler-ice isto model the aerodynamic characteristics on an iced airfoil. The flexible-height optimization successfullyyields a reattachment point similar to the original iced-airfoil case, while the reattachment point is pushedfurther downstream for the fixed-height optimization. This has a direct influence on the c p distribution,see Fig. 14. In fact, the c p distribution of the fixed-height case has four times larger deviation from thetarget distribution, i.e. r in Eq. (15), compared to the flexible-height case. It is also observed that forthe fixed-height optimization, the optimal parameters q opt are located at the borders of the consideredadmissible space. This means that the global optimum could change value if the admissible space wouldbe different. Recall that the admissible range of the parameters are limited due to the geometrical andmeshing constraints. In contrast, for the flexible-height optimization, the optimum q opt is found withinthe admissible space and leads to a better agreement with the result of the original iced airfoil. Therefore,despite being used in the literature, keeping the heights in the spoiler-ice model fixed is not a sufficientlyaccurate assumption, and considering also the height as an design parameter leads to higher fidelity in theresults.The obtained c p distribution from the flexible-height optimization is almost identical to the target dis-tribution within the optimization range except on the pressure side around x/c = 0 . 2. The slight deviationsare due to the fact that the original ice on the pressure side is of the streamwise type, which is basicallymore difficult to be represented by a spoiler-ice model.17 SecondaryVortex PrimaryVortex (a) SecondaryVortex PrimaryVortex (b) SecondaryVortex PrimaryVortex (c) Figure 13: Comparison of the flow structures around the leading edge. The line integral convolution (LIC) technique [10] hasbeen used to visualize the velocity field along with the contours of the streamwise velocity for (a) the original iced airfoil, (b)the optimized spoiler ice with fixed heights, and (c) the optimized spoiler ice with flexible heights. . . . . . . x/c − − c p Downstreamof the ice Original iceSpoiler-ice: fixed heightSpoiler-ice: flexible height Figure 14: Comparison between the distributions of the pressure coefficient obtained for the original ice profile (shown inFig. 12), and the optimized spoiler-ice models. The corresponding flow fields are represented in Fig. 13. The solid and dashedlines specify the suction and pressure sides of the airfoil, respectively. The black lines represent the target distribution, c p,t ,which is also shown in Fig. 12(b). Note that the optimization range is 0 . < x/c < 1, which specifies a region downstream ofthe original ice. The BO has also shown to be an efficient optimizer even for this applied CFD case. The coupling to theCFD code and meshing software are quite straight-forward, and can be done non-intrusively. Note that noderivatives of the target function are needed, which makes the present approach suitable for general CFDapplications. 7. Summary and conclusions The Bayesian optimization based on Gaussian process regression (BO-GPR) has been applied to differentCFD problems ranging from purely academic to industrially relevant setups, using state-of-the-art simula-tion methods. The diversity of the examples with different numbers of design parameters helps to betterunderstand the applicability and performance of the BO-GPR approach which has not been frequently usedin CFD and turbulent flows despite being a primary choice in other fields, e.g. data sciences. The use of theBO-GPR is straightforward noting that the Bayesian optimization is among the gradient-free approachesand hence only requires forward evaluation of the quantities of interest (QoIs), i.e. running a CFD code.In the BO, a sequentially-updating set of samples is taken from the space of the design parameters, suchthat the sequence converges to the global optimum after a finite number of iterations. A main advantageof the BO-GPR is its versatility, meaning the algorithm is non-intrusive and can be utilized with any CFDcode and for any number of design parameters. The use of a GPR surrogate provides the possibility ofestimating the confidence in the values of the cost function over the admissible space of the parameters.What is needed to create the optimization loop is to automatize the whole process through appropriate linksbetween different scripts and codes. The software developed in the present work is based on GPy [24] and GPyOpt [63] libraries and will be made available online.Although the optimization problems studied in the present study contain as many as 8 design parameters,the maximum number of iterations to ensure obtaining a reliable global optimum is no more than 90.This number should be compared with the number of function evaluations necessary e.g. for adjoint-basedoptimizations. More interestingly, by changing the number of design parameters from 2 to 8, the number ofrequired flow simulations does not significantly increase. Therefore, the computational effort for BO-GPRis observed to be affordable for many practical applications, especially if the RANS (Reynolds-averagedNavier–Stokes) approach is employed for the simulation of turbulence. However, similar to other gradient-free approaches, the BO-GPR may become inefficient for a large number of design parameters. As shownfor instance by Wu et al. [71], the efficiency of BO-GPR can be improved by adding sensitivity information19rom the gradients of the cost function computed by adjoint methods. Application and analysis of suchmodifications will be considered in the future extensions of the present study.We demonstrate the BO-GPR approach on two comparably simple problems in Sections 3 and 4 for ananalytical test problem and a lid-driven cavity. In both cases, our approach could find the global optimum,as compared to previous literature results of the same cases. The optimization required about 20 iterationsfor the four-dimensional problem, exploring the complete parameter space in an efficient manner. It can beobserved that the Bayesian technique quickly focuses on the regions where the potential optima are located.In addition, the uncertainty information about the optima can be retrieved from the approach.A more physically relevant problem is considered on the example of the shape optimization of theupper wall of a two-dimensional channel. Here, the goal is to obtain a given target pressure-gradient (PG)distributions in a turbulent boundary layer developing over a flat wall, which is achieved by adjustingthe location of the upper wall in a variable-height channel flow. The flow is simulated using the open-source finite-volume-based solver OpenFOAM [70] adopting the RANS approach to model the turbulence.Accurate results are obtained for the target pressure gradient regardless of being constant or varying in thestreamwise direction. Different numbers of design parameters 2, 4, and 8 are considered for which 25, 50,and 80 simulations are respectively needed to find the corresponding global optimum. The results of theoptimum configuration for the constant-PG flows show an excellent agreement with the desired pressure-gradient distribution, and as consequence also with the correlations and experimental data reported in theliterature for integral quantities of pressure-gradient boundary layers. In addition, promising results for thenon-constant target PG distribution suggests that the BO-GPR approach can be applied to various practicalapplications including the design of inserts in wind tunnels. This is of special interest noting that the currentapproach for such designs is based on trial and errors, with limited accuracy and robustness.Finally, we also apply the Bayesian optimization to an industrially-relevant case, in an effort to showthat the approach can readily be used in an applied context. The goal is to optimize the configuration ofthe spoiler ice which is a simplified yet useful model for the actual ice profiles formed on the airfoils in coldconditions. The objective is to match, within a small margin of error, the resulting pressure distribution ofan airfoil with modeled ice with that of the same airfoil but with actual ice profile. The RANS simulationsusing OpenFOAM for design parameters of dimension 5 and 7 are examined for which the maximum numberof iterations in the optimization is found to be equal to 50 and 90, respectively. The encouraging resultsof the optimization prove the applicability of the BO-GPR framework to a relatively complex parameteroptimization and also re-confirm the validity of the spoiler-ice models. It is also observed that keeping theheights in the spoiler-ice model fixed, as it is convenient in the literature, would lead to inaccurate flowsimulations. Moreover, although the spoiler-ice concept has been mainly used in the literature for modelingthe horn ice, see Refs. [45, 44, 46], the present study reveals that the model can also be applicable to thestreamwise ice, conditioned on adopting the optimal parameters.The present study can be extended in several ways. Our future attempts will include the use of scale-resolving approaches for simulation of the turbulent flows, instead of RANS which has been adopted here.In this way, even the effect of time-averaging on transient data can be considered. In addition, largernumbers of design parameters could be inquired, to identify the impact of dimensionality on the completeoptimization problem. Acknowledgments YM acknowledges the Keio-KTH double degree program and the financial support from the NSK Schol-arship Foundation. SR acknowledges the financial support from the Linn´e FLOW Centre at KTH andthe EXCELLERAT project which has received funding from the European Union’s Horizon 2020 researchand innovation programme under grant agreement No 823691. PS, NT and SR also acknowledge fundingby the Knut and Alice Wallenberg Foundation via the KAW Academy Fellow programme. KF acknowl-edges the financial support from the Japan Society for the Promotion of Science (KAKENHI grant number:18H03758). The simulations in Section 6 were performed on the resources provided by the Swedish NationalInfrastructure for Computing (SNIC) at PDC (KTH Royal Institute of Technology).20 eferences [1] ANSYS Inc. ICEM CFD User Manual , 2011. Release 14.0.[2] T. Back. Evolutionary Algorithms in Theory and Practice: Evolution Strategies, Evolutionary Programming, Genetic Al-gorithms . Oxford University Press, 1996. ISBN 9780195356700. URL https://books.google.de/books?id=htJHI1UrL7IC .[3] M. Berggren. Numerical solution of a flow-control problem: Vorticity reduction by dynamic boundary ac-tion. SIAM Journal on Scientific Computing , 19(3):829–860, May 1998. doi: 10.1137/s1064827595294678. URL https://doi.org/10.1137/s1064827595294678 .[4] H. Bijl, D. Lucor, S. Mishra, and C. Schwab, editors. Uncertainty Quantification in ComputationalFluid Dynamics . Springer International Publishing, 2013. doi: 10.1007/978-3-319-00885-1. URL https://doi.org/10.1007/978-3-319-00885-1 .[5] A. Bobke, R. Vinuesa, R. ¨Orl¨u, and P. Schlatter. History effects and near equilibrium in adverse-pressure-gradientturbulent boundary layers. Journal of Fluid Mechanics , 820:667–692, May 2017. doi: 10.1017/jfm.2017.236. URL https://doi.org/10.1017/jfm.2017.236 .[6] V. B. L. Boppana and J. S. B. Gajjar. Global flow instability in a lid-driven cavity. Interna-tional Journal for Numerical Methods in Fluids , 62(8):827–853, 2010. doi: 10.1002/fld.2040. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.2040 .[7] G. Box and N. Draper. Empirical Model-Building and Response Surfaces . Wiley Series in Probability and Statistics.Wiley, 1987. ISBN 9780471810339. URL https://books.google.se/books?id=QO2dDRufJEAC .[8] M. Bragg, A. Broeren, and L. Blumenthal. Iced-airfoil aerodynamics. Progress in Aerospace Sciences , 41(5):323–362, July2005. doi: 10.1016/j.paerosci.2005.07.001. URL https://doi.org/10.1016/j.paerosci.2005.07.001 .[9] A. D. Bull. Convergence rates of efficient global optimization algorithms. J. Mach. Learn. Res. , 12(null):2879–2904, Nov.2011. ISSN 1532-4435.[10] B. Cabral and L. C. Leedom. Imaging vector fields using line integral convolution. In Proceedings of the 20thAnnual Conference on Computer Graphics and Interactive Techniques , SIGGRAPH ’93, page 263–270, New York,NY, USA, 1993. Association for Computing Machinery. ISBN 0897916018. doi: 10.1145/166117.166151. URL https://doi.org/10.1145/166117.166151 .[11] H.-S. Chung and J. Alonso. Using gradients to construct cokriging approximation models for high-dimensional designoptimization problems . 2002. doi: 10.2514/6.2002-317. URL https://arc.aiaa.org/doi/abs/10.2514/6.2002-317 .[12] F. H. Clauser. Turbulent boundary layers in adverse pressure gradients. Journal of the Aeronautical Sciences , 21(2):91–108, Feb. 1954. doi: 10.2514/8.2938. URL https://doi.org/10.2514/8.2938 .[13] F. H. Clauser. The turbulent boundary layer. In Advances in Applied Mechanics , pages 1–51. Elsevier, 1956. doi:10.1016/s0065-2156(08)70370-3. URL https://doi.org/10.1016/s0065-2156(08)70370-3 .[14] S. Deck, F. Gand, V. Brunet, and S. B. Khelil. High-fidelity simulations of unsteady civil aircraft aerodynamics: stakesand perspectives. Application of zonal detached eddy simulation. Philosophical Transactions of the Royal Society A:Mathematical, Physical and Engineering Sciences , 372(2022):20130325, Aug. 2014. doi: 10.1098/rsta.2013.0325. URL https://doi.org/10.1098/rsta.2013.0325 .[15] A. M. DeGennaro, C. W. Rowley, and L. Martinelli. Uncertainty quantification for airfoil icing using poly-nomial chaos expansions. Journal of Aircraft , 52(5):1404–1411, 2015. doi: 10.2514/1.C032698. URL https://doi.org/10.2514/1.C032698 .[16] M. O. Deville, P. F. Fischer, and E. H. Mund. High-Order Methods for Incompressible Fluid Flow . Cambridge Monographson Applied and Computational Mathematics. Cambridge University Press, 2002. doi: 10.1017/CBO9780511546792.[17] C. B. Dilgen, S. B. Dilgen, D. R. Fuhrman, O. Sigmund, and B. S. Lazarov. Topology op-timization of turbulent flows. Computer Methods in Applied Mechanics and Engineering , 331:363 – 393, 2018. ISSN 0045-7825. doi: https://doi.org/10.1016/j.cma.2017.11.029. URL .[18] Y. Fan, W. Li, M. Atzori, R. Pozuelo, P. Schlatter, and R. Vinuesa. Decomposition of the mean friction drag inadverse-pressure-gradient turbulent boundary layers. Phys. Rev. Fluids , 5:114608, Nov 2020. doi: 10.1103/PhysRevFlu-ids.5.114608. URL https://link.aps.org/doi/10.1103/PhysRevFluids.5.114608 .[19] P. F. Fischer, J. W. Lottes, and S. G. Kerkemeier. NEK5000: Open source spectral element CFD solver. Available at: http://nek5000.mcs.anl.gov , 2008. URL http://nek5000.mcs.anl.gov .[20] P. Frazier. A tutorial on Bayesian optimization. ArXiv , abs/1807.02811, 2018.[21] S. Gantasala, N. Tabatabaei, M. Cervantes, and J.-O. Aidanp¨a¨a. Numerical investigation of the aeroelastic behavior ofa wind turbine with iced blades. Energies , 12(12):2422, Jun 2019. ISSN 1996-1073. doi: 10.3390/en12122422. URL http://dx.doi.org/10.3390/en12122422 .[22] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processingfacilities. International Journal for Numerical Methods in Engineering , 79(11):1309–1331, 2009. doi: 10.1002/nme.2579.URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2579 .[23] P. W. Goldberg, C. K. I. Williams, and C. M. Bishop. Regression with input-dependent noise: A Gaussian processtreatment. In Proceedings of the 1997 Conference on Advances in Neural Information Processing Systems 10 , NIPS ’97,page 493–499, Cambridge, MA, USA, 1998. MIT Press. ISBN 0262100762.[24] GPy. GPy: A gaussian process framework in python. http://github.com/SheffieldML/GPy , since 2012.[25] R. B. Gramacy. Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences . ChapmanHall/CRC, Boca Raton, Florida, 2020. URL http://bobby.gramacy.com/surrogates/ .[26] M. Gunzburger. Adjoint equation-based methods for control problems in incompressible, viscous flows. low, Turbulence and Combustion , 65(3/4):249–272, 2000. doi: 10.1023/a:1011455900396. URL https://doi.org/10.1023/a:1011455900396 .[27] M. Hinze, A. Walther, and J. Sternberg. An optimal memory-reduced procedure for calculating adjoints of the instationaryNavier-Stokes equations. Optimal Control Applications and Methods , 27(1):19–40, 2006. doi: 10.1002/oca.771. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/oca.771 .[28] D. R. Jones, M. Schonlau, and W. J. Welch. Efficient global optimization of expensive black-box func-tions. Journal of Global Optimization , 13(4):455–492, 1998. doi: 10.1023/a:1008306431147. URL https://doi.org/10.1023/a:1008306431147 .[29] R. R. Kerswell, C. C. T. Pringle, and A. P. Willis. An optimization approach for analysing nonlinear stability with transitionto turbulence in fluids as an exemplar. Reports on Progress in Physics , 77(8):085901, aug 2014. doi: 10.1088/0034-4885/77/8/085901.[30] J. Larson, M. Menickelly, and S. M. Wild. Derivative-free optimization methods. Acta Numerica , 28:287–404, 2019. doi:10.1017/S0962492919000060.[31] J. Laurenceau, M. Meaux, M. Montagnac, and P. Sagaut. Comparison of gradient-based and gradient-enhanced response-surface-based optimizers. AIAA Journal , 48(5):981–994, 2010. doi: 10.2514/1.45331. URL https://doi.org/10.2514/1.45331 .[32] L. Laurent, R. L. Riche, B. Soulier, and P.-A. Boucard. An overview of gradient-enhanced metamodels with applications. Archives of Computational Methods in Engineering , 26(1):61–106, July 2017. doi: 10.1007/s11831-017-9226-3. URL https://doi.org/10.1007/s11831-017-9226-3 .[33] B. Lindgren and A. V. Johansson. Evaluation of the flow quality in the MTL wind-tunnel. Report TRITA-MEK 2002:13,Dept. of Mechanics, KTH, Sweden, 2002.[34] O. A. Mahfoze, A. Moody, A. Wynn, R. D. Whalley, and S. Laizet. Reducing the skin-friction drag of a turbulent boundary-layer flow with low-amplitude wall-normal blowing within a Bayesian optimization framework. Phys. Rev. Fluids , 4:094601,Sep 2019. doi: 10.1103/PhysRevFluids.4.094601. URL https://link.aps.org/doi/10.1103/PhysRevFluids.4.094601 .[35] F. Menter, M. Kuntz, and R. Langtry. Ten years of industrial experience with the SST turbulence model. Heat and MassTransfer , 4, 01 2003.[36] J. Moˇckus. On Bayesian methods for seeking the extremum. In G. I. Marchuk, editor, Optimization Techniques IFIPTechnical Conference Novosibirsk, July 1–7, 1974 , pages 400–404, Berlin, Heidelberg, 1975. Springer Berlin Heidelberg.ISBN 978-3-540-37497-8.[37] B. Mohammadi and O. Pironneau. Shape optimization in fluid mechanics. Annual Reviewof Fluid Mechanics , 36(1):255–279, 2004. doi: 10.1146/annurev.fluid.36.050802.121926. URL https://doi.org/10.1146/annurev.fluid.36.050802.121926 .[38] Y. Morita. Robust optimisation of the pressure-gradient distribution in turbulent boundary layers. Master’s thesis, KTHRoyal Institute of Technology, Sweden, 2020.[39] Y. Nabae, K. Kawai, and K. Fukagata. Bayesian optimization of traveling wave-like wall deformation for turbulentdrag reduction. In The Ninth JSME-KSME Thermal and Fluids Engineering Conference (TFEC9) , Oct. 2017. PaperTFEC9-1045.[40] T. Nakazawa. Increasing the critical Reynolds number by maximizing energy dissipation problem. In A. Segalini, editor, Proceedings of the 5th International Conference on Jets, Wakes and Separated Flows (ICJWSF2015) , pages 613–620,Cham, 2016. Springer International Publishing. ISBN 978-3-319-30602-5.[41] J. Nocedal and S. Wright. Numerical Optimization . Springer New York, 2006. doi: 10.1007/978-0-387-40065-5. URL https://doi.org/10.1007/978-0-387-40065-5 .[42] N. Offermans, A. Peplinski, O. Marin, and P. Schlatter. Adaptive mesh refinement for steady flows in Nek5000. Computers & Fluids , 197:104352, 2020. ISSN 0045-7930. doi: https://doi.org/10.1016/j.compfluid.2019.104352. URL .[43] T. A. Oliver, N. Malaya, R. Ulerich, and R. D. Moser. Estimating uncertainties in statistics computed from direct numericalsimulation. Physics of Fluids , 26(3):035101, 2014. doi: 10.1063/1.4866813. URL https://doi.org/10.1063/1.4866813 .[44] M. Papadakis, S. Alansatan, and M. Seltmann. Experimental study of simulated ice shapes on a NACA 0011 airfoil ,pages 1–29. 1999. doi: 10.2514/6.1999-96. URL https://arc.aiaa.org/doi/abs/10.2514/6.1999-96 .[45] M. Papadakis, S. Alansatan, and S.-C. Wong. Aerodynamic characteristics of a symmetric NACA section with simulatedice shapes , pages 1–36. 2000. doi: 10.2514/6.2000-98. URL https://arc.aiaa.org/doi/abs/10.2514/6.2000-98 .[46] M. Papadakis, B. G. Laflin, G. Youssef, and T. Ratvasky. Aerodynamic scaling experiments with simulated ice accretions ,pages 1–42. 2001. doi: 10.2514/6.2001-833. URL https://arc.aiaa.org/doi/abs/10.2514/6.2001-833 .[47] S. V. Patankar. Numerical heat transfer and fluid flow . Series on Computational Methods in Mechanics and ThermalScience. Hemisphere Publishing Corporation (CRC Press, Taylor & Francis Group), 1980. ISBN 978-0891165224. URL .[48] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning (Adaptive Computation and MachineLearning) . The MIT Press, 2005. ISBN 026218253X.[49] S. Rezaeiravesh and M. Liefvendahl. Effect of grid resolution on large eddy simulation of wall-bounded turbulence. Physicsof Fluids , 30(5):055106, May 2018. doi: 10.1063/1.5025131. URL https://doi.org/10.1063/1.5025131 .[50] S. Rezaeiravesh, R. Vinuesa, and P. Schlatter. An uncertainty-quantification framework for assessing accuracy, sen-sitivity, and robustness in computational fluid dynamics. arXiv:2007.07071 , 2020. doi: arXiv:2007.07071. URL https://arxiv.org/abs/2007.07071 .[51] P. Sagaut, S. Deck, and M. Terracol. Multiscale and Multiresolution Approaches in Turbulence: LES, DES and Hy-brid RANS/LES Methods: Applications and Guidelines . Imperial College Press, 2013. ISBN 9781848169876. URL ttps://books.google.se/books?id=MzW6CgAAQBAJ .[52] C. Sanmiguel Vila, R. Vinuesa, S. Discetti, A. Ianiro, P. Schlatter, and R. ¨Orl¨u. Experimental realisation of near-equilibrium adverse-pressure-gradient turbulent boundary layers. Experimental Thermal and Fluid Science , 112:109975,Apr. 2020. doi: 10.1016/j.expthermflusci.2019.109975. URL https://doi.org/10.1016/j.expthermflusci.2019.109975 .[53] P. Schlatter and R. ¨Orl¨u. Assessment of direct numerical simulation data of turbulent boundary layers. Journal of Fluid Me-chanics , 659:116–126, July 2010. doi: 10.1017/s0022112010003113. URL https://doi.org/10.1017/s0022112010003113 .[54] A. Scotto Di Perrotolo. A theoretical framework for Bayesian optimization convergence. Master’s thesis, KTH, Optimiza-tion and Systems Theory, 2018.[55] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. de Freitas. Taking the human out of the loop: A review of Bayesianoptimization. Proceedings of the IEEE , 104(1):148–175, Jan 2016. ISSN 0018-9219. doi: 10.1109/JPROC.2015.2494218.[56] J. P. Slotnick, A. Khodadoust, J. J. Alonso, D. L. Darmofal, W. D. Gropp, E. A. Lurie, and D. J. Mavriplis. CFD vision2030 study: A path to revolutionary computational aerosciences, 2014.[57] R. C. Smith. Uncertainty Quantification: Theory, Implementation, and Applications . Society for Industrial and AppliedMathematics, Philadelphia, PA, USA, 2013. ISBN 161197321X, 9781611973211.[58] N. Tabatabaei, M. Cervantes, C. Trivedi, and J.-O. Aidanp¨a¨a. Numerical study of aerodynamic characteristics of asymmetric NACA section with simulated ice shapes. Journal of Physics: Conference Series , 753:022055, sep 2016. doi:10.1088/1742-6596/753/2/022055. URL https://doi.org/10.1088/1742-6596/753/2/022055 .[59] N. Tabatabaei, M. Cervantes, and C. Trivedi. Time-dependent effects of glaze ice on the aerodynamic characteristics ofan airfoil. International Journal of Rotating Machinery , 2018, 2018.[60] N. Tabatabaei, M. Raisee, and M. J. Cervantes. Uncertainty Quantification of Aerodynamic Icing Losses in Wind TurbineWith Polynomial Chaos Expansion. Journal of Energy Resources Technology , 141(5), 04 2019. ISSN 0195-0738. doi:10.1115/1.4042732. URL https://doi.org/10.1115/1.4042732 . 051210.[61] C. Talnikar, P. Blonigan, J. Bodart, and Q. Wang. Parallel optimization for large eddy simulations. Proceedings of theCTR Summer Program , page 315, 2014.[62] ´A. Tanarro, R. Vinuesa, and P. Schlatter. Effect of adverse pressure gradients on turbulent wing boundary layers. Journalof Fluid Mechanics , 883:A8, 2020. doi: 10.1017/jfm.2019.838.[63] The GPyOpt authors. GPyOpt: A Bayesian Optimization framework in Python. http://github.com/SheffieldML/GPyOpt , 2016.[64] D. Th´evenin and G. Janiga, editors. Optimization and Computational Fluid Dynamics . Springer Berlin Heidelberg, 2008.doi: 10.1007/978-3-540-72153-6. URL https://doi.org/10.1007/978-3-540-72153-6 .[65] F. A. C. Viana, T. W. Simpson, V. Balabanov, and V. Toropov. Special section on multidisciplinary design optimization:Metamodeling in multidisciplinary design optimization: How far have we really come? AIAA Journal , 52(4):670–690,Apr. 2014. doi: 10.2514/1.j052375. URL https://doi.org/10.2514/1.j052375 .[66] R. Vinuesa, P. H. Rozier, P. Schlatter, and H. M. Nagib. Experiments and computations of localized pressuregradients with different history effects. AIAA Journal , 52(2):368–384, Feb. 2014. doi: 10.2514/1.j052516. URL https://doi.org/10.2514/1.j052516 .[67] R. Vinuesa, A. Bobke, R. ¨Orl¨u, and P. Schlatter. On determining characteristic length scales in pressure-gradient turbulent boundary layers. Physics of Fluids , 28(5):055101, May 2016. doi: 10.1063/1.4947532. URL https://doi.org/10.1063/1.4947532 .[68] R. Vinuesa, R. ¨Orl¨u, C. S. Vila, A. Ianiro, S. Discetti, and P. Schlatter. Revisiting history effects in adverse-pressure-gradient turbulent boundary layers. Flow, Turbulence and Combustion , 99(3-4):565–587, Aug. 2017. doi: 10.1007/s10494-017-9845-7. URL https://doi.org/10.1007/s10494-017-9845-7 .[69] J. Weatheritt and R. Sandberg. A novel evolutionary algorithm applied to algebraic mod-ifications of the RANS stress–strain relationship. Journal of Computational Physics , 325:22–37, 2016. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2016.08.015. URL .[70] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby. A tensorial approach to computational continuum mechan-ics using object-oriented techniques. Computers in Physics , 12(6):620–631, 1998. doi: 10.1063/1.168744. URL https://aip.scitation.org/doi/abs/10.1063/1.168744 .[71] J. Wu, M. Poloczek, A. G. Wilson, and P. Frazier. Bayesian optimization with gradients. In I. Guyon,U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Ad-vances in Neural Information Processing Systems 30 , pages 5267–5278. Curran Associates, Inc., 2017. URL http://papers.nips.cc/paper/7111-bayesian-optimization-with-gradients.pdf .[72] H. Xiao and P. Cinnella. Quantification of model uncertainty in RANS simulations: A review. Progress inAerospace Sciences , 108:1–31, 2019. ISSN 0376-0421. doi: https://doi.org/10.1016/j.paerosci.2018.10.001. URL ..