Extension of the frequency-domain pFFT method for wave structure interaction in finite depth

To analyze wave interaction with a large scale body in the frequency domain, a precorrected Fast Fourier Transform (pFFT) method has been proposed for infinite depth problems with the deep water Green function, as it can form a matrix with Toeplitz and Hankel properties. In this paper, a method is proposed to decompose the finite depth Green function into two terms, which can form matrices with the Toeplitz and a Hankel properties respectively. Then, a pFFT method for finite depth problems is developed. Based on the pFFT method, a numerical code pFFT-HOBEM is developed with the discretization of high order elements. The model is validated, and examinations on the computing efficiency and memory requirement of the new method have also been carried out. It shows that the new method has the same advantages as that for infinite depth.


Introduction
Boundary element method (BEM) has been widely used in the analysis of wave-structure interaction. Compared with other numerical methods, such as the finite element method (FEM) or the finite volume method (FVM), it has the advantage to reduce the number of dimensions of problem domain by one and lower the difficulty in meshing and requirement of storage as well. When using the free-surface Green function, meshing is only required on the body surface, making further improvement in meshing and storage saving. Therefore, in the linear frequency domain analysis, the BEM with the free-surface Green function has become a good choice.
However, the matrix equation of the BEM is a dense system. To solve the equation by Gauss elimination method, the computational time is and the memory space required to store the whole matrix is , where N is the number of unknowns. If an iterative method is employed, the computational time can be reduced to . Both those two methods lead to a difficulty that requires large memory space and powerful computational capacity when handling large-scale problems.
A couple of accelerated methods have been developed to reduce the computational time and memory requirement for the traditional BEM based on iterative methods. The common idea is to systematically approximate the influence of the kernel between source and field points which are sufficiently far apart relative to the panel size, and to reduce the time and required memory for matrix-vector products. One of the suitable approaches for such processes, where the kernel is oscillatory, is the precorrected Fast Fourier Transform (pFFT) method (Newman and Lee, 2002). In this method, the input geometry is overlaid with a uniform rectangular grid, and the far-field influence of the distributed singularities on the body surface is represented by the influence from the point singularities on the grids. The representation allows the FFT to be involved for the matrix-vector production with considerable efficiency in evaluating the far-field influence accurately. The computational effort for the product is , where N g is the number of grid nodes. This implies a reduction in the computational cost in comparison with traditional BEM especially when N increases to a significant high level, because N g is smaller than N essentially. Furthermore, there is no need to store the matrix during the whole process in the pFFT method, saving large quantities of memory space.
The method was first employed in the field of electrostatic and electro-dynamic analysis with Rankine source by Phillips and White (1997), and was subsequently extended to wave-structure analysis with the free surface Green function by Korsmeyer et al. (1999), Kring et al. (2000) and Dai et al. (2004). The above works are on the basis of a constant panel method. Jiang et al. (2012) developed the pFFT method with higher-order boundary element method, not only obtaining the same or higher accuracy but also further reducing computational time and required memory.
However, so far, previous research of the pFFT method for wave interaction with bodies has been only for the case of infinite depth, because the two parts of the free-surface Green function in infinite depth, including Rankine source and an infinite integral, can be separately used in the calculation of three-dimensional discrete convolution accelerated by the FFT. These two parts will form Toeplitz or Hankel matrix in each direction for matrix-vector product. But for the case of finite depth, the infinite integral itself does not suit the features of neither Toeplitz nor Hankel matrix. In this paper, a decomposition method for the free surface Green function in finite water is proposed to overcome this difficulty, and a pFFT method for finite depth is established. An effective method is also given for the computation of the Green function components numerically, and then a numerical code is implemented with the discretization of higher order elements.
In Section 2, the mathematical model of the linear frequency domain analysis for wave-structure interaction is briefly given, in which the pFFT method is applied to accelerate the computation on the basis of that. In Section 3, we focus on the method of decomposing the two parts of the free-surface Green function in finite depth in the step of convolution. Except that, other steps of the pFFT method are introduced briefly. We adopt the same idea of far-field approximation and near-field correction in finite depth so that they are implemented uniformly with the case of infinite depth. In Section 4, the proposed method is employed in the computation of a truncated cylinder and a hemisphere in different water depths, which are compared with an existing HOBEM model to demonstrate the accuracy of the pFFT-HOBEM model and the effectiveness of the decomposition. Efficiency and memory required are also examined in this section.

Φ
Assuming that the fluid is incompressible and inviscid and the flow is irrotational, the fluid velocity can be represented as the gradient of the velocity potential which satisfies the Laplace equation: For linear problems in the frequency domain, the velocity potential can be expressed in a complex spatial form as: ω where is the wave frequency. The complex spatial poten-ϕ (x, y, z) tial satisfies the Laplace equation in the fluid domain Ω and the linearized boundary conditions on the free surface, the body surface and the plain sea bed: v n where d is the water depth, g, the gravitational acceleration, and , the normal velocity of a point on the body surface.
ϕ (x, y, z) For the problem of wave interaction with a structure, the velocity potential can be decomposed into three components: where is the incident potential, , the diffraction potential, and , the radiation potential generated from the body motion. Except for the incident waves, all the other wave components satisfy the Sommerfeld condition in the far field. ϕ r Associated with motions of a body with six degrees of freedom, the radiation potential can be expressed as follows: where refer to the amplitudes of the body translation, refer to the amplitudes of body rotation. is defined as the radiation potentials generated by unit body motion in the j-th freedom. By combining the unit radiation potential with boundary condition (5), the boundary condition on the body surface can be represented in a uniform way as follows: ∂ϕ where is the unit normal vector on the body surface which is positive when it points outwards the fluid domain, , and is the position vector of the rotation center. Furthermore, note that and ϕ 7 = ϕ d , and the body boundary conditions can be put together in the following form: Substituting the above velocity potentials and the free surface Green function , which satisfies the boundary conditions on the plain sea bed, the linearized free surface and the Sommerfeld condition into the second Green's theorem, we can obtain the boundary integral equation as follows after application of the above boundary conditions: where α is the solid angle coefficient associated with the geometry of body surface. The Green function in finite depth was given by John (1950): where is the horizontal distance between a field point and a source point, is the wave number in infinite depth, and the integration tour is around the singularity underneath. For infinite depth, the Green function can be written as follows: (13) This integral equation can be solved by a boundary element method (Teng and Taylor, 1995) to obtain the potentials at the nodes on the body surface. Then, the exciting forces, added mass and radiation damping can be calculated by integration of the dynamic pressure over the body surface as follows:

Precorrected fast Fourier transform method in finite depth
After discretization of Eq. (11) by BEM, a set of linear algebraic equations can be formed in the following form: where N is the number of nodes for higher order boundary element method (HOBEM).
To set up the matrix A is a process of operations, and it needs memory space for storage of matrix A as the matrix is dense. To solve this equation, it requires operations when a direct method is employed, and operations when an indirect method is employed. Therefore, when tackling a very large scale problem, such as the analysis for wave interaction with a Very Large Floating Body (VLFB) or a Mobile Offshore Bases (MOB), the capacity of a normal computer may not meet the demand for memory space and computational time even though an iterative method is employed.
The pFFT method can reduce memory requirement to and reduce the computation operation to , where N g is the number of grids which is smaller than N. In the pFFT method the computation of potentials due to the source distribution over the body surface is substituted by the influence due to point singularities on the grids and computed by a FFT-accelerated convolution, except the near-field influence which is evaluated by the conventional BEM directly. In summary, the pFFT method can improve computation efficiency and save required memory with little accuracy loss.
To implement the process mentioned above, the pFFT method should be divided into a series of steps: grid set-up, projection, convolution, interpolation and precorrection, which are shown in Figs. 2 and 3. Each step is labeled by a  serial number, and the shadow area is the domain where the precorrection should be carried out when a source point is in the center cube. Note that we adopt the same procedure for finite water depth as the case for infinite depth. Next, we give a brief description of all other procedures but focus on the decomposition method for finite water depth in the step of convolution.
(1) Grid set-up Overlay the whole body geometry with a uniform rectangular grid, and an array of small cubes are obtained with amount of where M is the division in each direction. These small cubes are referred as cells. Then all elements on the body surface are sorted into cells in which they are located. Point singularities are set at half the spacing of the vertices, and an array of 3×3×3 point singularities in each cell is applied to obtain high accuracy in the matrix-vector product. We should mention that the size of the panel elements is associated with the geometry of the body surface and the wave length of the incident waves while the arrangement of point singularities is not. The whole process is represented in Fig. 2. (2) Projection Based on the established grid, an equivalent set of point singularities can be obtained from the set of the element singularity distributions with the projection operators. These are matrices deduced from a collocation problem for each cell that matches the potential due to the grid singularities with the potential due to the singularity distributions on the element at a set of test points on a sphere which encases the entire problem domain. In this step, the operator is based on the fundamental solution of the Laplace equation, but not the free-surface Green function (Kring et al., 2000). The step can be represented in Fig. 4. (3) Convolution This step is to compute the potential at a grid point due to the free-surface Green functions at the other grid points by FFT-accelerated convolution. One of the key principles for applying the pFFT method is whether the matrix-vector product is Toeplitz or Hankel system or not.
In detail, the matrix-vector product can be expressed as an equation as follows: where is the notation of the field points and the source points, q is the strength of singularity at the point, represents a function with the form of independent variable that forms Toepltz matrix along the x, y and z directions, finally leading to a triply nested Toeplitz matrix. Similarly, represents a function that forms the Toepltz matrix along the x and y directions while forming the Hankel matrix along the z direction, which is also a triply nested matrix.
For infinite depth, the Green function is defined in Eq. (13) and the two components are displayed as follows: For the finite depth Green function defined in Eq. (12), the first two terms, Rankine source and its image about the sea bed can be represented as and respectively, but the third term of the infinite integral about z does not meet any features of the two types of matrices. To tackle this difficulty, the free-surface Green function for finite depth can be rewritten in the following form: by substituting the multiplication of hyperbolic cosines. Thus, the Green function can be written in the form of: Notice that the function is an even function of variable t.
has the same feature of that will form the Toepltz matrixes along the x, y and z axes, and has the same feature of that will form the Toepltz matrixes along the x and y axes while form the Hankel matrix along the z axis. For the BEM model, an efficient numerical code for the free surface Green function must be set up in advance TENG Bin, SONG Zhi-jie China Ocean Eng., 2017, Vol. 31, No. 3, P. 322-329 325 F (R, z − ζ) (Newman, 1985). This is the same for the new components, and . Here, a more easily operable way is proposed for the computation of and by application of the established code for the Green Function G, instead of developing new numerical codes for them again. The two components can be computed by the following equations: In terms of the number of singularities N g , there are unknowns in total, N g unknowns in the Toepltz matrix and unknowns in the Hankel matrix. The number of multiplications for the decomposition is , which is quite a small cost because it is only implemented once in the step of projection and can be utilized in every circulation of iterative steps afterwards.
It should be noticed that the above equations do not work well all the time, and special treatment must be done. One case is when a field point and a source point are at the same position, the Green function cannot be calculated. For this case, we may give the Green function component G T or G H an arbitrary value, such as 0.0, temporarily, and substitute it later by the result from the near field direct calculation in the step of pre-correction. Another case is that when the two points are in close proximity, it will lead to the inaccuracy in the calculation of the new Green function components. For this case, the same procedure as for the two points coinciding will be applied to deal with it.
(4) Interpolation This step is to obtain the potential at each node in the elements of the body surface from the potentials at the grid singularities according to test points on the sphere, where the interpolation operators are essentially the transpose of the projection operators. Note that interpolation is for nodes in an element rather than element distributions. The operator is also based on the Laplace equation, which is actually the transpose of the projection operator for nodes. (

5) Precorrection
The approximation by the FFT on the grid is not accurate enough for near-field influences. In this step, one cell and the others sharing vertices with it are considered as the domain of precorrection, where the influences of elements in this domain obtained from the grid-based approximation should be subtracted and then directly computed by the BEM.

Numerical results
Based on the above processes, a pFFT-HOBEM program is developed for wave interaction with bodies in finite water depth. The program is validated, and the examination on its efficiency is also carried out. In this program, the GCR iterative method is employed to resolve the linear al-gebraic equations. All the computations are carried out in a computer with Inter core i7-4790 (3.6 GHz and 4 cores) and 16 GB of RAM.

Validation
(1) Truncated cylinder A truncated cylinder of radius a=1 and draft T/a=1 in water depths of d/a=3, 4 and 5 is computed to validate the present pFFT method. The rotation center is defined at the point (0, 0, 0). The mesh scheme is shown in Fig. 5, in which there are 50 divisions in the azimuthal direction, 10 divisions in the radial direction on the bottom and 10 divisions in the vertical direction. The total number of elements on the body surface is 50×10+50×10=1000, and the number of nodes is 3001. The non-dimensional wave number range is and the discreteness is . The comparisons on the surge force and the heave added mass by the present pFFT-HOBEM and a conventional HOBEM with direct solvers are shown in Figs. 6-8. It can be seen that the results from the two methods agree very well at all computation frequencies. It indicates that the present pFFT-HOBEM model is correct and accurate.
(2) Hemisphere For further validation, a hemisphere of radius a=1.0 in water depths of d/a=2.0, 3.0 and 5.0 is considered. The rotation center is defined at the point (0, 0, 0). A mesh with 30 divisions in the horizontal azimuthal direction and 30 divisions in the vertical circular direction is applied in the computation as shown in Fig. 9. Thus, the total amount of elements is , and the number of nodes is 2701. The comparisons on the surge force and the heave added mass by the present pFFT-HOBEM and a conventional HOBEM with direct solvers are shown in Figs. 10-12. It can be seen that the results from the both models agree very well. It proves again that the present pFFT-HOBEM model is correct and accurate.

Efficiency and memory usage
To evaluate the efficiency and the usage of memory space of the present pFFT-HOBEM model, a couple of cases with different amount of nodes on the body surface are calculated. The truncated cylinder in Section 4.1 is again chosen for the examination. In the examination the wave ka = 1.0 number and the water depth are kept as and d/a = 3.0 for all the cases. The comparison on the computational time is shown in Fig. 13.
It can be seen that with the increase of the node number, the computational time of the pFFT-HOBEM increases more slowly than those of the conventional HOBEMs with direct and iterative solvers. When the number is smaller than 3000, the pFFT-HOBEM and the conventional HOBEM models are at the same level on the computational time. However, when the number of nodes is larger than 3000, the pFFT-HOBEM model uses much less computational time than those of the conventional HOBEM's. With 9000 nodes, the pFFT-HOBEM only uses less than 50% time of the direct-HOBEM and 70% time of the iterative-HOBEM. Therefore, the pFFT is especially suitable for the analysis with large quantity of unknowns.    TENG Bin, SONG Zhi-jie China Ocean Eng., 2017, Vol. 31, No. 3, P. 322-329 327 The comparison on the memory usages is shown in Table 2. It can be seen that the memory usage of the pFFT-HOBEM increases almost linearly with the increase of the node number, while that of the conventional HOBEM increases almost in a quadratic way. When the node number is 1501, the memory usage of the pFFT-HOBEM is about 74.23% of that of the conventional HOBEM; When the node number is 9001, the memory usage of the pFFT-HOBEM is about 11.21% of that of the conventional HOBEM. It means that in handling a problem with a very large quantity of unknowns, the pFFT is a good choice to be used rather than the conventional HOBEM.