# Development of a Three-Dimensional Fully Nonlinear Potential Numerical Wave Tank for a Heaving Buoy Wave Energy Converter

Authors: Sung-Jae Kim; Weoncheol Koo
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The hydrodynamic performance of a vertical cylindrical heaving buoy-type floating wave energy converter under large-amplitude wave conditions was calculated. For this study, a three-dimensional fully nonlinear potential-flow numerical wave tank (3D-FN-PNWT) was developed. The 3D-FN-PNWT was based on the boundary element method with Rankine panels. Using the mixed Eulerian–Lagrangian (MEL) method for water particle movement, nonlinear waves were produced in the PNWT. The PNWT can calculate the wave forces acting on the buoy accurately using an acceleration potential approach. The constant panels and least-square gradient reconstruction method were applied to regridding of computational boundaries. An artificial damping zone was employed to satisfy the open-sea conditions at the end free surface boundaries. The diffraction and radiation problems were solved, and their solutions were confirmed by a comparison with previous studies. The interaction of the incident wave, floating body, and power take-off (PTO) behavior was examined in the time domain using the developed 3D-FN-PNWT. From comparison, the difference between the conventional linear analysis and the nonlinear analysis in large-amplitude waves was examined.


---

## Body

## 1. Introduction

Wave-body interactions are an essential part in analyzing the dynamical response of floating structures. Generally, wave nonlinearity increases with increasing ocean wave amplitude. Under high-wave conditions, various nonlinear phenomena of floating body motions can also occur. In particular, it is important to simulate a wave energy converter (WEC) considering wave-body interactions because the WEC can produce considerable energy from relatively high waves. Owing to the recent advances in computational fluid dynamics (CFD), it has become possible to develop viscous numerical wave tanks. On the contrary, a significantly large grid and considerable computational effort are required to account for wave-body interactions. The potential numerical wave tank (PNWT) technique, however, has a relatively short computation time and is capable of simulating a range of nonlinear phenomena of free surface behaviors with floating bodies.The PNWT technique can describe nonlinear waves based on the mixed Eulerian–Lagrangian (MEL) method proposed by Longuet-Higgins and Cokelet [1]. This method first analyzes the potential flow in the fluid domain from an Eulerian point of view and tracks the movement of the boundary nodes on the free surface and the buoy surface from a Lagrangian point of view. Based on this, nonlinear wave phenomena, wave-body interactions, and additional resistance of the ship can be applied to many floating body simulation problems.Many studies have examined wave-body interactions using the PNWT technique since the 1980s. The research can be classified into the following three categories: (i) diffraction problem for a buoy, (ii) radiation problem for a buoy, and (iii) freely floating problem of a buoy.For the diffraction problem, Yang and Ertekin [2] conducted a 3D-PNWT technique to estimate the wave loads acting on a fixed body. They considered the Sommerfeld [3] boundary condition as the radiation boundary condition and a fourth-order Adams–Bashforth integration method for a time-marching scheme. They used the Stokes second-order waves as the incident waves. Boo and Kim [4] also solved the diffraction problem for a cylinder in regular waves using a nonlinear PNWT technique. Based on this, Boo [5] described the diffracted waves caused by fixed structures under linear and nonlinear irregular wave conditions. For this, a high-order boundary element method and artificial damping zone for the external conditions were applied. Ferrant et al. [6] compared the results of the 2nd order diffraction approach and fully nonlinear analysis for a fixed cylindrical body. They confirmed that fully nonlinear analysis results were more accurate than the 2nd order results. Bai and Eatock Taylor [7] considered the 3D-PNWT based on the higher order boundary element method by dividing the computational domain for computational efficiency. They solved the diffraction problems in nonlinear regular and focused waves. By solving the diffraction problem, wave scattering due to an obstruction (buoy or bottom) and excitation body force can be predicted.For radiation problems, Ng and Issacson [8] and Issacson and Ng [9] solved the radiation problem in the time domain by applying Taylor series expansions and Stokes 2nd order perturbation approach. Wu and Eatock Taylor [10] applied fully nonlinear waves to a two-dimensional PNWT using the finite element method. Maiti and Sen [11] solved the radiation using two-dimensional cross sections of forced oscillating single-stranded and catamaran hulls by 2D-FN-PNWT analysis. Koo and Kim [12] solved the nonlinear waves generated by a wedge-type wave maker using a similar technique. Datta and Sen [13] applied the B-spline method to represent the body boundaries more clearly. Using the 3D-FN-PNWT, they solved the radiation and diffraction problems for a hemisphere and the Wigley hull. Bai and Eatock Taylor [14] also solved the same problem for a cylindrical body based on the higher order boundary element method. Zhou et al. [15] compared their results of various motion modes with those reported by Bai and Eatock Taylor [14]. Recently, Oh et al. [16] solved the radiation problem for a wave maker by applying the desingularized indirect PNWT. Feng et al. [17] solved the linear diffraction and radiation problems for a hemispheric buoy by calculating the integral of the Rankine source on the free surface panel.Koo and Kim [18, 19] examined the nonlinear hydrodynamic responses of various shaped buoys using the 2D-FN-PNWT technique with constant panels. They used Tanizawa’s acceleration potential method and mode decomposition method. Abbasnia and Ghiasi [20] examined the responses of two fully submerged cylinders using the 2D-FN-PNWT technique. A nonuniform rational B-spline (NURBS) was used to calculate the higher spatial derivatives of the boundaries of the body and free surface. Based on this, Abbasnia et al. [21, 22] carried out numerical analysis of a two-dimensional fully submerged body.Bai and Eatock Taylor [23] developed a three-dimensional high-order boundary element PNWT based on their previous studies [7, 14] to analyze the hydrodynamic performance of a cylindrical buoy and a flared vertical cylinder. The indirect method from Wu and Eatock Taylor [24] and the LU decomposition method for an inverse matrix calculation were adopted. Zhou et al. [25, 26] examined the response of a flared vertical buoy with the Stokes 5th order wave and solitary wave using the 3D-PNWT with higher order elements. They used the disturbed potential method proposed by Ferrant et al. [6]. Abbasnia and Soares [27] also examined the motion of the Wigley hull using the PNWT technique with the 3D NURBS method. In addition, the PNWT technique was used for wave resistance and wave pattern analysis for a ship with a forward speed [28, 29].Some studies of wave energy converters using the PNWT technology have been attempted. Lee et al. [30] and Kim et al. [31] conducted a numerical and experimental study of a backward bent duct buoy- (BBDB-) type floating wave energy converter by applying the 2D-FN-PNWT. However, wave energy converters can generate power mainly in the nonlinear condition (high wave amplitude and large body motion), so various nonlinear effects must be considered. As far as the authors know, no public literature has been available to numerically analyze a floating-type wave energy converter using the 3D-FN-PNWT to consider the coupling effect between wave-body-PTO (power take-off) systems.In this study, numerical modeling of a vertical cylindrical buoy was developed using a 3D-FN-PNWT technique. The constant panel method (CPM) was adopted. The CPM requires more elements for an accurate interpretation than the high-dimensional panel method, but the definition of the elements and calculation points are simple and numerically stable using the appropriate simulation model. The mixed Eulerian–Lagrangian (MEL) method was applied to simulate the nonlinear water particle motion effectively. The acceleration potential approach using the indirect method was used to calculate the total force and displacement of the body. To describe the body boundary condition in an acceleration field, the simplified formulation from Letournal et al. [32] was adopted. To develop the PNWT based on the CPM, the least-square gradient reconstruction method and the inverse distance weighting method were adopted. The least-square gradient reconstruction method, which is a well-known technique in CFD, was used to calculate the spatial derivatives of the elements in the CPM. In addition, the inverse distance weighting method was adopted for regridding of whole computational meshes.The developed 3D-FN-PNWT technique can be used for a range of hydrodynamic problems, such as wave generation and propagation, diffraction problem, radiation problem, and floating body dynamics. The wave-body-PTO interactions can also be modeled by appropriate numerical modeling of the PTO system. A hydraulic-cylinder-type PTO system, which can be applied mainly to a moving-body-type point absorber, was adopted in this study. Kim et al. [33] provided details and verification of the numerical modeling of a PTO system.An analysis of a fully nonlinear wave-body-PTO interaction was performed. The nonlinear effect was examined by comparing the numerical results on various wave steepness conditions.


## 2. Mathematical Formulation

### 2.1. Rankine Panel Method

In the hydrodynamic analysis of ships and ocean structures, the effect of viscosity is generally not significant and the analysis has been carried out under the assumption of potential flow, which is an inviscid, irrotational, and incompressible fluid. The governing equation for the potential fluid domain is the Laplace equation. Using Green’s 2nd identity, the governing equation can be transformed into the boundary integral equation as follows:(1)αϕ=∫SG∂ϕ∂n−ϕ∂G∂ndS,where α is the solid angle and G is the Green function, which can be expressed as(2)G=1r,where r is the distance between the source point and the field point on the boundary. The influence matrix using the Green function in the CPM is constructed according to r because the potential value on the element does not change. The method reported by Hess and Smith [34] can be applied when r is small. When r is large, however, the multipole expansion method is used [35, 36]. The source potential and dipole potential by the multipole expansion method can be expressed as follows:(3)∬SGdS=∑m=0∞∑n=0∞−1m+nm!n!Imn∂m+n∂xm∂ynx2+y2+z2−1/2,∬S∂G∂ndS=∑m=0∞∑n=0∞−1m+n+1m!n!Imnz∂m+n∂xm∂ynx2+y2+z2−3/2,where Imn denotes the moment of the panel about the origin. The subscripts m and n mean the orders for the x- and y-axes, respectively. The computational time is reduced considerably by calculating the moment of each element in advance and performing the influence-matrix calculation. All calculations were conducted by CPU E5-2620 2.40 GHz without a parallel computation. In the case without the multipole expansion method, the method reported by Hess and Smith [34] was applied to calculate the influence matrix, irrespective of the magnitude of r. The computation time of the influence matrix reduced the time by a factor of four compared to the case where the multipole extension method was not considered.Figure1 gives an overview of a 3D-FN-PNWT. The numerical simulation of the NWT is a type of mixed boundary value problem with various boundary conditions. The computational fluid domain was surrounded by the incident wave boundary, free surface boundary, body boundary, side wall and back wall, and bottom boundaries. The image method was used for the bottom boundary, and the computational nodes were not required on the bottom. The Stokes 2nd order wave potential profile was applied to the incident wave boundary condition, and the angle of incident waves was set to zero:(4)∂ϕ∂n=nx∂ϕ∂x=nxgAkωcoshkz+hcoshkhcoskx−ωtlinear+nx38A2ω2kcosh2kz+hsinh4khcos2kx−ωtStokes2nd,onSin,where g, A, k, ω, and h are the gravitational acceleration, wave amplitude, wave number, wave frequency, and water depth, respectively. n→=nxi→+nyj→+nzk→ is the normal vector on the boundary surface. The kinematic and dynamic boundary conditions in the following equations were satisfied on the free surface boundary [18]:(5)δηδt=∂ϕ∂z−∇ϕ−v→⋅∇η,onSF,(6)δϕδt=−gη−12∇ϕ2+v→⋅∇ϕ,onSF,where v→ denotes the velocity of the calculation node on the surface and δ/δt and ∂/∂t denote the material derivative (=∂/∂t+u⋅∇) and the partial derivative, respectively. Depending on how the node velocity is expressed, there are two types of time integration for the free surface boundary conditions, such as the full-Lagrangian method and semi-Lagrangian method. The full-Lagrangian method expresses the velocity of a computational node the same as the velocity of a water particle on the free surface (v→=ξ˙,ζ˙,η˙). In the semi-Lagrangian method, the node velocity moves only in the vertical direction of the water particles on the free surface (v→=0,0,η˙). The full-Lagrangian method can describe nonlinear waves directly in terms of the velocity of water particles in the x, y, and z directions. On the contrary, it is necessary to rearrange the positions of the calculation points in the x and y directions. For the semi-Lagrangian method, rearranging the positions of the computational nodes in the x direction at the free surface boundaries is unnecessary because the position of the node moves only vertically. In this study, the semi-Lagrangian method was applied, and the nonlinear free surface boundary condition can be expressed as(7)δηδt=∂ϕ∂z−∂ϕ∂x∂η∂x−∂ϕ∂y∂η∂y,onSF,(8)δϕδt=−gη−12∇ϕ2+δηδt∂ϕ∂z,onSF.Figure 1
Overview of a three-dimensional potential numerical wave tank.The impervious boundary condition given in the following condition was applied to the side and back walls:(9)∂ϕ∂n=0,onSW.The buoy boundary condition is also expressed as an impermeable boundary condition (equation (9)) when the diffraction problem is solved. When the buoy is a floating body, the boundary condition can be expressed by the following equation in the radiation problem:(10)∂ϕ∂n=V→⋅n→,onSB,where V→ and n→ denote the buoy velocity and normal vector of each boundary element, respectively.


### 2.2. Artificial Damping Zone

To realize open-sea conditions, an artificial damping zone was installed on the free surface boundary near the end wall. The artificial attenuation scheme has been applied widely to various types of waves, such as nonlinear waves, linear waves, regular waves, irregular waves, and multidirectional waves in the numerical simulation. Many studies on forced wave damping have been carried out since the 1980s [18, 37–41]. In particular, Kim et al. [42] compared and analyzed a range of artificial damping terms and ramp functions using a 3D-FN-PNWT technique. In this study, the ϕn−η-type artificial damping term and the low slope at the front of the damping zone (equation (13)) were used, as reported by Kim et al. [42]. Equations (11) and (12) show the kinematic and dynamic free surface boundary conditions with a damping term for each boundary condition to prevent the reflected waves due to the end wall boundary. The ramp function used in the damping terms is given in equation (13) [42]:(11)δηδt=∂ϕ∂z−∂ϕ∂x∂η∂x−∂ϕ∂y∂η∂y−μE2η,onSFE,(12)δϕδt=−gη−12∇ϕ2+δηδt∂ϕ∂z−μE1∂ϕ∂n,onSFE,(13)μE1x=μE101−cosπ2x−xedld,where μE1 and μE2 are the end artificial damping coefficients (μE2=kμE1) and SFE, μE10, ld, and xed denote the free surface boundary with the end artificial damping zone (see Figure 1), magnitude of the end artificial damping coefficient, length of the damping zone, and starting location of the damping zone, respectively. In this study, the length of the artificial damping zone was set to one wavelength. In addition, the frontal artificial damping zone was also applied in front of the input boundary. This can prevent the re-reflected waves caused by the side wall, body, and radiated waves. On the contrary, in order to not affect the incident wave, the damping term was applied with the reference values (ϕn∗ and η∗) calculated by the theoretical solution of the incident wave. The following equations show the boundary conditions for the frontal artificial damping zone [18]:(14)δηδt=∂ϕ∂z−∂ϕ∂x∂η∂x−∂ϕ∂y∂η∂y−μF2η−η∗,onSFF,(15)δϕδt=−gη−12∇ϕ2+δηδt∂ϕ∂z−μF1∂ϕ∂n−∂ϕ∂n∗,onSFF,where μF1 and μF2 denote the frontal damping coefficients and SFF is the free surface boundary with the frontal artificial damping zone (see Figure 1). The length of the damping zone and ramp function are the same as those of the end damping zone:(16)μF1x=μF101−cosπ2xfd−xld,where μF10 and xfd denote the frontal artificial damping coefficient and the starting point of the frontal damping zone.To avoid reflection from the side wall, the side damping zone was also installed (equations (17) and (18)). The side damping zone was designed to operate from the moment the incident wave reached at a nodal point, x. Equation (19) expresses the time at which the incident wave reached at point x. Equation (20) shows the switch function to apply the damping term according to the arrival point of the wave. The artificial damping coefficient is gradually increased by the temporal ramp function from the moment the incident wave reaches a certain point:(17)δϕδt=−gη−12∇ϕ2+v→⋅∇ϕ−κ1μS1∂ϕ∂n−∂ϕ∂n∗,(18)δηδt=∂ϕ∂z−∇ϕ−v→⋅∇η−κ1μS2η−η∗,(19)TS0=xcg,(20)κ1=0,fort<TS0,1−cosπt−TS02TRamp,forTS0≤t,1,forTS1<t,where μS1 and μS2 are the side artificial damping coefficients and κ1, TRamp, TS0, and TS1 are the switch function, the duration of the ramp function, the time at which the incident wave is reached, and the time when the damping terms are fully applied, respectively. In addition, the spatial ramp function was applied as(21)μS1x=μS101−cosπ2y−ysdld,where μS10 and ysd denote the frontal artificial damping coefficient and the starting point of the side damping zone.


### 2.3. Least-Square Gradient Reconstruction Method

The Green–Gauss method and least-square gradient reconstruction method are typical tools for calculating the gradient of a certain property for an unstructured grid. The least-square gradient reconstruction method can express the specific physical quantity “Qi” at the neighboring calculation point “i” based on a specific calculation point “0” using a Taylor series transformation [43]:(22)Qi=Q0+∂Q∂x0Δxi+∂Q∂y0Δyi+∂Q∂z0Δzi+∂2Q∂x2012Δxi2+∂2Q∂y2012Δyi2+∂2Q∂x2012Δyi2+∂2Q∂x∂y0ΔxΔy+∂2Q∂y∂z0ΔyΔz+∂2Q∂x∂z0ΔxΔz,where Δxi, Δyi, and Δzi are the difference of the x, y, and z coordinates between a specific calculation point “0” and the neighboring calculation point “i,” respectively. In this study, the wave height at the free surface and the velocity potential at the free surface and buoy boundaries were considered in the physical quantities, “Q.” The weight term was applied using the distance between the specific calculation point “0” and the neighboring point “i,” and the difference between the left side and the right side of equation (22) can be expressed as(23)fwi,∂Q∂x0,∂Q∂y0,∂Q∂z0,∂2Q∂x20,∂2Q∂y20,∂2Q∂z20,∂2Q∂x∂y0,∂2Q∂y∂z0,∂2Q∂x∂z0=∑i=1NwiΔQi−∂Q∂x0Δxi+∂Q∂y0Δyi+∂Q∂z0Δzi+∂2Q∂x2012Δxi2+∂2Q∂y2012Δyi2+∂2Q∂z2012Δzi2+∂2Q∂x∂y0ΔxΔy+∂2Q∂y∂z0ΔyΔz+∂2Q∂x∂z0ΔxΔz2,where wi denotes the weight for the specific calculation point “0” of the ith neighboring calculation point, and it can be expressed by the reciprocal of the fifth power of the distance from the ith calculation point to the specific calculation point. More details are given in Appendix.


### 2.4. Acceleration Potential Approach

To obtain the pressure and force acting on the buoy due to wave loading, it is essential to estimate the time derivative of the velocity potential (ϕt=∂ϕ/∂t) accurately. Therefore, the acceleration potential approach has been used [18, 44]. This approach solved the Laplace equation with boundary conditions in the acceleration potential field, and the time derivative of the velocity potential was obtained as a solution. The following equation shows the definition of the acceleration term [18, 44]:(24)a→=δ∇ϕδt=∂∇ϕ∂t+∇ϕ⋅∇∇ϕ=∇∂ϕ∂t+12∇ϕ2=∇Φ,where Φ denotes the acceleration potential. On the contrary, this term cannot satisfy the Laplace equation. Hence, the time derivative of the velocity potential can be obtained directly from the following equation to satisfy the Laplace equation:(25)ϕt=∂ϕ∂t=Φ−12∇ϕ2.The body boundary condition is difficult to describe in the acceleration potential field because of the nonlinear termϕt. The expression reported by Letournel et al. [32] was adopted to express the body boundary condition, which was derived from the approach of Tanizawa [44] and Cointe [45]. The following equations show the body boundary condition in the acceleration field. Koo and Kim [18] reported the detailed boundary conditions presented here:(26)∂ϕt∂n=n→⋅a→0+ω˙→×r→+qB,(27)qB=ω→⋅s→1∂ϕ∂s2−2v→⋅s→2−ω→⋅s→2∂ϕ∂s1−2v→⋅s→1+kn1v→⋅s→1∂ϕ∂s1−v→⋅s→1+kn2v→⋅s→2∂ϕ∂s2−v→⋅s→2+v→⋅n→∂2ϕ∂s12+∂2ϕ∂s22+kn1+kn2∂ϕ∂n,where a0, ω˙→, and r→ denote the translate acceleration of the buoy, angular acceleration of the buoy, and local coordinate vector on the body boundary, s→1 and s→2 are the tangent vectors of the quadratic panel, and kn1 and kn2 are the radius of gyration of each tangent vector.To calculate the buoy force and displacement using the acceleration potential method, the indirect calculation method was applied together. To calculate the time derivative of velocity potential on a floating body surface, the body boundary condition should be separated because it contains the body acceleration. So, the time derivative of velocity potential can be expressed as(28)∂ϕ∂t=∑i=16aiφi+φ7,where ai and φi denote the body acceleration and i-th mode component of the acceleration potential from the i-th mode motion. Modes with “i” less than 4 mean translational motion, and modes over 4 mean rotational motion. φ7 is the acceleration potential component from the velocity field. The indirect method is a method of calculating φ7 by applying the virtual velocity potential (ψi) satisfying the governing equation (equation (29)) with the Green 2nd identity (equation (30)):(29)∇2ψi=0,(30)∬Sφ7∂ψi∂n−ψi∂φ7∂nds=0.Substituting the boundary condition into equation (30) yields(31)∬SBφ7nids=∬SBψiqBds−∬SF∂ψi∂n−gη−12∇ϕ2ds+∬Sinψi∂ϕt∂nds.This method cannot calculate the exact hydrodynamic pressure acting on each element, but it is a quick way to calculate the force acting on the buoy. The following equations express the pressure and total force acting on the buoy, respectively:(32)P=−ρ∑i=16aiφi+φ7+gz+12∇ϕ2,(33)mai=Fi=∫SBPnids−Wδi3−FPTO,where W and FPTO are the buoy weight and PTO force, respectively.


### 2.5. Inverse Distance Weighting

For regridding of meshes, the inverse distance weighting method (IDW) was adopted. It is one of the conventional methods for multivariate interpolation with a known scattered quantity of points. This is a method applied to the geographic computation field and has recently been applied to the CFD analysis [46].(34)Zp=∑i=1nZiwi∑i=1nwi,where Zi and wi are the known quantity of the ith known calculation point and the weighting of the ith known point for the target point and Zp is the unknown quantity for the target point. The weighting was calculated by the reciprocal of the fourth power of the distance from the ith calculation point to the target point.For regridding of meshes, free surface calculation nodes are relocated first by using the MEL method. Then, IDW is carried out to redefine the free surface calculation panels. Based on the free surface panels and buoy displacement, the other nodes and panels are reestablished.


### 2.6. Additional Numerical Schemes

A range of numerical techniques have been applied for stable numerical simulations. To reduce the transient mode caused by an abrupt change due to wave generation or body motion, the ramp function shown in equation (35) was applied to the incident wave boundary condition. For ramping time, 4T (four wave periods) was adopted.(35)rt=1,fort>TRamp,1−cosπ2tTRamp,fort≤TRamp.A nonphysical saw tooth instability can occur when simulating the nonlinear waves on the free surface. The Chebyshev five-point smoothing technique and B-spline technique were used as a smoothing technique to avoid the instability phenomenon [1, 18, 47].The positions of the calculation nodes on the boundary in the fully nonlinear numerical simulation should be rearranged according to their physical positions at every time step. The influence function for the computational domain was then varied at every time step. To convert the inverse matrix of the influence function for entire boundary elements, the generalized minimal residual (GMRES) method was used at every time step [48]. The GMRES method is an iterative solver for a linear algebraic computation. The right preconditioning was also applied to minimize the residuals quickly.


### 2.7. PTO Numerical Modeling

The numerical model of a hydraulic PTO system was developed as expressed in equation (32). This PTO system is generally used in a heaving buoy-type WEC ([33, 49]).(36)FPTO=−minGz˙i,signz˙iScΔp,where Sc and Δp are the sectional area of the hydraulic cylinder and difference between HP (high pressure) and LP (low pressure), respectively, and G is the numerical slope to reduce numerical errors. The value is set to five times the excitation force on the body that can be obtained from the diffraction problem.


## 2.1. Rankine Panel Method

In the hydrodynamic analysis of ships and ocean structures, the effect of viscosity is generally not significant and the analysis has been carried out under the assumption of potential flow, which is an inviscid, irrotational, and incompressible fluid. The governing equation for the potential fluid domain is the Laplace equation. Using Green’s 2nd identity, the governing equation can be transformed into the boundary integral equation as follows:(1)αϕ=∫SG∂ϕ∂n−ϕ∂G∂ndS,where α is the solid angle and G is the Green function, which can be expressed as(2)G=1r,where r is the distance between the source point and the field point on the boundary. The influence matrix using the Green function in the CPM is constructed according to r because the potential value on the element does not change. The method reported by Hess and Smith [34] can be applied when r is small. When r is large, however, the multipole expansion method is used [35, 36]. The source potential and dipole potential by the multipole expansion method can be expressed as follows:(3)∬SGdS=∑m=0∞∑n=0∞−1m+nm!n!Imn∂m+n∂xm∂ynx2+y2+z2−1/2,∬S∂G∂ndS=∑m=0∞∑n=0∞−1m+n+1m!n!Imnz∂m+n∂xm∂ynx2+y2+z2−3/2,where Imn denotes the moment of the panel about the origin. The subscripts m and n mean the orders for the x- and y-axes, respectively. The computational time is reduced considerably by calculating the moment of each element in advance and performing the influence-matrix calculation. All calculations were conducted by CPU E5-2620 2.40 GHz without a parallel computation. In the case without the multipole expansion method, the method reported by Hess and Smith [34] was applied to calculate the influence matrix, irrespective of the magnitude of r. The computation time of the influence matrix reduced the time by a factor of four compared to the case where the multipole extension method was not considered.Figure1 gives an overview of a 3D-FN-PNWT. The numerical simulation of the NWT is a type of mixed boundary value problem with various boundary conditions. The computational fluid domain was surrounded by the incident wave boundary, free surface boundary, body boundary, side wall and back wall, and bottom boundaries. The image method was used for the bottom boundary, and the computational nodes were not required on the bottom. The Stokes 2nd order wave potential profile was applied to the incident wave boundary condition, and the angle of incident waves was set to zero:(4)∂ϕ∂n=nx∂ϕ∂x=nxgAkωcoshkz+hcoshkhcoskx−ωtlinear+nx38A2ω2kcosh2kz+hsinh4khcos2kx−ωtStokes2nd,onSin,where g, A, k, ω, and h are the gravitational acceleration, wave amplitude, wave number, wave frequency, and water depth, respectively. n→=nxi→+nyj→+nzk→ is the normal vector on the boundary surface. The kinematic and dynamic boundary conditions in the following equations were satisfied on the free surface boundary [18]:(5)δηδt=∂ϕ∂z−∇ϕ−v→⋅∇η,onSF,(6)δϕδt=−gη−12∇ϕ2+v→⋅∇ϕ,onSF,where v→ denotes the velocity of the calculation node on the surface and δ/δt and ∂/∂t denote the material derivative (=∂/∂t+u⋅∇) and the partial derivative, respectively. Depending on how the node velocity is expressed, there are two types of time integration for the free surface boundary conditions, such as the full-Lagrangian method and semi-Lagrangian method. The full-Lagrangian method expresses the velocity of a computational node the same as the velocity of a water particle on the free surface (v→=ξ˙,ζ˙,η˙). In the semi-Lagrangian method, the node velocity moves only in the vertical direction of the water particles on the free surface (v→=0,0,η˙). The full-Lagrangian method can describe nonlinear waves directly in terms of the velocity of water particles in the x, y, and z directions. On the contrary, it is necessary to rearrange the positions of the calculation points in the x and y directions. For the semi-Lagrangian method, rearranging the positions of the computational nodes in the x direction at the free surface boundaries is unnecessary because the position of the node moves only vertically. In this study, the semi-Lagrangian method was applied, and the nonlinear free surface boundary condition can be expressed as(7)δηδt=∂ϕ∂z−∂ϕ∂x∂η∂x−∂ϕ∂y∂η∂y,onSF,(8)δϕδt=−gη−12∇ϕ2+δηδt∂ϕ∂z,onSF.Figure 1
Overview of a three-dimensional potential numerical wave tank.The impervious boundary condition given in the following condition was applied to the side and back walls:(9)∂ϕ∂n=0,onSW.The buoy boundary condition is also expressed as an impermeable boundary condition (equation (9)) when the diffraction problem is solved. When the buoy is a floating body, the boundary condition can be expressed by the following equation in the radiation problem:(10)∂ϕ∂n=V→⋅n→,onSB,where V→ and n→ denote the buoy velocity and normal vector of each boundary element, respectively.


## 2.2. Artificial Damping Zone

To realize open-sea conditions, an artificial damping zone was installed on the free surface boundary near the end wall. The artificial attenuation scheme has been applied widely to various types of waves, such as nonlinear waves, linear waves, regular waves, irregular waves, and multidirectional waves in the numerical simulation. Many studies on forced wave damping have been carried out since the 1980s [18, 37–41]. In particular, Kim et al. [42] compared and analyzed a range of artificial damping terms and ramp functions using a 3D-FN-PNWT technique. In this study, the ϕn−η-type artificial damping term and the low slope at the front of the damping zone (equation (13)) were used, as reported by Kim et al. [42]. Equations (11) and (12) show the kinematic and dynamic free surface boundary conditions with a damping term for each boundary condition to prevent the reflected waves due to the end wall boundary. The ramp function used in the damping terms is given in equation (13) [42]:(11)δηδt=∂ϕ∂z−∂ϕ∂x∂η∂x−∂ϕ∂y∂η∂y−μE2η,onSFE,(12)δϕδt=−gη−12∇ϕ2+δηδt∂ϕ∂z−μE1∂ϕ∂n,onSFE,(13)μE1x=μE101−cosπ2x−xedld,where μE1 and μE2 are the end artificial damping coefficients (μE2=kμE1) and SFE, μE10, ld, and xed denote the free surface boundary with the end artificial damping zone (see Figure 1), magnitude of the end artificial damping coefficient, length of the damping zone, and starting location of the damping zone, respectively. In this study, the length of the artificial damping zone was set to one wavelength. In addition, the frontal artificial damping zone was also applied in front of the input boundary. This can prevent the re-reflected waves caused by the side wall, body, and radiated waves. On the contrary, in order to not affect the incident wave, the damping term was applied with the reference values (ϕn∗ and η∗) calculated by the theoretical solution of the incident wave. The following equations show the boundary conditions for the frontal artificial damping zone [18]:(14)δηδt=∂ϕ∂z−∂ϕ∂x∂η∂x−∂ϕ∂y∂η∂y−μF2η−η∗,onSFF,(15)δϕδt=−gη−12∇ϕ2+δηδt∂ϕ∂z−μF1∂ϕ∂n−∂ϕ∂n∗,onSFF,where μF1 and μF2 denote the frontal damping coefficients and SFF is the free surface boundary with the frontal artificial damping zone (see Figure 1). The length of the damping zone and ramp function are the same as those of the end damping zone:(16)μF1x=μF101−cosπ2xfd−xld,where μF10 and xfd denote the frontal artificial damping coefficient and the starting point of the frontal damping zone.To avoid reflection from the side wall, the side damping zone was also installed (equations (17) and (18)). The side damping zone was designed to operate from the moment the incident wave reached at a nodal point, x. Equation (19) expresses the time at which the incident wave reached at point x. Equation (20) shows the switch function to apply the damping term according to the arrival point of the wave. The artificial damping coefficient is gradually increased by the temporal ramp function from the moment the incident wave reaches a certain point:(17)δϕδt=−gη−12∇ϕ2+v→⋅∇ϕ−κ1μS1∂ϕ∂n−∂ϕ∂n∗,(18)δηδt=∂ϕ∂z−∇ϕ−v→⋅∇η−κ1μS2η−η∗,(19)TS0=xcg,(20)κ1=0,fort<TS0,1−cosπt−TS02TRamp,forTS0≤t,1,forTS1<t,where μS1 and μS2 are the side artificial damping coefficients and κ1, TRamp, TS0, and TS1 are the switch function, the duration of the ramp function, the time at which the incident wave is reached, and the time when the damping terms are fully applied, respectively. In addition, the spatial ramp function was applied as(21)μS1x=μS101−cosπ2y−ysdld,where μS10 and ysd denote the frontal artificial damping coefficient and the starting point of the side damping zone.


## 2.3. Least-Square Gradient Reconstruction Method

The Green–Gauss method and least-square gradient reconstruction method are typical tools for calculating the gradient of a certain property for an unstructured grid. The least-square gradient reconstruction method can express the specific physical quantity “Qi” at the neighboring calculation point “i” based on a specific calculation point “0” using a Taylor series transformation [43]:(22)Qi=Q0+∂Q∂x0Δxi+∂Q∂y0Δyi+∂Q∂z0Δzi+∂2Q∂x2012Δxi2+∂2Q∂y2012Δyi2+∂2Q∂x2012Δyi2+∂2Q∂x∂y0ΔxΔy+∂2Q∂y∂z0ΔyΔz+∂2Q∂x∂z0ΔxΔz,where Δxi, Δyi, and Δzi are the difference of the x, y, and z coordinates between a specific calculation point “0” and the neighboring calculation point “i,” respectively. In this study, the wave height at the free surface and the velocity potential at the free surface and buoy boundaries were considered in the physical quantities, “Q.” The weight term was applied using the distance between the specific calculation point “0” and the neighboring point “i,” and the difference between the left side and the right side of equation (22) can be expressed as(23)fwi,∂Q∂x0,∂Q∂y0,∂Q∂z0,∂2Q∂x20,∂2Q∂y20,∂2Q∂z20,∂2Q∂x∂y0,∂2Q∂y∂z0,∂2Q∂x∂z0=∑i=1NwiΔQi−∂Q∂x0Δxi+∂Q∂y0Δyi+∂Q∂z0Δzi+∂2Q∂x2012Δxi2+∂2Q∂y2012Δyi2+∂2Q∂z2012Δzi2+∂2Q∂x∂y0ΔxΔy+∂2Q∂y∂z0ΔyΔz+∂2Q∂x∂z0ΔxΔz2,where wi denotes the weight for the specific calculation point “0” of the ith neighboring calculation point, and it can be expressed by the reciprocal of the fifth power of the distance from the ith calculation point to the specific calculation point. More details are given in Appendix.


## 2.4. Acceleration Potential Approach

To obtain the pressure and force acting on the buoy due to wave loading, it is essential to estimate the time derivative of the velocity potential (ϕt=∂ϕ/∂t) accurately. Therefore, the acceleration potential approach has been used [18, 44]. This approach solved the Laplace equation with boundary conditions in the acceleration potential field, and the time derivative of the velocity potential was obtained as a solution. The following equation shows the definition of the acceleration term [18, 44]:(24)a→=δ∇ϕδt=∂∇ϕ∂t+∇ϕ⋅∇∇ϕ=∇∂ϕ∂t+12∇ϕ2=∇Φ,where Φ denotes the acceleration potential. On the contrary, this term cannot satisfy the Laplace equation. Hence, the time derivative of the velocity potential can be obtained directly from the following equation to satisfy the Laplace equation:(25)ϕt=∂ϕ∂t=Φ−12∇ϕ2.The body boundary condition is difficult to describe in the acceleration potential field because of the nonlinear termϕt. The expression reported by Letournel et al. [32] was adopted to express the body boundary condition, which was derived from the approach of Tanizawa [44] and Cointe [45]. The following equations show the body boundary condition in the acceleration field. Koo and Kim [18] reported the detailed boundary conditions presented here:(26)∂ϕt∂n=n→⋅a→0+ω˙→×r→+qB,(27)qB=ω→⋅s→1∂ϕ∂s2−2v→⋅s→2−ω→⋅s→2∂ϕ∂s1−2v→⋅s→1+kn1v→⋅s→1∂ϕ∂s1−v→⋅s→1+kn2v→⋅s→2∂ϕ∂s2−v→⋅s→2+v→⋅n→∂2ϕ∂s12+∂2ϕ∂s22+kn1+kn2∂ϕ∂n,where a0, ω˙→, and r→ denote the translate acceleration of the buoy, angular acceleration of the buoy, and local coordinate vector on the body boundary, s→1 and s→2 are the tangent vectors of the quadratic panel, and kn1 and kn2 are the radius of gyration of each tangent vector.To calculate the buoy force and displacement using the acceleration potential method, the indirect calculation method was applied together. To calculate the time derivative of velocity potential on a floating body surface, the body boundary condition should be separated because it contains the body acceleration. So, the time derivative of velocity potential can be expressed as(28)∂ϕ∂t=∑i=16aiφi+φ7,where ai and φi denote the body acceleration and i-th mode component of the acceleration potential from the i-th mode motion. Modes with “i” less than 4 mean translational motion, and modes over 4 mean rotational motion. φ7 is the acceleration potential component from the velocity field. The indirect method is a method of calculating φ7 by applying the virtual velocity potential (ψi) satisfying the governing equation (equation (29)) with the Green 2nd identity (equation (30)):(29)∇2ψi=0,(30)∬Sφ7∂ψi∂n−ψi∂φ7∂nds=0.Substituting the boundary condition into equation (30) yields(31)∬SBφ7nids=∬SBψiqBds−∬SF∂ψi∂n−gη−12∇ϕ2ds+∬Sinψi∂ϕt∂nds.This method cannot calculate the exact hydrodynamic pressure acting on each element, but it is a quick way to calculate the force acting on the buoy. The following equations express the pressure and total force acting on the buoy, respectively:(32)P=−ρ∑i=16aiφi+φ7+gz+12∇ϕ2,(33)mai=Fi=∫SBPnids−Wδi3−FPTO,where W and FPTO are the buoy weight and PTO force, respectively.


## 2.5. Inverse Distance Weighting

For regridding of meshes, the inverse distance weighting method (IDW) was adopted. It is one of the conventional methods for multivariate interpolation with a known scattered quantity of points. This is a method applied to the geographic computation field and has recently been applied to the CFD analysis [46].(34)Zp=∑i=1nZiwi∑i=1nwi,where Zi and wi are the known quantity of the ith known calculation point and the weighting of the ith known point for the target point and Zp is the unknown quantity for the target point. The weighting was calculated by the reciprocal of the fourth power of the distance from the ith calculation point to the target point.For regridding of meshes, free surface calculation nodes are relocated first by using the MEL method. Then, IDW is carried out to redefine the free surface calculation panels. Based on the free surface panels and buoy displacement, the other nodes and panels are reestablished.


## 2.6. Additional Numerical Schemes

A range of numerical techniques have been applied for stable numerical simulations. To reduce the transient mode caused by an abrupt change due to wave generation or body motion, the ramp function shown in equation (35) was applied to the incident wave boundary condition. For ramping time, 4T (four wave periods) was adopted.(35)rt=1,fort>TRamp,1−cosπ2tTRamp,fort≤TRamp.A nonphysical saw tooth instability can occur when simulating the nonlinear waves on the free surface. The Chebyshev five-point smoothing technique and B-spline technique were used as a smoothing technique to avoid the instability phenomenon [1, 18, 47].The positions of the calculation nodes on the boundary in the fully nonlinear numerical simulation should be rearranged according to their physical positions at every time step. The influence function for the computational domain was then varied at every time step. To convert the inverse matrix of the influence function for entire boundary elements, the generalized minimal residual (GMRES) method was used at every time step [48]. The GMRES method is an iterative solver for a linear algebraic computation. The right preconditioning was also applied to minimize the residuals quickly.


## 2.7. PTO Numerical Modeling

The numerical model of a hydraulic PTO system was developed as expressed in equation (32). This PTO system is generally used in a heaving buoy-type WEC ([33, 49]).(36)FPTO=−minGz˙i,signz˙iScΔp,where Sc and Δp are the sectional area of the hydraulic cylinder and difference between HP (high pressure) and LP (low pressure), respectively, and G is the numerical slope to reduce numerical errors. The value is set to five times the excitation force on the body that can be obtained from the diffraction problem.


## 3. Numerical Results

### 3.1. Wave Generation and Propagation

Numerical simulations of wave generation and propagation were performed using the 3D-FN-PNWT technique. The dimensions of the 3D-FN-PNWT were six wavelengths in the longitudinal direction, one wavelength in the width direction, and 0.5 wavelengths in the vertical direction. The incident wave satisfied the deep-water condition because the water depth was set to 0.5 wavelengths in the vertical direction. The artificial damping zone, in which the length was two wavelengths, was installed near the end wall along the wave propagation direction.For input waves, the Stokes 2nd order waves can be expressed from the theoretical solution for an arbitrary depth as [50](37)ηt=Acos−ωt+kA3−tanh2kh4tanh3khcos−2ωt.Figure2 shows the time series of free surface elevations from numerical simulations and theoretical solutions. The input wave conditions were as follows: the wave period of three seconds, the wave height of 0.4 m, and the water depth of 3 m. The wave elevations from the 3D-FN-PNWT simulation agreed well with the theoretical Stokes 2nd order waves. Because the wave steepness (H/λ) was 0.0315 (≈1/32), the mean wave elevation was greater than zero because of wave nonlinearity.Figure 2
Comparison of the time series of wave elevations from the numerical and theoretical results (num: numerical results; theory: theoretical results).Figure3 shows the time series of the free surface elevations on various measured points and various wave steepness conditions. Wave elevations measuring all computational free surface points (four wavelengths) were the same without any reflections from the end wall. The amplitudes of the waves decreased at a point larger than four times the wavelength at which the artificial damping zone started. These results show that the artificial damping condition was well operated and the numerical wave tank satisfied the open-sea conditions. In Figure 3(b), when the wave steepness is small (H/λ=1/60), a linear wave with the same crest height and trough height is generated. On the contrary, as wave steepness increases, the crest wave height increases and is greater than the trough height. This is because the higher order term of the wave increases with increasing wave steepness.Figure 3
Time series of wave elevations on various measurement points (a) and various wave steepness conditions (b).
(a)
(b)


### 3.2. Diffraction Problem

To solve the diffraction problem of wave-body interactions, the Stokes 2nd order waves as the input condition were generated and propagated to the constrained body. The acceleration potential method was used to calculate the excitation force acting on the body. The dimensions of the computational domain were four wavelengths in the x direction, two wavelengths in the y direction, and 0.5 wavelengths in the z direction. The length of the artificial damping zone was set to one wavelength. Figure 4 shows meshes to present a cylindrical body and free surface around the body. The computational domain was modeled as a half model with the x-axis symmetry. All meshes in the whole calculational fluid domain were expressed using quadrilateral panels. The free surface meshes around the body were constructed radially from the body meshes. The body was a vertical cylinder with a draft and radius of 1 m each.Figure 4
Meshes to present a cylindrical body and free surface.Figure5 shows the time series of horizontal and vertical forces acting on a stationary vertical circular cylinder by solving the diffraction problem. The horizontal and vertical excitation forces were divided with ρgR2A, where R is the radius of the cylindrical buoy. As like a physical wave tank experiment, the hydrodynamic forces acting on the buoy occur after three times the wave period, which is the time when the incoming wave reaches the buoy. When the simulation time exceeds 8 times the wave period, the time series of both hydrodynamic forces are in a steady state. From this, the amplitude of horizontal and vertical forces can be obtained and compared with the results of the frequency-domain analysis. Figure 6 shows the convergence tests for the number of nodes per wavelength and interval of the time step. The vertical excitation forces on various kR conditions were compared. As the number of nodes increases and the interval of the time step decreases, the results of hydrodynamic analysis are converged. Therefore, optimal simulation conditions are adopted in this study (the number of nodes per wavelength is 20, and the interval of the time step is T/64).Figure 5
Time series of horizontal and vertical forces acting on a stationary vertical circular cylinder from the diffraction problem (kR = 0.44,A = 0.1 m, and R = d = 1.0 m under the deep-water condition).Figure 6
Convergence test for the number of nodes per wavelength (a) and time step (b) (R = d = 1.0 m under the deep-water condition).
(a)
(b)Figure7 compares the vertical and horizontal excitation forces acting on the vertical cylinder with the results of WAMIT based on the wave’s Green function. The draft and radius of the cylinder were 1 m each. The numerical results of the linear time-domain analysis (Time (Lin)) and fully nonlinear time-domain analysis using the 3D-FN-PNWT (Time (Non-Lin)) were compared. Linear time-domain analysis used linear boundary conditions, and the computational nodes were assumed to be in the form of the first harmonic function. The linear and nonlinear analysis results were in good agreement with the frequency-domain results because the amplitude of the input wave was small and satisfied the linear condition.Figure 7
Comparison of the excitation forces for a truncated vertical cylinder on various numerical analysis methods (R = d = 1 m and A = 0.1 m under the deep-water condition).Figure8 compares the horizontal excitation force on a bottom-mounted circular cylinder for various wave steepness (H/λ=0.01and0.05). In this case, the water depth was set to 1 m, which is the same draft of the buoy. Using the second weakly nonlinear analysis, Bai and Teng [51] solved the same problem, and the results of the present simulation were in good agreement.Figure 8
Comparison of the horizontal excitation force for a bottom-mounted vertical cylinder on various wave steepness conditions (R = h = d = 1 m).


### 3.3. Radiation Problem

The radiation problem was conducted using the developed 3D-FN-PNWT technique. The radiation problem applies a forced harmonic motion to the buoy to measure the radiation force on the body. The added mass and radiation damping coefficients can be obtained from the radiation force (equations (35) and (36)).(38)maij=1πξ0ω∫0TFitsinωtdt+cijω2,(39)bij=−1πξ0∫0TFitcosωtdt,where ξ0 is the amplitude of the forced oscillating motion. More detailed information can be found in the report by Koo and Kim [12]. To solve this problem, the frontal, side, and end artificial damping zones were applied. Because there were no incident waves, both the frontal and end damping schemes were the same and the direction of the ramping function was different. Other specifications of the developed 3D-FN-PNWT were the same as those of the diffraction problem.Figure9 shows the time series of the displacement, velocity, and acceleration of a buoy and the radiated vertical force acting on the buoy in the radiation problem. In this case, the buoy was forcibly moved as a sine curve. To reduce the unstable numerical phenomenon, the amplitude of buoy displacement was gradually increased by the ramp function. The ramping duration was set to four times the buoy oscillating period. So, the transient motion components and horizontal force became steady-state results after five times the oscillating period. Using steady-state results, all other results of the radiation problem were compared with frequency-domain results.Figure 9
(a) Time series of the displacement, velocity, and acceleration of a buoy and (b) radiated vertical force acting on the buoy from the radiation problem (kR = 0.44 andξ0 = 0.1 m under the deep-water condition).Figure10 shows the added mass and radiation damping coefficients for the heave motion of a vertical cylinder buoy at various oscillation periods. The calculation results of the 3D-FN-PNWT were compared with the frequency-domain and linear time-domain results. The frequency-domain results were obtained using WAMIT, and the time-domain results were obtained using the linear 3D-PNWT. When the forced motion amplitude was small (H = 0.1 m) in the radiation problem, the nonlinear results agreed very well with both the linear frequency- and time-domain results because of small-amplitude radiated waves. However, it can be seen that the added mass and damping coefficients decrease slightly as the forced motion amplitude increases.Figure 10
Comparison of the added mass and radiation damping coefficients of a vertical cylinder buoy on various calculation methods (R = d = 1 m and h = 3 m).Figure11 compares the second-order heave radiation forces with those of other studies. The computational model had a radius, draft, and water depth of 1.0 m, 0.5 m, and 1.5 m, respectively. Teng [51] used the high-order boundary element method in the frequency domain. Bai [53] performed time-domain analysis using the B-spline-based boundary element method. In addition, Shao and Faltinsen [54] solved the same problem by applying weakly nonlinear hydrodynamic analysis. In the present 3D-FN-PNWT, the fast Fourier transform was used to extract the second-order wave component from the time series of the total wave forces. The present results agreed with Bai’s and Shao and Faltinsen’s results. The second-order wave component of the radiation force increases with increasing kR.Figure 11
Comparison of the 2nd order force components of a radiation force on various methods (R = 1 m, R/d = 2, h/R = 1.5, and ξ0/R = 0.05).


### 3.4. Heave-Only Allowed Problem

The analysis of the heave-only allowed body was performed using the 3D-FN-PNWT. All motion modes were restrained except for the heave motion. This case can be used to analyze the heaving buoy-type point absorber (HPA). The Stokes 2nd order wave was applied to the input boundary. The time-varying buoy force and buoy displacement were obtained using the acceleration potential approach and indirect method. The artificial damping zone was applied, as in the case of the diffraction problem.Figure12 shows the time history of the buoy displacement and velocity in two different initial conditions. Figure 12(a) presents results of the free decay test, and Figure 12(b) presents results of the velocity decay test. Both time results were compared with the results of the second-order expansion method by Bai and Teng [51]. Both results were in good agreement. Figure 13 compares the heave RAOs for various calculation analyses. The case of the 0.05 m incident wave amplitude, which is a relatively linear wave condition, was compared. The results of the frequency domain, linear time domain, and fully nonlinear time domain showed good agreement. Figure 14 shows the comparison of heave RAOs between the present analysis and the previous study of Kim et al. [55]. They conducted numerical analysis based on linear wave theory and experimental study for a vertical cylindrical heaving buoy. The results of the present study are closer to the experimental data than the linear ones because the nonlinear effects of the heaving buoy motion are taken into account.Figure 12
Time series of the buoy vertical displacement of a circular cylindrical buoy when initial displacement and velocity are given (A = 0.1 m, R = d = 1.0 m, and h = 3 d). (a) Free decay test. (b) Velocity decay test.
(a)
(b)Figure 13
Comparison of heave RAOs on various analysis methods (R = d = 1 m, h = 3 m, and A = 0.05 m).Figure 14
Comparison of heave RAOs on various analysis approaches (R = 0.2 m and d = 0.52 m under the deep-water condition).


### 3.5. Wave-Body-PTO Interaction

Using the 3D-FN-PNWT technique, a wave-body-PTO interaction was analyzed with hydraulic PTO modeling in the body system. To examine the effects of fully nonlinear analysis, linear and nonlinear analyses were carried out for small and large incident wave heights, respectively. The hydraulic PTO system was modeled, as shown in equation (33) [33]. The inner diameter of the hydraulic cylinder was fixed to 0.04 m, and the diameter of the cylinder rod was 0.018 m.Figure15 shows the time series of the buoy displacement, velocity, and PTO force acting on the buoy and the instantaneous generated power. In this case, the pressure difference of the hydraulic cylinder and hydraulic circuit was chosen as 12 bar. The whole simulation time was set to 30 times the wave period. As the buoy displacement occurs, the PTO force and wave power generate simultaneously. When the simulation time is longer than 13 times the wave period, the steady-state results occur on all time series.Figure 15
Time series of buoy displacement (a), velocity (b), PTO force (c), and generated power (d) (kR = 0.72,A = 0.3 m, and Δp = 12 bar).
(a)
(b)
(c)
(d)Figure16 compares the heave RAOs and time-averaged generated power on various wave steepness conditions. In this case, the optimal PTO force (FPTO/Fz = 0.4) was applied as given by Kim et al. [56]. In the high-kR condition (short period), the effect of wave steepness is very small. However, heave RAO and generated power decrease in the long-period condition (T/T0 is greater than 1). In particular, the effect of wave steepness is significant at the resonant wave period.Figure 16
Comparison of heave RAOs (a) and time-averaged generated power (b) on various wave steepness conditions (under the deep-water condition).
(a)
(b)Table1 shows high-order frequency components of the vertical displacement of the buoy for various wave steepness conditions. To separate zeroth-, first-, and second-order components, FFT (fast Fourier transform) was carried out. As wave steepness increases, the first component of the buoy displacement decreases and zeroth and second components increase. In particular, comparing the case of H/λ = 1/50 and H/λ = 1/20 in the resonance area (T/T0 = 1.0), the first-order component of the high-steepness case decreases 11% and mean and second-order components increase 62% and 114%, respectively. Therefore, as wave steepness increases, the generated power decreases by 10% because of the increase of the nonlinear wave effect.Table 1
Frequency components of the buoy displacement for various wave steepness conditions.
Z (nω)
H/λ
T/T0 = 0.8
T/T0 = 1.0
T/T0 = 1.4
1/20
1/30
1/40
1/50
1/20
1/30
1/40
1/50
1/20
1/30
1/40
1/50
0th
0.057
0.040
0.030
0.024
0.090
0.082
0.067
0.055
0.051
0.037
0.028
0.023
1st
0.236
0.247
0.251
0.252
1.371
1.487
1.525
1.542
0.826
0.849
0.859
0.863
2nd
0.002
0.002
0.001
0.001
0.065
0.048
0.037
0.030
0.104
0.073
0.056
0.045


## 3.1. Wave Generation and Propagation

Numerical simulations of wave generation and propagation were performed using the 3D-FN-PNWT technique. The dimensions of the 3D-FN-PNWT were six wavelengths in the longitudinal direction, one wavelength in the width direction, and 0.5 wavelengths in the vertical direction. The incident wave satisfied the deep-water condition because the water depth was set to 0.5 wavelengths in the vertical direction. The artificial damping zone, in which the length was two wavelengths, was installed near the end wall along the wave propagation direction.For input waves, the Stokes 2nd order waves can be expressed from the theoretical solution for an arbitrary depth as [50](37)ηt=Acos−ωt+kA3−tanh2kh4tanh3khcos−2ωt.Figure2 shows the time series of free surface elevations from numerical simulations and theoretical solutions. The input wave conditions were as follows: the wave period of three seconds, the wave height of 0.4 m, and the water depth of 3 m. The wave elevations from the 3D-FN-PNWT simulation agreed well with the theoretical Stokes 2nd order waves. Because the wave steepness (H/λ) was 0.0315 (≈1/32), the mean wave elevation was greater than zero because of wave nonlinearity.Figure 2
Comparison of the time series of wave elevations from the numerical and theoretical results (num: numerical results; theory: theoretical results).Figure3 shows the time series of the free surface elevations on various measured points and various wave steepness conditions. Wave elevations measuring all computational free surface points (four wavelengths) were the same without any reflections from the end wall. The amplitudes of the waves decreased at a point larger than four times the wavelength at which the artificial damping zone started. These results show that the artificial damping condition was well operated and the numerical wave tank satisfied the open-sea conditions. In Figure 3(b), when the wave steepness is small (H/λ=1/60), a linear wave with the same crest height and trough height is generated. On the contrary, as wave steepness increases, the crest wave height increases and is greater than the trough height. This is because the higher order term of the wave increases with increasing wave steepness.Figure 3
Time series of wave elevations on various measurement points (a) and various wave steepness conditions (b).
(a)
(b)


## 3.2. Diffraction Problem

To solve the diffraction problem of wave-body interactions, the Stokes 2nd order waves as the input condition were generated and propagated to the constrained body. The acceleration potential method was used to calculate the excitation force acting on the body. The dimensions of the computational domain were four wavelengths in the x direction, two wavelengths in the y direction, and 0.5 wavelengths in the z direction. The length of the artificial damping zone was set to one wavelength. Figure 4 shows meshes to present a cylindrical body and free surface around the body. The computational domain was modeled as a half model with the x-axis symmetry. All meshes in the whole calculational fluid domain were expressed using quadrilateral panels. The free surface meshes around the body were constructed radially from the body meshes. The body was a vertical cylinder with a draft and radius of 1 m each.Figure 4
Meshes to present a cylindrical body and free surface.Figure5 shows the time series of horizontal and vertical forces acting on a stationary vertical circular cylinder by solving the diffraction problem. The horizontal and vertical excitation forces were divided with ρgR2A, where R is the radius of the cylindrical buoy. As like a physical wave tank experiment, the hydrodynamic forces acting on the buoy occur after three times the wave period, which is the time when the incoming wave reaches the buoy. When the simulation time exceeds 8 times the wave period, the time series of both hydrodynamic forces are in a steady state. From this, the amplitude of horizontal and vertical forces can be obtained and compared with the results of the frequency-domain analysis. Figure 6 shows the convergence tests for the number of nodes per wavelength and interval of the time step. The vertical excitation forces on various kR conditions were compared. As the number of nodes increases and the interval of the time step decreases, the results of hydrodynamic analysis are converged. Therefore, optimal simulation conditions are adopted in this study (the number of nodes per wavelength is 20, and the interval of the time step is T/64).Figure 5
Time series of horizontal and vertical forces acting on a stationary vertical circular cylinder from the diffraction problem (kR = 0.44,A = 0.1 m, and R = d = 1.0 m under the deep-water condition).Figure 6
Convergence test for the number of nodes per wavelength (a) and time step (b) (R = d = 1.0 m under the deep-water condition).
(a)
(b)Figure7 compares the vertical and horizontal excitation forces acting on the vertical cylinder with the results of WAMIT based on the wave’s Green function. The draft and radius of the cylinder were 1 m each. The numerical results of the linear time-domain analysis (Time (Lin)) and fully nonlinear time-domain analysis using the 3D-FN-PNWT (Time (Non-Lin)) were compared. Linear time-domain analysis used linear boundary conditions, and the computational nodes were assumed to be in the form of the first harmonic function. The linear and nonlinear analysis results were in good agreement with the frequency-domain results because the amplitude of the input wave was small and satisfied the linear condition.Figure 7
Comparison of the excitation forces for a truncated vertical cylinder on various numerical analysis methods (R = d = 1 m and A = 0.1 m under the deep-water condition).Figure8 compares the horizontal excitation force on a bottom-mounted circular cylinder for various wave steepness (H/λ=0.01and0.05). In this case, the water depth was set to 1 m, which is the same draft of the buoy. Using the second weakly nonlinear analysis, Bai and Teng [51] solved the same problem, and the results of the present simulation were in good agreement.Figure 8
Comparison of the horizontal excitation force for a bottom-mounted vertical cylinder on various wave steepness conditions (R = h = d = 1 m).


## 3.3. Radiation Problem

The radiation problem was conducted using the developed 3D-FN-PNWT technique. The radiation problem applies a forced harmonic motion to the buoy to measure the radiation force on the body. The added mass and radiation damping coefficients can be obtained from the radiation force (equations (35) and (36)).(38)maij=1πξ0ω∫0TFitsinωtdt+cijω2,(39)bij=−1πξ0∫0TFitcosωtdt,where ξ0 is the amplitude of the forced oscillating motion. More detailed information can be found in the report by Koo and Kim [12]. To solve this problem, the frontal, side, and end artificial damping zones were applied. Because there were no incident waves, both the frontal and end damping schemes were the same and the direction of the ramping function was different. Other specifications of the developed 3D-FN-PNWT were the same as those of the diffraction problem.Figure9 shows the time series of the displacement, velocity, and acceleration of a buoy and the radiated vertical force acting on the buoy in the radiation problem. In this case, the buoy was forcibly moved as a sine curve. To reduce the unstable numerical phenomenon, the amplitude of buoy displacement was gradually increased by the ramp function. The ramping duration was set to four times the buoy oscillating period. So, the transient motion components and horizontal force became steady-state results after five times the oscillating period. Using steady-state results, all other results of the radiation problem were compared with frequency-domain results.Figure 9
(a) Time series of the displacement, velocity, and acceleration of a buoy and (b) radiated vertical force acting on the buoy from the radiation problem (kR = 0.44 andξ0 = 0.1 m under the deep-water condition).Figure10 shows the added mass and radiation damping coefficients for the heave motion of a vertical cylinder buoy at various oscillation periods. The calculation results of the 3D-FN-PNWT were compared with the frequency-domain and linear time-domain results. The frequency-domain results were obtained using WAMIT, and the time-domain results were obtained using the linear 3D-PNWT. When the forced motion amplitude was small (H = 0.1 m) in the radiation problem, the nonlinear results agreed very well with both the linear frequency- and time-domain results because of small-amplitude radiated waves. However, it can be seen that the added mass and damping coefficients decrease slightly as the forced motion amplitude increases.Figure 10
Comparison of the added mass and radiation damping coefficients of a vertical cylinder buoy on various calculation methods (R = d = 1 m and h = 3 m).Figure11 compares the second-order heave radiation forces with those of other studies. The computational model had a radius, draft, and water depth of 1.0 m, 0.5 m, and 1.5 m, respectively. Teng [51] used the high-order boundary element method in the frequency domain. Bai [53] performed time-domain analysis using the B-spline-based boundary element method. In addition, Shao and Faltinsen [54] solved the same problem by applying weakly nonlinear hydrodynamic analysis. In the present 3D-FN-PNWT, the fast Fourier transform was used to extract the second-order wave component from the time series of the total wave forces. The present results agreed with Bai’s and Shao and Faltinsen’s results. The second-order wave component of the radiation force increases with increasing kR.Figure 11
Comparison of the 2nd order force components of a radiation force on various methods (R = 1 m, R/d = 2, h/R = 1.5, and ξ0/R = 0.05).


## 3.4. Heave-Only Allowed Problem

The analysis of the heave-only allowed body was performed using the 3D-FN-PNWT. All motion modes were restrained except for the heave motion. This case can be used to analyze the heaving buoy-type point absorber (HPA). The Stokes 2nd order wave was applied to the input boundary. The time-varying buoy force and buoy displacement were obtained using the acceleration potential approach and indirect method. The artificial damping zone was applied, as in the case of the diffraction problem.Figure12 shows the time history of the buoy displacement and velocity in two different initial conditions. Figure 12(a) presents results of the free decay test, and Figure 12(b) presents results of the velocity decay test. Both time results were compared with the results of the second-order expansion method by Bai and Teng [51]. Both results were in good agreement. Figure 13 compares the heave RAOs for various calculation analyses. The case of the 0.05 m incident wave amplitude, which is a relatively linear wave condition, was compared. The results of the frequency domain, linear time domain, and fully nonlinear time domain showed good agreement. Figure 14 shows the comparison of heave RAOs between the present analysis and the previous study of Kim et al. [55]. They conducted numerical analysis based on linear wave theory and experimental study for a vertical cylindrical heaving buoy. The results of the present study are closer to the experimental data than the linear ones because the nonlinear effects of the heaving buoy motion are taken into account.Figure 12
Time series of the buoy vertical displacement of a circular cylindrical buoy when initial displacement and velocity are given (A = 0.1 m, R = d = 1.0 m, and h = 3 d). (a) Free decay test. (b) Velocity decay test.
(a)
(b)Figure 13
Comparison of heave RAOs on various analysis methods (R = d = 1 m, h = 3 m, and A = 0.05 m).Figure 14
Comparison of heave RAOs on various analysis approaches (R = 0.2 m and d = 0.52 m under the deep-water condition).


## 3.5. Wave-Body-PTO Interaction

Using the 3D-FN-PNWT technique, a wave-body-PTO interaction was analyzed with hydraulic PTO modeling in the body system. To examine the effects of fully nonlinear analysis, linear and nonlinear analyses were carried out for small and large incident wave heights, respectively. The hydraulic PTO system was modeled, as shown in equation (33) [33]. The inner diameter of the hydraulic cylinder was fixed to 0.04 m, and the diameter of the cylinder rod was 0.018 m.Figure15 shows the time series of the buoy displacement, velocity, and PTO force acting on the buoy and the instantaneous generated power. In this case, the pressure difference of the hydraulic cylinder and hydraulic circuit was chosen as 12 bar. The whole simulation time was set to 30 times the wave period. As the buoy displacement occurs, the PTO force and wave power generate simultaneously. When the simulation time is longer than 13 times the wave period, the steady-state results occur on all time series.Figure 15
Time series of buoy displacement (a), velocity (b), PTO force (c), and generated power (d) (kR = 0.72,A = 0.3 m, and Δp = 12 bar).
(a)
(b)
(c)
(d)Figure16 compares the heave RAOs and time-averaged generated power on various wave steepness conditions. In this case, the optimal PTO force (FPTO/Fz = 0.4) was applied as given by Kim et al. [56]. In the high-kR condition (short period), the effect of wave steepness is very small. However, heave RAO and generated power decrease in the long-period condition (T/T0 is greater than 1). In particular, the effect of wave steepness is significant at the resonant wave period.Figure 16
Comparison of heave RAOs (a) and time-averaged generated power (b) on various wave steepness conditions (under the deep-water condition).
(a)
(b)Table1 shows high-order frequency components of the vertical displacement of the buoy for various wave steepness conditions. To separate zeroth-, first-, and second-order components, FFT (fast Fourier transform) was carried out. As wave steepness increases, the first component of the buoy displacement decreases and zeroth and second components increase. In particular, comparing the case of H/λ = 1/50 and H/λ = 1/20 in the resonance area (T/T0 = 1.0), the first-order component of the high-steepness case decreases 11% and mean and second-order components increase 62% and 114%, respectively. Therefore, as wave steepness increases, the generated power decreases by 10% because of the increase of the nonlinear wave effect.Table 1
Frequency components of the buoy displacement for various wave steepness conditions.
Z (nω)
H/λ
T/T0 = 0.8
T/T0 = 1.0
T/T0 = 1.4
1/20
1/30
1/40
1/50
1/20
1/30
1/40
1/50
1/20
1/30
1/40
1/50
0th
0.057
0.040
0.030
0.024
0.090
0.082
0.067
0.055
0.051
0.037
0.028
0.023
1st
0.236
0.247
0.251
0.252
1.371
1.487
1.525
1.542
0.826
0.849
0.859
0.863
2nd
0.002
0.002
0.001
0.001
0.065
0.048
0.037
0.030
0.104
0.073
0.056
0.045


## 4. Conclusion

Various hydrodynamic analyses of the vertical cylinder were carried out using a three-dimensional fully nonlinear potential numerical wave tank (3D-FN-PNWT). The developed 3D-FN-PNWT was based on the boundary element method with constant panels. The least-square gradient reconstruction method was applied to calculate the spatial gradient of the boundary elements. Time-varying instantaneous forces on the body with nonlinear waves were calculated using the MEL method and the acceleration potential approach. To realize open-sea conditions, artificial damping zones, such as the front, end, and side areas of the domain, were imposed to avoid the reflection waves. For verification, wave generation and propagation and wave diffraction and radiation problems were solved and compared with the results of previous studies. In particular, the second-order diffraction wave components were compared with those of other studies.Numerical modeling of the wave-body-PTO interaction was performed using the developed 3D-FN-PNWT technique. The heave RAOs of linear analysis and fully nonlinear analysis were compared for small- and high-amplitude wave conditions. Both the linear and fully nonlinear results were in good agreement in the case of a small-amplitude wave because of the linear assumed input. In contrast to the case of a small-amplitude wave, large differences between the linear and fully nonlinear results occurred, particularly at the resonance period. When the optimal PTO condition was applied, the heave RAO and time-averaged generated power were reduced by approximately 10% because of the nonlinear wave effect. Therefore, the hydrodynamic performance of wave-body-PTO simulation in the linear analysis is overestimated compared to that in fully nonlinear analysis, which may be important for the design of a wave energy converter.



---

*Source: 5163597-2019-09-19.xml*