Modelling a Heaving Point-Absorber with a Closed-Loop Control System Using the DualSPHysics Code

: The present work addresses the need for an efﬁcient, versatile, accurate and open-source numerical tool to be used during the design stage of wave energy converters (WECs). The device considered here is the heaving point-absorber developed and tested by Sandia National Laboratories. The smoothed particle hydrodynamics (SPH) method, as implemented in DualSPHysics, is proposed since its meshless approach presents some important advantages when simulating ﬂoating devices. The dynamics of the power take-off system are also modelled by coupling DualSPHysics with the multi-physics library Project Chrono. A satisfactory matching between experimental and numerical results is obtained for: (i) the heave response of the device when forced via its actuator; (ii) the vertical forces acting on the ﬁxed device under regular waves and; (iii) the heave response of the WEC under the action of both regular waves and the actuator force. This proves the ability of the numerical approach proposed to simulate accurately the ﬂuid–structure interaction along with the WEC’s closed-loop control system. In addition, radiation models built from the experimental and WAMIT results are compared with DualSPHysics by plotting the intrinsic impedance in the frequency domain, showing that the SPH method can be also employed for system identiﬁcation.


Introduction
A wide variety of wave energy converter (WEC) designs are currently being developed, but very few are generating electricity commercially (see, e.g., [1,2]). According to the JRC Ocean Energy Status Report (2016 edition) [3], the highest percentage of research and development efforts have been focused on the technology of the point absorbers. These devices consist of a floating buoy whose oscillating motion, heaving and/or pitching, is converted into electricity by means of a power take-off (PTO) system. Modeling of these WECs should account for device kinematics, hydrodynamics, electromechanics, and closedloop control systems. The incorporation of a control system is of particular importance since, in order to extract power in an efficient manner, the PTO unit must be controlled. Furthermore, a significant body of research has shown the benefits of applying more advanced control techniques [4,5].
A number of approaches are typically used to perform design analyses of point absorbers. Physical experiments are widely used and usually play a key role in concept exploration and model validation/identification. Some works focused on the wave tank testing of point-absorbers are [6][7][8][9]. Nevertheless, physical experiments have high costs associated with the construction, deployment and maintenance of the prototypes and facilities. This, along with the exponential growth of the computational resources in the last decades, has boosted the use of numerical methods during the design stages of WECs. Different reviews highlight the main advances of the numerical modelling of WECs achieved during the last decades [10][11][12][13][14]. Furthermore, ref. [15] reviewed several numerical methods focused on point absorbers and concluded that advanced computational tools allow simulation of the device's hydrodynamic response in the time domain as well as integration of different control strategies.
Potential flow models are most often used to reproduce the response of a pointabsorbers under operational sea states (see, e.g., [16][17][18]). These models apply the boundary element method (BEM), based on a discretization of the device boundaries, to solve the wave-body interaction effects in the frequency domain. This solution can be then integrated using the Cummins approach [19] to obtain the time-domain solution. Potential flow codes, such as WAMIT developed by the Massachusetts Institute of Technology [20] or NEMOH developed by the Ecole Centrale de Nantes [21], are based on linear wave theory and small amplitude motion approximation and assume the fluid to be incompressible, inviscid and its motion to be irrotational. Nevertheless, in order to maximize power absorption, point-absorbers are designed to work near resonance conditions, where large motions are common. Therefore, besides the drawback of neglecting viscosity and turbulence, the linear approaches assumed in the potential flow solvers make them unsuitable to reproduce the usual nonlinear dynamics of WECs.
Higher fidelity modeling of WECs is performed using Computational Fluid Dynamics (CFD) methods that do not require any of the previous simplifications but are more time consuming and complex. They solve the Navier-Stokes equations, which model the behavior of the fluid, by numerically discretizing space and time. CFD methods can be classified as mesh-based or meshless. The mesh-based models most frequently used in fluid dynamics consider the fluid as a continuum discretized in control volumes (finite volume method). They require expensive mesh generation and present severe technical challenges to capturing the free surface as well as the nonlinearities within rapidly changing geometries. Some research works where the finite volume method is applied to the analysis of a single-body point-absorber are [22][23][24][25]. When simulating a multi-body device, an ad-hoc approach to model the kinematics is often employed [26]. Among the different meshless models recently developed, the smoothed particle hydrodynamics (SPH) method has become the most popular one for hydrodynamic problems in the presence of freesurface flows [27,28]. SPH is a meshfree method where the fluid is discretized by using discrete sample points, named particles, that move according to their local velocity and carry the physical properties of the fluid with them. The physical properties values of each particle are updated at each time step by an interpolation of the surrounding particles' quantities. Although the SPH method has not achieved the level of mathematical and computational maturity of some mesh-based models, such as the finite volume method, it presents several advantages when modeling free-surface flows with wave-structure interaction: (i) no additional algorithms are required for free-surface detection; (ii) efficient treatment of large deformations; (iii) rapidly moving complex geometries are handled in a straightforward way; (iv) coefficient discontinuities and singular forces are easily implemented. The pioneering works of [29,30] applied the SPH method to oscillating wave surge devices. In the case of point absorbers [31], compared its hydrodynamic response computed with SPH and with a finite volume method, while [32,33] studied its interaction with extreme waves using SPH methods.
DualSPHysics is a SPH-based open-source software developed by the Universidade de Vigo (Spain), The University of Manchester (UK), University of Parma (Italy), Instituto Tecnico de Lisboa (Portugal), New Jersey Institute of Technology (USA) and Universitat Politecnica de Catalunya (Spain). It is considered to be one of the most efficient SPH solvers available [34]. In addition to its Central Processing Unit (CPU) version, it can also be executed on Graphics Processing Unit (GPU) cards, which allows it to have powerful computing capabilities in a personal computer. The coupling of DualSPHysics with the Energies 2021, 14, 760 3 of 20 multi-physics library Project Chrono, developed by the University of Wisconsin-Madison (USA) and University of Parma (Italy) [35], enables the implementation of the mechanical restrictions needed to reproduce the complex mechanism of the PTO systems. DualSPHysics has already been applied to simulate different WECs; [36] studied the interaction of regular waves with a floating oscillating water column and [37] satisfactorily simulated an oscillating wave surge converter, including validation with experimental data. Recently, ref. [38] proved the capability of DualSPHysics to carry out the efficiency and survivability analysis of a heaving point absorber.
In the present work, DualSPHysics coupled with Project Chrono is used to simulate the main experiments conducted when designing a heaving point-absorber under working conditions. The device under study is the 1/17th-scale model designed and built by SANDIA National Laboratories. Different tests are conducted here to prove the code's ability to simulate: (i) radiation tests, where the WEC is forced by means of its actuator in calm water and, (ii) static and dynamic response tests, where the effect of incoming waves on the WEC is considered. Time-domain results from the SPH simulations are compared with the experimental ones. In the case of the radiation model, which is built in the frequency domain, results obtained with DualSPHysics and with potential flow-based solver WAMIT are compared with the experimental model, showing the strengths and weaknesses of each model. The paper is organized as follows: Section 1 provides an introduction and the state-of-the-art of the topic, Section 2 describes in detail the numerical model, Section 3 presents the experimental campaign, Section 4 shows the results of the radiation tests, Section 5 presents the results of the static and dynamic response tests, and Section 6 synthesizes the main conclusions of the work.

Smoothed Particle Hydrodynamics
The meshless numerical approach defined by the smoothed particle hydrodynamics (SPH) method is described here, in particular the DualSPHysics code will be used to perform the simulations. In the present section, the SPH equations implemented in the DualSPHysis code are described, the coupling with Project Chrono is presented and, finally, special attention is paid to the high-performance computing possibilities of the software.

Governing Equations
In SPH, the fluid is discretized into a set of particles where the physical quantities (position, velocity, density, and pressure) are approximated as an interpolation of the values of the neighboring particles. The contribution of each particle depends on the inter-particle distance and the characteristic smoothing length (h) defined by a kernel function (W). The kernel function must fulfil several properties [39], such as positivity inside its compact support, normalization and monotonically decreasing with distance. In the DualSPHysics code, a quintic Wendland kernel [40] with a compact support of radius 2h is applied.
a being the particle where the physical quantities are being calculated and b the surrounding particles, the system of governing equations can be written in the discrete SPH formalism as where t is the time, r is the position, v is the velocity, p is the pressure, ρ is the density, m is the mass, c is the numerical speed of sound, and g is the gravitational acceleration. The artificial viscosity (Π ab ) proposed in [39] will be used and the density diffusion term proposed by [41] is applied in the present simulations using δ = 0.1. Note that the previous three equations have four unknowns: r, v, p and ρ. The system can be closed by solving a Poisson-like equation, where the fluid is considered fully incompressible (ISPH), or by adding an equation of state, where the fluid is treated as weakly compressible (WSPH). In the case of DualSPHysics, the fluid pressure at a certain particle is obtained as a function of its density, according to Tait's equation of state: γ = 7 being the polytropic constant [42], and ρ 0 = 1000 kg·m −3 the reference density of the fluid. The speed of sound, c, is artificially lowered to keep a reasonable time step value (calculated according to the Courant-Friedrichs-Lewy condition) throughout the simulation, while also guaranteeing that density variations are lower than 1%.

Boundary Conditions and Fluid-Driven Objects
Solid boundary conditions have been implemented in different ways in the SPH models, although there is no unanimity yet about the best approach [43]. In the Dual-SPHysics code, solids (such as tanks, floating objects, wavemakers, etc.) are discretized by a set of boundary particles that differ from the fluid ones. The Dynamic Boundary Conditions (DBCs) [44] are boundary particles that satisfy the same equations as fluid particles; however, they do not move according to the forces exerted on them. When a fluid particle approaches a boundary particle and the distance becomes smaller than the interaction distance (2h), the density of the affected boundary particles increases, resulting in a pressure increase. In turn, this results in a repulsive force being exerted on the fluid particle due to the pressure term in the momentum equation (2). DBCs have been successfully applied in some coastal engineering problems due to their capability of discretizing complex 3-D geometries into a set of boundary particles, as presented in the study of wave-runup of armor block breakwaters [45,46]. Nevertheless, DBCs present some drawbacks such as an over dissipation that leads to unphysical large boundary layers. A modification of DBCs (the so called mDBCs) has been presented in [47]. Within this implementation, the boundary particles are arranged in the same way as in the original DBCs, but with the boundary interface located half the initial interparticle distance (dp) away from the innermost layer of boundary particles. For each boundary particle a ghost node is projected into the fluid across that boundary interface. Fluid properties are then computed at the ghost node through a corrected SPH approximation proposed by [48] and finally mirrored back to the boundary particles. This new approach is applied in the simulations of this work since it provides a more accurate and smoother pressure field, resulting in a reduction in the unphysically large boundary layer, as already shown in [48].
On the other hand, fluid-driven objects can be simulated in DualSPHysics in a straightforward way. First, the net force per unit of mass, f , of each floating particle q is computed as the summation of the contributions of all surrounding fluid particles b: where the interactions between particles q and b are solved following the mDBCs approach.
Since the floating object is considered as rigid, its motion can be obtained by solving the basic equations of the rigid body dynamics: Energies 2021, 14, 760

of 20
M being the total mass of the object, V the linear velocity, I the moment of inertia, Ω the rotational velocity and R the center of mass. By integrating Equations (6) and (7) in time, the values of V and Ω at the beginning of the next time step can be predicted. The velocity of each particle of the floating object is given by: Finally, each floating particle q is moved to the position obtained when integrating Equation (8) in time. The work in [49] proved that this approach conserves linear and angular momentum. The research in [50] proved the accuracy of DualSPHysics in simulating fluid-driven objects by comparing the motions of a freely floating box under nonlinear regular waves with experimental data.

Coupling with Project Chrono
As introduced before, the Project Chrono library is coupled with DualSPHysics to solve the fluid-structure interaction problem with mechanical constraints applied on rigid bodies. In this way, the user can define dynamic properties such as friction and restitution coefficients for collisions, as well as restitution forces from springs and damper systems. The last ones are of interest in this work since the controller of the PTO system will be numerically defined as a spring-damper system installed in the WEC. In addition, this coupling procedure allows implementation of a closed-loop system by defining the PTO force as a function of the WEC's position and velocity. A more complete description and validation of the coupling between DualSPHysics and Chrono can be found in [51].

High-Performance Computing
The need to accelerate CFD codes is more demanded in the case of the SPH models due to the Lagrangian nature of the solver. One of the reasons of the high computational cost of the SPH codes is that during the particle interactions, each particle will interact with hundreds of neighboring particles at each time step. As opposed to mesh-based methods, in SPH the computational nodes (particles) are allowed to move, i.e., the particle of interest and its surrounding particles move. Therefore, the neighbors of the particle where the physical quantities are being computed are different each time step. A neighbor search needs to be implemented, increasing the computational workload [52]. For all these reasons the parallelization of SPH models is a more arduous task than in the case of mesh-based codes where the nodes are fixed in time.
DualSPHysics is developed to run on two different hardware architectures: the Central Processing Unit (CPU) and Graphics Processing Unit (GPU). The GPU execution is performed using the NVIDIA's programming framework Compute Unified Device Architecture (CUDA). Particle interactions can be implemented on the GPU for only one particle using one execution CUDA thread to compute the force resulting from the interaction with all its neighbors [53]. In general, all the sequential tasks and operations that involve a loop over all particles are performed using the parallel architecture of the GPU cores [54]. Therefore, DualSPHysics can perform large simulations by using the computational power of the GPUs, but it can be also executed on workstations without GPUs by using the current multi-core CPU version of the code (with OpenMP).

Experimental Campaign
The WEC device studied here is a roughly 1/17th-scale model designed by Sandia National Laboratories [55]. It is schematically depicted in Figure 1 and its most relevant physical parameters are summarized in Table 1. This device was designed to work in three degrees-of-freedom: heave, surge and pitch. However, for cases considered in this study, the surge and pitch motions are locked and the device is allowed to move in heave only. During the experiments, the effect of the PTO system is reproduced by means of an actuator located on top of the float, as shown in Figure 1. Therefore, the total rigid-body mass is the sum of the masses of the float and the actuator slider, as noted in Table 1. The vertical position of the float is measured by the sinusoidal encoder of the actuator, while the vertical force is measured by load cells. More information about the experimental data acquisition can be found in [55][56][57].
is the sum of the masses of the float and the actuator slider, as noted in Table 1. The vertical position of the float is measured by the sinusoidal encoder of the actuator, while the vertical force is measured by load cells. More information about the experimental data acquisition can be found in [55][56][57]. Physical tests were conducted in the Maneuvering and Sea Keeping (MASK) basin located in the Naval Surface Warfare Center, Carderock Division. As shown in Figure 2, the MASK basin is 110 m long, 73 m wide and 6.1 m deep, except for a 15.2 m wide trench parallel to the long side of the basin (dashed line), which has a depth of 10.7 m. Waves were sent diagonally across the basin. The MASK basin is spanned by a 115 m bridge (not shown in Figure 2 for the sake of simplicity), which can move, thanks to a rail system, to the center of the basin width as well as rotate up to 45° from the centerline. The MASK basin carriage, where the models and instrumentation are mounted, is suspended beneath the bridge, and is permitted to travel along the rails. The basin is equipped with a wavemaker system consisting of 216 paddles: 108 paddles along the long side of the basin, 60  Physical tests were conducted in the Maneuvering and Sea Keeping (MASK) basin located in the Naval Surface Warfare Center, Carderock Division. As shown in Figure 2, the MASK basin is 110 m long, 73 m wide and 6.1 m deep, except for a 15.2 m wide trench parallel to the long side of the basin (dashed line), which has a depth of 10.7 m. Waves were sent diagonally across the basin. The MASK basin is spanned by a 115 m bridge (not shown in Figure 2 for the sake of simplicity), which can move, thanks to a rail system, to the center of the basin width as well as rotate up to 45 • from the centerline. The MASK basin carriage, where the models and instrumentation are mounted, is suspended beneath the bridge, and is permitted to travel along the rails. The basin is equipped with a wavemaker system consisting of 216 paddles: 108 paddles along the long side of the basin, 60 paddles in a ninety-degree arc, and 48 paddles along the short edge of the basin. Two artificial beaches, with a 12 • slope, are used to avoid wave reflections from the two edges opposite to the wavemakers. Note that the numerical tanks that will be used for the simulations in the next sections aim to reproduce the same wave-WEC interactions as in the wave basin ( Figure 2). Several reductions in the tank dimensions were applied to reduce the computational runtime, however, due to the geometry of the float, it is necessary to keep all the simulations in 3D. paddles in a ninety-degree arc, and 48 paddles along the short edge of the basin. Two artificial beaches, with a 12° slope, are used to avoid wave reflections from the two edges opposite to the wavemakers. Note that the numerical tanks that will be used for the simulations in the next sections aim to reproduce the same wave-WEC interactions as in the wave basin ( Figure 2). Several reductions in the tank dimensions were applied to reduce the computational runtime, however, due to the geometry of the float, it is necessary to keep all the simulations in 3D. In this research, data from three different experimental campaigns were used: the MASK1, MASK2B and MASK3 campaigns. All of them were conducted by SANDIA National Laboratories between the years 2016 and 2019 and are publicly available at [55], [56] and [57].

Radiation Test
Radiation tests consist of forcing the WEC in otherwise calm water with a known force in the vertical direction via its actuator. In the same fashion as the physical tests, the WEC's motion is numerically restricted such that the float can only move in heave. The results shown in this section consider first a monochromatic force signal (F) and then a multisine force signal as inputs. The block diagram in Figure 3 shows that the radiation FRF (frequency-response-function) is the inverse of the intrinsic impedance, Zi, while the output obtained from the radiation test is the heaving velocity, v, which can be calculated as A suitable numerical tank was designed to perform the radiation tests (the top view is shown in Figure 4). Since the goal of the numerical tank is to mimic open-sea conditions, it is necessary to avoid reflection from any of the boundaries. For this reason, the tank has In this research, data from three different experimental campaigns were used: the MASK1, MASK2B and MASK3 campaigns. All of them were conducted by SANDIA National Laboratories between the years 2016 and 2019 and are publicly available at [55][56][57].

Radiation Test
Radiation tests consist of forcing the WEC in otherwise calm water with a known force in the vertical direction via its actuator. In the same fashion as the physical tests, the WEC's motion is numerically restricted such that the float can only move in heave. The results shown in this section consider first a monochromatic force signal (F) and then a multisine force signal as inputs. The block diagram in Figure 3 shows that the radiation FRF (frequency-response-function) is the inverse of the intrinsic impedance, Z i , while the output obtained from the radiation test is the heaving velocity, v, which can be calculated as paddles in a ninety-degree arc, and 48 paddles along the short edge of the basin. Two artificial beaches, with a 12° slope, are used to avoid wave reflections from the two edges opposite to the wavemakers. Note that the numerical tanks that will be used for the simulations in the next sections aim to reproduce the same wave-WEC interactions as in the wave basin ( Figure 2). Several reductions in the tank dimensions were applied to reduce the computational runtime, however, due to the geometry of the float, it is necessary to keep all the simulations in 3D. In this research, data from three different experimental campaigns were used: the MASK1, MASK2B and MASK3 campaigns. All of them were conducted by SANDIA National Laboratories between the years 2016 and 2019 and are publicly available at [55], [56] and [57].

Radiation Test
Radiation tests consist of forcing the WEC in otherwise calm water with a known force in the vertical direction via its actuator. In the same fashion as the physical tests, the WEC's motion is numerically restricted such that the float can only move in heave. The results shown in this section consider first a monochromatic force signal (F) and then a multisine force signal as inputs. The block diagram in Figure 3 shows that the radiation FRF (frequency-response-function) is the inverse of the intrinsic impedance, Zi, while the output obtained from the radiation test is the heaving velocity, v, which can be calculated as A suitable numerical tank was designed to perform the radiation tests (the top view is shown in Figure 4). Since the goal of the numerical tank is to mimic open-sea conditions, it is necessary to avoid reflection from any of the boundaries. For this reason, the tank has A suitable numerical tank was designed to perform the radiation tests (the top view is shown in Figure 4). Since the goal of the numerical tank is to mimic open-sea conditions, it is necessary to avoid reflection from any of the boundaries. For this reason, the tank has a depth of 1.50 m and the lateral 7.20 m long boundaries form a square with the axisymmetric WEC at its center. Furthermore, periodic conditions [58] are applied to the lateral boundaries and numerical damping is applied to the lateral and bottom water particles, affecting those particles located less than 0.60 m away from the laterals and/or less than 0.20 m away from the bottom. The numerical damping system [59] employs a quadratic decay function to gradually reduce the velocity of fluid particles within the damping zone at each time step according to their location.
Energies 2021, 14, 760 8 of 20 a depth of 1.50 m and the lateral 7.20 m long boundaries form a square with the axisymmetric WEC at its center. Furthermore, periodic conditions [58] are applied to the lateral boundaries and numerical damping is applied to the lateral and bottom water particles, affecting those particles located less than 0.60 m away from the laterals and/or less than 0.20 m away from the bottom. The numerical damping system [59] employs a quadratic decay function to gradually reduce the velocity of fluid particles within the damping zone at each time step according to their location.

Monochromatic Force Signal
The monochromatic force input considered here is the one from Case130 of the MASK3 experimental campaign [57]: a sinusoidal signal with an amplitude of 800 N and a frequency of 0.60 Hz, which corresponds to a period of 1.67 s.
Using the same numerical tank described in Figure 4, different simulations for several resolutions are executed. In the SPH simulations, the resolution is defined by the initial interparticle distance, dp, which is employed to generate the particles that will initially discretize the numerical domain. In this case, three values of the initial interparticle distance are chosen: dp = 3.0, 2.0, 1.5 cm.
The WEC's heave response, in terms of position (z) and velocity (v), for each resolution is shown in Figure 5 and compared with the corresponding experimental time series. It can be noted that the device reaches a steady state in which it oscillates harmonically with the frequency of the actuator force signal. Figure 5 also shows that the heave position and velocity time series satisfy the π/2 rad theoretical phase lag expected between them. Qualitatively, the matching between the numerical and experimental results observed in Figure 5 is satisfactory, both in amplitude and phase, for all resolutions.

Monochromatic Force Signal
The monochromatic force input considered here is the one from Case130 of the MASK3 experimental campaign [57]: a sinusoidal signal with an amplitude of 800 N and a frequency of 0.60 Hz, which corresponds to a period of 1.67 s.
Using the same numerical tank described in Figure 4, different simulations for several resolutions are executed. In the SPH simulations, the resolution is defined by the initial interparticle distance, dp, which is employed to generate the particles that will initially discretize the numerical domain. In this case, three values of the initial interparticle distance are chosen: dp = 3.0, 2.0, 1.5 cm.
The WEC's heave response, in terms of position (z) and velocity (v), for each resolution is shown in Figure 5 and compared with the corresponding experimental time series. It can be noted that the device reaches a steady state in which it oscillates harmonically with the frequency of the actuator force signal. Figure 5 also shows that the heave position and velocity time series satisfy the π/2 rad theoretical phase lag expected between them. Qualitatively, the matching between the numerical and experimental results observed in Figure 5 is satisfactory, both in amplitude and phase, for all resolutions. To quantify the accuracy of the results, the non-dimensional error estimator named index of agreement, d1, is employed. According to [60], it can be defined in a general way as: To quantify the accuracy of the results, the non-dimensional error estimator named index of agreement, d 1 , is employed. According to [60], it can be defined in a general way as: N being the total number of records of the variable, C and E its numerical and experimental values, respectively, and the overbar indicates average. Equation (10) shows that d 1 is bounded between zero and one, where one represents perfect agreement between the two signals. Table 2 shows the index of agreement, total number of particles and runtime for each simulation. The computational times shown in Table 2 correspond to 26 s of physical time simulated in a GeForce RTX 2080 GPU card. Note that the lower the dp, the more particles involved in the simulation and the longer the execution time. The values of d 1 confirm that results are highly accurate for the three values of dp and only a slight improvement is achieved when the two finest resolutions are used.

Polychromatic Force Signal
In this section, a multisine force signal is applied to the WEC in order to obtain the intrinsic impedance of the system and its resonant frequency. This allows the performance an interesting comparison in the frequency domain between results from DualSPHysics, physical experiments and potential-flow code WAMIT. Figure 6 shows that the input force, which covers a range of frequencies from 0.25 to 1.00 Hz, is a pink multisine signal, since its amplitude (and therefore, its power spectral density) decreases with frequency. To set the resolution and length of the input force signal, a compromise between accuracy and computational cost must be achieved. Table 2 shows that the accuracy is slightly improved at a reasonable computational cost when using dp = 2.0 cm, and Figure 7a shows that the magnitude of the force applied is, especially at certain instants, significantly lower than in the monochromatic case. Therefore, dp = 2.0 cm is used here. To keep the computational times lower than a week (using a GeForce RTX 2080 GPU card), a single period of 100 s of the multisine signal is employed.

Polychromatic Force Signal
In this section, a multisine force signal is applied to the WEC in order to obtain the intrinsic impedance of the system and its resonant frequency. This allows the performance an interesting comparison in the frequency domain between results from DualSPHysics, physical experiments and potential-flow code WAMIT. Figure 6 shows that the input force, which covers a range of frequencies from 0.25 to 1.00 Hz, is a pink multisine signal, since its amplitude (and therefore, its power spectral density) decreases with frequency. To set the resolution and length of the input force signal, a compromise between accuracy and computational cost must be achieved. Table 2 shows that the accuracy is slightly improved at a reasonable computational cost when using dp = 2.0 cm, and Figure 7a shows that the magnitude of the force applied is, especially at certain instants, significantly lower than in the monochromatic case. Therefore, dp = 2.0 cm is used here. To keep the computational times lower than a week (using a GeForce RTX 2080 GPU card), a single period of 100 s of the multisine signal is employed.       Figure 7a, is applied to the device. Since the input signal is periodic (even though only one period is considered), the force and velocity can be described by their complex quantities in the frequency-domain, denoted asF andv, respectively. In this way, the block diagram in Figure 3 and Equation (9) show that the intrinsic impedance can be obtained as: whereF andv are calculated from their time series (Figure 7) using the discrete Fourier transform via the Fast Fourier Transform (FFT). The same procedure is used to obtain the intrinsic impedance from the experimental force input signal and heave velocity output [61]. More specifically, data from Case063 of the MASK2B experimental campaign [56] are considered here, where a single period of a pink multisine signal is also used but, in this case, with a total period of 300 s (instead of the 100 s period used in the SPH simulation). This difference between the numerical and experimental periods is important because the lower the period, the fewer frequency components considered, i.e., lower frequency resolution. This is because the frequency sampling is inversely proportional to the period thus, in this case, the frequency sampling is 1/100 Hz for numerical results and 1/300 Hz for experimental results.
The intrinsic impedance is a frequency dependent complex-valued variable that accounts for the effects of inertia, radiation, viscosity and hydrostatic stiffness. According to [62], the intrinsic impedance, Z i (ω), can be modelled as: where the viscous force has been assumed to be linear with velocity via a friction coefficient, b f , the radiation force is defined by the radiation damping, b r , and the added mass, m a , the hydrostatic force depends on the hydrostatic coefficient, k hs , and the inertia force simply depends on the mass of the device, M. It is important to note that Equations (11) and (12) allow attainment of the resonance frequency of the device, i.e., the frequency at which resonance is achieved. As stated in [62], resonance occurs when the heave velocity is in phase with the excitation force. Thus, the resonance frequency is the one that guarantees that Im(Z i ) = 0 Ns/m (which directly implies ∠Z i = 0 rad). Figure 8 depicts the magnitude (top left), phase (bottom left), real part (top right) and imaginary part (bottom right) of the intrinsic impedance in the frequency domain. The results of the DualSPHysics simulation (in blue) are compared with the experimental ones (in orange) [61] and with the ones calculated using the potential-flow method WAMIT (in black). Both signals, the experimental and SPH one, have some noise; this noise being more relevant in the latter due to the lower frequency resolution employed, as previously explained. To avoid this undesirable but inevitable effect, a second order Savitzky-Golay filter is used to smooth the SPH and experimental results (dashed lines in Figure 8). Leaving aside the plot of Re(Z i ), Figure 8 demonstrates a high level of agreement between the SPH, experimental and WAMIT results. The resonance frequency, obtained as the frequency at which the imaginary part (or, equivalently, the phase) of the intrinsic impedance crosses the x-axis, is around 0.65 Hz for all cases. Regarding the plot of the real part of Z i shown in Figure 8, the difference between the experimental and WAMIT curves is roughly constant with frequency, revealing an offset caused by the fact that WAMIT is neglecting the viscous effects, which have a significant impact on Re(Z i ), as shown in Equation (12). This mismatch was investigated in depth in [61]. The real part of the intrinsic impedance obtained with SPH nicely matches the experimental data at high frequencies; however, when low frequencies are considered, Re(Z i ) is highly underestimated. This means that when the frequency is high enough, the SPH simulation is able to reproduce the physics of the problem correctly, but when the frequency is low, there is a phenomenon, most likely related with viscosity, that is not being adequately captured. the physics of the problem correctly, but when the frequency is low, there is a phenomenon, most likely related with viscosity, that is not being adequately captured.

Static and Dynamic Response Test
In this section, the effect of incoming waves on the WEC is considered. First, a static excitation test is conducted by keeping the device fixed and calculating the forces exerted by regular waves of height H = 0.09 m and period T = 1.54 s (which corresponds to a wavelength L = 3.70 m in deep-water conditions). Then, during the dynamic response test the device is allowed to oscillate, and its heave response is calculated as a result of the incoming waves and the actuator force. The regular waves have, in this case, a wave height of H = 0.11 m and period of T = 1.58 s (corresponding, in deep water, to a wavelength of L = 3.90 m).
The numerical domain employed in the simulations is common for the two aforementioned cases since the wave conditions are very similar. A detailed diagram of the top

Static and Dynamic Response Test
In this section, the effect of incoming waves on the WEC is considered. First, a static excitation test is conducted by keeping the device fixed and calculating the forces exerted by regular waves of height H = 0.09 m and period T = 1.54 s (which corresponds to a wavelength L = 3.70 m in deep-water conditions). Then, during the dynamic response test the device is allowed to oscillate, and its heave response is calculated as a result of the incoming waves and the actuator force. The regular waves have, in this case, a wave height of H = 0.11 m and period of T = 1.58 s (corresponding, in deep water, to a wavelength of L = 3.90 m).
The numerical domain employed in the simulations is common for the two aforementioned cases since the wave conditions are very similar. A detailed diagram of the top and lateral views of the numerical domain is shown in Figure 9. Although the numerical tank is very different from the experimental basin to reduce the computational runtime, it has been carefully designed such that the wave conditions acting on the device are the same. The reduction in the water depth from 6.10 m to 2.00 m is such that the deep-water condition present in the experimental campaign is conserved and, therefore, the wavelength remains constant. The correct wave propagation is guaranteed by placing the device one wavelength away from the wavemaker (considering the case with longer waves). The wavemaker consists of 17 independent pistons, each one equipped with an active wave absorption system [59] that absorbs the waves radiated by the device, guaranteeing an adequate wave generation throughout the entire simulation. To avoid lateral reflection, the width of the numerical domain is set to 5.23 m (roughly three times the buoy diameter). In addition, periodic boundary conditions are also applied in the lateral limits, and the numerical damping mechanism described in Section 4 is used in the 0.40-m wide zones closer to the lateral boundaries (to damp the transversal velocities). Finally, a 1:2 steep beach together with numerical damping (to damp the longitudinal velocities) prevents any possible reflection downstream of the device. The overall absorption capabilities of the beach and numerical damping used together here are above 96% for both wave conditions, which means that over 96% of the incident wave energy is being dissipated. the beach and numerical damping used together here are above 96% for both wave conditions, which means that over 96% of the incident wave energy is being dissipated.

Static Response Test
The physical test Case004 from the MASK1 campaign [55] is replicated here numerically using the DualSPHysics code. It consists of obtaining the forces exerted by regular waves of height H = 0.09 m and period T = 1.54 s (which corresponds to L = 3.70 m) on a fixed WEC. Denoting the free-surface elevation with η, its block diagram is shown in Figure 10.

Static Response Test
The physical test Case004 from the MASK1 campaign [55] is replicated here numerically using the DualSPHysics code. It consists of obtaining the forces exerted by regular waves of height H = 0.09 m and period T = 1.54 s (which corresponds to L = 3.70 m) on a fixed WEC. Denoting the free-surface elevation with η, its block diagram is shown in Figure 10.

Static Response Test
The physical test Case004 from the MASK1 campaign [55] is replicated here numerically using the DualSPHysics code. It consists of obtaining the forces exerted by regular waves of height H = 0.09 m and period T = 1.54 s (which corresponds to L = 3.70 m) on a fixed WEC. Denoting the free-surface elevation with η, its block diagram is shown in Figure 10.  Figure 11 shows the experimental and numerical time series of the vertical force acting on the buoy for three different resolutions, namely dp = 3.0, 2.0 and 1.5 cm. In this figure, the exact instantaneous values of the experimental force are depicted with black crosses while the corresponding fitted curve, generated using a fourth order Savitzky-Golay function, is represented by a black dashed line. As expected, the obtained responses, both numerical and experimentally, are sinusoidal and their periods are equal to the wave period. The figure shows, qualitatively, that the results obtained with any of the resolutions employed provide a good matching, in terms of amplitude and phase, with the experimental curve. Note, nevertheless, that the force output obtained with dp = 3.0 cm presents an important amount of noise during its peaks, while finer resolutions do not. The index of agreement, d1, defined in Equation (10) is used here as a quantitive measure  Figure 11 shows the experimental and numerical time series of the vertical force acting on the buoy for three different resolutions, namely dp = 3.0, 2.0 and 1.5 cm. In this figure, the exact instantaneous values of the experimental force are depicted with black crosses while the corresponding fitted curve, generated using a fourth order Savitzky-Golay function, is represented by a black dashed line. As expected, the obtained responses, both numerical and experimentally, are sinusoidal and their periods are equal to the wave period. The figure shows, qualitatively, that the results obtained with any of the resolutions employed provide a good matching, in terms of amplitude and phase, with the experimental curve. Note, nevertheless, that the force output obtained with dp = 3.0 cm presents an important amount of noise during its peaks, while finer resolutions do not. The index of agreement, d 1 , defined in Equation (10) is used here as a quantitive measure of the total error in each case. Its values, shown in Table 3, certify the good matching observed in Figure 11, as well as the fact that the accuracy slightly improves when decreasing dp. Table 3 also includes the runtimes of the execution of 26 s of physical time using the same GPU card indicated before.  Table 3, certify the good matching observed in Figure 11, as well as the fact that the accuracy slightly improves when decreasing dp. Table 3 also includes the runtimes of the execution of 26 s of physical time using the same GPU card indicated before. Figure 11. Numerical and experimental time series of the force for three values of dp.

Dynamic Response Test
During the dynamic response test, the WEC moves in heave under the action of waves and the actuator force, Fa. This actuator force may be either predefined externally (open-loop system) or defined during the test by means of the feedback provided by the PTO (closed-loop system). Both possibilities, depicted in Figure 12 via block diagrams, are considered in this section based on Case149 of the MASK3 experimental campaign [57]. During the experiment, the actuator force applied to the WEC is the reaction of the PTO force, i.e., Fa = − FPTO; FPTO being modelled as a proportional-integral (PI) controller that depends on the device's heave position, z, and velocity, v [63]: where kPTO and bPTO are, respectively, the stiffness and damping coefficient of the PTO system or, equivalently, the integral and proportional gains of the PI controller.

Dynamic Response Test
During the dynamic response test, the WEC moves in heave under the action of waves and the actuator force, F a . This actuator force may be either predefined externally (open-loop system) or defined during the test by means of the feedback provided by the PTO (closed-loop system). Both possibilities, depicted in Figure 12 via block diagrams, are considered in this section based on Case149 of the MASK3 experimental campaign [57]. During the experiment, the actuator force applied to the WEC is the reaction of the PTO force, i.e., F a = − F PTO ; F PTO being modelled as a proportional-integral (PI) controller that depends on the device's heave position, z, and velocity, v [63]: where k PTO and b PTO are, respectively, the stiffness and damping coefficient of the PTO system or, equivalently, the integral and proportional gains of the PI controller.
waves and the actuator force, Fa. This actuator force may be either predefined externally (open-loop system) or defined during the test by means of the feedback provided by the PTO (closed-loop system). Both possibilities, depicted in Figure 12 via block diagrams, are considered in this section based on Case149 of the MASK3 experimental campaign [57]. During the experiment, the actuator force applied to the WEC is the reaction of the PTO force, i.e., Fa = − FPTO; FPTO being modelled as a proportional-integral (PI) controller that depends on the device's heave position, z, and velocity, v [63]: where kPTO and bPTO are, respectively, the stiffness and damping coefficient of the PTO system or, equivalently, the integral and proportional gains of the PI controller.
(a) (b) To search the space of controller gains for the optimal tuning around a predicted value, the proportional and integral gains vary every 45 s during the physical experiment. In this section, however, a time window where the gains remain constant is selected from the physical experiment [57]. In particular, t ∈ (271,316) is considered, where the actuator To search the space of controller gains for the optimal tuning around a predicted value, the proportional and integral gains vary every 45 s during the physical experiment. In this section, however, a time window where the gains remain constant is selected from the physical experiment [57]. In particular, t ∈ (271,316) is considered, where the actuator force (shown in Figure 13) is defined by k PTO = 3527 N/m and b PTO = −1754 Ns/m. The wave conditions remain constant throughout the entire experiment: regular waves of height H = 0.11 m and period T = 1.58 s, which qualifies as a deep-water wave with a wavelength of L = 3.90 m. Figure 13 shows that, as expected, the frequency of the actuator force is coincident with the wave frequency and that, after the transient observed during the first period, its amplitude is approximately constant around 160 N. This initial transient is avoided in the simulations by considering the time series from t = 273 s.  Figure 13 shows that, as expected, the frequency of the actuator force is coincident with the wave frequency and that, after the transient observed during the first period, its amplitude is approximately constant around 160 N. This initial transient is avoided in the simulations by considering the time series from t = 273 s. The numerical tank shown in Figure 9 is employed in all the simulations in this section. From the information available in Tables 2 and 3, it is concluded that dp = 2 cm provides a good compromise between accuracy and computational cost, while avoiding any possible noise in the results.

Open-loop System
The time series of the experimental actuator force ( Figure 13) is applied here to the numerical device as a predefined input signal in DualSPHysics. The WEC's heave response is the result of the interaction between the externally-imposed actuator force and the action of regular waves. Therefore, the phase between the wave elevation and the actuator force has a very important effect on the results. Even though it is not possible to know with certainty the wave elevation at the WEC location, the heave position z is available and its phase with respect to Fa can be effectively used instead. Three different phases are simulated simply by delaying the time at which the external actuator force is applied in each simulation (cases A, B and C). The actuator force and heave position time series for the three numerical cases are shown, along with the experimental ones, in Figure 14.
Regarding the numerical results, the actuator force signal is shifting in phase while the heave position is not, causing a phase difference between them of about π rad in case A, 3π/4 rad in case B, and π/2 rad in case C. The experimental phase between Fa and z, observed in the bottom plot of Figure 14, is very similar to the one obtained with case C, i.e., approximately π/2 rad. The numerical tank shown in Figure 9 is employed in all the simulations in this section. From the information available in Tables 2 and 3, it is concluded that dp = 2 cm provides a good compromise between accuracy and computational cost, while avoiding any possible noise in the results.

Open-Loop System
The time series of the experimental actuator force ( Figure 13) is applied here to the numerical device as a predefined input signal in DualSPHysics. The WEC's heave response is the result of the interaction between the externally-imposed actuator force and the action of regular waves. Therefore, the phase between the wave elevation and the actuator force has a very important effect on the results. Even though it is not possible to know with certainty the wave elevation at the WEC location, the heave position z is available and its phase with respect to F a can be effectively used instead. Three different phases are simulated simply by delaying the time at which the external actuator force is applied in each simulation (cases A, B and C). The actuator force and heave position time series for the three numerical cases are shown, along with the experimental ones, in Figure 14. Regarding the numerical results, the actuator force signal is shifting in phase while the heave position is not, causing a phase difference between them of about π rad in case A, 3π/4 rad in case B, and π/2 rad in case C. The experimental phase between F a and z, observed in the bottom plot of Figure 14, is very similar to the one obtained with case C, i.e., approximately π/2 rad.   Figure 15 compares the heave position and velocity time series of each numerical case with the experimental results. The three cases overestimate the amplitudes of z and v but, as the phase lag between F a and z approaches the experimental one ( Figure 14), so do the amplitudes of the heave position and velocity. However, the phase difference is not known before running the simulation, thus it becomes a trial-and-error process. This approach, besides being inefficient, is also unrealistic because the actuator force is very rarely given by a predefined signal, but rather by its PTO force.

Closed-Loop System
The actuator force is applied to the WEC in the physical experiment [57] by means of a closed-loop system as the reaction of the PTO force. The PTO force depends on the device's heave response according to the coefficients k PTO and b PTO , as described in Equation (13). Project Chrono allows numerical implementation of this feedback loop simply by setting k PTO = 3527 N/m and b PTO = −1754 Ns/m. Therefore, thanks to the coupling of DualSPHysics with Project Chrono it is possible to simulate the dynamic response test with its corresponding closed-loop system using the numerical tank defined in Figure 9 and the wave conditions previously described (H = 0.11 m and T = 1.58 s). Figure 16 shows the numerical (solid lines) and experimental (dashed lines) time series of the actuator force and heave position. The numerical phase difference between F a and z now matches the experimental one. This is automatically accomplished because the actuator force is being calculated at each time step as a function of the instantaneous heave position and velocity, instead of being pre-imposed externally.
Energies 2021, 14, x FOR PEER REVIEW 17 of 22 Figure 15 compares the heave position and velocity time series of each numerical case with the experimental results. The three cases overestimate the amplitudes of z and v but, as the phase lag between Fa and z approaches the experimental one (Figure 14), so do the amplitudes of the heave position and velocity. However, the phase difference is not known before running the simulation, thus it becomes a trial-and-error process. This approach, besides being inefficient, is also unrealistic because the actuator force is very rarely given by a predefined signal, but rather by its PTO force.

Closed-loop System
The actuator force is applied to the WEC in the physical experiment [57] by means of a closed-loop system as the reaction of the PTO force. The PTO force depends on the device's heave response according to the coefficients kPTO and bPTO, as described in Equation (13). Project Chrono allows numerical implementation of this feedback loop simply by setting kPTO = 3527 N/m and bPTO = −1754 Ns/m. Therefore, thanks to the coupling of Du-alSPHysics with Project Chrono it is possible to simulate the dynamic response test with its corresponding closed-loop system using the numerical tank defined in Figure 9 and the wave conditions previously described (H = 0.11 m and T = 1.58 s). Figure 16 shows the numerical (solid lines) and experimental (dashed lines) time series of the actuator force and heave position. The numerical phase difference between Fa and z now matches the experimental one. This is automatically accomplished because the actuator force is being calculated at each time step as a function of the instantaneous heave position and velocity, instead of being pre-imposed externally.  Figure 17 shows the heave position and velocity obtained numerically and experimentally. The matching in terms of phase is nearly perfect, while in terms of amplitude it is slightly overestimated. If the index of agreement (Equation (10)) is used to quantify the error, a value of d1 = 0.92 is obtained for both z and v when using a closed-loop system, which is significantly higher than the one obtained with the most accurate case when using an open-loop system (d1 = 0.74). By comparing results in Figure 17 with the ones in Figure 15, it can be also seen that the amplitude overestimation obtained now with a closed-loop system is lower than in any of the three cases with an open-loop system. Note as well that the amplitude of the heave movement is of approximately 0.02 m, which is the value of the dp used, meaning that small motions can be correctly simulated even with relatively coarse resolutions. Therefore, this numerical simulation satisfactorily reproduces the response of the WEC equipped with a proportional-integral controller.  Figure 17 shows the heave position and velocity obtained numerically and experimentally. The matching in terms of phase is nearly perfect, while in terms of amplitude it is slightly overestimated. If the index of agreement (Equation (10)) is used to quantify the error, a value of d 1 = 0.92 is obtained for both z and v when using a closed-loop system, which is significantly higher than the one obtained with the most accurate case when using an open-loop system (d 1 = 0.74). By comparing results in Figure 17 with the ones in Figure 15, it can be also seen that the amplitude overestimation obtained now with a closed-loop system is lower than in any of the three cases with an open-loop system. Note as well that the amplitude of the heave movement is of approximately 0.02 m, which is the value of the dp used, meaning that small motions can be correctly simulated even with relatively coarse resolutions. Therefore, this numerical simulation satisfactorily reproduces the response of the WEC equipped with a proportional-integral controller.

Conclusions
This work aims to prove the capability of the DualSPHysics code coupled with the multiphysics Project Chrono library to reproduce the behavior of a point-absorber under different conditions, including the closed-loop controller defined by the PTO system. By comparing numerical results with physical tests conducted by Sandia National Laboratories, DualSPHysics has proved to be a valuable numerical tool to conduct experiments that are fundamental during the design stage of point absorbers. In this work, results were only compared in terms of the forces acting on the WEC and its motions; however, Du-alSPHysics also has the ability to easily calculate other parameters (such as pressure, particles trajectories or the vorticity field) that can be useful when modelling and optimizing WECs in future works.
The radiation models built from the WAMIT and DualSPHysics simulations were compared with the experimental model by plotting the intrinsic impedance calculated in each case in the frequency domain. It was found that the SPH results match the experimental ones except for the real part of Zi at low frequencies, that is significantly underestimated. On the other hand, WAMIT also provides good results except for Re(Zi), which is underestimated at all frequencies. This offset, which is roughly constant with frequency, is due to the fact that the potential-flow solver neglects viscosity. Thus, the most efficient way to numerically obtain the radiation model is to use WAMIT (or another potentialflow solver) but adding a correction that takes into account the viscous effects. The proven ability of the SPH method to reproduce the viscous effects at high frequencies allows for an estimation of the viscous correction needed for the potential-flow solver. Therefore, a fast and reliable procedure for system identification can be set up by combining the results obtained with a CFD simulation and a potential solver.
It has been shown that the heave response of the device forced via its actuator in calmed water can be accurately obtained even using coarse resolutions with DualSPHys-

Conclusions
This work aims to prove the capability of the DualSPHysics code coupled with the multiphysics Project Chrono library to reproduce the behavior of a point-absorber under different conditions, including the closed-loop controller defined by the PTO system. By comparing numerical results with physical tests conducted by Sandia National Laboratories, DualSPHysics has proved to be a valuable numerical tool to conduct experiments that are fundamental during the design stage of point absorbers. In this work, results were only compared in terms of the forces acting on the WEC and its motions; however, Dual-SPHysics also has the ability to easily calculate other parameters (such as pressure, particles trajectories or the vorticity field) that can be useful when modelling and optimizing WECs in future works.
The radiation models built from the WAMIT and DualSPHysics simulations were compared with the experimental model by plotting the intrinsic impedance calculated in each case in the frequency domain. It was found that the SPH results match the experimental ones except for the real part of Z i at low frequencies, that is significantly underestimated. On the other hand, WAMIT also provides good results except for Re(Z i ), which is underestimated at all frequencies. This offset, which is roughly constant with frequency, is due to the fact that the potential-flow solver neglects viscosity. Thus, the most efficient way to numerically obtain the radiation model is to use WAMIT (or another potential-flow solver) but adding a correction that takes into account the viscous effects. The proven ability of the SPH method to reproduce the viscous effects at high frequencies allows for an estimation of the viscous correction needed for the potential-flow solver. Therefore, a fast and reliable procedure for system identification can be set up by combining the results obtained with a CFD simulation and a potential solver.
It has been shown that the heave response of the device forced via its actuator in calmed water can be accurately obtained even using coarse resolutions with DualSPHysics. The forces exerted by regular waves on the fixed WEC have also been correctly calculated. Finally, during the dynamic response test (when both regular waves and the actuator force are applied simultaneously), it has been proven that the closed-loop system needs to be considered in order to reproduce the nature of the physical experiments. Thanks to the coupling of DualSPHysics with Project Chrono, it is possible to model the PTO force that defines the feedback loop. In this way, the experimental phase between the two inputs is automatically satisfied and the numerical heave response matches satisfactorily the experimental one. Funding: Funded by Xunta de Galicia (Spain) under project ED431C 2017/64 "Programa de Consolidación e Estructuración de Unidades de Investigación Competitivas (Grupos de Referencia Competitiva)" cofunded by European Regional Development Fund (ERDF). Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.