Influence of the Drag Force on the Average Absorbed Power of Heaving Wave Energy Converters Using Smoothed Particle Hydrodynamics

Abstract

In this paper, we investigated how the added mass, the hydrodynamic damping and the drag coefficient of a Wave Energy Converter (WEC) can be calculated using DualSPHysics. DualSPHysics is a software application that applies the Smoothed Particle Hydrodynamics (SPH) method, a Lagrangian meshless method used in a growing range of applications within the field of Computational Fluid Dynamics (CFD). Furthermore, the effect of the drag force on the WEC’s motion and average absorbed power is analyzed. Particularly under controlled conditions and in the resonance region, the drag force becomes significant and can greatly reduce the average absorbed power of a heaving point absorber. Once the drag coefficient has been determined, it is used in a modified equation of motion in the frequency domain, taking into account the effect of the drag force. Three different methods were compared for the calculation of the average absorbed power: linear potential flow theory, linear potential flow theory modified to take the drag force into account and DualSPHysics. This comparison showed the considerable effect of the drag force in the resonance region. Calculations of the drag coefficient were carried out for three point absorber WECs: one spherical WEC and two cylindrical WECs. Simulations in regular waves were performed for one cylindrical WEC with two different power take-off (PTO) systems: a linear damping and a Coulomb damping PTO system. The Coulomb damping PTO system was added in the numerical coupling between DualSPHysics and Project Chrono. Furthermore, we considered the optimal PTO system damping coefficient taking the effect of the drag force into account.

1. Introduction

Wave energy is a potential source of clean electricity that can make a significant contribution to the de-carbonization of the world’s electricity supply. A wide variety of devices to harvest wave energy, known as wave energy converters (WECs), have been designed over recent decades [1]. However, none of them have yet achieved the level of technological development needed to be economically viable [2]. The type of WEC studied in the current research is the heaving point absorber, one of the most investigated types of WECs [3]. These devices typically consist of a floating buoy moved by the waves and connected to a PTO system, which converts the float’s movement into electricity. In order to assess the performance and survivability of wave energy converters, which are necessary for exploiting wave energy, the related wave-induced hydrodynamic forces and WEC motions have to be investigated. Physical experiments are widely used; nevertheless, they have very high costs and require a high level of expertise, and scaling may become an important issue in some cases. On the other hand, numerical methods have become very popular in recent years [4], mainly due to the unprecedented growth of the computational resources available. A complete review on the numerical methods used to simulate the hydrodynamic response of point absorbers can be found in [5].

Most studies concerning WECs employ potential flow theory based on the linearized form of the Navier–Stokes equations. Applying linear potential flow theory allows the numerical modeling of WECs in the time or frequency domain [6], enabling fast calculations of the WEC’s motion. In this research, the open-source linear potential flow software NEMOH [7] is used for the calculation of the hydrodynamic coefficients and the motion of the WEC. Furthermore, this method needs to assume small amplitude oscillations of the WEC and the fluid to be incompressible and inviscid and have an irrotational motion. This results in an underestimation of the wave-induced forces on WECs under highly-nonlinear sea states [8]. Previous simplifications are too restrictive when modeling WECs; therefore, it may be more appropriate to employ higher fidelity models, generally known as CFD (Computational Fluid Dynamics) methods. The most commonly used CFD methods in hydrodynamics are mesh-based. Different point absorbers have been studied numerically using these methods in [9,10]. Despite being very robust mathematically and computationally, mesh-based methods face important challenges when capturing the free surface and the rapidly-evolving nonlinearities. Due to their ability to overcome these drawbacks, meshless CFD methods have gained attention in recent years, with the Smoothed Particle Hydrodynamics (SPH) method being the most widely used [11,12,13]. In contrast with the mesh-based methods, in SPH, the fluid is discretized in a series of points, named particles, that move with the velocity calculated from the Navier–Stokes equations carrying all physical properties with them [14]. The meshless Lagrangian formulation makes the SPH method a very interesting alternative when simulating free-surface flows with a wave–structure interaction, such as the case investigated in this work. In [15], the hydrodynamic response of a point absorber was computed using SPH and a finite volume solver. In [16,17], the authors exploited the capabilities of the SPH to model nonlinearities to study the interaction between a point absorber and extreme waves. In this research, the open-source software DualSPHysics [18] (available at www.dual.sphysics.org) is employed to obtain the hydrodynamic and drag coefficients of different point absorbers. Furthermore, DualSPHysics has been recently coupled [19] with Project Chrono [20], enabling the effect of the PTO system to be included in the simulations. This software has proven to be a valuable tool in the modeling of the wave–structure interaction in general and of floating WECs in particular (see [17,21,22,23,24]).

Estimating the hydrodynamic coefficients of floating bodies has been done before with CFD [25,26,27] and with SPH [28]. The drag coefficient has also been estimated using SPH in [29,30,31], but these studies were limited only to a fixed body in a current. In this research, three point absorbers with different float geometries are considered: one spherical and two cylindrical. The drag coefficient for each is calculated with DualSPHysics and validated with results from previous CFD simulations in the case of the sphere [32] or with experimental data in the cylindrical cases [33,34]. The linear model (in the frequency domain) has also been extended in order to take into account the effect of the drag force by writing the drag force term from the Morison equation in its Fourier series. Finally, the motion and absorbed power of a cylindrical heaving WEC are calculated using the original linear model, the extended linear model and the SPH method. These results are obtained considering two different kinds of PTO systems: a linear damper and a Coulomb damper. In order to model the Coulomb damper PTO system with DualSPHysics, the DualSPHysics–Project Chrono coupling was extended.

The paper is organized as follows: Section 2 describes the basic theoretical principles of SPH and its implementation in DualSPHysics; Section 3 provides a description of the methodology followed to obtain the hydrodynamic and drag coefficients and the modified equation of motion that includes the effect of the drag force; Section 4 outlines the numerical setup employed and the different test cases; Section 5 discusses the main results obtained using the original linear theory, the extended linear theory and DualSPHysics; finally, Section 6 presents an overview of the conclusions.

2. Smoothed Particle Hydrodynamics—DualSPHysics

DualSPHysics applies the SPH method, which is a meshless Lagrangian method used in a growing range of applications within the field of CFD. In SPH, the fluid is discretized in a set of particles for which the position, velocity, density and pressure are computed by solving the Navier–Stokes equations and by the interpolation of the values of neighboring particles. The contribution of each of the neighboring particles depends on the inter-particle distance and on the kernel function W, which has an area of influence defined by the smoothing length h (see [14,18,23,35]).

2.1. Governing Equations

The Navier–Stokes equations written in their SPH notation are solved each timestep for each of the particles. The momentum and continuity equations are given in Equations (1) and (2), respectively. In the following equations, the physical properties of particle a are calculated, with b representing each of its neighboring particles:

$$\frac{d{\mathit{v}}_{\mathit{a}}}{dt}=-\sum _b{m}_b\left(\frac{P_b+P_a}{{\rho }_b{\rho }_a}+{\Pi }_{ab}\right){\nabla }_a{W}_{ab}+\mathit{g}$$

(1)

$$\frac{d{\rho }_a}{dt}=\sum _b{m}_b({\mathit{v}}_{\mathit{a}}-{\mathit{v}}_{\mathit{b}}){\nabla }_a{W}_{ab}+\delta h{c}_a\sum _b{\psi }_{ab}·{\nabla }_a{W}_{ab}\frac{m_b}{{\rho }_b}$$

(2)

where m is mass, $\mathit{v}$ is velocity, P is pressure, $\rho$ is density and $\mathit{g}$ is the gravity acceleration. $W_{ab}$ represents the kernel function and depends on the distance between particles a and b. In this work, a Quintic kernel [36] is applied, since this type of kernel is well suited for general free-surface problems (see [19]). One of the main advantages of the Quintic kernel is that the tensile instability that appears using other kernels, such as the Cubic Spline, is avoided using the kernel adopted in this work. More information can be found in [37]. The diffusion term introduced on the right-hand side of the continuity equation (Equation (2)) acts as a numerical noise filter, thereby improving the numerical stability and smoothing the density and pressure field (see [38,39,40]). In this paper, the recently proposed density diffusion term of [41] is applied, since it has been proven to produce more accurate results for the pressure field near boundaries while keeping the computational cost limited. The term ${\psi }_{ab}$ takes the following form:

$${\psi }_{ab}=2({\rho }_b^D-{\rho }_a^D)\frac{{\mathit{r}}_{\mathit{ab}}}{\Vert {\mathit{r}}_{\mathit{ab}}\Vert }$$

(3)

where ${\mathit{r}}_{\mathit{ab}}={\mathit{r}}_{\mathit{a}}-{\mathit{r}}_{\mathit{b}}$ with ${\mathit{r}}_{\mathit{k}}$ being the position of particle k and ${\rho }^D$ being the dynamic density, equal to the difference of the total (${\rho }^T$) and hydrostatic (${\rho }^H$) density: ${\rho }^D={\rho }^T-{\rho }^H$.

${\Pi }_{ab}$ represents the artificial viscosity term, proposed in [14]:

$${\Pi }_{ab}=\left\{\begin{array}{cc}\frac{-\alpha \overline{c_{ab}}{\mu }_{ab}}{\overline{{\rho }_{ab}}}\hfill & \mathrm{if}{\mathit{v}}_{\mathit{ab}}·{\mathit{r}}_{\mathit{ab}}<0\hfill \\ 0\hfill & \mathrm{if}{\mathit{v}}_{\mathit{ab}}·{\mathit{r}}_{\mathit{ab}}>0\hfill \end{array}\right.$$

(4)

where ${\mathit{v}}_{\mathit{ab}}={\mathit{v}}_{\mathit{a}}-{\mathit{v}}_{\mathit{b}}$, ${\mathit{v}}_{\mathit{k}}$ is the velocity of particle k, $m{u}_{ab}=\frac{h{\mathit{v}}_{\mathit{ab}}·{\mathit{r}}_{\mathit{ab}}}{r_{ab}^2+{\eta }^2}$, $\overline{c_{ab}}=0.5(c_a+c_b)$ is the mean speed of sound, ${\eta }^2=0.1h^2$ and $\alpha$ is a coefficient tuned for proper dissipation (see [18]). The artificial viscosity term, ${\Pi }_{ab}$, is added in the momentum equation based on the Neumann–Richtmeyer artificial viscosity with the aim of reducing oscillations and stabilizing the SPH scheme, following the work in [14]. DualSPHysics code uses this numerical viscosity; however, more physical treatments need to be implemented in the future in order to properly identify the turbulence by using the K-Epsilon model or Sub-Particle Scale model. Throughout this paper, $\alpha$ is set to be equal to 0.01 (as in [21,42]) or 0.00 (in the case of inviscid simulations). Since the fluid is weakly-compressible in DualSPHysics, an equation of state is used to calculate the fluid pressure as a function of the density instead of solving a Poisson-like equation. The speed of sound $c_s$ is artificially lowered such that a reasonable time step can be employed while ensuring that density variations are kept lower than 1% during the simulation.

$$P_a=\frac{c_s^2{\rho }_0}{\gamma }\left[{\left(\frac{{\rho }_a}{{\rho }_0}\right)}^{\gamma }-1\right]$$

(5)

where $\gamma =7$ is the polytropic constant and ${\rho }_0$ is the reference fluid density.

2.2. Boundary Conditions

In DualSPHysics, the boundaries are described by a set of particles for which the same equations (Equations (1) and (2)) as used for fluid particles are solved. However, the particles belonging to the boundary do not move according to the forces acting on them: boundary particles remain either fixed or move according to a predefined movement. These boundary conditions are called Dynamic Boundary Conditions (DBC) and have the advantage of being able to deal with complex geometries and being computationally efficient. However, due to excessive repulsive forces near the boundary between a structure and the fluid, a gap appears of the order of magnitude of the smoothing length h (see [43]). Furthermore, at this same boundary, unphysical values of pressure are observed. These issues are dealt with in the modified DBC (mDBC) implementation recently added to DualSPHysics, as described in [43]. For each boundary particle of mDBC, a ghost node is mirrored into the fluid. Boundary particles then receive the properties of the fluid at the position of the ghost node. The density is calculated using a first-order consistent SPH interpolation, as proposed in [44]. The mDBC method has proven to show results with more realistic physical values for pressure at the boundaries and a significant reduction in the size of the gap between the fluid and boundary without a significant extra computational cost.

2.3. Floating Bodies

DualSPHyics has the capability to accurately simulate fluid-driven objects, as described and validated in [23,45], and is used extensively in this paper. The force per unit mass acting on one boundary particle k of the floating body is calculated by summing the contributing forces exerted on this boundary particle k by fluid particles a within the compact support of the kernel:

$${\mathit{f}}_k=\sum _a{\mathit{f}}_{ka}$$

(6)

where ${\mathit{f}}_k$ is the force per unit mass acting on boundary particle k of the floating body and ${\mathit{f}}_{ka}$ is the force per unit mass acting on boundary particle k exerted by fluid particle a, calculated with Equation (1) [46]. Once the force acting on the floating body is computed, the body’s motion can be determined, assuming it is rigid:

$$M\frac{d\mathit{V}}{dt}=\sum _k{m}_k{\mathit{f}}_k$$

(7)

$$I\frac{d\mathsf{\Omega }}{dt}=\sum _k{m}_k({\mathit{r}}_k-\mathit{R})\times {\mathit{f}}_k$$

(8)

where M is the body’s total mass, I is the moment of inertia, $\mathit{V}$ is the velocity, $\mathsf{\Omega }$ is the rotational velocity, $\mathit{R}$ is the centre of mass, $m_k$ is the mass of floating boundary particle k and ${\mathit{r}}_k$ is the position of floating particle k. Equations (7) and (8) are integrated over time in order to compute the velocity ${\mathit{v}}_k$ of the floating particle:

$${\mathit{v}}_k=\mathit{V}+\mathsf{\Omega }\times ({\mathit{r}}_k-\mathit{R})$$

(9)

2.4. Coupling Dualsphysics—Project Chrono

Project Chrono is an open-source software package that enables the numerical modeling of mechanical constraints and collisions between objects [20]. It has recently been successfully coupled to DualSPHysics [19,47] and is used in this case for the modeling of the PTO system of a WEC.

Linear damping was already implemented in the coupling, which is applied here for the simulation with a linear damping PTO system (see Equation (10)).

$$F_{PTO,l}\left(t\right)=-B_{PTO,l}v\left(t\right)$$

(10)

where $B_{PTO,l}$ is the linear PTO system damping coefficient and $v\left(t\right)$ is the WEC’s heave velocity. The DualSPHysics–Project Chrono model is here extended with the Coulomb damping model, applying the force described in Equation (11), where $B_{PTO,c}$ is the PTO damping coefficient for a Coulomb damping PTO system.

$$F_{PTO,c}\left(t\right)=-B_{PTO,c}sign\left(v\left(t\right)\right)$$

(11)

where $B_{PTO,c}$ is the Coulomb damping coefficient.

3. Methodology

In this study, both linear potential flow theory as well as the SPH method are applied. Therefore, this section begins with a brief overview of the linear potential flow theory. Next, we describe how the hydrodynamic coefficients can be computed using DualSPHysics. Finally, we describe how the equation of motion based on linear potential flow theory was modified to take into account the effect of the drag force. A schematic overview of the applied theories and equations used in the current study for the calculation of the motion and the average absorbed power of the studied WECs is given in Figure 1.

In the calculation of the hydrodynamic coefficients and the equation of motion in the present study, we assume linear potential flow theory [6], which is a subset of linear wave theory that allows the fluid velocity, $v_w$, to be expressed as the gradient of the time dependent velocity potential $\Phi$ (Equation (12)).

$$v_w=\nabla \Phi$$

(12)

Potential flow theory makes some essential assumptions: (i) the fluid is inviscid, (ii) the fluid is incompressible and (iii) the flow is irrotational [3]. Furthermore, it is assumed that the motion amplitude of the WEC is much smaller than the wavelength. It is assumed that all time-dependent variables oscillate with the same wave angular frequency $\omega$, enabling the separation of the time dependence from the time-independent velocity potential $\Phi$,

$$\Phi (x,y,z,t)=\Re \left\{\widehat{\Phi }(x,y,z)e^{i\omega t}\right\}$$

(13)

where $\widehat{\Phi }$ is the complex amplitude of the velocity potential. Due to application of the principle of superposition, linear potential flow theory allows the separation of the total velocity potential into the following components (Equation (14)):

$${\Phi }_t(x,y,z)={\Phi }_i+{\Phi }_d+\sum _i^6{\Phi }_r$$

(14)

where ${\Phi }_t$ is the total velocity potential, ${\Phi }_i$ is the incident wave velocity potential, ${\Phi }_d$ is the diffracted wave velocity potential and ${\sum }_i^6{\Phi }_r$ is the sum of the radiated wave velocity potentials for each degree of freedom (DoF) of the WEC. In the current study, only one DoF for each WEC is modeled: namely the heave motion. Since the oscillatory motion of the waves and the WEC is considered to be harmonic, all forces and displacements can be decomposed into their spatial and temporal dependencies [6]. Therefore, the WEC’s heave displacement can be written as

$$\xi \left(t\right)=\Re \left\{\widehat{\xi }\left(\omega \right)e^{i\omega t}\right\}$$

(15)

where $\widehat{\xi }$ is the complex amplitude of the WEC’s heave displacement. Applying Newton’s second law of motion and writing all terms in its complex amplitudes results in the following equation (see [6]):

$$-{\omega }^{2}m\widehat{\xi }={\widehat{F}}_e+{\widehat{F}}_r+{\widehat{F}}_{hs}+{\widehat{F}}_{PTO}$$

(16)

where ${\widehat{F}}_e$ is the complex amplitude of the excitation force, ${\widehat{F}}_r$ is the complex amplitude of the radiation force, ${\widehat{F}}_{hs}$ is the complex amplitude of the hydrostatic force and ${\widehat{F}}_{PTO}$ is the complex amplitude of the PTO system force. It is here assumed that the PTO system is a linear damping PTO system, following Equation (10). Equation (16) can be further rewritten, following the same approach as in [6]:

$$-{\omega }^{2}m\widehat{\xi }={\widehat{F}}_e-i\omega B_{33}\widehat{\xi }+{\omega }^2{A}_{33}\widehat{\xi }-K_{33}\widehat{\xi }-i\omega B_{PTO,l}\widehat{\xi }$$

(17)

where $B_{33}$ is the heave component of the hydrodynamic damping, $A_{33}$ is the heave component of the added mass, $K_{33}$ is the hydrostatic spring stiffness or hydrostatic coefficient equal to $K_{33}=\rho g{A}_d$ with $A_d$ being the cross-sectional area at the undisturbed sea level, and $B_{PTO,l}$ is the PTO system damping coefficient for a linear PTO system without a spring coefficient. From Equation (17), the complex amplitude of the device’s heave motion can be obtained,

$$\widehat{\xi }=\frac{{\widehat{F}}_e}{-{\omega }^2(m+A_{33})+K_{33}+i\omega (B_{33}+B_{PTO,l})}$$

(18)

Equation (18) is further extended to include the effect of the drag force in Section 3.2. Equation (18) can be used to calculate the average absorbed power $P_{av}$ [6]:

$$P_{av}=\frac{1}{T}{\int }_0^T{B}_{PTO,l}v{\left(t\right)}^{2}dt=\frac{1}{2}{B}_{PTO,l}{\omega }^2{\left|\widehat{\xi }\right|}^2$$

(19)

An optimal value for $B_{PTO,l}$ can be calculated, leading to the maximal average absorbed power [48]:

$$B_{PTO,l,opt}=\sqrt{B_{33}^2+{\left(\omega (m+A_{33})-\frac{K_{33}}{\omega }\right)}^2}$$

(20)

3.1. Determination of Hydrodynamic Coefficients Using DualSPHysics

The added mass, hydrodynamic damping and drag coefficient are estimated in DualSPHysics by running a forced motion test. The test consists of the WEC heaving according to a fixed sinusoidal motion in still water. The heave motion and heave velocity of the WEC can be described as follows:

$$z\left(t\right)=asin\left(\omega t\right)$$

(21)

$$v\left(t\right)=a\omega cos\left(\omega t\right)$$

(22)

where a is the heave amplitude of the WEC with a forced heave motion.

Following potential flow theory, the total vertical force $F_z$ acting on the heaving WEC, as the sum of the hydrodynamic, hydrostatic, drag and PTO system force, is expressed as

$$F_z\left(t\right)=(A_{33}{\omega }^2-\rho g{A}_d)asin\left(\omega t\right)-B_{33}\omega acos\left(\omega t\right)+F_v\left(t\right)$$

(23)

where $F_v\left(t\right)$ is the viscous force, written as the semi-empirical Morison equation [49,50]:

$$F_v\left(t\right)=-\frac{1}{2}\rho A_d{C}_{d}v\left(t\right)\left|v\left(t\right)\right|=-\frac{1}{2}\rho A_d{C}_d{\left(\omega a\right)}^{2}cos\left(\omega t\right)\left|cos\left(\omega t\right)\right|$$

(24)

where $C_d$ is the drag coefficient. The viscous force $F_v\left(t\right)$ can be approximated by writing out its Fourier series and retaining only the first term: Equation (25). It is preferred to retain only the first term since this allows the modification of the equation of motion depending on only one single frequency. This is justified since the second appearing term of the drag force’s Fourier series has a magnitude five times lower than the first frequency component. The same approach was used in [51].

$$F_v\left(t\right)\approx -\frac{4}{3\pi }\rho A_d{C}_d{\left(\omega a\right)}^{2}cos\left(\omega t\right)$$

(25)

Now, the force $F_z\left(t\right)$ from Equation (23) can be written as a sum of a sine and a cosine, containing only one frequency component:

$$F_z\left(t\right)=(A_{33}{\omega }^2-\rho g{A}_d)asin\left(\omega t\right)+(-B_{33}\omega a-\frac{4}{3\pi }\rho A_d{C}_d{\left(\omega a\right)}^2)cos\left(\omega t\right)$$

(26)

One way of estimating the hydrodynamic coefficients is to apply the least-square method when comparing the theoretical total force with Equation (23) with the total force acting on the WEC, calculated with DualSPHysics, $F_{SPH}\left(t\right)$. However, since multiple coefficients have to be estimated, we chose to apply Fourier analysis. A similar approach was followed in [25,26,28,52]. The coefficients of the Fourier series of the force $F_z\left(t\right)$ (Equation (26)) acting on the heaving WEC can be used to calculate the hydrodynamic coefficients.

$$b_1=\frac{2}{T}{\int }_0^T{F}_{SPH}\left(t\right)sin\left(\omega t\right)$$

(27)

$$A_{33}=\frac{b_1+\rho g{A}_{d}a}{{\omega }^{2}a}$$

(28)

Equation (28) assumes that the cross-sectional area of the WEC, $A_d$, is constant. This is true for a cylinder, but not for, e.g., a sphere, where the varying cross-sectional area causes non-linear static Froude–Krylov forces (see [32,53]). The two main non-linearities affecting point absorber WECs in the resonance region are the drag forces and the non-linear Froude–Krylov forces, the latter being a purely geometrical effect [32,53,54]. Therefore, the heave amplitude of the spherical WEC was kept low in order to limit these non-linear Froude–Krylov forces.

$$a_1=\frac{2}{T}{\int }_0^T{F}_{SPH}\left(t\right)cos\left(\omega t\right)$$

(29)

$$B_{33}=-\frac{a_1}{\omega a}$$

(30)

Equation (30) is only applicable when the fluid is inviscid. The first simulations were performed in an inviscid fluid by setting the artificial viscosity coefficient in DualSPHysics equal to zero: $\alpha =0.0$. This does not lead to numerical instabilities of the water particles at the water surface since the density-diffusion term is applied [41]. Once the inviscid simulations are finished, the derived hydrodynamic damping coefficient $B_{33}$ can be used in viscous simulations to estimate the drag coefficient $C_d$:

$$C_d=-\frac{a_1+B_{33}\omega a}{\frac{4}{3\pi }\rho A_d{\omega }^2{a}^2}$$

(31)

The test cases discussed in this paper all have a high Reynolds number ($Re$), meaning that the drag coefficient mainly depends on the pressure distribution surrounding the WEC and not on the distribution of the wall shear stresses (see [55] (Chapter 9)). In order to obtain an accurate pressure distribution, it is essential to apply the mDBC method instead of DBC. It is described in [53] that the total drag can be separated into a scale-independent and a scale-dependent (or viscous) drag component. Viscous effects depend on the Reynolds number and do not scale with Froude scaling. On the other hand, there is form drag and vortex-induced drag, which is independent of scale (Figure 2). Since artificial viscosity is used in the DualSPHysics simulations throughout this paper, it is expected that the viscous drag component is not modeled correctly. However, it follows from [53] that the viscous drag term only covers a small part of the total drag (1–4%), whereas the vortex-induced drag covers 18–19% of the total drag. Vortex shedding is generated by high pressure and velocity gradients close to regions of high wall curvature [53]. Vorticity is modeled in DualSPHysics since the software supports flow rotationality and there is numerical (artificial) viscosity. The presence of these vortices is proven later in the present study (see Figure 3). Vortices surrounding the heaving WEC disrupt the dynamic pressure field, causing a loss in heave amplitude. Furthermore, in [4], it was concluded that pressure or form drag is the main contributor of drag forces, caused by flow separation and vortex shedding, whereas skin friction drag was considered negligible.

3.2. Derivation of the Equation of Motion and the Optimal Damping Coefficient Including the Effect of Drag Forces

After the drag coefficient has been determined with DualSPHysics, it can be used in potential flow theory. The aim of this section is to derive a modified equation of motion including the effect of the drag force. Furthermore, a modified equation for the value of $B_{PTO,l,opt}$ is derived, since it has been proven that the drag coefficient has an influence on the optimal linear PTO damping coefficient of a heaving WEC (see [34,56]).

3.2.1. Linear Damping PTO System

The PTO force from a linear damping PTO system is written as follows in the time domain:

$$F_{PTO,l}\left(t\right)=-B_{PTO,l}v\left(t\right)$$

(32)

The complex amplitude of the linear damping PTO system force is given by

$${\widehat{F}}_{PTO,l}\left(\omega \right)=-i\omega B_{PTO,l}\widehat{\xi }$$

(33)

In order to obtain a modified equation of motion and a new estimate for the optimal damping coefficient taking into account the effect of the drag force, defined in Equation (24), this force has to be written in the frequency domain. In a first step, it is assumed the WEC is oscillating in still water. The drag force contains a quadratic term, which can be approximated as follows by calculating its Fourier series and retaining only the first frequency component (see also [51]):

$$v\left(t\right)=\Re \left\{i\omega \widehat{\xi }{e}^{i\omega t}\right\}$$

(34)

$$\Rightarrow v\left(t\right)=-\omega |\widehat{\xi }|sin\left(\omega t\right)$$

(35)

$$-v\left(t\right)·\left|v\left(t\right)\right|\approx -\frac{8}{3\pi }\left(\omega \right|\widehat{\xi }\left|\right)v\left(t\right)$$

(36)

$$\Rightarrow -v\left(t\right)·\left|v\left(t\right)\right|\approx \Re \left\{-\frac{8}{3\pi }\left(\omega \right|\widehat{\xi }\left|\right)i\omega \widehat{\xi }{e}^{i\omega t}\right\}$$

(37)

Filling in the result of Equation (37) in the drag force gives

$$F_v\left(t\right)=-\frac{1}{2}\rho A_d{C}_{d}v\left(t\right)·\left|v\left(t\right)\right|$$

(38)

$$\Rightarrow F_v\left(t\right)=\Re \left\{-\frac{4}{3\pi }\rho A_d{C}_d\left(\omega \right|\widehat{\xi }\left|\right)i\omega \widehat{\xi }{e}^{i\omega t}\right\}$$

(39)

Equation (39) results in the following expression for the complex amplitude of the drag force:

$${\widehat{F}}_v=-\frac{4}{3\pi }\rho A_d{C}_d\left(\omega \right|\widehat{\xi }\left|\right)i\omega \widehat{\xi }$$

(40)

$$\Rightarrow {\widehat{F}}_v=-{\widehat{c}}_{f}i\omega \widehat{\xi }$$

(41)

In this case, ${\widehat{c}}_f$ is a real number; this will not be the case when the fluid is no longer still. Equation (41) can now used in the equation of motion:

$$-{\omega }^{2}m\widehat{\xi }={\widehat{F}}_e+{\widehat{F}}_r+{\widehat{F}}_{hs}+{\widehat{F}}_{PTO}+{\widehat{F}}_v$$

(42)

$$\Rightarrow -{\omega }^{2}m\widehat{\xi }={\widehat{F}}_e+{\omega }^2{A}_{33}\widehat{\xi }-i\omega B_{33}\widehat{\xi }-K_{33}\widehat{\xi }-i\omega B_{PTO,l}\widehat{\xi }-i\omega {\widehat{c}}_f\widehat{\xi }$$

(43)

$$\widehat{\xi }=\frac{\widehat{F_e}}{(-{\omega }^2(m+A_{33})+K_{33})+i\omega (B_{33}+B_{PTO,l}+{\widehat{c}}_f)}$$

(44)

Equation (38) is, however, only true if the water surrounding the WEC is standing still, which is no longer the case when the WEC is moving in waves. In the formula of the drag force, the relative velocity difference $v^{\prime }\left(t\right)=v\left(t\right)-v_0\left(t\right)$ has to be used, which is the difference between the WEC’s velocity $v\left(t\right)$ and the vertical velocity of the surrounding water particles $v_0\left(t\right)$ (see Equation (45)).

$$v_0\left(t\right)=-\frac{H\pi }{T}sin\left(\omega t\right)e^{-k{z}_b}$$

(45)

where H is the wave height, k is the wave number and $z_b$ is the depth at which the vertical fluid particle velocity is calculated, taken as half of the draft in current research.

The relative velocity difference is written as a function of the WEC’s velocity $v\left(t\right)$ in order to simplify further results:

$$v^{\prime }\left(t\right)=\Re \left\{{\alpha }_{v}i\omega \widehat{\xi }{e}^{i{\varphi }_{\alpha }}{e}^{i\omega t}\right\}$$

(46)

where ${\alpha }_v$ is the ratio of the amplitude of $v^{\prime }\left(t\right)$ to the amplitude of $v\left(t\right)$ and ${\varphi }_{\alpha }$ is the phaseshift between $v^{\prime }\left(t\right)$ and $v\left(t\right)$. Harmonic decomposition is applied in a similar manner as in Equation (36):

$$-v^{\prime }\left(t\right)·|v^{\prime }\left(t\right)|\approx -\frac{8}{3\pi }({\alpha }_v\omega |\widehat{\xi }\left|\right)·v^{\prime }\left(t\right)$$

(47)

$$\Rightarrow -v^{\prime }\left(t\right)·\left|v^{\prime }\left(t\right)\right|\approx \Re \left\{-\frac{8}{3\pi }({\alpha }_v\omega |\widehat{\xi }\left|\right){\alpha }_v{e}^{i{\varphi }_{\alpha }}i\omega \widehat{\xi }{e}^{i\omega t}\right\}$$

(48)

The Equation (48) is now used to represent the viscous force in the time and frequency domain.

$$F_v\left(t\right)=-\frac{1}{2}\rho A_d{C}_d{v}^{\prime }\left(t\right)\left|v^{\prime }\left(t\right)\right|=\Re \left\{-\frac{4}{3\pi }\rho A_d{C}_d{\alpha }_v^2\left(\omega \right|\widehat{\xi }\left|\right)e^{i{\varphi }_{\alpha }}i\omega \widehat{\xi }{e}^{i\omega t}\right\}$$

(49)

$$\Rightarrow \widehat{F_v}=-\frac{4}{3\pi }\rho A_d{C}_d{\alpha }_v^2\left(\omega \right|\widehat{\xi }\left|\right)e^{i{\varphi }_{\alpha }}i\omega \widehat{\xi }$$

(50)

$$\Rightarrow \widehat{F_v}=-{\widehat{c}}_f·i\omega \widehat{\xi }$$

(51)

with

$${\widehat{c}}_f=\frac{4}{3\pi }\rho A_d{C}_d{\alpha }_v^2\left(\omega \right|\widehat{\xi }\left|\right)e^{i{\varphi }_{\alpha }}$$

(52)

Equation (50) is the complex amplitude of the viscous force written as a function of the complex amplitude of the WEC’s heave motion $\widehat{\xi }$. This force can be introduced in the equation of motion of the WEC:

$$\widehat{\xi }=\frac{\widehat{F_e}}{(-{\omega }^2(m+A_{33})+K_{33}-\Im \left({\widehat{c}}_f\right)\omega +i\omega (B_{33}+B_{PTO,l}+\Re \left({\widehat{c}}_f\right)))}$$

(53)

Note that ${\alpha }_v$ and ${\varphi }_{\alpha }$ are not known in advance since they depend on $\widehat{\xi }$, so an iterative procedure has to be followed to solve Equation (53).

Once a solution for $\widehat{\xi }$ is found, ${\widehat{c}}_f$ can be calculated from Equation (52), and an updated value of the optimal damping coefficient of the heaving WEC can be calculated:

$$B_{PTO,l,opt}=\sqrt{{(B_{33}+\Re \left({\widehat{c}}_f\right))}^2+{\left(\omega (m+A_{33})-\frac{K_{33}-\Im \left({\widehat{c}}_f\right)\omega }{\omega }\right)}^2}$$

(54)

3.2.2. Coulomb Damping PTO System

In case a hydraulic PTO system is applied, the PTO system force can be approximated as a Coulomb damping system [48]:

$$F_{PTO,c}\left(t\right)=-B_{PTO,c}sign\left(v\left(t\right)\right)$$

(55)

where $B_{PTO,c}$ is the Coulomb damping coefficient. The first frequency component of the Fourier series of this PTO system force can be used as approximation in the frequency domain (see Equation (56)). ${\widehat{F}}_{PTO,c}$ is written in such a way that it resembles the complex amplitude of the linear PTO system force (see Equation (33)):

$${\widehat{F}}_{PTO,c}=-\frac{4}{\pi \omega |\widehat{\xi }|}{B}_{PTO,c}i\omega \widehat{\xi }$$

(56)

The second frequency component, at $3\omega$, has an amplitude of one-third of the first frequency component. However, it can be proven that this second-frequency component (and all other higher-order frequency components) does not change the average absorbed power of the WEC (see also later in Section 5). Therefore, only the first frequency component of $F_{PTO,c}\left(t\right)$ is retained below. This now means that $F_{PTO,c}\left(t\right)$ can be written as the linear damping PTO system force with an equivalent damping coefficient $B_{PTO,l,eq}$:

$$B_{PTO,l,eq}=\frac{4}{\pi \omega |\widehat{\xi }|}{B}_{PTO,c}$$

(57)

This $B_{PTO,l,eq}$ is then used in Equation (53). Since $B_{PTO,l,eq}$ depends on the WEC’s heave amplitude $|\widehat{\xi }|$, an iterative approach is followed to solve the Equation (53). The value for $B_{PTO,c}$ leading to the maximum average absorbed power $P_{av}$, $B_{PTO,c,opt}$, can then be estimated as

$$B_{PTO,c,opt}=\frac{\pi }{4}\omega \left|\widehat{\xi }\right|·B_{PTO,l,opt}$$

(58)

with $B_{PTO,l,opt}$ from Equation (54).

4. Test Cases and Numerical Setup

The methods used to determine $C_d$ described in Section 3.1 are applied to three different WECs:

• A heaving sphere with a diameter of 5 m, as studied in [32]. For this WEC, the added mass and hydrodynamic damping were also estimated wtih DualSPHysics.
• A cylindrical WEC with a diameter of 0.5 m and draft of 0.11 m, as studied in [17,56], hereafter referred to as “cylinder1”.
• A cylindrical WEC with a diameter of 0.3 m and draft of 0.28 m, as studied in [34,57], hereafter referred to as ’cylinder2’. This WEC is studied in more detail by analyzing its response amplitude operator (RAO) and comparing it with the experimental results found in [57]. The WEC is also simulated with two kinds of PTO systems (a linear damping and a Coulomb damping PTO system) and the average absorbed power is compared for a range of damping coefficients.

Spherical and cylindrical shapes of the WECs were chosen since these are the shapes that are regularly used for heaving point absorbers. A spherical WEC was chosen since this shape induces extra non-linear forces, named non-linear Froude—Krylov forces, due to their non-uniform cross-sectional area [54]. Two cylindrical shapes were chosen since a flat cylinder provides a balance between the power absorption, WEC bandwidth and material cost considerations (see [58]), whereas for a slender cylinder, the drag force is expected to be more significant [34].

An overview of important settings and parameters applied during the DualSPHysics simulations is listed below.

• Simulations were carried out with two different types of boundary conditions: Dynamic Boundary Conditions (DBC) and modified Dynamic Boundary Conditions (mDBC), as described in Section 2.2. Both SPH results were compared to the theoretical force calculated with Equation (23), with the hydrodynamic coefficients from NEMOH and $C_d=0.45$ as in [32]. It is clear from Figure 4 that mDBC gave significantly better results compared to DBC. When mDBC was applied, the repulsive forces exerted by the boundary particles of the sphere were much smaller than when applying DBC, resulting in a smaller gap between the sphere and the fluid.
• Artificial viscosity was applied with an artificial viscosity coefficient $\alpha =0.01$. Artificial viscosity was introduced into SPH in [14] and was used primarily due to its simplicity [18]. It was stated in [59] that this artificial viscosity corresponds to an equivalent kinematic viscosity of $\frac{15}{112}\alpha c_{s0}h$ (in 2D), which is generally much higher than the real kinematic viscosity of water $\nu$ = 1 × 10${}^{-6}$ m²/s. Therefore, one way of reducing the numerical dissipation caused by artificial viscosity is by lowering $\alpha$; however, this was not preferred since the value of $\alpha =0.01$ has been proven to give the best results in the validation of wave flumes to study the wave propagation and wave loadings exerted onto coastal structures [21,42] and is also the value used when simulating the WEC in regular waves. Only in the case where the hydrodynamic coefficients $A_{33}$ and $B_{33}$ are computed is $\alpha$ set to be equal to zero, as described in Section 3.1.
• The initial speed of sound was set to $c_{s0}=15\sqrt{gd}$, with d being the depth of the numerical wave basin. It was found that convergence was reached with a lower resolution when the speed of sound $c_{s0}$ was decreased. This can be related to the influence of $c_{s0}$ on the viscosity: it is stated in [59] that the equivalent kinematic viscosity associated with the artificial viscous term has the form $\frac{15}{112}\alpha c_{s0}h$. Further decreasing $c_{s0}$ leads to overly large timesteps and therefore less accurate results.
• A convergence test was done by varying the interparticle distance $dp$ and studying the resulting hydrodynamic coefficients $A_{33}$, $B_{33}$ and the drag coefficient $C_d$, calculated with Equations (28), (30) and (31), respectively. The hydrodynamic coefficients were compared to results from potential flow theory obtained with NEMOH and the drag coefficient was compared to results from previous experimental or numerical tests from [32,33,34]. The results of these convergence tests are described in Section 5.
• The domain size of the basin was set to be large enough to avoid interaction with side walls (see Figure 5). Sloped sidewalls were provided as well as numerical damping layers, with the aim of reducing side wall reflection.

5. Results and Discussion

5.1. Estimation of the Hydrodynamic Coefficients and of the Drag Coefficient

Before calculating the drag coefficient $C_d$, wer briefly checked whether DualSPHysics was able to accurately compute the hydrodynamic coefficients corresponding to the added mass $A_{33}$ and the hydrodynamic damping $B_{33}$. This computation was carried out for the heaving sphere described in Section 4 by simulating this heaving sphere in an inviscid fluid ($\alpha =0.0$) in DualSPHysics. Equations (28) and (30) were applied in order to calculate $A_{33}$ and $B_{33}$. The obtained values were compared to the values found by NEMOH, which applies linear potential flow theory. In order to obtain results comparable to linear theory, the heave amplitude a of the sphere was kept rather low, at $a=0.25$ m. Finally, the results of $A_{33}$ and $B_{33}$ are displayed in Figure 6, showing that DualSPHysics was able to accurately compute $A_{33}$ and $B_{33}$, apart from some discrepancies in the low-frequency region.

In the following simulations, carried out to compute the drag coefficient $C_d$, $\alpha$ was set to be equal to 0.01 as suggested in [21]. As discussed in Section 3.1, DualSPHysics supports flow rotationality and adds artificial viscosity, resulting in vortex shedding. Vortex shedding is expected to be present in regions with a high wall curvature, such as the corners of a heaving cylindrical WEC, and results in a significant part of the total drag force. Therefore, a brief qualitative analysis was carried out to investigate the presence of vorticity in DualSPHysics simulations. The cylindrical WEC with a diameter if 0.5 m and draft of 0.11 m—cylinder1—was simulated in DualSPHysics in regular waves with $H=0.16$ m, $T=1.5$ s and $B_{PTO,l}=1100$ Ns/m, in a similar way as was done in an experimental study in [56]. The high value of $B_{PTO,l}$ was chosen in order to induce a large phase difference and thus high relative velocities and high vortex-induced drag. A qualitative comparison between the numerical and experimental model is displayed in Figure 3, showing the magnitude of vorticity calculated in DualSPHysics and the flow patterns surrounding the WEC during experiments. It is clear that DualSPHysics shows high vorticity at the sides and corners of the heaving WEC, similar to the results in the experimental case studied in [56].

The presence of vortex shedding around the heaving cylindrical WEC in DualSPHysics simulations implied that there was a drag force acting on the WEC. In the remainder of this section, we determine whether this drag force was accurately modeled in DualSPHysics by calculating the drag coefficient $C_d$, applying Equation (31). This $C_d$ was calculated for the three WECs described in Section 4: one spherical WEC and two cylindrical WECs. In all cases, the values of $A_{33}$ and $B_{33}$ used in Equation (31) were those obtained from NEMOH, which is the same approach as followed in [32]. This means that $C_d$ does not only account for the viscous non-linear effects, but also for other non-linear effects, as also described in [29]. Figure 7 shows the results for the calculation of $C_d$ using the approach described in Section 3.1, for different ratios of $D/dp$, with D being the WEC’s diameter and $dp$ the interparticle distance. The convergence analyses are performed for the sphere studied in [32] with $a=1.5$ m and $T=9$ s, the cylinder studied in [17,56], with $a=0.045$ m and $T=1.5$ s (cylinder1 in Figure 7) and the cylinder studied in [34,60], with $a=0.1$ m and $T=1.2$ s (cylinder2 in Figure 7). The $C_d$ value obtained from previous numerical tests or experiments performed in [32,33,34] is displayed as well; it is clear that the result of DualSPHysics converges towards a slightly higher value of $C_d$ than obtained in previous numerical or experimental tests. One possible explanation for why the DualSPHysics simulations slightly overestimate the value of $C_d$ is the use of the artificial viscosity value with $\alpha =0.01$, resulting in excessive equivalent kinematic viscosity as described in Section 4. In future work, a more physical viscosity and turbulence model will be implemented in DualSPHysics. The $C_d$ value obtained for cylinder1 simulated in DualSPHysics was compared to the $C_d$ obtained from the results described in [33]; i.e., $C_d=1.5$. In [33,61], experimental and numerical tests were carried out for a heaving vertical cylinder with similar values for the Keulegan–Carpenter number $KC$ and the Reynolds number $Re$. The $C_d$ value obtained for cylinder2 in DualSPHysics was compared to the $C_d$ value determined experimentally in [34]; i.e., $C_d=1.4$. In [34,60], the drag coefficient $C_d$ was determined by performing one or multiple heave decay tests for the cylinder. However, this approach was not followed here since decay tests in DualSPHysics require a very high resolution for decent accuracy (see [62]). The test for cylinder2 performed in DualSPHysics was carried out at a period of $T=1.2$ s and a heave amplitude $a=0.1$ m, resulting in values for $KC$ and $Re$ of the same order of magnitude as those in the free decay tests, which is important since $C_d$ depends on both $KC$ and $Re$ [63].

An overview of the heave amplitude a and the heave period T applied during the DualSPHysics tests is displayed in

Table 1. Heave amplitude a and heave period T applied during the forced oscillation with DualSPHysics of a spherical WEC and two cylindrical WECs, as well as the obtained drag coefficient $C_d$.

	Spherical WEC	Cylinder1	Cylinder2
T [s]	9	1.5	1.2
a [m]	1.5	0.045	0.1
${\mathit{C}}_{\mathit{d}}$ [-]	0.78	1.65	1.50

, as well as the obtained value of $C_d$.

5.2. Cylindrical WEC in Regular Waves

Once the $C_d$ values were determined, a study was carried out on the WEC’s behavior in regular waves. This study was carried out only for cylinder2, since experimental validation for this WEC was available from [34]. The modified equation of motion including the effect of the drag force, Equation (53), was used for the calculation of the heave motion and the average absorbed power. The $C_d$ used in Equation (53) equalled 1.5, determined with DualSPHysics (see Table 1). The results of the heave motion and average absorbed power were then compared with the results obtained from simulations in DualSPHysics.

5.2.1. Undamped Heaving WEC

Before studying the heaving WEC with a PTO system, the response amplitude operator (RAO) of the undamped heaving WEC was calculated. This RAO could be calculated by comparing different modeling techniques: (i) linear potential flow theory—Equation (18), (ii) linear potential flow theory with a correction term taking into account the effect of the drag force—Equation (53), (iii) SPH simulations using DualSPHysics and (iv) experimental results using data from [34,57]. This RAO was calculated for the undamped case ($B_{PTO,l}=0.0$ Ns/m) with $H=0.08$ m as was done in the experiments. The dimensions of the numerical wave basin modeled in DualSPHysics differed for each wave frequency, since the WEC should be placed at least one wavelength away from the piston [17]. An example of the numerical wave basin for a specific wave frequency is given in Figure 8. This numerical wave basin was provided with a beach and numerical damping (introduced in [21]) in order to reduce wave reflection. Furthermore, periodic boundaries were provided at the sides, reducing reflection from the side walls. Figure 9 shows that the SPH simulations fit well with the experimental results and that the linear theory with a correction term for the drag force significantly increased accuracy compared to the conventional linear potential flow theory. It is clear that applying linear potential flow theory without the inclusion of the drag force greatly overestimates the RAO, at least for the WEC analyzed in the present study.

5.2.2. Linear Damping PTO System

In the next step, simulations with cylinder2 were carried out with a linear damping PTO system and the average absorbed power was calculated in three different ways: (i) linear potential flow theory, (ii) linear potential flow theory with $C_d$ = 1.5 in Equation (53) and (iii) with DualSPHysics. The results are displayed in Figure 10. Similar to the results was found in [34], the optimal PTO system damping coefficient $B_{PTO,l,opt}$ was significantly larger than $B_{33}$ ($\frac{B_{PTO,l,opt}}{B_{33}}\approx 7$). The average absorbed power was significantly lower when including the drag force, due to the high relative velocities in the resonance region leading to considerable drag forces. It is clear from Figure 10 that the average absorbed power calculated with DualSPHysics lay close to the results obtained from the modified equation of motion with $C_d$ = 1.5. The lower values for average absorbed power obtained with DualSPHysics could be due to the use of only one constant value for $C_d$ in the equation of motion and thus in the calculation of the average absorbed power, while $C_d$ actually increased with decreasing amplitude [32]. At $B_{PTO,l}=0.0$ Ns/m, the WEC’s heave amplitude was approximately $0.10$ m whereas at $B_{PTO,l}=35$ Ns/m the WEC’s heave amplitude was only $0.05$ m.

From Figure 10, the optimal linear PTO damping coefficient $B_{PTO,l,opt}$ can be found. This value was then used to compute the velocity of the WEC in regular waves with $H=0.15$ m, $T=1.2$ s and compare the different modeling techniques—see Figure 11. Figure 11 shows that the velocity calculated with linear potential flow theory with $C_d=1.5$ lay close to the velocity calculated with DualSPHysics, while linear potential flow theory with $C_d=0.0$ overestimates the velocity.

5.2.3. Coulomb Damping PTO System

The same procedure was repeated with a Coulomb damping PTO system: the average absorbed power was calculated with (i) linear potential flow theory, (ii) linear potential flow theory with $C_d$ = 1.5 and (iii) with DualSPHysics. In the case in which linear potential flow theory was applied, Equation (56) was applied for the PTO system force in the frequency domain. For the DualSPHysics calculation, the PTO system force with Coulomb damping, expressed in Equation (55), was implemented in Project Chrono and applied in the DualSPHysics–Project Chrono coupling. The results are displayed in Figure 10. From Figure 10, it is clear that the calculations using Equation (56) gave results close to the DualSPHysics results, meaning that the approximation of the Coulomb damping PTO system force as its first frequency component is valid.

From Figure 10, the optimal Coulomb damping coefficient $B_{PTO,c,opt}$ could be found. This value was then used to compute the velocity of the WEC in regular waves with $H=0.15$ m, $T=1.2$ s and compare the different modeling techniques (see Figure 11). Figure 11 shows that the velocity calculated with linear potential flow theory with $C_d=1.5$ lay close to the velocity calculated with DualSPHysics, while linear potential flow theory with $C_d=0.0$ overestimated the velocity. Figure 11 also shows that the velocity computed with DualSPHysics was no longer a sine wave since the WEC was latched briefly when its velocity reaches zero. This is an inherent property of the Coulomb damping PTO system force $F_{PTO,c}$ which is a square wave and latches the WEC shortly when it reaches its maximum or minimum displacement [48,64]. This behavior was not visible in the linear potential flow theory, with $C_d=0.0$ or $C_d=1.5$, since these theories simplify $F_{PTO,c}$ as a sine wave.

6. Conclusions

In this study, we investigated how the hydrodynamic coefficients of a floating body (added mass and hydrodynamic damping) can be determined with DualSPHysics when using appropriate settings (see Section 4). The artificial viscosity coefficient $\alpha$ should be set to be equal to zero and the heave amplitude of the moving WEC should be kept sufficiently low in order to allow a fair comparison with hydrodynamic coefficients calculated with NEMOH. It was found that the phenomenon of vortex shedding, responsible for causing a significant part of the drag force, is present in DualSPHysics simulations of a cylindrical WEC (see Figure 3). A significant and novel result of the current study is that the drag coefficient $C_d$ of heaving WECs can be determined with DualSPHysics. This is not only useful for WECs but for floating offshore structures in general. For the test cases studied in this research, the interparticle distance $dp$ was found to be required to be lower than $D/50$ with D being the WEC’s diameter in order to obtain accurate results. mDBC should be applied a boundary condition in order to achieve accurate results with reasonable resolution. The results were validated with experimental data obtained from [33,34].

Once an accurate value of $C_d$ is found, the effect of the drag force can be included in the equation of motion of the WEC (Equation (53)) in the frequency domain. In the current study, it was shown that this Equation (53) allowed accurate calculation of the WEC’s heave motion, with results similar to those calculated with DualSPHysics, validated with experimental data from [34]. Furthermore, an updated value of the optimal linear PTO system damping coefficient $B_{PTO,l,opt}$ can be calculated with Equation (54). The average absorbed power taking into account the effect of the drag force was calculated and compared to results from DualSPHysics simulations. The results show that the average absorbed power near resonance is significantly lower when including the drag force and that the $B_{PTO,l,opt}$ obtained from the modified equation of motion taking into account the drag force lies significantly higher than the $B_{PTO,l,opt}$ obtained from the equation of motion without the inclusion of the drag force. These findings are confirmed by simulations of the heaving WECs in DualSPHysics. Future research could further extend the modified equation of motion by applying it to WECs that oscillate with one or more than one degree of freedom, such as pitching or surging WECs.

Besides a linear damping PTO system, a Coulomb damping PTO system has been included in the DualSPHysics–Project Chrono coupling. An approximation of the Coulomb damping PTO system force has also been included in the equation of motion, and the results are compared with DualSPHysics. A major outcome of the current study is that the modified equation of motion including the effect of the drag force gives results similar to those obtained from DualSPHysics for heaving WECs with a linear damping or Coulomb damping PTO system.