K A Belibassakis

G Arnaud

V Rey

J Touboul

Propagation and scattering of waves by dense arrays of impenetrable cylinders in a waveguide

Keywords:

Coupled modal-BEM for propagation and scattering by dense arrays of impenetrable cylinders in a waveguide.

• Comparison of modal-BEM with simplified model based on Foldy-Lax theory and experimental data.

• Efficient calculation of reflection and transmission properties, supporting design and optimization studies.

• Present models can be used for similar scattering problems in acoustic and electromagnetic waveguides.

Introduction

Porous structures are often studied in the context of shoreline protection, see, e.g., [START_REF] Sulisz | Wave reflection and transmission at permeable breakwaters of arbitrary cross section[END_REF][START_REF] Yu | Wave motion through porous structures[END_REF]. A particular kind is artificial coastal cellular breakwaters which contain openings of very narrow width, typically much less than one wavelength; see, e.g., [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: Interference process and specific surface effects on damping[END_REF][START_REF] Koraim | Mathematical study for analyzing caisson breakwater supported by two rows of piles[END_REF][START_REF] Manam | Scattering of water waves by vertical porous barriers: An analytical approach[END_REF][START_REF] Somervell | Hydrodynamic characteristics of vertical cellular breakwater[END_REF] and the references cited there. The latter structures permit part of flow and energy to propagate to the lee side, and are exploitable for dissipation of incident wave energy and improvement of water quality in a basin. Furthermore, wave transformation and attenuation in mangrove forests and marine vegetation is also a subject of interest in coastal engineering, taking into account their role in supporting local economy and stabilizing and protecting the coastal zone. In many studies mangrove trunks and roots are treated as vertical cylindrical elements located in the water column, [START_REF] Massel | Surface wave propagation in mangrove forests[END_REF]; see also [START_REF] Chen | Numerical simulations of wave propagation over a vegetated platform[END_REF]. Particular arrangements used in experiments and numerical models are regular arrays of vertical cylinders; see e.g., [START_REF] Augustin | Laboratory and numerical studies of wave damping by emergent and near-emergent wetland vegetation[END_REF][START_REF] Lowe | Oscillatory flow through submerged canopies: 1. Velocity structure[END_REF][START_REF] Mei | Long waves through emergent coastal vegetation[END_REF][START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF]. In the case when the scatterer diameter is very small compared to the characteristic wavelength, typically of the order of a few percent, the structure could be macroscopically considered as porous medium. For waves of small steepness the problem is treated in the framework of linear theory. Also, if the length of the porous structure is small, viscous dissipation can be approximately omitted.

For protection structures composed by dense arrays of vertical cylinders, laboratory physical models have been recently used to study and characterize the complex hydrodynamic processes involved during wave propagation within such porous media; see [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: Interference process and specific surface effects on damping[END_REF]. Furthermore, sloshing experiments in a rectangular tank, completely filled with vertical cylinders undergoing forced horizontal motion around the natural frequency of the first sloshing mode, have been used by Molin and Remy [START_REF] Molin | An investigation of the resonant sloshing motion in a rectangular tank with multiple vertical cylinders[END_REF] to identify the dispersion equation of linear waves propagating through such dense arrays of vertical cylinders; see also [START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF]. However, when considering the same porous structures with gaps, as e.g. in periodic arrangements (as in the case of segmented porous breakwaters) or during tests in a rectangular wave tank where the structure occupies the middle region, as studied by Arnaud [START_REF] Arnaud | Houle a la Cote: Propagation, Impacts Et Ouvrages Innovants[END_REF], additional phenomena enter into play, like scattering and interaction with the tank walls, which forms a waveguide, and the appearance of critical frequencies.

When the vertical cylinders in the examined finite-width protection structure extent from the top (free surface) to the bottom surface and the depth is constant, the problem can be treated by the Helmholtz equation on the horizontal plane, supplemented with additional conditions at the sides of the tank or for accounting for the effects of periodic boundaries in the case of segmented porous breakwater. The latter constitutes a scattering problem by many impenetrable bodies in a waveguide. Similar problems are frequently encountered in acoustics, elasticity and electromagnetism; see, e.g., [START_REF] Cai | Acoustical scattering by arrays of cylinders in waveguides[END_REF][START_REF] Martin | Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles[END_REF][START_REF] Thierry | mu-diff: a Matlab toolbox for multiple scattering problems by disks[END_REF] and the references cited there. Also, similar applications have been studied by Li and Mei [START_REF] Li | Bragg scattering by a line array of small cylinders in a waveguide. Part 1. Linear aspects[END_REF] in the case of array of cylinders in a long channel, where multiple scattering can excite Bragg resonances if the frequency falls within one of the band gaps.

In this work a coupled numerical scheme, based on a modal expansions, in conjunction with boundary element technique, is developed for treating propagation and scattering by finite structures consisted of regular arrangements of impenetrable cylinders inside a waveguide. Numerical results are presented and discussed concerning the wave details both inside and outside the porous structure, illustrating its applicability operating as breakwater unit in a periodic arrangement. Subsequently, focusing on the reflection and transmission properties far from the structure, a simplified model is derived, based on Foldy-Lax theory; see, e.g., Martin [START_REF] Martin | Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles[END_REF]. The latter model provides an equivalent refraction index for the wave scattering by the porous structure located strictly inside the waveguide; see e.g., [START_REF] Cassier | Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: Mathematical justification of the Foldy-Lax model[END_REF][START_REF] Ahmad | The equivalent refraction index for the acoustic scattering by many small obstacles: with error estimates[END_REF]. Results obtained by the present fully coupled and simplified methods are compared against experimental measurements collected in the wave tank of SeaTech of University of Toulon, France, showing good agreement. The present models support further design and optimization studies concerning the examined porous structures, and can be used for studying similar scattering problems in acoustic and electromagnetic waveguides.

Mathematical formulation

We consider small-amplitude harmonic water waves of angular frequency ω, propagating in the x-direction in a channel or a tank of constant water depth h, with vertical side walls located at y = ±d; see Fig. 1. In this environment the waves are scattered by a structure of length L and breadth B, located at the center of the domain. The structure consists of multiple bottom mounted vertical cylinders of radius a = D/2 in a symmetric arrangement, however, the extension of the present method to treat more general arrangements is straightforward. The above configuration models also the interaction of normally incident waves with a y-periodic arrangement of such structures of breadth B and gaps 2d -B between them. Furthermore, we will consider that the radius of the cylinder(s) is small with respect to the incident wavelength ka ≪ 1,

where k denotes the characteristic wavenumber ω 2

= kg tanh (kh) .

(

) 1
We consider the multiple scattering problem in the above waveguide excited by any superposition of the incoming propagating modes by

Φ + (x, y, z; t) = Re ( Φ + (x, y, z) exp (-iωt) ) , (2)
where i = √ -1 and the complex potential Φ + (x, y, z) satisfies the Laplace equation and the no-entrance boundary conditions at the bottom z = -h and the side walls y = ±d. Since the boundaries are vertical surfaces, evanescent modes corresponding to the imaginary roots of the dispersion relation are not excited and the three-dimensional incident wave potential is given by

Φ + (x, y, z) = cosh (k (z + h)) cosh (kh) Np ∑ n=1 A + n Y n (y) exp (ik xn x) . (3)
Considering the above incident field to be symmetric with respect to the centerline (y = 0) and to satisfy the no-entrance boundary conditions on the side walls (y = ±d), the transverse eigenfunctions and eigenvalues appearing in the above expansion are given by

Y n (y) = W n cos ( (n -1) π y d ) , and k xn = √ k 2 - π 2 (n -1) 2 d 2 , n = 1, 2, 3, . . ., (4)
where the constants W n can be fixed by normalization, and thus,

W 1 = 1/

√

2d and

W n = 1/ √ d, n = 2, 3, . . . . (5)
We easily derive from the above relations the number of propagating wavenumbers

N p = 1 + [kd/π ] , (6)
where [•] denotes the integer part. It is also seen from the above relations, that for water wavenumbers k (ω) = (n -1) π /d, n = 2, 3, 4.., which correspond to the following frequencies

ω n = √ ((n -1) gπ /d) tanh ((n -1) πh/d), n = 2, 3 . . . , (7)
a new mode, with index n, is added to the set of propagating terms. It is seen from Eqs. ( 3) and ( 4), that at these frequencies the new mode behaves like a stationary field that could possibly develop a resonance in the waveguide. Furthermore, using

Eq. ( 7), for n = 1, we obtain the following frequency

ω 1 = √ (gπ/d) tanh (πh/d), (8)
which identifies the upper limit of frequency for possible excitation of the waveguide only by the first (n = 1) mode corresponding to plane wave. Thus, for ω < ω 1 the incident field in the waveguide consists only by the parallel wave

Φ + (x, y, z) = A + 1 cosh (k (z + h)) cosh (kh) exp (ikx) , (9)
while for higher frequencies, ω > ω 1 , the incident field Φ + may contain in general more propagating modes (Eq. ( 3)).

Using the method of separation of variables complete representations of the wave field are derived in the form of normal mode series; see e.g., [START_REF] Jensen | Computational Ocean Acoustics[END_REF]. In the upwave region (left semi-infinite strip) the wave field consists of incident, reflected and evanescent components (generated from the scatterer), with amplitudes

A (1) n , Φ (1) (x, y, z) = cosh (k (z + h)) cosh (kh) ⎛ ⎝ Np ∑ n=1 A + n Y n (z) exp (ik xn x) + ∞ ∑ n=1
A (1) n Y n (y) exp (-

ik xn x) ⎞ ⎠ , (10a)
while in the downwave region (right half-strip) the field is

Φ (3) (x, y, z) = cosh (k (z + h)) cosh (kh) ∞ ∑ n=1 A (3) n Y n (y) exp (ik xn x) , (10b)
and consists only of transmitted and evanescent modes, with amplitudes A

(3)

n . Since the incident wave field contains only the propagating water-wave modes and the solid boundaries (both side walls of the channel and the multiple scatterers) are vertical surfaces, the evanescent water-wave modes characterized by vertical structure cos (k n (z + h)) /cos (k n h) and wavenumbers k n , n = 1, 2, .. corresponding to the imaginary roots of the dispersion relation (1) on the complex plane, are not excited. The total three-dimensional field in the considered waveguide contains both the incident and the diffracted component due to structure composed by many vertical cylinders, Φ (x, y, z) = Φ + (x, y, z) + Ψ (x, y, z), and is represented by

Φ (x, y, z) = ϕ (x, y) cosh (k (z + h)) cosh (kh) , (11a)
where ϕ (x, y) = ϕ + (x, y) + ψ (x, y) and

ϕ + (x, y) = Np ∑ n=1 A + n Y n (y) exp (ik xn x) . (11b)
Thus, the problem is reduced to the Helmholtz equation on the horizontal plane with respect to the surface wave potential ϕ (x, y) (or the free surface elevation η (x, y)) satisfying

( ∇ 2 + k 2 ) ϕ (x, y) = 0, in D = {-∞ < x < ∞, -d ≤ y ≤ d } , (12a)
∂ n ϕ (x, y) = 0, (x, y) ∈ ∂D C = ∪ i=1,N C ∂D i , (12b)
∂ y ϕ (x, y) = 0, y = ±d, (12c) ϕ
(x, y) ∼ Np ∑ n=1 A + n Y n (y) exp (ik xn x) + reflected waves for x → -∞, -d ≤ y ≤ d, (12d) ϕ (x, y) ∼ transmitted waves for x → -∞, -d ≤ y ≤ d, ( 12e
)
where ∇ = ( ∂ x , ∂ y )

is the surface gradient, n indicates the normal vector on the boundary of the N C cylinders (∂D i ) directed to the exterior of the fluid domain and ∂ n ϕ (x, y) the corresponding normal derivative of the wave potential on the horizontal plane.

Some properties of the scattering field ψ can be revealed by using the integral representation

ψ (x, y) ≈ ∫ β (x s ) G (x|x s ) dx s ,
involving the Green's function G (x|x s ) in the considered waveguide (also termed as a line source), which is given by (see [START_REF] Jensen | Computational Ocean Acoustics[END_REF]):

G (x|x 0 ) = G (x, y|x s , y s ) = i 2 ∑ n=1 Y n (y) Y n (y s ) exp (ik xn |x -x s |) k xn
, where x = (x, y) denotes the field point and x s = (x s , y s ) the corresponding source with density β, respectively. From the above it is seen that at the critical frequencies specified by Eq. ( 9), the x-wavenumbers k xn = 0 (see Eq. ( 4)), and if these modes are excited, the scattering field could possibly become singular resonating the waveguide. Before proceeding to the reformulation of the problem, Eqs. [START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF], as a matching boundary value problem based on the fundamental solution of the Helmholtz equation (or the free-space Green's function), it is worth mentioning here that the above complex Green's function is slowly convergent, especially near the singularity as the field point (x, y) approaches the source point (x s , y s ).

This fact leads to increased computational complexity and cost, especially concerning the calculation of the BEM influence coefficients. On the other hand, the evaluation of the latter quantities using the free-space Green's function is much easier and faster, at the cost of increased discretization required due to the treatment of additional boundaries (the side walls and the matching interfaces) that need to be considered.

Using Eqs. [START_REF] Lowe | Oscillatory flow through submerged canopies: 1. Velocity structure[END_REF], the above multiple scattering problem in the infinite strip -∞ < x < ∞, -d ≤ y ≤ d, is reformulated as a matching-boundary value problem in the bounded subdomain -ℓ < x < ℓ, -d ≤ y ≤ d containing strictly into its interior the N c cylindrical scatterers, as follows:

( ∇ 2 + k 2 ) ϕ (2) (x, y) = 0, in D (2) {-ℓ < x < ℓ , -d ≤ y ≤ d} , (13a)
∂ n ϕ (2) (x, y) = 0, (x, y) ∈ ∂D C , (13b)

∂ y ϕ (2) (x, y) = 0, -ℓ < x < ℓ , y = ±d, (13c)
ϕ (2) (x, y) = ϕ (1) (x, y) and

∂ x ϕ (2) (x, y) = ∂ x ϕ (1) (x, y) , x = -ℓ, -d ≤ y ≤ d, (13d) ϕ (2) (x, y) = ϕ (3) (x, y)
and

∂ x ϕ (2) (x, y) = ∂ x ϕ (1) (x, y) , x = ℓ, -d ≤ y ≤ d, (13e)
where x = ±ℓ denotes the position of the interfaces shown also by using dotted lines in Fig. 1(a), and

ϕ (1) (x, y) = Np ∑ n=1 A + n Y n (y) exp (ik xn x) + ∞ ∑ n=1 A (1) n Y n (y) exp (-ik xn x) , (14a)
ϕ (3) (x, y) = ∞ ∑ n=1 A (3) n Y n (y) exp (ik xn x) . (14b)
We note that all terms in the above infinite series with index greater than N p present exponential decay (always faster for larger n) because the wavenumbers k xn defined from Eq. ( 4) are imaginary always increasing in modulus. In the sequel the series appearing in the above equations are truncated keeping only a finite number N m of terms (N m > N p ) that are sufficient for numerical convergence.

Given the amplitudes of the incident modes

{ A + n , n = 1, 2, . . . , N p } , the problem is solved in the box-shaped bounded subregion (-ℓ < x < ℓ, -d ≤ y ≤ d)
, finding the unknown wave potential ϕ (2) and the sets of unknown coefficients

A (1)
n , A (3) n , n = 1, 2, . . . . In the following, a hybrid modal-BEM will be presented for treating the above problem. For simplicity in the presentation, from now on the wave in the middle region ϕ (2) (x, y) will be simply denoted as ϕ (x, y) and the corresponding domain as D.

Moreover, energy flux conservation is easily derived by using Green's theorem for the field in the waveguide, in conjunction with the wall conditions Eqs. (12c), from which we finally obtain

Im {∫ y=d y=-d ϕ * (x, y) dϕ (x, y) dx dy } = const, (15a)
for any transverse section x = const of the waveguide outside the scatterer, where ϕ * denotes the complex conjugate of the potential. Applying the above equation, in conjunction with the representations [START_REF] Arnaud | Houle a la Cote: Propagation, Impacts Et Ouvrages Innovants[END_REF] in the far upwave region x ≪ -ℓ and in the far downwave region x ≫ -ℓ, and using the fact that the functions ∥Y n ∥ = 1, the following relation is derived:

Np ∑ n=1 k xn ( ⏐ ⏐ A + n ⏐ ⏐ 2 - ⏐ ⏐ A (1) n ⏐ ⏐ 2 ) = Np ∑ n=1 k xn ⏐ ⏐ A (3) n ⏐ ⏐ 2 . ( 15b
)
The above relation will be used to check satisfaction of energy conservation from the calculated values of the reflection and transmission amplitudes, as obtained by the present modal-BEM that will be described in detail in the next sections.

The matching-boundary value problem

In the bounded (middle) subregion the wave field is represented by a single layer boundary integral representation based on simple sources, as follows

ϕ (x) = ∫ ∂D σ (x s ) G (x|x s ) ds (x s ) , x ∈ D ∪ ∂D, (16a)
where x = (x, y) denotes the field point, x s = (x, y) ∈ ∂D the source point on the boundary, and the involved function G (x|x s ) is the fundamental solution of Helmholtz equation on the plane, which should not to be confused with the Green's function G (x|x s ) of the previous section. The fundamental solution is defined by means of the Hankel function of zero order and first kind as follows Using the jump relations of the single layer we obtain from Eqs. ( 16) the following integral representations concerning the derivatives of the potential

G (x|x s ) = i 4 H (1) 0 (k |x -x s |) . (16b)
∇ϕ = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∇ x ∫ ∂D σ (x s ) G (x|x s ) ds (x s ) = ∫ ∂D σ (x s ) ∇ x G (x|x s ) ds (x s ) , for x ∈ D σ (x) n 2 + ∫ ∂D σ (x s ) ∇ x G (x|x s ) ds (x s ) , for x ∈ ∂D (17)
Then by calculating the wave velocity for points on the boundaryx ∈ ∂D we obtain the representation of the normal derivative of the potential ∂ϕ ∂n

⏐ ⏐ ⏐ ⏐ ∂D ≡ n∇ϕ = σ (x) 2 + ∫ ∂D σ (x s ) ∂ n G (x|x s ) ds (x s ) , x ∈ ∂D, (18)
where n denotes the unit normal vector on the boundary (directed outwards the domain D), as shown in Fig. 2, and

∂ n G (x|x s ) = n ∇ x G (x|x s ).
Then, the Neumann boundary condition on the solid boundary of the scatterers and the side walls is

σ (x) 2 + ∫ ∂D σ (x s ) ∂ n G (x|x s ) ds (x s ) = 0 , for x ∈ ∂D C and -ℓ < x < ℓ , y = ±d. ( 19
)
The formulation is completed by including the matching condition at the two interfaces separating the middle domain D from the two semi-infinite strips. To this aim using the truncated versions of Eqs. ( 14) keeping only a finite number N m of terms N m > N p , and the boundary integral representations ( 16) and [START_REF] Thierry | mu-diff: a Matlab toolbox for multiple scattering problems by disks[END_REF], we obtain for points x on the left interface {x = -ℓ, -d ≤ y ≤ d } the following equation expressing the matching of the potential

ϕ| x=-ℓ ≡ ∫ ∂D σ (x s ) G (x|x s ) ds (x s ) = Np ∑ n=1 A + n Y n (y) exp (ik xn ℓ) + Nm ∑ n=1
A (1) n Y n (y) exp (-ik xn ℓ) , [START_REF] Ahmad | The equivalent refraction index for the acoustic scattering by many small obstacles: with error estimates[END_REF] and as concerns the matching of the normal derivative ∂ϕ ∂n (1) n Y n (y) exp (-ik xn ℓ) .

⏐ ⏐ ⏐ ⏐ x=-ℓ ≡ σ (x) 2 + ∫ ∂D σ (x s ) ∂ n G (x|x s ) ds (x s ) = - Np ∑ n=1 ik xn A + n Y n (y) exp (ik xn ℓ) + Nm ∑ n=1 ik xn A
A (1) n = exp (ik xn ℓ) (⟨ ϕ| x=-ℓ , Y n ⟩ -A + n exp (ik xn ℓ) ) , n = 1, 2, 3 . . . , (22)
where ⟨f (y) , g (y)⟩ = ∫ y=ℓ y=-ℓ f (y) g (y) dy denotes the inner product in -d ≤ y ≤ d.

Using the above result in Eq. ( 21) we obtain the final form of the non-local boundary condition expressing the matching of the potential on the left interface x = -ℓ:

∂ϕ ∂n ⏐ ⏐ ⏐ ⏐ x=-ℓ - Nm ∑ n=1 ik xn ⟨ ϕ| x=-ℓ , Y n ⟩ Y n (y) = -2ik xn A + n exp (ik xn ℓ) Y n (y) , -d ≤ y ≤ d, (23)
where A + n is the n-amplitude of the incident mode constituting the forcing of the problem. Quite similarly, for points on the right interface {x = ℓ, -d ≤ y ≤ d } we obtain the following result concerning the coefficients of the expansion of the potential in the right half strip

A (3) n = exp (-ik xn ℓ) ⟨ ϕ| x=ℓ , Y n ⟩ , n = 1, 2, 3 . . . . (24)
and using this result we obtain the following non-local condition, expressing the matching of the potential on the right interfacex = ℓ:

∂ϕ ∂n ⏐ ⏐ ⏐ ⏐ x=ℓ - Nm ∑ n=1 ik xn ⟨ ϕ| x=ℓ , Y n ⟩ Y n (y) = 0, -d ≤ y ≤ d ( 25
)
where

ϕ| x=-ℓ ≡ ∫ ∂D σ (x s ) G (x|x s ) ds (x s )
and ∂ϕ ∂n

⏐ ⏐ ⏐ ⏐ x=-ℓ ≡ σ (x) 2 + ∫ ∂D σ (x s ) ∂ n G (x|x s ) ds (x s ) considered for points {x = ℓ, -d ≤ y ≤ d}.
The examined problem possesses a unique solution, with the possible exception of a set of discrete frequencies constituting the spectrum of irregular frequencies connected with the single layer representation used; see, e.g., [START_REF] Marburg | Treating the phenomenon of irregular frequencies[END_REF][START_REF] Martin | Identification of irregular frequencies in simple direct integral equation methods for scattering by homogeneous inclusions[END_REF] and the references cited there. Also, the possible existence of trapped modes at specific frequencies could cause non-uniqueness of the problem (see, e.g., [START_REF] Mciver | Embedded trapped modes for obstacles in two-dimensional waveguides[END_REF]). In the present work we focus on the numerical solution of the above system in the case of regular frequencies only, leaving the study of the interesting problem concerning irregular frequencies and singular conditions in future extensions.

The modal Boundary Element Method

The problem is solved by applying a Boundary Element Method to calculate the integral operators, in conjunction with a collocation technique. The external boundary of the middle domain consists of four parts (walls, incidence and transmission sides) and the boundaries of discs in the arrangement, numbered consecutively as indicated in Fig. 2. Moreover, linear elements are used to discretize the various parts of the boundary, ensuring continuity of geometry, in conjunction with a piecewise constant approximation of the source distribution σ (x) , x ∈ ∂D. Thus, the unknowns are the value of the density in each element σ j , j = 1, 2, . . . , N el , where N el is the total number of elements used to discretize all parts of the boundary, and the equations are satisfied at specific control points x j ∈ ∂D, selected to be the centers of the boundary elements, j = 1, 2, . . . , N el ; see Fig. 2. In this way, the integrals representing the values of the potential and its normal derivative on the boundary involved in the boundary and matching conditions Eqs. ( 19), ( 23) and [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF], are approximated by ϕ (

x j ) = N el ∑ n=1 σ n F jn , ∂ n ϕ ( x j ) = N el ∑ n=1 σ n U jn . (26)
In the above relations, the matrices F jn and U jn are defined by the following integrals defining induced quantities from the n-linear element to the j-control point, respectively,

F jn = ∫ (n) H (1) 0 ( k ⏐ ⏐ x j -x (s) ⏐ ⏐ ) ds, ( 27
)
and

U jn = δ jn 2 -k ∫ (n) H (1) 1 ( k ⏐ ⏐ x j -x (s) ⏐ ⏐ ) n x n ( x j -x (s) ) + n y n ( y j -y (s) ) ⏐ ⏐ x j -x (s) ⏐ ⏐ 2 ds. ( 28
)
In the above equations x (s) = (x (s) , y (s)) denotes the parametric equation of the n-element with s denoting the physical length, n = (

n x j , n y j )
is the outer unit normal vector at the center of the j-element, and H

(1)

1 is the Hankel function of first kind and first order. For j ̸ = n the integrals in Eqs. ( 27) and ( 28) are regular and, for small nondimensional wavenumbers kδs, where δs is a representative value of the meshsize on the boundary (element length), and small distances k ⏐ ⏐ x j -x n ⏐ ⏐ , they can be calculated by usual numerical integration rules, and in the present work an adaptive Simpson's rule is applied. Difficulties could arise for higher frequencies, which is the case in similar problems in acoustic and electromagnetic waveguides, where these integrals are more efficiently calculated by means of Filon's rule; see, e.g. [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF]. For large distances k

⏐ ⏐ x j -x n ⏐ ⏐ ≫ 1, the
above integrals are more rapidly and efficiently calculated by means of multipole expansions of the kernels; see, e.g. [START_REF] Liu | Fast Multipole Boundary Element Method -Theory and Applications in Engineering[END_REF]. On the other hand, for j = n the integral in Eq. ( 27), expressing the induced potential from a linear element to its center, is weakly singular. These integrals contain logarithmic singularity and are conveniently treated by transformation of the integration variable. In particular, simple transformations of the form u 2

= (ss mid ), where s mid the coordinate of the middle of the element, bringing the integration points closer to the singularity, are sufficient in order to reduce the integrals to regular ones and apply similar as before numerical rules. Finally, we notice that for j = n the singular integral in Eq. ( 28) exists in the sense of Cauchy principal value. However, it becomes simply zero due to symmetry properties of the kernel about the point of singularity. The above facts result in a very fast calculation of the matrices F jn and U jn even for complex configurations discretized by using large number of elements.

On the basis of the above, the boundary condition Eq. ( 19) on the surface of the scatterers ∂D C and on the side walls y = ±d, becomes

N el ∑ n=1 σ n U jn = 0, for j ∈ {J (∂D c ) , J (y = ±d) } , (29)
where J (∂D c ) , J (y = ±d) denote the corresponding set of indices of the panels on these parts of the boundary ∂D. The infinite series appearing in Eqs. ( 23) and [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF], respectively, are truncated keeping only a finite number of terms N m > N p , exploiting the fact that all terms above N p present exponential decay (always faster for larger n) because the wavenumbers k xn defined from Eq. ( 4) are imaginary always increasing in modulus. Furthermore, using the BEM representations [START_REF] Liu | Fast Multipole Boundary Element Method -Theory and Applications in Engineering[END_REF] in Eqs. ( 23) and ( 25) we obtain:

N el ∑ n=1 σ n U jn -i Nm ∑ n=1 k xn Y n ( y j ) N el ∑ p=1 σ p F mp ⟨H m , Y n ⟩ = -2ik xn ( δ nm A + m ) exp (ik xn ℓ) Y n ( y j
) , p, m, and j ∈ J (x = -ℓ) ,

N el ∑ n=1 σ n U jn -i Nm ∑ n=1 k xn Y n ( y j ) N el ∑ p=1 σ p F mp ⟨H m , Y n ⟩ = 0, p, m, and j ∈ J (x = ℓ) , (30)
where H m (y) denotes the unit step function with support on the m-element of the vertical interfaces. Clearly, the above equations materialize the discrete form of the matching-boundary conditions by the present modal-BEM formulation, and the linear system is solved with respect to the discrete complex unknowns {σ n , n = 1, 2, . . . , N el }.

Validation in the case of scattering of parallel incident waves by hard cylinder

In order to illustrate the applicability and efficiency of the present method, we first consider the simple case of scattering of plane harmonic waves by a single hard cylindrical scatterer in infinite domain, corresponding to the limiting case of waveguide of infinite width (d → ∞). This will help us also to estimate the number of panels on each cylinder that is needed for convergence of the results. In this case analytical solution is available, which is obtained by separation of variables in cylindrical coordinates.

Assuming that the circular scatterer of radius a is located at the center of the domain (x, y) = (0, 0), the analytical solution is given by (see, e.g., [START_REF] Martin | Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles[END_REF][START_REF] Mei | The Applied Dynamics of Ocean Surface Waves[END_REF]), ϕ (r, θ) = exp (ikx) + ψ (r, θ) ,

(32a) ψ (r, θ) = ∞ ∑ m=0 A m cos (mθ) H (1) m (kr) , with A m = -ε m i m J ′ m (ka) Φ (1) m (ka) , (32b)
where r = ( 3 for N el = 31 and 61 elements on the surface of the impermeable cylinder and compared against the analytical solution. In this example, we consider plane water waves of frequency f = 1.5 Hz propagating in constant depth h = 0.23 m, and thus, the shallowness parameter is h/λ = 0.34 corresponding to intermediate water depth. The waves are scattered by a circular, bottom mounted cylinder of radius a = 0.5 m, which is comparable to the dimensions of the porous structure examined in this work as presented in more detail in the next section. Thus, the nondimensional wavenumber in the case of Fig. 3 is ka = 4.66. We clearly observe the excellent convergence of the present BEM scheme to the analytical solution.

x 2 + y 2 ) 1/2 , θ =

Numerical results for the many cylinder structure in the waveguide

In recent research by Arnaud et al. [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: Interference process and specific surface effects on damping[END_REF] wave propagation through arrangements of vertical cylinder arrays is investigated, that is examined for possible applications as permeable breakwaters and protective structures. In particular, laboratory physical models have been systematically tested in [START_REF] Arnaud | Houle a la Cote: Propagation, Impacts Et Ouvrages Innovants[END_REF] to characterize the complex hydrodynamic processes involved during wave propagation and to study the effect of specific surface or porosity on the wave energy dissipation; see Figs. 4 and5. Several arrangements with different cylinder diameters have been considered in order to study the effects of the specific surface while keeping the porosity constant. In particular, three cylinder diameters are considered in [START_REF] Arnaud | Houle a la Cote: Propagation, Impacts Et Ouvrages Innovants[END_REF], with D = 0.020 m, 0.032 m and 0.050 m, arranged in blocks of length L = 0.3 m and breadth B = 1.2 m. The cylinders are regularly disposed along two perpendicular axes forming a 45 deg angle with the longitudinal axis. Porosity in the case of the examined arrangements of vertical cylinders attain the same value, given in terms of diameter D and number N C of cylinders as follows:

γ = 1 -N C π a 2

LB

, where a = D/2.

(33)

The structure of Fig. 4 with cylinders of diameter D = 0.032 m and the indicated spacing is modeled by arranging 17 rows of 4 cylinders in x direction and, between them, 16 rows of 5 cylinders, i.e. a total number Nc = 148 cylinders, as shown in Fig. 6.

Using the above values in Eq. ( 33) it results in γ = 0.67 which approximates well the reported value in the experiments [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: Interference process and specific surface effects on damping[END_REF][START_REF] Arnaud | Houle a la Cote: Propagation, Impacts Et Ouvrages Innovants[END_REF]. Although the above value of porosity is relatively low, the term 'dense arrays of vertical cylinders' is also used herein to characterize this type of structures; see also [START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF]. The above structure operating as porous segmented breakwater is tested in the wave tank of Seatech shown in Fig. 5, in water depth h = 0.23 m. The tank is 10 m long and 2.6 m wide, and the porous breakwater is located in the middle of the tank so that the side gaps between the structure and the walls of the tank are equal to 0.7 m. For this configuration and water depth the critical frequencies estimated from Eq. ( 7) are listed in Table 1 below.

In the following, we consider forcing only by the first incident mode (plane incident waves) in the tank, and results

obtained by the present modal-BEM are shown in Fig. 6 for three frequencies f = 1 Hz, 1.5 Hz and 1.8 Hz, for which experiments were performed in the tank. Present calculations are based on total M = 152 boundary sections modeling the waveguide and all the cylinders in the porous structure. More specifically, 50 elements are used along the side walls and 64 elements in the entrance and exit boundaries. Also, 30 elements are used to discretize each one of the 148 cylinders in the arrangement, resulting in a total number of 4668 elements, which is found enough for numerical convergence. In the colorplots of Fig. 6 (left column) the real part and the modulus of the wave field, at each one of the above frequencies, is shown. In all cases we observe that the field is resolved with high accuracy, both outside and inside the arrangement of vertical cylinders composing the porous structure. In the right subplots of Fig. 6 the computed results along the section y = 0.1 m (at 10 cm from the centerline of the tank) are shown. The real part of the wave field is plotted by using thin solid line and the imaginary part by dashed line, respectively. Also, the modulus of the wave field, which is proportional to the distribution of the amplitude of the free-surface elevation along the same longitudinal section, is shown by using thick lines. In these plots, the limits of the structure are indicated by vertical dashed lines. We observe that at all frequencies the structure allows much part of the flow and energy to be transmitted to the downwave side. In particular, as concerns energy transmission, the examined structure seems to be more efficient at lower frequencies. Furthermore, although in cases (a)

and (b) the frequency (f = 1 Hz and 1.5 Hz) is close to the second and fourth critical value (see Table 1), respectively, still no special effect is reflected in the results, indicating that the possible resonances are very local.

For the examined configuration, more results concerning the reflection and transmission of waves are presented in Fig. 7, characterizing also the energy exchange between the first incident mode (with amplitude A + 0 ) and the higher modes that can be excited by the scattering structure in the tank, as the frequency exceeds the critical values of Table 1. In this figure the nondimensional coefficients associated with the reflection and the transmission modes are plotted, as calculated by the present modal-BEM. These coefficients are defined by the corresponding amplitudes in the representations, Eqs. (14a) and (14b), respectively, as follows

R n = ⏐ ⏐ A (1) n ⏐ ⏐ /A + 0 , and T n = ⏐ ⏐ A (3) n ⏐ ⏐ /A + 0 , n = 1, 2, 3, . . . . (34)
In particular, the first three coefficients, containing most of the energy, are plotted against the nondimensional wavenumber kh = 2π h/λ, for frequencies from 0.8 to 2.1 Hz, corresponding to wave conditions ranging from intermediate water depth to deep water. We observe in Fig. 7, similarly as in the case of finite underwater steps or submerged breakwaters, that the behavior of transmission and reflection coefficients presents alternating maxima and minima at specific frequencies. Moreover, as the frequency of the incident wave increases, higher modes are excited, carrying energy in the transmission region. We also remark here that the shown results, obtained by using the above discretization, are found to be consistent with the energy conservation relation Eq. (15b), presenting error less than 1%, which diminishes with further increase of the number of boundary elements.

In the following, we will consider the above as reference data to compare with results obtained by a simplified model and experimental measurements in the tank. The computational cost of the simplified model is significantly less that the above more complete modal-BEM solution, facilitating its systematic application to investigate the reflection and transmission properties of such porous structures in more complex environment and supporting further design and optimization studies. Comparative data for computational requirements and time will be provided at the end of Section 5.2.

A simplified coupled-mode model on the horizontal plane

In the case of the scattering of time-harmonic waves by large number of densely arranged small obstacles in a homogeneous medium, where the size of the bodies is quite small compared with the wavelength, the computation cost of the numerical simulation of the field by boundary element methods, in conjunction with a Dirichlet-to-Neumann conditions, can become excessively high, making difficult its use in optimization studies and application to scattering problems by many such structures in realistic environments. In this section an approximate, simplified model, with very low computational requirements, is presented and tested for the solution of the problem. The simplified model is based on the Helmholtz equation on the horizontal plane, and treats the porous structure as an inclusion in the waveguide, characterized by an effective index of refraction k eff . The latter is estimated as discussed in more detail in Appendix by an adapted version Foldy-Lax theory (see, e.g., [START_REF] Martin | Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles[END_REF]). Thus, the governing equation is

( ∇ 2 + κ 2 ) ϕ (x, y) = 0, (35)
where the wavenumber parameter in the water subregion is the standard one κ 2 = k 2 (as obtained from Eq. ( 1) for the given frequency and water depth), and in the porous subregion this parameter is approximated by an effective value as follows,

κ 2 = k 2 eff = γ k 2 , ( 36
)
where γ is the porosity coefficient (defined by Eq. ( 33)).

The solution of the Helmholtz equation ( 35) is obtained by a domain decomposition technique. The problem is solved taking into account the incidence-reflection condition at the entrance and the transmission condition at the exit of the waveguide, the side-wall boundary conditions, and the matching conditions on the interfaces separating the porous from the water subregions, shown by using dashed lines in Fig. 8. The usual interface conditions (see, e.g., [START_REF] Jensen | Computational Ocean Acoustics[END_REF]) concern the normal derivative ∂ϕ/∂n and the value of the wave field ϕ (x, y) on the two sides of the interface. In the present work we use the following conditions ∂ϕ ∂n which expresses the interface condition for the flow rate through the interface (see, e.g., [START_REF] Sulisz | Wave reflection and transmission at permeable breakwaters of arbitrary cross section[END_REF], Eq. 2.14), in conjunction with ϕ| water = S ϕ| porous , (38) which expresses the interface condition for the pressure across the same interface, where S is an empirically estimated coefficient (see, e.g., [START_REF] Sulisz | Wave reflection and transmission at permeable breakwaters of arbitrary cross section[END_REF], Eq. 2.13). It is remarked here that in several works concerning the propagation of water waves through or over porous media, see e.g., [START_REF] Sulisz | Wave reflection and transmission at permeable breakwaters of arbitrary cross section[END_REF][START_REF] Yu | Wave motion through porous structures[END_REF][START_REF] Gu | Gravity waves over porous bottom[END_REF], the model by Sollitt and Cross [START_REF] Sollitt | Wave transmission through permeable breakwaters[END_REF] is used to describe the porous flow. In the latter model an extra term is included in the momentum equation representing the effect of added resistance caused by the virtual mass of the grains in the porous medium, and when applied to the present configuration it results to a modification of the water-wave dispersion relation, as follows

⏐ ⏐ ⏐ ⏐ water = γ ∂ϕ ∂n ⏐ ⏐ ⏐ ⏐ porous , (37)
ω 2 S = k p g tanh ( k p h ) , (39)
from which the wavenumber k p associated with the propagating mode in the porous subregion can be calculated. In the latter works S = 1 + C M (1 -γ ) /γ is described as a coefficient representing the effect of inertia in the porous flow. The inertia coefficient C M for circular discs, especially in the frequency band and corresponding Reynolds numbers (wave inertial to viscous forces) considered here, takes values near unity. Thus, in the present case S ≈ 1/γ >1 and the wavenumber in the porous subregion derived from Eq. (39) would be k p > k. Furthermore, theoretical and experimental analysis involving sloshing experiments by Molin & Remy (2016) [START_REF] Molin | An investigation of the resonant sloshing motion in a rectangular tank with multiple vertical cylinders[END_REF], shows that the following modified form of the dispersion relation better describes the wave properties in the porous region

ω 2 √ S = k p g tanh ( k p h/ √ S ) , (40)
which results in k p = k/ √ γ , and thus also k p > k. Therefore, for S > 1, both the above dispersion relations Eqs. (39) and (40), produce values for the wavenumber k p that are greater than the corresponding value k in the water region (for the same frequency), and the appropriate use of the dispersion relation to represent the waves over the porous region is thoroughly discussed in [START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF]. On the other hand, in the present simplified model the wave parameter in the porous subregion is smaller than the water wavenumber, k eff = √ γ k < k. However, it should be stressed here that k eff is not representative of the water wave behavior over the porous structure, but rather it provides us with an effective index of refraction for modeling the scattering effect of the porous structure treated as an inclusion on the basis of applications based on the Helmholtz equation on the horizontal plane. In the sequel a coupled-mode system will be presented for solving the above problem Eqs. ( 35)-( 38) and approximating the performance of the porous structure with very small computation cost.

Coupled-mode system

On the basis of coupled mode theory, the problem is approximately solved on the horizontal plane by decomposing the domain into three parts denoted by (I), (II) and (III) in Fig. 8, where (I) and (III) are incident and transmission water subregions in the waveguide and the field is provided by the expansions (14a) and (14b), respectively. In the middle subdomain (II), which includes the porous medium modeling the structure, a similar modal expansion is introduced, as follows ϕ (2)

(x, y) = ∑ n=1 ϕ n (x) Y (2) n (y) , -ℓ ≤ x ≤ ℓ , -d ≤ y ≤ d p and d p ≤ y ≤ d , (41)
where the new functions

{ Y (2) n (y) , n = 1, 2, . . . . } are symmetric, Y (2)
n (y) = Y (2) n (-y), in -d ≤ y ≤ 0 and 0 ≤ y ≤ d, respectively, and are obtained as the solution of the following eigenvalue problem in 0 ≤ y ≤ d :

d 2 Y (2) n (y) dy 2 + ( k 2 (y) -k 2 n ) Y (2) n (y) = 0, in 0 < y < d, (42a)
Y (2) n

( y = d p + 0 ) = SY (2) n ( y = d p -0 ) , (42c) dY (2)
n

( y = d p + 0 ) dy = γ dY (2) n ( y = d p -0 ) dy , (42d)
where

k 2 (y) = k 2 for d > y > d p and k 2 (y) = k 2 eff for 0 < y < d p . (42e)
The eigenfunctions associated with the above problem are given by

Y (2) n (y) = W cos (√ k 2 eff -k 2 n y ) , in 0 < y < d p , (43a)
Y (2) n (y) = W cos ( √ k 2 -k 2 n d p ) cos ( √ k 2 -k 2 n ( d -d p ) ) cos ( √ k 2 -k 2 n (d -y) ) , in d > y > d p , (43b)
and the constant W can be again fixed by normalization. In the sequel the above eigenfunctions will be simply denoted

as Y n (y) = Y (2) n (y) /   Y (2)
n   , and since they are orthogonal it holds ⟨Y n , Y m ⟩ = δ nm , where δ nm denotes Kronecker's delta. Moreover, the corresponding eigenvalues {k n , n = 1, 2, . . .} are determined by the roots of the following equation

γ √ k 2 eff -k 2 n sin (√ k 2 eff -k 2 n d p ) cos ( √ k 2 -k 2 n ( d -d p ) ) + S √ k 2 -k 2 n sin ( √ k 2 -k 2 n ( d -d p ) ) cos (√ k 2 eff -k 2 n d p ) = 0. (44)
In Fig. 9 the first four vertical eigenfunctions Y n (y), in 0 < y < d = 1.3 m, are plotted, for frequency f = 1.8 Hz and water depth h = 0.23 m. In particular three cases of porosity coefficients are examined: (a) γ = S = 1 in uniform layer (e.g., in water regions I and III), (b) for γ = 0.67 and S = 1/γ , and (c) γ = 0.67 and S = 1. In all cases the interface is located at

Comparison of simplified model against the modal-BEM solution and experimental data

In order to illustrate the applicability of the present simplified model for evaluating the performance of porous breakwater structures composed by dense arrangements of vertical columns, in this section we present comparisons with solutions obtained by the modal-BEM and laboratory data. In particular, the experimental arrangement tested in the SeaTech wave tank, in water depth 0.23 m, shown in Fig. 5 By comparing the plots of Fig. 10 with the corresponding ones in the left column of Fig. 6, we conclude that the calculated wave patterns are compatible and thus, the present simplified model is able to provide useful information concerning the wave details in the water region outside the porous structure in the waveguide, and especially, its reflection and transmission properties. This is further indicated by considering the distributions of the calculated wave amplitude along the longitudinal section y = 0.1 m (at a distance 10 cm from the centerline of the tank), for the two frequencies f = 1 Hz and f = 1.8 Hz, which are presented in Fig. 11. In this figure results obtained by the present simplified model (shown by using solid lines) are plotted against the full modal-BEM solution (shown by using thick red lines), and are further compared with the experimental data indicated by using circles. In general, the results compare well with the full solution, especially at the front side of the porous breakwater dominated by reflection. An exception is observed in the case of the higher frequency f = 1.8 Hz in the downwave region just after the porous structure, where the reduced-order model locally overestimates the transmitted wave amplitude.

For the studied configuration, it is worth mentioning here that the results of the present modal-BEM are more accurate, however, they have been obtained by using a discretization of Nel = 4668 boundary elements. Compared to the modal-BEM data, the dimension of the coupled-mode system is only 2.8% and the required computation time is only 0.01%. Although the results of the coupled-mode simplified model are not as accurate, the compression of the computation requirement supports its systematic application to investigate the performance of the porous structure and further design and optimization studies of various arrangements of many such structures in more realistic environment.

Conclusions

In this work a numerical scheme based on modal expansions and Boundary Element Method is developed for treating the problem of wave scattering by dense array of multiple vertical cylinders, extended from the top to the bottom of the water column. The vertical cylinders are regularly arranged in a block and tested in a tank, modeling a porous segmented breakwater. Numerical results are shown and discussed concerning the wave details both inside and outside the arrangement, and the applicability of such structures operating as breakwater units. Focusing on the reflection and transmission properties and on the scattered field outside the porous structure, a simplified model is also derived based on an adaptation of Foldy-Lax theory. The latter provides us an equivalent index of refraction for the wave scattering by many small obstacles located in a part of the waveguide. Results obtained by the present models are compared against experimental measurements collected in the wave tank, showing relatively good agreement. The reduced-order model permits a rapid and efficient calculation of the reflection and transmission properties of the examined structures by significantly reducing the computational cost, and further supports design and optimization studies. Extensions to include non-symmetric fields in the waveguide, which can be excited by non symmetric porous structures and/or arbitrary cylinder arrangements can be straightforward developed. Finally, the present analysis may find useful applications to similar scattering problems in acoustic and electromagnetic waveguides.

Acknowledgment

Support of Kostas Belibassakis to visit Université de Toulon is greatly acknowledged.

Appendix. An approximate model of the porous medium

In order to derive an approximate model for the scattering of time-harmonic waves by large number of densely arranged small obstacles in a homogeneous medium, we consider here the work by Foldy [START_REF] Foldy | The multiple scattering of waves. I. General theory of isotropic scattering by randomly distributed scatterers[END_REF] concerning multiple isotropic scattering in a medium which contains randomly distributed small scatterers, which was further extended by Lax [START_REF] Lax | Multiple scattering of waves[END_REF][START_REF] Lax | Multiple scattering of waves. II. The effective field in dense systems[END_REF]. Foldy-Lax theory is based on the assumption that the scattered field is approximated by the waves emitted by point sources placed at the centers of the obstacles, with density proportional to the wave field (see also [START_REF] Martin | Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles[END_REF] Chap. 8). We start here by considering the following decomposition of the total wave field into incident ϕ inc (x) = ϕ inc (x, y) and scattering field ϕ scat (x, y) ϕ (x, y) = ϕ inc (x, y) + ϕ scat (x, y) , where ϕ scat (x, y) = qn 0 ∫ Ω ϕ (x s ) H (1) 0 (k |x -x s |) dx s ,

Fig. 1 .

 1 Fig. 1. Schematic presentation of the waveguide containing a porous scatterer constituted by many vertical cylinders in the middle. (a) Top view, (b) vertical section through the channel centerline.

Fig. 2 .

 2 Fig.2. Hybrid modal-BEM formulation and discretization of various sections of the boundary in the middle subdomain of the waveguide. Side walls are the boundary parts 1 and 3, and the incident and transmission interfaces are the parts 4 and 2, respectively. The cylindrical scatterers are numbered as parts 5, 6, . . . , M. Also, the outward normal vector on various parts of the boundary is indicated.

( 21 )

 21 Using the orthogonality and completeness properties of the basis functions {Y n (y) , n = 1, 2, . . .} in the interval -d ≤ y ≤ d, the coefficients { A (1)n , n = 1, 2 . . .}are expressed from Eq. (20) as follows

Fig. 3 .

 3 Fig. 3. (a) Waves scattered by a hard cylinder (ka = 4.66) as calculated by the present BEM using N el = 31 elements. (b) Real part of the potential on the surface of the cylinder using N el = 31 (circles), N el = 61 (crosses) against the analytical solution (solid lines). (c) Same as before for the imaginary part of the potential.

Fig. 4 .

 4 Fig. 4. Porous breakwater constructed by regular arrangement of vertical cylinders.

  2 tan -1 (y/(x + r)) are the cylindrical coordinates and J ′ m (u) = dJm(u) du and Φ(1) m (u) = dH (1) m (u) du , and ε m is the Neumann symbol (ε m = 1, m = 0, and ε m = 2, m ≥ 1). Indicative results obtained by the present BEM are shown in Fig.

Fig. 5 .

 5 Fig. 5. Tests in the wave tank of Seatech (Univ. Toulon) with the porous breakwater unit located in the middle of the tank.

Fig. 6 .

 6 Fig. 6. Calculated wave field (only half symmetric part) in the waveguide modeling the tank experiments in water depth h = 0.23 m, at frequencies (a) f = 1 Hz, (b) f = 1.5 Hz, (c) f = 1.8 Hz, excited by the first incident mode. The right column subplots indicate computed results along the section y = 0.1 m from the centerline. Real part of the field is indicated by thin solid line, imaginary part by dashed line, and modulus with thick lines, respectively. The limits of the structure are indicated by vertical dashed lines.

Fig. 7 .

 7 Fig. 7. Transmission coefficients (T n , shown by using black lines) and reflection coefficients (R n , shown by using red lines) of the examined porous structure against nondimensional wavenumber kh, concerning the first three modes. Results for the first mode are indicated by using solid lines and for the rest by dashed lines, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .

 8 Fig. 8. Domain decomposition of the plane waveguide including the porous structure (porous region) in the middle part shown by thick dashed lines.

Fig. 9 .

 9 Fig. 9. First four vertical eigenfunctions Y n (y), in 0 < y < d = 1.3 m, in the case of f = 1.8 Hz, h = 0.23 m, for: (a) γ = S = 1 in uniform layer (e.g., in water regions I and III), (b) in the middle subregion (II), for γ = 0.67 and S = 1/γ , and (c) γ = 0.67 and S = 1, with the interface located at y = d p = 0.6 m, as indicated by horizontal dashed line. The lower part (P) denotes porous medium and the upper part (W) water, respectively.

  dY

  y = 0 and y = d,

  is considered, for regular incident waves with frequency in the range 1 Hz-1.8 Hz. Results obtained by the simplified model using γ = 0.67, corresponding to the porosity of the structure of Fig. 4, are presented in Fig. 10. More specifically, in the top subplots the real part of the calculated wave field is shown and in the lower subplot the modulus of the field, respectively, for the three plane incident wave frequencies: (a) 1 Hz, (b) 1.5 Hz and (c) 1.8 Hz. The colorplots concerning the modulus of the wave field are also representative of the spatial distribution of the amplitude of the free surface elevation in the domain normalized with respect to the incident wave. The numerical results shown have been obtained by the simplified model, Eqs. (25), using totally 25 modes.

Fig. 10 .

 10 Fig. 10. Calculated wave field as obtained by the present simplified model in water depth h = 0.23 m, for three frequencies. For each frequency, in the top subplots the real part of the field is shown and in the lower subplots its modulus (or free surface elevation amplitude), respectively. The porous subregion is indicated by dashed lines.

Table 1

 1 Critical frequencies of Seatech tank, for water depth h = 0.23 m and distance between centerline and side walls d = 1.3 m.

	n	1	2	3	4	5	6
	f n (Hz)	0.55	0.98	1.29	1.53	1.72	1.89

y = d p = 0.6 m, and is indicated by using horizontal dashed lines. In the last two cases of Fig. 9(b) and (c) the discontinuities concerning the eigenfunctions and its y-derivative are evident at the interface separating the lower part corresponding to the porous (P) medium and the upper part corresponding the water (W) region, respectively. Consequently, in the present simplified model, the problem is solved by finding the modal amplitudes ϕ n (x) in order to satisfy Helmholtz equation in subregion (II) ∇ 2 ϕ (2) + k 2 (y) ϕ (2) = 0, (45) in conjunction with the following matching conditions ∂ϕ (1) ∂x = N (y) ∂ϕ (2) ∂x and ϕ (1) = F (y) ϕ (2)

∂ϕ (3) ∂x = N (y) ∂ϕ (2) ∂x and ϕ (1) = F (y) ϕ (2)

where ϕ (1) and ϕ (3) are given by the modal expansions, Eqs. ( 14), in subregions (I) and (III), respectively. Also, the functions N (y) and F (y) are determined from parameters γ and S as follows:

and

and

Subsequently, by projecting Eq. (45), at each position x, on the transverse basis

the following system of horizontal equations is derived with respect to the (unknown) modal amplitudes

where a prime denotes differentiation with respect to x, and the coefficients of the system are defined by

Moreover, the following boundary conditions are derived from the matching conditions, Eqs. (46),

where the coupling matrices A mn and C mn are defined by

The above problem can be semi-analytically treated by introducing the following decomposition of the solution into incident and reflected waves

where k n are obtained from Eq. (44). Using the above equation and the corresponding one concerning the x-derivative

the system of equations, Eq. (48), is automatically satisfied and the boundary-matching conditions, Eqs. (50), result in the following algebraic system with respect to the unknown amplitudes α n and β n :

where E n = exp (2ik n ℓ). The system equations ( 54) is finally solved by truncating the series equation (41) keeping a number of terms sufficient for the convergence of the results. which is a modification of the version based on discrete point sources for ϕ scat (x, y) (as e.g., described in [16], Eq. 8.50) to a continuous distribution of the same singularities. In the above equation Ω denotes the region occupied by the densely distributed scatterers. In the case of plane harmonic incident waves propagating along the x-direction in infinite medium ϕ inc (x, y) = exp (ikx). Moreover, q is a proportionality coefficient to be determined, and n 0 denotes the number of scatterers per unit volume in Ω. For the present structure, treated as an inclusion, n 0 is expressed in terms of porosity as:

. Using the properties of the fundamental solution G (x|x s ) of the Helmholtz equation on the plane (defined by Eq. ( 16b)), we have

where δ (xx s ) denotes Dirac's delta function. Subsequently, by applying the Helmholtz operator to Eq. (A.1) we obtain:

) ,

where k eff represents an effective index of refraction associated with the porous medium, which could behave as dissipative with characteristic coefficient controlled by the imaginary part of k eff .

In order to obtain an approximation for the coefficient q, we focus on the vicinity of each one scatterer -circular disc of radius a, which is assumed to be located well in the interior of the subregion Ω, and consider a cell ω of the porous medium containing the scatterer. Next, we split the integral defining the scattering field in Eq. (A.1) in two parts as follows:

Without loss of generality, we select the center of cell ω to be located at the origin of the coordinate system. Also, since only one disc is included in ω, it holds n 0 ∫ ω dx s = 1. Taking into account the singular behavior of the function H (1) 0 (k |x -x s |) for x ∈ ω, the contribution of the second term in the right-hand side of Eq.(A.4) to the scattering field is considered to be dominant. In addition, the total field under the integral is approximated by the incident one ϕ ≈ϕ inc , and thus, locally we obtain the following approximation, for x = (x, y) ∈ ω (and x s = 0), ϕ scat (x, y) ≈ ϕ inc qH (1) 0 (k |x|) = qH (1) 0 (k |x|) .

(A.5)

Comparing the above result with the leading order term in the Fourier-Bessel series (Eq. ( 32)) representing the scattering field by a single circular disc, i.e. ϕ scat (x, y) ≈ A 0 H (1) 0 (k |x|) , where A 0 = -J 1 (ka) /H (1) 1 (ka) ,

we finally obtain the approximation q ≈ A 0 .

(A.7)

Moreover, for small nondimensional wavenumbers ka ≪ 1, the coefficient A 0 becomes A 0 ≈ -iπ 4 (ka) 2 , and using this result in Eq. (A.7) we obtain q ≈ -iπ 4 (ka) 2 . Substituting the latter value in the expression for k eff given by Eq. (A.3) we finally obtain the following estimation concerning the effective index of refraction of the porous medium in terms of the porosity γ ,

(A.8)