TECHNICAL FIELD

The invention relates to design of an aerodynamic component, such as an aerofoil wing, wind turbine blade, etc. In particular, the invention relates to automatically designing a shape of an aerodynamic component that is optimised against a plurality of objectives, e.g. performance characteristics when the aerodynamic component is in a fluid flow.

BACKGROUND

The design of complex aerodynamic structures or components such as aircraft wings or wind turbine blades typically requires taking into account various, often conflicting, attributes and conditions. In particular, it may be desired to determine an optimal design - e.g. geometric shape/dimensions - of such aerodynamic structures with respect to a number of parameters, e.g. maximising lift, minimising drag, etc. In this way, aerodynamic structure design may be regarded as a multi-objective optimisation problem. Such optimisation will also typically be performed with respect to a number of constraints, e.g. maximum allowable structure weight, maximum/minimum dimensions in various directions, etc.

Furthermore, it is typically the case when designing aerodynamic structures that explicit forms of one or more of the objectives to be optimised are not known. This means that standard optimisation techniques, that may rely on knowledge of expressions of the objectives as a function of various variables, e.g. derivative-based techniques, cannot be used. Even if an explicit functional expression is known for a particular objective, evaluation of the objective may be time consuming to an unduly onerous extent. For instance, evaluation of the objective may involve executing a complex computational / numerical simulation, or performing a time-consuming physical experimentation. Many optimisation techniques rely on evaluating objectives at a relatively large number of points, thereby limiting their practical effectiveness for expensive objectives.

The design of aerodynamic components often makes use of Computational Fluid Dynamics (CFD) simulations. Such simulations tend to be numerically complex and computationally expensive. As such, performing CFD simulations for lots of different designs (e.g. shapes) of an aerodynamic structure to determine an optimal design - or near optimal design - by trial and error is generally not feasible.

It is against this background to which the present invention is set.

SUMMARY OF THE INVENTION

According to an aspect of the invention there is provided a method of designing an aerodynamic component. The method comprises computer-implemented steps. The method comprises: defining a plurality of variables whose values represent a geometric shape of the aerodynamic component to be designed; defining a plurality of objectives for the aerodynamic component to be designed, wherein each objective is a desired value or desired range of values of a performance characteristic of the aerodynamic component; defining a plurality of constraints to be satisfied by the aerodynamic component to be designed, the plurality of constraints including a constraint on values that each of the respective plurality of variables can take; and, determining a plurality of initial sets of values of the variables, each initial set representing a different geometric shape of the aerodynamic component.

The method comprises, for each of the plurality of initial sets of values of the variables: executing a computational fluid dynamics (CFD) solver to simulate fluid flow around the aerodynamic component of geometric shape defined by the respective initial set of variable values to obtain values of one or more CFD solver outputs indicative of performance characteristics of the aerodynamic components in the fluid flow; and, determining a value of each of the plurality of objectives based on the obtained one or more CFD solver outputs to obtain a pair of the respective initial set of values of the variables and a corresponding set of values of the objectives. The method comprises constructing a surrogate model describing a relationship between the plurality of variables and the plurality of objectives, the surrogate model being constructed based on the plurality of pairs of initial sets of variable and objective values; and, solving the surrogate model subject to the defined plurality of constraints to determine Pareto optimal solutions of the surrogate model.

The method comprises iteratively performing steps of: determining one or more further sets of values of the variables each representing a different geometric shape of the aerodynamic component; for each of the plurality of further sets of values of the variables: executing the CFD solver to simulate fluid flow around the aerodynamic component of geometric shape defined by the respective further set of variable values to obtain values of the one or more CFD solver outputs; determining a value of each of the plurality of objectives based on the obtained one or more CFD solver outputs to obtain a pair of the respective further set of values of the variables and a corresponding set of values of the objectives; updating the constructed surrogate model to be based on the plurality of pairs of initial sets of variable and objective values and the one or more pairs of further sets of variable and objective values; and, solving the updated surrogate model subject to the defined plurality of constraints to determine updated Pareto optimal solutions of the updated surrogate model, until a defined stopping criterion has been satisfied, wherein values of the plurality of variables of the designed aerodynamic component are selected to correspond to values of the plurality of variables of a solution on an approximated Pareto front of the updated surrogate model, the Pareto front being approximated based on the determined updated Pareto optimal solutions.

The aerodynamic component may be an aerofoil. Optionally, the aerofoil may be a wing.

The plurality of objectives may include one or more of: maximising lift exhibited by the aerodynamic component in a defined fluid flow; minimising drag exhibited by the aerodynamic component in the defined fluid flow; and, minimising mass of the aerodynamic component.

The plurality of variables may include a position of one or more control points on a surface of the aerodynamic component.

Each control point may be at a respective defined distance along a chord line of the aerodynamic component. The value of each control point may be indicative of the position of the control point in a vertical direction. Optionally, the vertical direction may be perpendicular to the chord line.

The plurality of variables may include one or more control points on a suction side of the aerodynamic component and one or more control points on a pressure side of the aerodynamic component.

Determining the plurality of initial sets of values of the variables may comprise performing Latin Hypercube Sampling of a decision space defined to include each combination of values of the plurality of variables that satisfies the plurality of constraints. Determining the plurality of initial sets of values of the variables may comprise performing aspect ratio scaled Latin Hypercube Sampling of the decision space.

The CFD solver may be an OpenFOAM simulation software tool.

The method may comprise defining one or more fixed parameter values for input to the CFD solver, representative of one or more initial conditions and/or boundary conditions. The CFD solver may be executed based on the defined fixed parameter values.

The method may comprise linking values of the variables to values of input variable parameters for input to the CFD solver, the CFD solver being executed based on the values of input variable parameters. The method may comprise linking values of output parameters obtained from the CFD solver to values of the objectives.

Constructing the surrogate model or updated surrogate model may comprise performing interpolation of the sets of variable and objective values. The interpolation may be performed by implementing a radial basis function method, implementing an automatic radial basis function method, or implementing a Gaussian process model.

Solving the surrogate model may comprise performing a scalarisation approach. The scalarisation approach may comprise a weighted-sum method or an e-constraint method.

Solving the surrogate model or updated surrogate model may comprise performing an evolutionary approach. The evolutionary approach may comprise an elitist non-dominated sorting genetic algorithm.

Solving the surrogate model or updated surrogate model may comprise performing an evolutionary approach. The evolutionary approach may comprise a multi-objective sequential quadratic programming algorithm.

Solving the surrogate model or updated surrogate model may comprise performing a Gaussian process regression method.

The defined stopping criterion may include that a defined maximum number of iterations of the iterative steps have been performed. The defined stopping criterion may include that a hypervolume improvement metric is within a defined distance of zero. Optionally, the hypervolume improvement metric may be Kalman filtered. The defined stopping criterion may include that a mutual domination metric is within a defined distance of zero.

The determination of the one or more further sets of values of the variables may be based on a determined value of a hypervolume improvement metric associated with the updated Pareto optimal solutions of the updated surrogate model.

The or each further set of values of the variables (i.e. evaluation points) may be determined to be the or each further set of values that maximises the hypervolume improvement metric.

If a determined value of a hypervolume metric associated with the updated Pareto optimal solutions of the previous iteration is below a defined threshold value, then the determination of the one or more further sets of values of the variables may be based on a stochastic search algorithm.

According to another aspect of the invention there is provided a method of manufacturing an aerodynamic component. The method comprises designing the aerodynamic component to be manufactured in accordance with the design method defined above. The method comprises manufacturing the designed aerodynamic component, comprising manufacturing the aerodynamic component to have the geometric shape defined by the values of the respective plurality of variables of the solution on the approximated Pareto front of the updated surrogate model.

Manufacturing the aerodynamic component may comprise forming a negative mould in which the designed aerodynamic component is to be formed.

According to another aspect of the invention there is provided an aerodynamic component designed in accordance with the design method defined above.

The aerodynamic component may be manufactured in accordance with the manufacturing method defined above.

BRIEF DESCRIPTION OF THE DRAWINGS Examples of the invention will now be described with reference to the accompanying drawings, in which:

Figure 1 shows the steps of a computer-implemented method in accordance with an aspect of the invention;

Figures 2(a) and 2(b) schematically illustrate how a degree of regularisation can be incorporated into an interpolation to account for noisy data in accordance with a step of the method of Figure 1 in some examples of the invention;

Figure 3 schematically illustrates an outline of an aerofoil or wing whose shape is to be designed in accordance with an example of the method of Figure 1 ;

Figure 4 shows a plot of objective values indicative of performance at different iterations of the aerofoil of Figure 3 designed according to the method of Figure 1;

Figures 5(a)-(f) show schematic plots of the shape of the aerofoil of Figure 3 corresponding each of six Pareto optimal solutions obtained at different iterations of the method of Figure 1 ;

Figure 6 shows a plot of all of the objective evaluation points at different iterations of an aerofoil whose shape is to be designed in accordance with another example of the method of Figure 1 ;

Figures 7(a)-(c) show schematic plots of the shape of the aerofoil designed in the example of Figure 6 at respective iterations of the method;

Figure 8 schematically illustrates an outline of an aerofoil or wing whose shape is to be designed in accordance with a further example of the method of Figure 1 ;

Figure 9(a) shows a mesh defined for use by a computational fluid dynamics (CFD) solver used to simulate performance of the aerofoil of Figure 8 as part of the method of Figure 1 , and Figure 9(b) shows a zoomed view of part of the mesh of Figure 9(a) as well as the aerofoil of Figure 8 to be simulated; Figures 10(a)-(f) show schematic plots of the shape of the aerofoil designed in the example of Figure 8 at respective iterations of the method; and,

Figure 11 shows a plot of values of an objective function, defined as part of the further example of Figure 8, against the number of iterations of the method of Figure 1 that have been performed.

DETAILED DESCRIPTION

Examples of the invention provide a method and system for automatically optimising the design of aerodynamic components, parts or structures with respect to multiple objectives and constraints, where evaluation of a particular design for an aerodynamic part involves executing a Computational Fluid Dynamics (CFD) simulation. In particular, the invention provides for determining an optimal design of an aerodynamic part in a manner that minimises, or significantly reduces, the number of CFD simulations that need to be performed relative to previous approaches, e.g. trial and error, to obtain a solution. The invention is therefore beneficial in that optimal designs for aerodynamic components are obtained, where this has not been achievable using previous approaches. The invention is also beneficial in that such optimal designs are achieved automatically, and in a timely manner. Further advantages associated with examples of the invention will become apparent, and are described, in the following.

There is provided an optimisation method and system that uses evaluated points of an objective function to determine which point(s) to evaluate next with the aim of optimising the function in the fewest evaluations. Specifically, in the present context an objective function evaluation involves performing or executing a (relatively expensive) CFD simulation for a particular design/shape of the aerodynamic part whose design is being optimised. The outputs from the simulation are, or can be used to calculate, the objective function evaluation at a particular point, i.e. for a particular set of parameter values defining the design/shape of the aerodynamic part. These objective function evaluations are used to approximate a model describing the objective function, the model then being used to determine the next points, i.e. next aerodynamic part / component design, to be evaluated (simulated using CFD). The selection of the next point(s) at which to evaluate the objective function, i.e. the aerodynamic part design for which to execute the CFD simulation, may be chosen with the aim of improving the model of the objective function (‘enrichment’) or with the aim of optimising the present approximation (‘convergence’). Examples of the invention may be implemented on a system, which may be in the form of any suitable computing device, for instance one or more functional units or modules implemented on one or more computer processors. Such functional units may be provided by suitable software running on any suitable computing substrate using conventional or custom processors and memory. The one or more functional units may use a common computing substrate (for example, they may run on the same server) or separate substrates, or one or both may themselves be distributed between multiple computing devices. A computer memory may store instructions for performing the methods performed by the system, and the processor(s) may execute the stored instructions to perform the method.

In accordance with the present invention, a step of a computational method to automatically design an aerodynamic component of optimal shape includes defining various aspects of the problem to be solved. The various aspects may be defined by a user. The method may involve defining or setting one or more parameters that define one or more aspects of the aerodynamic part to be designed. Such parameters may be fixed or constant irrespective of the eventual design of the aerodynamic part. For instance, in an example in which an aerodynamic component in the form of an aerofoil or wing is to be designed, a chord length of the aerofoil may be fixed, i.e. a length between a leading edge and a trailing edge of the aerofoil.

The method involves defining a set of a plurality of (control) variables whose values are used to represent a geometric shape of the aerodynamic component to be designed. An aim of the method is to determine values of the variables in the set that result in an aerodynamic part of optimal shape, i.e. that optimises various objectives (discussed below). For instance, in an example in which an aerofoil or wing is to be designed, a plurality of (control) points on the surface of the aerofoil may be defined as variables. In particular, the vertical position of each of the control points may be varied, and the specific vertical position of the control points (along with one or more fixed parameters) define the specific geometric shape of the aerofoil. As a purely illustrative example, a total of four control points on an aerofoil surface may be defined, two on a suction side of the aerofoil and two on a pressure side of the aerofoil. For instance, the two control points on each side of the aerofoil may be located at particular distances along a (defined or fixed) chord length of the aerofoil, e.g. 25% and 75% along the chord length. The determination of the specific shape of the aerofoil based on the values / positions of the control variables / points may be performed in any suitable way, e.g. by a suitably defined interpolation. For instance, the aerofoil shape / profile may be obtained by a thin plate spline radial basis function interpolation of the control points (where the x coordinate is the particular distance I position along the chord length, and the y coordinate is the determined vertical value at said position). It will be understood that any suitable number of variables (control points) may be defined.

The method involves defining a plurality of objectives to be achieved or exhibited by the aerodynamic component to be designed. The objectives relate to desired physical or performance characteristics of the aerodynamic component. In particular, the objectives may be that the designed aerodynamic component achieves or exhibits a desired value, or is within a desired range of values, of various specified physical or performance characteristics. It may be that different objectives to be achieved by the aerodynamic component conflict with one another. In an example in which an aerofoil is to be designed, one objective may be that it is desired to maximise lift exhibited by the aerofoil. Another objective may be that it is desired to minimise drag exhibited by the aerofoil. It will be understood that any suitable number of objectives to be achieved may be defined.

The method involves defining various constraints to be satisfied by the aerodynamic component to be designed. The constraints may for instance relate to physical constraints placed on the aerodynamic component. In particular, the constraints may include constraints on the values of the various variables that define the shape of the aerodynamic component. In the above-mentioned example of an aerofoil, the variables are the positions (e.g. vertical positions in a suitably defined coordinate system) of a number of (control) points on the suction and pressure surfaces of the aerofoil. In such an example, the constraints may include upper and lower (vertical) values or positions that each of the control points may be.

A computational simulation tool or solver is used to simulate performance of the aerodynamic component in appropriate fluid flow. The simulation tool is in particular a Computational Fluid Dynamics (CFD) tool or solver. For instance, in the case of an aerofoil the CFD solver may simulate airflow around an aerofoil, where the aerofoil is at a defined inclination angle relative to the direction of airflow. It will be understood that any suitable aerodynamic component and any suitable flow characteristics (type of fluid, speed of flow, laminar/turbulent flow, etc.) may be simulated using the CFD solver. The inputs to the CFD solver include inputs based on the variables and other fixed parameters that define the shape of the aerodynamic part to be simulated. Other inputs to the CFD solver include initial and boundary conditions, e.g. no-slip condition at interface between fluid flow and aerodynamic part surface, uniform flow field in the far field from the aerodynamic part, etc. The outputs of the CFD solver are, or are used to calculate, various aspects of the aerodynamic part’s performance, e.g. aerodynamic performance, including performance characteristics relating to the various objectives that it is desired for the aerodynamic part to achieve or exhibit. As mentioned above, such aerodynamic performance characteristics may include an amount of lift and/or an amount of drag exhibited by the aerodynamic component whose performance in the presence of an airflow is being simulated by the CFD solver.

In examples described herein, the CFD solver that is used is the OpenFOAM (Open- source Field Operation And Manipulation) simulation software tool. OpenFOAM is a numerical solver for determining solutions to various continuum mechanics problems, in particular CFD-based problems. OpenFOAM also includes various pre-processing and post-processing functionality. It will be understood that different CFD simulation tools I numerical solvers may be used in different examples of the invention.

A relationship between the shape of an aerodynamic component - as defined at least in part by the values of the defined variables - and the performance of the aerodynamic component- with respect to at least some of the defined objectives - is not known a priori, either by way of an explicit functional relationship or otherwise. Examples of the invention involve approximating a model that describes a relationship between the variables and objectives, where the approximated model may then be used to select variable values that correspond to an aerodynamic component that exhibits optimal performance relative to the defined objectives. Examples of the invention in particular use surrogate methods to approximate the relationship between variables and objectives, whereby the objectives - or objective function - is approximated by some surrogate function using currently known evaluations of the objectives.

An issue arises in that each execution of the CFD solver - i.e. simulating performance of a single design of the aerodynamic performance - is both computationally and time intensive. This means that it is not feasible to adopt a trial-and-error approach by modifying the design of the aerodynamic part - by adjusting the variable values - and then testing the performance - by executing the CFD solver, for an innumerate number of different designs. Examples of the invention therefore guide the design of the aerodynamic component using the approximated model. Unlike in some previous approaches, however, in which approximated models may be used to guide the design of components, e.g. mechanical components, in a manner that improves performance relative to objectives, examples of the invention do not simply guide design of aerodynamic parts of improved performance but rather of optimal performance with respect to defined objectives. This is discussed in greater detail below.

In accordance with examples of a method of the invention, for the design of a particular aerodynamic component with defined variables, constraints and objectives, an initial approximation of the model may be determined. In order to determine the initial approximation, an initial set of training data is needed. The training data includes inputs and corresponding outputs for the model. In particular, the training data includes sets of values of the plurality of defined variables and corresponding values of the defined objectives for each set of variable values.

Such sets of variable values and corresponding objective values will generally not be available a priori. As such, the CFD solver needs to be executed a plurality of times, i.e. for a plurality of different sets of variable values each corresponding to a respective shape of the aerodynamic component, to determine corresponding objective values for each of the sets of variable values. In general, it may be the case that the greater the number of sets of variable values and corresponding objective values that are available for the model approximation, the more accurate the approximated model will be. However, as it is expensive to execute the CFD solver, then it desirable to minimise the number of CFD solver executions that are performed to generate the initial or training set, while obtaining a training set that provides the most accurate approximation of the model describing the relationship between the variables and the objectives.

In accordance with examples of the invention, a step of the method may therefore involve selecting or determining the sets of variable values used as input to execute the CFD solver to generate the training or initial set. These selected or determined sets of variable values may be referred to as sampling points. One desirable feature of the sampling is that it is non-local. Otherwise, the objective function(s) of interest would likely only be well approximated in some local region of the overall region of possible variable value sets, and poorly approximated outside of this local region. Expressed differently, it may be desired that no two sampling points are too close to one other. Also, it is desired that the sampling covers the overall region or decision space to generate a good global approximation to the objective(s) of interest. As mentioned, for reasons of computational time and expense, it is desired that relatively few sampling points are needed.

In some examples, (direct) random sampling may be implemented to obtain a set of sampling points in the domain of interest, i.e. across the various possible combinations of values of the input variables. This approach benefits from being relatively simple, both conceptually and to implement. A prescribed number of points may be sampled from the domain of interest, i.e. decision space, with respect to a prescribed probability distribution. One such probability distribution is a uniform distribution, where each point in the decision space has an equal chance of being selected. In order to guard against clustering behaviour of points, which may leave relatively large parts of a decision space unsampled (and which becomes more of an issue for larger-dimensional problems), some examples may employ low-discrepancy (or quasi-random) sequences when performing sampling. A low-discrepancy sequence is a sequence of points that well distribute themselves over a domain, but maintain some stochastic behaviour. A low-discrepancy sequence may be open or closed. A low-discrepancy sequence is open if it is initially constructed on the assumption that it is desired to continue adding points while still maintaining a uniform-like distribution over the sampling domain. A low-discrepancy sequence is closed if it is a sequence of points that is optimised for a given sample size. For instance, a standard lattice I grid of points is an example of a closed low-discrepancy sequence.

Examples of the invention may implement a Latin Hypercube Sampling (LHS) technique to select the sampling points. The LHS technique guards against selected points becoming clustered. The LHS technique benefits from being relatively simple to implement / compute, and that its sampling is more uniformly distributed across a domain than direct sampling. Additional distribution criteria may be adopted when using an LHS technique in order to encourage uniform sampling with good coverage across a decision space, particularly for larger-dimensional problems. It will be understood by the skilled person how to implement an LHS approach to select sampling points from a prescribed domain of interest. Furthermore, an aspect ratio scaled LHS method may be used to select the sampling points, which is suitable for non-square decision spaces.

Examples of the invention may implement a maximin sampling technique to select the sampling points. This approach aims to ensure that a minimum distance between any two points is maximised, so as to provide sampling that uniformly covers the decision space. The maximin sampling approach may be defined as follows. Let x _j denote a series of samples of the (decision) space and denote the set of all N sample points by S _N (Ω ). It is desired to find the optimal configuration of points that forms a solution to the problem: where for x,y e (1, p. denotes the rescaled L ² norm between two points.

The maximin approach benefits from obtaining an even coverage of the decision space. To address issues relating to the potential computational cost associated with this approach, an approximation to the optimal sampling of the maximin technique may be implemented for improved performance. One approximation approach may be as follows. Randomly generate a decision vector x ₁ ∈ Ω . Randomly generate s decision vectors y _i ∈ (1, 1 = 1, ... , s and calculate the distances Then, letting Repeat until N decision vectors have been chosen.

Examples of the invention may implement a Halton sequence, which provides an efficient approach for providing good coverage of the decision space. Examples of the invention may implement Sobol sampling, which can be beneficial for larger-dimensional problems and is relatively cheap to compute. The skilled person would understand how to implement Halton sequence and Sobol sampling approaches.

The particular sampling method adopted may depend on dimensionality of the particular problem under consideration and on time constraints. Sobol sampling may be attractive as it may provide the most uniform sampling over the decision space, without being excessively more time consuming than other approaches such as direct sampling. In some examples, at least some of the sampling points may be user-defined or user-selected.

Before a surrogate model of the objective function may be determined, constructed or approximated, the objective values at each of the sampling points need to be determined. This involves executing the CFD solver at each of the selected sampling points. That is, for each sampling point, the CFD solver is executed to simulate an aerodynamic part of geometric shape defined by the values of the (input) variables corresponding to the respective sampling point, to determine the various physical or performance characteristics of said aerodynamic part in order to determine the objective values corresponding to the particular variable values.

A surrogate model may then be determined, constructed or approximated based on the determined set of decision-objective pairs S _N x f(S _N ) where S _N (Ω) = x ₁ , ...,x _N is the set of sample points and f(S _N ) = f(x ₁ ) ,..., f(x _N ) is the collection of evaluates. One option for constructing the surrogate model is via the use of trigonometric interpolation (or Fourier interpolation). In the present context of expensive-to-compute black-box functions in which relatively little training data, i.e. decision-objective pairs, is available, trigonometric interpolation may suffer from inaccuracies with relatively large variations between sampling points.

Examples of the invention may implement a Radial Basis Function (RBF) method to perform the interpolation to construct the surrogate model. The skilled person would understand how to implement an RBF method. There are a variety of different basis functions that may be chosen for this approach, such as linear functions, cubic functions, thin plate spline functions, Gaussian functions, and quadratic functions. The radial basis method is highly accurate for problems with low numbers of sample points, as is the case in the present context. Different choices of the radial basis function can produce different interpolation, and an appropriate function can be chosen to minimise the error to the objective function that it is desired to interpolate. In particular, for an unknown objective function as in the present context, cubic and thin plate spline basis functions have been found to provide good approximations.

Examples of the invention may implement an automatic RBF method to perform the interpolation to construct the surrogate model. In the automatic RBF method, the best or most appropriate radial basis function to perform the interpolation is chosen automatically, based on previous data. The automatic RBF method may be performed as follows. Consider N > 1 sampling points S = (x^ ..., x _N ) with which it is desired to construct a surrogate model. From this sample, choose a randomly generated index k e {1, ..., N} and denote by S ^x = {x1 _. ..,x _N }\{x _k } the set of sampling points not including x _k . Now, for each implemented radial basis function, i.e. the radial basis function cp is the thin plate spline, linear, cubic, etc., denote these by φi for i = 1, ...,s where s denotes the number of different basis functions. For each i = 1, ...,s, construct an RBF surrogate model based on the sample points S ^x . Then, choose the radial basis function φj to generate a surrogate model based on the full set of sample points S according to i.e. the basis function which best predicts the function of interest with less sampling points.

Examples of the invention may implement a Kriging method (or a stochastic / Gaussian process model) to perform the interpolation. Its application to multi-objective functions may be performed component wise. Heuristically, a Kriging interpolation model includes two components: a first component that is deterministic and captures the general trend of the provided data, e.g. via some regression; and, a second component that is a stochastic variable and measures the error between the simple regression model and the true function. The skilled person would understand how to implement the Kriging interpolation method. Fast and efficient implementation of approximations to the Kriging method are available when large numbers of sampling points are being used, to address the possible slowdown in execution of the Kriging method in such circumstances.

In accordance with the invention, a step of a method may then include solving the (constructed) surrogate model, i.e. solving a multi-objective problem in the present context. Solving a multi-objective problem means determining the Pareto front of solutions. As is known in the art, the Pareto front is composed of the Pareto optimal solutions of the optimisation problem under consideration, where a solution is Pareto optimal if it is nondominated. The Pareto front is typically difficult to compute, and may only be approximated via numerical methods. In general, in multi-objective optimisation it is not possible to find or enumerate all of the points on a Pareto front. Hence, when solving a multi-objective problem, it may be the aim to seek the best discrete representation of the Pareto front.

Different approaches for solving the multi-objective problem of the surrogate model may be implemented. In some examples, one or more scalarisation approaches may be utilised for this purpose. Scalarisation approaches encompass a set of methods that aim to reduce multi-objective optimisation problems to single objective optimisation problems. Often, to find a discrete representation of the Pareto front, a sequence of single objective optimisation problems are solved. Typically, these methods abide by the following rules: a solution to a corresponding single objective optimisation problem should provide a Pareto optimal solution; and, any point in the set of all Pareto optimal solutions should be attainable through solving a certain scalarisation of the multi-objective problem. To obtain a good representation of the Pareto front, a scalarised problem should be solved repeatedly by varying the parameters of the scalarisation. Various techniques known in the art have been developed for determining scalarisation parameters.

One scalarisation approach is a weighted sum method in which, given a sequence of non- negative weights w ∈R ^m such that an auxiliary scalarisation function where x e fl, may be constructed. To obtain up to N Pareto optimal solutions, the weights may be varied N times via w for each j = 1, ... ,N For each weight , it is required to solve: subject to g _j (x) ≤ 0,j = 1, ... ,p and h _l (x) - 0, l = 1, ... ,q

The solution to the single objective scalarised problem may be proven to be weakly Pareto optimal, and furthermore Pareto optimal if w _i > 0 for all i = 1, ... ,k or if the solution is unique. In the weighted sum method, any Pareto optimal solution may be found by altering the weights only if the problem is convex, i.e. the Pareto front is convex.

Another scalarisation approach is an ε-constraint method, in which one of the objective functions, say f _l of the multi-objective optimisation problems is selected to be optimised, whilst the others are converted into inequality constraints with variable bounds given by a sequence of points ε ∈ R ^m . To obtain Pareto optimal solutions, the bounds ε are varied and the following problem is solved: subject to f _i (x) ≤ for all j = 1,

I, g _j (x) < 0,j = 1, ...,p and h _l (x) = 0, 1 = 1, ..., q. The solution to this problem may be proven to always be weakly Pareto optimal, and indeed strongly Pareto optimal if and only if it solves the problem for each l = 1, where ε _j = fj(x) for j = 1, ... , m, j # I. In addition, a unique solution may be proven to be strongly Pareto optimal for any upper bounds ε _j . Unlike the weighted sum method, any Pareto optimal solution may be obtained via the e-constraint method, even if the Pareto set is not convex.

In examples of the invention, one or more evolutionary approaches may be implemented for solving the multi-objective problem of the surrogate model. In an evolutionary approach, an iteratively evolving population of solutions is used to approximate a Pareto optimal set. Such approaches benefit from not needing derivate information, meaning that they are relatively simple to implement, can be made computationally efficient via parallelisation, and can be applied to many different types of problem. In general, evolutionary algorithms (EAs), or genetic algorithms, typically begin their search with a population of solutions within some feasible region. After an initial selection, EAs enter an iterative step of updating their current population via the use of standard mutation procedures, e.g. selection, crossover, mutation and elite-preservation, to select a new population that reflects all of the best parts of the previous iteration. This may be performed in conjunction with a stochastic step so as to explore more of the search space and reduce the chances of becoming stuck on local optima. The iteration step stops once one or more pre-specified stopping criteria are met. Heuristically, EAs may be regarded as following a similar approach to that of natural selection.

Examples of the invention may implement an EA in the form of an elitist non-dominated sorting genetic algorithm (NSGA-II). The NSGA-II algorithm follows the general outline of an EA with a modified mating and survival process. In particular, it uses a so-called ‘elitist principle’, i.e. where only the elites of a given population are given the chance to be carried onto the next generation. Furthermore, to preserve diversity amongst solutions along the Pareto front, i.e. a good quality solution, it keeps an eye on crowding distance between solutions. The NSGA-II algorithm emphasises non-dominated solutions.

For each iteration t, the NSGA-II algorithm proceeds as follows. Feed into the algorithm a parent population P _t of size N and produce an offspring population Q _t , also of size N, using standard genetic operations. Combine the populations Q _t and P _t to create a new population R _t of size 2N. Perform a non-dominated sorting procedure on R _t , i.e. extract from R _t a series of ranked non-dominated fronts, named Fi. The final population is now filled via an elitist approach, i.e. the elite population (fronts Fi) from the non-dominated sorting are automatically transferred into the new population until the new population reaches a desired size (in this case, size N). For any front F _k that would take the total population size of P _t+1 above N, a crowding distance sorting of points is performed in that front, where only the points in F _k with the greatest crowding distance (i.e. those that are least crowded) are included in the final population P _t+1 until the desired size of N being reached, with the remaining points being discarded. NSGA-II exhibits strong performance and speed.

In examples of the invention, a so-called MOSQP (Multi-Objective Sequential Quadratic Programming) algorithm may be implemented to solve the surrogate model. The MOSQP algorithm beneficially addresses the low rate of convergence and low quality solutions obtained in some other approaches. The MOSQP algorithm solves a constrained multi- objective problem of the form subject to g _j (x) ≤ 0,j = 1, ...,p and h _l (x) = 0, l = 1, ... , q via a two-stage process, where denotes the decision variable space, f :Ω → R ^m is a vector-valued objective function with m individual objectives, g _j :Ω fl → R, h _l :Ω → R for j = 1, ...,p and l = 1, ...,q respectively, are real-valued inequality and equality constraint functions.

The MOSQP algorithm is an iterative algorithm with a two-stage process. Within each iteration, MOSQP performs an enrichment (to improve diversification) and convergence (to improve dominance of solutions) step. In particular, enrichment aims at computing a high quality approximation of the Pareto front by increasing the number of Pareto optimal points in its current iteration, and convergence aims at driving the previously computed points towards the true Pareto front. In each of these two stages, MOSQP performs a quadratic programming (QP) step to obtain search directions for either enrichment of the Pareto front, or its convergence towards optimality. The steps of the MOSQP algorithm are outlined in the following pseudocode.

As noted from the MOSQP algorithm, a candidate point may be one of three states at any one time: ‘not active’; ‘active: stopped’; and, ‘active: not stopped’. A ‘not active’ point is a point that is either dominated, or unreliable, and will not be considered further. An ‘active: stopped’ point is a point that is feasible for further improvement, but has been performed upon in the current iteration. An ‘active: not stopped’ point is a point that is feasible for further improvement, and has not been performed upon in the current iteration.

Also, for any σ > 0, there is defined a merit function Φ and Φ , for i = 1, ..., m in the algorithm as follows:

Φ and Φ _i are typically non-differentiable, although directional derivatives Φ '(x, σ; d) at any x and direction d 0 exist. This function combines the objective and constraints (like a barrier function).

The cleanup procedure mentioned in the MOSQP algorithm refers to a procedure in which points are ranked according to a crowding distance algorithm from the least crowded to the most crowded (ignoring boundary points of the Pareto front) and deactivate (i.e. turn into non-active points) points that are most crowded until the number of active points reaches n. The non-dominated procedure mentioned in the MOSQP algorithm refers to a procedure in which we return the points which are non-dominated within the input set. It may be desired to provide a good spread of solutions along the Pareto front with a view to providing future improvements of the surrogate problem. The MOSQP algorithm is beneficial as it demonstrates a good spread of solutions along the Pareto front relative to some other solvers.

Examples of the invention may include using a Gaussian process regression (GPR) method to solve a constrained black box problem. However, issues may arise regarding feasibility regions for solutions. Let lb _i ≤ x _i ≤ ub _i for i ∈ 1, ...,n, where n is the number of variables, f denotes the black box objective function and g ₁ , ... g _m are the black box constraints, which are of the form g _t ≤ 0. A GPR interpolates a set of sample points across the entire feasible region and can provide confidence intervals for the interpolation. The quantity gpr(x, ci) can be defined as the value of the prediction x with confidence interval ci, with gpr _f ,gpr _g1 , ...,gpr _gm being defined as the GPR of the objective function and constraints, respectively.

After sampling the initial points, it may be found that none of the points are feasible with respect to the defined black box constraints. In order to find a feasible point, the following problem may be solved: min |ci|, subject to gpr _{g j} (x, ci _j ) ≤ 0 ∀j ∈ 1, ..., m, lb _i ≤ x _i ≤ ub _i

∀i ∈ 1, ..., n, ci _min < cij < ci _max ∀j ∈ 1, ..., m, Here, the objective is to keep the magnitude of the confidence interval across each constraint to a minimum. The first constraint implies that for a given ci there must exist an x which is feasible for each gpr _g .. As such, the solution x* to this problem may be regarded as the most likely feasible point. The latter two constraints ensure that x and ci remain within their respective bounds.

When optimising a function f it is desired to obtain the global optimum, but also to understand f s behaviour across the whole feasible region. In a sense, these may be regarded as a single objective as an understanding of f’s values across the entire feasible region is needed in order that there is confidence that the optimum solution has been found. A cyclic method or approach may be defined to satisfy both of these aims, in particular according to the following algorithm, where ciList indicates the various confidence intervals to be considered and ciListLength indicates the number of different confidence intervals to be considered. In processes involving Gaussian process regression or Gutmann methods, for instance, such a cyclic strategy assists in not only optimising the function but also to get a better understanding of the function. At one end of the cycle, points will be generated which are ‘far away’ from previously sampled points. As these points will be in areas that have not previously been visited then a better (overall) understanding of the function is obtained. At the other end of the cycle, the focus is on achieving the optimal point, which may happen to be close to a previously sampled point. Such a point may not dramatically improve the interpolation, but it will give a solution near a local optimal point.

The problem (2a) is solved subject to (2b) and (2c) repeatedly until some stopping criteria has been met. In particular, the GPR of the objective function is minimised, with a given confidence interval, whilst ensuring that the GPR for the constraints, with the same confidence interval, are feasible. When ci = 0, there is no confidence interval, and so the optimal point is being determined based on the current interpolation at that point. When ci > 0, the point that is providing the smallest objective function is being determined, given that at the current iteration there is a ci% chance that the objective could produce this value and there is also a ci% that the constraints could be feasible.

By cycling through various values of ci, the two aims mentioned above are being achieved. As mentioned, when ci = 0 the current interpolations are being optimised, and thus the optimal solution is being sought. For larger values of ci, points are being explored that achieve optimal solutions with a high confidence interval. With GPRs, the confidence regions are largest between (evaluated) points that are farthest apart and/or where the gradient of the interpolation changes. Therefore, for larger values of ci areas are being explored that are either previously relatively unexplored or where the values of the objective function changes significantly. One or more performance metrics may be used to analyse the quality of the determined approximate solutions, i.e. how good the approximation to the surrogate model is. Factors of interest when analysing the quality of solutions may include the so-called closeness of the approximated Pareto front, and the distribution of points along the approximated Pareto front. For instance, a more uniform distribution of points along the Pareto front may typically be regarded as desirable. Various metrics are known in the art and will be familiar to the skilled person, including spread metrics, a hypervolume metric, etc.

The approach to solve the multi-objective problems in the present context is iterative in nature, with the aim of improving the approximation to a surrogate model as more iterations are performed. For such an approach, one or more stopping criteria for the approximation algorithm will typically be defined, with the aim being that the algorithm stops once a sufficient solution - which may be user-defined - has been obtained, without performing unnecessary further computations (iterations) that would waste time. One stopping criterion that may be implemented is that the algorithm stops after a defined maximum number of iterations have been performed. Other stopping criteria that may be implemented include the hypervolume improvement metric and/or the mutual domination metric, both of which are known in the art. Both of these metrics converge to zero as the level of improvement to the approximation reduces with the addition of more points. However, these indicators are non-monotonic, meaning that, on their own, they may not necessarily be ideal for use as stopping criteria. A Kalman filter may be used to assist with this issue, in particular to determine when a solution has stagnated around zero. Specifically, the Kalman filter can provide a global stopping criteria, which may use each of the previous solutions obtained by the algorithm when evaluating whether an indicator has stagnated around zero.

Once the surrogate model has been solved to obtain the Pareto front of the surrogate model, a determination as to which point(s) - i.e. which set(s) of values of the defined variables - to select next for evaluation (using the CFD solver) needs to be performed. The overall aim is to optimise the approximation of the Pareto front obtained from the surrogate model, while minimising the number of expensive evaluations - i.e. CFD solver simulations - that are needed. Selection of which points to evaluate - i.e. which sets of input variable values (corresponding to particular aerodynamic components) to provide as input to the CFD solver to simulate performance (i.e. obtain objective values) needs to be made carefully. In particular, points are chosen to gain as much useful information as possible to maximise improvement of the approximation of the Pareto front of the objective function. The selection of points to be evaluated may be a balance between enrichment - i.e. improving the current approximation of the surrogate model - and convergence - i.e. exploiting the current approximation of the surrogate model to optimise the objectives. It may be that earlier iterations of the procedure are more weighted towards enrichment to improve the initial approximation of the surrogate model, and subsequent or later iterations of the procedure are more weighted to convergence to try to obtain an optimal solution of the (better approximated) model. The method may include determining a value of a performance indicator or metric that assigns a real number to an approximated solution of the surrogate model, i.e. an approximated Pareto set of the surrogate model, that is indicative of the performance of the algorithm. Such a performance indicator that may be used is a hypervolume metric.

The procedure to select the next points for evaluation may follow a multistep process, and may include methods or techniques for improving the performance profile of the algorithm. Examples of the invention use a measure of hypervolume improvement as a metric for performance improvement. In particular, the point that is selected to be the next point for evaluation may be the point that maximises the increase in the size of the hypervolume (based on the surrogate model). This provides a balanced improvement metric and promotes solutions that converge to the true Pareto front of the original problem, and solutions that provide a good spread of solutions along the Pareto front. Indeed, the largest hypervolume is obtained when the entire Pareto front is obtained. Also, if a sufficient hypervolume - according to an appropriately defined measure - is not obtained over the most (suitably defined number of) recent iterations of the black-box algorithm, then a stochastic search algorithm may be used to select the next point(s) for evaluation so as not to become stuck at local optima and to encourage exploration of the decision space.

As mentioned above, objectives for a given set of variables may be evaluated using a CFD software tool, e.g. OpenFOAM. An interface between an application implementing the multi-objective optimisation and the CFD software allows for the variable values in the optimisation method to be implemented in the CFD software for executing a particular simulation, and the outputs from the simulation in the CFD software can be interpreted or determined as the objectives in the optimisation method. The objectives and constraints that may be defined for the particular aerodynamic component design problem under consideration can come from various external applications or combinations of applications. For instance, an objective may be defined or constructed such that its output depends on the output of various (possibly different) applications. For instance, a singular application may be an OpenFOAM simulation with certain specified parameters from which application outputs may be extracted that can be stated as an objective. When setting up a particular optimisation problem for designing an aerodynamic part, the variables of the optimisation problem (e.g. certain points on the pressure and suction sides of an aerofoil) are defined or specified as mentioned above, and then these are linked to application inputs. That is, the defined variables for the optimisation problem are linked to quantities or parameters that are defined as inputs for the CFD solver. Furthermore, the outputs obtained from the CFD solver (or other / additional applications) are then linked to defined objectives to be optimised in the optimisation problem. These can be updated and extracted after each iteration of the algorithm.

The iterative nature of the present invention is beneficial in that it allows an optimisation algorithm/routine and a CFD simulation work in conjunction with one another. In particular, this approach means that the number of computationally expensive CFD simulations is kept to a minimum, and that the results of the CFD simulations are validated as soon as the optimisation routine predicts the next sampling points (as the CFD solver is executed at each iteration). At each iteration, a point from the Pareto front is generated based on the (current) surrogate model, with the CFD simulation then being performed at said point. The output obtained from the CFD solver is then used as (further) training data to update I retrain the surrogate model to then find another point on the Pareto front. With this approach, the present invention advantageously updates the surrogate model multiple times with the aim of determining multiple points on the (true) Pareto front. This contrasts with some previous approaches, which rely heavily on an initial surrogate model, which is improved by a relatively high number of samples. Such previous approaches are therefore not only more time and computationally expensive - by virtue of the numerous samples required - but it also means that the optimisation routine would be less informed of an optimal solution region compared to the present invention. Expressed differently, use of a single surrogate model in previous approaches means that the surrogate model needs to be of a high standard to begin with, which in turn means that a relatively high number of expensive CFD simulations needs to be performed. In contrast, in the present invention a relatively small number of CFD simulations are performed initially, with further (small numbers of) CFD simulations being performed at each iteration in order to improve the surrogate model. This approach beneficially means that performing too many / unnecessary expensive CFD simulations is avoided, as the surrogate model is analysed at each iteration. Examples of the present invention are advantageous in that they select points that maximise the size of hypervolume (which are then used to update the surrogate model). This contrasts with some previous approaches in which angles are used to predict the Pareto front. The approach of the present invention in this regard is beneficial for several reasons. Firstly, a better representation of the objective space is obtained, meaning a more diverse set of solutions is provided. Secondly, a better trade-off between objectives is obtained, meaning that solution quality is improved. Thirdly, a better guided search is provided as this approach assists the optimisation algorithm with the search process.

The approach of the present invention is advantageous in that it can be used to optimise aerodynamic components of complex shape, including 3D shapes. This is because the CFD solver in the present approach implements ‘full’ or ‘high-fidelity’ surrogate models (with computational saving coming from maximising the hypervolume and retraining the model in an iterative manner, as described above). This contrasts with some previous approaches, in which an interpolation model is trained on points using ‘low-fidelity’ CFD simulations. While such previous approaches can reduce computational time for training, low-fidelity CFD simulation is unsuitable for complex aerodynamic component shapes, and for optimisation problems with a medium to high number of variables to be optimised (e.g. 10-40+ variables) or 3D shapes.

Figure 1 summarises the steps of a method 10 in accordance with examples of the invention. Steps of the method 10 may be computer-implemented, i.e. performed by one or more computing devices or computer processors. At step 101 of the method 10, a plurality of variables whose values represent a geometric shape of the aerodynamic component to be designed is defined. Also at step 101, a plurality of objectives for the aerodynamic component to be designed is defined. Each objective is a desired value or desired range of values (which could include simple being greater/less than a defined threshold value, or could include maximising/minimising a parameter) of a performance characteristic of the aerodynamic component, e.g. lift, drag, etc. Also at step 101, a plurality of constraints to be satisfied by the aerodynamic component to be designed is defined. These include a constraint on values that each of the respective plurality of variables can take.

At step 102, the method 10 includes determining a plurality of initial sets of values of the variables, each initial set representing a different geometric shape of the aerodynamic component. For instance, the initial sets of variable values may be selected using Latin Hypercube Sampling or aspect ratio scaled Latin Hypercube Sampling of the decision space.

At step 103, the method 100 involves, for each of the plurality of initial sets of values of the variables, executing a computational fluid dynamics (CFD) solver, e.g. OpenFOAM, to simulate fluid flow around the aerodynamic component of geometric shape defined by the respective initial set of variable values to obtain values of one or more CFD solver outputs indicative of performance characteristics of the aerodynamic components in the fluid flow. This step also includes, for each of the plurality of initial sets of values of the variables, determining a value of each of the plurality of objectives based on the obtained one or more CFD solver outputs to obtain a pair of the respective initial set of values of the variables and a corresponding set of values of the objectives.

At step 104 of the method 10, a surrogate model describing a relationship between the plurality of variables and the plurality of objectives is constructed. In particular, the surrogate model is constructed based on the plurality of pairs of initial sets of variable and objective values. For instance, a Radial Basis Function (RBF) method may be implemented to perform an interpolation to construct the surrogate model. At step 105, the method 10 involves solving the surrogate model subject to the defined plurality of constraints to determine Pareto optimal solutions of the surrogate model. This may involve using a scalarisation approach (e.g. a weighted-sum method or an e-constraint method), an evolutionary approach (e.g. an elitist non-dominated sorting genetic algorithm or a multi-objective sequential quadratic programming algorithm), or a Gaussian process regression method.

At step 106, the method 10 involves selecting or determining one or more further sets of values of the variables each representing a different geometric shape of the aerodynamic component. At step 107, the method 10 involves, for each of the plurality of further sets of values of the variables, executing the CFD solver to simulate fluid flow around the aerodynamic component of geometric shape defined by the respective further set of variable values to obtain values of the one or more CFD solver outputs. This step also involves determining, for each of the plurality of further sets of values of the variables, a value of each of the plurality of objectives based on the obtained one or more CFD solver outputs to obtain a pair of the respective further set of values of the variables and a corresponding set of values of the objectives. At step 108, the method 10 involves updating the constructed surrogate model to be based on the plurality of pairs of initial sets of variable and objective values and the one or more pairs of further sets of variable and objective values. At step 109, the method 10 involves solving the updated surrogate model subject to the defined plurality of constraints to determine updated Pareto optimal solutions of the updated surrogate model.

Steps 106 to 109 are performed iteratively until a defined stopping criterion has been satisfied. Values of the plurality of variables of the designed aerodynamic component are selected to correspond to values of the plurality of variables of a solution on an approximated Pareto front of the updated surrogate model, the Pareto front being approximated based on the determined updated Pareto optimal solutions.

Examples of the invention benefit from being able to deal with noisy data relating to the defined objectives and constraints. For non-noisy data, simulations of black box objectives and constraints may be performed at specific points to create an interpolation for each objective and constraint. These interpolations are then used to determine the next point the simulation(s) is to be performed at. With these interpolations, it is ensured that the interpolations return the results of the simulation for any simulated points. For noisy objectives and/or constraints, simulations may be performed that do not necessarily reproduce the same outputs for the same inputs. Current interpolation methods will fail as they cannot pass through two different values at the same point. Also, if there are two points relatively close to each other, then if they have opposing amounts of noise the interpolations can become less accurate. Figures 2(a) and 2(b) schematically illustrate how a degree of regularisation can be incorporated into the interpolation to account for noisy data. In particular, Figure 2(a) shows a radial basis function (RBF) interpolation 20 of observed points / data where the observed data is assumed to be exact (not noisy). On the other hand, Figure 2(b) shows an RBF interpolation 21 of the same observed data where the interpolation assumes that the observed data includes some noise. As is seen, the interpolation 20 in Figure 2(a) is much ‘bumpier’ than the interpolation 21 in Figure 2(b), and it is likely that the ‘smoother’ interpolation 21 of Figure 2(b) reflects the function more accurately. In order to regularise the interpolation - to account for assumed noise - rather than forcing the interpolation to pass through the observed data, it is instead forced to be within a certain range of the observed data. Increasing the size of this range will ‘flatten’ the interpolation, whereas decreasing the size of this range will result in a move back towards the ‘non-noisy’ interpolations. In some examples, the error can be user-defined. For instance, in some CFD simulations the lift and drag outputs oscillate, so the error would correspond to the size of the oscillations. However, in many cases the size or amount of the error is not given. In such cases, an error size or amount needs to be assigned to each data point. To do this, a multiobjective problem is solved in which the aim is to minimise the amount of error for each data point while also minimising the ‘bumpiness’ of the interpolation. A range of error values will be provided to choose from.

Minimising bumpiness may be performed according to: where x _t are the previous sample points, is a target objective value, s _n is the current interpolation, and μ _n (y) is defined as where A, A _n come from the interpolation. In particular, where and where 0 is the interpolation kernel. Also, P _ij = P _j (x _i ), where p is a polynomial used in the interpolation. Finally, where m is the order of the polynomial p.

Examples of the invention benefit from being able to handle both black box and white box objectives and constraints. That is, the optimisation methods of the invention can handle a wide variety of constrained black and white box optimisation problems that may include black and white box constraints. This includes being able to handle mixed box optimisation problems that include a mix of black and white box objectives and/or constraints. A black box function is a function whose behaviour is to be observed (and approximated) based on inputs and outputs to the function, where no explicit expression describing the behaviour of the function is known and/or where a relational expression for the function is prohibitively expensive to compute. By white box is meant that an expression for the function of interest is known and is relatively easy to compute.

Previous methods may only be capable of handling black box constraints. In particular, these may be embedded into the objective function acting as a penalty. Example methods of the invention can handle both white and black box constraints simultaneously, as appropriate, by not creating a surrogate model for the white box constraints. In a problem of mixed box type, exploration of the white box function feasible region is not performed I needed, with instead exploration of the black box function feasible region performed I focussed on.

An example implementing a method in accordance with the invention is now described. In particular, the described example relates to designing an aerofoil - i.e. the shape of the aerofoil - that is optimised according to a plurality of objectives. Specifically, the objectives include both a black box objective and a white box objective.

Figure 3 schematically illustrates an outline of an aerofoil or wing 30. The aerofoil has boundary thickness τ. A first objective of the aerofoil to be designed is that its shape is optimised to maximise a ratio between its coefficient of lift and drag. A second objective is that the mass of the aerofoil is minimised. This is a multi-objective problem with two objectives. The first objective is a black box objective, f ₁ (x) = C _l (x)/C _d (x), where C _l (x) is a coefficient of lift, C _d (x) is a coefficient of drag, and x is a decision variable representing the shape of the aerofoil. The second objective is a white box objective, f ₂ (x) = m(x), where m(x) is the mass of the aerofoil. The mass of the aerofoil is treated as a white box objective in the described example as it may be approximated via the analytic equation m(x) = P(x)Tp, where P(x) denotes the perimeter of the aerofoil induced by decision variable x, and T and p are assumed to be constant parameters depicting the thickness of the aerofoil boundary and the density of the material forming the aerofoil, respectively. As m(x) is simply a scalar multiple of P(x) in this example, then for brevity it is assumed that f ₂ (x) = P(x)

In the described example, the surrogate models for the black box objective are constructed using radial basis function interpolation. A stopping criterion of a total of 1000 (expensive) evaluations with respect to the first objective f _r was defined. Evaluations of the objective in this example were performed using OpenFOAM to simulate behaviour of the aerofoil in the presence of an airflow. An approximation to the Pareto front - i.e. the set of approximate optimal solutions for the multi-objective problem - is found. Initial points were selected using Latin Hypercube Sampling, with 200 points being selected. A multi- objective sequential quadratic programming algorithm was used for the multi-objective optimisation. The metric C _l /C _d was used as a metric through OpenFOAM’s result output, while surface area was a white box problem. The fluid setup in OpenFOAM was: density equal to 1.225, viscosity equal to 1.802e-05 and freestream velocity equal to 35. The particular CFD solver used was RAS SpalartAllmaras.

Figure 4 illustrates a plot of results 40 obtained for the described example. In particular, Figure 4 shows a plot of values of the second objective f ₂ against values of the first objective fa for the aerofoils of each of the iterations (where each point / dot represents the solution obtained at each iteration). In the described example, six Pareto optimal points I solutions were obtained, indicated as black dots and each labelled 401. Table 1 below indicates details of the approximated values of each of the points/solutions of the obtained Pareto set, including the iteration number of the solution as well as the associated values of the objectives fa and fa. Figures 5(a)-(f) show schematic plots of the shape of the aerofoil corresponding each of the respective six Pareto optimal solutions obtained in the described example.

In the described example, some of the obtained Pareto front solutions exhibit vortex shedding I boundary layer separation. To account for this behaviour, it may be desired to constrain the variation of lift over the simulation period within the CFD solver. In particular, this may be introduced into the above problem by enforcing an inequality constraint of the form

(i.e. a black box constraint) where C _l,var (x) denotes the total variation of the coefficient of lift over the CFD simulation period, and denotes a constant value that acts as an upper bound on the amount of variation allowed.

Another example implementing a method in accordance with the invention is now described. In particular, the shape of an aerofoil or wing is optimised using the multi- objective optimisation techniques described above in combination with the OpenFOAM CFD solver. The optimisation problem may be expressed as where f ₁ (x) = Cl(x)/Cd(x) and f ₂ (x) = log(ClVariation(x)), where 12 denotes the space of parameters that control the shape of the aerofoil. It is desired to minimise the logarithm of ClVariation as experiments have shown that this function can scale rapidly, and the logarithm provides a way to rescale this parameter appropriately for the multi-objective solver.

The CFD solver needs to interface with the multi-objective solver. The OpenFOAM application interface takes in the following parameters. A geometry that it is desired to simulate within OpenFOAM. In the described example, this is in the form of a ‘.str file representative of an initial shape of the aerofoil. The interface takes in a design box (with respect to the geometry), i.e. the region of the geometry which it is desired to optimise over. The interface takes in various other OpenFOAM-specific parameters, including a number of computer processors to run simulations with.

The geometry specified in the described example is indicated in the pseudocode below. As shown, this includes parameters specifying a maximum camber of the aerofoil, a camber distance, a maximum thickness of the aerofoil, naca points (i.e. the number of points on the wing spline used to make the . stl file, where a larger number of points gives a smoother wing), a wing span of the aerofoil, a position of the aerofoil in the simulation region, and a rotation angle of the aerofoil.

Given the specified geometry, it is desired to optimise the shape of the aerofoil. With respect to this geometry, a design box is constructed that encompasses the imported aerofoil. The design box is permitted to undergo free form deformation. The specific design box parameters set for the described example are indicated in the pseudocode below.

To setup the imported application, the above-defined geometry and design box parameters are provided along with other parameters needed by OpenFOAM to run. This can include the total / maximum number of iterations to be performed, fluid density, fluid viscosity, free stream velocity of air, and number of processors to perform the simulation. The specific parameters in the described example are specified in the below pseudocode.

The optimisation problem is then constructed. In the described example, no constraints are defined and so a constraints variable field is set as empty. It is desired to obtain the optimal shape of the aerofoil, and so the variables that can be altered are set to be elements of the design box (where the size of the search space has been defined upon creation of the design box).

The OpenFOAM application outputs parameters relating to the lift Cl, the drag Cd, and also their variation (in this case computed as the total variation of the values of the last 400 seconds of the simulation), CIVaviation, CdV aviation. The objectives are constructed by setting the objective outputs to be Cl) Cd and CIV aviation, respectively, by accessing the OpenFOAM application outputs.

The multi-objective method defined in the described example has a stopping criterion of 300 algorithm iterations and 50 initial points obtained by Latin Hypercube sampling. An output, e g. in the form of an output file from the multi-objective method solver, including a list of all of the evaluation inputs and outputs of the objectives is obtained. Figure 6 shows a plot 60 of f ₂ against for all of the evaluation points obtained in the described example, and in particular indicates the points of the approximated Pareto front, specifically indicated as dark / black dots and labelled 601. Pareto optimal solutions were obtained at iterations 51 , 127, 156, 181 , 200, 249, 250, 278, 291 , 299, and 336. Figure 7 shows example outputs of a computed aerofoil obtained in the described example. In particular, Figures 7(a), 7(b) and 7(c) respectively show the shape of the aerofoil obtained at iterations 51 (Cl/Cd = 12.98, log(ClVariation) = —21.81), 200 (Cl/Cd = 34.99, log(ClVariation) = —20.29), and 249 (Cl/Cd = 39.62, log(ClVariation) = -14.70).

A further example illustrating a method of the invention is described. With reference to Figure 8, in this example a stationary aerofoil or wing 80 with chord length c = 1 is placed in a uniform flow field. Figure 8 illustrates the geometry dimensions and boundary conditions for the problem. The horizontal velocity component u _x is set as u _m = 1 at the farfield, with the assumption that the flow is uniform at a distance away from the aerofoil. The viscosity and density are set as = 0.01 and p = 1, respectively, such that the Reynolds number obtained is Re = 100. The CFD simulations are performed in OpenFOAM, and the boundary conditions chosen to achieve this are as indicated in Table 2.

The mesh 90 used in OpenFOAM in the described example consists of roughly 30460 quadrilateral elements and is illustrated in Figure 9(a), with Figure 9(b) showing a zoomed or focussed view 901 of the mesh 90 with the aerofoil 80 at the centre. The objective f is set such that the ratio of a coefficient of lift C _L and a coefficient of drag C _D is minimised, i.e. minimise f = C _L /C _D . In the described example, these coefficients are given by the following expressions: where F _L (t) and F _D (t) respectively denote the aerodynamic forces associated with lift and drag. The variables considered in the optimisation are the vertical positions of two control points on either side of the aerofoil (suction and pressure) which are located at 25% and 75% of the chord length, respectively. The upper and lower bounds of these four variables are set as var ₄ → [—0.5,0.025], var ₂ → [—0.5,0.025], var ₃ → [0.025,0.5], var ₄ → [0.025,0.5].

Figures 10(a)-(f) schematically illustrate six of the iterations solved during the optimisation procedure, namely, iterations 1, 5, 27, 40, 48 and 62. It is clear that varying the position of the four control points - labelled as a first pair 1001a, b and a second pair 1002a,b - makes a significant difference to the design of the aerofoil 80. Note that the other points indicated in Figures 10(a)-(f) that define the shape of the aerofoil 80 are fixed between iterations. Iteration 62 was deemed to provide the best objective. The pressure side of the aerofoil obtained at iteration 62 is the largest, while the suction side is the lowest, thereby producing the largest C _L . This, in combination with the small drag coefficient, results in the optimal objective function. Figure 11 shows a plot 1101 of the value of the objective function f as the number of iterations increases. It may be seen that a converged solution is obtained after around iteration 50, where the objective function value is approximately 0.9.

Many modifications may be made to the described examples without departing from the scope of the appended claims.

Although examples of the invention have been described with reference to an aerodynamic component in the form of an aerofoil or wing, such as an aircraft wing, automotive vehicle rear wing I spoiler, etc., it will be understood that examples of the invention are applicable to the design of other types of aerodynamic component, e.g. wind turbine blade, helicopter blade, etc. Also, although examples of the invention have been described with reference to the design of two-dimensional aerodynamic components, it will be understood that the methods described herein are equally applicable to the design of three-dimensional aerodynamic components.