**1. Introduction**

Developing efficient, cost competitive and survivable wave energy converters (WECs), has proven, in the last years, to be a challenging task [1–3]. Fluid dynamic analyses are fundamental at every design stage to evaluate the behaviour of the WEC [4], so different methods for simulating the fluid forces and the interaction between waves and semi-submerged bodies have been developed, with increasing levels of fidelity and associated computational cost. Solving the complete set of Navier–Stokes Equations (NSE) for offshore engineering applications is expensive and not advised because the computational cost is too high compared to the accuracy desired [5]. Therefore, the NSE are simplified to obtain a linear equivalent potential flow theory (PF), where the solution is obtained through linearization, thanks to assumptions of inviscid fluid and irrotational, incompressible flow. In this way, one of the most complex set of equations is transformed into a Laplace problem. Under such hypotheses, the Boundary Element Method (BEM) can be implemented, which is a numerical approach for the solution of linear partial differential equations, expressed as integral equations. BEM represents a valid alternative to Reynolds Average Navier–Stokes (RANS) simulation: Being three (or more) orders of magnitude faster [6], BEM is used in many industrial and research applications, especially when computational resources are limited and accuracy is not of primary importance. Linear approaches produce satisfactory results when modelling linear waves and the WEC is moving around the operative point. In fact, non linearities of fluid forces are less important when oscillations

are small and smooth. The success of linear models is justified by their many perks, which include, but not limited to, extreme simplicity, compliance with the superposition principle and great flexibility that makes them compatible with a vast array of mathematical tools. The linearization is always referred to the operating point. In many applications this is a reasonable assumption, since the working condition is not far respect a specific setpoint; in contrast, in the case of wave energy conversion, the objective is to drive the system as far away from equilibrium as possible, in order to maximize the produced energy. This is likely to excite nonlinear dynamics, resulting in non-representative linear models, which usually underestimate energy losses and overestimate the productivity of the WEC [7].

Two main assumptions must be satisfied when using BEMs:


The better these hypotheses are met, the higher the accuracy of the results obtained with the linear model [8]. Although linear PF methods have been used successfully in many offshore engineering applications, the linearizing assumptions are often challenged by realistic WEC behaviour, where large amplitude motions result in viscous drag, flow separation, vortex shedding and other nonlinear hydrodynamic effects [5]. In order to improve the linear time domain model of a WEC, described by the Cummins' equation [9], it is possible to include a viscous damping term as an external force, without significantly increment the computational time required to run the model. In this paper, a notional un-coupled linear and quadratic dependence of the viscous force on the velocity is assumed, hence considering two additive terms. It follows that two viscous drag coefficients should be identified, one for the linear contribution and one for the quadratic contribution. The addition of such non-linear damping force is necessary to avoid divergence in simulations under certain conditions: for example, the model of a forced, undamped system, is likely to explode into very large, non-realistic motion close to resonance conditions. Diverging behaviour can happen without a proper extinction coefficient and is not possible to model some of the conditions the WEC will face at sea [10].

A way to determine damping coefficients is to conduct an experimental campaign, analysing the free decay motion of the WEC [11,12]. For each degree of freedom (DOF), the WEC is released from a set of initial conditions in order to highlight when non-linearities become important; this test is usually performed to identify the roll damping for ships [13]. In order to consider only the significant kinematics of the free decay, the analysed time trace usually covers a few oscillations before the transient elapses; for different starting conditions, the resulting angular acceleration, velocity and rotation angle are post-processed to estimate damping parameters. However, experimental campaigns [14] are expensive and usually not affordable in the early design stage of a WEC. Such a critical obstacle can be overcome by using the use of a CFD Numerical Wave Tank (CNWT). CFD is a Hi-Fi approach to the problem that makes the inclusion of all non-linearities possible. However, CFD is so computationally expensive that it cannot be used for extensive analysis or study of operative conditions: Its most common use is when large non-linearities are expected (survivability, large motion) or for tuning and validating lower-fidelity models, or investigating behaviour of related devices such as tuning tanks [15]. Furthermore, turbulence is a critical feature to handle. There are mainly three techniques to deal with it: Direct Numerical Simulations (DNS), in such simulations, all of the motions contained in the flow are resolved, Large Eddy Simulations (LES) in which NSEs are filtered over the space, finally RANS in which flow quantities are split into the mean value and the perturbation. The resulting equations of LES and RANS are essentially identical, the difference is that in LES most of the turbulence energy is resolved [16]. The URANS approach is chosen because it is the best cost-efficient method considering accuracy and solver time (which could increase easily by several orders of magnitude from RANS to LES and DNS).

The aim of this paper is to define a CFD setup able to emulate a real tank experiment and to describe a methodology to identify the viscous damping coefficients along the two rotational DOFs of a WEC, namely pitch and roll. In particular, the WEC is a prismatic pitching platform that extracts energy from

the pitching motion. On the one hand, viscous losses in pitch are detrimental for power extraction, since they increase the dissipation in the DOF exploited to harvest energy; on the other hand, roll response to oblique or short-crested waves is also detrimental, since it absorbs part of the incoming energy, hence decreasing the energy available to the pitch DOF. Therefore, it is important to obtain a reliable prediction of both pitch and roll responses. The goal of this paper is to increase the fidelity of models based on the potential flow theory by adding the effect of viscous non-linearities. It is worthwhile to remark that most of the literature in the wave energy field is concerned with the identification of viscous losses in the DOF used for energy extraction only (usually a smooth surfaces), while ancillary DOFs are often neglected (with sharp edges in the case study herein considered). Hence, this paper also purports to discuss specific challenges and intrinsic differences between different DOFs, despite using the same test-case and identification technique. Results carried out from simulations are consistent with expectations based on the shape of the hull: the linear term of the damping force prevails in the pitch free decay analysis while the quadratic term, that models the viscous non-linearities, is predominant in the roll DOF.

This paper is organized as follows. After a brief presentation of the device considered in Section 2, the mathematical model of the methodology presented in this paper is discussed in Section 3. The CFD simulation and its setup is explained in Section 4 with emphasis on damping techniques to avoid wave reflections and overset mesh. In the final part of CFD section the mesh convergence study is presented and discussed. Results of simulations are shown in Section 5, firstly as regards CFD scenes, then kinematics and finally the application of methodology discussed in the mathematical section. Ultimately the conclusions of Section 6 involve the analysis of the previously shown outcomes and a last post-processing phase by solving the equation of motion in the time domain, including the non-linear damping effect and the optimization of coefficient identification through the goodness of fit.

#### **2. Device**

The present case considers the WEC shown in Figure 1. Such device, so-called Inertial Sea Wave Energy Converter (ISWEC), is a pitching gyroscope-based WEC, designed for the Mediterranean Sea. ISWEC absorbs wave energy through a floating hull, that drives a gyroscope within contained by means of reacting inertial effects. In particular, waves force the hull into pitching motion, which is transmitted to the gyroscopic system located on a platform inside the floater; the reaction is due to the gyroscopic effects on the spinning flywheel, which induces a torque that is transformed by the power take off (PTO) into electrical energy.

**Figure 1.** Schematic representation of the Inertial Sea Wave Energy Converter (ISWEC) device hull, gyroscope, power take off (PTO) and reference system with respect to the incoming wave.

The hull considered in this work has length *LH* equals to 10m and width *WH* = 5 m.

*J. Mar. Sci. Eng.* **2020**, *8*, 355

#### **3. Mathematical Model and Identification Method**

The time domain linear model of a floating WEC is based on the Cummins' equation [9].

$$(M + A\_{\infty})\ddot{X}(t) + \int\_{0}^{t} h\_{r}(t - \tau)\dot{X}(t)d\tau + KX(t) = F\_{w}(t) \tag{1}$$

where the terms of Equation (1) are:


The convolution term represents, along with *A*∞, the effect of wave radiation in an ideal fluid and it is often referred to as fluid memory effect, because it considers the energy of the radiated waves due to the past motion of the body.

The hydrodynamic parameters in Equation (1) are computed via BEM-based software (e.g., NEMOH [17] and ANSYS Aqwa [18]). It is worth noting that in Equation (1) there are no linearity limitations on the external forces, hence the model is potentially capable to deal with non-linear contributions. Viscous effects are modelled via the addition of a nonlinear term to the right-hand side of Equation (1), divided into linear and quadratic dependence on the pitching velocity ˙ *δ*.

The identification of such damping parameter is conducted using the logarithmic-decrement decay method [11,12,19], which considers the rate of decay of oscillation starting from a non-equilibrium position. During a free decay test, the system oscillates, for each DOF, at its single natural frequency, thus it is possible to determine both the natural period and the hydrodynamic damping. The non-linear equation for free decay pitch motion is:

$$(I\_{\mathbb{FS}} + A(\omega\_\delta))\ddot{\delta}(t) + F\_{\text{vis}, \mathbb{FS}}^{\text{NL}} + K\_{\mathbb{FS}}\delta(t) = 0 \tag{2}$$

where the subscript 55 refers to the pitching DOF, *δ* is the pitch angle and the frequency-dependent added mass *A* is calculated at the natural frequency of the system (*ωδ*). *F*NL vis,55 is the nonlinear viscous force, defined as

$$F\_{\rm vis,55}^{\rm NL} = B\_{\rm 55,1} \delta + B\_{\rm 55,2} \delta |\delta| \tag{3}$$

where *B*55,1 is the coefficient for linear dissipation, and *B*55,2 is the coefficient for quadratic dissipation. Note that *B*55,1 accounts for both linear viscous damping and radiation damping. Equation (3) can be directly included in a non-linear time domain model, with the coefficients *B*55,1 and *B*55,2 identified as the ones that minimize the error between a fidelity-benchmark (CFD, for instance) and the time-domain model. However, it is often crucial to preserve the linear structure of the lower-fidelity model, for example to retain the ability to simulate the model in the frequency domain. Therefore, it is useful to linearize Equation (3), while preserving the ability to discriminate between linear and quadratic contributions. According to [20], it is possible to assume that for each half-cycle, the oscillation is reasonably sinusoidal. Under such an assumption, the non-linear (quadratic) term of Equation (3) can be linearized using the Fourier series expansion:

$$
\delta |\delta| = \frac{8}{3\pi} \omega\_\delta \delta\_\mathbf{i} \delta \tag{4}
$$

where *δ<sup>i</sup>* and *ωδ* are the amplitude and the frequency of oscillation of the *i*-th oscillation cycle, respectively. Consequently, the linearized damping force can be defined as:

$$F\_{\rm vis,55}^{\rm L} = B\_{\rm 55,tot} \delta = \left( B\_{\rm 55,1} + \frac{8}{3\pi} \omega\_\delta \delta\_i B\_{\rm 55,2} \right) \delta \tag{5}$$

*J. Mar. Sci. Eng.* **2020**, *8*, 355

With the inclusion of *F*<sup>L</sup> vis,55 and dividing Equation (2) by the equivalent inertia (*I*<sup>55</sup> + *A*(*ωδ*)), it is possible to obtain the non-dimensional linear equation of motion in the canonical form:

$$
\delta \, \delta + 2\kappa\_{cq}\delta + \omega\_\delta^2 \delta = 0 \tag{6}
$$

where

$$
\alpha\_{eq} = \alpha + \frac{4}{3\pi} \omega\_\delta \delta\_i \beta \tag{7}
$$

with *α* = *<sup>B</sup>*55,1 <sup>2</sup>(*I*55+*A*(*ωδ*)) and *<sup>β</sup>* <sup>=</sup> *<sup>B</sup>*55,2 (*I*55+*A*(*ωδ*)). Equation (6) describes a linear underdamped system, for which the envelope curve of a decay starting from *δ*<sup>0</sup> is defined as:

$$
\delta = \delta\_0 e^{a\_{\text{eq}}t} \tag{8}
$$

Applying Equation (8) for two consecutive peaks *i* and *i* + 1 of the decay curve, it is possible to calculate the logarithmic decay:

$$\frac{\delta\_{i}}{\delta\_{i+1}} = e^{a\_{eq}(t\_{i+1} - t\_i)}\tag{9}$$

Thus, identifying the equivalent linear extinction coefficient:

$$\alpha\_{\epsilon q} = \frac{1}{t\_{i+1} - t\_i} \ln \left( \frac{|\delta\_i|}{|\delta\_{i+1}|} \right) \approx \mathfrak{a} + \frac{4}{3\pi} \omega\_\delta \delta\_{mcm,i} \beta \tag{10}$$

where

$$
\delta\_{\text{mean},i} = \frac{|\delta\_i| + |\delta\_{i+1}|}{2} \tag{11}
$$

and the damped natural pitch frequency can be calculated as:

$$
\omega\_{\delta\_0} = \sqrt{\omega\_{\delta}^2 + \kappa\_{eq}^2} \tag{12}
$$

This method allows to calculate the extinction curve of the *αeq* as a function of the mean amplitude *δmean* for each cycle. To exploit all the available experimental data and obtain a more accurate estimation, it is possible to calculate *αeq* and *δmean* for each oscillation cycle, considering both maxima and minima peaks and grouping the informations of all the decay tests for the specific DOF. The linear regression fit is then performed on the calculated points, with respect to the curve expression:

$$
\alpha\_{t\emptyset} = a\delta\_{mean} + b \tag{13}
$$

which allows to identify the linear and quadratic extinction coefficients:

$$a \equiv b$$

$$\beta \equiv \frac{3\pi}{4\omega\_{\delta}} a$$

#### **4. CFD Model**

The analysis has been carried out on the commercial software Star-CCM+ from Siemens [21] and a Unsteady Reynold Average Navier–Stokes (URANS) approach has been used. Particular attention should be paid to the setup of the CFD simulation environment, since different fluid interfaces, fluid-structure interaction and large solid body motion present numerical problems of challenging resolution. Under such conditions, reliable results can only be achieved by an appropriate and comprehensive model definition. The following subsections detail the most relevant aspects to take under consideration when performing pitch and roll free decay tests of a floating pitching platform as the one considered in this paper.

#### *4.1. Domain, Boundary Conditions and Damping Length*

Due to its oscillation, the floater generates radiated waves that propagate towards the boundaries of the domain, effectively dissipating energy. In order to correctly represent the amount of energy that is subtracted from the floater due to wave generation, it is crucial to capture the wave pattern generated by the hull during the oscillations. In this regard, the domain needs to be large enough to contain the most energetic part of the wave pattern, freely propagating in each direction. Considering *LH* and *WH*, respectively length and width of the hull, domain overall length is 20 *LH* and its width 10 *WH*. Since the wave are generated by the hull of length *LH*, depth of domain is set to 2.5 *LH* in order to respect the deep water approximation Depth/Wave Length > 0.5.

The choice of the dimensions of the domain also depends on wave reflection phenomena, which are a critical problem in both real and numerical tanks. In fact, reflected waves modify the wave field and ultimately affect the dynamics of the floater; Furthermore, it is difficult, if not impossible, to accurately separate the contribution of radiated and reflected waves, leading to improper measures, in case of real experiments, and to numerical noise and simulation crash, in case of numerical tanks. A possible solution consists in adding a damping zone at the boundaries of the domain in order to decrement the energy and the height of waves before they reach the frontier of the domain. Therefore, damping zones hinder the generation of reflected waves, leaving fluid domain perturbed mainly by the radiated waves for the whole duration of the significant decay.

The damping zone should have an appropriate length (at least twice the wavelength [21]) and intensity of the damping [22]. Damping waves is possible by introducing an additional resistance term *Sz* as a source into the momentum equation relative to the vertical motion *w*, according to [22]:

$$S\_z = \rho \left( f\_1 + f\_2 |w| \right) \frac{e^k - 1}{e - 1} w \tag{16}$$

$$k = \left(\frac{\mathbf{x} - \mathbf{x}\_{sd}}{\mathbf{x}\_{cd} - \mathbf{x}\_{sd}}\right)^{n} \tag{17}$$

where:


In this particular case, considering a wavelength (*λ*) comparable to the hull's length, considering the results of Peric et al. [23] and after a trial-and-error evaluation stage, based on estimation of the simulation's residuals, the following values have been imposed:

$$f\_1 = 0; \qquad f\_2 = \frac{10^2 \pi}{\lambda}; \qquad n = 2 \tag{18}$$

Bearing this in mind, the waves reflection phenomenon has been addressed imposing a damping factor *k* in Equation (17) within a distance of 20 m from the boundaries as shown in Figure 2, that represents a top-view of the sea surface of the simulation.

**Figure 2.** VOF Damping Factor *k*.

The boundaries conditions are shown in Figure 3. With respect to the boundaries named *top* and *bottom*, they are set as *velocity inlet*, because a unique phase, air and water, respectively, with zero-velocity is imposed. At *side*, *right* and *left*, the boundary conditions are set as *pressure outlet*, and the pressure is forced to be the atmospheric one for the air section and the hydrostatic one for the water part.

200 400 600 800 1000 1200 1400

Being the domain a function of wavelength, a parametric approach has been used: For the different free decay tests and different initial conditions the wave produced are quite various and so the domain and mesh refinements, to maximize efficiency of the simulation.

#### *4.2. Mesh Generation*

100

200

300

400

500

600

700

In CFD there are two ways to deal with moving bodies: Mesh morphing/remeshing techniques are best when the movement of the body is very small, while they are too computationally expensive and potentially unstable with large displacements of the floater, since the mesh morphing task becomes a significant portion of the overall computational burden. Alternatively, when large motions are expected, overset methods should be preferred [24,25]. Thus the overset approach has been chosen for this study and the domain is split in two different overlapped parts: background and overset. The overset is the one moving with the body, while the background is fixed. With the overset approach,

no remesh operation is needed because the mesh never deforms nor its quality decreases, since it moves with the body and remains unaltered. The overset and background separate regions exchange information due to an overlapping area consisting of some cells, called donor and acceptor, where the solution is interpolated trough Chimera Interpolation techniques according to [26,27]. Equations are solved separately in these different areas and bound together only at the interface between the two zones.

While the dimensions of the background are function of the wavelength, the dimensions of the overset region are constant for all the tests because they are referred to the hull and not to the wave pattern. Another important volumetric refinement is the one around the overset which ensures that, in the zone where information are exchanged, cells belonging to overset and background have the same volume, in order to reduce discretization and interpolation errors. When dealing with trimmed block mesh, it is is important to bear in mind that cells can increase their size only by double.
