Highly Accurate Experimental Heave Decay Tests with a Floating Sphere: A Public Benchmark Dataset for Model Validation of Fluid–Structure Interaction

: Highly accurate and precise heave decay tests on a sphere with a diameter of 300 mm were completed in a meticulously designed test setup in the wave basin in the Ocean and Coastal Engineering Laboratory at Aalborg University, Denmark. The tests were dedicated to providing a rigorous benchmark dataset for numerical model validation. The sphere was ballasted to half submergence, thereby ﬂoating with the waterline at the equator when at rest in calm water. Heave decay tests were conducted, wherein the sphere was held stationary and dropped from three drop heights: a small drop height, which can be considered a linear case, a moderately nonlinear case, and a highly nonlinear case with a drop height from a position where the whole sphere was initially above the water. The precision of the heave decay time series was calculated from random and systematic standard uncertainties. At a 95% conﬁdence level, uncertainties were found to be very low—on average only about 0.3% of the respective drop heights. Physical parameters of the test setup and associated uncertainties were quantiﬁed. A test case was formulated that closely represents the physical tests, enabling the reader to do his/her own numerical tests. The paper includes a comparison of the physical test results to the results from several independent numerical models based on linear potential ﬂow, fully nonlinear potential ﬂow, and the Reynolds-averaged Navier–Stokes (RANS) equations. A high correlation between physical and numerical test results is shown. The physical test results are very suitable for numerical model validation and are public as a benchmark dataset.


Introduction
Numerical models with complex fluid-structure interactions are often developed to simulate motions of floating bodies in the ocean, which can be applied to assess the performances of wave energy devices; see, e.g., [1,2]. Despite the complexity of such models, the discretization and assumptions needed to formulate the numerical model mathematically inevitably introduce errors, for many of which the influences are unknown. Engineers may struggle to identify whether linear wave theory can be applied with sufficient accuracy or more advanced computational fluid dynamics (CFD) methods should be used. Physical tests of high accuracy and reproducibility are paramount for validation and calibration purposes when using such advanced methods; see, e.g., [3,4].
The International Energy Agency Technology Collaboration Programme for Ocean Energy Systems (OES) has initiated the OES Wave Energy Converters Modelling Verification and Validation working group (formerly OES Task 10). Here, multiple research institutions and R&D companies from 12 countries collaborate with the focus on the development of numerical models for simulating wave energy converters (WECs) [5]. A floating sphere was chosen as a practical representation of a simple wave energy convertor buoy, and numerical modelling of the decay of a sphere was completed as an initial test case [6][7][8]. The resulting simulations from the different members showed widespread simulation results, which highlighted the need for knowing the true, real-world results for the considered test case together with the associated measurement uncertainties. In order to validate and calibrate numerical models, a high-quality benchmark dataset was needed. Such datasets were lacking, so during a Danish-granted EUDP project [9] a sphere model was built, and tests were performed in the wave basin in the Ocean and Coastal Engineering Laboratory at Aalborg University in Denmark. The test design, namely, the release mechanism and the construction of the sphere, was optimized through several stages to mitigate sources of uncertainties. A 300 mm diameter aluminum sphere model with changeable ballasts-see Figure 1-was chosen as the most practical and accurate representation of a sphere for physical heave decay tests dedicated to producing a highly accurate benchmark dataset. The benchmark dataset is publicly available in the supplementary material of the present paper; see Appendix A. The iterations in the design and construction process of the physical test setup are described in [10], which is also included as supplementary material in the Descriptions folder. In [10], the tests are referred to as the Kramer Sphere Cases. A new test case was formulated to accurately represent the performed tests and allow for numerical replications for model validation against the benchmark. Three different drop heights were investigated. The aim of the present study was to estimate the precision and accuracy of the physical decay tests using uncertainty analysis and comparisons to state-of-the-art hydrodynamic numerical models for all three drop heights. Using this approach, the applicability of the benchmark dataset to validation of numerical modelling of the presented test case is accounted for. The presented uncertainty analysis is based on the ASME Performance Test Code Test Uncertainty [11], which is in accordance with the A new test case was formulated to accurately represent the performed tests and allow for numerical replications for model validation against the benchmark. Three different drop heights were investigated. The aim of the present study was to estimate the precision and accuracy of the physical decay tests using uncertainty analysis and comparisons to state-of-the-art hydrodynamic numerical models for all three drop heights. Using this approach, the applicability of the benchmark dataset to validation of numerical modelling of the presented test case is accounted for. The presented uncertainty analysis is based on the ASME Performance Test Code Test Uncertainty [11], which is in accordance with Energies 2021, 14, 269 3 of 36 the methodologies and nomenclature of the ISO/IEC Guide 98-3 Guide to Expression of Uncertainty in Measurement (GUM) [12], but contains a more technical treatment.
In Section 1.1, the test case is presented. All physical parameters are given to mimic the setup of the conducted heave decay tests. The reader can set up his/her own numerical model based on the information given herein, and thereafter apply the generated benchmark dataset for comparison/validation. Dedicated measurements of certain physical parameters, such as air pressure and viscosity, are not included in the test case. These are instead considered in the uncertainty analysis in Section 3.
The test case was given to participants of the OES working group, who independently formulated numerical models to simulate the decay tests utilizing miscellaneous modelling approaches. In the order of descending fidelity, these models included finite volume method (FVM) 3D unsteady Reynolds-averaged Navier-Stokes (URANS) models, boundary element method (BEM) fully nonlinear potential flow (FNPF) models, and BEM linear potential flow (LPF) models. The utilized numerical modelling approaches are presented in Section 1.2.

The Test Case
Consider an ideal sphere with a diameter D and a mass m. In a local Cartesian coordinate system, the origin coincides with the geometrical center of the sphere and with the z-axis vertical oriented upwards. The center of gravity is CoG. The local acceleration due to gravity is g.
The sphere floats between an air and a water phase, when at rest (equilibrium). The water phase has the density ρ w , while the density of air is disregarded. A fixed global Cartesian coordinate system is defined from the still water level; the xy-plane coincides with the plane of the free water surface, and the z-axis is vertical oriented upwards towards the air phase; see Figure 2. The sphere is half-submerged when at rest, and with the CoG on the z-axis (underneath the center of buoyancy), the local and global coordinate system axes will coincide when the sphere is at rest; see Figure 2. The seabed is horizontal with a depth of d = 3D.
Energies 2020, 13, x FOR PEER REVIEW 3 of 36 methodologies and nomenclature of the ISO/IEC Guide 98-3 Guide to Expression of Uncertainty in Measurement (GUM) [12], but contains a more technical treatment. In Section 1.1, the test case is presented. All physical parameters are given to mimic the setup of the conducted heave decay tests. The reader can set up his/her own numerical model based on the information given herein, and thereafter apply the generated benchmark dataset for comparison/validation. Dedicated measurements of certain physical parameters, such as air pressure and viscosity, are not included in the test case. These are instead considered in the uncertainty analysis in Section 3.
The test case was given to participants of the OES working group, who independently formulated numerical models to simulate the decay tests utilizing miscellaneous modelling approaches. In the order of descending fidelity, these models included finite volume method (FVM) 3D unsteady Reynolds-averaged Navier-Stokes (URANS) models, boundary element method (BEM) fully nonlinear potential flow (FNPF) models, and BEM linear potential flow (LPF) models. The utilized numerical modelling approaches are presented in Section 1.2.

The Test Case
Consider an ideal sphere with a diameter and a mass . In a local Cartesian coordinate system, the origin coincides with the geometrical center of the sphere and with the -axis vertical oriented upwards. The center of gravity is . The local acceleration due to gravity is .
The sphere floats between an air and a water phase, when at rest (equilibrium). The water phase has the density , while the density of air is disregarded. A fixed global Cartesian coordinate system is defined from the still water level; the xy-plane coincides with the plane of the free water surface, and the z-axis is vertical oriented upwards towards the air phase; see Figure 2. The sphere is half-submerged when at rest, and with the on the z-axis (underneath the center of buoyancy), the local and global coordinate system axes will coincide when the sphere is at rest; see Figure 2. The seabed is horizontal with a depth of = 3 . Initial conditions of zero velocity and zero acceleration are applied in all test setups. Under the assumption of a rigid body, the sphere has six degrees of freedom (DoF). Translations relative to the rest condition in the directions of the local , , and -axes are defined as surge , sway , and heave , respectively. Rotations relative to the rest condition around the local , , and -axes are defined as roll , pitch , and yaw , respectively. Three initial test setups are investigated with displacements of the sphere in positive heave given by the drop height = 0.1 , 0.3 , 0.5 ; see Figure 3. Initial conditions of zero velocity and zero acceleration are applied in all test setups. Under the assumption of a rigid body, the sphere has six degrees of freedom (DoF). Translations relative to the rest condition in the directions of the local x, y, and z-axes are defined as surge x 1 , sway x 2 , and heave x 3 , respectively. Rotations relative to the rest condition around the local x, y, and z-axes are defined as roll x 4 , pitch x 5 , and yaw x 6 , respectively. Three initial test setups are investigated with displacements of the sphere in positive heave given by the drop height H 0 = {0.1D, 0.3D, 0.5D}; see Figure 3.  The sphere is released, and around eight natural periods in heave should be captured for comparison to the benchmark dataset. The physical parameters of the test case are presented in Table 1. The utilized initial conditions match those of previous tests carried out under the OES working group.

Numerical Modelling Blind Tests of the Test Case
Participants of the OES working group independently developed numerical models to simulate the test case presented in Section 1.1 and to compare results against the benchmark. Only the governing physical parameters of the test case, given in Table 1, were shared with the participants, and the numerical modelling of the test case was thus carried out in blind without any shared information on domain geometry, resolution, turbulence modelling, etc. Various types of numerical models were developed by the participants. The specifications of the numerical model developed by each participant are presented in Appendix B. In general, three categories of numerical models were used: (i) FVM-based Reynolds-averaged Navier-Stokes (RANS) models, (ii) BEM-based fully nonlinear potential flow (FNPF) models, and (iii) BEM-based LPF models. These are introduced in the following subsections.
An analytical solution of the Navier-Stokes (N-S) equations would yield an exact model of the fluid flow of any Newtonian fluid, such as water. In their most general form, the N-S equations are the formulation of conservation of mass, momentum, and energy into a set of nonlinear partial differential equations. Currently, no analytical solutions to the N-S equations exist, but several numerical solutions have been established, introducing various simplifying assumptions and levels of inaccuracies. In general, decreasing the complexity of the mathematical problem by simplifying assumptions will yield less accurate numerical models, but increase the computational efficiency creating more feasible models. The influences of the errors introduced by the numerical model are strongly casespecific, and no generic model with a perfect balance of accuracy and efficiency is currently available.

RANS Models
Within high-fidelity CFD modelling of WECs, RANS models have become the model of choice [13]. The RANS equations are based on Reynolds decomposition and ensembleaveraging of the N-S equations. This reformulation of the N-S equations introduces a term referred to as the Reynolds stress, which accounts for the contribution of turbulent fluctuations to the fluid momentum. Turbulence structures are not resolved in RANS models, and thus computational effort is significantly decreased relative to, e.g., direct numerical simulations (DNS). Larger unsteady mean flow structures are captured from The sphere is released, and around eight natural periods in heave should be captured for comparison to the benchmark dataset. The physical parameters of the test case are presented in Table 1. The utilized initial conditions match those of previous tests carried out under the OES working group.

Numerical Modelling Blind Tests of the Test Case
Participants of the OES working group independently developed numerical models to simulate the test case presented in Section 1.1 and to compare results against the benchmark. Only the governing physical parameters of the test case, given in Table 1, were shared with the participants, and the numerical modelling of the test case was thus carried out in blind without any shared information on domain geometry, resolution, turbulence modelling, etc. Various types of numerical models were developed by the participants. The specifications of the numerical model developed by each participant are presented in Appendix B. In general, three categories of numerical models were used: (i) FVMbased Reynolds-averaged Navier-Stokes (RANS) models, (ii) BEM-based fully nonlinear potential flow (FNPF) models, and (iii) BEM-based LPF models. These are introduced in the following subsections.
An analytical solution of the Navier-Stokes (N-S) equations would yield an exact model of the fluid flow of any Newtonian fluid, such as water. In their most general form, the N-S equations are the formulation of conservation of mass, momentum, and energy into a set of nonlinear partial differential equations. Currently, no analytical solutions to the N-S equations exist, but several numerical solutions have been established, introducing various simplifying assumptions and levels of inaccuracies. In general, decreasing the complexity of the mathematical problem by simplifying assumptions will yield less accurate numerical models, but increase the computational efficiency creating more feasible models. The influences of the errors introduced by the numerical model are strongly case-specific, and no generic model with a perfect balance of accuracy and efficiency is currently available.

RANS Models
Within high-fidelity CFD modelling of WECs, RANS models have become the model of choice [13]. The RANS equations are based on Reynolds decomposition and ensembleaveraging of the N-S equations. This reformulation of the N-S equations introduces a term referred to as the Reynolds stress, which accounts for the contribution of turbulent fluctuations to the fluid momentum. Turbulence structures are not resolved in RANS models, and thus computational effort is significantly decreased relative to, e.g., direct numerical simulations (DNS). Larger unsteady mean flow structures are captured from the unsteady RANS (URANS) formulation (see, e.g., [14]), to the extent allowed by the Energies 2021, 14, 269 5 of 36 temporal resolution. In the present paper, URANS models are developed from the opensource framework of OpenFOAM (versions 5.0, 7, and v1912) [15] and the commercial code StarCCM+ 13.06 [16]. The numerical models utilize the FVM to discretize the RANS equations. The interface between the two fluid phases is tracked by a volume of fluid (VOF) advection scheme; see, e.g., [17]. The models further assume incompressible, isothermal, immiscible flows.

FNPF Models
In the FNPF category of CFD models, further assumptions for the second-order nonlinear N-S equations are made; i.e., the fluid domain is assumed inviscid and irrotational, thereby introducing potential flow theory, which reduces the governing equations of the fluid domain to Laplace's equation [18]. The boundaries of the fluid domain evolve in time, to be able to capture finite-amplitude waves and have a time-varying wetted body surface. The boundary conditions of the fluid domain are fully nonlinear in the sense that the velocity potential satisfies the nonlinear kinematic and dynamic boundary conditions at the free surface. No-flow boundary conditions are satisfied at solid boundaries [19]. In this study, the FNPF commercial code SHIPFLOW-Motions 6 [20] was applied. Here, a mixed Eulerian and Lagrangian (MEL) scheme [21] is utilized to capture the nonlinear free surface. The positions of free surface particles are then tracked in time in a Lagrangian representation of the flow problem, allowing for the advection of mesh nodes [22]. A rigid six-DoF model is included to update the position of the wetted surface at each time step.

LPF Models
At the low-fidelity end of CFD models to simulate WECs are the LPF models, which despite rather gross assumptions of linearity in both the governing equation (Laplace) and the boundary conditions, produce useful simulations for engineering purposes and indeed are very time-efficient; see, e.g., [23]. The dynamic response of marine structures is commonly analyzed in the frequency domain using LPF theory [23][24][25][26]. Time-domain models are based on hydrodynamic coefficients solved in the frequency domain and inserted into the Cummins equation [27,28]; see Appendix C for further information. In the present paper, hydrodynamic coefficients are calculated in the frequency domain from the BEM-based LPF software WAMIT [29]. Five models of various levels of accuracy are considered. The LPF0 model is based on the solution to a traditional one-DoF massspring-damper system with constant hydrodynamic coefficients; i.e., the added mass, the hydrodynamic damping, and the hydrostatic stiffness are merely evaluated at a single frequency (damped natural frequency). Furthermore, the draft-dependency is disregarded in the calculation of the hydrodynamic coefficients, in which the sphere is considered static at the neutrally buoyant position (submergence to the equator). The LPF1-4 models are based on the Cummins equation, allowing the description of arbitrary motions (multiple frequencies) rather than a regular motion (single frequency). For LPF1, the hydrodynamic coefficients in the frequency domain are calculated for the neutrally buoyant position and are assumed as linear. Various levels of nonlinearities (draft-dependencies) are added as extension of each other to LPF2, 3, and 4: Respectively, the hydrostatic stiffness, the added mass at infinite frequency, and the convolution part of the radiation force are nonlinearized. The utilized LPF models are thoroughly presented in Appendix C.

Materials and Experimental Setup
In the present section, the materials and setup of the physical heave decay tests conducted at Aalborg University are presented. Four repetitions were carried out for each drop height.

The Sphere Model
The sphere model was constructed using computer numerical control (CNC) machining of two aluminum blocks into two hemisphere shells of equal outer radii. A thread was cut internally at the equator of the sphere to be able to assemble and disassemble the two hemisphere shells; see Figures 4 and 5a. A thin rubber gasket was installed to seal the model when assembled; see Figure 5b. The sphere was designed with an adjustable internal ballast system. A thread was tapped internally at the bottom of the model to fix ballast weights; see section view A-A in Figure 4. cut internally at the equator of the sphere to be able to assemble and disassemble the two hemisphere shells; see Figures 4 and 5a. A thin rubber gasket was installed to seal the model when assembled; see Figure 5b. The sphere was designed with an adjustable internal ballast system. A thread was tapped internally at the bottom of the model to fix ballast weights; see section view A-A in Figure 4.
Additional threads were tapped externally at the top and bottom of the sphere to allow attachment of lines for decay tests and future tests, including mooring and power take-off (PTO). For line attachment to the sphere model, custom-made M8 nuts were used; see Figure 6a. In the presented tests, a line was merely mounted to the top of the sphere to displace it in the positive z-direction as the initial condition. A nut was installed at the bottom external thread with a cover of polyvinyl chloride (PVC) tape; see Figure 6b,c. The sphere model was marked with thin lines to have a reference system of and , as also seen in Figure 6b,c.  The sphere model was constructed using computer numerical control (CNC) machining of two aluminum blocks into two hemisphere shells of equal outer radii. A thread was cut internally at the equator of the sphere to be able to assemble and disassemble the two hemisphere shells; see Figures 4 and 5a. A thin rubber gasket was installed to seal the model when assembled; see Figure 5b. The sphere was designed with an adjustable internal ballast system. A thread was tapped internally at the bottom of the model to fix ballast weights; see section view A-A in Figure 4.
Additional threads were tapped externally at the top and bottom of the sphere to allow attachment of lines for decay tests and future tests, including mooring and power take-off (PTO). For line attachment to the sphere model, custom-made M8 nuts were used; see Figure 6a. In the presented tests, a line was merely mounted to the top of the sphere to displace it in the positive z-direction as the initial condition. A nut was installed at the bottom external thread with a cover of polyvinyl chloride (PVC) tape; see Figure 6b,c. The sphere model was marked with thin lines to have a reference system of and , as also seen in Figure 6b,c.  Additional threads were tapped externally at the top and bottom of the sphere to allow attachment of lines for decay tests and future tests, including mooring and power take-off (PTO). For line attachment to the sphere model, custom-made M8 nuts were used; see Figure 6a. In the presented tests, a line was merely mounted to the top of the sphere to displace it in the positive z-direction as the initial condition. A nut was installed at the bottom external thread with a cover of polyvinyl chloride (PVC) tape; see Figure 6b,c. The sphere model was marked with thin lines to have a reference system of x and y, as also seen in Figure 6b,c.
An optical 3D motion capture system was utilized to track four reflective markers installed on top of the model. In order to minimize the reflections from the model itself, the upper hemisphere shell was painted matte black. Ballast weights were CNC machined from stainless steel and mounted internally at the bottom of the lower hemisphere; see Figure 7. An optical 3D motion capture system was utilized to track four reflective markers installed on top of the model. In order to minimize the reflections from the model itself, the upper hemisphere shell was painted matte black. Ballast weights were CNC machined from stainless steel and mounted internally at the bottom of the lower hemisphere; see Figure 7. The machined components (i.e., the hemisphere shells and the ballast weights), were constructed with a precision of 0.1 mm. The dimensions of the additional components (i.e., nuts and reflective markers), were known with the same precision. The weight of each of the individual components of the sphere model was measured on precision scales with a precision of 0.1 g. A 3D computer-aided design (CAD) drawing of the sphere model was created in which densities were ascribed to the individual components from the measured weights. The total mass, total center of gravity (in the local coordinate system defined in  An optical 3D motion capture system was utilized to track four reflective markers installed on top of the model. In order to minimize the reflections from the model itself, the upper hemisphere shell was painted matte black. Ballast weights were CNC machined from stainless steel and mounted internally at the bottom of the lower hemisphere; see Figure 7. The machined components (i.e., the hemisphere shells and the ballast weights), were constructed with a precision of 0.1 mm. The dimensions of the additional components (i.e., nuts and reflective markers), were known with the same precision. The weight of each of the individual components of the sphere model was measured on precision scales with a precision of 0.1 g. A 3D computer-aided design (CAD) drawing of the sphere model was created in which densities were ascribed to the individual components from the measured weights. The total mass, total center of gravity (in the local coordinate system defined in The machined components (i.e., the hemisphere shells and the ballast weights), were constructed with a precision of 0.1 mm. The dimensions of the additional components (i.e., nuts and reflective markers), were known with the same precision. The weight of each of the individual components of the sphere model was measured on precision scales with a precision of 0.1 g. A 3D computer-aided design (CAD) drawing of the sphere model was created in which densities were ascribed to the individual components from the measured weights. The total mass, total center of gravity (in the local coordinate system defined in Section 1.1), and total moments of inertia of the sphere model installed with ballast to generate half-submergence are given in Table 2. In the Supplementary Materials under the Descriptions folder, the dimensions, weights, and centers of gravity are given for all individual components.

Experimental Setup and Equipment
The decay tests were carried out in the wave basin in the Ocean and Coastal Engineering Laboratory at Aalborg University in Denmark. The wave basin measured 13.00 × 8.44 m, and a water depth of 900 mm was used for all tests. The wave basin had vertical wavemaker pistons and vertical passive wave absorber elements installed. The wavemaker pistons were inactive during the tests. The sphere model was released in the middle of the basin; see Figure 8. A camera was mounted for documentation purposes, and three wave gauges were installed to measure the radiated waves from the decays and reflected waves; see Figures 8 and 9. Wave gauge data were collected, partly to assess reflections, and partly to analyze radiated waves in further work. The position of the sphere model was tracked by a Qualisys Motion Capture System; four Oqus7+ cameras at 300 fps with invisible, infrared strobes were mounted in the air phase, pointing towards the model; see Figures 9 and 10.
Energies 2020, 13, x FOR PEER REVIEW 8 of 36 Section 1.1), and total moments of inertia of the sphere model installed with ballast to generate half-submergence are given in Table 2. In the supplementary material under the Descriptions folder, the dimensions, weights, and centers of gravity are given for all individual components.

Experimental Setup and Equipment
The decay tests were carried out in the wave basin in the Ocean and Coastal Engineering Laboratory at Aalborg University in Denmark. The wave basin measured 13.00 × 8.44 m, and a water depth of 900 mm was used for all tests. The wave basin had vertical wavemaker pistons and vertical passive wave absorber elements installed. The wavemaker pistons were inactive during the tests. The sphere model was released in the middle of the basin; see Figure 8. A camera was mounted for documentation purposes, and three wave gauges were installed to measure the radiated waves from the decays and reflected waves; see Figures 8 and 9. Wave gauge data were collected, partly to assess reflections, and partly to analyze radiated waves in further work. The position of the sphere model was tracked by a Qualisys Motion Capture System; four Oqus7+ cameras at 300 fps with invisible, infrared strobes were mounted in the air phase, pointing towards the model; see   The release of the sphere model was initiated by a mechanical system consisting of a pushrod and a small electrical actuator; see Figure 11. A line was mounted to the top of the sphere model at the one end and to a small nut at the other end. The nut was supported by the pushrod preceding the initialization of the tests. A trigger signal was sent to the actuator which displaced the pushrod backwards (towards the actuator), thereby removing the support of the sphere model. The release time was measured by highspeed cameras (960 fps) to less than 1/960 s [10]. The line connecting the sphere model to the pushrod was a Suffix ® 832 line with 8 braided fibers and 32 weaves per inch (thickness 0.30 mm, weight 0.18 g/m).  The release of the sphere model was initiated by a mechanical system consisting of a pushrod and a small electrical actuator; see Figure 11. A line was mounted to the top of the sphere model at the one end and to a small nut at the other end. The nut was supported by the pushrod preceding the initialization of the tests. A trigger signal was sent to the   The release of the sphere model was initiated by a mechanical system consisting of a pushrod and a small electrical actuator; see Figure 11. A line was mounted to the top of the sphere model at the one end and to a small nut at the other end. The nut was supported by the pushrod preceding the initialization of the tests. A trigger signal was sent to the  actuator which displaced the pushrod backwards (towards the actuator), thereby removing the support of the sphere model. The release time was measured by highspeed cameras (960 fps) to less than 1/960 s [10]. The line connecting the sphere model to the pushrod was a Suffix ® 832 line with 8 braided fibers and 32 weaves per inch (thickness 0.30 mm, weight 0.18 g/m). The sphere model was displaced in positive heave to approximately match the test case drop heights , as given in Table 1. The sphere model was kept at a given drop height, until the model and the free water surface were at rest; see Figure 12. The initial calmness of the sphere model (measured drop heights, velocities, and accelerations) and the free water surface are quantified in Section 3.  The sphere model was displaced in positive heave to approximately match the test case drop heights H 0 , as given in Table 1. The sphere model was kept at a given drop height, until the model and the free water surface were at rest; see Figure 12. The initial calmness of the sphere model (measured drop heights, velocities, and accelerations) and the free water surface are quantified in Section 3.
Energies 2020, 13, x FOR PEER REVIEW 10 of 36 actuator which displaced the pushrod backwards (towards the actuator), thereby removing the support of the sphere model. The release time was measured by highspeed cameras (960 fps) to less than 1/960 s [10]. The line connecting the sphere model to the pushrod was a Suffix ® 832 line with 8 braided fibers and 32 weaves per inch (thickness 0.30 mm, weight 0.18 g/m). The sphere model was displaced in positive heave to approximately match the test case drop heights , as given in Table 1. The sphere model was kept at a given drop height, until the model and the free water surface were at rest; see Figure 12. The initial calmness of the sphere model (measured drop heights, velocities, and accelerations) and the free water surface are quantified in Section 3.

Results
The measured heave decay time series and the associated systematic and random uncertainties are accounted for in the present section. Furthermore, deviations between the ideal test case and the physical tests are quantified and considered. Heave x 3 of the sphere was measured as the displacement of the sphere in the global z-axis. The influences of rotations in roll and pitch on the heave measurements of the sphere model were included in the uncertainty analysis.

Decay Measurements and Expanded Uncertainty
The measured heave decay time series are presented for the three investigated drop heights in Figure 13. To mitigate the effect of small variations in the drop height between the repetitions, the heave decay time series were normalized with the respective measured drop heights H 0,m . Time was normalized with the damped natural period in heave T e0 = 0.76 s; see Appendix C.

Results
The measured heave decay time series and the associated systematic and random uncertainties are accounted for in the present section. Furthermore, deviations between the ideal test case and the physical tests are quantified and considered. Heave of the sphere was measured as the displacement of the sphere in the global -axis. The influences of rotations in roll and pitch on the heave measurements of the sphere model were included in the uncertainty analysis.

Decay Measurements and Expanded Uncertainty
The measured heave decay time series are presented for the three investigated drop heights in Figure 13. To mitigate the effect of small variations in the drop height between the repetitions, the heave decay time series were normalized with the respective measured drop heights , . Time was normalized with the damped natural period in heave = 0.76 s; see Appendix C. The measured heave decay time series included with 95% confidence intervals (CIs) around the sample mean are presented for each of the investigated drop heights in Figure  14. To be able to distinguish the 95% CIs from the sample mean, a zoom of the first trough is included in Figure 14. Both the normalized and raw heave decay time series can be found in larger formats in Appendix D, where the 95% CIs are upscaled to be able to visualize the time-dependency of the CIs. The 95% CIs were calculated from the Taylor series method (TSM) in accordance with the recommendations in [11]. The calculation of both the random and systematic uncertainties in the physical heave decay tests are described in the present section.
The time-dependent, two-sided 95% CI on the sample mean ( ) was established from expanding the combined standard uncertainty ( ) by the value , following the Student's distribution [30]. refers to the confidence level and is the number of degrees of freedom (not to be confused with the previously introduced rigid body motions, but rather the independent variables in the calculation of ) given by = − 1, with being the number of repetitions.
where ( ) is referred to as the expanded uncertainty, and . , = 3.182 [30]. The combined standard uncertainty ( ) was calculated as the root-sum-square of the random standard uncertainty and the systematic standard uncertainty as per TSM [11]: The random standard uncertainty of the sample mean was directly calculated from the sample standard deviations at each instant of time (ISO Type A) as The measured heave decay time series included with 95% confidence intervals (CIs) around the sample mean are presented for each of the investigated drop heights in Figure 14.
To be able to distinguish the 95% CIs from the sample mean, a zoom of the first trough is included in Figure 14. Both the normalized and raw heave decay time series can be found in larger formats in Appendix D, where the 95% CIs are upscaled to be able to visualize the time-dependency of the CIs. The 95% CIs were calculated from the Taylor series method (TSM) in accordance with the recommendations in [11]. The calculation of both the random and systematic uncertainties in the physical heave decay tests are described in the present section.
The time-dependent, two-sided 95% CI on the sample mean X(t) was established from expanding the combined standard uncertainty u X 3 (t) by the value t Cv following the Student's t distribution [30]. C refers to the confidence level and v is the number of degrees of freedom (not to be confused with the previously introduced rigid body motions, but rather the independent variables in the calculation of u X ) given by v = N − 1 with N being the number of repetitions.
where U X (t) is referred to as the expanded uncertainty, and t 0.95,3 = 3.182 [30]. The combined standard uncertainty u X (t) was calculated as the root-sum-square of the random standard uncertainty s X and the systematic standard uncertainty b X as per TSM [11]: The random standard uncertainty of the sample mean was directly calculated from the sample standard deviations at each instant of time (ISO Type A) as The systematic standard uncertainty was calculated as the root-sum-square of the elemental systematic standard deviations; see Table 3. The quantification of the elemental systematic standard deviations is accounted for later in the present section.
The systematic standard uncertainty was calculated as the root-sum-square of the elemental systematic standard deviations; see Table 3. The quantification of the elemental systematic standard deviations is accounted for later in the present section. Figure 14. Normalized heave decay time series and 95% CIs with a zoom of the first trough. Figure 14. Normalized heave decay time series and 95% CIs with a zoom of the first trough. By multiplying the expanded uncertainty time series U X (t) for each drop height with the respective averaged measured drop heights H 0,m , each expanded uncertainty (with a confidence level of 95%) was given a physical dimension (length in mm); see Figure 15.  the respective averaged measured drop heights , , each expanded uncertainty (with a confidence level of 95%) was given a physical dimension (length in mm); see Figure 15. The precision of the motion capture system (incl. calibration) was assessed from displacements of the sphere model in heave with high-precision blocks 50.0 mm in height. By comparing position time series, the systematic standard uncertainty of the motion capture system setup was found to be 0.01 mm (ISO Type A); refer to [10] for further information. The systematic standard uncertainty introduced by vibrations of the bridge (reference frame for the motion capture system) after release of the sphere model was conservatively assessed through a simple supported beam analogy to be less than 0.1 mm (ISO Type B). The systematic standard uncertainties from the deflections of the support rods of the reflective markers were estimated from the magnitude of the change in acceleration of the decaying sphere from time zero to the first trough in the heave time series (~16.5 m/s 2 for = 0.5 ), which is in the same order of magnitude as , allowing the deflection to be assessed by including the weight of an additional reflective marker. Conservatively, the systematic standard uncertainties introduced from deflections in the global -direction of the support rods of the reflective markers were included as 0.1 mm (ISO Type B) for = 0.5 . The systematic standard uncertainties for the lower drop heights were linearly scaled down.
Rotations in roll and pitch resulted in small deviations between the measured heave of the sphere model (global coordinate system) and the actual heave, as the reflective markers were placed at a certain distance from the center of rotation (305 mm on average). The motions in heave resulting from the time-dependent roll and pitch were calculated, and the systematic standard uncertainty on the measured heave was found by the rootsum-square (ISO Type A). The maximum measured rotation in pitch or roll is 0.5°-see Figure 16-corresponding to an approximately 0.01 mm decrease of the global -coordinate of the reflective markers. The precision of the motion capture system (incl. calibration) was assessed from displacements of the sphere model in heave with high-precision blocks 50.0 mm in height. By comparing position time series, the systematic standard uncertainty of the motion capture system setup was found to be 0.01 mm (ISO Type A); refer to [10] for further information. The systematic standard uncertainty introduced by vibrations of the bridge (reference frame for the motion capture system) after release of the sphere model was conservatively assessed through a simple supported beam analogy to be less than 0.1 mm (ISO Type B). The systematic standard uncertainties from the deflections of the support rods of the reflective markers were estimated from the magnitude of the change in acceleration of the decaying sphere from time zero to the first trough in the heave time series (~16.5 m/s 2 for H 0 = 0.5D), which is in the same order of magnitude as g, allowing the deflection to be assessed by including the weight of an additional reflective marker. Conservatively, the systematic standard uncertainties introduced from deflections in the global z-direction of the support rods of the reflective markers were included as 0.1 mm (ISO Type B) for H 0 = 0.5D. The systematic standard uncertainties for the lower drop heights were linearly scaled down.
Rotations in roll and pitch resulted in small deviations between the measured heave of the sphere model (global coordinate system) and the actual heave, as the reflective markers were placed at a certain distance from the center of rotation (305 mm on average). The motions in heave resulting from the time-dependent roll and pitch were calculated, and the systematic standard uncertainty on the measured heave was found by the root-sum-square (ISO Type A). The maximum measured rotation in pitch or roll is 0.5 • -see Figure 16 -corresponding to an approximately 0.01 mm decrease of the global z-coordinate of the reflective markers. The mean values of the expanded uncertainty time series for 0 < t/T e0 < 8 multiplied with H 0,m are 0.44, 0.24, and 0.09 mm for the target drop heights of 0.5D, 0.3D, and 0.1D, respectively, which correspond to about 0.3% of the drop height for all cases.

Six DoF Motions
In Figure 16, time series of the six-DoF rigid body motions of the sphere model measured from the optical motion capture system are presented for H 0 = 0.5D. The measured six-DoF motions for H 0 = {0.1D, 0.3D} are presented in Appendix D. The influences on the heave measurements from roll and pitch of the sphere model were included in the uncertainty analysis; see Table 3. The mean values of the expanded uncertainty time series for 0 < ⁄ < 8 multiplied with , are 0.44, 0.24, and 0.09 mm for the target drop heights of 0.5 , 0.3 , and 0.1 , respectively, which correspond to about 0.3% of the drop height for all cases.

Six DoF Motions
In Figure 16, time series of the six-DoF rigid body motions of the sphere model measured from the optical motion capture system are presented for = 0.5 . The measured six-DoF motions for = 0.1 , 0.3 are presented in Appendix D. The influences on the heave measurements from roll and pitch of the sphere model were included in the uncertainty analysis; see Table 3.

Initial Calmness of the Sphere Model
The test case imposes zero velocity and zero acceleration as initial conditions on the sphere. To investigate the initial calmness of the sphere model, the heave (position) time series and time derivatives preceding the drop (i.e., for −0.3 < / < 0), were assessed; see Figure 17. The position time series were subtracted with the respective measured drop heights to get zero as reference value. A moving average filter with a size of 21 samples was utilized to filter the acceleration time series.

Initial Calmness of the Sphere Model
The test case imposes zero velocity and zero acceleration as initial conditions on the sphere. To investigate the initial calmness of the sphere model, the heave (position) time series and time derivatives preceding the drop (i.e., for −0.3 < t/T e0 < 0), were assessed; see Figure 17. The position time series were subtracted with the respective measured drop heights to get zero as reference value. A moving average filter with a size of 21 samples was utilized to filter the acceleration time series.

Frequency Content
The three normalized heave decay time series (Figure 13) with 0 < / < 8 were converted to a periodic signal by mirroring about / = 0; see Figure 18a. The one-sided spectral densities were calculated through FFT analysis; see Figure 18b.

Frequency Content
The three normalized heave decay time series (Figure 13) with 0 < t/T e0 < 8 were converted to a periodic signal by mirroring about t/T e0 = 0; see Figure 18a. The one-sided spectral densities were calculated through FFT analysis; see Figure 18b.

Reflections and Initial Calmness of the Water Phase
The measured surface elevation time series at the three wave gauges can be seen for the highest drop height (four repetitions) in Figure 19. Reflective walls (wave maker) were at 4.22 m from the sphere model location; see Figure 8. A radiated wave needed to travel to the reflective wall and back (i.e., 2 ⋅ 4.22 = 8.44 m), before reaching the sphere model. The time = 8.44/ , where is the celerity of a linear wave with period , is included in Figure 19. Reflected waves propagated past the locations of wave gauges 1, 2, and 3 for around 2.0, 1.3, and 0.7 periods before , respectively. Decay time series presented up to / = 8 are not under the influence of reflections from waves with the period ; see Figure 19. This can be considered a conservative estimate, as the main wave front of radiated waves would have propagated with the group velocity rather than the phase velocity. The measured surface elevations from the other drop heights are included in Appendix D. The initial calmness of the free water surface was assessed by the surface elevation time series prior to the release of the sphere model; see Figure 20.

Reflections and Initial Calmness of the Water Phase
The measured surface elevation time series at the three wave gauges can be seen for the highest drop height (four repetitions) in Figure 19. Reflective walls (wave maker) were at 4.22 m from the sphere model location; see Figure 8. A radiated wave needed to travel to the reflective wall and back (i.e., 2·4.22 = 8.44 m), before reaching the sphere model. The time t r0 = 8.44/c, where c is the celerity of a linear wave with period T e0 , is included in Figure 19. Reflected waves propagated past the locations of wave gauges 1, 2, and 3 for around 2.0, 1.3, and 0.7 periods before t r0 , respectively. Decay time series presented up to t/T e0 = 8 are not under the influence of reflections from waves with the period T e0 ; see Figure 19. This can be considered a conservative estimate, as the main wave front of radiated waves would have propagated with the group velocity rather than the phase velocity. The measured surface elevations from the other drop heights are included in Appendix D.

Reflections and Initial Calmness of the Water Phase
The measured surface elevation time series at the three wave gauges can be seen for the highest drop height (four repetitions) in Figure 19. Reflective walls (wave maker) were at 4.22 m from the sphere model location; see Figure 8. A radiated wave needed to travel to the reflective wall and back (i.e., 2 ⋅ 4.22 = 8.44 m), before reaching the sphere model. The time = 8.44/ , where is the celerity of a linear wave with period , is included in Figure 19. Reflected waves propagated past the locations of wave gauges 1, 2, and 3 for around 2.0, 1.3, and 0.7 periods before , respectively. Decay time series presented up to / = 8 are not under the influence of reflections from waves with the period ; see Figure 19. This can be considered a conservative estimate, as the main wave front of radiated waves would have propagated with the group velocity rather than the phase velocity. The measured surface elevations from the other drop heights are included in Appendix D.   The mean and standard uncertainty of the surface elevation time series for all repetitions and wave gauges over −1 < / < 0 are both 0.0000 m.

Uncertainties of Physical Parameters
The values and standard uncertainties of the physical parameters from the test case in Section 1.1 are presented for the physical tests in Table 4. Standard uncertainties were calculated from the sample standard deviations; see Equation (3). Physical parameters not included in the test case-the influences of which were found to vary insignificantly between indoor laboratories of about 20 °C-are also included in Table 4 to easily be available to the reader (for inclusion in high-fidelity numerical models).

Comparison of Decay Measurements to Numerical Modelling Blind Tests
In the present section, the numerical heave decay time series are presented that were obtained from the numerical models of the test case by modelling approaches of various fidelity, as introduced in Section 1.2, and with the properties outlined in Appendixes A and B. Comparisons of the full time series for all drop heights are shown in Figure 21. In The mean and standard uncertainty of the surface elevation time series for all repetitions and wave gauges over −1 < t/T e0 < 0 are both 0.0000 m.

Uncertainties of Physical Parameters
The values and standard uncertainties of the physical parameters from the test case in Section 1.1 are presented for the physical tests in Table 4. Standard uncertainties were calculated from the sample standard deviations; see Equation (3). Physical parameters not included in the test case-the influences of which were found to vary insignificantly between indoor laboratories of about 20 • C-are also included in Table 4 to easily be available to the reader (for inclusion in high-fidelity numerical models).

Comparison of Decay Measurements to Numerical Modelling Blind Tests
In the present section, the numerical heave decay time series are presented that were obtained from the numerical models of the test case by modelling approaches of various fidelity, as introduced in Section 1.2, and with the properties outlined in Appendices A and B. Comparisons of the full time series for all drop heights are shown in Figure 21. In Figure 22 the initiation of the decay for H 0 = 0.5D is shown. The first trough and crest of the decay time series are shown in Figures 23 and 24, respectively. In Figure 25, the comparison of decay time series is shown merely for the numerical models of higher fidelity, i.e., FNPF and RANS models.
Energies 2020, 13, x FOR PEER REVIEW 18 of 36 Figure 22 the initiation of the decay for = 0.5 is shown. The first trough and crest of the decay time series are shown in Figures 23 and 24, respectively. In Figure 25, the comparison of decay time series is shown merely for the numerical models of higher fidelity, i.e., FNPF and RANS models.

Discussion
The measured heave decay time series are seen in Figure 14. The repeatability between the test repetitions is very high for each of the three drop heights, as seen both from Figure 14 and from the random standard uncertainties of the heave decay time series. On average, these are around 0.07, 0.03, and 0.01 mm, respectively, for the three drop heights in descending order, corresponding to less than 0.1% of the initial respective drop heights. However, the random standard uncertainty is largely time-dependent, and the maxima are factors of 4-7 times larger than the average. In general, the random uncertainty decreases when the sphere model decreases in speed and vice versa. This is both visible over time and over the three investigated drop heights. Over time, two maxima (in magnitude) are expected in the speed time series per natural period, and these maxima damp out over time (to less than 10% of the first maxima after ~5 ); see Figure 17b. This broadly correlates with the time-variation of the expanded uncertainty in Figure 15, for which the timevariation is governed by the random uncertainty (over the systematic). Over the three drop heights, the random uncertainty decreases with the drop height, where obviously the sphere model will oscillate with lower speeds for lower drop heights; see Figures 15 and 16b. The observations of dependence between the random uncertainties and the speed of the sphere model are ascribed to marker-image-shape-distortions increased by higher relative speeds between the optical motion capture system and the test specimen, as reported in [32].
Apart from the systematic uncertainty modelled from the influences of roll and pitch on the heave measurements, the systematic uncertainty is modelled as a time-invariant. The systematic standard uncertainty stemming from the roll and pitch time series does not exceed 0.02 mm, and as the total systematic standard uncertainty is taken as the rootsum-square of elemental systematic standard uncertainties of significantly higher values, the total systematic uncertainty is practically modelled as a time-invariant. As the random uncertainty largely is dictated by the sphere model speed (equal to zero twice per natural period), the dominating nature of the time-varying uncertainty is alternately systematic and random. As the sphere model damps out, it will eventually be dominated by the systematic uncertainties, seen as the offsets in Figure 15. The reader should note that systematic uncertainties are not directly modelled from the test measurements as with random uncertainties, but rather on estimates and engineering judgment. This is indicated by the ISO Type categorization in Section 3; see [11] for further information.
In Figure 13, the normalized heave decay time series for the three drop heights can be seen relative to one another. Most notably, for increasing drop heights, the initial damped natural period in heave increases. This is in accordance with the spectra shown in Figure 18, where the peak in the spectrum for the highest drop height is shifted to a slightly lower frequency.

Discussion
The measured heave decay time series are seen in Figure 14. The repeatability between the test repetitions is very high for each of the three drop heights, as seen both from Figure 14 and from the random standard uncertainties of the heave decay time series. On average, these are around 0.07, 0.03, and 0.01 mm, respectively, for the three drop heights in descending order, corresponding to less than 0.1% of the initial respective drop heights. However, the random standard uncertainty is largely time-dependent, and the maxima are factors of 4-7 times larger than the average. In general, the random uncertainty decreases when the sphere model decreases in speed and vice versa. This is both visible over time and over the three investigated drop heights. Over time, two maxima (in magnitude) are expected in the speed time series per natural period, and these maxima damp out over time (to less than 10% of the first maxima after~5T e0 ); see Figure 17b. This broadly correlates with the time-variation of the expanded uncertainty in Figure 15, for which the time-variation is governed by the random uncertainty (over the systematic). Over the three drop heights, the random uncertainty decreases with the drop height, where obviously the sphere model will oscillate with lower speeds for lower drop heights; see Figures 15 and 16b. The observations of dependence between the random uncertainties and the speed of the sphere model are ascribed to marker-image-shape-distortions increased by higher relative speeds between the optical motion capture system and the test specimen, as reported in [32].
Apart from the systematic uncertainty modelled from the influences of roll and pitch on the heave measurements, the systematic uncertainty is modelled as a time-invariant. The systematic standard uncertainty stemming from the roll and pitch time series does not exceed 0.02 mm, and as the total systematic standard uncertainty is taken as the root-sumsquare of elemental systematic standard uncertainties of significantly higher values, the total systematic uncertainty is practically modelled as a time-invariant. As the random uncertainty largely is dictated by the sphere model speed (equal to zero twice per natural period), the dominating nature of the time-varying uncertainty is alternately systematic and random. As the sphere model damps out, it will eventually be dominated by the systematic uncertainties, seen as the offsets in Figure 15. The reader should note that systematic uncertainties are not directly modelled from the test measurements as with random uncertainties, but rather on estimates and engineering judgment. This is indicated by the ISO Type categorization in Section 3; see [11] for further information.
In Figure 13, the normalized heave decay time series for the three drop heights can be seen relative to one another. Most notably, for increasing drop heights, the initial damped natural period in heave increases. This is in accordance with the spectra shown in Figure 18, where the peak in the spectrum for the highest drop height is shifted to a slightly lower frequency.
The ideal heave decay tests described as the test case in Section 1.1 only allow oscillations in heave (one-DoF system). Naturally, imperfections will activate additional DoF, which under the assumption of rigid body motions are quantified in Figure 16 for H 0 = 0.5D. As reflective markers are mounted on the upper hemisphere of the sphere model, rotations in pitch and/or roll influence the measurement of the position of the sphere model in the global coordinate system. These influences were accounted for in the uncertainty analysis; see Table 3. Slight drifts occur in surge, sway, and yaw during the decay. The drifts have a negligible influence on the heave decays.
The physical parameters from the test case are listed in Table 4; associated standard uncertainties and values of additional physical parameters are not given for the test case. The values given in Table 4 quantify the certainty with which the governing physical parameters of the test setup are known. All physical parameters from the test case comply very well with the values given in Table 4. The relative deviations between the measured drop heights are the largest, but are basically without influence on the presented results, since normalizing with respect to the measured drop height in each repetition practically eliminates deviations between repetitions. The initial calmness of the sphere model and water phase are analyzed from time series preceding the drop; see Figures 16 and 20. Both the sphere model and the water phase are considered completely calm for practical applications.

Comparison to Numerical Modelling Blind Tests
Numerical models have successfully been formulated to represent the test case presented in Section 1.1. The majority of the numerical models depict the heave decay time series from the physical tests very well; see Figure 21. The largest deviations between physical and numerical tests occur for the LPF models, where the deviations are more pronounced for higher drop heights. This was expected, as nonlinearities increasingly govern the heave decay as the drop height is increased. The LPF0 and LPF1 models, introduced in Appendix C, have a significant negative phase shift within the first natural period relative to the physical tests and the models of higher fidelity; see . As a result of the phase shift, large deviations from the 95% CI from the physical tests of around 50 mm (i.e., 33% of H 0 ), occur for the LPF0 and LPF1 models at H 0 = 0.5D. Not considering the phase shifts, but merely the magnitudes of troughs and crests, the LPF0 and LPF1 models, respectively, deviate with around 12-13 and 1-5 mm (i.e., 9% and 1-3% of H 0 ) at the first trough and crest; see Figures 23 and 24. The LPF0 model oscillates with the damped natural frequency of a one-DoF spring-mass-damper system with constant hydrodynamic coefficients, and thus is not capable of including broader frequency contents, which may explain the larger phase shifts for larger drop heights; see Figure 21. The linearization of the hydrostatic force in the LPF1 model spuriously increases the acceleration, as discussed in Appendix C. As the drop height is decreased, the heave decay will oscillate with T e0 and the assumption of linear hydrostatics will become more accurate. Consequently, the LPF0 and LPF1 models become increasingly accurate in both amplitude and phase for lower drop heights; see Figure 21. The inclusion of nonlinear hydrostatics in the LPF2 and 3 models significantly reduces the phase shifts; see, e.g., Figure 23. The constant a ∞ 33 term in the LPF2 model, however, spuriously delays the decay at initiation-see Figure 22-and in general increases the deviation from the physical tests when the sphere is displaced from its rest condition at which the constant a ∞ 33 term is evaluated; see Figures 21 and 24. Only including the draft-dependency of the a ∞ 33 term in the radiation force as in the LPF3 model (see Appendix C) introduces large deviations at the first trough at H 0 = 0.5D; see Figure 23. The inclusion of draft-dependency of the convolution part of the radiation force, as done in the LPF4 model (refer to Appendix C for further information), does not yield more accurate results. Despite the large deviations at the first trough, the LPF3 model captures all subsequent crests and troughs in the H 0 = 0.5D case with an accuracy close to those of the RANS models, and is thus significantly more accurate than the LPF2 model with constant a ∞ 33 . At H 0 = 0.1D, the LPF2 and 3 models perform with maximum deviations of around 1 mm, which are comparable to the deviations of the models of higher fidelity.
The FNPF and RANS models deviate with less than 1 mm for H 0 = 0.1D, corresponding to 3% of H 0 . At the first trough, the models FNPF1, RANS1 and RANS5 lie within the 95% CI of the physical measurements, while the RANS2 and RANS4 models deviate with less than 0.3 mm (i.e., less than 1% of H 0 ). Deviations at the first trough have the same order of magnitude for H 0 = 0.3D, whereas at H 0 = 0.5D, the deviations increase to around 1-3 mm (i.e., 1-2% of H 0 ), with the exception of the RANS2 and RANS3 models, which are actually within the (narrow) 95% CI. The kinematics, and thus velocity gradients, are largest within the first natural period, leading to high demands on the near-wall meshing and treatment (mesh morphing, wall functions, etc.) in the RANS models. However, from Figure 25, there is a general tendency of the largest deviations to occur at 1 < t/T e0 < 4 (even when taking into account the decrease of the CI width; see Figure 15). Assuming the time-error of the motion capture system to be negligible, the reasoning behind the tendency of largest deviations to not occur during the first natural period is two-fold: (i) in a RANS model, errors from the numerical discretization and iterations accumulate, and (ii) turbulence increases over the first periods and when the sphere changes direction. The former includes numerical errors of turbulence parameters if calculated in a turbulence model, while the latter refers to the increase of the complexity of the water phase over time (emergence of high-frequency perturbations of the free surface and sub-grid vortices) and how model errors of either not including a turbulence model (laminar simulations) or the inaccuracies associated with a given model thus become more pronounced with time. The deviations tend to reduce for 4 < t/T e0 which is ascribed to the low amplitudes themselves rather than an increase in the accuracy, as the continued increase in the phase shifts (up to around 0.04 s, i.e., 0.05T e0 ) also suggests. An increased accuracy from inclusion of a turbulence model (k-omega-SST) can be seen by comparing the RANS2 and RANS3 models in Figure 25.
Troughs and crests for the RANS models are calculated with deviations of maximally 1 mm, 2 mm, and 4 mm, respectively, for the three drop heights in ascending order. This corresponds to deviations up to 3% of H 0 . The FNPF model has similar deviations for the two lowest drop heights, while the deviations at H 0 = 0.5D are up to 8 mm or 5% of H 0 . For H 0 = 0.5D, the maximum of deviations at troughs and crests are an order of magnitude higher for the LPF models than the RANS models, which indicates the potential pitfalls of LPF models for large-amplitude motions.

Conclusions
A sphere model was constructed to accurately represent the formulated test case. Physical parameters of the test setup were quantified, and associated uncertainties were generally found to be low. The precision of the physical test results is very high and was quantified by time-varying systematic and random uncertainties of the heave time series. At a 95% confidence level, the uncertainties were on average 0.09, 0.24, and 0.44 mm for the target drop heights in ascending order, corresponding to about 0.3% of the respective drop heights. The uncertainty of the optical motion capture system increased with larger velocities of the test specimen, and for the largest drop height the uncertainty was less than 1.5 mm, corresponding to less than 1% of the drop height.
Strong correlations were found between the physical test results and the results from independent numerical modelling blind tests for LPF, FNPF, and RANS models, ranged with increasing fidelities. At the lowest drop height, the deviations are less than 1 mm for all models, which corresponds to less than 3% of the drop height (disregarding the regular motion model LPF0). Deviations of the LPF models increase for higher drop heights. The performance of the FNPF model is in general better than the LPF models, but deviations are larger than those of the RANS models for the highest drop height. RANS models produce heave decay time series with deviations of 0-4 mm at troughs and crests for the highest drop height, which correspond to 0-3% of the drop height. Deviations are smaller for the lower drop heights. It should be mentioned that the results from the RANS models have a larger spread than the physical results, and various models are outside of the 95% CIs at various periods during the decay. The comparison of the numerical and physical test results suggests that the LPF and partly the FNPF models should be used with care in applications with motions of very large amplitudes, whereas the RANS models, if proper convergence is reached, are capable of producing accurate results for all drop heights.
The high correlations of multiple independent numerical modelling blind tests with the physical tests demonstrate the use of the test case and the physical test results in validating numerical models. Taking this into account, together with the high repeatability and quantified uncertainties of the physical tests, the measured heave decay time series of the sphere model provide a highly accurate solution to the test case, and are thus highly appropriate for numerical model validation. The heave decay time series are made public as a benchmark dataset in the Supplementary Materials of the present paper.
It is the intention of the authors to perform further tests in the future, including motion of the sphere model in waves with PTO and motions in multiple DoF. If the reader is interested in following the future work, he/she is encouraged to become a member of the international working group by contacting the coordinator of the OES modelling task, Kim Nielsen (please request his contact details from the authors of this paper).

Supplementary Materials:
The benchmark dataset of the physical heave decay tests is publicly available from the supplementary material of the present paper online at https://www.mdpi.com/ 1996-1073/14/2/269/s1 and at the OES webpage [5]. In addition, all numerical modelling blind tests of the test case are available. Refer to Appendix A for detailed information about the contents of the supplementary material. Acknowledgments: The authors would like to thank the laboratory technicians at the Ocean and Coastal Engineering Laboratory at Aalborg University for assistance with the sphere model design and construction. A special thanks to Flemming Buus Bendixen and Sintex for the design, construction, and tests of a magnetic release mechanism, which eventually was not included in the final test design.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The Supplementary Materials to the present paper is structured under the folder Datafile with subfolders Descriptions, Experimental results, and Numerical results; see Figure A1. The folder Descriptions includes technical descriptions of the sphere model and the test setup (referred to in Sections 1 and 2). The folders Experimental results and Numerical results contain the results from the heave decay tests performed physically and numerically, respectively. Eleven numerical modelling approaches were performed on the test case, and thus eleven subdirectories are located under Numerical results; see Figure A1. For further information on the specifications of the numerical models, refer to Appendix B.
The results are given as text-files with columns containing time t test setup (referred to in Sections 1 and 2). The folders Experimental results and Numerical results contain the results from the heave decay tests performed physically and numerically, respectively. Eleven numerical modelling approaches were performed on the test case, and thus eleven subdirectories are located under Numerical results; see Figure A1. For further information on the specifications of the numerical models, refer to Appendix B.
The results are given as text-files with columns containing time [s] and heave [m]; see Figure A2.   cally, respectively. Eleven numerical modelling approaches were performed on the test case, and thus eleven subdirectories are located under Numerical results; see Figure A1. For further information on the specifications of the numerical models, refer to Appendix B.
The results are given as text-files with columns containing time [s] and heave [m]; see Figure A2.

Appendix B
As explained in Section 1.2, three categories of numerical models have been applied to the test case: (i) Reynolds-averaged Navier-Stokes (RANS) models, (ii) fully nonlinear potential flow (FNPF) models based on the boundary element method (BEM), and (iii) linear potential flow (LPF) models based on BEM.  ..
where f r,conv is the convolution part of the radiation force-i.e., WAMIT directly outputs the infinite frequency added mass coefficient a ∞ 33 , and the radiation impulse response functions (IRF) is calculated based on the damping coefficients: For a strictly linear model the coefficients are found for the structure located at rest at its statically neutrally buoyant position in the water. The results of such a model are given in the LPF1 model. However, one may try to extend the linear case by introducing nonlinear coefficients. When doing this the effects of the motion of the structure (i.e., the draft of the sphere) are included, but the water surface is considered calm. The easiest and most common first step is to include nonlinear buoyancy, which is done in LPF2. Further, the draft dependency of a ∞ 33 is included in LPF3, and finally, in addition, the radiation convolution function is included in LPF4. The models are outlined in Table A3. Measurements were performed using a force sensor which was connected to the mooring line. Two tests were performed, one test where the sphere was slowly lifted out of the water and the sensor was mounted at the mooring line going upward, and another test where the sphere was slowly submerged into the water and in this case the sensor was mounted under the water at a mooring line going downward. Simultaneous position and force measurements were recorded; see Figure A3. It is seen that the nonlinear Equation (A1) represents the measurements accurately, whereas the linear Equation (A3) is about 50% off when the sphere is fully submerged (x 3 /D = −0.5) or just lifted out of the water (x 3 /D = 0.5). Equation (A1) is utilized in the models with a nonlinear implementation of the hydrostatic force; i.e., LPF2, 3, and 4. (a total of 300 functions). The radiation impulse function to be used at a particular time step during the simulation was thus pieced together of the radiation impulse functions corresponding to the drafts of previous time history. Linear interpolation in the functions was used to get the values corresponding to the actual drafts.
Energies 2020, 13, x FOR PEER REVIEW 29 of 36 corresponding to the drafts of previous time history. Linear interpolation in the functions was used to get the values corresponding to the actual drafts.

Appendix C.7. Comparison of the LPF1-4 Models
A comparison of the simulation results from the LPF1-4 models with various levels of nonlinearities is particularly interesting for the tests conducted with the highest drop height-i.e., = 0.5 . These are shown in Figure A6 for the first two natural periods in heave. For these tests, the initial buoyancy force on the sphere is zero, as the draft is zero. The LPF1 model, however, under-predicts the initial downward hydrostatic force, see Figure A6, since in the linearized hydrostatics assumption, Equation (A1), the buoyancy of a cylinder with a radius and a height equal to the spherical radius is subtracted from the rest condition at = 0 (zero hydrostatic force). In the LPF2 model, the initial downward acceleration of the sphere is over-predicted due to the inclusion of a constant added mass term (the added mass should ideally be zero at initiation). The LPF1 model weighs out this error by the former mentioned error induced by the subtraction of the buoyancy of the cylinder, where it ideally should be the buoyancy of half a sphere. The volume of a cylinder is 1.5 times the volume of a sphere, causing the under-predicted hydrostatic force to exactly balance out the extra added mass ( , = 0.5 ) at initiation. Hence, the LPF1 model accelerates by at initiation, as is the case with the models LPF3 and 4, where the added mass at infinite frequency is calculated as a function of the draft. Regarding the convolution part of the radiation force, the LPF4 model is predicting a different force time series with higher frequency content. Consequently, the LPF4 model has a different response in the heave decay when compared to the LPF3 model. Not including any nonlinearities as in LPF1 model or only including nonlinear hydrostatics, as in the LPF2 model, produces large deviations from the more accurately formulated models with draft-dependent radiation forces implemented; see Figure A6. It is stressed that the comparison to physical tests or numerical models of higher fidelity is needed to evaluate the accuracy of any of the LPF models; see Sections 3 and 4. Appendix C.7. Comparison of the LPF1-4 Models A comparison of the simulation results from the LPF1-4 models with various levels of nonlinearities is particularly interesting for the tests conducted with the highest drop height-i.e., H 0 = 0.5D. These are shown in Figure A6 for the first two natural periods in heave. For these tests, the initial buoyancy force on the sphere is zero, as the draft is zero. The LPF1 model, however, under-predicts the initial downward hydrostatic force, see Figure A6, since in the linearized hydrostatics assumption, Equation (A1), the buoyancy of a cylinder with a radius and a height equal to the spherical radius is subtracted from the rest condition at x 3 = 0 (zero hydrostatic force). In the LPF2 model, the initial downward acceleration of the sphere is over-predicted due to the inclusion of a constant added mass term (the added mass should ideally be zero at initiation). The LPF1 model weighs out this error by the former mentioned error induced by the subtraction of the buoyancy of the cylinder, where it ideally should be the buoyancy of half a sphere. The volume of a cylinder is 1.5 times the volume of a sphere, causing the under-predicted hydrostatic force to exactly balance out the extra added mass (a ∞ 33,LPF1 = 0.5m) at initiation. Hence, the LPF1 model accelerates by g at initiation, as is the case with the models LPF3 and 4, where the added mass at infinite frequency is calculated as a function of the draft. Regarding the convolution part of the radiation force, the LPF4 model is predicting a different force time series with higher frequency content. Consequently, the LPF4 model has a different response in the heave decay when compared to the LPF3 model.
Not including any nonlinearities as in LPF1 model or only including nonlinear hydrostatics, as in the LPF2 model, produces large deviations from the more accurately formulated models with draft-dependent radiation forces implemented; see Figure A6. It is stressed that the comparison to physical tests or numerical models of higher fidelity is needed to evaluate the accuracy of any of the LPF models; see Sections 3 and 4.

Appendix D
Raw and normalized heave decay time series are presented in Figures A7 and A8, respectively. The measured surface elevation time series for the drop heights = 0.1 and = 0.3 are presented in Figure A9. The locations of wave gauges can be seen in Figure 8. The measured motions in all six DoF for the drop heights = 0.1 and = 0.3 are presented in Figures A10 and A11, respectively.        Figure A11. Time series of the measured motions in six DoF for H 0 = 0.3D.