Figure 2 shows the numerical solutions for the free-surface elevation at probe 5 (the focus location), along with the physical measurement, for the empty tank test 2BT2 (the intermediate steepness case). Also plotted is the variance density spectra of the results,

${\psi }_{\eta }\left(f\right)$

, computed over an interrogation window of 35.3–50.3 s. The linearised methods (LPT, WECSim 1 and 2 and non-linear Froude–Krylov) have not been included because these are assumed to have ‘perfect’ reproduction of the surface elevation at this position (the ‘target’ position).

Figure 2. Surface elevation at the target position, probe 5, from the empty-tank test of the mid-steepness wave, 2BT2: (a) time series, and (b) variance density spectra (computed over an interrogation window of 35.3–50.3 s). Physical measurements are plotted using a black dotted line; numerical submissions, from all participants, are shown using coloured lines

As the participants had access to the empty tank data, and the situation is greatly simplified, convergence of the calculated surface elevation is expected (at least between similar methods). In general, the agreement between the numerical models and the experimental data is good, and there is less variation across the numerical solutions compared with the series 3 results (Ransley et al., 2020) (although the participants are different and some of those with the least good results in series 3 have not participated in series 2); however, some of the results again demonstrate that there can still be significant discrepancies, even in this ‘visible’ case, depending on the implementation strategy. When viewed in frequency space, the same observations can be made, however; the frequency-domain analysis appears to have filtered out the obvious discrepancies in the OpenFOAM (waves2Foam) result, which might only be high-frequency distortions, and has highlighted some curious issues underlying the OpenFOAM (source-term) result, which could be a consequence of the wave generation strategy used. In addition to this, it has been noticed that the quality of the reproduction varies with respect to the position interrogated and at other probe locations there is more spread in the results. Ignoring the linearised methods, which demonstrate significant discrepancies away from the target location (due to the non-linear propagation of these waves not being described well by linear wave theory), the variations observed away from the target location highlight a clear issue with modern wave generation techniques; many of the methods that appear to perform the best use an iterative method to ‘tune’ their wave generation to achieve the desired surface elevation at the target location; however, if the result is not equally as good at other positions, can one be sure that the wave, including its underlying kinematics, has truly been reproduced as well as it appears to have been at the target location? Further, these methods tend to require the iterative procedure to be performed for every new case, which can increase the required computational resource considerably, particularly when using an already computational expensive method like CFD, and renders each simulation somewhat more difficult to reproduce by others. This issue is somewhat analogous to the ‘transfer functions’ used to control physical wavemakers which tend to be bespoke adjustments to the wavemaker displacements to achieve the desired result (typically only surface elevation is considered) at a particular position in the physical basin. This practice, in general, renders physical experiments as ‘facility-specific’ and, analogously, iterative wave tuning renders numerical simulations as ‘method-specific’ (albeit somewhat more transparent compared to the physical analogy). These observations, and the anticipated correlation between the quality of incident wave reproduction and that of the structural response, demonstrate a clear need for standardisation in numerical wave generation practices and also represent an additional complexity in judging the predictive capability of numerical models. There is no obvious trend in the quality of the reproductions as a function of wave steepness.

Figure 3 shows the predicted heave displacement of geometry 1's CoM when subject to the intermediate steepness wave, 2BT2, as well as the corresponding variance density spectra. In general, the predicted heave displacement is reasonably good across all model fidelities. Possibly because, in these cases, heave motion is dominated by inertia and restoring (rather than hydrodynamic) forces; however, there is a noticeable spread in the results (±10% of maximum heave) with the LPT model over-predicting and the WEC-Sim and Froude–Krylov methods under-predicting the amplitude of heave displacement; the NS solvers typically perform slightly better, so there is some evidence to suggest the quality of the reproduction is a function of model fidelity (i.e., some non-linear effects might be present). Despite this, as was the case in Series 3 (Ransley et al., 2020), the quality of the heave displacement prediction appears to correlate well with that of the surface elevation reproduction (for the non-linear methods which model the wave propagation). In frequency space, the non-linear/propagation models again have the same trends as in the surface elevation; the LPT method has an excellent main peak but significantly overestimates the higher frequency content; the Froude–Krylov method also has noticeable disparities at high frequencies (as well as underestimating the main peak); and there is considerable variation between the two WEC-Sim results (particularly at high frequencies), perhaps suggesting that the estimate of the viscous drag and non-linear forcing terms is important. In general, there is no obvious difference in the quality of the heave prediction with respect to the two geometries, although, for Geometry 2, the WEC-Sim 2 implementation appears to have some issues with premature oscillation of the structure early in the time series. Again, there is no obvious trend in the quality of the reproduction as a function of wave steepness.

Figure 3. Heave displacement of geometry 1's CoM when subject to the mid-steepness wave, 2BT2: (a) time series, and (b) variance density spectra (computed over an interrogation window of 35.3–50.3 s). Physical measurements are plotted using a black dotted line; numerical submissions, from all participants, are shown using coloured lines

Figure 4 shows the predicted surge displacement of Geometry 1's CoM when subject to the intermediate steepness wave, 2BT2, along with the corresponding variance density spectra. As in series 3, there is considerably more spread in the surge results compared to heave, particularly after the main wave crest has passed and the structure displays significant ‘drift’ motion. Again, in general, the NS solvers and hybrid methods perform best and, out of these, those which reproduced the incident wave the best also appear to reproduce the surge motion the best; however, the two NS methods utilising overset meshing (STARCCM+ and OpenFOAM (overset)) do consistently predict premature surge motion that is not seen in the other methods or the physical data. Considering the linearised methods, the Froude–Krylov method does predict drift motion but has considerable discrepancies (possible due to assumed linear wave dispersion or exclusion of a viscous drag estimate); as expected, the pure LPT method does not predict any drift motion because no non-linear effects are included; there is a clear difference between the two WEC-Sim results: WEC-Sim 2, which does not consider non-linear forcing and has a more approximate strategy to estimate the viscous drag (Hughes et al., 2020), is missing the drift motion completely; whereas, WEC-Sim 1, which includes weakly non-linear restoring and Froude–Krylov forcing (based on the instantaneous body position and wave elevation) and utilises a NS solver to estimate the viscous drag coefficients (van Rij et al., 2020), gives a remarkably good prediction for the drift motion, comparable with the least good NS solvers. Despite the quality of the WEC-Sim 1 prediction, all of the linearised methods used here neglect forces from the second-order difference frequency in the spectrum (available in the NREL FAST software (Duarte et al., 2014)) as well as those from radiation and drag damping, which were derived in the frequency domain by Roald et al. (2013), for example, and shown to be important in predicting the mean drift of floating structures (Stansby et al., 2019). These observations suggest that the inclusion of non-linear forcing terms, particularly those associated with the instantaneous body position (and, perhaps, more accurate estimates of viscous drag) are necessary to predict the surge motion in these focused wave cases. In frequency space, the variation in the prediction of low-frequency drift motion is evident and the correlation with the incident wave prediction appears to be lost at low frequencies. Again, there is no obvious difference in the quality of the reproduction with respect to the two geometries; however, qualitatively, there is an apparent reduction in the predictive capability as a function of wave steepness (particularly for the linearised models but also generally).

Figure 4. Surge displacement of geometry 1's CoM when subject to the mid-steepness wave, 2BT2: (a) time series; and (b) variance density spectra (computed over an interrogation window of 35.3–50.3 s). Physical measurements are plotted using a black dotted line; numerical submissions, from all participants, are shown using coloured lines

As in the CCP-WSI Blind Test Series 3 (Ransley et al., 2020), the mooring load in these cases is dominated by vertical motion – that is, heave. Consequently, the quality of the reproductions is similar to that of surface elevation and heave response (with some influence from surge motion) and there is no obvious difference in the quality of the reproduction as a function of either the wave steepness or the structure's geometric complexity. Curiously, all the models predict a greater peak in the variance density of the mooring load relative to the experiment, compared to that for the heave displacement.

Dynamic time warping (DTW) is a method for measuring the similarity between two data signals that vary in time. DTW was originally designed for automatic speech recognition (Sakoe and Chiba, 1978) and, in contrast to the linear alignment of signals by cross-correlation, it is used to find the optimal global alignment between two time series by non-linearly ‘warping’ the time axis of the signals until their dissimilarity is minimised. This process effectively removes, from the similarity measure, the influence of both ‘shifts’ in the time dimension – that is, phase discrepancies, and the ‘speeds’ of the two time series – that is, frequency discrepancies, offering, what is believed to be, an independent measure of the amplitude discrepancy between two time series. In addition, a ‘warping path’ is generated containing the information on how to translate, compress and expand the signals so that similar features are matched. DTW has been used in many fields, from analysis of electrocardiogram (ECG) measurements (Huang and Kinsner, 2002) to fault detection in waste-water treatment (Jun, 2011) and sewer flow monitoring (Dürrenmatt et al., 2013), but, despite this, to the authors’ knowledge, DTW has not been used to quantify the predictive capability of numerical models, particularly in the field of offshore and coastal engineering.

In this study, the Matlab function, dtw, is used. Figure 7 shows an example of the results from a DTW analysis of one of the blind heave submissions (a qualitatively accurate one in this case). Compared to the original data, the ‘warped’ data have excellent alignment in time, with only small differences now present due only to amplitude discrepancies. Also given in Figure 7, is the minimised Euclidean distance parameter between the warped experimental and numerical data, defined as

$$d\left(\mathcal{X},\mathcal{Y}\right)=\sum _{w_k}{D}_{n,m}$$

1

where

$D\in {\mathbb{R}}^{N\times M}$

is the distance matrix computed between all points in the two sequences,

$\mathcal{X}=\left(x_1,x_2,\dots ,x_n,\dots ,x_N\right)$

and

$\mathcal{Y}=\left(y_1,y_2,\dots ,y_n,\dots ,y_M\right)$

, given

$$D_{n,m}=\sqrt{{(x_n-y_m)}^2}$$

2

and

$W=\left(w_1,w_2,\dots ,w_k,\dots ,w_K\right)$

is the warping path, computed as a sequence of consecutive matrix elements defining a mapping between

$X$

and

$Y$

with the

$k$

th element being

$w_k=(n,m{)}_k$

, that minimises

$d$

(Dürrenmatt et al., 2013).

The warp path, for the example in Figure 7, shows noticeable compression (gradients

$<1$

), and expansion (gradients

$>1$

) is needed at the beginning and the end of the interrogation window, to align the numerical data with the experimental time series; but during the middle period the alignment is near perfect (gradient

$=1$

). The deviations, from perfect alignment, at the beginning of the warp path might be due to low amplitudes (in heave displacement in this case), and small Euclidean distances, which can lead to unrealistic warping paths when the input signals are noisy or erroneous (Dürrenmatt et al., 2013). The deviations towards the end of the warp path might indicate the moment when inconsistencies in the reflected waves contaminate the signal. In addition, the total length and the gradients of the warp path may offer measures of similarity that are independent of discrepancies in amplitude and, therefore, an opportunity to uncover trends, similar to those observed for heave, in the prediction of other degrees of freedom like pitch and surge.

Figure 7. An example of DTW analysis showing: (a) the original experimental heave displacement and one of the numerical submissions, from case 2BT2 (geometry 1), along with their warped signals and the minimised Euclidean distance between them; and (b) the warp path compared to perfect alignment

Figure 8 shows the Euclidean distances between the submitted numerical data and the physical measurements. The distances are calculated over an interrogation window from 35.3 to 50.3 s, from data interpolated onto a fixed sample frequency of 128 Hz, and are normalised by the standard deviation of the physical data over the same period (to achieve some independence with respect to the window size (Ransley et al., 2020)). Pitch results from the LPT method have been excluded due to an unresolved issue with the submitted data. Figure 8(a) shows that, when using this measure, the capability of the NS/hybrid codes, to predict the amplitudes of heave motion, correlates strongly with their ability to reproduce the incident wave amplitudes. This reinforces suggestions that the heave, in these cases, is dominated by inertia and restoring (rather than hydrodynamic) forces. Compared to NRMS (Figure 6), the use of Euclidean distance appears to separate the NS/hybrid methods into two distinct groups, and the trend in the data suggests a normalised Euclidean distance of zero is achievable, for heave displacement, provided perfect reproduction of the incident wave amplitudes. This suggests that the apparent limit in the capability to predict heave motion, found using NRMS (Ransley et al., 2020), is due to phase or frequency content discrepancies and the better performing group of NS methods have relatively higher discrepancies of this type compared to the poorer performing group. Figure 8(a) also shows that, when using Euclidean distance (as opposed to RMS), the capability of the linearised methods to predict heave motion is apparently improved, relative to the other methods, suggesting these methods suffer more from discrepancies in frequency and phase, than amplitude, when predicting heave motion. This observation demonstrates the benefit of being able to specify the wave elevation at the position of the structure precisely because the majority of the linearised predictions of the heave amplitudes are on par with the stronger group of NS/hybrid methods (and noticeably better than the weaker NS methods). Figure 8(b) shows the OpenFOAM (waves2Foam) method has a relatively poor prediction of the surge amplitudes, compared to the heave (and, therefore incident wave) amplitudes, and the OpenFOAM (source-term) method's prediction of surge amplitudes is stronger relative to heave. Other than this, in general, those NS/hybrid methods that predict the heave (and, therefore incident wave) amplitudes the best tend to predict the surge amplitudes the best (but the trend is far less pronounced compared to the heave–wave relationship). For the linearised methods, the spread in the predicted surge amplitudes is large: WEC-Sim 2 and the non-linear Froude–Krylov method have very poor predictions for the surge amplitudes of geometry 2 and geometry 1, respectively; but WEC-Sim 1 has a remarkably low Euclidean distance for surge (even noticeably lower than that of the STARCCM+ method used to tune its viscous drag coefficients). Figure 8(c) shows that the spread in the predicted pitch angle amplitudes (plotted on a log scale) is the greatest of all the degrees of freedom considered with the majority of the linearised methods giving very poor predictions, particularly for geometry 2. The WEC-Sim 1 implementation is, again, the exception with Euclidean distances for pitch comparable to the stronger NS/hybrid methods. These observations suggest that, for cases with unbroken waves interacting with simple floating systems, lower fidelity methods, like WEC-Sim, can demonstrate predictive capabilities similar to high-fidelity NS solvers, but only if improved estimates for the non-linear terms are found beforehand. For NS/hybrid methods, the quality of the predicted pitch amplitudes tends to follow that of the heave (and, therefore incident wave) amplitudes, but, like surge, the trend is less pronounced compared with the relationship between heave prediction and incident wave reproduction. Curiously, considering its apparently superior performance in other degrees of freedom, the PIC method has abnormally poor predictions of the pitch amplitudes for cases with geometry 2; the OpenFOAM (source term) and STARCCM + methods also have greater pitch Euclidean distances for geometry 2, compared to geometry 1, but the OpenFOAM (waves2Foam) method has the opposite trend. Last, Figure 8(d) shows the Euclidean distances for pitch against those for surge; clearly (with the exception of WEC-SIM 1), the linearised methods are the weakest at predicting the amplitudes of pitch and surge motion; for the rest of the methods, there does not appear to be a clear trend between the Euclidean distances in pitch and surge, but, considering both degrees of freedom together, the PIC (geometry 1), OpenFOAM (source term) (geometry 1) and hybrid methods (along with the WEC-Sim 1 implementation) appear superior to the rest of the NS solutions, which seem to be limited by a compromise between a strong performance in either pitch or surge (but not both).

Figure 8. Normalised Euclidean distance (relative to the experimental data) for: (a) heave displacement against surface elevation at wave probe 5 (during empty tank tests); (b) heave against surge displacement; (c) pitch angle against heave displacement; and (d) pitch angle against surge displacement. Data colour coded by underlying theory/method (red = NS solvers, green = LPT, cyan = hybrid, magenta = PIC); filled markers = Geometry 1; open markers = Geometry 2; marker sizes scaled according to the wave steepness,

$kA$

As mentioned earlier, selecting the right numerical model for any particular application is considered to be a compromise between the model fidelity and the required computational effort to gain a solution. Figure 9 shows the normalised Euclidean distance for both heave and surge displacement against the CPU effort required to generate the solutions. CPU effort has been defined as the execution time of the numerical solver (in seconds), multiplied by both the number of cores used and their speed in GHz then divided by the simulated time in seconds. For methods requiring ‘preprocessing’ steps to complete the simulations – for example, tuning of the incident waves, the CPU effort required for these has been added when these steps must be repeated for every wave case (when the preprocessing only needs to be performed once per geometry for example, calculation of hydrodynamic coefficients, this effort has not been included). It should be noted that participants were not asked to minimise the CPU effort as part of this study. The WECSim 1 results have been plotted excluding the CPU effort required to estimate the viscous drag coefficients which, for the actual submitted results, required the STARCCM + results (for each case individually) to have been found prior. This technically means the WECSim 1 CPU effort is the greatest (six orders of magnitude higher than that plotted); however, this has been considered to be misleading because typically the viscous drag coefficients for each geometry would be estimated more economically (rather than for each incident wave case individually). Following the blind test, the use of a single set of drag coefficients (for each geometry) was found to have a fairly minimal effect on the results (van Rij et al., 2020).

Figure 9. Normalised Euclidean distance (relative to the experimental data) for: (a) heave and (b) surge displacement against CPU effort. Data colour coded by underlying theory/method (red = NS solvers, green = LPT, cyan = hybrid, magenta = PIC); filled markers = Geometry 1; open markers =  Geometry 2; marker sizes scaled according to the wave steepness,

$kA$

. NB: WECSim 1 plotted in black as preprocessing effort has been excluded (technical CPU effort (including ‘atypical’ preprocessing)

$\sim {10}^8$

)

As expected, the methods based on linear theory require significantly less CPU effort than the NS solvers (up to eight orders of magnitude less) and for heave displacement, at least, there is no significant improvement in the Euclidean distance to warrant the additional computing cost (for the cases considered here). Considering the NS solvers, the two hybrid methods do demonstrate significant savings in CPU effort, with the exception of the OpenFOAM (source term) method, which exploits the symmetry of the test cases and is therefore noticeably more efficient than the other pure NS solvers. As noted in series 1 (Ransley et al., 2019), there is a large range in the required CPU effort for the NS solvers. It is suspected that this is mainly due to the specific implementations – that is, mesh/domain design, applied by the individual operators (rather than the efficiency of the underlying code) and, although participants were not asked to minimise the CPU effort, this highlights another opportunity for best-practice procedures to provide value.