A Fully Nonlinear HOBEM with the Domain Decomposition Method for Simulation of Wave Propagation and Diffraction

A higher-order boundary element method (HOBEM) for simulating the fully nonlinear regular wave propagation and diffraction around a fixed vertical circular cylinder is investigated. The domain decomposition method with continuity conditions enforced on the interfaces between the adjacent sub-domains is implemented for reducing the computational cost. By adjusting the algorithm of iterative procedure on the interfaces, four types of coupling strategies are established, that is, Dirchlet/Dirchlet-Neumman/Neumman (D/D-N/N), Dirchlet-Neumman (D-N), Neumman-Dirchlet (N-D) and Mixed Dirchlet-Neumman/Neumman-Dirchlet (Mixed D-N/N-D). Numerical simulations indicate that the domain decomposition methods can provide accurate results compared with that of the single domain method. According to the comparisons of computational efficiency, the D/D-N/N coupling strategy is recommended for the wave propagation problem. As for the wave-body interaction problem, the Mixed D-N/N-D coupling strategy can obtain the highest computational efficiency.


Introduction
As the offshore industry has advanced towards deeper water and harsher environment, the increasingly efficient and accurate numerical models for highly nonlinear free surface wave motions, as well as their interactions with offshore structures, attract more attention. One of the most widely used time domain methods is the Mixed Eulerian Lagrangian (MEL) time stepping technique (Longuet-Higgins and Cokelet, 1976). Based on this method, the fully nonlinear free surface boundary condition is tracked in the Lagrangian frame, and the flow field equation can be solved in the Eulerian description at every time step. This method has been applied in several fully nonlinear wave propagation and interaction problems. Ma et al. (2001aMa et al. ( , 2001b and Kim et al. (2006) developed a numerical wave tank based on the finite element method (FEM). The alternative boundary element method (BEM) was used by Boo (2002), Ferrant et al. (2003) and Bai and Eatock Taylor (2009) for simulating the three-dimensional wave diffraction. Hu et al. (2002), Bai and Eatock Taylor (2006) and Zhou et al. (2013) computed the radiated free surface transformation by the three-dimensional surface piercing structures. More complex structures with more extreme hydrodynamic behavior can also be simulated by the BEM above. Bai et al. (2014) performed a fully nonlinear analysis of nearing-trapping phenomenon around four cylinders. Feng and Bai (2015) studied the wave resonance problem in a narrow gap between two barges using fully nonlinear BEM model. Numerical results have been validated according to the comparison with the experimental data in this work, and hence the numerical challenge can be accepted by the fully nonlinear BEM model.
As is well-known, the traditional BEMs for solving the mixed boundary value problem are computationally expensive because the dense matrixes are generated in the linear matrix equation system. For the problem with N boundary element nodes, it requires at least O(N 2 ) operations, which is prohibitive when N exceeds several thousands. It leads to that the computations of the three dimensional nonlinear simulation are very demanding, especially the complex structures and/or large-scale problems. Various attempts have been carried out to reduce the computational cost, including the accelerated algorithms and domain decomposition method. Typical accelerated algorithms include the Fast Multiple method (Teng and Cou, 2006) and Precorrected-FFT method (Jiang et al., 2012;Teng and Song, 2017). Domain decomposition method has been adopted in fully nonlinear wave and structure interaction problem by Bai and Eatock Taylor (2009). The basic idea is to divide the computational region into many sub-domains, and to apply the boundary element method to each individual sub-domain, leading to far fewer unknowns compared with the original region. Thus, the computational cost per sub-domain is much reduced. Very rough estimate is that the computational effort will be decreased from O(N 2 ) to N s ×O((N/N s ) 2 ), where N s is the number of sub-domains. The key issue in the domain decomposition method is to satisfy suitable continuity requirements between adjacent sub-domains. When the wave flume is divided by multi-subdomains region, various coupling strategies will be generated if different matching schemes are utilized in different common surfaces. It in fact will significantly affect the computational efficiency. Zhang et al. (2013Zhang et al. ( , 2014 developed an efficient heterogeneous flow model combining incompressible viscous flow and potential flow models, where a D-N matching scheme (Bjørstad and Widlund, 1986;Bramble et al., 1986) is adopted between the above two sub-domains. Bai and Eatock Taylor (2009) applied a D/D-N/N coupling strategy (de Haas and Zandbergen, 1996) in a fully nonlinear numerical wave flume, in which multi-subdomains dividing can be performed.
This paper is a further extension of the numerical method presented by Bai and Eatock Taylor (2007). The higherorder boundary element method (HOBEM) using the domain decomposition method with various coupling strategies is adopted to enhance the computational efficiency. The computational efficiency of different coupling strategies is compared and emphasized, which is the most important aim of this paper.

Mathematical formulation and numerical algorithm
ϕ In order to simulate the three-dimensional wave-body interaction problem, a right-handed Cartesian coordinate system Oxyz is defined in Fig. 1, where the origin O is fixed on the mean water surface and z-axis points vertically upwards. Based on the potential flow theory, the fluid flow can be formulated in terms of a velocity potential (x, y, z, t) satisfying Laplace's equation (1) Ω in the fluid domain , and is also subject to various boundary conditions on all surfaces S of the fluid domain.
The boundary integral equation for the potential flow problem can be obtained through Green's second identity, where C(x 0 ) is the solid angle at the field point x 0 , and n is the unit normal vector from the source point x. The solid angle C(x 0 ) can be computed by (3)

ξ η
In this work, the boundary integral equation is discretized by quadratic isoparametric elements. After introducing the shape functions N j ( , ) of the higher-order boundary element method (HOBEM), the discretized equation of Eq.
(2) can be written as: where, N p and N n are the numbers of elements on the Dirichlet (free surface) and Neumann (body surface) boundaries, respectively. The terms in the coefficient matrix A and the right-hand side vector B can refer to the pioneering work in Bai and Eatock Taylor (2006) and Bai and Eatock Taylor (2009).

Domain decomposition method
3.1 Coupling strategies In Eq. (4), A is a dense matrix, which requires O(N 3 ) operations by the direct method (such as the LU decomposition). If the discretization of Eq. (4) is sufficiently well conditioned, the iterative method (such as the GMRES procedure) requiring O(N 2 ) operations can be utilized. However, both of the O(N 3 ) and O(N 2 ) operations are time consuming if the number of unknowns is too large. In this study, the domain decomposition method is adopted to reduce the computational cost of fully nonlinear BEM model. The essential idea is to separate the computational domain into For the purpose of illustration, the computational domain with five sub-domains is depicted in Fig. 2, where the interfaces between the sub-domains are defined as . Different to the other boundaries of Domain , there is no prescribed boundary conditions on the interfaces, that is, both of the potential and its normal derivative are not known a priori. Therefore, an initial Dirichlet condition or Neumann condition needs to be imposed on the interfaces to satisfy the solution conditions in the sub-domains. As an example, if the potential on the interface equals the exact solution of the Laplace's equation in the computational domain , the normal derivatives of the potential and in the sub-domains and will be continuous. In the same way, the exact solution of on the interface will also lead to a continuous between the sub-domains and . The stop criterion on the interface for judging the convergence of the iterative solution can be written as follows: where and are the tolerances of Potential and its normal derivative , and n i and n j point out of the sub-domains and , respectively. In numerical simulations, the stop criterion in Eq. (6) needs to be satisfied on all the interfaces in Domain within every time step.
In the present domain decomposition method, the types of initial conditions (Dirichlet or Neumann) and the algorithm of iterative procedure on the interfaces, which are defined as the coupling strategy in this paper, can significantly affect the accuracy and efficiency of the simulations. A pioneering D/D-N/N coupling strategy in Fig. 2a has been established by Bai and Eatock Taylor (2007), which is referred to as Method 1 in what follows. In this work, many more coupling strategies, including D-N, N-D, and Mixed D/N-N/D coupling strategies, are constructed, and the detailed procedures of them are introduced as follows.

D-N coupling strategy
Before the description, two basic concepts adopted in this paper should be clarified, which are matching scheme and coupling strategy. The matching scheme refers to the iterative procedure one interface between two sub-domains, while the coupling strategy is focused on what matching schemes are utilized in all the interfaces. A typical example  could be cited that the D-N matching scheme is adopted on the interface in Fig. 2b, while the D-N coupling strategy means all the interfaces in the computational domain are operated by the D-N matching scheme. In the solution of elliptic problems, the D-N matching scheme has been conducted by Bjørstad and Widlund (1986) and Bramble et al. (1986), successfully. It is adopted in the present fully nonlinear BEM method as follows, where between the subdomains and is discussed.
(2) Solve Laplace's equation with the Dirichlet condition (k) on in . This yields on in .
(3) Formulate a Neumann condition on for as assign the computed solutions from , (5) Calculate the maximum error (k+1) on 12 based on the potentials (k) and (k) in and , (8) ε (6) If (k+1) satisfies the prescribed stop criterion, then exit the loop, and the final approximate solutions on the interface are ) .
(9) (7) Go to Step 2 and repeat the procedure with k = k + 1. As illustrated in this figure, all the interfaces between the sub-domains are operated by D-N matching schemes, correspondingly the D-N coupling strategy is constructed, which is defined as Method 2 in this paper.

N-D coupling strategy
The N-D coupling strategy, which is defined as Method 3 in this work, indicates that all the interfaces between the adjacent sub-domains are operated by N-D matching schemes. As shown in Fig. 2c, the similar process with the D-N matching scheme can be observed in the procedure of N-D matching scheme. The main discrepancy between them is the N-D matching scheme beginning with the initial Neumann condition on in . Laplace's equa- tion with Neumann condition is solved in , and the solution (k) is yielded on the interface . Correspondingly, the Dirichlet condition can be imposed on in as follows: (11) ε n If (k+1) satisfies the prescribed stop criterion, the final approximate solutions on the interface can be computed as follows: λ where the coefficient = 0.5 is adopted in Eq. (12).

Mixed D-N/N-D coupling strategy
The Mixed D/N-N/D coupling strategy, defined as Method 4, equipped with both the D-N and N-D matching schemes in the interfaces between the adjacent sub-domains. A typical example can be observed in Fig. 2d, where the interfaces and are operated by D-N matching scheme, while the interfaces and are operated by N-D matching scheme. The further reason and advantage of this coupling strategy will be discussed in the next sub-section.

Numerical implementations
Meshes of a fixed vertical circular cylinder in a numerical wave flume are generated in Fig. 3, where three sub-domains are involved. In numerical implementations, structured quadrilateral meshes are utilized on the vertical surfaces, and unstructured triangular meshes are adopted on the free surface generated by the Delaunay triangulation method. A double or triple node is employed on the intersections between the free surface and vertical surfaces. According to the multiple-node technique, the intersection lines can be determined conveniently without additional treatment. Only the nodes on the free surface need to be updated and the oth- er multiple nodes on the solid surfaces can be allowed to move with it.
In the D-N and N-D matching schemes discussed above, Laplace's equation with Dirichlet or Neumann conditions on the interfaces need to be solved in the sub-domains at each step. Besides, another Dirichlet condition is specified on the free surface, leading to some double nodes are distributed on the intersection line between the interface and free surface. These double nodes share the same position and the same boundary value, and therefore the corresponding coefficients in Eq. (4) will be exactly the same. This will result in a singular algebraic equation system. The semi-discontinuous elements allocated in the top layers of elements at the interfaces can be adopted to deal with the difficulty. γ = 0.6 Fig. 4 shows the semi-discontinuous element on the top layer of the interface and the adjacent continuous element on the free surface. The double nodes on the intersection line between the free surface and interface are the reason of singularity by the Dirichlet conditions. The semi-discontinuous element approach locates the new collocation points inside the element. According to the different positions of the new collocation nodes x i , the different coefficients in the matrix of Eq. (4) can be generated, and therefore the singularity can be avoided in the algebraic equation system. In numerical simulations, the discontinuity coefficient is adopted in the local intrinsic coordinates for specifying the positions of the new collocation nodes. The relevant modifications, such as the values of Dirichlet conditions and the triangular polar-coordinate transformations on the new collocation points, are needed to be conducted in the operations of coupling strategies. The semi-discontinuous element approach can also be found in Sun et al. (2008) for removing the irregular frequencies and Bai and Eatock Taylor (2007) for performing the D/D-N/N coupling strategy.
Although the semi-discontinuous element approach could avoid the singularity, the condition number of the matrix equation is still a little poor, leading to the slow convergence when using an iterative scheme, such as the GMRES procedure. Therefore, the LU decomposition pro- cedure have to be adopted to obtain the solutions, which is faster than the iterative scheme for the poorly conditioned matrix equations. Numerical simulations suggested that the computational effort is acceptable for the wave propagating problem as the number of unknowns is small in each subdomain. However, it is a little time consuming for the wavebody interaction problem because more elements have to be utilized in the sub-domain with the body. This is the main reason for constructing the Mixed D/N-N/D coupling strategy in this paper. As shown in Fig. 2d, the D-N and N-D matching schemes are adopted on the interfaces and , respectively, leading to Neumann boundary conditions appear on both the vertical common surfaces in the body sub-domain . On this occasion, there is no singularity in the body sub-domain and the GMRES procedure can be used to solve the matrix equation. Consequently, the computational cost can be reduced by adjusting the D-N and N-D coupling scheme in the Mixed D/N-N/D coupling strategy.
In this work, the numerical simulations start from the still state, which means that the zero potential and zero velocity are imposed on all the surfaces, including the interfaces between the sub-domains. The initial values of Dirichlet or Neumann conditions are set to be =0 or on all the interfaces at the first time step. By solving the boundary integral equation and conducting the transfer of information between the adjacent sub-domains, the unknowns vector B in Eq. (4) can be obtained by different coupling strategies. The standard fourth-order Runge-Kutta scheme (RK4) is adopted for updating the geometries and potentials on the free surface. On the interfaces, the obtained potentials or velocities at this time step are utilized as the initial values at the next time step, which can allow the new beginning of the coupling strategies. During the above processes, the stop criterion in Eq. (6) needs to be satisfied on all the interfaces within every time step.

Numerical results
In this section, four types of coupling strategies mentioned above are adopted for simulating the wave propagating in the wave flume with/without a diffracting body, where the computational efficiencies of them are focused.

Wave propagating in a wave flume
ρ In the presentation of present results, the fluid density , gravitational acceleration ɡ and the water depth d are set to unity. The other parameters such as the time and distances are taken to be non-dimensionalized, effectively. As shown in Fig. 3, only half of the numerical wave flume needs to be considered according to the symmetric character. The total length and half-width of the wave flume are taken as 15.7 and 0.31, respectively. A monochromatic wave is generated by a piston-like wave maker situated at x = -7.7 at the left end of the wave flume undergoing the following motion, Fig. 4. Semi-discontinuous element on the top layer of the common surface.
where S 0 (t) and U 0 (t) are the displacement and velocity of the piston-like wave maker, respectively; and are the motion amplitude and frequency, respectively. A numerical zone based on Bai and Eatock Taylor (2006) is situated for avoiding the wave reflection near the right hand end of the wave flume, which starts at a distance of one incident wavelength measured from the right end. In this simulation, the 10 sub-domains of equal size scheme are adopted to divide the whole computational domain.
The convergence of the computation by various coupling strategies is examined. The propagating wave generated by the wave maker moving at and with various stop criterion and in Eqs. (8) and (11) is computed. There are 6 elements in the vertical direction on the vertical surfaces. The 10×3 intervals on the boundary of the free surface are adopted in each subdomain. So that the total elements and nodes are 2200 and 6840, respectively. Note that the mesh convergence adopted here has been validated in Bai and Eatock Taylor (2007). In numerical simulation, 20 time steps are chosen to perform the calculation ε ϕ ε n ε ϕ ε n for each wave period. Fig. 5 shows the time history of the wave elevation at x=-3.12 by various stop criteria and .
The results indicate that the stop criterion =1×10 -4 and =5×10 -4 can give converged calculations, and it can be taken as the baseline of the numerical simulations. Then, the comparisons between different coupling strategies and single domain results are shown in Fig. 6. It can be observed that the domain decomposition methods agree well with those in a single domain, which suggests that the domain decomposition approaches presented here can obtain the accurate results.
Furthermore, the efficiency of the domain decomposition method with various coupling strategies is investigated for the same case of and . The meshes in the single domain is adopted for the basic mesh scheme in numerical simulations, where the equal numbers of elements on the free surface and solid walls are adopted as the cases with various sub-domains. The total numbers of elements and nodes in the whole domain with various sub-domains are different because of the existence of the interfaces. Six kinds of sub-domains numbers are considered for the efficiency test, as shown in Table 1, in which they are  Fig. 7 also suggest that the CPU time tends to increase if more sub-domains are adopted in the present domain decomposition method. The rapid increase of the iterative times is the major reason for the results, in which the transfer of information between the interfaces approach to significantly large when too many sub-domains are utilized. Furthermore, comparing the different coupling strategies, we can see that Method 1 gives rise to the fewest CPU time for the wave propagation problem.

Wave propagating in a wave flume
The wave diffraction around a bottom mounted vertical cylinder in a rectangular wave flume is simulated in this section. The total length and half-width of the numerical wave flume are 11.3 and 0.62, respectively, while the wave maker is situated at x=-7.0. Compared with the wave propagating problem in pervious sub-section, a doubled width of the flume is adopted for reducing the wave reflection by the flume side walls. Fine mesh density near the body surface is used in the simulation, which are 16 elements around one half of the vertical cylinder. As for the other surfaces, largely the same meshes with the wave propagating problem are equipped, including the free surface, solid surfaces and common surfaces.
The time histories of the wave forces on the cylinder in the x-direction are given in Fig. 8. The incident waves are generated by the piston-like wave maker with the amplitude at the frequency . In this simulation, the domain decomposition method is performed based on 9 subdomains dividing. All the numerical results by various coupling strategies can agree well with those obtained in single domain, indicating the accuracy of the present domain decomposition approaches. Table 2 tabulates the mesh resolution for various subdomain schemes. Since the existence of the cylinder, more elements have to be adopted in the body sub-domain. In this table, they are (Elements or nodes on each sub-domains without cylinder)/(Elements or nodes on each sub-domains with cylinder)/(Elements or nodes on total domain). The comparisons of CPU time for these dividing schemes by various coupling strategies are shown in Fig. 9, in which the time for single domain results are chosen for normalization. It can be seen that the normalized CPU time is even larger than 1 for Methods 1-3 when the number of the subdomain is N s = 3. It is due to the LU decomposition procedure has to be used in the body sub-domain, which is time consuming since the unknowns is 2272. With the increase of sub-domains number, the computational time gives rise to small, implying the advantage of domain decomposition method. Moreover, because the GMRES procedure is always adopted in the body subdomain, Method 4, the Mixed D/N-N/D coupling strategy, can give rise to the fewest computer time  in all methods. It is recommended to adoption when the body exists.

Discussion
The domain decomposition method in this work can further cooperate with parallel computation. According to the allocation of the thread to each sub-domain, the MPI parallel computation can be carried out in the domain decomposition method. Theoretically, all of the four types of coupling strategies can extend to parallel computation. In the D-N, N-D, and Mixed D/N-N/D coupling strategies, the boundary integral equation in the downstream sub-domain has to be established after solving Laplace's equation in the upstream sub-domain. For example, in Fig. 2b, the Neumann condition on the interface in the sub-domain has to be obtained according to the solution in the subdomain . In other words, the computation in the sub-domain must be performed after the solution in the sub-domain . The threads allocated to the sub-domains will be executed in order, one after another. It in fact is not the real MPI Parallel because the threads cannot be operated at the same time. An alternative method for performing the real MPI parallel is that an initial guess needs to be imposed on the interfaces in the downstream sub-domains. In the D/D-N/N coupling strategy, the boundary values on the interfaces in the sub-domain is independent of the adjacent subdomain. The threads allocated to the sub-domains can be executed at the same time. Therefore, the D/D-N/N coupling strategy has the most advantage for further extending to MPI parallel.

Conclusion
Numerical simulations for fully nonlinear regular wave propagation and diffraction around a fixed vertical circular cylinder in a wave flume are carried out in this study. The higher-order boundary element methods using the domain decomposition method with various coupling strategies are adopted for increasing the computational efficiency. Numerical investigations suggested that the domain decomposition method can reduce the computational cost, significantly. The D/D-N/N coupling strategy can obtain the fewest CPU time when considering the wave propagation without body problem in the present results. If the wave-body interaction problem is simulated, the mixed D-N/N-D method is recommended for higher efficiency.