A two-dimensional boundary element method with generating absorbing boundary condition for floating bodies of arbitrary shape in the frequency domain

Abstract

A two-dimensional (2D) boundary element method is developed for the rapid assessment of the hydrodynamic performance of floating structures in waves. The boundary element method is based on potential flow and has panels along all boundaries of the fluid domain – not only along the boundary of the floater – to make the extension to second order feasible. Panels along all boundaries requires the development of generating absorbing boundary conditions for use at radiation boundaries to send incident waves into the domain while absorbing waves originating from the floating body at the same boundary, at the same time. The model is verified by means of conservation of energy of a heaving wave energy converter, and by means of the propagation of second-order waves. The performance in terms of conservation of energy with 12 panels per wave length is good, the generating absorbing boundary condition works according to expectation and the second-order wave propagation corresponds to theory.

1. Background

It is convenient to have a rapid assessment method for the performance of floating structures in free surface waves. Rapid assessment was shown to capture much of the physical behaviour of wave interaction with a non-moving wall in Bos et al. [2]. While Xu and Wellens [14,15] introduced rapid assessment methods for flexible floating structures that need to be large compared to the wave lenghth, we aim to develop a method for floating rigid structures of any shape without a restriction in size. Boundary element methods have an order fewer unknowns compared to field methods (Computational Fluid Dynamics, CFD), making them attractive in terms of computational cost.

A boundary element method (BEM) discretizes the problem at the boundaries to solve the fluid domain. Only potential flow is considered. Potential flow theory is based on two fundamental equations, the conservation of mass and the conservation of momentum, that need to be solved in the fluid domain. The equation for conservation of momentum stems from an equation of motion for the fluid.

Boundary conditions need to bo added to the fundamental equations. One is that the water particles can not enter the body, so there is a no-penetration boundary condition on the body: the velocity of the water particles is the same as the velocity of the body. A similar boundary condition applies to free surface and bottom. The side walls of the fluid domain need to let incident waves in and let radiated and diffracted waves out. The body has its own equations of motion which completes the system of equations.

The fundamental equations with boundary conditions can be solved in different ways. There is the option to solve the system fully analytically. This has been done for a cylinder by Dean [3] and Ursell [11], Ursell and Dean [12] up to first order. The mathematics required to solve the system are extensive. Ogilvie [9] extended the approach to include second order forces on a stationary cylinder, and on a submerged cylinder that was free to respond. The analytical approaches, especially to second order, have only be applied to a limited set of geometries that are cylindrical of elliptical in shape.

Our chosen method for solving the system of equations for more complex shapes is by discretization with panels along a boundary of interest. Sources are used to satisfy the boundary condition at hand at the position of each panel. This results in a system of equations with an equation for each panel that solves for its source strength. The total potential that describes the flow in the entire fluid domain can be constructed from the combination of all source strengths. Many examples of panel methods exist.

Frank [4] used sources with Green’s functions to solve the boundary conditions at the free surface and at infinity. An advantage of using Green’s functions is that only the body surface itself needs to be approximated with panels. A more involved example comes from Kim, Kim and Lucas [7] and Kim and Yue [6]. They make use of a Green’s function to solve the second order diffraction problem for the pressures and body motions in three dimensions. The software WAMIT uses the methods from this work to calculate the second order flow around floating bodies. An advantage of these methods is that the number of unknowns is low (only the body); a disadvantage is that the mathematics are complex and that extension to third order is even more complex.

A different approach is taken by Yeung [16] and Ballast [1]. Instead of using a Green’s function to solve for the boundary conditions, simpler sources are used at the positions of the panels. The boundary conditions at the free surface and at infinity are solved by having panels on the free surface and including radiation boundary conditions at the side walls of the fluid domain. This increases the number of panels, but reduces the complexity of the calculations. The method of Yeung [16] is suitable for two and three dimensions, but first order and limited to the radiation problem only (no incident waves and diffraction).

Our objective is to develop a fast method with radiation and diffraction in the frequency domain that can be extended to second order. To make this possible, the approach is to place panels along all boundaries of the fluid domain and not only along the boundary of the floating object as is commonly done (Frank [4]). The approach implies that radiation conditions, also called absorbing boundary conditions (ABC), are imposed on the side walls of the fluid domain for the waves that leave the domain in the radiation problem. For the diffraction problem it is required that diffracted waves leave the domain over the same boundary and at the same time as incident waves enter the fluid domain. Following Wellens and Borsboom [13], a generating absorbing boundary condition (GABC) has been derived for the diffraction problem.

Table 1 shows how the research in this article relates to existing literature with boundary element methods. The table compares methods based on a number of characteristics that are important for what we would like to achieve. The first characteristics say whether the methods are 2D, 3D or both. Some methods feature incident waves, others do not. ‘Simple sources’ describes whether a method uses a Green’s function which takes care of the boundary conditions of the fluid domain or not. Whether a method works in the frequency domain or in the time domain is an indication of how fast it is. ‘Coupling’ indicates whether incident waves and radiated waves can enter and leave the domain simultaneously. And ‘second order’ states whether a method is first or second order. Our method is ready for second order: it features second-order wave propagation, but does not yet solve the fluid structure interaction with the body to second order. In summary, the method presented in this article is the only frequency domain method with absorbing boundary conditions and an extension towards second order.

Table 1

State of the art BEM

	2D	3D	Incoming waves	Simple sources	Frequency domain	(G)ABC	Second-order
Kim and Yue [6]	No	Yes	Yes	No	No	No	Yes
Yeung [16]	Yes	Yes	No	Yes	Yes	No	No
Frank [4]	Yes	No	Yes	No	Yes	No	No
Ballast [1]	No	Yes	Yes	Yes	No	Yes	Yes
This research	Yes	No	Yes	Yes	Yes	Yes	Ready

The next chapter presents the mathematical formulation of the problem. Section 3 describes how the governing equations have been discretized and included in the method. The method is evaluated in Section 4 with the first-order response of a floating wave energy converter and with second-order wave propagation. The conclusions are summarized in the final chapter.

2. Mathematical formulation

2.1. Fundamental equations and boundary conditions

A potential Φ is defined such that the gradient of the potential gives the velocity vector, $\nabla \mathrm{\Phi }=\overrightarrow{\mathit{u}}$ . Conservation of mass then becomes $$\begin{array}{}\text{(1)}& \hfill {\nabla }^2\mathrm{\Phi }=0.& \end{array}$$

Conservation of momentum is obtained through the Bernoulli equation $$\begin{array}{}\text{(2)}& \hfill {\displaystyle \frac{\mathit{p}}{\mathit{\rho }}}=-{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{t}}}-{\displaystyle \frac{1}{2}}|\nabla \mathrm{\Phi }{|}^2-\mathrm{g}\mathit{y}.& \end{array}$$ Here, p is the pressure, t is time, g is the acceleration of gravity and y is the vertical coordinate that is defined as positive in opposite direction to the force of gravity.

A fluid domain is given in Fig. 1. It features the right-handed axis system with y, and x in the direction orthogonal to y. The total potential Φ is subdivided into ${\mathrm{\Phi }}_{\mathit{I}}$ , ${\mathrm{\Phi }}_{\mathit{D}}$ and ${\mathrm{\Phi }}_{\mathit{R}}$ . The explanation of these symbols and the symbols for identifying the boundaries in Fig. 1 is given in Table 2.

Fig. 1.

A fluid domain.

Table 2

List of symbols for components of the total potential and the boundaries of the fluid domain

Abbreviation	Physical meaning
${\mathrm{\Phi }}_{\mathit{I}}$	Incident wave
${\mathrm{\Phi }}_{\mathit{D}}$	Diffracted wave
${\mathrm{\Phi }}_{\mathit{R}}$	Radiated wave
${\mathit{S}}_{\mathit{B}}$	Body surface
${\mathit{S}}_{\mathit{F}}$	Free surface
${\mathit{S}}_{\mathit{R}}$	Radiation boundary
${\mathit{B}}_0$	Bottom boundary

At the bottom boundary ${\mathit{B}}_0$ , with position $\mathit{y}=-\mathit{h}$ and h being the water depth, a no-penetration boundary condition is solved $$\begin{array}{}\text{(3)}& \hfill {\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{y}}}=0\text{at }\mathit{y}=-\mathit{h}.& \end{array}$$

A kinematic condition is solved at free surface ${\mathit{S}}_{\mathit{F}}$ $$\begin{array}{}\text{(4)}& \hfill {\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{t}}}+\nabla \mathrm{\Phi }·\nabla \mathit{\eta }={\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{y}}}\text{at }\mathit{y}=\mathit{\eta }\mathrm{,}& \end{array}$$ with η the free surface. The kinematic condition states that fluid particles at the position of the free surface follow the motion of the free surface.

A highlight of this paper is the generation of waves at radiation boundaries ${\mathit{S}}_{\mathit{R}}$ . To achieve this, the well-known Sommerfeld radiation condition, also known as absorbing boundary condition (ABC) is changed in to generating absorbing boundary condition (GABC) as in Wellens and Borsboom [13]. The Sommerfeld condition is based on a propagating wave in a certain direction with phase velocity c. At the left radiation boundary, waves that propagate to the left are absorbed; at the right radiation boundary waves that propagate to the right are absorbed. For a wave propagating to the right, at the right-hand radiation boundary, the Sommerfeld condition reads $$\begin{array}{}\text{(5)}& \hfill \mathrm{(}{\displaystyle \frac{\partial }{\partial \mathit{t}}}+\mathit{c}{\displaystyle \frac{\partial }{\partial \mathit{x}}}\mathrm{)}\mathrm{\Phi }=0.& \end{array}$$

The wave speed c is defined as $\mathit{\omega }\mathrm{/}\mathit{k}$ , with ω the wave frequency and k the wave number. In the frequency domain, for a wave mode with argument $\mathrm{(}\mathit{\omega }\mathit{t}-\mathit{k}\mathit{x}\mathrm{)}$ , Eq. (5) simplifies to $$\begin{array}{}\text{(6)}& \hfill \mathrm{(}{\displaystyle \frac{\partial }{\partial \mathit{x}}}+\mathrm{i}\mathit{k}\mathrm{)}\mathrm{\Phi }=0.& \end{array}$$

This is also the expression used by Yeung [16] for his radiation problem. For the diffraction problem, the left-hand radiation boundary is also a generation boundary for the incident waves propagating to the right. Following Wellens and Borsboom [13], who derived a GABC for irregular waves in the time domain, the right-hand side of Eq. (6) is given a non-zero value equal to the operators on the left-hand side applied to a wave mode with argument $\mathrm{(}\mathit{\omega }\mathit{t}+\mathit{k}\mathit{x}\mathrm{)}$ , yielding a GABC for the frequency domain $$\begin{array}{}\text{(7)}& \hfill \mathrm{(}{\displaystyle \frac{\partial }{\partial \mathit{x}}}-\mathrm{i}\mathit{k}\mathrm{)}\mathrm{\Phi }=\mathrm{(}{\displaystyle \frac{\partial }{\partial \mathit{x}}}-\mathrm{i}\mathit{k}\mathrm{)}{\mathrm{\Phi }}_{\mathit{I}}=-2\mathrm{i}\mathit{k}{\mathrm{\Phi }}_{\mathit{I}}& \end{array}$$

The performance of boundary condition Eq. (7) is demonstrated in Section 4.

2.2. Linear wave theory

It is assumed that the waves are not steep, and are long enough to neglect surface tension. For linear potential flow theory all equations are linearized around $\mathit{y}=0$ , meaning higher than first-order terms are neglected. The kinematic boundary condition then becomes $$\begin{array}{}\text{(8)}& \hfill {\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{t}}}={\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{y}}}\text{at }\mathit{y}=0.& \end{array}$$

The atmospheric pressure at the free surface is assumed to be zero $$\begin{array}{}\text{(9)}& \hfill \mathit{p}=0\text{at }\mathit{y}=\mathit{\eta }.& \end{array}$$

Linearizing the Bernoulli equation, Eq. (2), and substituting Eq. (9), leads to the dynamic free surface boundary condition $$\begin{array}{}\text{(10)}& \hfill {\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{t}}}+\mathrm{g}\mathit{\eta }=0\text{at }\mathit{y}=0.& \end{array}$$

Linear potential flow theory is also known as Airy wave theory. The potential function of an Airy wave $$\begin{array}{}\text{(11)}& \hfill \mathit{\eta }=\mathit{A}sin\mathrm{(}\mathit{\omega }\mathit{t}-\mathit{k}\mathit{x}\mathrm{)}& \end{array}$$ that satisfies the fundamental equations and the linearized boundary conditions is $$\begin{array}{}\text{(12)}& \hfill \mathrm{\Phi }={\displaystyle \frac{\mathit{A}\mathrm{g}}{\mathit{\omega }}}{\displaystyle \frac{cosh\mathrm{(}\mathit{k}\mathrm{(}\mathit{h}+\mathit{y}\mathrm{)}\mathrm{)}}{cosh\mathrm{(}\mathit{k}\mathit{h}\mathrm{)}}}cos\mathrm{(}\mathit{\omega }\mathit{t}-\mathit{k}\mathit{x}\mathrm{)}\mathrm{,}& \end{array}$$ with A the amplitude of the Airy wave.

The characteristic equation of the system of equations is the dispersion relation $$\begin{array}{}\text{(13)}& \hfill {\mathit{\omega }}^2=\mathit{k}\mathrm{g}tanh\mathrm{(}\mathit{k}\mathit{h}\mathrm{)}.& \end{array}$$

It relates the wave frequency and the wave number. Because $\mathit{c}=\mathit{\omega }\mathrm{/}\mathit{k}$ , the dispersion relation implies that waves of different frequency propagate at different phase velocities c. The wave length is found from the wave number as $\mathit{L}=2\mathit{\pi }\mathrm{/}\mathit{k}$ and the wave period is found from the wave frequency as $\mathit{T}=2\mathit{\pi }\mathrm{/}\mathit{\omega }$ .

2.3. Second-order wave theory

In order to solve the non-linear wave propagation problem, it is common to use perturbation theory. This was first done by Stokes [10] for free surface gravity waves. It was repeated by Newman [8] for waves in infinite water depth. With perturbation theory, a function or solution variable is approximated by a series of terms of increasing order $$\begin{array}{}\text{(14)}& \hfill \mathit{\xi }={\mathit{\xi }}^{\mathrm{(}0\mathrm{)}}+\mathit{\epsilon }{\mathit{\xi }}^{\mathrm{(}1\mathrm{)}}+{\mathit{\epsilon }}^2{\mathit{\xi }}^{\mathrm{(}2\mathrm{)}}+\cdots & \end{array}$$

Here, ${\mathit{\xi }}^{\mathrm{(}\mathit{i}\mathrm{)}}$ are the expansion terms in the perturbation series approximating solution variable ξ. ε is a small term (for waves it is the wave steepness), so that higher-order terms become smaller as the order increases. An approximate solution is found by truncating the series. From now on, approximated variables are given a superscript $\mathrm{(}1\mathrm{)}$ or $\mathrm{(}2\mathrm{)}$ to indicate the order of approximation.

The second-order boundary condition (Stokes [10]) is referred to as the quadratic forcing function (QFF) $$\begin{array}{}\text{(15)}& \hfill -4\mathit{\omega }{\mathrm{\Phi }}^{\mathrm{(}2\mathrm{)}}+\mathrm{g}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\mathrm{(}2\mathrm{)}}}{\partial \mathit{y}}}=\text{QFF}\mathrm{,}& \end{array}$$ with $$\begin{array}{ll}\hfill \text{QFF}=& -{\displaystyle \frac{\partial }{\partial \mathit{t}}}\mathrm{(}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\mathrm{(}1\mathrm{)}}}{\partial \mathit{x}}}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\mathrm{(}1\mathrm{)}}}{\partial \mathit{x}}}+{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\mathrm{(}1\mathrm{)}}}{\partial \mathit{y}}}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\mathrm{(}1\mathrm{)}}}{\partial \mathit{y}}}\mathrm{)}\\ \text{(16)}& \hfill & +{\displaystyle \frac{1}{\mathrm{g}}}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\mathrm{(}1\mathrm{)}}}{\partial \mathit{t}}}{\displaystyle \frac{\partial }{\partial \mathit{y}}}\mathrm{(}{\displaystyle \frac{{\partial }^2{\mathrm{\Phi }}^{\mathrm{(}1\mathrm{)}}}{\partial {\mathit{t}}^2}}+\mathrm{g}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\mathrm{(}1\mathrm{)}}}{\partial \mathit{y}}}\mathrm{)}.\end{array}$$

An index to find the non-linearity for waves in finite water depth is the Ursell number $$\begin{array}{}\text{(17)}& \hfill {\mathit{N}}_{\mathrm{Ursell}}={\displaystyle \frac{\mathit{A}\mathrm{/}\mathit{L}}{{\mathrm{(}\mathit{h}\mathrm{/}\mathit{L}\mathrm{)}}^3}}={\displaystyle \frac{\mathit{A}{\mathit{L}}^2}{{\mathit{h}}^3}}.& \end{array}$$ It is defined as the steepness/(relative depth)³. The higher the Ursell number, the more non-linear a wave is going to behave. The second-order potential that satisfies the second-order boundary condition at the free surface (Stokes [10]) is $$\begin{array}{}\text{(18)}& \hfill {\mathrm{\Phi }}^{\mathrm{(}2\mathrm{)}}={\displaystyle \frac{3{\mathit{A}}^2\mathit{\omega }}{8}}{\displaystyle \frac{cosh\mathrm{(}2\mathit{k}\mathrm{(}\mathit{h}+\mathit{y}\mathrm{)}\mathrm{)}}{{sinh}^4\mathrm{(}\mathit{k}\mathit{h}\mathrm{)}}}sin\mathrm{(}2\mathrm{(}\mathit{\omega }\mathit{t}-\mathit{k}\mathit{x}\mathrm{)}\mathrm{)}& \end{array}$$

2.4. Pressures and forces

To find the pressures in the domain use is made of the linearized Bernoulli equation $$\begin{array}{}\text{(19)}& \hfill {\mathit{p}}^{\mathrm{(}1\mathrm{)}}=-\mathit{\rho }{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\mathrm{(}1\mathrm{)}}}{\partial \mathit{t}}}-\mathit{\rho }\mathrm{g}\mathit{y}.& \end{array}$$

The forces on the body are the integration of the pressures along the body surface ${\mathit{S}}_{\mathit{b}}$ , see Fig. 1 $$\begin{array}{}\text{(20)}& \hfill {\mathit{F}}^{\mathrm{(}1\mathrm{)}}={{\displaystyle \int }}_{{\mathit{S}}_{\mathit{B}}^{\mathrm{(}0\mathrm{)}}}{\mathit{p}}^{\mathrm{(}1\mathrm{)}}\mathrm{d}\mathit{s}.& \end{array}$$

The buoyancy force on the floating body is found differently. The change in buoyancy in the heave direction is caught by the hydrostatic term in the Bernoulli equation $\mathit{\rho }\mathit{g}\mathit{y}$ , in which y is replaced by the vertical position a of the floater. In the presence of waves, the pressure depends on the relative position of the structure with respect to the free surface elevation $\mathit{\eta }-\mathit{a}$ . The pressure is then integrated over the waterplane area ${\mathit{S}}_{\mathit{A}}$ $$\begin{array}{}\text{(21)}& \hfill {\mathit{F}}_{\mathit{b}}={{\displaystyle \int }}_{{\mathit{S}}_{\mathit{A}}}\mathit{\rho }\mathrm{g}\mathrm{(}{\mathit{\eta }}^{\mathrm{(}1\mathrm{)}}-{\mathit{a}}^{\mathrm{(}1\mathrm{)}}\mathrm{)}\mathrm{d}\mathit{s}.& \end{array}$$

This expression for the buoyancy force is accurate for wall-sided structures and small motions. For non-wall-sided structures it is the zeroth-order approximation. Other motion directions require a similar evaluation.

2.5. Motions

The motions of the floating object are found from an equation of motion for the body. For the vertical direction $$\begin{array}{}\text{(22)}& \hfill {\mathit{\zeta }}_{\mathit{y}}={\displaystyle \frac{{\mathit{F}}_{\mathit{y}}^{\mathrm{(}1\mathrm{)}}}{-{\mathit{\omega }}^2\mathrm{(}\mathit{M}+\mathit{\mu }\mathrm{)}-\mathrm{i}\mathit{\omega }\mathit{\lambda }+\mathit{K}}}\mathrm{,}& \end{array}$$ in which ${\mathit{\zeta }}_{\mathit{y}}$ is the response amplitude operator (RAO), defined as the ratio between the body motion amplitude $\hat{\mathit{a}}$ and the wave amplitude A. M is the body mass, μ the added mass, λ the damping and K the coefficient of the buoyancy. The added mass and damping of the floating body are found from a radiation simulation. The total force on the floater is found from the sum of forces found from a radiation simulation and a diffraction simulation.

2.6. Power prediction & optimal damping

When the floating body is equipped with a generator to harvest wave energy, a power take-off (PTO) device, the power produced can be estimated using the parameters from a BEM simulation. Assume that ${\mathit{\lambda }}_{\mathrm{PTO}}$ is the damping coefficient of the PTO. The produced power as a function of time then becomes $$\begin{array}{}\text{(23)}& \hfill {\mathit{P}}_{\mathit{A}}\mathrm{(}\mathit{t}\mathrm{)}={\mathit{\lambda }}_{\text{PTO}}\dot{\mathit{a}}{\mathrm{(}\mathit{t}\mathrm{)}}^2.& \end{array}$$

Velocity $\dot{\mathit{a}}\mathrm{(}\mathit{t}\mathrm{)}$ is the time derivative of motion $\mathit{a}\mathrm{(}\mathit{t}\mathrm{)}$ of the body, with $\hat{\mathit{a}}$ its amplitude.

The average power in one wave period T then becomes $$\begin{array}{}\text{(24a)}& \hfill & {\bar{\mathit{P}}}_{\mathit{A}}={\displaystyle \frac{1}{\mathit{T}}}{{\displaystyle \int }}_0^{\mathit{T}}{\mathit{P}}_{\mathit{d}}\mathrm{(}\mathit{t}\mathrm{)}\mathrm{d}\mathit{t}\text{(24b)}& \hfill & {\bar{\mathit{P}}}_{\mathit{A}}={\displaystyle \frac{1}{\mathit{T}}}{{\displaystyle \int }}_0^{\mathit{T}}{\mathit{\lambda }}_{\text{PTO}}{\mathit{\omega }}^2{\hat{\mathit{a}}}^2{sin}^2\mathrm{(}\mathit{\omega }\mathit{t}\mathrm{)}\mathrm{d}\mathit{t}\text{(24c)}& \hfill & {\bar{\mathit{P}}}_{\mathit{A}}={\displaystyle \frac{1}{2}}{\mathit{\lambda }}_{\text{PTO}}{\mathit{\omega }}^2{\hat{\mathit{a}}}^2.\end{array}$$

Using RAO ${\mathit{\zeta }}_{\mathit{y}}$ with the additional damping ${\mathit{\lambda }}_{\text{PTO}}$ of the generator, an expression for the average power in terms of the damping is found $$\begin{array}{}\text{(25a)}& \hfill & {\mathit{\zeta }}_{\mathit{y}}={\displaystyle \frac{{\mathit{F}}_{\mathit{y}}}{-{\mathit{\omega }}^2\mathrm{(}\mathit{M}+\mathit{\mu }\mathrm{)}-\mathrm{i}\mathit{\omega }\mathrm{(}\mathit{\lambda }+{\mathit{\lambda }}_{\text{PTO}}\mathrm{)}+\mathit{K}}}\text{(25b)}& \hfill & {\bar{\mathit{P}}}_{\mathit{A}}={\displaystyle \frac{1}{2}}{\mathit{\lambda }}_{\text{PTO}}{\mathit{\omega }}^2{\mathit{\zeta }}_{\mathit{y}}^2{\mathit{A}}^2\text{(25c)}& \hfill & {\bar{\mathit{P}}}_{\mathit{A}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\mathit{\lambda }}_{\text{PTO}}{\mathit{\omega }}^2{\mathit{A}}^2{\mathit{F}}_{\mathit{y}}^2}{{\mathrm{(}\mathit{K}-{\mathit{\omega }}^2\mathrm{(}\mathit{M}+\mathit{\mu }\mathrm{)}\mathrm{)}}^2+{\mathrm{(}\mathit{\omega }\mathrm{(}\mathit{\lambda }+{\mathit{\lambda }}_{\text{PTO}}\mathrm{)}\mathrm{)}}^2}}.\end{array}$$

The optimal damping coefficient ${\mathit{\lambda }}_{\mathrm{opt}}$ of the generator is then found from the expression obtained after taking the derivative of the average power with respect to ${\mathit{\lambda }}_{\text{PTO}}$ and finding its root(s) $$\begin{array}{}\text{(26a)}& \hfill & {\displaystyle \frac{\partial {\bar{\mathit{P}}}_{\mathit{A}}}{\partial {\mathit{\lambda }}_{\text{PTO}}}}=0\text{(26b)}& \hfill & {\displaystyle \frac{1}{2}}{\displaystyle \frac{{\mathit{\omega }}^2{\mathit{\zeta }}_{\mathit{y}}^2{\mathit{F}}_{\mathit{y}}^2\mathrm{(}{\mathrm{(}\mathit{K}-{\mathit{\omega }}^2\mathrm{(}\mathit{M}+\mathit{\mu }\mathrm{)}\mathrm{)}}^2+{\mathit{\omega }}^2\mathrm{(}{\mathit{\lambda }}^2-{\mathit{\lambda }}_{\text{PTO}}^2\mathrm{)}\mathrm{)}}{{\mathrm{(}{\mathrm{(}\mathit{K}-{\mathit{\omega }}^2\mathrm{(}\mathit{M}+\mathit{\mu }\mathrm{)}\mathrm{)}}^2+{\mathrm{(}\mathit{\omega }\mathrm{(}\mathit{\lambda }+{\mathit{\lambda }}_{\text{PTO}}\mathrm{)}\mathrm{)}}^2\mathrm{)}}^2}}=0\text{(26c)}& \hfill & {\mathit{\lambda }}_{\text{opt}}={\mathit{\lambda }}_{\text{PTO}}=\sqrt{\mathrm{(}\mathit{K}\mathrm{/}\mathit{\omega }-\mathit{\omega }{\mathrm{(}\mathit{M}+\mathit{\mu }\mathrm{)}}^2+{\mathit{\lambda }}^2}\end{array}$$

2.7. Conservation of energy

Conservation of energy is used to verify the simulations in Section 4. As energy enters and leaves the domain through the radiation boundaries ${\mathit{S}}_{\mathit{R}}$ in Fig. 1, it is convenient to define dimensionless numbers in terms of the directions of the wave components encountered at those boundaries: the incident wave component, the reflected wave component and the transmitted wave component. Reflection is the combination of diffracted and radiated waves that propagate in opposite direction to the incident wave. Transmission is the combination of incident, diffracted and radiated waves that propagate in the direction of the incident wave. In equations $$\begin{array}{}\text{(27a)}& \hfill & \mathit{R}={\displaystyle \frac{{\mathrm{\Phi }}_{\mathit{D}}-{\mathit{\zeta }}_{\mathit{y}}{\mathrm{\Phi }}_{\mathit{R}}}{{\mathrm{\Phi }}_{\mathit{I}}}}\text{(27b)}& \hfill & \mathit{T}={\displaystyle \frac{{\mathrm{\Phi }}_{\mathit{D}}+{\mathrm{\Phi }}_{\mathit{I}}-{\mathit{\zeta }}_{\mathit{y}}{\mathrm{\Phi }}_{\mathit{R}}}{{\mathrm{\Phi }}_{\mathit{I}}}}.\end{array}$$

Energy is conserved when the energy flux, or power, that enters the domain, is equal to the power that leaves the domain plus the power that is absorbed by the generator $$\begin{array}{}\text{(28)}& \hfill {\mathit{P}}_{\mathit{I}}={\mathit{P}}_{\mathit{R}}+{\mathit{P}}_{\mathit{T}}+{\mathit{P}}_{\mathit{A}}.& \end{array}$$

Because the incident, reflected and transmitted wave power are related to the incident wave amplitude, the reflected wave amplitude and the transmitted wave amplitude as $$\begin{array}{}\text{(29)}& \hfill {\mathit{P}}_{\mathit{I}\mathrm{,}\mathit{R}\mathrm{,}\mathit{T}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathit{\rho }\mathrm{g}{\mathit{A}}_{\mathit{I}\mathrm{,}\mathit{R}\mathrm{,}\mathit{T}}^2\mathit{\omega }}{\mathit{k}}}\mathrm{,}& \end{array}$$ and because $$\begin{array}{}\text{(30)}& \hfill {\mathit{P}}_{\mathit{A}}={\displaystyle \frac{1}{2}}{\mathit{\lambda }}_{\text{PTO}}{\mathit{\omega }}^2{\mathit{\zeta }}_{\mathit{y}}^2{\mathit{A}}_{\mathit{I}}^2\mathrm{,}& \end{array}$$ energy conservation leads to the following relation $$\begin{array}{}\text{(31)}& \hfill {\mathit{R}}^2+{\mathit{T}}^2+{\mathit{\eta }}_{\mathit{h}}=1\mathrm{,}& \end{array}$$ in which ${\mathit{\eta }}_{\mathit{h}}={\mathit{P}}_{\mathit{A}}\mathrm{/}{\mathit{P}}_{\mathit{I}}$ is the efficiency of the power take-off system.

3. Model description

A boundary element method has been adopted to find the hydrodynamic properties of a floating body with an arbitrary geometry and to solve the complete system of motion with both the diffracted as the radiated waves (leading to reflected and transmitted wave components). The choice is made for a method with panels along all boundaries and simple sources instead of Green’s functions, as this will make the extension to second order more feasible. The linear problem is split in a diffraction and a radiation problem, which, when combined, describe the complete equation of motion of the floater. An example of a domain with panels on all boundaries is shown in Fig. 2. The center of each panel features a collocation point.

Fig. 2.

Panel layout.

3.1. Discretization

The approximate solution is achieved by solving a matrix equation for all boundaries simultaneously. Since this is partly a mixed type boundary condition, the potential and the normal velocity must be defined on all panels. The condition is solved for each collocation point (observation point) in the center of each respective panel. Collocation points are points in the fluid along the boundary of the domain, not exactly on the boundary but with a small distance to the boundary to prevent singularities in the integration, see Fig. 3.

Fig. 3.

Panel definitions.

Source panels are used to find the solution to the boundary problem, with sources spread out over the whole panel. The coordinates of the panels are given in the form of a complex vector z $$\begin{array}{}\text{(32)}& \hfill \mathbf{z}=\mathbf{x}+\mathrm{i}\mathbf{y}.& \end{array}$$

Using z, a linear function is created that describes the path s from the start of the panel at ${\mathit{z}}_{\mathit{a}}$ to the end of the panel at ${\mathit{z}}_{\mathit{b}}$ $$\begin{array}{}\text{(33)}& \hfill {\mathbf{z}}_1\mathrm{(}\mathit{s}\mathrm{)}={\mathbf{z}}_{\mathit{a}}+\mathit{s}\mathrm{(}{\mathbf{z}}_{\mathit{b}}-{\mathbf{z}}_{\mathit{a}}\mathrm{)}& \end{array}$$

The influence matrix ${\mathbf{P}}_{[\mathit{m}\mathrm{,}\mathit{n}]}$ describes the influence of the potential of panel m on panel n. The influence matrix ${\mathbf{Q}}_{[\mathit{m}\mathrm{,}\mathit{n}]}$ is used to find the influence of panel m on panel n with regard to the normal velocity of the panel (pointing into the fluid domain) $$\begin{array}{}\text{(34a)}& \hfill & {\mathbf{P}}_{[\mathit{n}\mathrm{,}\mathit{m}]}={{\displaystyle \int }}_{\text{panel}}{\displaystyle \frac{{\mathbf{l}}_{[\mathit{n}]}}{2\mathit{\pi }}}ln\mathrm{(}{\mathbf{z}}_{\mathit{c}[\mathit{m}]}-{\mathbf{z}}_{1[\mathit{n}]}\mathrm{(}\mathit{s}\mathrm{)}\mathrm{)}\mathrm{d}\mathit{s}\text{(34b)}& \hfill & {\mathbf{Q}}_{[\mathit{n}\mathrm{,}\mathit{m}]}={{\displaystyle \int }}_{\text{panel}}{\displaystyle \frac{{\mathbf{l}}_{[\mathit{n}]}}{2\mathit{\pi }}}{\displaystyle \frac{\partial }{\partial {\mathbf{n}}_{\mathit{p}}}}ln\mathrm{(}{\mathbf{z}}_{\mathit{c}[\mathit{m}]}-{\mathbf{z}}_{1[\mathit{n}]}\mathrm{(}\mathit{s}\mathrm{)}\mathrm{)}\mathrm{d}\mathit{s}\end{array}$$ Here, l is a vector containing the lengths of panels.

The influence matrices can pose a challenge. A numerical approach is taken to find the potential value influence matrix ${\mathit{P}}_{[\mathit{m}\mathrm{,}\mathit{n}]}$ and the normal velocity influence matrix ${\mathbf{Q}}_{[\mathit{m}\mathrm{,}\mathit{n}]}$ . The matrices are found by numerical integration of sources along each panel. The challenge arises when $\mathit{m}=\mathit{n}$ . At this panel, the collocation point ${\mathit{z}}_{\mathit{c}}$ is very close to the panel surface. If the integration step is not small enough, integration errors will occur on this panel as the distance between the panel and the collocation point is critical in the result. To ensure a good result, the integration step size has been reduced until the solution converges.

Linear boundary condition matrices

With these descriptions of the normal velocities and potential in the entire domain, the boundary conditions can be written in a discretized form. To solve for all boundaries simultaneously, the problem must be brought into a single matrix equation. $$\begin{array}{}\text{(35)}& \hfill {\mathbf{C}}_{[\mathit{n}\mathrm{,}\mathit{m}]}^{\mathrm{(}1\mathrm{)}}=\left\{\begin{array}{cc}{\mathbf{Q}}_{[\mathit{n}\mathrm{,}\mathit{m}]}-\mathrm{i}\mathit{k}{\mathbf{P}}_{[\mathit{n}\mathrm{,}\mathit{m}]}\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{L}}\hfill \\ {\mathbf{Q}}_{[\mathit{n}\mathrm{,}\mathit{m}]}\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{b}}\hfill \\ {\mathbf{Q}}_{[\mathit{n}\mathrm{,}\mathit{m}]}+\mathit{k}{\mathit{P}}_{[\mathit{n}\mathrm{,}\mathit{m}]}\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{F}}\hfill \\ {\mathbf{Q}}_{[\mathit{n}\mathrm{,}\mathit{m}]}-\mathrm{i}\mathit{k}{\mathit{P}}_{[\mathit{n}\mathrm{,}\mathit{m}]}\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{R}}\hfill \\ {\mathbf{Q}}_{[\mathit{n}\mathrm{,}\mathit{m}]}\mathrm{,}\hfill & \mathit{n}\in {\mathit{B}}_0\hfill \end{array}\right.& \end{array}$$

This matrix C contains all influence matrices that need to be solved for. Another matrix, describing the known behaviour of the fluid at its boundaries is also determined. This matrix is different for the diffraction problem and the radiation problem.

Diffraction problem: $$\begin{array}{}\text{(36)}& \hfill {\mathbf{K}}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}=\left\{\begin{array}{cc}-2\mathrm{i}\mathit{k}{\mathbf{\Phi }}^{\mathrm{(}1\mathrm{)}}\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{L}}\hfill \\ 0\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{b}}\hfill \\ 0\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{F}}\hfill \\ 0\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{R}}\hfill \\ 0\mathrm{,}\hfill & \mathit{n}\in {\mathit{B}}_0\hfill \end{array}\right.& \end{array}$$

Radiation problem: $$\begin{array}{}\text{(37)}& \hfill {\mathbf{K}}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}=\left\{\begin{array}{cc}0\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{L}}\hfill \\ {\dot{\mathbf{a}}}_{\mathit{n}[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{b}}\hfill \\ 0\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{F}}\hfill \\ 0\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{R}}\hfill \\ 0\mathrm{,}\hfill & \mathit{n}\in {\mathit{B}}_0\hfill \end{array}\right.& \end{array}$$ with ${\dot{\mathbf{a}}}_{\mathit{n}[\mathit{n}]}$ the normal body velocities $$\begin{array}{}\text{(38)}& \hfill {\dot{\mathbf{a}}}_{\mathit{n}[\mathit{n}]}={\mathbf{n}}_{\mathit{b}[\mathit{n}]}·{\mathbf{u}}_{\mathit{b}[\mathit{n}]}.& \end{array}$$

The source strengths ${\mathbf{q}}_{[\mathit{n}]}$ are found as the solution of $$\begin{array}{}\text{(39)}& \hfill {\mathbf{C}}_{[\mathit{n}\mathrm{,}\mathit{m}]}^{\mathrm{(}1\mathrm{)}}{\mathbf{q}}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}={\mathit{K}}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}.& \end{array}$$

And from the resulting source strengths, the potential and the free surface values are reconstructed $$\begin{array}{}\text{(40a)}& \hfill & {\mathbf{\Phi }}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}={\mathbf{P}}_{[\mathit{n}\mathrm{,}\mathit{m}]}^{\mathrm{(}1\mathrm{)}}{\mathbf{q}}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}\text{(40b)}& \hfill & {\mathit{\eta }}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}={\displaystyle \frac{-\mathrm{i}\mathit{\omega }}{\mathit{g}}}{\mathbf{\Phi }}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}\text{for }\mathit{n}\in {\mathit{S}}_{\mathit{F}}\end{array}$$

Second-order boundary condition matrices

To find the second-order results, a new influence matrix G is created. This matrix replaces the ${\mathbf{Q}}_{[\mathit{m}\mathrm{,}\mathit{n}]}$ matrix in finding the derivatives of the potential in y. $$\begin{array}{}\text{(41)}& \hfill {\mathbf{G}}_{[\mathit{n}\mathrm{,}\mathit{m}]}={\displaystyle \frac{\partial }{\partial \mathit{y}}}{\mathbf{P}}_{[\mathit{n}\mathrm{,}\mathit{m}]}=2\mathit{k}tanh\mathrm{(}2\mathit{k}\mathrm{(}\mathit{h}+\mathit{y}\mathrm{)}\mathrm{)}{\mathbf{P}}_{[\mathit{n}\mathrm{,}\mathit{m}]}& \end{array}$$

The $tanh\mathrm{(}2\mathit{k}\mathrm{(}\mathit{h}+\mathit{y}\mathrm{)}\mathrm{)}$ term originates from the transformation of the $cosh\mathrm{(}2\mathit{k}\mathrm{(}\mathit{h}+\mathit{y}\mathrm{)}\mathrm{)}$ to a $sinh\mathrm{(}2\mathit{k}\mathrm{(}\mathit{h}+\mathit{y}\mathrm{)}\mathrm{)}$ when taking the partial derivative in y. This matrix can only be used for a propagating wave without a floating body in the domain, because this use of the derivative is not valid for evanescent waves near the floating body. The complete influence matrix becomes $$\begin{array}{}\text{(42)}& \hfill {\mathbf{C}}_{[\mathit{n}\mathrm{,}\mathit{m}]}^{\mathrm{(}2\mathrm{)}}=\left\{\begin{array}{cc}\mathrm{g}{\mathbf{G}}_{[\mathit{n}\mathrm{,}\mathit{m}]}-4{\mathit{\omega }}^2{\mathbf{P}}_{[\mathit{n}\mathrm{,}\mathit{m}]}\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{F}}\hfill \\ {\mathbf{G}}_{[\mathit{n}\mathrm{,}\mathit{m}]}\mathrm{,}\hfill & \mathit{n}\in {\mathit{B}}_0\hfill \end{array}\right.& \end{array}$$

Since the presented solution does not contain a body, the right-hand side is given by $$\begin{array}{}\text{(43)}& \hfill {\mathbf{K}}_{[\mathit{n}]}^{\mathrm{(}2\mathrm{)}}=\left\{\begin{array}{cc}{\text{QFF}}_{[\mathit{n}]}\mathrm{,}\hfill & \mathit{n}\in {\mathit{S}}_{\mathit{F}}\hfill \\ 0\mathrm{,}\hfill & \mathit{n}\in {\mathit{B}}_0\mathrm{,}\hfill \end{array}\right.& \end{array}$$ after which the source strengths are solved from $$\begin{array}{}\text{(44)}& \hfill {\mathbf{C}}_{[\mathit{n}\mathrm{,}\mathit{m}]}^{\mathrm{(}2\mathrm{)}}{\mathbf{q}}_{[\mathit{n}]}^{\mathrm{(}2\mathrm{)}}={\mathbf{K}}_{[\mathit{n}]}^{\mathrm{(}2\mathrm{)}}.& \end{array}$$

4. Results

The implemention of the method is verified by means of first-order simulations with a floating cylinder featuring a power take-off device (PTO) and second-order simulations of a propagating free surface wave.

4.1. First-order results with a floating body

The vertical motion of a floating cylinder with a circular cross section and radius $\mathit{R}=1$ m in waves is considered. The domain has a size of 25m on either side of the cylinder and a depth of 100 m. The panel size was taken to be 0.52 m, so that twelve panels fit the circumference of the cylinder and 100 panels fit the free surface in horizontal direction on either side of the cylinder. In depth direction, the size of the panels grows exponentially with each panel in downward direction being 2% larger than the panel above it. Diffraction problems with multiple wave frequencies are considered with unit wave height. Radiation problems are simulated for those same frequencies with unit vertical velocity of the floater. Figure 4 shows free surface of the diffraction simulation and the radiation simulation for a wave propagating to the right with a wave frequency of 2.3 rad/s. The diffraction problem features the incident wave component (propagating to the right) and the reflected (diffracted) wave component (propagating to the left) on the left of the cylinder, forming a partial standing wave system. Without a functional GABC at the left boundary of the domain, the wave height to the left of the cylinder would not have a finite value. The right side of the cylinder features a transmitted wave component, which is the combination of the incident wave component and the diffracted wave component. In the radiation problem, the wave components on either side propagate away from the heaving cylinder. The wave component on the left of the cylinder adds to the reflection, the wave component on the right adds to the transmission.

Fig. 4.

Wave profiles with body for the diffraction problem and for the radiation problem.

Diffraction force

The vertical force on the cylinder is found from the discrete potential that is reconstructed from the vector of source strengths $$\begin{array}{}\text{(45)}& \hfill {\mathit{F}}_{\mathit{y}}=\mathrm{i}\mathit{\rho }\mathit{\omega }{\displaystyle \sum }_{\mathit{n}}{\mathbf{\Phi }}_{[\mathit{n}]}^{\mathrm{(}1\mathrm{)}}{\mathbf{s}}_{[\mathit{n}]}{\mathbf{n}}_{\mathit{b}\mathrm{,}\mathit{y}[\mathit{n}]}\text{for }\mathit{n}\in {\mathit{S}}_{\mathit{b}}\mathrm{,}& \end{array}$$ in which ${\mathbf{s}}_{[\mathit{n}]}$ is the vector with panel sizes and ${\mathbf{n}}_{\mathit{b}\mathrm{,}\mathit{y}[\mathit{n}]}\mathit{j}$ the projection of the panel’s normal vector on the y-axis.

For the diffraction problem, the vertical force is shown in Fig. 5. The panel on the left of the figure shows the real part $Re\mathrm{(}{\mathit{F}}_{\mathit{y}}\mathrm{)}$ and imaginary part $Im\mathrm{(}{\mathit{F}}_{\mathit{y}}\mathrm{)}$ of the vertical force. The panel on the right of the figure shows the force amplitude $|{\mathit{F}}_{\mathit{y}}|$ and the phase ${tan}^{-1}\mathrm{(}Im\mathrm{(}{\mathit{F}}_{\mathit{y}}\mathrm{)}\mathrm{/}Re\mathrm{(}{\mathit{F}}_{\mathit{y}}\mathrm{)}$ of the force with respect to the wave.

Fig. 5.

Diffraction force on a cylinder in heave.

Radiation force

For the radiation problem, the vertical force is determined in the same way as in Eq. (45). From that force, the added mass coefficient μ is found from the part of the force in phase with the motion $$\begin{array}{}\text{(46)}& \hfill \mathit{\mu }=-{\displaystyle \frac{Re\mathrm{(}{\mathit{F}}_{\mathit{y}}\mathrm{)}}{{\mathit{\omega }}^2\hat{\mathit{a}}}}.& \end{array}$$

The damping coefficient λ is found as the part of the vertical force out of phase with the motion $$\begin{array}{}\text{(47)}& \hfill \mathit{\lambda }={\displaystyle \frac{Im\mathrm{(}{\mathit{F}}_{\mathit{y}}\mathrm{)}}{\mathit{\omega }\hat{\mathit{a}}}}.& \end{array}$$

The added mass and damping for the floating cylinder are shown as a function of frequency in Fig. 6.

Fig. 6.

μ and λ for a cylinder in heave.

When the added mass and damping are made non-dimensional in a way similar to Yeung [16] as ${\mathit{\mu }}_{\mathit{n}\mathit{d}}=\mathit{\mu }\mathrm{/}0.5\mathit{\rho }{\mathit{R}}^2$ and ${\mathit{\lambda }}_{\mathit{n}\mathit{d}}=\mathit{\lambda }\mathrm{/}0.5\mathit{\rho }\mathit{\omega }{\mathit{R}}^2$ , their values can be compared to existing literature. The non-dimensional values for the added mass and damping are shown in Fig. 7, where they are compared with the results of Frank [4].

Fig. 7.

Verification hydrodynamic coefficients.

Equation of motion

The forces from the diffraction problem and the coefficients from the radiation problem are combined with the buoyancy of the floating cylinder in Eq. (21) and the mass $\mathit{M}=0.5\mathit{\rho }\mathit{\pi }{\mathit{R}}^2$ to determine the RAO according to Eq. (25a). At first, a value of ${\mathit{\lambda }}_{\mathrm{PTO}}=0$ is used as if no PTO device is present. The RAO as a function of one over the wave length, made non-dimensional with the diameter of the cylinder, is shown in Fig. 8, where it is compared to the results of Zhao et al. [17] and Isaacson and Baldwin [5]. The differences with Zhao et al. [17] are small. They are just as small as the differences between Zhao et al. [17] and Isaacson and Baldwin [5] and due to discretization errors.

Fig. 8.

RAO results.

With a PTO device, and a PTO damping ${\mathit{\lambda }}_{\mathrm{PTO}}={\mathit{\lambda }}_{\mathrm{opt}}$ , see Eq. (26c), the RAO – now as a function of frequency – becomes as it is depicted in Fig. 9. The additional damping of the PTO reduces the motion response of the floater. In the same figure that contains the RAO with PTO damping, the efficiency ${\mathit{\eta }}_{\mathit{h}}$ of the PTO device is plotted. One can tell that, as expected from the derivation in Eq. (26c), the efficiency peaks near the natural frequency of the floater.

Fig. 9.

RAO and efficiency of ideally damped heaving cylinder.

Energy conservation

The results of the heaving cylinder with ideal damping according to the equations in Section 2.7 are plotted in Fig. 10 and Fig. 11. Figure 10 shows the radiation coefficient R and the transmission coefficient (T) from the simulation with incoming waves, radiation and diffraction. Energy is conserved when ${\mathit{R}}^2+{\mathit{T}}^2=1$ . The figure demonstrates that this is very nearly the case, apart from a region of frequencies near the natural frequency of the floater. Near that region, the radiated waves have the most influence on the reflection coefficient and the transmission coefficient. Because there is a discretization error in the representation of the wave length (and hence phase velocity of the wave), the effect of the discretization error on R, T and ${\mathit{R}}^2+{\mathit{T}}^2$ is largest in the region near the natural frequency of the floater.

Fig. 10.

Radiation and transmission coefficients and energy conservation.

The results with power-takeoff (PTO) device on the floater are shown in Fig. 11. It gives the reflection coefficient, the transmission coefficient, the efficiency of the power-takeoff and the total energy flux to and from the system. As can be seen, hardly any energy is lost to the system, apart from the region of frequencies near the natural frequency of the floater. This verifies the results for first-order simulations with a PTO present.

Fig. 11.

Energy conservation with PTO.

4.2. Second-order results for a propagating wave

Panels were placed along all boundaries of the simulation domain – instead of only along the boundary of the floater as in many BEM implementations (Frank [4]) – so that second-order simulations would be possible. A second-order propagating wave is simulated. The wave height was 2 m, the wave period was 5.6 s. The water depth was chosen to be 100 m. The panel size was 1 m so that approximately 100 panels were used per wave length. The size of the panels in vertical direction starts out at 1 m, but increases exponentially in vertical direction with each panel further down 2% larger than the panel above it. The linear and second order results for a propagating wave generated with a GABC on the left and an ABC on the right are shown in Fig. 12. These are compared to the analytical relations given above. It is shown that the first and second-order components of the total surface elevation are a good match with the analytical solution.

Fig. 12.

First and second-order wave surface elevation.

5. Conclusions

A novel boundary element method is developed with panels at all domain boundaries and a straightforward Green’s function so that extension to second order is feasible. Simulations are performed in the frequency domain. The diffraction problem with incident waves is solved in one part of the method, the radiation problem in another. A generating and absorbing boundary condition for the side walls of the domain was derived and implemented to make the combination of incident and diffraction waves at the same boundary possible. The following conclusions are found:

Added mass and damping found from the radiation problem are nearly the same as those of Frank [4].

The response amplitude operator found from the combination of two simulations, one with incident waves and diffraction, and one with radiation, closely resembles those in earlier literature.

Without a power take-off device, energy is nearly conserved in the combined incoming wave, diffraction and radiation problem. Energy is not conserved perfectly due to discretization errors.

With a power take-off device, tuned for optimal damping at the natural frequency of the floating body, the energy lost to the floating system is nearly equal to the energy taken up by the power-takeoff device. The small difference is the result of discretization errors.

Following the verification with second-order wave propagation, the model is suitable for extension to second order.

In first-order and second-order simulations, the GABC and ABC at the left and right radiation boundaries function according to theory.

Funding

The authors report no funding.

Conflict of Interest

Peter Wellens is an Editorial Board Members of this journal, but was not involved in the peer-review process nor had access to any information regarding its peer-review.

Author contributions

Maarten Gabriel: Conceptualization, Methodology, Formal analysis, Investigation, Visualization, Writing – original draft. Peter Wellens: Conceptualization, Writing – review and editing, Supervision