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

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 2^nd order diffraction approach and fully nonlinear analysis for a fixed cylindrical body. They confirmed that fully nonlinear analysis results were more accurate than the 2^nd 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 2^nd 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 5^th 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 2^nd identity, the governing equation can be transformed into the boundary integral equation as follows:$$\begin{array}{}\text{(1)}& \alpha \varphi ={\int }_S\left(G\frac{\partial \varphi }{\partial n}-\varphi \frac{\partial G}{\partial n}\right)dS,\end{array}$$where $\alpha$ is the solid angle and $G$ is the Green function, which can be expressed as$$\begin{array}{}\text{(2)}& G=\frac{1}{r},\end{array}$$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:$$\begin{array}{c}\text{(3)}& {\iint }_{S}GdS=\sum _{m=0}^{\infty }{\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^{m+n}}{m!n!}{I}_{mn}\frac{{\partial }^{m+n}}{\partial x^m\partial y^n}\left(x^2+y^2+z^2\right)}^{-1/2},{\iint }_S\frac{\partial G}{\partial n}dS=\sum _{m=0}^{\infty }{\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^{m+n+1}}{m!n!}{I}_{mn}z\frac{{\partial }^{m+n}}{\partial x^m\partial y^n}\left(x^2+y^2+z^2\right)}^{-3/2},\end{array}$$where $I_{mn}$ 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.

Figure 1 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 2^nd order wave potential profile was applied to the incident wave boundary condition, and the angle of incident waves was set to zero:$$\begin{array}{c}\text{(4)}& \frac{\partial \varphi }{\partial n}=n_x\frac{\partial \varphi }{\partial x}=n_x\frac{gAk}{\omega }\frac{\cosh k\left(z+h\right)}{\cosh kh}\cos \left(kx-\omega t\right)\left(\text{linear}\right)+n_x\frac{3}{8}{A}^2\omega \left(2k\right)\frac{\cosh 2k\left(z+h\right)}{{\sinh }^{4}kh}\cos 2\left(kx-\omega t\right)\left(\text{Stokes\hspace{0.17em}}2\text{nd}\right),\\ \text{on\hspace{0.17em}}{S}_{\text{in}},\end{array}$$where $g$, $A$, $k$, $\omega$, and $h$ are the gravitational acceleration, wave amplitude, wave number, wave frequency, and water depth, respectively. $\overrightarrow{n}=n_x\overrightarrow{i}+n_y\overrightarrow{j}+n_z\overrightarrow{k}$ 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]:$$\begin{array}{}\text{(5)}& \frac{\delta \eta }{\delta t}=\frac{\partial \varphi }{\partial z}-\left(\nabla \varphi -\overrightarrow{v}\right)\cdot \nabla \eta ,\hspace{1em}\text{on\hspace{0.17em}}{S}_F,\text{(6)}& \frac{\delta \varphi }{\delta t}=-g\eta -\frac{1}{2}{\left|\nabla \varphi \right|}^2+\overrightarrow{v}\cdot \nabla \varphi ,\hspace{1em}\text{on\hspace{0.17em}}{S}_F,\end{array}$$where $\overrightarrow{v}$ denotes the velocity of the calculation node on the surface and $\delta /\delta t$ and $\partial /\partial t$ denote the material derivative ($=\left(\partial /\partial t\right)+u\cdot \nabla$) 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 ($\overrightarrow{v}=\left(\dot{\xi },\dot{\zeta },\dot{\eta }\right)$). In the semi-Lagrangian method, the node velocity moves only in the vertical direction of the water particles on the free surface ($\overrightarrow{v}=\left(\mathrm{0,0},\dot{\eta }\right)$). 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$$\begin{array}{}\text{(7)}& \frac{\delta \eta }{\delta t}=\frac{\partial \varphi }{\partial z}-\frac{\partial \varphi }{\partial x}\frac{\partial \eta }{\partial x}-\frac{\partial \varphi }{\partial y}\frac{\partial \eta }{\partial y},\hspace{1em}\text{on\hspace{0.17em}}{S}_F,\text{(8)}& \frac{\delta \varphi }{\delta t}=-g\eta -\frac{1}{2}{\left|\nabla \varphi \right|}^2+\frac{\delta \eta }{\delta t}\frac{\partial \varphi }{\partial z},\hspace{1em}\text{on\hspace{0.17em}}{S}_F.\end{array}$$

[figure omitted; refer to PDF]

The impervious boundary condition given in the following condition was applied to the side and back walls:$$\begin{array}{}\text{(9)}& \frac{\partial \varphi }{\partial n}=0,\hspace{1em}\text{on\hspace{0.17em}}{S}_W.\end{array}$$

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:$$\begin{array}{}\text{(10)}& \frac{\partial \varphi }{\partial n}=\overrightarrow{V}\cdot \overrightarrow{n},\hspace{1em}\text{on\hspace{0.17em}}{S}_B,\end{array}$$where $\overrightarrow{V}$ and $\overrightarrow{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 ${\varphi }_n-\eta$-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]:$$\begin{array}{}\text{(11)}& \frac{\delta \eta }{\delta t}=\frac{\partial \varphi }{\partial z}-\frac{\partial \varphi }{\partial x}\frac{\partial \eta }{\partial x}-\frac{\partial \varphi }{\partial y}\frac{\partial \eta }{\partial y}-{\mu }_{E2}\eta ,\hspace{1em}\text{on\hspace{0.17em}}{S}_{FE},\text{(12)}& \frac{\delta \varphi }{\delta t}=-g\eta -\frac{1}{2}{\left|\nabla \varphi \right|}^2+\frac{\delta \eta }{\delta t}\frac{\partial \varphi }{\partial z}-{\mu }_{E1}\frac{\partial \varphi }{\partial n}\text{,}\hspace{1em}\text{on\hspace{0.17em}}{S}_{FE},\text{(13)}& {\mu }_{E1}\left(x\right)={\mu }_{E10}\left(1-\cos \left(\frac{\pi }{2}\frac{x-x_{\text{ed}}}{l_{\text{d}}}\right)\right),\end{array}$$where ${\mu }_{E1}$ and ${\mu }_{E2}$ are the end artificial damping coefficients (${\mu }_{E2}=k{\mu }_{E1}$) and $S_{FE}$, ${\mu }_{E10}$, $l_{\text{d}}$, and $x_{\text{ed}}$ 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 (${\varphi }_n^{\ast }$ and ${\eta }^{\ast }$) calculated by the theoretical solution of the incident wave. The following equations show the boundary conditions for the frontal artificial damping zone [18]:$$\begin{array}{}\text{(14)}& \frac{\delta \eta }{\delta t}=\frac{\partial \varphi }{\partial z}-\frac{\partial \varphi }{\partial x}\frac{\partial \eta }{\partial x}-\frac{\partial \varphi }{\partial y}\frac{\partial \eta }{\partial y}-{\mu }_{F2}\left(\eta -{\eta }^{\ast }\right),\hspace{1em}\text{on\hspace{0.17em}}{S}_{FF},\text{(15)}& \frac{\delta \varphi }{\delta t}=-g\eta -\frac{1}{2}{\left|\nabla \varphi \right|}^2+\frac{\delta \eta }{\delta t}\frac{\partial \varphi }{\partial z}-{\mu }_{F1}\left(\frac{\partial \varphi }{\partial n}-{\frac{\partial \varphi }{\partial n}}^{\ast }\right),\hspace{1em}\text{on\hspace{0.17em}}{S}_{FF},\end{array}$$where ${\mu }_{F1}$ and ${\mu }_{F2}$ denote the frontal damping coefficients and $S_{FF}$ 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:$$\begin{array}{}\text{(16)}& {\mu }_{F1}\left(x\right)={\mu }_{F10}\left(1-\cos \left(\frac{\pi }{2}\frac{x_{\text{fd}}-x}{l_{\text{d}}}\right)\right),\end{array}$$where ${\mu }_{F10}$ and $x_{\text{fd}}$ 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:$$\begin{array}{}\text{(17)}& \frac{\delta \varphi }{\delta t}=-g\eta -\frac{1}{2}{\left|\nabla \varphi \right|}^2+\overrightarrow{v}\cdot \nabla \varphi -{\kappa }_1{\mu }_{S1}\left\{\frac{\partial \varphi }{\partial n}-{\left(\frac{\partial \varphi }{\partial n}\right)}^{\ast }\right\},\text{(18)}& \frac{\delta \eta }{\delta t}=\frac{\partial \varphi }{\partial z}-\left(\nabla \varphi -\overrightarrow{v}\right)\cdot \nabla \eta -{\kappa }_1{\mu }_{S2}\left(\eta -{\eta }^{\ast }\right),\text{(19)}& T_{S0}=\frac{x}{c_g},\text{(20)}& {\kappa }_1=\left\{\begin{array}{cc}0,\hfill & \text{for\hspace{0.17em}}t<T_{S0},\hfill \\ 1-\cos \frac{\pi \left(t-T_{S0}\right)}{2T_{\text{Ramp}}},\hfill & \text{for\hspace{0.17em}}{T}_{S0}\le t,\hfill \\ 1,\hfill & \text{for\hspace{0.17em}}{T}_{S1}<t,\hfill \end{array}\right.\end{array}$$where ${\mu }_{S1}$ and ${\mu }_{S2}$ are the side artificial damping coefficients and ${\kappa }_1$, $T_{\text{Ramp}}$, $T_{S0}$, and $T_{S1}$ 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$$\begin{array}{}\text{(21)}& {\mu }_{S1}\left(x\right)={\mu }_{S10}\left(1-\cos \left(\frac{\pi }{2}\frac{y-y_{\text{sd}}}{l_{\text{d}}}\right)\right),\end{array}$$where ${\mu }_{S10}$ and $y_{\text{sd}}$ 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 “$Q_i$” at the neighboring calculation point “$i$” based on a specific calculation point “$0$” using a Taylor series transformation [43]:$$\begin{array}{}\text{(22)}& Q_i=Q_0+{\left.\frac{\partial Q}{\partial x}\right|}_0\Delta x_i+{\left.\frac{\partial Q}{\partial y}\right|}_0\Delta y_i+{\left.\frac{\partial Q}{\partial z}\right|}_0\Delta z_i+{\left.\frac{{\partial }^{2}Q}{\partial x^2}\right|}_0\frac{1}{2}\Delta x_i^2+{\left.\frac{{\partial }^{2}Q}{\partial y^2}\right|}_0\frac{1}{2}\Delta y_i^2+{\left.\frac{{\partial }^{2}Q}{\partial x^2}\right|}_0\frac{1}{2}\Delta y_i^2+{\left.\frac{{\partial }^{2}Q}{\partial x\partial y}\right|}_0\Delta x\Delta y+{\left.\frac{{\partial }^{2}Q}{\partial y\partial z}\right|}_0\Delta y\Delta z+{\left.\frac{{\partial }^{2}Q}{\partial x\partial z}\right|}_0\Delta x\Delta z,\end{array}$$where $\Delta x_i$, $\Delta y_i$, and $\Delta z_i$ 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$$\begin{array}{}\text{(23)}& f\left(w_i,{\left.\frac{\partial Q}{\partial x}\right|}_0,{\left.\frac{\partial Q}{\partial y}\right|}_0,{\left.\frac{\partial Q}{\partial z}\right|}_0,{\left.\frac{{\partial }^{2}Q}{\partial x^2}\right|}_0,{\left.\frac{{\partial }^{2}Q}{\partial y^2}\right|}_0,{\left.\frac{{\partial }^{2}Q}{\partial z^2}\right|}_0,{\left.\frac{{\partial }^{2}Q}{\partial x\partial y}\right|}_0,{\left.\frac{{\partial }^{2}Q}{\partial y\partial z}\right|}_0,{\left.\frac{{\partial }^{2}Q}{\partial x\partial z}\right|}_0\right)={\sum _{i=1}^N{w}_i\left[\Delta Q_i-\left({\left.\frac{\partial Q}{\partial x}\right|}_0\Delta x_i+{\left.\frac{\partial Q}{\partial y}\right|}_0\Delta y_i+{\left.\frac{\partial Q}{\partial z}\right|}_0\Delta z_i+{\left.\frac{{\partial }^{2}Q}{\partial x^2}\right|}_0\frac{1}{2}\Delta x_i^2+{\left.\frac{{\partial }^{2}Q}{\partial y^2}\right|}_0\frac{1}{2}\Delta y_i^2+{\left.\frac{{\partial }^{2}Q}{\partial z^2}\right|}_0\frac{1}{2}\Delta z_i^2+{\left.\frac{{\partial }^{2}Q}{\partial x\partial y}\right|}_0\Delta x\Delta y+{\left.\frac{{\partial }^{2}Q}{\partial y\partial z}\right|}_0\Delta y\Delta z+{\left.\frac{{\partial }^{2}Q}{\partial x\partial z}\right|}_0\Delta x\Delta z\right)\right]}^2,\end{array}$$where $w_i$ denotes the weight for the specific calculation point “$0$” of the i^th neighboring calculation point, and it can be expressed by the reciprocal of the fifth power of the distance from the i^th 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 (${\varphi }_t=\partial \varphi /\partial 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]:$$\begin{array}{}\text{(24)}& \overrightarrow{a}=\frac{\delta \nabla \varphi }{\delta t}=\frac{\partial \nabla \varphi }{\partial t}+\left(\nabla \varphi \cdot \nabla \right)\nabla \varphi =\nabla \left(\frac{\partial \varphi }{\partial t}+\frac{1}{2}{\left|\nabla \varphi \right|}^2\right)=\nabla \Phi ,\end{array}$$where $\Phi$ 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:$$\begin{array}{}\text{(25)}& {\varphi }_t=\frac{\partial \varphi }{\partial t}=\Phi -\frac{1}{2}{\left|\nabla \varphi \right|}^2.\end{array}$$

The body boundary condition is difficult to describe in the acceleration potential field because of the nonlinear term ${\varphi }_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:$$\begin{array}{}\text{(26)}& \frac{\partial {\varphi }_t}{\partial n}=\overrightarrow{n}\cdot \left({\overrightarrow{a}}_0+\overrightarrow{\dot{\omega }}\times \overrightarrow{r}\right)+q_B,\text{(27)}& q_B=\left(\overrightarrow{\omega }\cdot {\overrightarrow{s}}_1\right)\left(\frac{\partial \varphi }{\partial s_2}-2\left(\overrightarrow{v}\cdot {\overrightarrow{s}}_2\right)\right)-\left(\overrightarrow{\omega }\cdot {\overrightarrow{s}}_2\right)\left(\frac{\partial \varphi }{\partial s_1}-2\left(\overrightarrow{v}\cdot {\overrightarrow{s}}_1\right)\right)+k_{n1}\left(\overrightarrow{v}\cdot {\overrightarrow{s}}_1\right)\left(\frac{\partial \varphi }{\partial s_1}-\left(\overrightarrow{v}\cdot {\overrightarrow{s}}_1\right)\right)+k_{n2}\left(\overrightarrow{v}\cdot {\overrightarrow{s}}_2\right)\left(\frac{\partial \varphi }{\partial s_2}-\left(\overrightarrow{v}\cdot {\overrightarrow{s}}_2\right)\right)+\left(\overrightarrow{v}\cdot \overrightarrow{n}\right)\left(\frac{{\partial }^2\varphi }{\partial s_1^2}+\frac{{\partial }^2\varphi }{\partial s_2^2}+\left(k_{n1}+k_{n2}\right)\frac{\partial \varphi }{\partial n}\right),\end{array}$$where $a_0$, $\overrightarrow{\dot{\omega }}$, and $\overrightarrow{r}$ denote the translate acceleration of the buoy, angular acceleration of the buoy, and local coordinate vector on the body boundary, ${\overrightarrow{s}}_1$ and ${\overrightarrow{s}}_2$ are the tangent vectors of the quadratic panel, and $k_{n1}$ and $k_{n2}$ 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$$\begin{array}{}\text{(28)}& \frac{\partial \varphi }{\partial t}=\sum _{i=1}^6{a}_i{\phi }_i+{\phi }_7,\end{array}$$where $a_i$ and ${\phi }_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. ${\phi }_7$ is the acceleration potential component from the velocity field. The indirect method is a method of calculating ${\phi }_7$ by applying the virtual velocity potential (${\psi }_i$) satisfying the governing equation (equation (29)) with the Green 2^nd identity (equation (30)):$$\begin{array}{}\text{(29)}& {\nabla }^2{\psi }_i=0,\text{(30)}& {\iint }_S\left({\phi }_7\frac{\partial {\psi }_i}{\partial n}-{\psi }_i\frac{\partial {\phi }_7}{\partial n}\right)ds=0.\end{array}$$

Substituting the boundary condition into equation (30) yields$$\begin{array}{}\text{(31)}& {\iint }_{S_B}{\phi }_7{n}_{i}ds={\iint }_{S_B}{\psi }_i\left(q_B\right)ds-{\iint }_{S_F}\frac{\partial {\psi }_i}{\partial n}\left(-g\eta -\frac{1}{2}{\left|\nabla \varphi \right|}^2\right)ds+{\iint }_{S_{\text{in}}}{\psi }_i\left(\frac{\partial {\varphi }_t}{\partial n}\right)ds.\end{array}$$

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:$$\begin{array}{}\text{(32)}& P=-\rho \left(\sum _{i=1}^6{a}_i{\phi }_i+{\phi }_7+gz+\frac{1}{2}{\left|\nabla \varphi \right|}^2\right),\text{(33)}& m{a}_i=F_i={\int }_{S_B}P{n}_{i}ds-W{\delta }_{i3}-F_{\text{PTO}},\end{array}$$where $W$ and $F_{\text{PTO}}$ 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].$$\begin{array}{}\text{(34)}& Z_{\text{p}}=\frac{{\sum }_{i=1}^n{Z}_i{w}_i}{{\sum }_{i=1}^n{w}_i},\end{array}$$where $Z_i$ and $w_i$ are the known quantity of the i^th known calculation point and the weighting of the i^th known point for the target point and $Z_{\text{p}}$ is the unknown quantity for the target point. The weighting was calculated by the reciprocal of the fourth power of the distance from the i^th 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.$$\begin{array}{}\text{(35)}& r\left(t\right)=\left\{\begin{array}{cc}1,& \text{for\hspace{0.17em}}t>T_{\text{Ramp}},\\ 1-\cos \left(\frac{\pi }{2}\frac{t}{T_{\text{Ramp}}}\right),& \text{for\hspace{0.17em}}t\le T_{\text{Ramp}}.\end{array}\right.\end{array}$$

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]).$$\begin{array}{}\text{(36)}& F_{\text{PTO}}=-\min \left(G{\dot{z}}_i,\text{sign}\left({\dot{z}}_i\right)S_{\text{c}}\Delta p\right),\end{array}$$where $S_{\text{c}}$ and $\Delta 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 2^nd order waves can be expressed from the theoretical solution for an arbitrary depth as [50]$$\begin{array}{}\text{(37)}& \eta \left(t\right)=A\left\{\cos \left(-\omega t\right)+kA\frac{3-{\tanh }^2\left(kh\right)}{4{\tanh }^3\left(kh\right)}\cos \left(-2\omega t\right)\right\}.\end{array}$$

Figure 2 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 2^nd order waves. Because the wave steepness ($H/\lambda$) was 0.0315 (≈1/32), the mean wave elevation was greater than zero because of wave nonlinearity.

[figure omitted; refer to PDF]

Figure 3 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/\lambda =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.

[figures omitted; refer to PDF]

3.2. Diffraction Problem

To solve the diffraction problem of wave-body interactions, the Stokes 2^nd 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 omitted; refer to PDF]

Figure 5 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 $\rho g{R}^{2}A$, 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 omitted; refer to PDF]

[figures omitted; refer to PDF]

Figure 7 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 omitted; refer to PDF]

Figure 8 compares the horizontal excitation force on a bottom-mounted circular cylinder for various wave steepness ($H/\lambda =0.01\text{and\hspace{0.17em}}0.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 omitted; refer to PDF]

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)).$$\begin{array}{}\text{(38)}& m_{aij}=\frac{1}{\pi {\xi }_0\omega }{\int }_0^T{F}_i\left(t\right)\sin \omega tdt+\frac{c_{ij}}{{\omega }^2},\text{(39)}& b_{ij}=-\frac{1}{\pi {\xi }_0}{\int }_0^T{F}_i\left(t\right)\cos \omega tdt,\end{array}$$where ${\xi }_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.

Figure 9 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 omitted; refer to PDF]

Figure 10 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 omitted; refer to PDF]

Figure 11 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 omitted; refer to PDF]

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 2^nd 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.

Figure 12 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.

[figures omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

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.

Figure 15 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.

[figures omitted; refer to PDF]

Figure 16 compares the heave RAOs and time-averaged generated power on various wave steepness conditions. In this case, the optimal PTO force (F_PTO/F_z = 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/T₀ is greater than 1). In particular, the effect of wave steepness is significant at the resonant wave period.

[figures omitted; refer to PDF]

Table 1 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/T₀ = 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.

Acknowledgments

This work was supported by Inha University Research Grant.

Appendix

Appendix

Least-Square Gradient Reconstruction Method

$$\begin{array}{c}\text{(A.1)}& \frac{\partial f}{\partial X_{0,x}}=\sum _{i=1}^{N}2w_i\Delta x_{i}F,\frac{\partial f}{\partial X_{0,y}}=\sum _{i=1}^{N}2w_i\Delta y_{i}F,\\ \frac{\partial f}{\partial X_{0,z}}=\sum _{i=1}^{N}2w_i\Delta z_{i}F,\\ \frac{\partial f}{\partial X_{0,xx}}=\sum _{i=1}^{N}2w_i\frac{1}{2}\Delta x_i^{2}F,\\ \frac{\partial f}{\partial X_{0,yy}}=\sum _{i=1}^{N}2w_i\frac{1}{2}\Delta y_i^{2}F,\\ \frac{\partial f}{\partial X_{0,zz}}=\sum _{i=1}^{N}2w_i\frac{1}{2}\Delta z_i^{2}F,\\ \frac{\partial f}{\partial X_{0,xy}}=\sum _{i=1}^{N}2w_i\Delta x_i\Delta y_{i}F,\\ \frac{\partial f}{\partial X_{0,yz}}=\sum _{i=1}^{N}2w_i\Delta y_i\Delta z_{i}F,\\ \frac{\partial f}{\partial X_{0,xz}}=\sum _{i=1}^{N}2w_i\Delta x_i\Delta z_{i}F,\end{array}$$ where the function “$F$” can be expressed as$$\begin{array}{}\text{(A.2)}& F=\Delta Q_i-\left({\left.\frac{\partial Q}{\partial x}\right|}_0\Delta x_i+{\left.\frac{\partial Q}{\partial y}\right|}_0\Delta y_i+{\left.\frac{\partial Q}{\partial z}\right|}_0\Delta z_i+{\left.\frac{{\partial }^{2}Q}{\partial x^2}\right|}_0\frac{1}{2}\Delta x_i^2+{\left.\frac{{\partial }^{2}Q}{\partial y^2}\right|}_0\frac{1}{2}\Delta y_i^2+{\left.\frac{{\partial }^{2}Q}{\partial z^2}\right|}_0\frac{1}{2}\Delta z_i^2+{\left.\frac{{\partial }^{2}Q}{\partial x\partial y}\right|}_0\Delta x\Delta y+{\left.\frac{{\partial }^{2}Q}{\partial y\partial z}\right|}_0\Delta y\Delta z+{\left.\frac{{\partial }^{2}Q}{\partial x\partial z}\right|}_0\Delta x\Delta z\right).\end{array}$$

Using this, a 9 × 9 linear algebraic system can be transformed into a matrix of Ax = b. The size of the A matrix is 9 × 9, and the size of the x matrix and b matrix is 9 × 1. The following equation shows the A matrix, x matrix, and b matrix:$$\begin{array}{c}\text{(A.3)}& A=\left(\begin{array}{ccccccccc}{A}_{11}& A_{12}& A_{13}& A_{14}& A_{15}& A_{16}& A_{17}& A_{18}& A_{19}\\ A_{21}& A_{22}& A_{23}& A_{24}& A_{25}& A_{26}& A_{27}& A_{28}& A_{29}\\ A_{31}& A_{32}& A_{33}& A_{34}& A_{35}& A_{36}& A_{37}& A_{38}& A_{39}\\ A_{41}& A_{42}& A_{43}& A_{44}& A_{45}& A_{46}& A_{47}& A_{48}& A_{49}\\ A_{51}& A_{52}& A_{53}& A_{54}& A_{55}& A_{56}& A_{57}& A_{58}& A_{59}\\ A_{61}& A_{62}& A_{63}& A_{64}& A_{65}& A_{66}& A_{67}& A_{68}& A_{69}\\ A_{71}& A_{72}& A_{73}& A_{74}& A_{75}& A_{76}& A_{77}& A_{78}& A_{79}\\ A_{81}& A_{82}& A_{83}& A_{84}& A_{85}& A_{86}& A_{87}& A_{88}& A_{89}\\ A_{91}& A_{92}& A_{93}& A_{94}& A_{95}& A_{96}& A_{97}& A_{98}& A_{99}\end{array}\right),x=\left(\begin{array}{c}{x}_{11}\\ x_{21}\\ x_{31}\\ x_{41}\\ x_{51}\\ x_{61}\\ x_{71}\\ x_{81}\\ x_{91}\end{array}\right),\\ b=\left(\begin{array}{c}{b}_{11}\\ b_{21}\\ b_{31}\\ b_{41}\\ b_{51}\\ b_{61}\\ b_{71}\\ b_{81}\\ b_{91}\end{array}\right),\end{array}$$where matrix A is a symmetric matrix. Thus, equation (A.4) is satisfied, and equation (A.5) means the value of each element of the A matrix:$$\begin{array}{c}\text{(A.4)}& A_{ij}=A_{ji},\text{(A.5)}& A_{11}=w_i\Delta x_i^2,A_{12}=w_i\Delta x_i\Delta y_i,\\ A_{13}=w_i\Delta x_i\Delta z_i,\\ A_{22}=w_i\Delta y_i^2,\\ A_{23}=w_i\Delta y_i\Delta z_i,\\ A_{33}=w_i\Delta z_i^2,\\ A_{14}=\frac{1}{2}{w}_i\Delta x_i^3,\\ A_{15}=\frac{1}{2}{w}_i\Delta x_i\Delta y_i^2,\\ A_{16}=\frac{1}{2}{w}_i\Delta x_i\Delta z_i^2,\\ A_{24}=\frac{1}{2}{w}_i\Delta x_i^2\Delta y_i,\\ A_{25}=\frac{1}{2}{w}_i\Delta y_i^2\Delta z_i,\\ A_{26}=\frac{1}{2}{w}_i\Delta y_i\Delta z_i^2,\\ A_{34}=\frac{1}{2}{w}_i\Delta x_i^2\Delta z_i,\\ A_{35}=\frac{1}{2}{w}_i\Delta x_i\Delta y_i^2,\\ A_{36}=\frac{1}{2}{w}_i\Delta z_i^3,\\ A_{17}=w_i\Delta x_i^2\Delta y_i,\\ A_{18}=w_i\Delta x_i\Delta y_i\Delta z_i,\\ A_{19}=w_i\Delta x_i^2\Delta z_i,\\ A_{27}=w_i\Delta x_i\Delta y_i^2,\\ A_{28}=w_i\Delta y_i^2\Delta z_i,\\ A_{29}=w_i\Delta x_i\Delta y_i\Delta z_i,\\ A_{37}=w_i\Delta x_i\Delta y_i\Delta z_i,\\ A_{38}=w_i\Delta y_i\Delta z_i^2,\\ A_{39}=w_i\Delta x_i\Delta z_i^2,\\ A_{44}=\frac{1}{4}{w}_i\Delta x_i^4,\\ A_{45}=\frac{1}{4}{w}_i\Delta x_i^2\Delta y_i^2,\\ A_{46}=\frac{1}{4}{w}_i\Delta x_i^2\Delta z_i^2,\\ A_{55}=\frac{1}{4}{w}_i\Delta y_i^4,\\ A_{56}=\frac{1}{4}{w}_i\Delta y_i^2\Delta z_i^2,\\ A_{66}=\frac{1}{4}{w}_i\Delta z_i^4,\\ A_{47}=\frac{1}{2}{w}_i\Delta x_i^3\Delta y_i,\\ A_{48}=\frac{1}{2}{w}_i\Delta x_i^2\Delta y_i\Delta z_i,\\ A_{49}=\frac{1}{2}{w}_i\Delta x_i^3\Delta z_i,\\ A_{57}=\frac{1}{2}{w}_i\Delta x_i\Delta y_i^3,\\ A_{58}=\frac{1}{2}{w}_i\Delta y_i^3\Delta z_i,\\ A_{59}=\frac{1}{2}{w}_i\Delta x_i\Delta y_i^2\Delta z_i,\\ A_{67}=\frac{1}{2}{w}_i\Delta x_i\Delta y_i\Delta z_i^2,\\ A_{68}=\frac{1}{2}{w}_i\Delta y_i\Delta z_i^3,\\ A_{69}=\frac{1}{2}{w}_i\Delta x_i\Delta z_i^3,\\ A_{77}=w_i\Delta x_i^2\Delta y_i^2,\\ A_{78}=w_i\Delta x_i\Delta y_i^2\Delta z_i,\\ A_{79}=w_i\Delta x_i^2\Delta y_i\Delta z_i,\\ A_{88}=w_i\Delta y_i^2\Delta z_i^2,\\ A_{89}=w_i\Delta x_i\Delta y_i\Delta z_i^2,\\ A_{99}=w_i\Delta x_i^2\Delta z_i^2.\end{array}$$

In addition, the values of each element of the x matrix and b matrix can be expressed as follows:$$\begin{array}{c}\text{(A.6)}& x_{11}=\frac{\partial Q}{\partial x},x_{21}=\frac{\partial Q}{\partial y},\\ x_{31}=\frac{\partial Q}{\partial z},\\ x_{41}=\frac{{\partial }^{2}Q}{\partial x^2},\\ x_{51}=\frac{{\partial }^{2}Q}{\partial y^2},\\ x_{61}=\frac{{\partial }^{2}Q}{\partial z^2},\\ x_{71}=\frac{{\partial }^{2}Q}{\partial x\partial y},\\ x_{81}=\frac{{\partial }^{2}Q}{\partial y\partial z},\\ x_{91}=\frac{{\partial }^{2}Q}{\partial x\partial z},\\ \text{(A.7)}& b_{11}=w_i\Delta x_i\Delta Q_i,b_{21}=w_i\Delta y_i\Delta Q_i,\\ b_{31}=w_i\Delta z_i\Delta Q_i,\\ b_{41}=\frac{1}{2}{w}_i\Delta x_i^2\Delta Q_i,\\ b_{51}=\frac{1}{2}{w}_i\Delta y_i^2\Delta Q_i,\\ b_{61}=\frac{1}{2}{w}_i\Delta z_i^2\Delta Q_i,\\ b_{71}=w_i\Delta x_i\Delta y_i\Delta Q_i,\\ b_{81}=w_i\Delta y_i\Delta z_i\Delta Q_i,\\ b_{91}=w_i\Delta x_i\Delta z_i\Delta Q_i.\end{array}$$

References

[1] M. S. Longuet-Higgins, E. D. Cokelet, "The deformation of steep surface waves on water. I: a numerical method of computation," Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 350 no. 1660,DOI: 10.1098/rspa.1976.0092, 1976.

[2] C. Yang, R. C. Ertekin, "Numerical simulation of nonlinear wave diffraction by a vertical cylinder," Journal of Offshore Mechanics and Arctic Engineering, vol. 114 no. 1,DOI: 10.1115/1.2919950, 1992.

[3] A. Sommerfeld, Partial Differential Equations in Physics, 1949.

[4] S. Y. Boo, C. H. Kim, "Fully nonlinear diffraction due to a vertical circular cylinder in a 3-D HOBEM numerical wave tank," Proceedings of the Sixth International Offshore and Polar Engineering Conference, .

[5] S. Y. Boo, "Linear and nonlinear irregular waves and forces in a numerical wave tank," Ocean Engineering, vol. 29 no. 5, pp. 475-493, DOI: 10.1016/s0029-8018(01)00055-5, 2002.

[6] P. Ferrant, D. L. Touzé, K. Pelletier, "Non-linear time-domain models for irregular wave diffraction about offshore structures," International Journal for Numerical Methods in Fluids, vol. 43 no. 10-11, pp. 1257-1277, DOI: 10.1002/fld.506, 2003.

[7] W. Bai, R. Eatock Taylor, "Numerical simulation of fully nonlinear regular and focused wave diffraction around a vertical cylinder using domain decomposition," Applied Ocean Research, vol. 29 no. 1-2, pp. 55-71, DOI: 10.1016/j.apor.2007.05.005, 2007.

[8] J. Y. T. Ng, M. Issacson, "Second-order wave interaction with two-dimensional floating bodies by a time-domain method," Applied Ocean Research, vol. 15 no. 2, pp. 95-105, DOI: 10.1016/0141-1187(93)90024-r, 1993.

[9] M. Issacson, J. Y. T. Ng, "Time-domain second-order wave interaction with three-dimensional floating bodies," International Journal of Offshore and Polar Engineering, vol. 5 no. 3, 1995.

[10] G. X. Wu, R. Eatock Taylor, "Time stepping solutions of the two-dimensional nonlinear wave radiation problem," Ocean Engineering, vol. 22 no. 8, pp. 785-798, 1995.

[11] S. Maiti, D. Sen, "Time-domain wave diffraction of two-dimensional single and twin hulls," Ocean Engineering, vol. 28 no. 6, pp. 639-665, DOI: 10.1016/s0029-8018(00)00025-1, 2001.

[12] W. C. Koo, M. H. Kim, "Numerical simulation of nonlinear wave and force generated by a wedge-shape wave maker," Ocean Engineering, vol. 33 no. 8-9, pp. 983-1006, DOI: 10.1016/j.oceaneng.2005.09.002, 2006.

[13] R. Datta, D. Sen, "A B-spline-based method for radiation and diffraction problems," Ocean Engineering, vol. 33 no. 17-18, pp. 2240-2259, DOI: 10.1016/j.oceaneng.2005.12.002, 2006.

[14] W. Bai, R. Eatock Taylor, "Higher-order boundary element simulation of fully nonlinear wave radiation by oscillating vertical cylinders," Applied Ocean Research, vol. 28 no. 4, pp. 247-265, DOI: 10.1016/j.apor.2006.12.001, 2006.

[15] B. Z. Zhou, D. Z. Ning, B. Teng, W. Bai, "Numerical investigation of wave radiation by a vertical cylinder using a fully nonlinear HOBEM," Ocean Engineering, vol. 70,DOI: 10.1016/j.oceaneng.2013.04.019, 2013.

[16] S. Oh, S. K. Cho, D. Jung, H. G. Sung, "Development and application of two-dimensional numerical wave tank using desingularized indirect boundary integral equation method," Journal of Ocean and Technology, vol. 32 no. 6, pp. 447-457, DOI: 10.26748/ksoe.2018.32.6.447, 2018.

[17] A. Feng, Y. You, H. Cai, "An improved Rankine source panel method for three dimensional water wave problems," International Journal of Naval Architecture and Ocean Engineering, vol. 11 no. 1, pp. 70-81, DOI: 10.1016/j.ijnaoe.2018.02.001, 2019.

[18] W. Koo, M.-H. Kim, "Freely floating-body simulation by a 2D fully nonlinear numerical wave tank," Ocean Engineering, vol. 31 no. 16, pp. 2011-2046, DOI: 10.1016/j.oceaneng.2004.05.003, 2004.

[19] W. C. Koo, M. H. Kim, "Fully nonlinear wave-body interactions with surface-piercing bodies," Ocean Engineering, vol. 34 no. 7, pp. 1000-1012, DOI: 10.1016/j.oceaneng.2006.04.009, 2007.

[20] A. Abbasnia, M. Ghiasi, "Simulation of nonlinear wave interaction with dual cylinders in numerical wave tank," Transactions of FAMENA, vol. 37 no. 1, pp. 35-48, 2013.

[21] A. Abbasnia, M. Ghaisi, J. Barandiaran, C. Guedes Soares, "An implicit model of a submerged horizontal cylinder oscillating about an off-centered axis as a wave energy converter," Renewable Energies Offsore, pp. 247-255, 2015.

[22] A. Abbasnia, M. Ghiasi, A. H. Abbasnia, "Irregular wave transmission on bottom bumps using fully nonlinear NURBS numerical wave tank," Engineering Analysis with Boundary Elements, vol. 82, pp. 130-140, DOI: 10.1016/j.enganabound.2017.06.008, 2017.

[23] W. Bai, R. Eatock Taylor, "Fully nonlinear simulation of wave interaction with fixed and floating flared structures," Ocean Engineering, vol. 36 no. 3-4, pp. 223-236, DOI: 10.1016/j.oceaneng.2008.11.003, 2009.

[24] G. Wu, R. Eatock Taylor, "The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies," Ocean Engineering, vol. 30, pp. 387-400, 2003.

[25] B. Z. Zhou, G. X. Wu, B. Teng, "Fully nonlinear wave interaction with freely floating non-wall-sided structures," Engineering Analysis with Boundary Elements, vol. 50, pp. 117-132, DOI: 10.1016/j.enganabound.2014.08.003, 2015.

[26] B. Z. Zhou, G. X. Wu, Q. C. Meng, "Interactions of fully nonlinear solitary wave with a freely floating vertical cylinder," Engineering Analysis with Boundary Elements, vol. 69, pp. 119-131, DOI: 10.1016/j.enganabound.2016.05.004, 2016.

[27] A. Abbasnia, C. Guedes Soares, "Exact evaluation of hydrodynamic loads on ships using NURBS surfaces and acceleration potential," Engineering Analysis with Boundary Elements, vol. 85,DOI: 10.1016/j.enganabound.2017.09.005, 2017.

[28] H. G. Sung, S. T. Grilli, "BEM computations of 3D fully nonlinear free surface flows caused by advancing surface disturbances," International Journal of Offshore and Polar Engineering, vol. 18 no. 4, pp. 292-301, 2008.

[29] H. Yan, Y. Liu, "An efficient high-order boundary element method for nonlinear wave–wave and wave-body interactions," Journal of Computational Physics, vol. 230 no. 2, pp. 402-424, DOI: 10.1016/j.jcp.2010.09.029, 2011.

[30] K.-R. Lee, W. Koo, M.-H. Kim, "Fully nonlinear time-domain simulation of a backward bent duct buoy floating wave energy converter using an acceleration potential method," International Journal of Naval Architecture and Ocean Engineering, vol. 5 no. 4, pp. 513-528, DOI: 10.2478/ijnaoe-2013-0150, 2013.

[31] S.-J. Kim, W. Koo, M.-H. Kim, "Nonlinear time-domain NWT simulations for two types of a backward bent duct buoy (BBDB) compared with 2D wave-tank experiments," Ocean Engineering, vol. 108, pp. 584-593, DOI: 10.1016/j.oceaneng.2015.08.038, 2015.

[32] L. Letournel, G. Ducrozet, A. Babarit, P. Ferrant, "Proof of the equivalence of Tanizawa–Berkvens’ and Cointe–van Daalen’s formulations for the time derivative of the velocity potential for non-linear potential flow solvers," Applied Ocean Research, vol. 63, pp. 184-199, DOI: 10.1016/j.apor.2017.01.010, 2017.

[33] S.-J. Kim, W. Koo, M.-J. Shin, "Numerical and experimental study on a hemispheric point-absorber-type wave energy converter with a hydraulic power take-off system," Renewable Energy, vol. 135, pp. 1260-1269, DOI: 10.1016/j.renene.2018.09.097, 2019.

[34] J. L. Hess, A. M. O. Smith, "Calculation of nonlifting potential flow about arbitrary three-dimensional bodies," 1962. Report no. E.S. 40622

[35] J. N. Newman, "Distributions of sources and normal dipoles over a quadrilateral panel," Journal of Engineering Mathematics, vol. 20 no. 2, pp. 113-216, DOI: 10.1007/bf00042771, 1986.

[36] S. Sen, "Efficient computation of influence coefficients in Rankine-panel based methods of computational ship hydrodynamics," Oceanic Engineering International, vol. 5 no. 1, pp. 16-30, 2000.

[37] G. R. Baker, D. J. Merion, S. A. Orazag, "Applications of a generalized vortex method to nonlinear free-surface flows," Proceedings of the Third International Conference on Numerical Ship Hydrodynamics, pp. 179-191, .

[38] R. Cointe, P. Geyer, B. King, B. Molin, M. Tramoni, "Nonlinear and linear motions of a rectangular barge in a perfect fluid," Proceedings of the 18th Symposium on Naval Hydrodynamics, pp. 85-99, .

[39] S. T. Grilli, J. Horillo, "Periodic wave shoaling over barred-beaches in a fully nonlinear numerical wave tank," Proceedings of the 8th International Offshore and Polar Engineering Conference, Montreal (ISOPE), vol. 3, pp. 294-300, .

[40] S. A. Hong, M. H. Kim, "Nonlinear wave forces on a stationary vertical cylinder by HOBEM-NWT," Proceedings of the Tenth International Offshore and Polar Engineering Conference, vol. 3, pp. 214-220, .

[41] M. Israeli, S. A. Orszag, "Approximation of radiation boundary conditions," Journal of Computational Physics, vol. 41 no. 1, pp. 115-135, DOI: 10.1016/0021-9991(81)90082-6, 1981.

[42] M. W. Kim, W. Koo, S. Y. Hong, "Numerical analysis of various artificial damping schemes in a three-dimensional numerical wave tank," Ocean Engineering, vol. 75, pp. 165-173, DOI: 10.1016/j.oceaneng.2013.10.012, 2014.

[43] E. Sozer, C. Brehm, C. C. Kiris, "Gradient calculation methods on arbitrary polyhedral unstructured meshes for cell-centered CFD solvers," Proceedings of the 52nd Aerospace Sciences Meeting, AIAA SciTech Forum, (AIAA 2014-1440),DOI: 10.2514/6.2014-1440, .

[44] K. Tanizawa, "A nonlinear simulation method of 3-D body motions in waves (1st report)," Journal of the Society of Naval Architects of Japan, vol. 1995 no. 178, pp. 179-191, DOI: 10.2534/jjasnaoe1968.1995.178_179, 1995.

[45] R. Cointe, "Nonlinear simulation of transient free surface flows," Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, pp. 85-99, .

[46] S. Sen, G. De Nayer, M. Breuer, "A fast and robust hybrid method for block-structured mesh deformation with emphasis on FSI-LES applications," International Journal for Numerical Methods in Engineering, vol. 111 no. 3, pp. 273-300, DOI: 10.1002/nme.5465, 2017.

[47] A. Abbasnia, M. Ghiasi, "Fully nonlinear wave interaction with an array of truncated barriers in three dimensional numerical wave tank," Engineering Analysis with Boundary Elements, vol. 58, pp. 79-85, DOI: 10.1016/j.enganabound.2015.03.015, 2015.

[48] Y. Saad, M. H. Schultz, "GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems," SIAM Journal of Science Statistical Computation, vol. 7 no. 3, pp. 856-869, DOI: 10.1137/0907058, 1986.

[49] A. F. D. O. Falcão, "Modelling and control of oscillating-body wave energy converters with hydraulic power take-off and gas accumulator," Ocean Engineering, vol. 34 no. 14-15, pp. 2021-2032, DOI: 10.1016/j.oceaneng.2007.02.006, 2007.

[50] M. W. Dingemans, Water Wave Propagation over Uneven Bottoms: Linear Wave Propagation, 1997.

[51] W. Bai, B. Teng, "Simulation of second-order wave interaction with fixed and floating structures in time domain," Ocean Engineering, vol. 74, pp. 168-177, DOI: 10.1016/j.oceaneng.2013.07.014, 2013.

[52] B. Teng, "Second order wave action on 3D floating body," Journal of Hydrodynamics, Series A, vol. 10 no. 3, pp. 316-327, 1995.

[53] W. Bai, "Non-linear wave interaction with arbitary 3-D bodies," . Ph.D. thesis, Dalian University of Technology, Dalian, China, 2001, in Chinese

[54] Y.-L. Shao, O. M. Faltinsen, "Use of body-fixed coordinate system in analysis of weakly nonlinear wave-body problems," Applied Ocean Research, vol. 32 no. 1, pp. 20-33, DOI: 10.1016/j.apor.2010.05.004, 2010.

[55] J. Kim, I. H. Cho, M. H. Kim, "Numerical calculation and experiment of a heaving-buoy wave energy converter with a latching control," Ocean System Engineering, vol. 9 no. 1,DOI: 10.12989/ose.2019.9.1.001, 2019.

[56] S. J. Kim, W. Koo, C. H. Jo, "Assessment of latching control for the hemispheric heaving buoy type point absorber with and without nonlinear Froude-Krylov force acting on the buoy," Proceedings of the ASME 38th OMAE 2019, .