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
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,
Conservation of momentum is obtained through the Bernoulli equation
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
Fig. 1.
Table 2
Abbreviation | Physical meaning |
Incident wave | |
Diffracted wave | |
Radiated wave | |
Body surface | |
Free surface | |
Radiation boundary | |
Bottom boundary |
At the bottom boundary
A kinematic condition is solved at free surface
A highlight of this paper is the generation of waves at radiation boundaries
The wave speed c is defined as
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
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
The atmospheric pressure at the free surface is assumed to be zero
Linearizing the Bernoulli equation, Eq. (2), and substituting Eq. (9), leads to the dynamic free surface boundary condition
Linear potential flow theory is also known as Airy wave theory. The potential function of an Airy wave
The characteristic equation of the system of equations is the dispersion relation
It relates the wave frequency and the wave number. Because
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
Here,
The second-order boundary condition (Stokes [10]) is referred to as the quadratic forcing function (QFF)
An index to find the non-linearity for waves in finite water depth is the Ursell number
2.4.Pressures and forces
To find the pressures in the domain use is made of the linearized Bernoulli equation
The forces on the body are the integration of the pressures along the body surface
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
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
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
Velocity
The average power in one wave period T then becomes
Using RAO
The optimal damping coefficient
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
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
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
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.
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.
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
Using z, a linear function is created that describes the path s from the start of the panel at
The influence matrix
The influence matrices can pose a challenge. A numerical approach is taken to find the potential value influence matrix
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.
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:
Radiation problem:
The source strengths
And from the resulting source strengths, the potential and the free surface values are reconstructed
Second-order boundary condition matrices
To find the second-order results, a new influence matrix G is created. This matrix replaces the
The
Since the presented solution does not contain a body, the right-hand side is given by
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
Fig. 4.
Diffraction force
The vertical force on the cylinder is found from the discrete potential that is reconstructed from the vector of source strengths
For the diffraction problem, the vertical force is shown in Fig. 5. The panel on the left of the figure shows the real part
Fig. 5.
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
The damping coefficient λ is found as the part of the vertical force out of phase with the motion
The added mass and damping for the floating cylinder are shown as a function of frequency in Fig. 6.
Fig. 6.
When the added mass and damping are made non-dimensional in a way similar to Yeung [16] as
Fig. 7.
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
Fig. 8.
With a PTO device, and a PTO damping
Fig. 9.
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
Fig. 10.
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.
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.
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:
1. Added mass and damping found from the radiation problem are nearly the same as those of Frank [4].
2. 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.
3. 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.
4. 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.
5. Following the verification with second-order wave propagation, the model is suitable for extension to second order.
6. 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
References
[1] | A. Ballast, Water Waves, Fixed Cylinders and floating Spheres, 2004, p. 302. ISBN 9064648956. |
[2] | R.W. Bos, M.v.d Eijk, J.H. Besten and P.R. Wellens, A reduced order model for FSI of tank walls subject to wave impacts during slosing, International Shipbuilding Progress 69 (2022). |
[3] | W.R. Dean, On the reflexion of surface waves by a submerged plane barrier, Mathematical Proceedings of the Cambridge Philosophical Society 41: (3) ((1948) ), 231–238. doi:10.1017/S030500410002260X. |
[4] | W. Frank, Oscillations of Cylinders in or Below the Free Surface of Deep Fluids, (1967) . |
[5] | M. Isaacson and J. Baldwin, Wave Propagation Past a Pile-Restrained Floating Breakwater 8(4) (1998), 1–5. |
[6] | M. Kim and D.K. Yue, The complete second-order diffraction solution for an axisymmetric body. Part 2. Bichromatic incident waves and body motions, Journal of Fluid Mechanics 211: ((1990) ), 557–593. doi:10.1017/S0022112090001690. |
[7] | Y. Kim, S.H. Kim and T. Lucas, Advanced Panel Method for Ship Wave Inviscid Flow Theory (SWIFT), (1989) . |
[8] | J.N. Newman, Marine Hydrodynamics – Waves and Wave Effects (Ch. 6), Marine hydrodynamics (1977), 237–327. ISBN 978-0-262-14026-3. |
[9] | T.F. Ogilvie, First- and second-order forces on a cylinder submerged under a free surface, Journal of Fluid Mechanics 16: (3) ((1963) ), 451–472. doi:10.1017/S0022112063000896. |
[10] | G.G. Stokes, On the theory of oscilatory waves, Mathematical and Physical Papers 1: ((1847) ), 225–228. ISBN 0521419662. doi:10.1017/CBO9780511702242.016. |
[11] | F. Ursell, On the Heaving Motion of a Circular Cylinder on the Surface of a Fluid, April 1948, 1949. |
[12] | F. Ursell and W.R. Dean, Surface waves on deep water in the presence of a submerged circular cylinder, Mathematical Proceedings of the Cambridge Philosophical Society 46: (01) ((1950) ), 153. doi:10.1017/S0305004100025573. |
[13] | P.R. Wellens and M. Borsboom, A generating and absorbing boundary condition for dispersive waves in detailed simulations of free-surface flow interaction with marine structures, Computers and Fluids (2019), 104387. doi:10.1016/j.compfluid.2019.104387. |
[14] | P. Xu and P.R. Wellens, Fully nonlinear hydroelastic modeling and analytic solution of large-scale floating photovoltaics in waves, Journal of Fluids and Structures 109: ((2022) ), 103446. doi:10.1016/j.jfluidstructs.2021.103446. |
[15] | P. Xu and P.R. Wellens, Theoretical analysis of nonlinear fluid–structure interaction between large-scale polymer offshore floating photovoltaics and waves, Ocean Engineering 249: ((2022) ), 110829. doi:10.1016/j.oceaneng.2022.110829. |
[16] | R.W. Yeung, A singularity-distribution method for free surface flow problems with an oscillating body, 1973. |
[17] | X. Zhao, D. Ning, C. Zhang and H. Kang, Hydrodynamic investigation of an oscillating buoy wave energy converter integrated into a pile-restrained floating breakwater, Energies 10(5) (2017). ISBN 8641184708. doi:10.3390/en10050712. |