Hydrodynamic Analysis of a Breakwater-Integrated Heaving-Buoy-Type Wave Energy Converter with an Optimal Artificial Damping Scheme

Abstract

A three-dimensional frequency-domain numerical wave tank (FR-NWT) based on the Rankine panel method was developed. An optimal artificial damping zone (ADZ) scheme was first applied to the FR-NWT to prevent reflection waves from the end walls. Parametric studies of ramp function shape with artificial damping coefficients and damping zone length were conducted to find a proper damping scheme for the frequency domain program. Applying both the Sommerfeld radiation condition and the ADZ scheme to the frequency domain program can reduce the length of the ADZ to less than one wavelength. The FR-NWT developed by the authors was used to calculate the hydrodynamic response of a hemispherical-heaving buoy wave energy converter (WEC) integrated with a seawall-type breakwater of infinite length. A linear power take-off system was used to calculate power generation of the WEC. The global motion of the WEC combined with the breakwater was up to 1.85 times greater than that of the WEC without the breakwater. Moreover, the capture width ratio of the WEC increased approximately 3.67 times more than that of the WEC without the breakwater.

1. Introduction

As the importance of renewable energy has increased due to the Kyoto Protocol and the Paris agreement in United Nations Framework Convention on Climate Change (UNFCCC), interest in the development and utilization of marine renewable energy has also increased [1]. Wave energy is one of the most abundant alternative energy sources [2]. Oscillating water column (OWC)-type wave energy converters (WECs), point-absorber-type WECs, and overtopping-type WECs are common forms of WECs. Among them, the point-absorber type can have the highest efficiency because it instantly converts the motion of a floating body into electrical energy.

Many studies have focused on the power efficiency of WECs [3,4,5]. One way to increase the power efficiency is that WECs could be integrated with other marine structures, such as breakwaters or seawalls [6]. By integrating WECs with the breakwater, the reflected or re-reflected waves from the breakwater can improve the performance of WECs. Additionally, more frequent extreme weather conditions in the future could increase flooding of seawall-type breakwaters, leading to catastrophic coastal flooding. Thus, various studies have been conducted to strengthen the resilience of coastal protection by placing structures such as pillars in front of the breakwater or changing the shape of the breakwater [7,8,9,10,11].

In general, various numerical analyses have been used to evaluate the hydrodynamic performance of WECs, including frequency-domain analysis, time-domain analysis, and spectral-domain analysis [6]. In particular, frequency-domain analysis based on linear potential theory has been used for the initial design of the floating body because the calculation is fast and efficient. The wave-body interaction problem is important for analyzing the dynamic response of floating structures, such as WECs.

Many studies on wave-body interaction problems have been conducted since the 1980s. The wave-body problem can usually be solved by numerical analysis using the free-surface wave Green function or Rankine source panel method, assuming potential flow [12]. The wave-body problem can be divided into a diffraction problem and a radiation problem. Solving the radiation problem requires a proper radiation condition to make an open sea condition. One of the well-known radiation boundary conditions is the Sommerfeld condition [13]. This boundary condition is easy to apply and does not require additional free surface meshes. Although this boundary condition has the advantage of not requiring an additional free surface mesh, it is ineffective under nonlinear wave conditions [14]. The reason is that the Sommerfeld radiation condition is a method of forcibly radiating from the end boundary using the information of propagated incident waves. Therefore, there is a limit to radiating nonlinear waves or scattered waves by objects at the end boundary.

In the Rankine panel approach used for time-domain numerical wave tank (NWT) development, it is important to establish far-field radiation conditions. In addition, because the time domain NWT has the same conditions as the physical experimental wave tank, it is essential to solve the reflection problem of the incident wave to satisfy the open sea condition [14]. Therefore, an artificial damping zone scheme was applied to satisfy the open sea conditions in the time-domain NWT. In this scheme, by adding artificial damping terms to the free surface boundary condition, all reflected waves except the incident wave were forcibly attenuated in the incident wave boundary. All waves are forcibly attenuated at the end-wall boundary. On the other hand, the artificial damping scheme requires additional free-surface panels because the length of the artificial damping zone requires at least one wavelength to properly attenuate reflected waves [14,15].

Many studies on the artificial damping system focusing on minimizing the damping length and maximizing damping performance have been reported. Several studies have examined the wave radiation boundary conditions in the time domain. The wave radiation problem was first solved by applying an artificial damping scheme [16]. In addition, a study was conducted by applying two damping terms ($\varphi$, $\eta$) to kinematic and dynamic free surface boundary conditions [17]. A dual damping system that combined the artificial damping zone with a piston-type wave absorber at the end boundary was proposed in [18].

An artificial damping scheme applying two damping terms ($\varphi$, $\eta$) only to the kinematic free surface boundary condition was developed [19]. In the two-dimensional fully nonlinear NWT development, two damping terms (${\varphi }_n$, ${\eta }_n$) were used for kinematic and dynamic free surface boundary conditions, respectively [20]. An optimal ramp function and length of the artificial damping zone applicable to a three-dimensional NWT were proposed [9]. Four different shapes of ramp function applied to the damping coefficient were introduced to evaluate the damping performance of the artificial damping zone. In addition, it was confirmed that the total wave energy in the computational fluid domain was preserved despite the propagation of incident wave propagation.

A three-dimensional, time-domain, fully nonlinear NWT was equipped with frontal, side, and end artificial damping zones [21,22]. In addition, a three-dimensional NWT combining two radiation boundary conditions (Sommerfeld condition and artificial damping zone (ADZ)) was developed [23,24]. The application of the combined conditions was very effective in removing waves from the sidewall and end boundaries.

Studies using the ADZ technique based on the Rankine panel method are used mainly for time-domain analysis. On the other hand, a two-dimensional frequency domain program with an artificial damping zone was developed [25]. This was the numerical simulation for the wave radiation problem of a floating body in two-layer fluid.

In this study, a three-dimensional frequency-domain numerical wave tank (FR-NWT) based on the Rankine panel method was developed to analyze the hydrodynamic performance of a heaving-buoy-type WEC integrated with a breakwater. Since this frequency-domain analysis model can solve the radiation and diffraction problems in the rectangular fluid domain, various wave–ocean structural interactions can be calculated.

The ADZ scheme, which is mainly used for the time-domain analysis, was first applied to the 3D frequency-domain, FR-NWT. In addition, to apply the ADZ of the optimal condition, the optimal length of ADZ, type of damping terms, and shape function of ramping function were estimated. The far-field radiation performance was compared by applying a combined radiation boundary condition (both ADZ and Sommerfeld radiation conditions).

Based on the developed FR-NWT, the hydrodynamic performance of a hemispherical-buoy-type WEC combined with a seawall-type breakwater of infinite length was evaluated. The interaction between the floating WEC and the breakwater was used to amplify the global motion and generate power under specific conditions.

The rest of this paper is organized as follows. Section 2 presents numerical schemes for developing the FR-NWT. Section 3 can be divided into two parts. The first part focuses on investigation of the radiation boundary condition suitable for application to the FR-NWT. The second part focuses on the motion amplification and power generation efficiency of the WEC integrated with the breakwater using the developed FR-NWT. Finally, Section 4 summarizes the results and future challenges.

2. Numerical Model Development

2.1. Boundary Value Problem

In the computational domain, a three-dimensional FR-NWT assumes that the fluid is potential flow, which is inviscid, irrotational, and incompressible. In this case, the governing equation is the Laplace equation (Equation (1)) using the velocity potential and continuity equation.

$${\nabla }^2\varphi =0$$

(1)

The Laplace equation can be transformed to the boundary-integral equation using the Green second identity. The boundary-integral equation can be written as follows:

$$\alpha {\varphi }_i={{\displaystyle \iint }}_{\Omega }\left(G_{ij}\frac{\partial {\varphi }_j}{\partial n}-{\varphi }_j\frac{\partial G_{ij}}{\partial n}\right)ds$$

(2)

where $\varphi$, $\alpha$, and $G_{ij}$ are the total velocity potential, solid angle, and Green function, respectively. The solid angle is 0.5 on the boundaries of the computational domain. Equation (3) shows the simple Rankine source as the Green function used in this study, and the image method is used to express the flat seabed and reduce the computational demands (Equation (4)).

$$G=\frac{1}{4\pi }\left(\frac{1}{R_1}+\frac{1}{R_2}\right)$$

(3)

$$\begin{gathered} R_1=\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2+{\left(z-z_0\right)}^2} \\ {R}_2=\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2+{\left(z+z_0+2h\right)}^2} \end{gathered}$$

(4)

where ($x$, $y$, $z$) and ($x_0$, $y_0$, $z_0$) are the coordinates of the field and source point, respectively; and $h$ is the water depth. The total velocity potential can be divided into the incident wave potential (${\varphi }_I$), radiation potential (${\varphi }_R$), and diffraction potential (${\varphi }_D$), as expressed in Equation (5). Incident wave potential can be calculated using Equation (6). Additionally, the radiation potential was decomposed conventionally as six degrees of freedom-of-motion mode components.

$$\varphi ={\varphi }_I+{\varphi }_R+{\varphi }_D$$

(5)

$${\varphi }_I=-i\frac{gH}{2\omega }\frac{\mathit{cosh}k\left(z+h\right)}{\mathit{cosh}kh}·e^{ikx}$$

(6)

$${\varphi }_R=i\omega {\displaystyle \sum }_{j=1}^6{\zeta }_j{\varphi }_j$$

(7)

where $\omega$, $k$, and $H$ are the wave frequency, wave number in the finite water depth, and wave height, respectively, and ${\zeta }_j$ denotes the complex amplitudes of the body motion in its six degrees of freedom. To obtain radiated wave potential and diffracted wave potential, the radiation problem and the diffraction problem must be solved, respectively. The radiation problem aims to analyze oscillating bodies in calm water. The diffraction problem evaluates the wave excitation forces acting on the stationary body in waves. Figure 1 presents an overview of the computational domain of 3D FR-NWT, where $S_B$, $S_{fs}$, $S_{Bt}$, and $S_w$ are the body boundary surface, free-surface boundary surface, bottom boundary surface, and sidewall boundary surface, respectively.

2.2. Boundary Condition

The body boundary condition should be applied differently depending on the hydrodynamic problem. In the radiation problem, the body boundary condition was described as the velocity of the rigid body ($V_j$) times the normal vector ($n_j$) in Equation (8).

$$\frac{\partial {\varphi }_j}{\partial n}=V_j·n_{j}on{S}_B \mathrm{for}\mathrm{the}\mathrm{radiation}\mathrm{problem}\mathrm{only}$$

(8)

where $j$ is modes of rigid body motion, e.g., $i$ = 1–3 for translational motion and $i$ = 4–6 for the rotational motion. Equation (9) expresses the body boundary condition for the diffraction problem.

$$\frac{\partial {\varphi }_D}{\partial n}=-\frac{\partial {\varphi }_I}{\partial n}on{S}_B\mathrm{for}\mathrm{the}\mathrm{diffraction}\mathrm{problem}\mathrm{only}$$

(9)

The linearized free-surface boundary condition, which is a combination form of the kinematic and dynamic free-surface boundary conditions, is expressed as Equation (10) [25].

$$\frac{\partial \varphi }{\partial z}-\frac{{\omega }^2}{\mathrm{g}}\varphi +C=0$$

(10)

where $g$ and $C$ denote the gravitational acceleration ($=9.81 \mathrm{m}/{\mathrm{s}}^2)$ and an artificial damping term if ADZ is applied, respectively. Section 3.2 provides details of the damping terms. There are three options for the side- and end-wall boundary condition in the frequency domain: the Sommerfeld radiation boundary condition, an impermeable boundary condition, and the wave-adapted wall condition. Conventionally, FR-NWT uses the Sommerfeld boundary condition (Equation (11)) for the open sea condition, as reported elsewhere [22]. $K$ is the wave number in infinite water depth, and R denotes the distance from any fixed point. In the case of adopting the ADZ method, the wall is expressed by the impermeable condition ($\partial \varphi /\partial n=0$, second option). The last is the same as the body boundary condition for the diffraction problem in Equation (12). This condition is used as the boundary condition for a fixed breakwater in the frequency-domain diffraction problem. Moreover, the image method is used to express the flat bottom:

$$\underset{R\to \infty }{\lim }\sqrt{R}\left( \frac{\partial \varphi }{\partial R}+ik\varphi \right)=0 on S_w$$

(11)

$$\frac{\partial {\varphi }_D}{\partial n}=-\frac{\partial {\varphi }_I}{\partial n}on{S}_w$$

(12)

Finally, the boundary-integral equations for the radiation and diffraction problems can be solved with those given boundary conditions. The radiated wave potential and diffracted wave potential can be obtained from those hydrodynamic problems. The HG matrix (influence matrix) was evaluated using the Hess and Smith [26] method and the multipole expansion method [1]. This study used the LU decomposition method as a matrix solver, as well as the back-substitution method.

2.3. Equation of Motion

Each hydrodynamic coefficient and wave excitation force/moment can be calculated by solving the boundary value problem. The radiation velocity potential was used to obtain the added mass or added moment of inertia (Equation (13)), as well as the radiation damping coefficient (Equation (14)).

$$A_{ij}=Re\{\rho \underset{S_B}{\overset{}{{\displaystyle \iint }}}{\varphi }_i{n}_j dS\}$$

(13)

$$B_{ij}=Im\{-\rho \omega \underset{S_B}{\overset{}{{\displaystyle \iint }}}{\varphi }_i{n}_j dS\}$$

(14)

where $A_{ij}$ and $B_{ij}$ denote the added mass or added moment of inertia and radiation damping coefficients, respectively. Furthermore, the wave excitation forces were evaluated based on the Bernoulli equation with the incident wave potential and the diffracted wave potential as in Equation (15).

$$F_{e{x}_j}=\underset{S_B}{\overset{}{{\displaystyle \iint }}}\left(-\rho \frac{\partial \left({\varphi }_I+{\varphi }_D\right)}{\partial t} \right)n_j dS=i\rho \omega \underset{S_B}{\overset{}{{\displaystyle \iint }}}\left({\varphi }_I+{\varphi }_D\right)n_j dS$$

(15)

The equation of motion, including the external force and the PTO force, was set according to Equation (16).

$$\left(m+A_{ij}\right)\ddot{{\zeta }_j}+B_{ij}\dot{{\zeta }_j}+C_{ij}{\zeta }_j=F_{e{x}_j}+F_{pto}$$

(16)

where $C_{ij}$ represents the restoring coefficients, which can be easily obtained by a theoretical solution. $F_{pto}$ is the additional applied force from the power take-off (PTO) system. The effect of the PTO system is usually presented as a linear damping system (Equation (17)). Equation (18) represents the mean power generated by the floating body velocity according to the incident wave frequency. In general, the linear generator or the hydraulic system was selected as the PTO system for the movable body-type WEC system [1,27,28,29]. On the other hand, the PTO force was not easy to consider in the frequency domain model because of its strong nonlinearity. In this calculation, an equivalent linear model was applied. In addition, the radiated wave elevation from the oscillating body can be evaluated as Equation (19).

$$F_{pto}=b_{pto}\dot{\zeta }$$

(17)

$$\overline{P_g}=\frac{1}{2}{b}_{pto}{\dot{\zeta }}^2$$

(18)

$${\eta }_{{\mathrm{R}}_{\mathrm{j}}}=-\frac{1}{\mathrm{g}}{\left(\frac{\partial {\varphi }_j}{\partial t}\right)}_{z=0}$$

(19)

where $b_{pto}$ is the additional damping coefficient from the PTO system.

3. Numerical Results and Discussion

3.1. Numerical Model and Convergence Test

Hydrodynamic analysis for a hemispherical heaving-buoy-type WEC system was performed based on the aforementioned FR-NWT, a Fortran-based in-house program for hydrodynamic analysis in the frequency domain. The WEC system employs a hemispheric heaving buoy as an actuator.

Table 1. Calculation conditions for Sommerfeld radiation condition.

Characteristics	Value	Unit
Wave frequency ($\omega$)	1.5–3.5	rad/s
Wave height (H)	1.0	m
Radius of a floating body (R)	1.0	m
Draft of a floating body (d)	1.0	m
Water depth (h)	20.0	m

lists the conditions for a floating body and wave conditions. Prior to performing the analysis of the WEC system, the shape of the body was fixed, and the unit wave height was used to verify the radiation boundary conditions to be applied to the FR-NWT. In particular, the water depth ($h$) was set to 20 m to neglect the seafloor effect in all applied wave frequency ranges. Figure 2 presents the mesh system of 3D-FR-NWT, which has a rectangular shape, to analyze the influence of a flat wall on a floating body motion. A floating body was located in the middle of the free surface. The size of the computational fluid domain was 3λ (wavelengths) × 3λ. The sensitivity of the number of free-surface nodes was examined by conducting convergence tests on the number of nodes per wavelength by applying from a sparse mesh system to a finer mesh system with the Sommerfeld radiation condition. Figure 3 compares the heave added mass and heave radiation damping coefficient on the number of collocation nodes per wavelength at $\omega$ = 3.2 rad/s. The results were compared with a frequency-domain hydrodynamic analysis tool, WAMIT [30]. The results were based on the wave Green function, which satisfied the radiation condition very well. When the number of nodes per wavelength was more than 25, the hydrodynamic coefficient converged, providing the same results as those reported with WAMIT. Twenty-five collocation nodes per free surface were used in the free surface for the remaining numerical simulations.

3.2. 3D-FR-NWT with Radiation Boundary Conditions

Based on the optimal mesh system, in this study, we applied three types of radiation boundary conditions (1: Sommerfeld radiation condition (SRC), 2: artificial damping zone (ADZ), 3: combined radiation condition between the SRC and ADZ) and compared their performance.

First, the SRC was applied only at the sidewall boundary surface. Figure 4 compares the hydrodynamic coefficients (added mass and radiation damping coefficient) and wave excitation vertical forces on the hemispheric heaving buoy. The results were in agreement with the WAMIT results, except for under low-frequency conditions. This is because low-frequency radiated waves (long waves) are reflected at the end boundary of the computational domain (see Figure 5). Hence, in the case of the rectangular calculation domain, if only the Sommerfeld radiation condition is applied, the low-frequency wave near the apex cannot be properly radiated on the end boundary. Therefore, when numerical modeling is performed in a rectangular wave tank, an artificial damping zone must be adopted to solve the wave radiation problem.

An artificial damping scheme was used on the free surface to prevent the reflected or re-reflected waves in the fluid domain. Figure 6 presents a plan view image of the free surface with artificial damping zones installed at all boundaries of 3D-FR-NWT. Two artificial damping terms were added to the linearized free surface boundary condition (see

Table 2. Free surface boundary conditions with two different types of artificial damping terms.

Damping Type	Free Surface Boundary Condition
$\eta -\varphi$ type	$\frac{{\partial }^2\varphi }{\partial t^2}+g\frac{\partial \varphi }{\partial z}+{\mu }_1\frac{\partial \varphi }{\partial t}+{\mu }_2\varphi =0$
${\varphi }_n-\eta$ type	$\frac{{\partial }^2\varphi }{\partial t^2}+g\frac{\partial \varphi }{\partial z}+{\mu }_1\frac{\partial \varphi }{\partial t}+{\mu }_2\frac{\partial \varphi }{\partial t}=0$

). The efficiency of the artificial damping technique can be increased by applying two damping terms instead of one to the free surface boundary condition [14]. The $\eta -\varphi$-type damping method substitutes the damping terms proportional to the wave elevation and velocity potential to the linearized free surface boundary condition. The ${\varphi }_n-\eta$-type uses ${\varphi }_n-$ and $\eta -$ damping terms for the free surface boundary condition based on a previous study [13]. The relationship between the two coefficients can be expressed as ${\mu }_1^2=4{\mu }_2$, which means that the frequency of the incident wave does not change when it enters the artificial damping zone.

Based on the optimal case of [14], Figure 7 compares three shape functions of the artificial damping coefficient for various angles ($\theta$). All shapes have a small initial slope when entering the damping zone so that the occurrence of the reflected waves at the initial point of ADZ was reduced. In the case of ‘shape 1’, the slope of damping coefficients increased rapidly compared to the cases of ‘shape 2’ and ‘shape 3’. Figure 8 shows the radiated wave amplitude to which the three ramp shapes were applied to compare the wave damping effect. The wave frequency was 1.5 rad/s. In the results using ‘shape 1’ and ‘shape 2’, some fluctuations were observed at the starting point of ADZ. This is the effect of the reflected wave at the front of ADZ because the initial slope of the shape function was relatively large. In the case of ‘shape 3′, however, fluctuation of the wave amplitude did not occur. This means that the waves were attenuated properly at the start point of the artificial damping zone ($x/\lambda =1.0$). Therefore, in this study, we adopted the ‘shape 3’ ramp function for stable numerical analysis. The ramp function to which the damping coefficient was applied is presented for each artificial damping zone, as expressed in Equations (20)–(22).

$$u\left(x,y\right)=\frac{{\mu }_0}{1-\cos \theta }[1-\cos (\frac{\left|y-L_{fsy}\right|}{L_{damp}}\mathsf{\theta })]on{S}_{A1}$$

(20)

$$u\left(x,y\right)=\frac{{\mu }_0}{1-\cos \theta }[1-\cos (\frac{\sqrt{{\left(x-L_{fsx}\right)}^2+{(y-L_{fsy})}^2}}{\frac{L_{damp}}{\theta }}\mathsf{\theta })]on{S}_{A2}$$

(21)

$$u\left(x,y\right)=\frac{{\mu }_0}{1-\cos \theta }[1-\cos (\frac{\left|x-L_{fsx}\right|}{L_{damp}} \theta )]on{S}_{A3}$$

(22)

where $\left(x,y\right)$ and $\left(L_{fsx},L_{fsy}\right)$ denote the x and y coordinates of collocation nodes in ADZ and the starting points of ADZ, respectively. $L_{damp}$ is the length of ADZ and ${\mu }_0$ is the target value of the damping coefficient. The optimal angle ($\theta$) is $\pi /2$.

Figure 9 compares the results of applying two different types of damping terms. Shape 3 ramp function and one wavelength-damping zone were applied. In the time-domain result of [14], both the $\eta -\varphi$-type and ${\varphi }_n-\eta$-type damping schemes exhibited stable damping performance when the damping zone was 1.5 times the wavelength. On the other hand, in FR-NWT, the damping effect was sufficient, even when the length of the damping zone was approximately one wavelength, and there was little difference according to the damping term type. Briefly, the ${\varphi }_n-\eta$-type showed better damping performance in time-domain analysis, but there was no significant difference in the frequency-domain problem. The ${\varphi }_n-\eta$-type was judged to be slightly better, but there was little difference between the two damping types.

Figure 10 presents the heave-induced radiated wave amplitude with an artificial damping zone. The radiation problem was performed by forcibly oscillating a heaving buoy. No nonphysical waves around corners of the computational fluid domain were observed in either frequency case, unlike Figure 5. Therefore, the artificial damping zone works well in the high-frequency and low-frequency regions. Therefore, FR-NWT adopted the ${\varphi }_n-\eta$-type damping terms with the shape 3 ramp function.

A parametric study of the damping coefficients was conducted. For this analysis, the radiation problem was performed with various damping coefficients in $\omega$ $=1.5 \mathrm{and} 3.0 \mathrm{rad}/\mathrm{s}$, as shown in Figure 11. The results were compared with those of WAMIT. When the artificial damping coefficient is more than 2.0, the heave-added mass and damping coefficient converged and agreed well with the WAMIT results. Another parametric study on the length of the ADZ was performed (Figure 12).

The combined radiation boundary condition (SRC and ADZ) was employed to reduce the length of the ADZ, and the corresponding results were compared with those obtained when only the ADZ was applied. When the length of the ADZ was longer than one wavelength, the results of FR-NWT showed good agreement with those of WAMIT, as suggested in previous studies [22,31,32].

The heave motion RAOs were compared in Figure 13a by applying the radiation boundary conditions, and their corresponding computational demands were investigated (Figure 13b). The combined boundary condition (both ADZ and SRC) was applied to reduce the length of the damping zone, as described in [23,24]. When the result using the ADZ alone ($L_{damp}=0.5\lambda )$ and the results using the combined condition ($L_{damp}=0.5\lambda )$ were compared, the result using the combined radiation condition was more consistent with WAMIT. As a result, the length of ADZ could be reduced to a half wavelength by applying the combined radiation condition.

The heave RAO under the combined radiation condition agreed well with the results of WAMIT, in contrast to the results obtained with only the SRC applied. Figure 13b shows the calculation time for the ADZ only, as well as the combined condition (ADZ and SRC). The computational time was reduced by decreasing the number of meshes in the domain because of the combined boundary condition. Reducing the calculation time by reducing the total number of meshes does not affect the accuracy of the calculation results. It is important to reduce the computational demand to perform the analysis over a wide frequency range. For this, the combined radiation condition was applied to the FR-NWT. Future calculations will be performed using the FR-NWT equipped with the above optimized radiation boundary conditions to simulate a heaving buoy-type WEC-integrated breakwater.

3.3. Wave Energy Converter Integrated with Breakwater

The dynamic response of the WEC in front of the wall-type breakwater was analyzed using the FR-NWT, to which the ADZ and SRC were applied together. Figure 14 shows an overview of the computational fluid domain, and Equations (8) and (12) were used for the radiation and diffraction boundary conditions of the breakwater, respectively. A single WEC integrated with a seawall-type breakwater of infinite length was considered.

In addition, no artificial damping zone was installed on the free surface near the breakwater to examine the effect of wave reflection. To verify this numerical model, the heave-motion RAO of a cylindrical floating WEC in the presence of a breakwater was estimated and compared with the results reported elsewhere [33] (Figure 15). The WEC had a diameter of ($D$) = 2.0 m and draft of ($d$) = 1.0 m. The distance between the breakwater and WEC was set to five times the floater’s diameter ($D$). The calculation results were agreed upon well in both cases with or without a breakwater. The interaction between the WEC and the breakwater amplified the heave RAOs at the frequency region lower than the natural frequency, although the amplitude of the RAOs was reduced at $\omega =$2.05–2.2 rad/s. Therefore, the distance between a floater and breakwater should be carefully selected to enhance the power take-off performance.

The hydrodynamic performance of a hemispherical heaving-buoy-type WEC integrated with a breakwater was investigated.

Table 3. Analysis conditions of a WEC integrated with a breakwater.

	Characteristics	Value	Unit
Body (Hemisphere)	Diameter (D)	2.0	m
	Draft (d)	1.0	m
	Mass (m)	2084	kg
Wave condition	Incident wave frequency （$\omega$)	1.2–4.0 $\left(\mathsf{\Delta }\omega =0.05\right)$	rad/s
	Wave height (H)	1.0	m
Wall	Nondimensional distance between wall and body (α/λ)	0.2–0.5	-

lists the calculation conditions of the WEC system. The distance between the body and the wall (α) was established by fixing the diffraction parameter ($D/\lambda$) and comparing the heave RAOs according to $\alpha /\lambda$ (Figure 16).

As $D/\lambda$ was smaller than 0.2, that is, the body is slender, $\alpha /\lambda$ with the minimum heave RAO of the WEC was 0.25, and $\alpha /\lambda$ with the maximum heave RAO was 0.5. It seems that standing waves are caused by waves reflected from the breakwater. Therefore, the motion amplification is greater when the WEC lies at half the incident wavelength. On the other hand, as $D/\lambda$ increases, $\alpha /\lambda$ with the maximum heave RAO of WEC decreases. Therefore, $\alpha$ must be determined according to the wavelength. If α becomes too large or small, it is too difficult to integrate with the breakwater; hence, only two cases were analyzed (case 1–2). Thus, $\alpha$ was fixed to 2.41 m in case 1 and 3.85 m in case 2.

Assuming that a hemispherical WEC has a radius of 1 m, the heave RAO of the WEC (see Figure 17) was obtained. Motion amplification occurred up to 1.87 times in the low-frequency region in case 1 (see

Table 4. Maximum motion increment of the WEC compared to the case of no breakwater.

Case Number	Motion Increment of WEC
	In Low-Frequency Range	In Resonance Range
Case 1	1.87	2.24
Case 2	1.70	2.04

). In addition, the maximum motion amplification occurred approximately 2.24 times near the resonance frequency region in the presence of a breakwater compared to the case of a WEC with no breakwater. This is caused by the interaction between the wall and the WEC. Additionally, it can be seen that the heave RAO converged to 2.0 in the low-frequency region.

The power generation of the WEC was calculated using the PTO force to determine the effect of motion amplification on power generation. In this study, the mean generated power was calculated using Equation (23) in [34], assuming that the PTO force is linearly proportional to the motion of the WEC.

$$P_{abs}=\frac{1}{2}{b}_{pto}{\omega }^2{\left|{\zeta }_3\right|}^2$$

(23)

$$\begin{gathered} CW=\frac{P_{abs}}{J} \\ CWR[\%]=\frac{CW}{B}\times 100 \end{gathered}$$

(24)

where $b_{pto}$ is the PTO damping coefficient. The capture width (CW) is defined as the ratio of the extracted wave power, $P_{abs}$ (in $W$), to the input wave energy, $J$(in $W/m$) [35]. Thus, CW is interpreted as the wave-front width completely absorbed by the WEC from the incoming wave. The capture width ratio (CWR) is the ratio of CW to the characteristic length of the WEC. Therefore, the CWR reflects the efficiency of the hydrodynamic absorption of a WEC. The characteristic length ($B$), which is generally used to calculate the CWR, is the width of the WEC [36]. In this study, $B$ is the same as the diameter of the WEC ($D$).

A case study was performed to determine the optimal PTO damping coefficient for each incident wave frequency. Figure 18 shows the contour of the heave RAO of the WEC according to $b_{pto}$ and $\omega$. As the PTO damping coefficient increases, the heave RAO of the WEC decreased significantly near the resonant frequency ($\omega$ = 3.2 rad/s). This is where the kinetic energy of the floating body is replaced with PTO damping (energy extraction). In the low-frequency region, the CWR increases as the PTO damping (extraction) coefficient increases (Figure 19). More incident wave energy is extracted at the lower-frequency wave region. With the breakwater, the CWR increases by approximately 3.67 times compared to without the breakwater. This is a significant improvement in the hydrodynamic efficiency of the WEC because of the breakwater. However, the CWR is very low in some frequency ranges ($\omega =2.0-2.8 \mathrm{rad}/\mathrm{s}$) compared to the WEC without the breakwater. This is because the wave reflected by the breakwater at $\omega =2.0-2.8 \mathrm{rad}/\mathrm{s}$ (when the WEC is located at the node of the standing wave) interferes with the movement of the WEC. This phenomenon has been reported in some studies [37,38]. Therefore, to increase the energy extraction efficiency of the WEC, it is necessary to design the WEC system to be integrated with the breakwater, considering the sea conditions of the installation area.

4. Conclusions

A three-dimensional frequency-domain numerical wave tank (FR-NWT) was developed for the analysis of a hemispherical heaving-buoy WEC integrated with a wall-type breakwater. The FR-NWT, a Fortran-based hydrodynamic in-house program, was based on the boundary element method with Rankine panels.

The ADZ scheme, which is mainly used for time-domain analysis programs, was first applied to the development of a 3D frequency-domain FR-NWT. In addition, to apply the ADZ of the optimal condition, the optimal length of ADZ, type of damping term, and shape function of ramping function were estimated. The ${\varphi }_n-\eta$-type artificial damping term showed better damping performance than the $\eta -\varphi$-type in the time-domain analysis in a previous study, but there was no significant difference between the two types in the FR-NWT.

The result using the ADZ scheme was more accurate than the result using the SRC as the radiation boundary condition in the rectangular computational domain of FR-NWT. Using the combined radiation condition (using both SRC and ADZ) for an effective far-field radiation condition, the length of ADZ could be reduced to 0.5 times the incident wavelength while maintaining the accuracy of the calculation. In addition, the calculation time could be shortened.

The dynamic response of a hemispherical heaving-buoy-type WEC integrated with a wall-type breakwater was analyzed using the developed FR-NWT, and the energy extraction efficiency was calculated. Linear PTO damping was used in this study. The motion amplification was greater when the WEC led at half the incident wavelength. Therefore, the distance between the WEC and the breakwater ($\alpha$) should be determined according to the wavelength of the target sea state.

According to various PTO damping coefficients, the heave RAO of WEC increased 1.87 times in the low-frequency region and 2.2 times in the resonance period. The CWR of the WEC increased as the motion of the WEC increased. The CWR of the WEC combined with the breakwater increased by 3.67 times compared to the case without the breakwater in a low-frequency area.

The presence of breakwaters can cause a negative effect on power generation in some frequency ranges. Thus, the energy extraction efficiency of a WEC can be improved by designing the floating body through integration with a breakwater, considering the sea states of the target installation area. The development of FR-NWT is a pre-emptive and essential study for the verification of various future studies, such as the WEC layout optimization and time-domain analysis, as well as the seafloor effect of floating bodies.