Meshless Method for Analysis of Permeable Breakwaters in the Proximity of A Vertical Wall

In the present work, the improved version of the meshless singular boundary method (ISBM) is developed for analyzing the performance of bottom standing submerged permeable breakwaters in regular normally incident waves and in the proximity of a vertical wall. Both single and dual prismatic breakwaters of rectangular and trapezoidal forms are examined. The physical problem is cast in terms of the Laplace equation governing an irrotational flow and incompressible fluid motion with appropriate mixed type boundary conditions, and solved numerically using the ISBM. To model the permeability of the breakwaters fully absorbing boundary conditions are assumed. Numerical results are presented in terms of hydrodynamic quantities of the reflection coefficients. These are firstly validated against the results of a multi-domain boundary element method (BEM) developed independently for a previous study. The agreement between the results of the two methods is excellent. The coefficients of reflection are then computed and discussed for a variety of structural conditions including the breakwaters height, width, spacing, and absorbing permeability. Effects of the proximity of the vertical plane wall are also investigated. The breakwater’s width is found to have only marginal effects compared with its height. Permeability tends to decrease the minimum reflections. These coefficients show periodic variations with the spacing relative to the wavelength. Trapezoidal breakwaters are found to be more cost-effective than the rectangular breakwaters. Dual breakwater systems are confirmed to perform much better than single structures.


Introduction
Submerged breakwaters are popular structures designed to safeguard coastal areas in situations where complete protection from waves is not required. When these structures are constructed far from the shore, they reduce wave energy by reflecting much of the incident waves and hence diminish wave transmission (Lorenzoni et al., 2016). In the situations where submerged breakwaters are built in front of coastal seawalls their primary function is to reduce wave reflections (Koley et al., 2015a). These structures can offer potentially cost-effective solutions. They can effectively create calmer areas in their leeward side (as in harbors and marinas), have the ability to decrease erosion and the sedi-ment transport, and do not completely inhibit fish passage and water circulation (Mangor et al., 2017).
To identify the conditions at which the maximum/minimum reflections would occur (Bragg resonant reflections), wave interactions with submerged breakwaters have been studied both experimentally and theoretically. Experimental studies include those of Dattatri et al. (1978), Hsu et al. (2002), Cho et al. (2004) and Liu et al. (2016). All these studies indicate that Bragg resonant reflections are improved by increasing both the breakwater's height and the number of breakwaters. Analytical solutions based on the eigenfunction expansion method and in the context of linear potential theory were also provided by a number of in-vestigators such as Abul-Azm (1994), Cho et al. (2004), Twu and Liu (2004), Koley et al. (2015a), Liu et al. (2016) and Zhao et al. (2017aZhao et al. ( , 2017b. Numerical studies based on the boundary element method (BEM) in the context on the linear wave theory have also been developed successfully by a number of investigators including Hsu and Wu (1998), Yueh and Chuang (2009), Koley et al. (2015aKoley et al. ( , 2015b and Zhao et al. (2017aZhao et al. ( , 2017b. In recent years, a number of meshless methods have also been exploited to study the transmission of waves over submerged breakwaters, such as the method of fundamental solutions (MFS) (Tsai and Young, 2011), the regularized meshless method (RMM) (Chen et al., 2011;Ouyang et al., 2016), and the singular boundary method (SBM) (Chen et al., 2014;Fu et al., 2015;. The results of these investigations are confirmed by comparison with previous experimental data, and/or analytical and BEM methods. Collocation meshless methods are advantageous as only scattered source nodes without any connectivity are required for the numerical calculations, eliminating the need for a mesh and singular numerical integration over the elements as in the traditional BEM. Collocation meshless methods are broadly classified into two categories, domaintype and boundary-type methods. In the first category source nodes are distributed over the whole domain to seek a solution such as in the method of particular solution (MPS) (Lin et al., 2017a), the localized method of approximate particular solution (LMAPS) (Lin et al., 2017b), and the multiple scale method of particular solutions (MSMPS) (Lin et al., 2018b). In the second category the source nodes are placed only on the domain boundaries to find a solution as in the MFS, the RMM, and the SBM.
The MFS uses single layer potentials as its kernel (basis) functions and a nonphysical (fictitious) boundary other than the physical one to avoid the singularity of the kernels. This is a major drawback of the method. Also for large scale problems the MFS may produce the ill-conditioned coefficient matrix which may lead to inaccurate results upon its solution. Lin et al. (2011) examined the ill-conditioning of the MFS and obtained more stable results by combination with a regularization technique. Lin et al. (2016aLin et al. ( , 2016b also considered the fast solution of dense systems in the MFS for several multi-dimensional wave and Helmholtz problems.  proposed a dual-level MFS and investigated its effectiveness for three-dimensional exterior high frequency acoustic problems. To avoid the idea of using a fictitious boundary in the MFS, Young et al. (2005) proposed the RMM that employs the same physical boundary and double layer potentials as kernel functions with the regularization technique of the subtracting and adding-back methods. However, the use of double layer potentials can affect the whole accuracy of the method. The SBM proposed by Chen and Wang (2010) and Chen and Fu (2010), does not require the fictitious boundary either but uses single layer potentials as its kernel functions and the concept of the origin intensity factors. In the traditional SBM the origin intensity factors are evaluated by an inverse interpolation technique for both the Dirichlet and Neumann boundary equations. To carry out this technique a cluster of sample nodes placed inside the physical domain is required. The solution accuracy may be sensitive to the location of the sample nodes. To overcome the shortcoming of the sample nodes Chen and Gu (2012), Gu et al. (2012) and Gu andChen (2013, 2014) proposed an improved version of the original SBM (ISBM). In this method the disingularisation of the origin intensity factors for the Neumann boundary equation is carried out by the regularization technique of subtracting and adding-back. For the Dirichlet boundary equation the origin intensity factors are evaluated either by an improved inverse interpolation technique that does not need the troublesome sample nodes, or by integration of the fundamental solution on line segments. The SBM has been recently applied to some challenging problems, including wave propagation analysis in periodic structures , seismic wave scattering (Lin et al., 2018a), acoustic wave problems (Li and Chen, 2018), and scattering of SH waves (Tang et al., 2018).
The SBM is a truly meshless boundary type collocation method which does not require a mesh or singular numerical integration over the elements and has good stability. It circumvents the fictitious boundary of the MFS, and has better accuracy than the RMM. The ISBM shares all the merits of the SBM, and in addition it does not require the sample nodes. The aim of the present paper is to develop a numerical model based on the ISBM and non-viscous flow theory to analyze numerically the performance of bottom standing single and dual submerged permeable breakwaters in regular normally incident waves and near a vertical wall. The breakwaters are of rectangular and trapezoidal forms. Permeability of the breakwaters is modeled with absorbing boundary conditions. Effects of several parameters relating to the geometry of the breakwaters including their permeability are examined.

Formulation of the problem
We consider in this study single and dual bottom standing submerged permeable breakwaters of rectangular and trapezoidal shapes in front of a vertical wall as shown in Fig. 1. For the sake of generality the method is developed for a system of dual trapezoidal breakwaters.
The present numerical solution is developed based on the linear (small amplitude) wave theory. The idealized geometry of the two-dimensional problem in a Cartesian system (x-y) is shown in Fig. 2. Regular normal waves of small amplitude a, wavelength L (measured between two consecutive crests) and period T, impinge from the left in the water of depth d. Assuming an irrotational flow and incompressible fluid motion, the problem is formulated using a velocity potential where Re denotes the real part, is the time independent spatial velocity potential, is the wave angular frequency, and t is the time. The wave number is the solution of the dispersion relation , where g is the gravitational acceleration. The wave field is totally specified if the two dimensional velocity potential is known.
The breakwaters are similar with the height h, bottom width w b , and top width w t . For rectangular breakwaters . The distance between the rear breakwater and vertical wall is . The dual breakwaters are separated by a distance measured from the centers of the breakwaters. The total fluid domain is divided in two regions as shown in Fig. 2. Region I at ̶ ∞ is the region where the waves are incoming (inflow). Region II is delimited by the walls of the breakwaters ( , , and for the front breakwater 1, and , , and for the back breakwater 2), the seabed rigid boundary , the free surface boundary , the reflection boundary of the inflow region, and the boundary of the vertical rigid wall. The spatial velocity potential satisfies the following conditions: (1) ∂φ ∂n = 0, y = 0 on the seabed rigid boundary Γ s ; where n is the normal to the boundary pointing out of the flow region.
For the permeability of the breakwaters we adopt absorbing boundary conditions, similar to those of Isaacson and Qu (1990), Chen et al. (2011), , and Zhao et al. (2017b). The front and back walls of the breakwaters (Boundaries and for the front breakwater, and and for the back breakwater) are treated as absorbing boundaries. The top walls of the breakwaters (Boundaries and ) are assumed to be rigid impermeable walls. These conditions are summarized as follows: (m=1,2,3,4) are the corresponding absorbing (permeability) parameters respectively of the boundaries , , and . G is a non-dimensional parameter which corresponds to the permeability of the absorbing layer of the medium covering the breakwater. In practice using absorbent materials (such as rubble material, rocks, interlocking armor units, …) along the breakwater front and/or back sides can partially/or fully dissipate incident wave energy depending on the absorption coefficients G as in a porous breakwater. In the limit when G=0, corresponding to an impermeable rigid wall, all wave energy is reflected back. In the other limit when G=1, corresponding to a fully absorbing wall, all wave energy is fully absorbed and dissipated within the porous medium. Further details can be found in Isaacson and Qu (1990), and Zhao et al. (2017a).
The reflection condition at the inflow region is expressed as: where is the incident velocity potential.
The reflection condition in the infinite strip problem is treated by transferring the far field potentials at a fictitious vertical boundary at a finite distance representing the left boundary of the fluid domain. The analytical expression at this boundary is given by: and A − where is an unknown complex coefficient to be determined. The disturbances are guaranteed to be out-going waves only (see for example Bakhti et al., 2017;Chioukh et al., 2017). The incident velocity potential is defined as: The special matching condition at the interfaces of the inflow region ensures a smooth transfer of the mass flow from one region to the next. Once the potentials are calculated by satisfying the reflection boundary condition of Eq. (6), they are matched to those of Eq. (7), and then the unknown coefficient is evaluated following the method of Yueh and Chuang (2012): where .
The reflection coefficient ( ) is determined from the following expression (see for example Bakhti et al., 2017;Chioukh et al., 2017): 3 Numerical solution by the ISBM For the numerical solution the total boundary of the whole computational domain is discretized as shown in Fig. 3 for single and double breakwaters, respectively.
In the ISBM, the nodal values of the potentials and their fluxes are expressed as the linear combinations of the fundamental solutions and their derivatives (Chen and Gu, 2012;Gu et al., 2012;Gu and Chen, 2013), where are unknown coefficients to be determined, and are respectively the collocation points and the source points , and N is the total number of points. are the essential boundary conditions (Dirichlet), are the natural boundary conditions (Neumann), and is the normal at the collocation point . The coefficients and are the source intensity factors corresponding respectively to the fundamental solution and its derivative.
is the fundamental solution of the 2D Laplace equation. It depends only on the Euclidean distance between the collocation points and the source points , i.e. , and is given together with its normal derivative as: and are the component of the normal at the collocation point .
The coefficients and are the diagonal elements of the ISBM interpolation matrices. They arise when the collocation points and source points coincide ( ). Direct evaluation of these coefficients is unfeasible because of the singularities inherent in the fundamental solution and its derivative. In this study the coefficients are evaluated simply by the integration of the fundamental solution on the line segments leading to simple analytical expression as in Brebbia and Dominguez (1992), and Gu and Chen (2014): q ii For the coefficients , a simple expression is derived by Gu et al. (2012), using a regularization process of subtracting and adding-back to remove singularities: where and are the half distances respectively between the collocation points and , and the source points and .
is the normal at the source point . The boundary conditions given by Eqs.
(2)-(6) are satisfied by a linear combination of Eqs. (11) and (12). The discretization process leads to: for the nodes (the free surface boundary) for the nodes (the reflection boundary at ) for the nodes + + + (the rigid boundaries) x i ∈Γb1 Γ b3 Γ b4 Γ b6 for the nodes + + + (the breakwaters absorbing boundaries) The resulting discretized Eqs. (16)-(19) are written in a more compact matrix form as: where N is the total number of the nodes on the whole domain boundary, e.g., , where , , , and are the numbers of the nodes respectively on the boundaries , , (rigid boundaries + + + ), and (breakwaters absorbing boundaries + + + ). The algebraic system of the equations expressed by Eq. (20) is solved numerically using a Gaussian elimination algorithm to yield the vector of unknowns . The potential and its derivative at the nodes are then computed using Eqs. (11) and (12).

Error analysis and validation of the numerical method
Error sensitivity and convergence of the numerical results with respect to the total number of the boundary collocation nodes are verified for single and dual breakwaters. Errors are calculated relative to the results of a multi-domain boundary element method (BEM) developed independently for another study (Bakhti et al., 2017;Chioukh et al., 2017). The relative error is expressed as: where C r is the reflection coefficient calculated from the present numerical method (ISBM), and is the reference reflection coefficient calculated from a converged solution of the BEM.  Table 1.
Relative errors versus the total number of the boundary collocation nodes are shown for the single breakwater in Fig. 4 and double breakwater in Fig. 5. It is clearly seen for all tests that errors are smaller for the smaller values of kd and decrease with the increase of the number of the nodes.
Further, to have an idea of the computational cost of the present method CPU times versus the number of the boundary nodes using ISBM and BEM are shown in Fig. 6 respectively for single and double breakwaters. It should be mentioned that the computer program was written in Fortran and run on a PC with an Intel Core i5-3230M Processor (@2.60 GHz) and 4GB RAM memory. The results show that the ISBM takes far less CPU time than the BEM.
To further test the validity of the present method, the numerical results of the reflection coefficients of two cases are verified against the results of the BEM. In all subsequent computations the whole boundary of the computational domain is discretized with 400 source nodes for single breakwaters and 600 source nodes for dual breakwaters. This guarantees that the computational errors are small and do not exceed 10 -2 (see . Similarly for the multi-domain BEM simulations the computational domain is discret- ized with 400 constant elements for single breakwaters and 600 constant elements for dual breakwaters. In both methods the vertical boundaries are selected following the criteria of  such that and respectively for single and dual breakwaters. For rectangular breakwaters it is adequate to exchange by w. The first test case examined is a bottom standing single porous rectangular breakwater for h/d=0.5, and w/d=0.5. Only the front face of the breakwater is made permeable, e.g. ( )=(1, 0). The structure is located in front of a vertical wall such that the normalized distance . Variations of the coefficients of C r with respect to kd using the ISBM are presented in Fig. 7 against those of the BEM. It is clearly shown that the results of both methods are in close agreement.
The second case examined is a permeable structure of bottom standing dual rectangular breakwaters, such that h/d=0.5, and w/d=0.5. For the combined structure only the front face of the front breakwater is made permeable. All other faces are made impermeable, e.g. ( )= (1, 0, 0, 0). The breakwaters are separated by a distance such that =1. The normalized interspacing between the back breakwater and vertical wall is set to . Results of the reflection coefficients C r of the structure are shown in Fig. 8 for both ISBM and BEM. It is evident that results of both methods are similar.

Single breakwaters
In Fig. 9, the variations of C r are shown versus kd for different values of h/d. The conditions for the calculations are such that ( )=(0.5, 0), and . The results in Fig. 9a for rectangular breakwaters with w/d=1 clearly show that C r is oscillatory with several peaks occurring at certain values of kd. At these peaks the reflection coefficients approach unity. The primary larger peaks occur at smaller values of kd. The maximum and minimum values of the curves are due to resonance between the incoming waves and those   reflecting from the back impermeable vertical wall. With the increase of kd, the troughs of the oscillating curves increase. The rising of the height of the porous breakwater has clear effect on C r . Increasing the breakwater height in the presence of the back impermeable results in smaller reflection coefficients indicating that high breakwaters function better in terms of coastal protection.
For a trapezoidal breakwater with a base of the same width as the rectangular breakwater ( =1) and top width half of the width of the base ( =0.5), the results in Fig. 9b are similar to some extent to those of the rectangular breakwater, but the spectrum of C r becomes slightly larger. The primary peaks tend to have full reflection; however the smaller values of C r slightly decrease compared with their counterparts in the rectangular breakwater. In addition, bear in mind that trapezoidal breakwaters require less construction materials than rectangular breakwaters, and this suggests that trapezoidal breakwaters are more cost-effective and constitute a better alternative in terms of shoreline protection.
In Fig. 10, the variations of C r are shown versus kd for different values of the breakwater widths. The calculations are for h/d=0.75, ( )=(0.5, 0), and . The results of C r in Fig. 10a for rectangular breakwaters do not show any prominent effects when the width to the breakwater increases especially in the smaller range of kd. The curves do not preserve the same periodicity when kd further increases. For a trapezoidal breakwater with a base of the same width as the rectangular breakwater and top width half of the width of the base ( =0.5), the results in Fig. 10b are fundamentally similar to those of the rectangular case in Fig. 10a. Fig. 11, the variations of the reflection coefficients C r versus kd are shown for h/d=0.75, and different values of ( ). The results in Fig. 11a for rectangular breakwaters with w/d=1, clearly show that when the permeability parameter increases, the minimum values of C r decrease. However the maximum peaks are not affected by the permeability. This implies that at the minimum values of C r more wave energy is dissipated, whereas at the maximum values of C r hardly any wave energy is dissipated. As discussed by Huang et al. (2011) and Koley et al. (2015a)    reduction in the wave reflection occurs by a combination of the wave energy dissipation due to the structure's permeability, and the phase change of the incident and reflected waves in the confined region between the permeable structure and vertical seawall. The change in phase of the incident and reflected waves leads to constructive and/or destructive interference of waves outside/within the confined region between the permeable structure and rigid seawall. The maximum in the wave reflection is caused by the constructive interference of the incident and reflected waves (in phase) by the permeable structure, whereas the minimum in the wave reflection is caused by the destructive interference of the incident and reflected waves (out of phase). Practically similar results are shown in Fig. 11b for trapezoidal breakwater with =1 and =0.5.
Effects of the back wall proximity on the reflection coefficients are shown in Fig. 12 where the variations of C r are plotted against for h/d=0.75, kd=1.5 and different values of . The permeability parameter is set to zero. For rectangular breakwaters with w/d=1, the results of C r in Fig. 12a show periodic variations with the increase of . When the permeability decreases, the minimum values decrease accordingly, but again the maximum values are not affected. It is noticed that the minimum values of C r are achieved at certain values of . This is due to the maximum out of phase wave resonance between the incident and reflected waves (in the confined region between the permeable breakwater and the vertical back wall (Huang et al., Koley et al., 2015a)). When the permeability increases, the reflection decreases even further. Therefore, the relative spacing must be carefully designed to acquire smaller reflection coefficients and hence better sheltering. For a trapezoidal breakwater with =1 and =0.5, the results of C r in Fig. 12b show similar trends, except that the minimums are slightly smaller than their counterparts in the rectangular breakwater. When the right permeability is used trapezoidal breakwaters can in fact deliver smaller reflections.

Double breakwaters
The double breakwaters in this section are considered to be of the same height of the rectangular breakwaters presented previously, but have only half their widths. To represent the permeability in this study only the front wall of the front breakwater is made permeable. All other walls are impermeable ( )=(0, 0, 0).
In Fig. 13, the variations of C r are shown versus kd for different values of h/d, =0.5, the normalized distance between breakwaters =3, and normalized distance between the back breakwater and vertical wall . The results in Fig. 13a for rectangular breakwaters with w/d=0.5 clearly show that the variations are periodic with the occurrence of a multitude of peaks. The larger peaks occur at the smaller values of kd. The minimum values of C r increase with the increase of kd. It is also evident that increasing the height of the breakwater decreases the minim- Fig. 11. Variations of C r versus kd for different (G 1 , G 2 ). Nadji CHIOUKH et al. China Ocean Eng., 2019, Vol. 33, No. 2, P. 148-159 155 For trapezoidal breakwaters with a base of the same width as the rectangular breakwaters ( =0.5) and top width half of the width of the base ( =0.5), the results in Fig. 13b show similar trends. The overall width of both breakwaters (2w/d=1 for rectangular and =1 for trapezoidal) is similar to the single breakwater of Fig. 9. However, the performance of using two breakwaters separated by a distance is seen to be much better than using a single breakwater. The widths of the spectrums around the peaks become narrower leading to the minimum reflections persisting longer.
Variations of C r versus kd for different values of the breakwater widths are shown in Fig. 14 for h/d=0.75, =3, , and =0.5. The results in Fig. 14a for rectangular breakwaters clearly show that the variations of C r are not affected much by the increase of w/d. This suggests that increasing the width of the breakwater does not improve the performance of the system. Similar trends are shown in Fig. 14b for trapezoidal breakwaters with a base of the same width as the rectangular breakwaters and =0.5. The spectrum is not much different from that of the rectangular breakwaters. Thus, it is advised to adopt trapezoidal breakwaters bearing in mind that the massive savings it would make in terms of construction time and costs. By comparing the results of Fig. 14 for dual breakwaters and results of Fig. 10 for single breakwaters it is again seen that using two breakwaters separated by a distance per-form much better than using a single breakwater, due to the minimum reflections which persist longer.
Effects of the permeability on the variations of C r versus kd are depicted in Fig. 15 for h/d=0.75, =3, , and different values of . The results in Fig. 15a are for rectangular breakwaters with w/d=0.5, and in Fig. 15b are for trapezoidal breakwaters with =0.5 and =0.5. Similar to the single breakwater results in Fig. 11, the occurrence of the minima in the wave reflection is due to the out of phase destructive interference between the incident and the reflected waves. Effects of increasing are seen to decrease the minimum values of C r even further, but the maximum values are not affected at all. In the dual breakwater case the minimum reflections persist longer implying the superiority of the dual systems even though their total width is similar to the width of the single breakwater in Fig.  11. The results of the trapezoidal breakwaters are quite similar to those of the rectangular breakwaters, suggesting that it would be more practical to adopt trapezoidal breakwaters.
Effects of varying the relative distance between the back breakwater and vertical seawall on C r are shown in Fig. 16 for h/d=0. 75, kd=1.5, , and different values of . The results in Fig. 16a are for rectangular breakwaters with w/d=0.5, and in Fig. 16b for trapezoidal breakwaters with =0.5 and =0.5. In Fig. 17 results of C r are depicted for the effects of increasing the relative distance between the breakwaters for . It is shown that all curves of C r vary periodically with the increasing values of or . The minimum values of C r occur at certain values of and . Again this is most probably caused by the destructive wave resonance in the space between the breakwaters and between the back breakwater and vertical seawall. Increasing the permeability tends to decrease the minimum values of C r . In practice, it is important to design a system with proper values of , and breakwater permeability to deliver minimum values of reflections.

Conclusions
The work carried out in this numerical study has employed the improved version of the original singular bound-ary method (ISBM) to evaluate the hydrodynamic quantities of the reflection of bottom standing submerged porous breakwaters in regular normally incident waves and in front of a vertical wall. Both single and dual prismatic breakwaters of rectangular and trapezoidal forms were examined. Correctness of the results of the present method, including accuracy and performance, were confirmed by comparison with the results of a multi-domain boundary element method (BEM) developed independently for a previous study. Effects of major design parameters, including the breakwaters height, width, spacing, and permeability were investigated for several wave conditions. The results demonstrate that in places where only partial protection from waves is required submerged breakwaters can be used successfully,   Nadji CHIOUKH et al. China Ocean Eng., 2019, Vol. 33, No. 2, P. 148-159 157 as they can substantially attenuate waves. Dual breakwaters were found to perform better than single breakwaters. The minimum reflections decrease with the increase of the height of the breakwaters, but the width was found to have a very limited effect. Increase of the permeability generally decreases the minimum reflections. These coefficients vary periodically with the interspacing (between the breakwaters and between the back breakwater and back vertical wall) relative to the wavelength. The results of rectangular and trapezoidal breakwaters are very similar. It is recommended to adopt trapezoidal breakwaters for shoreline protection bearing in mind the added massive savings it would make in terms of construction time and costs.