Response of Point-Absorbing Wave Energy Conversion System in 50-Years Return Period Extreme Focused Waves

Abstract

This work evaluates the survivability of a point-absorbing wave energy converter at sea states along and inside the 50-year environmental contour for a selected-site in North Sea, by utilizing CFD simulations. Focused wave groups based on NewWave theory are used to model the extreme waves. The numerical breaking waves have been previously predicted by the analytical breaking criterion, showing that the latter provides an accurate estimate for the breaking state. The forces on key components of the device and the system’s dynamics are studied and compared. Slamming loads are identified in the interaction with extreme waves, particularly with breaking waves, and compared with the analytical formulas for slamming estimation as suggested by industrial standards. Considering the extreme wave characteristics, the accompanied phenomena and the resulting WEC’s response, this work contributes to the identification of the design-waves given the environmental contour of the selected site. The top-left side of the contour is identified as the more critical area as it consists of steep and high waves inducing significant nonlinear phenomena, resulting in high loads.

1. Introduction

The offshore renewable energy sector has significant potential to contribute to the decarbonisation of the global energy systems. It encompasses a number of clean energy technologies, including bottom-fixed wind turbines, floating offshore wind and ocean energy technologies. Whereas bottom-fixed offshore wind has reached a commercial viability and is being installed at a larger scale, floating offshore wind, wave and tidal energy technologies are still in the development phase [1]. The challenges for ocean energy technologies include providing reliability, survivability and viability [2].

To ensure that the offshore energy systems are reliable and survive during harsh wind and wave conditions, the response of the systems under extreme conditions need to be modeled and well understood, and design loads must be determined for installation at different sites and for different safety factors. Any offshore structure requires safe and economic design which is closely related with the accurate description of extreme wave–structure interaction. This is the topic of the current paper, with the focus on wave energy converters (WECs) and extreme wave conditions. Due to the small size and highly dynamical motion of WECs, established results from traditional offshore industry cannot be readily translated to the wave energy sector. In addition, since wave energy has not converged around one or a few technology concepts, the design loads identified for one WEC are not necessarily reliable for other wave energy concepts.

A WEC system needs to be designed in order to survive extreme environmental conditions. Uncertainty around the offshore reliability includes the uncertainty in determining the extreme wave conditions. For wave energy systems, the largest loads do not always occur in the highest waves [3,4,5], and peak loads may be due to the motion history of the device and certain combinations of wave height and wave period [6]. International standards [7,8,9] suggest the best practices for determining the site-specific environmental conditions and the design load cases for marine structures, such as offshore floating wind platforms, wave energy converters, mooring systems, and offshore oil platforms.

Statistical methods are implemented and the n-year sea states are identified as the extreme wave conditions in an offshore location, e.g., n = 50 years. The characterization of the n-year return period sea states is a challenging process and entails large statistical uncertainty. Univariate or bivariate extreme value analysis methods are utilized; Peak-over-threshold (POT) is univariate method for the n-year wave definition, while environmental contour is a bivariate approach considering the joint distribution of wave height and period to characterize the environmental conditions. The available univariate and bivariate methods for extreme sea states characterization are reviewed in [10]. High quality data sources of environmental conditions are required for the extreme value analysis. The methodologies and challenges to obtain these data are further discussed in [11].

The inverse first-order reliability method (IFORM) is a common methodology that has been utilized for the environmental contour generation at a specific site [12], the response of the structure along the contour can then be evaluated using numerical and experimental methods. A 50-year return period is recommended for wave energy design [13]. In [14], different methods for obtaining the environmental contours were compared for different sites in the North Sea and Baltic Sea. The IFORM was investigated implementing several joint distributions and the resulting environmental contours were compared. Principal component analysis (PCA) method and POT method were also studied in [14]. The best fit was found using IFORM method utilizing a hybrid distribution, and the resulting 50-year environmental contour for Dowsing site in the UK is used in this paper.

The modeling of extreme wave events is a non-trivial task. The detailed simulation of irregular waves demands high computational cost, while the representation of extreme waves as regular waves is not a realistic approach. A focused wave group is an advantageous approach which has been utilized by [6,15,16,17,18] for the investigation extreme waves interaction with offshore installations.

Nonlinear phenomena accompanying the extreme wave–structure interaction introduce a high uncertainty level. Breaking waves represent a significant challenge to the design of any offshore structure, being a predominant source of nonlinearity, influencing the ability to predict the system’s dynamics. Determining the breaking wave is only half of the problem, the impact load on the structure is the important part. The consequences of breaking waves on WECs are further analysed and discussed in [19]. Research on the breaking wave loading is still very much on-going. The highly nonlinear nature of the phenomenon makes calculation difficult to assess. The effects of breaking waves on monopiles and offshore wind platform applications have been widely investigated, yet limited studies have been conducted for the impact on WECs [16,20,21,22].

In harsh sea conditions, WEC systems can be objected to the so-called slamming leading to significant structural response affecting the ultimate limit state design and the fatigue limit state design of WECs. For milder waves, the loading is relatively well understood using Morison’s equation [23]; however, for the case of extreme waves, the direct use of Morison’s equations will not provide accurate predictions for slamming loads. Numerical tools used in wave energy industry usually lack a function to include slamming loads in the simulations; therefore, the application of slamming loads in the design process is limited. Computational fluid dynamics (CFD) is a possible solution especially when the numerical models have been validated via experimental testings.

As reviewed above, there are still large uncertainties related to the impact of extreme waves on wave energy converters. This work approaches several of the above mentioned critical topics. An open question is about the proper selection of the design-waves for WEC applications. In the present study, the Dowsing site outside the UK is selected, and the extreme sea states are chosen from a 50-year environmental contour, previously identified in [14]. This could be a suitable location for offshore installation in terms of water depth and distance from the shore, but the impact of extreme waves on offshore energy installations at this site has not been investigated before. Another issue is the numerical representation of the extreme waves as a function of reduced computational cost but high accuracy. Here, the focused wave technique based on NewWave theory is utilized for the extreme wave reproduction in the CFD software OpenFOAM. The breaking behaviour of the selected extreme sea states is evaluated based on analytical formulas and CFD simulations. The findings are compared and discussed.

The WEC of consideration in this paper is a point-absorbing WEC based on Uppsala University WEC concept. The WEC consists of a surface buoy connected to a direct-driven linear generator. The limited stroke length of the generator may cause additional extreme loads on the system, which is of particular analogue here. The wave impact on the WEC dynamics and the key components (connection line, generator end-stop spring, and buoy) is studied, and the design choice of the end-stop spring evaluated. The effects of breaking and non-breaking waves on Uppsala University WEC are studied, and their impact on the device is compared. Another critical issue is the wave steepness and here it is evaluated in a realistic approach. The impact of wave steepness on the system is studied as a function of both wave height and wave period. Last but not least, the slamming loads as calculated from the CFD simulations are identified and compared with the analytical methods for the estimation of slamming loads in offshore structures.

The rest of the paper is organized as follows: the selection of extreme wave conditions and the related phenomena are presented in Section 2.1. The NewWave theory for the focused waves generation is briefly described in Section 2.2, and the WEC system is explained in Section 2.3. The set-up of the numerical simulations is presented in Section 2.4. The first results related to the numerical wave propagation in an empty numerical wave tank are presented and compared with the analytical solution in Section 3.1, the WEC dynamics and the forces on critical WEC subsystems are analysed in Section 3.2, and the effects of maximum induced loads are presented further in Section 3.3. In Section 4, the results of the present study are further discussed and suggestions for future work are provided. In Section 5, the main conclusions are drawn.

2. Materials and Methods

2.1. Extreme Wave Conditions

For unmanned offshore renewable energy installations, 50-year extreme wave contours are recommended as suitable design envelopes [7,8,9], which is an approach that has been taken by many authors [24,25,26,27].

In [14], the 50-year environmental contours were generated for the four different sites shown in Figure 1. The present work utilizes the 50-year environmental contour for the Dowsing site, outside the UK, due to the water depth and the distance from the shore. Six sea states along and inside the contour have been selected as shown in Figure 2. The characteristics of each sea state are listed in

Table 1. Characteristics of the sea states 5 to 9, along the Dowsing environmental contour. $H_s$ and $T_p$ are the significant wave height and peak period, respectively. The maximum theoretical crest amplitude of the focused wave, $A_{th}$, is calculated following the New Wave theory. The wave height of the focus wave, H, is defined as the height between the maximum wave crest and the minimum wave trough. $L_p$ is the peak wave length, the parameter $k_p{A}_{th}$ measures the steepness of the waves, where $k_p$ corresponds to the peak frequency, and $KC$ is the Keulegan–Carpenter number.

$\mathit{Sea}\mathit{State}$	${\mathit{H}}_{\mathit{s}}\left[\mathit{m}\right]$	${\mathit{T}}_{\mathit{p}}\left[\mathit{s}\right]$	${\mathit{A}}_{\mathit{th}}\left[\mathit{m}\right]$	$\mathit{H}\left[\mathit{m}\right]$	${\mathit{L}}_{\mathit{p}}\left[\mathit{m}\right]$	${\mathit{k}}_{\mathit{p}}{\mathit{A}}_{\mathit{th}}$	$\mathit{KC}$
5a	5.5	9.0	5.1	9.0	108	0.29	10.2
5b	3.7	9.0	3.4	6.0	108	0.20	6.9
6	6.5	11.5	6.0	10.4	150	0.25	12.0
7	7.1	13.0	6.6	11.5	174	0.24	13.2
8	5.4	14.0	5.0	8.7	190	0.17	10.0
9	3.7	17.5	3.4	6.0	244	0.09	6.9

.

In Table 1, the characteristics of the investigated sea states are summarized. In this study, the focused wave technique based on the NewWave theory is used to reproduce the extreme wave for the given sea state and compute its maximum crest amplitude [28]. The maximum crest amplitude of the focused wave is calculated according to New Wave theory. The wave height H refers to the height between the maximum wave crest and minimum wave trough and is used for the assessment of breaking criteria. The last column of the Table 1 provides information about the Keulegan–Carpenter number, which determines a relation between the inertia and the drag forces acting on the structure and is defined as:

$$KC=\frac{2\pi {\eta }_{max}}{D}$$

(1)

where ${\eta }_{max}$ is the maximum surface elevation at the position of the buoy, obtained in the numerical wave tank with the absence of the buoy, and D is the buoy diameter. For $KC<3$, the inertial force dominates, but, as the KC number increases, the contribution of the drag force increases as well. In the present study, different KC numbers are obtained for the examined sea states, ranging from 6.9 to 13.2 and indicating that the drag term cannot be neglected. These KC values are based on the analytical solution for the wave elevation. As it is shown later, the KC number is affected by the nonlinear phenomena; thus, its value changes.

2.1.1. Breaking Waves

Wave breaking is a natural process involving the conversion of wave energy into turbulent kinetic energy leading to a violent transformation of the free surface and massive hydrodynamic loads are exerted on marine structures [29]. During the wave propagation, once the wave amplitude reaches a critical level, the wave becomes unstable and eventually breaks. In particular, the water particle velocity at the tip of the crest exceeds the wave celerity, hence the upper part of the wave moves faster than the wave itself, so it gets ahead of the main part of the wave and breaks. Breaking waves occur at certain sites, mainly depending on the water depth, the wave length, the wave steepness and the seabed slope. When the water depth is larger than half the wave length ($d/L$ > 1/2), the waves corresponding to this water depth are deep water waves since the wave motion is unaffected by the seabed and their celerity depends on the wave length. As the wave length is larger than the water depth ($d/L$ < 1/20), the waves at this water depth are shallow water waves, the wave motion interferes with the seabed, and their celerity is mainly determined by the water depth. For waves in intermediate water depth (1/20 < $d/L$ < 1/2), the wave motion is partially affected by the seabed and the celerity depends on both wave length and water depth [30].

The wave breaking is a critical phenomenon since it gives rise to energy dissipation and the available energy for harnessing from WEC systems is reduced. There are three basic types of breaking waves: spilling, plunging and surging. Spilling and plunging waves are the most commonly encountered by WECs. The breaker type can be identified based on the Iribarren number [31]:

$$\xi =\frac{s}{\sqrt{H_o/L_o}}$$

(2)

where s is the seabed slope and $H_o$ and $L_o$ are the wave height and wave length in deep water, respectively. Normally, spilling breaker is for $\xi <0.5$, plunging breaker for $0.5<\xi <3.3$ and surging breaker for $\xi >3.3$.

The process of spilling breaking waves, as shown in Figure 3, begins with the formation of a small zone of the forward side of the crest where the water surface becomes rough and the formation of turbulent air-water mixture on the forward face of the wave resulting in whitecap which continuously entrains air riding the forward face of the wave. In plunging waves, a large-scale jet projected out from the forward face of the crest, as a well-defined curl forward with considerable energy which induces entrapment of air between the falling-jet and the forward face of the water, leading to bubbles formation and splash up.

The breaking waves can be identified using breaking criteria. An extensive review of the existing criteria for the inception of breaking waves was conducted by [32]. The breaking criteria can be classified into four types: the McCowan [33] type, the Miche [34] type, the Goda [35] type and the Munk [36] type. In [32], it is stated that Goda’s and Miche’s criteria are advantageous in general. Goda’s formula is suggested when the beach slope is milder than 1/10, whereas Miche’s formula is fairly good for very steep slopes but lacks accuracy for ordinary beach slope.

In the present study, the breaking wave criterion is calculated using the Goda’s criterion [35] for zero slope. Similar studies [37] have utilized Goda’s criterion for the investigation of breaking focused waves group:

$$\frac{H_b}{g{T}_p^2}=A[1-exp(-1.5\pi \frac{h}{L_p})]\frac{tanh\left(kh\right)}{2\pi }$$

(3)

where the $H_b$ is the breaking wave height, $T_p$ and $L_p$ are the peak wave period and peak wave length, k is the wave number and h is the water depth. The parameter A ranges from 0.12 to 0.18 for irregular waves [37]. In Figure 4, Goda’s criterion is depicted (Equation (3)) and the examined sea states are located providing a first estimate regarding their breaking behaviour. Looking at Figure 4, it can be seen that the sea states 5a, 6, and 7 are more prone to break in 22 m water depth compared to the sea states 5b, 8 and 9. However, the sea states 5a, 6 and 7 are included in the marginal breaking region.

2.1.2. Slamming Loads

Slamming loads are caused due to impact pressures which are high in magnitude and short in duration (some milliseconds) compared with pressures exerted by non-breaking waves. Point-absorber WEC type is prone to slamming loads, the buoy is exposed to significant impact pressures when the waves are propelling the buoy upwards and then falling back in the water. This is so-called bottom or vertical or straight slamming and it may occur when buoys with small draft are exposed to energetic wave conditions. The slamming is also result of breaking waves and this is the so-called breaking wave slamming. During the breaking phase, the wave hits the device inducing high loads on the WEC. It is considered that plunging waves are mainly the reason for the impulsive slamming loads [32]. However, Refs. [38,39] consider that slamming events are connected to breaking waves regardless of the breaking type.

The 2D model of Wagner is one of the first theories for the calculation of peak pressures due to slamming effects [40]. Later, Wagner’s theory extended to rigid cylinders and axisymmetric bodies considering the bottom slamming effects [41,42]. Recently, Wagner’s theory has been further extended to slamming loads on vertical and inclined cylinders due to breaking waves [43].

The Morison equation [23] is a good approximation for the quasi-static forces acting on a slender cylinder under the experience of ocean waves as the sum of the inertia force $F_M$ and the drag force $F_D$. The Morison equation for an oscillating cylinder is described as follows:

$$F=F_D+F_M=\frac{1}{2}\rho C_{d}D(u-\dot{X})|u-\dot{X}|+\rho C_m\pi \frac{D^2}{4}(\dot{u}-\ddot{X})$$

(4)

where $\rho$ is the water density, D is the diameter of the cylinder, u is the velocity of the water particle and $\dot{u}$ is the flow acceleration, X is the cylinder displacement, $\dot{X}$ is the cylinder velocity and $\ddot{X}$ is the cylinder acceleration. $C_d$ and $C_m$ are the drag and inertia coefficients, respectively. The inertia force is the sum of the Froude–Krylov forces, $\rho V(\dot{u}-\ddot{X})$, and the hydrodynamic mass force $\rho C_{a}V(\dot{u}-\ddot{X})$. $C_m=1+C_a$, where $C_a$ is the added mass coefficient.

Extreme waves and breaking waves induce very high impact forces on a structure with an extremely short duration. A dynamic component indicating the impact (slamming) force, $F_I$, must be included into Equation (4). The dynamic part is associated with the water impact on the cylinder and the new formulation is:

$$F=F_D+F_M+F_I$$

(5)

Therefore, the wave forces acting on the structure could be divided into the quasi-static part and the dynamic component. The quasi-static part varies in time according to the water propagation. The impact force acts for a time instant and is known as the slamming force. In coastal engineering, the impact force acting on piles due to breaking wave impact is equally distributed along the impact area and is usually calculated according to [43]:

$$F_I=\lambda {\eta }_b{C}_s\rho R{V}^2(1-\frac{V}{R}t)$$

(6)

where the term $\lambda {\eta }_b$ is the height of the impact area as shown in Figure 5. The force changes over the time t. It is observed that at $t=0$ the impact force has its maximum, i.e., at the beginning of the impact. The $C_s$ is the slamming coefficient, according to [44] $C_s=\pi$, whereas, according to [40], $C_s=2\pi$. In [32], several theories are listed suggesting empirical values for $C_s$, but there is still open discussion about the right value.

2.2. NewWave Theory

In this study, NewWave theory [28] is applied for the generation of focused waves. Spectrum wave components adjust their phase so that they focus at a specified time and location achieving energy concentration, and thus maximum wave elevation. The free surface profile is the superposition of N linear wave components, calculated as:

$$\eta (x,t)=\sum _{n=1}^N{\alpha }_{n}cos(k_n(x-x_f)-{\omega }_n(t-t_f))$$

(7)

where ${\alpha }_n$, $k_n$, ${\omega }_n$ are the wave component amplitude, the wave number and the angular frequency, respectively. $x_f$ is the focus position and $t_f$ is the focus time. The wave amplitude of each wave component is calculated as:

$${\alpha }_n=\frac{A_{th}{S}_n\left(\omega \right)\Delta {\omega }_n}{{\sum }_{n=1}^N{S}_n\left(\omega \right)\Delta {\omega }_n}$$

(8)

where $A_{th}$ is the theoretical crest amplitude and $S_n\left(\omega \right)$ is the wave spectrum; in this study, JONSWAP spectrum [45] is used. For the NewWave theory, the amplitude $A_{th}$ of the largest N waves is given by:

$$A_{th}=\sqrt{2m_{o}ln\left(N\right)}$$

(9)

where $m_o$ is the zeroth moment of spectrum. For the calculation of $A_{th}$, $N=1000$, since, in a 3 h sea state, there are approximately 1000 waves according to [46].

The linear horizontal and vertical velocities U and v are given by:

$$u(x,z,t)=\sum _{n=1}^N\frac{{\alpha }_{n}g{k}_n}{{\omega }_n}\frac{cosh{k}_n(z+d)}{cosh\left(k_{n}d\right)}cos(k_n(x-x_f)-{\omega }_n(t-t_f))$$

(10)

$$v(x,z,t)=\sum _{n=1}^N\frac{{\alpha }_{n}g{k}_n}{{\omega }_n}\frac{sinh{k}_n(z+d)}{cosh\left(k_{n}d\right)}sin(k_n(x-x_f)-{\omega }_n(t-t_f))$$

(11)

where z is the vertical coordinate system measured upwards from the mean water level, d is the water depth.

2.3. Uppsala University WEC

Uppsala University’s WEC has a vertical cylindrical buoy that is 3.4 m in diameter and 2.12 m in height. The mass of the buoy is 5736 kg, and the center of gravity is 0.24 m below the still water level. The translator of the PTO system pulls the buoy down due to the translator weight, and the buoy has an initial draft of 1.3 m. The mass of the translator is 6240 kg. The buoy is connected to the PTO through a connection line. Figure 6 shows the PTO system of Uppsala University’s WEC, which is a direct-driven linear generator that has only vertical internal movement.

The forces on the buoy consist of the hydrodynamic forces due to wave–structure interaction, the weight of the buoy, and the force from the connection line. Due to the connection line, the motion of the buoy and the translator are coupled and the forces from the PTO are transferred to the buoy. These forces are the translator weight, the friction damping, and the force from the upper end-stop spring. Once the translator exceeds the upper free stroke length of 1.2 m and hits the upper end-stop spring of the PTO, a force from the spring, $F_{spring}$, acts downward. When the upper spring is fully compressed, the connection wire acts as an elastic mooring cable, which adds further restraint force, $F_{add}$, in the system. The spring stiffness is 776 kN/m and has a compression length of 0.6 m, whereas the connection line stiffness is 300 kN/m. The friction damping force is up to 59 kN, varying according to the generator velocity. Once the translator exceeds the lower free-stroke length, it rests at the lower position of the generator hull, then the connection line slacks and the motion of the buoy becomes uncoupled. The lower end-stop spring force is not considered here.

When the line is not slack, the motion of the buoy and the translator are coupled through the connection line. The restraints from the PTO acting on the buoy are expressed via the restraint force, $F_{line}$, which is directed along a vector from the center of the buoy’s bottom surface to a fixed anchor point, as seen in Figure 7. The total force acting on the buoy is the sum of the hydrodynamic forces, the restraint force, and the weight of the buoy.

In OpenFOAM, the above-described motion is applied by an in-house developed restraint described by the following equations:

$$\begin{array}{c}{m}_b\ddot{r}\left(t\right)=\int \int p\widehat{n}dS+F_{\mathrm{line}}-m_b\overline{g}\end{array}$$

(12)

$$\begin{array}{c}{F}_{\mathrm{line}}={\delta }_{down}(-m_{t}g+F_{\mathrm{damping}}-F_{\mathrm{spring}}-F_{\mathrm{add}})\widehat{r}\end{array}$$

(13)

where $m_b$ is the mass of the buoy, and r is the position vector of the buoy. The $F_{\mathrm{line}}$ should always point downwards, otherwise it is set equal to zero. The friction damping force from the generator is:

$$F_{\mathrm{damping}}=-\frac{{\mu }_{\mathrm{fr}}}{max\left(\right|\dot{r}|,u_{\mathrm{thres}})}\dot{r}$$

(14)

with ${\mu }_{\mathrm{fr}}=59$ kN and $u_{\mathrm{thres}}$ is the threshold velocity. The spring force $F_{spring}$ and the additional elastic mooring force $F_{add}$ are nonzero only when the translator position exceeds the free and total stroke lengths, respectively:

$$F_{\mathrm{spring}}=\left\{\begin{array}{cc}{k}_{\mathrm{spring}}(r-r_{\mathrm{rest}})\hfill & \mathrm{if}(r-r_{\mathrm{rest}})>1.2\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(15)

$$F_{add}=\left\{\begin{array}{cc}{k}_{\mathrm{mooring}}(r-r_{\mathrm{rest}})\hfill & \mathrm{if}(r-r_{\mathrm{rest}})>1.8\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(16)

where $k_{spring}$ and $k_{mooring}$ are the spring and mooring line stiffness, respectively, $r_{rest}$ is the line length when the buoy is resting when the water level is still, the length of the upward/downward free-stroke length is equal to 1.2 m, and the total stroke length is equal to 1.8, as seen in Figure 6. When the translator exceeds the downward stroke length, the line slacks. This is implemented through the ${\delta }_{down}$:

$${\delta }_{\mathrm{down}}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}(r-r_{\mathrm{rest}})>-1.2\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(17)

2.4. Numerical Simulation

2.4.1. Governing Equations

The fluid motion is described by Navier–Stokes equations (NSE) expressing the conservation of mass, momentum and energy, written in a differential form [48]. For Newtonian, incompressible ($\rho =constant$) and isothermal fluid, the NSE are expressed as:

$$\nabla ·\mathbf{u}=0$$

(18)

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla ·\left(\rho \mathbf{u}\mathbf{u}\right)-\nabla ·(\mu \nabla \mathbf{u})-(\nabla \mathbf{u})·\nabla \mu =-\nabla p+F_b$$

(19)

where $\mathbf{u}$ is the three-dimensional velocity field, $\mu$ is the fluid dynamic viscosity, p is the pressure, and $F_b$ is the external body forces. Equation (18) is referred to as the continuity equation and Equation (19) is referred to as the momentum equation.

Numerical solvers implement the Reynolds-averaging approach in order to solve the NSEs. The flow quantities (like velocity $\mathbf{u}$) are approximated as the sum of the time-averaged mean flow ($\overline{\mathbf{u}}$) and the fluctuating part (${\mathbf{u}}^{\prime }$), leading to the Reynolds-averaged-Navier–Stokes (RANS) approach in which a new term; the Reynolds stress is introduced and additional turbulence models are needed to make the equations closed [49].

In the present analysis, the CFD software OpenFOAMv1906 was utilized. The RANS for two incompressible fluids, i.e., air and water, are solved using the cell-centered finite volume method. The free surface is approximated through the volume of fluid method (VOF) as presented by [50]. The phase fraction, $\alpha$, is used to indicate the mixture between air ($\alpha =0$) and water ($\alpha =1$). The local properties ($\rho$ and $\mu$) are then calculated from:

$$\rho =\alpha {\rho }_{water}+(1-\alpha ){\rho }_{air}$$

(20)

$$\mu =\alpha {\mu }_{water}+(1-\alpha ){\mu }_{air}$$

(21)

The conservation of the phase fraction is essential and an additional equation should be added to describe the motion of the phases:

$$\frac{\partial \alpha }{\partial t}+\nabla ·\left(\mathbf{u}\alpha \right)+\nabla ·\left({\mathbf{u}}_c\alpha (1-\alpha )\right)=0$$

(22)

where the last term is an artificial compression term and ${\mathbf{u}}_c=min[c_{\alpha }\left|\mathbf{u}\right|,max\left(\right|\mathbf{u}\left|\right)]$, where $c_{\alpha }=1$. The MULES algorithm is implemented to ensure that $\alpha$ remains between 0 and 1 and the interface compression scheme is used to maintain sharp interfaces during the CFD simulation.

The VOF-RANS is an important feature of numerical solvers, especially for wave energy applications in which wave overtopping and wave breaking phenomena occur. Although supported by the VOF as such, the accurate modelling of breaking waves is a difficult task. Turbulent effects become important to model the energy dissipation in the breaking [51]. According to [52], the $k-\omega$ SST turbulence model has better performance for the prediction of the free surface elevation and is recommended for breaking waves simulations in terms of accuracy and numerical efficiency. The $k-\omega$ SST turbulence model introduces two new transport equations for the turbulent kinetic energy k, and specific dissipation rate $\omega$ [53].

The overInterDyMFoam solver [54] is used for the simulations as it is designed to solve RANS for two incompressible flows keeping track of the surface elevation using the VOF. The pressure–velocity coupling using the PIMPLE (PISO-SIMPLE) algorithm which incorporates the phase fraction equation to derive a new pressure equation.

2.4.2. Computational Domain

In the wave propagation direction (x-axis), the length of the computational domain is defined as a function of the peak wave length; the length from the inlet to the focal position (where the buoy is located and the wave group is expected to focus) is $1\lambda$, and the distance downstream the buoy up to the outlet is equal to $2\lambda$. The computational domain height is twice the water depth for a total of 44 m, and the computational domain width is set at 100 m.

In order to identify the numerical wave tank (NWT) length from the buoy position up to the outlet, a study considering several lengths was conducted. In the empty NWT, the downstream length changes, four lengths were examined: $1\lambda$, $2\lambda$, $3\lambda$, $4\lambda$, and two different sea states; the sea state with the highest steepness (sea state 5a) and the sea state with the longest wave (sea state 9). For the identification of the NWT’s width, a similar study takes place for the empty NWT. Three different widths were examined: 50 m, 100 m and 170 m. For both studies, the surface elevation was measured at the position where the buoy will be placed. For the study of length identification as reference length is considered the case with $4\lambda$ and for the study of width, the reference width is considered to be 170 m. Then, root-mean-square-error (RMSE) analysis follows and the results are presented in

Table 2. Evaluation of NWT’s length downstream the body and NWT’s width. In terms of the length, four lengths have been considered: $1\lambda$, $2\lambda$, $3\lambda$, $4\lambda$. For the width, three cases have been considered: 50 m, 100 m, 170 m. The RMSE value for the surface elevation at the buoy’s position is calculated, having as reference the case for $4\lambda$ and 170 m for the evaluation of the length and width, respectively.

$\mathit{Width}\left[\mathit{m}\right]$	${\mathit{RMSE}}_{5\mathit{a},\mathit{w}}[\%]$	$\mathit{Downstream}\mathit{Length}$	${\mathit{RMSE}}_{5\mathit{a},\mathit{L}}[\%]$	${\mathit{RMSE}}_{9,\mathit{L}}[\%]$
50	0.7	$1\lambda$	2.25	0.25
100	0.5	$2\lambda$	0.41	0.019
-	-	$3\lambda$	0.046	0.0003

. The downstream length selected for the present study is $2\lambda$ as both sea states provides a very small error. The width is selected to be equal to 100 m.

The computational domain consists of four mesh regions: (i) the background mesh, (ii) a refinement region close to the water surface with height $3A_{th}$ and one level higher discretization compared to the background mesh, (iii) an even higher refinement region close to the water surface with height $2A_{th}$ and two levels higher discretization compared to background mesh, (iv) the buoy is surrounded by the overset mesh which has rectangular shape and dimensions $6.5R\times 3R\times 4.7R$. The resolution of the overset region is very close to the spatial discretization of the region (iii). The role of the refinements’ regions close to the water surface is to properly capture the wave propagation. A cell aspect ratio of 1 is used in the working region. Specifically, in the x-direction, the rectangular cells are utilized up to $0.5\lambda$ downstream from the buoy. The refinement regions (ii) and (iii) consist of rectangular cells as well. At the same time, the computational cost should be kept low; therefore, the mesh stretching technique is applied. In particular, in the x-direction, the cells are stretched gradually from $0.5\lambda$ downstream from the buoy to the end of the NWT, and the aspect ratio increases up to 4. In the y-direction, the cells start to stretch gradually at a distance $6R$ far from the buoy until the aspect ratio gets equal to 3. In the z-direction, stretching is applied in the cells of the background mesh, reaching an aspect ratio equal to 2.

The NWT consists of six patches and appropriate boundary conditions are assigned to each one. In

Table 3. Boundary conditions of the computational domain.

Variables	Inlet/Outlet	Seabed	Atmosphere	SideWalls
alpha.water	waveAlpha */zeroGradient	zeroGradient	inletOutlet	zeroGradient
Velocity	waveVelocity *	slip	pressureInletOutletVelocity	slip
Pressure	fixedFluxPressure	zeroGradient	totalPressure	fixedFluxPressure
k	zeroGradient	kqRWallFunction	inletOutlet	slip
$\omega$	zeroGradient	omegaWallfunction	inletOutlet	slip

* IHFOAM toolbox.

, the boundary conditions utilized in the present study are listed. In terms of the $k-\omega$$SST$ turbulence model, the initial values are calculated based on the incoming wave celerity, and $y+$ treatment was utilized. The approximation of the turbulence in the near-wall region wall functions were employed with the first cell-centered placed in the log-law region ensuring the accuracy of the result. According to [55], $y+\in [30,200].$

For the present simulations, the overInterDyMFoam solver was utilized, which has as a base the interFoam solver but is adapted for the overset mesh dynamic mesh method. For extreme wave–structure interaction simulations, the overset mesh dynamic mesh method is recommended [56], since the mesh instabilities due to large body response are eliminated. The motion of the floating body is solved by the standard sixDoFRigidBodyMotion solver. The so-called PIMPLE algorithm is used in the numerical model. The variable time step is adjusted to the Courant–Friedrichs–Lewy (CFL) condition and maximum Courant number is set to 0.5.

2.4.3. Wave Generation

IHFOAM is a newly developed 3D numerical toolbox specially designed for offshore engineering applications, allowing the two-phase solution [57]. Its core is based on OpenFOAM. The specific boundary conditions, applying the static boundary method, allow the generation of any wave type in the 3D domain, from the most simple regular waves to complex irregular directional sea states. One of the most favorable characteristics of IHFOAM is the reduced computational cost compared to other wave generation libraries [58]. IHFOAM has been utilized in similar studies [59,60,61].

The present study considers focused wave groups, and their shape is determined from the JONSWAP wave spectrum using the NewWave theory of [28]. For wave generation, the wave height, wave period and phase of each one of the components are defined at the inlet boundary. For the wave absorption, IHFOAM employs an active wave absorption method in which the fixed outlet boundary cancels the incoming waves by applying the opposite wave velocity based on Shallow Water Absorption. In the present study, long waves are reproduced, and, according to [62], spurious long waves are better absorbed by the active wave absorption, making the IHFOAM toolbox a proper choice for the present work. Following the analytical NewWave theory, the wave groups are designed to focus where the buoy is located.

2.4.4. Number of Wave Components

The number of wave components is critical for the accurate reproduction of focused waves. In this section, a study for the number of wave components determination is conducted. For the evaluation process, the sea state with the highest steepness is selected since it constitutes the most challenging case. The criterion of evaluation is the maximum wave height, $H_{max}$, which is the sum of the maximum crest and the shallowest trough. Six cases were considered with wave components: 60, 80, 100, 120, 156, 210. The CFD solution is compared with the analytic solution according to NewWave theory, $H_{max,anal.}$. In Figure 8, the ratio, $H_{max,anal.}/H_{max}$, as a function of the number of wave components, N, is presented. A large number of wave components for wave spectrum discretization results in more accurate replication of the surface elevation and wave group velocity at the inlet [63]. However, increasing the number of wave components results in higher computational cost. Implementing a tolerance of 2%, 156 wave components are used for the numerical simulations. In order to reduce the computational time at the inlet boundary, linear superposition is used for the surface elevation reconstruction as described in Section 2.2.

2.4.5. Mesh Convergence Study

Grid resolution and convergence study were conducted in order to obtain the final spatial discretization of the computational domain. The study uses the least-squared approach according to [64], and a 2D computational domain was utilized to simulate the steepest wave (sea state 5a). Four resolution cases were examined: Coarse, Medium, Fine, Finest. The cells close to the surface (refinement region iii) were refined by a factor $\sqrt{2}$. The spatial discretization uncertainty value for the maximum wave crest and deepest wave trough are provided in

Table 4. Spatial convergence study for the steepest wave (sea state 5a). The spatial discretizetion is evaluated through the numerical uncertainty U$[\%]$, for the maximum wave crest $U_{MaxCrest}$ and the deepest wave trough, $U_{MaxTrough}$.

$\mathit{Discretization}$  $\mathit{Level}$	$\Delta \mathit{z}$$\left[\mathit{m}\right]$	${\mathit{A}}_{\mathbf{th}}/\Delta \mathit{z}$	${\mathit{U}}_{\mathit{MaxCrest}}$$[\%]$	${\mathit{U}}_{\mathit{MaxTrough}}$$[\%]$
Coarse	0.425	12	2.25	0.82
Medium	0.30	17	0.41	0.73
Fine	0.213	24	0.046	0.61
Finest	0.15	34	0.046	0.56

as well as the discretization levels. The uncertainty values for these quantities are less than 1%, except for the maximum wave crest for the coarse case. Therefore, the ‘Medium’ discretization level is applied for all the simulations in the present study by adjusting the cell size $\Delta z$ in order to keep the ratio $A_{th}/\Delta z=17$.

3. Results

3.1. Extreme Focused Waves Propagation

Figure 9 depicts the numerical propagation of transient focused waves for the six selected sea states and compares with the linear NewWave profile. For the highly nonlinear sea states, the wave crest is higher and narrower, while the wave troughs are shallower. Simultaneously, with the increase of wave nonlinearity, there is a shift in real focal time (sea states 6 and 7). A wave is considered focused at the position where the wave time history is symmetric, i.e., the troughs on either side of the main crest have the same magnitude. The symmetry of waves aside the peak crest is lost (sea state 5a) as the wave nonlinearity increases. The sea states with highly nonlinear behavior are sea states 5a, 6 and 7. According to the breaking criterion (Section 2.1.1), sea states 5a, 6 and 7 are predicted to break, and this condition is confirmed by the CFD simulations. Figure 3 shows the wave propagation for the sea state 6. The vertical dashed black line defines the theoretical focal position where the buoy will be placed later. At $t-t_f=-2s$, the wave breaking process initiates, at $t-t_f=0s$, the evolution of the breaking process continues, and, 4 s later, the focused wave group is completely breaking. Sea state 8 is also characterized by nonlinear phenomena, and the maximum crest amplitude presents a steep and narrow profile and the wave is close to the breaking region. Sea states 5b and 9 are in better agreement with the analytic solution, but nonlinearities still exist. As it has already been stated in [63] and confirmed by this work, the simulation of breaking or near breaking focused waves is a difficult process. The breaking wave criterion, as given in Section 2.1.1 and depicted in Table 4, predicts well the wave breaking condition for the examined sea states since the numerical simulations show that, at sea states 5a, 6 and 7, spilling breaking waves are formulated. However, the breaking criterion considers the maximum wave height ($H_{max,th}$) as it is estimated from the NewWave theory. According to the CFD solution, the real wave steepness ($k_p{A}_{real}$) and real maximum wave height ($H_{max,real}$) are listed in

Table 5. Steepness ($k_{p}A$) and maximum wave height ($H_{max}$) for the focused wave group as calculated by the NewWave theory and the numerical solution (CFD).

$\mathit{Sea}$  $\mathit{State}$	${\mathit{k}}_{\mathit{p}}{\mathit{A}}_{\mathit{real}}$	${\mathit{k}}_{\mathit{p}}{\mathit{A}}_{\mathit{th}}$	${\mathit{H}}_{\mathit{max},\mathit{real}}$	${\mathit{H}}_{\mathit{max},\mathit{th}}$	${\mathit{KC}}_{\mathit{real}}$	${\mathit{KC}}_{\mathit{th}}$
5a	0.36	0.29	9.4	9.0	12.4	10.2
5b	0.22	0.20	6.0	6.0	7.6	6.9
6	0.27	0.25	9.0	10.4	12.8	12.0
7	0.26	0.24	10.3	11.5	14.5	13.2
8	0.17	0.17	7.2	8.7	10.4	10.0
9	0.09	0.09	5.8	6.0	7.2	6.9

and compared with the theoretical values. It can be recalled that the maximum wave height in a focused wave group is the sum of the shallowest wave trough and the preceding wave crest.

In the present study, the focused waves are in the category of intermediate waves, as they fall inside the range $1/20<d/L_p<1/2$. For intermediate waves, the wave velocity is affected both by the wave length and water depth. The existence of the seabed affects the shape of the longer waves (sea states 8 and 9) more as well as their kinematics. In particular, the shallower troughs can be explained also due to the intermediate water depth effects.

Extreme waves are much more non-sinusoidal than the swells, and Fourier analysis is implemented for the identification of the harmonics. Figure 10 shows the numerical wave frequency spectra at focal point obtained by applying Fast Fourier Transform at the time series of Figure 9. The effects of different types of wave nonlinearity is obvious. The linear component is the largest, but, for the breaking or near breaking waves (sea states 5a, 6, 7, 8), the nonlinear components are significant. In particular, 2nd order harmonics are immediately on either side of the main spectral peak, and even a small triple frequency term can easily be recognized. Despite the fact that the 3rd order harmonics are small, according to [63], they are the reason for the shift of focal position and time. The 2nd order harmonics with high frequency are responsible for the increase of the wave steepness. For sea states 5b and 9, the 2nd order harmonics are smaller but still exist and express the contribution of the seabed. In contrast, the lower frequency harmonics correspond to amplitude decrease beneath the focused wave group.

3.2. Response of the Wave Energy Converter

Figure 11, Figure 12 and Figure 13 present the history of WEC’s response at the six sea states. Looking at Figure 11, it can be noticed that the heave response of the buoy follows the focused wave propagation, and the maximum heave displacement is recorded just before the focal time. Low frequency oscillations appear during the heave response time series which are related to the change in surge and pitch direction, i.e., from positive surge/pitch to negative and visa versa. That is to say, once the surge displacement changes direction, the pitch motion changes also direction, and unsteadiness appears in the heave response expressed by small oscillatory motion. In Figure 14, the maximum range of the heave displacement is evaluated using the colormap approach. The maximum range of the heave response is the sum of the minimum heave and the preceding maximum heave displacement as recorded in Figure 11. From the colormap, it is concluded that the maximum range of heave motion appears at sea state 7, which is the sea state with the maximum wave height. In general, the maximum range of heave response is greater for the higher sea states (sea state 5a, 6, 7). Here, it is worth mentioning that the evaluation of maximum heave range is important for point-absorbing WEC technologies since it is a parameter related to the power performance.

Looking at Figure 12, the range of surge motion increases for the sea states with the higher wave height. The maximum range of surge displacement appears close to the focal time, i.e., before the main wave crest, the buoy appears in the maximum left-side horizontal displacement while, after the main wave crest, the maximum right-side horizontal displacement is recorded. The colormap of surge motion in Figure 14 shows the maximum range of surge response for the examined sea states along the 50-year environmental contour. However, when analyzing the surge response during the full simulation time, it can be noticed that, after the focal time, the steep waves tend to push the buoy downstream from its initial position, and the surge motion is shrunk, whereas, in longer waves, the buoy is moving horizontally around its initial position having a greater range of motion. The pitch motion (Figure 13), which is considered as positive in the anticlockwise direction, has a great range for all the examined cases.

In Figure 15, the heave RAO and surge RAO of the WEC are depicted, as a function of the wave steepness. The ${\mathrm{RAO}}_z$ and ${\mathrm{RAO}}_x$, with z and x corresponding to heave and surge respectively, are given from the following formulas:

$${\mathrm{RAO}}_z=\frac{\mathrm{maximum}\mathrm{heave}\mathrm{response}\mathrm{range}}{\mathrm{maximum}\mathrm{wave}\mathrm{height}}$$

(23)

$${\mathrm{RAO}}_x=\frac{\mathrm{maximum}\mathrm{surge}\mathrm{response}\mathrm{range}}{\mathrm{maximum}\mathrm{wave}\mathrm{height}}$$

(24)

where the maximum heave response range is the sum of the maximum heave displacement and the minimum preceding heave displacement of the buoy, the maximum surge response range is the sum of the maximum surge and the minimum preceding surge displacement of the buoy, and the maximum wave height is the sum of the maximum wave crest and the shallowest preceding wave trough. From Figure 15, it can be seen that the optimal steepness ($k_{p}Ath$) is 0.17 in terms of heave, $HeaveRAO$. For greater steepness, the heave, $HeaveRAO$, has lower values and, for sea states steeper than 0.25, the heave, $HeaveRAO$, is dropping significantly. However, for very small steepness ($k_{p}Ath=0.09$), the heave, $HeaveRAO$, is also reduced. The surge, $SurgeRAO$, is reduced almost linearly with the increasing wave steepness showing that the surge motion is prone to long waves rather than high waves. These results should be understood in the context of the examined WEC system: a point-absorbing buoy connected to a linear generator with limited stroke length. When the translator inside the generator hits the upper endstop, it cannot move further, and, by increasing the steepness, the heave motion is not increased further, leading to decreased RAOs for sea states with high waves. However, the excitation force acting on the WEC system due to the increase of the wave height and steepness is increased as well, and this is reflected by the fact that the force on the connection line and the end-stop spring appear to be greater values for the steep and high waves, as will be discussed in detail at the end of this section.

Figure 16 shows the numerical solution for the line force time series. A plateau at around 120 kN appears in all the sea states declaring that most of the simulation time the generator moves within the free-stroke length and the connection line force is equal to the sum of the translator weight and the friction damping force. However, there are some time instances in which several peaks are recorded. Comparing Figure 16 with Figure 17, it can be seen that the force peaks appear when the translator exceeds the upper free stroke-length (red dashed line declares the free stroke length region) and the end-stop spring placed on the top of the generator starts to be compressed. Therefore, the additional force is due to the upper end-stop spring compression, and it is equal to the product of the upper spring stiffness (776 kN/m) and the length of spring compression. Looking at Figure 11, the maximum heave response of the buoy coincides with the peaks in the line force and translator displacement. All the maximum values occur just before the focal time because the buoy interacts earlier with the propagating focused wave group, as the model is displaced away from its initial position. The wave measured at the focal position differs from the impact wave.

As shown in Figure 16, the peak line force occurs at sea states 5a and 7, slightly higher than 400 kN. The additional force ($400-120=280$ kN) corresponds to the additional spring force due to 0.36 m spring compression, as can be verified from the translator displacement in Figure 17. Sea state 5a is the steepest and sea state 7 has the highest wave height; therefore, not only the wave height, but also the steepness is critical for the maximum loads on the connection line. On the other side, sea state 5b and 9 have the same wave height but with sea state 5b being steeper. Comparing sea state 5b and 9, the maximum peak is 245 kN and 280 kN, respectively. The slightly higher peak in sea state 9 accounts for a further instantaneous displacement of the translator resulting in further compression of the upper end-stop spring. However, for sea state 5b, a small peak also appears prior to the main peak, but both last for a few seconds. Comparing the two sea states with the same amplitude but with different steepness, it can be identified that the steepness is mainly related to the number of the peaks, whereas the higher peak force relies on the instantaneous response due to longer wave–structure interaction. Sea state 6 presents also high peaks, but the maximum peak does not appear at $t-t_f=0$ as in the other sea states; here, the maximum appeared approximately 10 s after. The reason is that the upper end-stop spring is compressed more after 10 s, resulting in higher upper end-stop spring force; thus, the total line force is increased.

The breaking waves effects reflect some unexpected peaks at the time series of the connection line force for the sea states 5a, 6 and 7. These peaks appear when the line force follows a decline behavior, a few seconds after the time instant $t-t_f=0$. Trying to explain these peaks, it is figured out that the breaking wave pushes the buoy further away from its initial position, thus the translator continues compressing the upper end-stop spring. At the same time, the buoy experiences high force from the generator pushing it back to the initial position, and the buoy’s velocity is significantly dropping for a few time intervals. During these seconds, the translator stays in the upper position compressing the spring, but the contribution of the friction damping force is negligible, resulting in the previously mentioned small peaks in the connection line force.

Once the lower free stroke length of the translator is exceeded (lower dashed red line in Figure 17), the connection line slacks, and the force is equal to zero. The translator stands on the bottom of the generator, and the motion of the buoy is not affected anymore by the presence of the generator. Looking carefully at Figure 16 and Figure 17, it can be seen that, during sea states 8 and 9, the force in the connection line is zero, but it does not correspond to exceedance of the lower free stroke-length of the translator. The zero value is justified by the fact that the force in the connection line expresses the sum of the forces from the generator (damping force, friction force, translator weight) and, at some time instances, especially when the translator is moving downwards, these forces do not affect the motion of the buoy, and the contribution of the connection line force in buoy’s motion is zero.

Figure 18 provides color depiction of the maximum forces in critical subsystems of Uppsala University WEC, such as the line connecting the buoy with the generator and the upper end-stop spring in the generator. The maximum force on both subsystems is recorded at sea state 5a, 6 and 7 making them the proper choice as design waves. However, an interesting phenomenon is related to the duration of the upper end-stop spring force. The maximum force is recorded in the sea states 5a, 6 and 7, yet for longer waves the compression of the upper end-stop spring lasts longer. For example, in sea state 7, the spring’s compression lasts for 5.5 s, whereas, in sea state 5a, the compression lasts only for 2.5 s. This is also an important factor that needs to be considered at the design stage.

3.3. Slamming Loads in Extreme Wave Conditions

Figure 19 shows the time series of the induced-pressure during the numerical simulations, as recorded from four probes mounted on the buoy. Probes 1 and 3 are located on the frontal side (also denoted impact side), and Probes 2 and 4 are located on the back side (or rear side) of the cylinder. Probes 1 and 2 are mounted 0.24 m above the still water level (SWL) and Probes 3 and 4 are 0.97 below the SWL. The induced-pressure records are presented for the sea states 5a, 5b, 7 and 9. The selection of the sea states is because sea state 5a has the highest wave steepness, sea state 7 has the highest wave height, and the sea states 5b and 9 have equal wave height, but sea state 5b is steeper. The probes mounted above the SWL (except sea state 5b) record higher peaks in the induced-pressure implying the sudden impact of water over the probes. The frontal probes indicate the contribution of the slamming load, which is more intense for the sea states 5a and 7. The induced-pressures the buoy experiences, $\Delta p_{CFD}$, at each sea state, as calculated from CFD simulations, are listed in

Table 6. Induced-pressure on the buoy due to the short impact of the wave–structure interaction as it is given from CFD simulations, $\Delta p_{CFD}$, and, as calculated from the Equation (25), $p_{induced}$. The relative velocity, v, as calculated from the CFD simulations, is also given. The deviation between the theoretical and the numerical pressure is provided.

Sea State	$\Delta {\mathit{p}}_{\mathit{CFD}}$ [kPa]	v [m/s]	${\mathit{p}}_{\mathit{induced}}$ [kPa]	$\mathit{Deviation}$ [%]
5a	80	6.2	99	23.7
5b	20	2.9	26.4	32
7	95	6.6	112	18
9	28	3.5	38.5	37.4

, and they are related to the slamming load. A high pressure area is also detected in the rear part of the cylinder caused by the subsequent back wash of the water around the buoy giving rise to the secondary load cycles. The pressure undulations observed in the rear part of the cylinder are also related to the secondary load cycles, and they are associated with separation and vortex generation at the back side of the cylinder.

According to Det Norske Veritas [7], the space averaged induced (slamming) pressure on a structure can be calculated from:

$$p_{induced}=\frac{1}{2}\rho C_{pa}{v}^2$$

(25)

where $\rho$ is the fluid density, $C_{pa}$ is the space average slamming pressure coefficient and v is the relative normal velocity between the wave surface and the structure. There are two types of slamming and both are observed in the present work; the bottom (straight or vertical) slamming, which is the result of the upward propelling of a structure and the fall back, and the breaking wave slamming, which occurs when a wave hits the structure at a certain moment in its breaking phase [20]. Breaking slamming is considered more critical [20], which is also confirmed in the present study. According to [7], for breaking slamming $C_{pa}=5.15$, and for bottom slamming $C_{pa}=2\pi$.

The wave breaking phenomenon occur at sea states 5a and 7, hence breaking slamming is detected. As can be seen in Figure 20, at sea state 5a, the wave breaks just after the interaction with the buoy, while at sea state 7 the wave breaks just in the frontal head of the buoy. Waves at sea states 5b and 9 do not break (Figure 21); however, high induced-pressure is recorded due to a bottom slamming effect. The theoretical induced-pressure, $p_{induced}$, is calculated according to Equation (25) and compared with the $\Delta p_{CFD}$ in Table 6. The relative velocity v is determined from the CFD simulations, and its values are presented in Table 6. It is obvious that there is a dependency on the impact–pressure on the square of the relative velocity between the hitting water mass and the buoy. The deviation between the $\Delta p_{CFD}$ and $p_{induced}$ ranges from $18\%$ to $37.4\%$. The analytical induced-pressure, $p_{induced}$, is overestimated, which is probably due to the value of selected $C_{pa}$.

Figure 22 shows the inline force acting on the buoy, which expresses the sum of the buoyancy force and the force due to the dynamic pressure integration around the buoy surface. The buoyancy force is approximately equal to 65 kN justifying the fact that the inline force does not start from zero. The contribution of the dynamic pressure varies in time in accordance with water surface elevation associated with the wave cycle. For all sea states, the contribution of the dynamic pressure to the inline force is approximately up to 190 kN. The spikes in the inline force are associated with the buoy interaction with the extreme waves, inducing extra impact forces on the buoy. Apart from the spikes, the inline force time histories are dominated by high frequency oscillations. For small cylinders, the frequency of oscillation is higher compared to larger cylinders [43]. The oscillations take place just after the spikes, and this is the so-called secondary load cycle, and declares the small amount of backward flow acting on the structure once the wave front passes the structure. The secondary load reaches a peak as long as the backward flow continues to act on the structure and, after some time, the flow vanishes and the secondary load cycle finishes. Looking at Figure 19 and Figure 22, it can be concluded that the pressure peaks measured by the frontal probes are associated with the peaks in the inline force. The main peak force, $F_{I,CFD}$, which is indicated by black arrow in Figure 22, occurs when the focused wave front hits the buoy, and its value is given in

Table 7. The slamming force ($F_{I,CFD}$) as calculated from the CFD simulations, the percentage contribution of slamming force ($AdditionalLoad$) in the total inline force, and the factor $\lambda$ identifying the impact area on the buoy.

Sea State	${\mathit{F}}_{\mathit{I},\mathit{CFD}}$ [kN]	$\mathit{Additional}\mathit{Load}$$[\%]$	$\mathit{\lambda }$$[-]$
5a	264	139	0.08
5b	90	47	0.11
7	235	124	0.06
9	115	60.5	0.10

. The additional load, $AdditionalLoad$, acting on the buoy due to slamming is within the range 47%–139%. The contribution of slamming is significant and particularly for breaking sea states; therefore, it cannot be neglected at the design stage.

The impact force is highly dependent on the impact area; the area of cylinder mainly affected by the induced pressure. For the point-absorbing WEC, Equation (6) can be expressed as:

$$F_I=\lambda A{p}_{induced}$$

(26)

where $\lambda$ identifies the impact area on cylinder, A is the total area of cylinder, $p_{induced}$ is the induced/slamming pressure. In order to identify the factor $\lambda$, the Equation (26) is utilized implementing the values for inline force and induced-pressure as calculated from the CFD; $F_I=F_{I,CFD}$ and $p_{induced}=\Delta p_{CFD}$. The value of $\lambda$ is listed in Table 7, showing that a small part of the total area of the cylindrical buoy experiences the impact.

In [65], the effects of slamming on a slender pile were investigated in conjunction with the KC number, and it was concluded that the slamming load increases with increasing KC number; however, slamming was not observed for small KC numbers ($KC<8.26$). In the present work, the KC as defined from the CFD simulations ($K{C}_{real}$) is listed in Table 5 and confirms partially the previous statement. For higher KC numbers, the slamming force is clearly observed; at sea states 5a and 7 ($KC>10$), the slamming force is obvious, while, at sea states 5b and 9 ($KC⩽7.6$), the contribution of slamming force is considerably lower. However, at sea states with high KC numbers ($KC>10$), the rate of slamming variation does not differ significantly, and the additional load due to the slamming force is almost the same for both sea states. Unlike the offshore wind industry applications, the slamming on point-absorbing WEC depends on the instantaneous position of the body, which is highly dependent on the wave steepness. Sea state 7 has higher KC number but lower wave steepness compared to sea state 5a, justifying the fact of the small variation in slamming force.

4. Discussion

The survivability of WECs in extreme waves is a critical issue faced by the developers. The selection of the extreme waves, their numerical reproduction, the response in the design load cases and the identification of highly nonlinear phenomena that accompany the extreme waves are important stages in the design process. The extreme wave conditions for an offshore application are characterized via extreme wave statistics, such as the 50-year sea states. According to International Electrotechnical Commission, the design load cases for WECs in extreme waves need to be defined by utilizing the 50-year sea states for evaluating the structural load responses [8,9].

In the present study, the ability of Uppsala University WEC has been assessed in the 50-year sea states in the Dowsing site, in the UK. The hydrodynamic loads are identified, and the WEC’s response shows that the upper end-stop spring of the linear generator is never fully compressed. The end-stop spring acts as a protection for the hull showing that the design choice for the spring is appropriate. The response in fundamental subsystems, such as the buoy, the connection line, and the upper end-stop spring, is investigated. The extreme waves are numerically reproduced as focused waves based on NewWave theory. The breaking wave impacts on the WEC are investigated based on high-fidelity numerical simulations.

Six sea states are considered, but three of them are considered as critical design load cases. From WECs’ response and the maximum loads on important subsystems, sea states 5a, 6 and 7 can be identified as critical. A drawn conclusion is related to the identification of the crucial environmental contour area. The sea states at the left-top-sided area of the environmental contour results in high load responses, and it is suggested for evaluation during the design process. The waves are steep and high, prone to breaking behaviour and peak loads.

The wave steepness is an important factor for the wave breaking condition; however, only the steepness is not enough for the prediction of the WEC’s response. Comparing two waves with the same height, the steepness determines the loads; in the current work, the comparison of the sea states 5b and 9 showed that the lower steepness case (sea state 9) resulted in higher loads. The wave is much longer than high, pushing the buoy in wider surge motion, inducing higher force in the connection line and the upper end-stop spring. For waves with the same period but different heights, the steeper wave results in higher loads as shown from the comparison between sea states 5a and 5b. However, for waves with different periods and heights, the wave steepness is not enough to predict the critical wave, hence the waves need to be evaluated carefully. As shown in Figure 18, the maximum loads in important WEC subsystems (connection line and upper end-stop spring) are identified at sea states 5a, 6 and 7.

The importance of CFD simulations at the design stage of a WEC system is highlighted in the present study. The identification of peak loads related to extreme waves, like the breaking and bottom slamming load, can be identified via high-fidelity numerical simulations. An important factor for the slamming is the relative velocity between the wave and the device, and the impact area on the buoy, which is dependent on the instantaneous buoy position. Comparing the bottom and breaking slamming loads, it is clear that the breaking slamming load is much more significant for the structure, resulting in 139% additional load on the structure. As the slamming is a highly nonlinear phenomenon with significant implications, the role of CFD simulations is fundamental.

As future work, the following is suggested:

1.

The validation of the CFD model will continue through comparison with experimental tests in which extreme and breaking waves will be reproduced.

2.

Lower-fidelity computational toolboxes play an important role in the design stage due to the reduced computational cost. Therefore, it is essential to verify the capabilities of these toolboxes to capture the critical phenomena enclosed to the nonlinear behavior of the wave–structure interaction, such as the breaking wave behaviour and the resulting slamming effects.

3.

The variation in WEC design could also provide some interesting information about its response in extreme wave conditions. For example, how the buoy geometry contributes and how the generator’s parameters (such as the damping or the end-stop spring) could affect the response.

5. Conclusions

Goda’s breaking criterion was implemented and an accurate prediction for the breaking wave behaviour was achieved as confirmed from the CFD simulations. Three out of the six examined sea states were breaking waves.

The selected sea states were numerically reproduced as focused waves, based on NewWave theory, in an empty NWT. To compare the numerically reproduced waves with the analytical solution, Fourier analysis was conducted, and the higher harmonics contribution was identified particularly for the sea states with breaking waves occurring.

The highest response loads on the WEC were recorded at the sea states on the top-left side of the environmental contour. These are sea states with large wave steepness and wave height with breaking behaviour. Both bottom slamming and breaking slamming were identified, with the latter exerting higher load on the buoy. A small part of the buoy’s area is subjected to the slamming loads, and high pressure acts for milliseconds mainly in the frontal side of the buoy. The induced pressure was computed both using CFD and an analytical method recommended by industry standards, and it was seen that the pressure was overestimated using the latter method.

The force in the connection line is maximized when the upper end-stop spring is compressed. As the purpose of the spring is to prevent the moving parts of the generator from hitting the hull, the fact that the spring was never fully compressed, during any of the extreme sea states, shows that the design capacity of the WEC system was never exceeded.