Multi-fidelity Kriging extrapolation together with CFD for the design of the cross-section of a falling lifeboat

1.1.Motivation

Automatic design procedures based on Computational Fluid Dynamics (CFD) simulations are becoming increasingly important in practical ship design [36]. For an efficient optimisation phase of such a design procedure, accurate high-dimensional surrogate modelling is essential [43]. A surrogate model provides a cheap substitute for an expensive (CFD) simulation, enabling us to feasibly approximate and explore the design space. To illustrate the effectiveness and promise of using surrogate models in naval optimisation routines, Scholcz and Veldhuis [37] show that Surrogate-Based Global optimisation (SBGO) of a tanker’s hull reduces the required computational time from two weeks to only a day when compared with a direct optimisation algorithm. This is important since Raven and Scholcz [34] state that the time available for hull shape optimisation with the goal of, for instance, wave resistance minimisation is usually only one to two weeks or even less.

However, SBGO suffers from the curse of dimensionality [27]. Here, the curse of dimensionality refers to the exponential growth of the design-space volume when the number of dimensions or design variables increases, thereby drastically increasing the number of (expensive) simulations required to cover this volume and reliably construct the surrogate model. In other words, if we increase the number of design variables we optimise over, we drastically need to increase the number of CFD simulations, which is costly or even infeasible. Simultaneously, high-dimensional design spaces naturally arise in naval ship design [37], making this curse of dimensionality a limiting problem for using SBGO in practice.

Multi-fidelity modelling can improve the efficiency of SBGO, by using the usual expensive and precise data together with less expensive but less accurate data. Results as early as that of Leifsson et al. [22] show that optimised hull designs can be obtained at over 94% lower computational cost when using multi-fidelity optimisation as compared to a direct high-fidelity optimisation algorithm. However, these results typically do not scale to more design variables.

We want to lessen the effects of the curse of dimensionality even further and enable ourselves to feasibly solve more complex problems of higher dimensionality, or to solve problems less costly and more accurately in general. Shan and Wang [38] survey strategies of solving high-dimensional expensive black-box problems. These strategies include problem decomposition, screening (e.g. Global Sensitivity Analysis), mapping (e.g. lossless projection of high-dimensional spaces onto lower-dimensional spaces; or multi-fidelity approaches that map a lower-fidelity onto a higher-fidelity model (such as considered in this work), and (adaptive) space reduction. From these, Shan and Wang [38] mark mapping as a promising approach and note that research does not take advantage of the characteristics of the underlying expensive functions. Wu and Wang [46] reiterate this conclusion, demonstrating that this research direction still holds relevance today.

In this work, we recognise that the grid convergence property of CFD solvers is currently an unused source of information that could further improve the performance of the surrogate model and the corresponding optimisation process. Grid convergence states that the simulation solution converges to the ‘true’ simulation solution as the grid is refined. CFD simulations reported by for instance Eijk and Wellens [40] show the clear presence of grid convergence. Therefore, this work explores a novel hierarchical multi-fidelity mapping method exploiting this grid convergence information to relieve the effects of the curse of dimensionality at the expensive high fidelity.

1.2.Method background

Kriging, also known as a form of Gaussian Process Regression [32], is a particularly popular surrogate modelling technique due to two reasons. First, Kriging provides the best linear unbiased prediction (BLUP) when the assumption of Gaussian correlations between data samples holds [35]. Second, Kriging provides an estimation of the mean squared error corresponding to that prediction [9].

Kriging receives its name from Krige [19,20], who used a precursor method to estimate gold concentrations based on only a few boreholes. It was further mathematically developed by Matheron [25] and has been popularised for use outside the field of geostatistics by Sacks [35]. Jones [14] and Forrester et al. [9] provide a gentle and intuitive introduction to (Ordinary) Kriging.

Kriging and other forms of Gaussian Process Regression are the most popular surrogate models for use in Bayesian optimisation [27], of which the Efficient Global Optimisation (EGO) algorithm [13] is one form. The EGO algorithm adaptively samples additional points based on the Expected Improvement criterion [13], which is formulated using the estimated mean squared error that Kriging provides. Using this criterion, EGO provably converges to the global optimum in an optimisation process [24], although it may exhaustively sample around local optima [8].

In contrast to Gaussian Process Regression, Kriging in principle is an interpolating formulation, which causes problems in the presence of noisy data. To solve this, Forrester et al. [7] show how to introduce regression terms to the Kriging formulation in coherence with the EGO method through a process called re-interpolation. Jones [14] discusses the relation of Kriging with respect to other surrogate modelling methods and explores possible improvements to EGO and the Expected Improvement criterion. As one of the most promising directions for further work, he hinted towards multi-fidelity simulation: using expensive high-precision simulations together with simulations of a lower precision but a much lower cost. This enables us to respectively balance exploitation and exploration of the search space and thereby cut costs while improving the quality of the surrogate.

Since then, Multi-Fidelity Kriging has received significant research effort, see for instance the review papers of Fernández-Godino et al. [10], Peherstorfer et al. [28], Viana et al. [43]. Kennedy and O’Hagan [15] mathematically formulated the fundamentals, while Forrester et al. [8] provides a more transparent overview. Gratiet and Garnier [21] show that the model of Kennedy and O’Hagan [15] can be equivalently but more cheaply formulated as recursively nested independent Kriging problems, one for each fidelity level. Perdikaris et al. [30] improve upon this work by generalising the linear autoregressive11 formulation of Gratiet and Garnier [21] to better accommodate non-linear-relations.

Korondi et al. [17], Meliani et al. [26] use the recursive formulation of Gratiet and Garnier [21] to introduce a multi-fidelity approach to the EGO algorithm and formulate level-selection criteria that balance the expected reduction of the prediction variance with the costs of sampling a given level. The combination of Gratiet and Garnier [21] together with the multi-fidelity SBGO method of Meliani et al. [26] will serve as a reference to the method proposed in this work since that combination will feature the so-called nested Design of Experiments (DoE), a sampling strategy that provides the most complete information on the relation between the fidelities [15].

Optimisation methods based on Multi-Fidelity Kriging hold promise in maritime engineering considering the following accounts. Raven [33], Raven and Scholcz [12] apply Multi-Fidelity Kriging in a ship hull form optimisation with a dense set of potential flow results and a coarse set of RANS data, and positively compared the result with the single-fidelity formulation. Korondi et al. [18] performs reliability-based optimisation of a ducted propeller using an EGO-like strategy and a multi-fidelity Kriging surrogate. Although performed under compressible flow conditions, the same techniques are interesting to incompressible flows. Liu et al. [23] uses viscous and potential flow calculations as the high- and low-fidelity respectively and shows the multi-fidelity Kriging approach can obtain a more optimal hydrodynamic hull form than the single-fidelity model alone.

Other than in these accounts, to the authors’ knowledge, Multi-Fidelity Kriging does not seem to have been applied within maritime engineering, presumably because the increased complexity of the computational setup versus the perceived potential benefit. There are still unresolved questions regarding when it is advantageous to use a multi-fidelity over the standard single-fidelity approach [11,17,39] in terms of accuracy, cost, and assumed data characteristics.

1.3.Objective

The main research objective that is investigated is how to leverage the grid convergence property of CFD solvers in order to reduce the number of expensive high-fidelity samples required for an accurate (multi-fidelity) Kriging surrogate. Using this property can accelerate surrogate-based global design optimisation in maritime engineering. CFD with a free surface description is used to computed the accelerations for stylized bow shapes of a falling life boat. The main assumption in reaching the objective is that value differences of the objective function that minimized the accelerations, implied by the grid convergence property of a CFD solver, are proportional from one point in the design space to another. To our knowledge the assumption has not been investigated in the context of Kriging or Gaussian Processes for the purpose of Bayesian Optimisation.

2.1.Basic formulation

We want to predict an unknown higher fidelity level based on as few samples as possible at that unknown level. To do this, we will estimate the training points of the Kriging model at an unknown level l+1, initially based upon only one actual high-fidelity sample at a location of interest to our optimisation. We retrieve an initial sample zl+1,x0 with x0=Xl[argmin(Zl)] and intuitively determine the multiplicative difference:

2.2.Incorporating noise and variance

Noise and variance are in principle present in each of the values used in Eq. (1).

Note again that sampling noise of numerical origin can be regressed and estimated by introducing an additional hyperparameter. We can then describe the noise using a Gaussian variance that is fully integrated with the Kriging model prediction variance.

To project this combination of noise and Kriging variance onto our estimation of the high-fidelity, we should thus incorporate them into Eq. (1). Noticing each z in Eq. (1) is actually the mean of a Gaussian, we can alternatively write the corresponding Gaussian as:

Let it be clear that we use each zp,x∗ as a training point for the Kriging model at level l+1. Effectively the variances σp,x∗ should be used as regression parameters. They can be used similarly to how we regress noise (see Forrester et al. [7]) by correcting them with the process variance σˆ. However, since the addition to the correlation matrix is not uniform (heteroscedastic instead of homoscedastic noise), we require an adapted Kriging implementation to accommodate this.

2.3.Sample weighing

Equation (1) retrieves exact predictions (predicted value equals the true value) if:

Of course, the above conditions in general will not be met exactly. Additionally, we might not know they are not met based on a single high-fidelity sample at location x0.

To detect such events, or simply (locally) improve our response surface, we must be able to retrieve and integrate additional high-fidelity samples. The predictions provided by each high-fidelity sample should therefore be weighed to obtain a weighted prediction zp,x∗. Adapting Eq. (8a) and (8b) to express this:

We will determine the coefficients ci using two metrics:

Before defining ci,distance, let us consider why the tuned kernel of the medium-fidelity, instead of the high-fidelity, is used and can be used. First, we know that these hyperparameters correctly encode the dimensional scales because these are the same for each fidelity. Furthermore, it is reasonable to assume that the medium-fidelity kernel correctly encodes the variability of the objective function, relative per dimension. Would we on the contrary tune the kernel using only the high-fidelity samples, in number as little as one, we generally have insufficient information to do so effectively and retrieve inaccurate if not nonsensical results. Therefore, using hyperparameters tuned using the medium-fidelity is a justified and reasonable option.

However, although capturing the dimensional scales and relative variability of the objective function per dimension, the medium-fidelity kernel does not necessarily provide an interpolating solution at the high-fidelity samples in the extrapolated surrogate. This is however an absolute requirement.

Therefore, we apply the elementary linear operations of row addition and row multiplication to the correlation matrix to retrieve an interpolating solution while retaining the correlation information. In this form, negative corrected ‘correlations’ indicate that a sample should provide no contribution to an extrapolated point, and are set to zero. Lastly, the columns (with the column consisting of each high-fidelity sample’s contribution to an extrapolated point) of this matrix are rescaled such that the sum equals one. The full procedure that defines ci,distance is described by Algorithm 1.

Algorithm 1

Distance weighing: scaling the correlation matrix

Note that by weighing the contribution of the high-fidelity samples in this way, we implicitly construct a difference model like Zd(x) of the multi-fidelity Kriging model. In effect, using this distance-weighing method we can partly resolve inherent non-linearities of the objective function by means of an increased sampling density in a similar fashion to regular multi-fidelity Kriging.

To define ci,variance, consider that if the variance Sf,i of Ef,i in Eq. (9a) and (9b) is large, then the fraction Ef,i upon which we base our prediction is not reliable and thus our extrapolation based on location xi is not reliable.

A clear example of when such high variances Sf,i occur is when all levels of some sampled location are extremely close to each other (and noise is involved). Small changes then potentially lead to large fraction differences or even changes in sign and so the resulting predictions are spurious and unreliable.

To resolve this, we define ci,variance using Eq. (11), which weighs the contribution of each xi to an extrapolated point using the negative exponential of the extrapolation variances σi,x∗, as defined by Eq. (8b), normalised with the minimum mean σi,x∗‾, where the mean is taken over all extrapolation locations x∗ and the minimum over corresponding samples xi:

Please note that the Gaussian-inspired use of a negative exponential is a functioning solution that fits the desired properties. However, other approaches might work as well and the usage of some scaling hyperparameter or constant could improve performance. These considerations are however left to subsequent work.

2.4.Method weighing: ‘Hybrid’ multi-fidelity surrogate modelling

When there is little convergence or difference between levels, even with little noise the prediction variance Sf,i of Ef,i in Eq. (9a) and (9b) might be large and thus our prediction unreliable. This is a fail-mode of the proposed method independent of non-linear relations between levels and thus should be presented with a separate solution. As it happens, since the differences between levels in such a case are small, intuition tells us we can regard for instance the medium fidelity to be indicative of the truth and thus of importance to the optimisation process. In other words: if a prediction based on sample xi is unreliable, that sample’s weighted contribution to the overall truth prediction as formulated in section 2.3 should be comprised of a different, more reliable solution.

To this end, we introduce prediction weight-factor wi and the ‘reliable solution’ metrics zˆl,x∗, σˆl,x∗ into Eq. (9a) and (9b):

Before defining wi in Eq. (12a) and (12b), we elaborate upon the difference between zˆl,x∗, σˆl,x∗ and zl,x∗, σl,x∗ as used previously. Suppose a scenario where a sample of the medium fidelity l lies between or close to samples of the high-fidelity l+1. Then even when w=0, it would not be reasonable to discard all weighted information and just take the value of the medium fidelity level. In fact, we would run into problems such as a perceived high noise in the high-fidelity Kriging model when samples of levels l and l+1 are very proximate as usually encountered during an optimisation process. This means our intuition as described at the start of this section would only be of benefit when using the initial DoE where samples are at a distance from each other by definition.

To solve this problem, we must weigh the contributions of the two fidelities automatically and optimally; which is the very objective of any multi-fidelity surrogate data fusion model. Therefore, we should use the Kriging mean and MSE estimate of a multi-fidelity method such as that of Gratiet and Garnier [21] to define zˆl,x∗ and σˆl,x∗, respectively.22

When using the medium-fidelity as an alternative surrogate, the proposed method gains the characteristic of being a hierarchical method. Giselle Fernández-Godino et al. [11] define a multi-fidelity hierarchical method in optimisation as a method that optimizes using the low-fidelity and refines the found extrema using the high-fidelity.

It is important to notice that Eq. (9a) and (9b) and the method of Gratiet and Garnier [21] are both in principle interpolating and exact formulations at the locations of high-fidelity samples, and thus the formulation of Eq. (12a) and (12b) will be too. The medium fidelity, however, inherently is not interpolating.

If the variance Sf,xi of Ef,xi in Eq. (9a) and (9b) is large, then the fraction Ef,xi upon which we base our prediction is not reliable and thus our prediction is not reliable. We do not want to weigh samples relative to each other but define an absolute measure of the reliability per sample. This absolute measure is captured in wi, which puts weight on either the prediction of the proposed method or that of a multi-fidelity method such as that of Gratiet and Garnier [21].

Consequentially, when Sf,i becomes higher with respect to Ef,i we want wi to diminish to indicate the unreliability of the extrapolation. Furthermore, a certain amount of noise is fully acceptable. To achieve this, the following methodology is used:

2.5.Level selection procedure

If the proposed method above is to be used within a multi-fidelity global optimisation routine, a method is required to select the level at which we will retrieve the sample at a given location x∗. Korondi et al. [17] and Meliani et al. [26] introduce such level selection algorithms, based on the work of Gratiet and Garnier [21]. Since the proposed multi-fidelity methodology is not formulated recursively, it requires a different approach.

Equation (12b) gives the total variance at x∗ of the prediction level, including method weighing. However, let us focus on only the variances of the proposed method, given by Eq. (9b). We notice that Ef,i and Sf,i are pre-calculated constants and rewrite Eq. (9b) to:

Although the above variances are calculated point-wise and do not take neighbouring points into account as the Kriging prediction variances do, this difference is of no consequence since we only require an indication of the reducible variance. In fact, the methods of Korondi et al. [17] and Meliani et al. [26] do not correspond with the Kriged variances either, nor do they take irreducible model noise into account. The presence of the second requirement, therefore, is enough to provide us with a reliable level selection methodology.

The level selection procedure is completed by selecting the level that provides the maximum variance reduction σred,x∗ at level l+1 scaled by the cost to sample [17,26]:

2.6.Alternative EI reinterpolation procedure

The Expected Improvement (EI) criterion is used by the (Multi-fidelity) Efficient Global Optimisation algorithm to select the next best location to sample to explore the design space. In the presence of noise, we would reinterpolate the surrogate [7], to ensure the EI converges to zero at sampled locations.

In our case however, we are in principle additionally involved with the extrapolated points of the proposed method. The associated extrapolation variances express the uncertainty of the extrapolation that should ideally remain informative even after reinterpolation. The reinterpolation routine of Forrester et al. [7] would render this information unused. Therefore, this section introduces a reinterpolation strategy that can accommodate both information sources and has the added benefit of being a less costly procedure by not requiring additional training steps. The procedure is given by Algorithm 2.

Algorithm 2

Alternative implicit EI reinterpolation procedure, fit for the proposed method

Algorithm 2 ensures that the next sampling location does not automatically default to the maximum EI location, but occasionally samples the location of the best extrapolation instead.

It does so either when the EI at the extrapolation location is higher than the maximum EI, or the difference between these two options normalised with the involved maximum expected cost ratio is lower than the criterion ϵEI.

At the best extrapolation, the lower fidelities were already sampled and so the high fidelity will be sampled next. The used cost ratio is therefore calculated as the cost of the high fidelity over the minimal cost of sampling at the location of maximum EI, which in this case is the cost of (nested) sampling the lower fidelities.

Algorithm 2 ensures that we do not sample the lower fidelities at infinitesimally small intervals, but sometimes sample the higher fidelity to lower the actual non-corrected Expected Improvement.

Compared to the original reinterpolation routine of Forrester et al. [7] the surrogate does not have to be retrained. The routine of Algorithm 2 is thereby about half the computational cost.

Algorithm 2 is in effect near identical to the original reinterpolation routine of Forrester et al. [7] if there are no extrapolations present. The Expected Improvement of samples is pushed to a maximum of 0. The primary difference is believed to be found in the evaluation of the Kriging variance, between samples, since the variance at the samples influences the process variance that is solved for during the hyperparameter optimisation. However, in that regard, the procedure proposed here should be regarded as the more accurate method that is more true to the Kriging model.

2.7.Implementation details

The proposed method is dependent on accurate noise estimates to provide results as desired. This observation involves some considerations in the method implementation:

Code for the proposed method is available at https://github.com/robertwenink/MFK-Extrapolation.

3.1.Case description

The in-house CFD solver that has been used is called EVA [40–42] and developed with a focus on extreme fluid-structure interaction (FSI). EVA is a multiphase flow solver based on the 2D Navier–Stokes equations for the motion of a mixed fluid. The method uses a fixed staggered Cartesian grid on which the free-surface between water and air is described using a geometric Volume-of-Fluid (VoF) method. The interface is represented by piecewise linear line segments and transported using a direction-split scheme.

Examples of extreme FSI are found when waves impact with structures [2,29,42,44,45]. Note that when waves are concerned, simulations benefit from Generating Absorbing Boundary Conditions [5,6], both in terms of accuracy and of computing time. Slamming of an object onto a free water surface constitutes another type of impact simulation, demonstrated by means of a falling cylinder on water in van der Eijk and Wellens [40].

A well-known 2D test case for such impact simulations is the ‘free-falling wedge’ which in our case mimics the collision with water of a free-fall lifeboat: a rescue vessel that drops off an offshore drilling platform into the sea in case of emergency. After the rescue vessel slams the water, constituting the moment with the highest vertical accelerations, flow separation and the generation of a cavity follows. EVA [40] has been validated for the whole entry procedure.

The deformation capacity of the falling object itself, as in Bos and Wellens [3], Xu and Wellens [47,48], is not taken into account, and the object is considered completely impermeable and not porous [44]. We will modify the scenario with a 2D falling wedge to represent a single-objective optimisation problem in which we vary the shape of the falling life boat’s nose or ‘bottom’ in the way that we have set up the simulations.

The optimisation objective is to minimize the maximum acceleration of the wedge as the main indicator for the immediate safety of its hypothetical passengers:

Fixed parameters in this problem will be:

The mass will be varied between different optimisations (but fixed per optimisation problem) since it is interesting to observe its effect. For instance, a vessel with low mass and thus low specific weight compared to the water might benefit from a more axe-shaped bottom, because buoyancy or added mass effects might de-accelerate the vessel before the primary impact. On the contrary, a vessel with a high mass and thus high specific weight would slow down much less and would benefit from avoiding anything close to flat-plate impacts.

Furthermore, from a design perspective, it makes more sense to fix the mass than for instance to fix the mass-buoyancy ratio or specific weight of the vessel. This is because the main dimensions, required load capacity, and structural weight are more or less already fixed.

We will consider the following two masses in these experiments, where one should note that the 2D CFD solver EVA accepts a mass per meter:

Since we have fixed the mass per design optimization in our experimental setup, the optimization objective acceleration is proportional to the surface-integrated body forces in the vertical y-direction.

3.2.NURBS parametrisation

The parameters which determine the shape of the wedge’s bottom spline are the parameters that define the domain of the surrogate model. We will use Non-Uniform Rational B-Splines (NURBS) to define this shape, see Piegl and Tiller [31] for extensive information. NURBS are implemented in the geomdl python package [1].

We impose some constraints on the shape endpoints:

Both horizontal 10-degree angle restrictions are in place to prevent situations too close to flat-plate impacts, which are numerically hard to simulate. Further, we consider an angle of 90 degrees vertically up at the keel to be unrealistic and not physically constructible and likewise restrict it to (90 −10) = 80 degrees.

We will limit the dimensionality of the design problem to 2. Note that normally by using NURBS, shapes such as the flare bow [50] can be created. Likewise, if we do not impose topological constraints on the vertical location of the control points, we can arrive at catamaran-like shapes. These would provide an interesting class of shapes involving air entrapment. However, proper optimisation requires a one-to-one shape parametrisation, meaning that each input x retrieves a unique shape. Simultaneously, we are limited to 2D shape variations due to the development status of EVA [40]. In 2D, we did not succeed in attaining a one-to-one parametrisation when including the above interesting variations while retaining the simpler but nonetheless important shape variations. Especially triangular bottom shapes were commonly constructed by non-unique and completely different inputs. We choose to retain only the simpler and more interpretable shape variations, defined purely by the angular constraints above.

Because we are interested in the primary impact response and not necessarily in any buoyancy effects after the bilge has submerged, we decrease the freeboard to only 0.5 m to decrease the numerical domain size.

Given a NURBS shape parameterisation, a resulting body is exported to an image. From this image, an SDF (Signed Distance Function) is created, providing the shortest distances to the body boundary where negative values are inside the body and positive values are outside the body. This SDF function is then used to project a body-fraction field upon the simulation grid corresponding to the shape of the provided NURBS image. Lastly, based on this fraction field the body is reconstructed in the grid using the Piecewise Linear Interface Calculation technique.

An example result of the mapping of a NURBS parameterisation to a reconstructed body in-grid is provided in Fig. 2.

Fig. 2.

NURBS input image (a) and reconstructed output (b) for X = [06975, 09443].

3.3.Definition of the numerical grid and the fidelity levels

The definition and coarseness of the numerical grid directly influences the accuracy and costs of the simulation and thereby determine the effective fidelity level. Therefore, first the grid definition is provided, followed by a description of how that grid is refined to create the required fidelities.

The grid definition has to conform to some requirements or preferences:

The standard grid definition in the computational domain resulting from these considerations is given by Table 1.

The definition of ‘refined’ grids and fidelity levels are simple but important, so will be provided separately here.

The refined grids or grid refinements adhere to the same rules as the standard grid definition provided earlier. Of particular importance is that the refinement does not break the grid symmetry around the rigid body. The standard grid definition provided above is from hereon referred to as ‘refinement factor 1’.

A refinement factor 2 takes the standard grid definition (refinement factor 1) and multiplies the number of cells within the focussed band of cells by a factor 2. In such a way we would retrieve a number of cells within the focussed grid area as given in Table 2.

The seemingly arbitrary usage of refinement steps of 0.5 originates from the fact that the given base grid definition uses a grid coarseness at which the simulation output seemed to be with reasonable levels of noise, provided an appropriate filtering routine and a reasonable cost level. Therefore, this seemed a good medium fidelity level, which is indexed with a 1 using 0-based indexing. The three fidelity levels required by the proposed method could now be constructed using for instance the refinement factors 0.5, 1 and 4.

Given this information, the refinement for the lower fidelity was already fixed too, at l0=0.5, due to it being the only refinement below refinement level 1 that ensures a symmetrical grid definition. Although having quite a coarse grid definition, this implies we could better expect convergence as desired while simultaneously providing us with challenging, noisy data.

The refinements used to define the fidelities for all experiments are:

These refinement levels, based upon the constant low mass scenario, are consequently applied regardless of what might be best for the other experimental scenarios. The chosen refinement level of the high fidelity is primarily based on the computation time of a CFD run. A CFD run at the lowest fidelity is completed after 1 to 3 minutes, while for the medium fidelity, this is 5 to 9 minutes. At high fidelity, a CFD run takes between 16 and 50 minutes,44 with the majority taking around 45 minutes each. CFD simulations at any refinement level higher than 2 would easily surpass a simulation time of an hour, which is deemed pragmatically infeasible during this exploratory phase of research.

3.4.CFD solver signal filtering

Having described the procedures of mapping the NURBS shape definitions onto the computational grid, we identify the need of filtering either the force signal or, equivalently, the acceleration signal retrieved by EVA such that we can retrieve a more stable objective value with a smaller error. Unfortunately, the force signal retrieved from EVA contains spikes, and these are of a purely numerical origin. Sources of these peaks are for instance the pressure peaks as described in Kleefsman et al. [16], or the PLIC reconstruction.

This custom composite filter consists of:

The results of this filtering procedure are visualized by means of a force signal in Fig. 3. In the figure, a comparison is made between the original signal and some other filters.

Fig. 3.

Comparison of filters. Zoom of the signal for the best design of the high mass scenario with m=15,816 [kg/m]. The filtered result is more stable and not biased towards the force spikes that have a purely numerical origin like the other filtering methods are. By disregarding these spikes, the filtered result more accurately follows the actual trend of the data.

3.5.Considerations of the simulation setup

Results will be compared to the surrogate modelling method of Gratiet and Garnier [21], extended with the level selection procedure of Meliani et al. [26] in order to be used in MF-SBGO. As stated above, these methods will jointly be referred to as ‘the reference method’. The reference method has a low-fidelity level and a high-fidelity level. Note that the method proposed here, compared to the reference method, uses an additional fidelity level. A lower fidelity level was added so that that becomes the low-fidelity level in our proposed method. The original fidelity levels in the reference method will hence be referred to as medium-fidelity and high-fidelity levels in the proposed method. The other main difference with the reference method is that the proposed method features the assumption of uniform proportional grid convergence throughout the design space.

To provide a comprehensive comparison, the following metrics will be employed to offer insights into both the surrogate and SBGO methodology in terms of the costs, performance, and other underlying factors:

Performance is measured by comparing the predicted value with the real sampled results by using the Root Mean Squared Error (RMSE) as done in Toal [39]:

Random number generators are used in many aspects of the code:

4.1.Optimisation metrics and evaluation of the main assumption

Optimisation metrics for all the cases are given by Table 3, with the corresponding best-optimised designs reported by Fig. 8. The table compares the outcomes for the proposed optimisation method with the reference method [21,26] for the two main scenarios with different masses. Variations of the methods are made by using a different number, either 1 or 3, of starting high-fidelity samples.

For both the high mass (m=15,816 [kg/m]) and the low mass (m=1,255 [kg/m]) scenarios, the initial DoE has been fully sampled at each fidelity and represented in Fig. 4 and Fig. 5. These figures also show the design that is near-optimal. For each of the designs, the force as a function of time is provided at the three different fidelity levels so that it can be investigated how the designs interact with the water, and what the effect of grid refinement is. Although the optimisations themselves do not have access to the full high-fidelity level, these plots provide context to the optimised designs and insight into the overall convergence behaviour of EVA and the problem at hand.

Fig. 4.

Convergence for the falling wedge experiment with m=15,816 [kg/m] on the initial DoE of size 20 and over the refinement levels 0.5, 1 and 2. The black line plots the optimisation objective over the fidelity levels while the coloured lines provide the full CFD simulation time traces. The design for each NURBS parameterisation x is shown in the background. The design at x = [0.5982, 0.0332] is near-optimal.

The proposed method relies on the convergence of the objective value over fidelity levels – as stated in the main assumption. This convergence was assumed to be uniform across the design space. The validity of this assumption is inspected before analysing the optimisation results. When we look at Fig. 4 and Fig. 5 it is seen that there is the anticipated convergence at badly performing designs (large acceleration), but almost no convergence at well-performing designs (low acceleration). The observation is that CFD simulations near the optimum design will encounter less strenuous circumstances with lower forces, pressures, fluid velocities and accelerations present. Because the flow variables for these designs have lower values overall near the optimum designs, good estimates of these flow variables are already achieved at lower grid resolution. Therefore, we observe less convergence of the objective value for better designs.

Fig. 5.

Convergence for the falling wedge experiment with m=1,255 [kg/m] on the initial DoE of size 20 and over the refinement levels 0.5, 1 and 2. The black line plots the optimisation objective over the fidelity levels while the coloured lines provide the full CFD simulation time traces. The design for each NURBS parameterisation x is shown in the background. The design at x = [0.8809, 0.5214] turns out to be ‘optimal’, i.e. causes the algorithm to stop because the stopping criteria have been met.

4.2.High mass case and the low mass case

The end-states of the surrogate of the proposed method and the surrogate of the reference method [21,26] are plotted and visually compared to the ‘truth’ surrogate that is based on a fully sampled DoE at the high-fidelity in Fig. 6 for the high mass case with m=15,816 [kg/m]. Figure 7 shows the same for the low mass case with m=1,255 [kg/m].

Fig. 6.

Visualization at the optimisation end of the high mass case with m=15,816 [kg/m]. On the left: the proposed method with 3 starting high-fidelity samples + MFK optimisation that was able to find a better design than all other setups. On the right: the reference method, which ends with a high NRMSE due to a wrongly build trend model for the additive differences. More high-fidelity samples in the initial DoE might decrease the issue.

Fig. 7.

Visualization at the optimisation end of the low mass case with m=1,255 [kg/m]. On the left: the proposed method with 1 starting high-fidelity sample + MFK optimisation that has all the best cost statistics and the best ending NRMSE. On the right: the reference method, which like the high mass case ends with a high NRMSE.

Figure 6 reveals the following for the high mass case. In testing the method weighing procedure, the medium fidelity is not used if the MFK solution is available, which is from 3 starting high fidelity samples and on. The problem with using the medium fidelity as the alternative surrogate is that it is not principally a solution interpolating with the high-fidelity. Although the high fidelity is forcibly added and there is in this case typically barely any difference between the high and medium fidelity in areas where the proposed method does not work well and so method weighing is used, differences over short distances could be relatively large and thus noise estimates overestimated.

For the case where instead the surrogate modelling method of Gratiet and Garnier [21] provides the alternative surrogate, we are ending with a higher NRMSE, similar to the experiments with the reference method [21,26], even though starting with a very accurate surrogate. When comparing the time traces of the NRMSE together with the weighing variables wi over time, we observe that as we sample more high-fidelity samples, there is a heightened probability of hitting ‘no convergence’ areas and thus likely falling back to the alternative method. At the same time, the reference surrogate provides increasingly bad results when sampling the area near the optimum. This was not expected given a Pearson correlation of 0.990 which normally implies the reference is effective [39].

Continuing on this, the simulation start of the proposed method with MFK as the alternative surrogate knows a very low NRMSE. This coincides with the reference method having good accuracy at the start too. The proposed method already applies method weighing at the simulation start, and it seems we retrieve a result that takes the best of both worlds. On top of that, this is the only optimisation run that retrieves a better design. This shows the idea of method weighing can be beneficial.

Fig. 8.

Best results of the wedge optimisations with (a) constant low mass (b) constant high mass. Next to the found optimal shape, the figures provide the body forces in the y-direction for the original acceleration signals and the signals after filtering. Although we minimise the maximum acceleration, the body forces in the y-direction are given because 1) constant body masses are used and so the translation to acceleration is the same for each point in the design space, and 2) this provides insight into how the body mass influences the involved forces whereas Tab. 3 does not.

Figure 8a shows the force on the optimal high mass case design as a function of time, and the optimal design itself. Contrary to expectations, the bottom end of the design is somewhat blunt, starting with a low gradient near the nose, after which the gradient of the cross-section increases before it becomes lower near the sides. Considering the hydrodynamics of water jet formation and jet seperation from the sides of the wedge, that is associated with the force on the wedge becoming lower, this design makes perfect sense: the design tries to avoid a single high force peak by postponing jet separation for as long as possible with a more constant force over time as the outcome.

The body forces of the best-optimised low mass design, given by Fig. 8b, are seen to be stable around the maximum value for most of the impact time, which indicates this is a very good design where we cannot expect to gain more by reducing force peaks. The design maintains a constant deadrise for most of the height of the bottom before tapering off to a more horizontal end at the chine. The peak of the low mass case force is nearly an order of magnitude lower than the high mass case force peak.

While interpreting the objective values of the low mass scenario, one should consider that due to the 2D nature of EVA used the wedge or ‘lifeboat’ effectively drops onto the body of water over its full length at once. This results in much higher forces and accelerations compared to the high mass scenario of m=15,816 [kg/m]. Again, as with the high-mass case, there is low to no convergence around the optimum design.

4.3.Cost-efficiency

The objective of the proposed method is also to lower the costs of building an accurate surrogate and optimise using it, by lowering the required number of samples at the expensive high fidelity. To assess if successful, Table 3 considers the following main metrics:

For the sampling cost, it should be taken into account that the ‘low’ fidelity in Table 3 consists of only the medium fidelity for the reference method [21,26] and both the low- and medium-fidelity for the proposed method. Next to this, to improve the proposed method’s stability, extra samples around the high-fidelity samples are acquired to improve the initial noise estimates per fidelity level. Therefore, any gains by the proposed method at the high fidelity need to offset the sampling costs of the complete lowest fidelity to perform better than the reference method. The proposed method is able to do this for the simulations starting with a single high-fidelity sample.

In terms of surrogate accuracy expressed in the NRMSE, the proposed method gives comparable results to the reference method [21,26], but cannot be said to perform better. It is found to be more stable during the course of EGO optimisation. Given the nature of the problem, the reference seems to need more high-fidelity samples in the initial surrogate to perform well.

Because of the instability of the reference method and the corresponding upward parabolic surrogate, the reference might decide earlier than the proposed method to stop the optimisation based on the Expected Improvement criterion. The reference method might therefore seem to be more efficient in costs than it deserves in these scenarios.

As for the optimisation result, only the proposed method has been able to find a better design, for the high-mass case. More investigation is required to determine whether this is due to a structural strength of the method.