V Rey
email: rey@univ-tln.fr

G Arnaud

J Touboul

K Belibassakis

Water wave scattering by dense or sparse arrays of surface-piercing bodies by integral matching method

Keywords:

In this work, integral matching method is used for the numerical study of water wave scattering by arrays of surface-piercing structures. Such a method, already used for many years for rectangular breakwaters, is applied to structures of various shapes, including arrays of vertical cylinders, for which the contour of vertical boundaries is discretized into steps. Finite width channels of constant depth are considered, which may also correspond to uneven periodic structures along the transverse direction. Applications to sparse or dense arrays of vertical cylinders are discussed. For sparse cylinders, first order and higher order Bragg resonance are observed in both the incoming wave direction and other directions for the higher frequencies. For dense cylinders, considered as porous media, dispersion relation for the wave within the cylinders arrays is discussed through interference process due to the finite length of the porous medium (or the cylinders array) along the incoming wave direction. It is found that the wavelength slightly depends on the cylinders arrangement at given porosity and specific surface.

Introduction

During the last decades, a huge effort has been made to better understand and model water waves propagation through inhomogeneous media and their interaction with ocean and coastal structures. Arrays of vertically emerging bodies have been widely studied in ocean engineering in the context of wave diffraction by cylinder arrays and wave loads on these structures (see e.g. [START_REF] Kagemoto | Interactions among multiple three-dimensional bodies in water waves: an exact algebraic method[END_REF][START_REF] Ohl | Water wave diffraction by a cylinder array[END_REF] and references cited there) and in coastal engineering for shore sheltering by detached offshore breakwaters of rectangular shape through wave scattering and dissipation [START_REF] Dalrymple | Wave diffraction through offshore breakwaters[END_REF][START_REF] Abul-Azm | Oblique wave diffraction by segmented offshore breakwaters[END_REF]. For environmental purpose, emergent wetland vegetation can be modelled by dense arrays of vertical cylinders across which wave propagates with significant damping ( [START_REF] Mei | Waves of intermediate length through an array of vertical cylinder[END_REF] and references cited). Such permeable bodies of finite extent which may also be used in coastal engineering for sheltering the shore since a part of the incoming energy may be reflected through interference process in addition to the dissipation process ( [START_REF] Mei | Long waves through emergent coastal vegetation[END_REF][START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] and references cited there for a more complete review).

For single or sparse vertical cylinders of circular cross section, analytical solutions exist and take advantage of the axial symmetry (see [START_REF] Mei | The Applied Dynamics of Ocean Surface Waves[END_REF]). For an array of cylinders, approaches from sparse cylinders for which plane wave approximation can be used (see e.g. [START_REF] Simon | Multiple scattering in arrays of axisymmetric wave-energy devices. Part 1. A matrix method using a plane-wave approximation[END_REF]) to denser arrays where interference between adjacent bodies have to be considered, various analytical or numerical solutions are proposed (see [START_REF] Walker | Wave diffraction from linear arrays of cylinders[END_REF] and references cited there). Linton and Evans [START_REF] Linton | The interaction of waves with arrays of vertical circular cylinders[END_REF] solved exactly the problem of wave scattering by arrays of vertical cylinders through the use of a simplified expression of the velocity potential in the vicinity of a cylinder. When dissipation effects are concerned for rather dense arrays of cylinders, semi-analytical solutions based on multi-scale anaysis can be used as proposed by Mei et al. [START_REF] Mei | Long waves through emergent coastal vegetation[END_REF]. In the latter works, dissipation is considered at the scale of the bodies, while both the diameter of the bodies and the distance between adjacent bodies remain small compared to wavelength. Since permeability remains large enough to allow wave to propagate across the arrays of cylinders, interference process and wave scattering may be observed, leading to diffraction depending on the cylinder diameter and spacing when wavelenth is of the order of cylinders spacings. For longer waves of wavelength of the order of the array of cylinders, reflection due to the finite extend of the array in the incoming wave direction may be observed [START_REF] Mei | Long waves through emergent coastal vegetation[END_REF][START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF]. Such interference processes are related to the wave wavelength. The dispersion relation in the presence of vertical cylinders was recently discussed by Molin et al. [START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF] using results from sloshing experiments within a tank at resonant conditions.

Vertical cylinders of rectangular cross section are often considered as an academic case to model breakwaters in coastal engineering. Wave field around a detached breakwater on a plane beach of constant slope is described in [START_REF] Horikawa | Nearshore Dynamics and Coastal Processes[END_REF]. Diffraction effects are significant in the case considered since the wave wavelength is of the order of the breakwater width. When considering breakwaters along shoreline, the gap between adjacent breakwaters is also a key parameter for the wave diffraction process, as shown by Dalrymple and Martin [START_REF] Dalrymple | Wave diffraction through offshore breakwaters[END_REF] for normal wave incidence and by Abul-Azm and Williams [START_REF] Abul-Azm | Oblique wave diffraction by segmented offshore breakwaters[END_REF] for oblique wave incidence for thin breakwaters and constant depth conditions. If breakwaters are of uniform geometry and gaps identical, the geometrical domain reduces to a channel of width equal to the periodicity of the breakwaters. Since boundaries are rectilinear, cartesian coordinates are employed, and an eigenfunction expansion method is used for the solution, based on matching conditions of continuity at the gap location for the wave, written as the sum of the incoming (or transmitted) wave and additional modes, which correspond to either scattered wave modes or evanescent modes. The number of scattered modes depends on the channel width. In the case of rectangular body within a channel, trapped modes were evidenced for a vertical cylinder of rectangular cross-section [START_REF] Evans | Trapped modes in open channels[END_REF], or for more generally for bodies of different shapes [START_REF] Newman | Trapped-wave modes of bodies in channels[END_REF]. Under the assumptions that the fluid is incompressible and inviscid, and that its motion is irrotational, a general expression for the velocity potential can then be defined in waters of constant depth in a cartesian frame, as proposed for instance by Dalrymple and Martin [START_REF] Dalrymple | Wave diffraction through offshore breakwaters[END_REF]. If wave propagates across successive rectangular domains of variable width, the solution for the potential can be obtained by means of pressure and flow matching conditions at each interface between successive domains though an integral matching method. Such a method is also well known for wave propagation in the presence of abrupt variations of the depth in the two-dimensional case. In the case of smooth bottom variations, this method is still valid for the calculation of both reflected and transmitted waves [START_REF] Guazzelli | Higher-order Bragg reflection of gravity surface waves by periodic beds[END_REF][START_REF] Rey | Propagation and local behavior of normally incident gravity waves over varying topography[END_REF] even if the velocity field is not valid near the boundaries discretized in successive steps but can be improved by use of an additive term in the potential expansion [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions[END_REF]. In the present study, integral matching method is used for the study of wave scattering by surface piercing cylinders with vertical walls and of any cross-section shape. Indeed, for smooth cross-sections shapes including circles, a discretization of the contour is done to obtain rectangular domains.

In Section 2, the theoretical model and the numerical formulation are presented. Results and comparisons to previous studies are then reported and discussed in Section 3. Concluding remarks are then given in Section 4.

Theoretical model and numerical formulation

Linear potential theory

Wave is assumed to propagate along a channel of constant depth h and variable width. Surface piercing structures can be of irregular shape in the horizontal plane, their shape does not depend on the vertical direction, and boundaries are then vertical walls as shown in Fig. 1. A coordinate systel is considered with origin at some point at the mean free surface and z axis pointing upwards. For smooth surface contours, as for instance vertical cylinders, discretization into steps will be done. For unbounded periodic structure arrangements along y-axis, conditions of symmetry will allow to further consider geometrical domains bounded along y-axis. In the following, the incoming wave propagates along x-axis, towards x > 0. Scattered waves propagate in the channel outwards the structure.

General expressions for the velocity potential within domains of finite width

In the linear problem, a general expression of the velocity potential in a cartesian frame Φ(x, y, z, t) for a regular wave of angular frequency ω propagating over a constant water depth h is of the form

= + + Φ x y z t αω k k z h F x y e ( , , , ) cosh[ ( )] sinh(kh) ( , ) cc iωt ( 1
)
where α is the amplitude of free surface elevation. Wavenumber k verifies the dispersion relation:

= ω gk tanh(kh) 2 (2)
Considering the fluid domain to be bounded along y-axis in the interval d m < y < d M , with d = d Md m the width of the channel, velocity potential (1) takes the form

= + + Φ x y z t k z h ϕ x y e ( , , , ) cosh[ ( )] ( , ) cc iωt ( 3
)
with the free surface potential ϕ(x, y) given by Since ΔΦ = 0, where

∑ = + = ∞ -- + + ϕ x y A e A e ψ y ( , ) [ ] ( ) n n x n x n 0 ik ik xn xn (4)
= ∂ ∂ ∂ ∂ ∂ ∂ ( ) Δ , , x y z 2 2 2 2 2 2 , = - k k k ( ) xn 2 yn 2 1/2 (7)
For n = 0, k xn = k, and k xn is purely imaginary for k yn > k. The number of propagating modes n prop is then given by

= + ⎡ ⎣ ⎤ ⎦ n π 1 Int kd prop (8)
with Int denoting the integer part. For propagating modes n along xaxis, k xn = k cosθ n and k yn = k sinθ n , where θ n is the angle of propagation with respect to the x-axis at given mode n. For narrow channels or long waves (

< - k π d d M m
), only the first mode propagates along x-axis, the other are evanescent along x-axis. For the numerical processing, general expression of the velocity potential is truncated at a given order n = P, the first (P + 1) modes are hence considered. For simplicity, P will be chosen the same for all the domains in the following calculations.

Equations of continuity for a single channel

For two successive domains (see Fig.

m m M M ϕ x ϕ x m m M M ϕ x m m M M 1 2 1 2 1 2 1 2 1 2 mj 1 2 1 2 Mj j 1 2 (9)
Functions ψ j,m , ψ j,n form an orthogonal system for the scalar product

∫ = f g f g . dy j j d d j j mj Mj (10)
Taking advantage of the orthogonal set of functions ψ j,m , interface conditions are written in an integral form, for any n.

for n = 0, …, P. In the following, we will only considerer domains of single channel thanks to the symmetries along y-axis. For multiple channels, one may write in the integral form pressure continuity condition along the boundary (in the y-direction) between two successive domains. The velocity condition is written along the whole channel width, through integral form taking into account both the velocity condition between successive domains and the absence of velocity normal to the boundaries. The method is similar to the method used in the two-dimensional vertical case for a submerged obstacle (see for instance, [START_REF] Rey | A note on the scattering of obliquely incident gravity waves by cylindrical obstacles in waters of finite depth[END_REF][START_REF] Rey | Forces and moment on a horizontal plate due to regular and irregular waves in the presence of current[END_REF]). Similar method has been used by Belibassakis et al. [START_REF] Belibassakis | Three-dimensional wave diffraction in the vicinity of openings in coastal structures[END_REF] to study tridimensional diffraction in the vicinity of openings in coastal structures.

∑ = + < < = -- - + + - + + ϕ x y A e A e ψ y x x x ( , ) [ ] ( ) for j n P n j x x n j x x n j j j ( ) 0
( ) ik ( ) ( ) ik ( ) ( ) 1 j j j xjn xj ( ) 1 (15)
where ± A n j ( ) are the 2(P + 1) unknown complex coefficients for each bounded domain D j . The 2N(P + 1) unknown coefficients are then numerically solved thanks to the 2N(P + 1) matching conditions written at each abscissa x j , j = 1, …, N derived from Eq. ( 13), ( 14) and (15).

Energy conservation

The mean energy flux across y0z plane is given by:

∫ ∫ ∫ = → → + =- E T pv n 1 . dtdydz t t t T y y x h η m M ( 16
)
where → n is a unit vector normal to y0z and = -

∂ ∂ p ρ Φ t , → → = ∂ ∂ v n . Φ x .
In the absence of dissipation, conservation of wave energy flux along the direction of propagation 0x takes the form:

∑ = + + ⎡ ⎣ ⎢ + ⎤ ⎦ ⎥ - + + = - + + k k A k k A k k A k k A k k A | | | | | | 1 2 | | | | x x x M M n n n M n M 0 (0) 0 (0) 2 0 (0) 0 (0) 2 0 ( ) 0 ( ) 2 1 1 xn (0) (0) 2 xn ( ) ( ) 2 prop (17)
Defining the following reflection and transmission coefficients for modes n [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions[END_REF] and assuming the same values for the boundaries y m and y M for both incident wave and transmitted wave domains, expression [START_REF] Rey | Propagation and local behavior of normally incident gravity waves over varying topography[END_REF] becomes,

= = + - + - R A A T A A | | a n d | | n n n n M (0) 0 (0) ( ) 0 (0)
∑ = + + + = - R T k k R T 1 1 2 [ ] n n n n 0 2 0 2 1 1 xn (0) 2 2 prop ( 19
)
which reduces to the classical energy conservation condition

= + R T 1 0 2 0
2 when only the first mode (n = 0) is propagative. Reflected energy coefficient for the first mode is defined as the ratio between the reflected and incident energy, it is given by

= E R R0 0 2 (20)
Reflected energy coefficient for mode n > 1 is , the reflected coefficient R n defined in ( 18) can be higher than 1. The total reflected energy coefficient E R is given by

= E R k k E 1 2 n I Rn 2 xn (0) (21)
∑ = = - E E R n n 0 1 Rn prop (22)
Energy flux conservation (Eq. ( 17)) is used as a control parameter for the numerical computations. The relative error is found under 10 -5 whatever the conditions.

Results and comparisons

The method can be used for vertical cylinders of any section shape. In the following, examples are given which correspond to previous studies.

Single body in a channel of finite width

In this section, we consider the case of a single body in a channel of finite width. Let us note that this configuration also corresponds to the case of periodic bodies normal to the incoming wave direction, as considered in the following subsection.

Breakwater of rectangular shape

We consider here the case of offshore breakwaters as presented by Dalrymple and Martin [START_REF] Dalrymple | Wave diffraction through offshore breakwaters[END_REF]. In their paper, breakwaters are considered as infinitely thin vertical barriers. Influence of the length of such structures was studied recently by Zhu and Xi [START_REF] Zhu | Hydrodynamic characteristics of offshore and pile breakwaters[END_REF]. Such an assumption allows them to compare variational methods and eigenfunction methods when only the normally incident waves are progressive. A schematic view of one spatial period along y-axis of the breakwaters is presented in Fig. 3. The gap between adjacent breakwaters is 2l, the spatial periodicity which also corresponds to the channel width is b. The width along x-axis of the breakwaters is Δx. This length Δx is considered as infinitely small in [START_REF] Dalrymple | Wave diffraction through offshore breakwaters[END_REF], while in the present work, we will use for calculations Δx = b/100. Results concerning the amplitude of reflection of each of the 4 first modes and the reflected energy either total or for the first mode are presented in Fig. 4 versus dimensionless width kb. Conditions considered here are l/b = 0.25 (case presented in [START_REF] Dalrymple | Wave diffraction through offshore breakwaters[END_REF], Fig. 7 for the first mode), water depth h = 0.3 m and the width b, which corresponds to half the periodicity of the breakwaters along y-axis, is 1 m. We can observe that for kb < π, only the first mode (n = 0) is propagative. For kb ⩾ nπ, mode n becomes propagative. Reflected energy E R0 corresponding to only the first mode (normal to x-axis) and the total reflected energy E R are presented respectively in dashed dot lines and dashed lines. After an almost full reflection for kb = π, reflected energy decreases and increases again up to kb = 2π. Smaller variations are observed for higher values of kb. In fact, the main part of the reflected energy flux is due to the first mode, whatever the frequency in the present configuration. When kb ≃ nπ, huge oscillation may appear in the transverse direction. We indeed observe amplitudes above one for modes 2-4. These oscillation amplitudes however remain finite, of relative amplitude 2.52, 22.68 and 1.89 with respect to the amplitude of the incoming wave, respectively for the second, the third and the fourth modes (values above the vertical axis upper limit in Fig. 4). Even if the reflected amplitude is locally huge for kb ≃ nπ, the corresponding energy flux remains much lower according to the wavenumber component along x-axis (see Eq. ( 21)). These huge surface amplitude oscillations for trapped waves at or near resonance are well known and studied for multiple legs offshore platforms (see e.g. [START_REF] Ohl | Water wave diffraction by a cylinder array[END_REF]).

Dimensionless surface wave amplitude fields with respect to the incident wave amplitude are presented in Fig. 5 4 by vertical dashed lines. We can observe a very low transmission for all four cases since in this frequency range, near total reflection is observed (see Fig. 4). For kb -π =-0.1, a single plane wave is observed far upwave the breakwater whereas the far field incident field is the superposition of modes 1 and 2 for kb -π = 0.1. For kb -π =-10 -8 , only mode 1 is propagative whereas the two first modes are propagative for kb -π = 10 -8 . However, we can observe similar surface wave fields since there is a smooth transition between an evanescent mode for mode 2 tending to zero (almost no decrease along x < 0), to a propagating mode in a direction closed to normal with respect to the incident wave. In both two cases, we can observe huge oscillations forced by the reasoning wave conditions along the y-axis, due to the periodicity of the breakwaters. All the findings are also observed in Fig. 6(a) for y = 0 and in Fig Figure 6(b) for y = 0.99, near the boundary or axis of symmetry.

Vertical cylinder of circular section

For cylinders of circular section, boundaries are discretized into steps. For simplicity, the discretization length along x-axis is kept constant (see Fig. 7). Tests of convergence with respect to both the number of modes n modes and the number of steps of discretization n steps of the cylinders are presented in Figs. three first modes are propagative in this configuration. A colourplot of the surface wave amplitude field is given in Fig. 8. We can observe a rapid convergence of the wave behaviour, especially in the far field. This rapid convergence is also observed in Fig. 9 through the behaviour of the wave amplitude along two different transects, (a) at y = 0 which corresponds to the axis of the vertical cylinder, and (b) at y = 0.6 which corresponds to a section next to the cylinder position. We can observe that even if the discretization of the cylinder is rough, far field effects are already well reproduced for n steps = 5. In order to better quantify this convergence, the relative error of the wave amplitude in the near field and of the reflection coefficient for the three propagating modes, representing the far field behaviour, are presented respectively in Fig. 10a andb. We define here the relative error ER i of a quantity φ i by

= - φ φ φ ER | | i n i n ( 23
)
where index i corresponds to a given number of steps n steps for the calculation, and n to the maximum number of steps considered for the test, here, n = 10. We can observe that a relative good convergence is reached for n steps = 5, with globally smaller relative errors for the reflection coefficients. Such a rapid convergence for the far field was also observed and discussed in [START_REF] Rey | Propagation and local behavior of normally incident gravity waves over varying topography[END_REF] for wave propagating over varying bottoms discretized into successive steps, even if perfect continuity for the horizontal velocity is hardly reached at the steps for steep slopes even when increasing the number of evanescent modes. In the present case, modes n, n > 1 are either propagating or evanescent, convergence is rapidly reached as shown in Figs. 8 and Figure 9. In fact practically a compromise has to be done between the precision of the results and the size of the system to be solved numerically. For far field calculations as the reflected coefficient, we verify that a good convergence of the results is reached at a given n modes and n steps .

Periodic cylinder arrays

In this section, examples are given for arrangements of cylinders of circular section as already studied in previous works. We consider hereafter the cases of sparse (a/L ≪ 1) or dense (a/L of order 1) arrays, where a = D/2 and L are the cylinder radius and the distance between the centres of adjacent rows of cylinders, respectively (see Fig. 11). Also, W is the width of the channel (transverse periodicity of the arrangement). Reflection coefficient is discussed versus frequency or dimensionless wavenumber.

Sparse array

In this first section, we consider the case of sparse cylinders as 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 their study, only the first mode is propagative in the range of frequency considered. We extend here the study to higher frequencies, for which scattering is observed not only for the reflected wave along the x-axis, direction of propagation of the incident wave, but also for scattered waves (oblique waves). We will consider here the same configuration as Fig. 9 of [START_REF] Li | Bragg scattering by a line array of small cylinders in a waveguide. Part 1. Linear aspects[END_REF], where they compare results from their asymptotic solution with those from a finite element method, for a/L = 0.10, L/W = 1, N = 10, a, L, W and (2N + 1) being the cylinder radius, the distance between two successive cylinders along the x-axis, the channel width and the number of cylinders along x-axis, respectively. A sketch of the configuration is shown in Fig. 11. Results with the present method concerning the reflected wave are calculated for N = 5 (a row of 11 cylinders) with a discretization step n steps = 5 for each of the cylinders. Results concerning the reflection coefficient corresponding to each of the propagating modes are presented in Fig. 12 for a truncation of the series at either the fourth or fifth mode. For the higher values of the frequency range (kL/π = 5), three modes are propagative. Convergence is more rapid at lower frequencies, for which only one mode is propagative. We can observe that convergence is globally reached for P = 4 (five modes) on the range of frequency considered. We can notice that since the channel width W = L, only the first mode (n = 0) is propagative for kL < π. As in [START_REF] Li | Bragg scattering by a line array of small cylinders in a waveguide. Part 1. Linear aspects[END_REF], we can observe for mode 1 a strong reflection near kL = π and kL = 2π. These two peaks correspond to respectively 1st and 2nd order Bragg resonances which correspond to resonant wavemedium interactions when the spatial periodicity of the medium is half the wave wavelength at first order (here kL = π), or a multiple of half the wave wavelength for higher orders (here kL = pπ, p > 1). Let us note that these peaks are slightly shifted towards low wavenumbers, due to the finite size of the cylinders (see also [START_REF] Li | Bragg scattering by a line array of small cylinders in a waveguide. Part 1. Linear aspects[END_REF]). High values of the reflection coefficient above one is observed for mode 2 and for mode 3, for kL = 2π and kL = 4π respectively. As for the case of the breakwaters, they correspond to transverse wave resonances due to the periodicity of the array of cylinders along y-axis, since the width W = L in the present case. Periodic shorter and smaller oscillations with respect to kL = π are also observed which correspond to the finite length of the row of cylinders along x-axis. Indeed, if we consider the length of the patch of cylinders, L row = (2N + 1)L, conditions of resonance are kL/π = p/(2N + 1), p > 1.

For kL/π ⩽ 1, 2N + 1 values of extrema are then found. Let us note that this number more or less corresponds to Fig. 9 of [START_REF] Li | Bragg scattering by a line array of small cylinders in a waveguide. Part 1. Linear aspects[END_REF], which however refers to 21 cylinders (N = 11) instead of 11 (N = 5) in the present case. Concerning modes 2 and 3, high reflection coefficients (amplitude with respect to the amplitude of the incoming mode 1, see Eq. ( 18)) are also observed. They are also due to resonance conditions for these modes. Since these modes are generated within the array of cylinders, resonance conditions, which lead to a strong reflection, correspond to opposite phase conditions at successive cylinders, for k xn , the component of k n along x-axis. k xn L = (2p + 1)π/2, p > 1 where n = 1 for mode 2 and n = 2 for mode 3. Reflection for modes 2 and 3 are presented in Fig. 13a andb respectively. Total energy reflected and energy reflection due to modes 1-3 (n = 0-2) are also presented. We can observe that the main energy reflected is due to these oblique modes at their resonance conditions. Near full wave reflection, especially through mode 2, is observed for k x1 L/π of about 1.5 (p = 1) and k x1 L/π of about 2.5 (p = 2). Let us note that the shift of the resonant peak towards lower wavenumbers increases as the wavenumber increases, which corresponds to an increase of the relative size of the cylinders with respect to the wave wavelength. Contrary to those of the resonant peaks for modes 1 and 2, the shapes of the resonant peaks for mode 3 exhibit a more irregular shape at resonance. We can observe for instance a huge reflection coefficient for mode 2 just under k x2 L/ π = 1.5, and a significant reflection for mode 1 above. The asymmetric shape of the reflection coefficient of the third mode near resonance may then be related to the evolution of the other modes near this value. We also observe in Fig. 13b that a sudden drop of the reflection coefficient for the third mode occurs in the neighbourhood of k x2 L/π = 1.5. Extensive numerical calculation has shown that this result is due to the convergence of the present method, which appears to become poor in the close vicinity of the above frequency, where the condition of the matrix system becomes near singular.

Dense array

We consider here the case of dense arrays of cylinders as considered in the framework of the thesis of [START_REF] Arnaud | Houle à la côte: Propagation, impacts et ouvrages innovants[END_REF] and presented in a two dimensional wave tank in [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF]. In this 2D study, dense array of staggered circular cylinders is considered as a porous medium, with porosity enough to allow the wave to propagate through and to dissipate slowly. Effects of both porosity and specific surface are discussed in the above paper. We here only consider wave propagation across the dense array of cylinders, neglecting the damping effects. The aim is to discuss the wave properties in terms of its wavelength (or celerity) within the porous media through the calculation of the reflected coefficient. Because of the finite length of the porous medium, interference process leads to a oscillating behaviour of the reflection coefficient versus frequency, as shown in [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF]. Since reflection was measured versus frequency, ranging from shallow water to deep water conditions, a wide range of wave conditions were then considered on a large range of wavelengths. Other types of experiments carried out at resonance for a given wavelength can also be considered to characterize the wave wavelength, as done recently by Molin et al. [START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF]. In [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF], the velocity potentials are calculated through integral matching by considering the whole arrays of cylinders as a porous medium. The dispersion relation is

= ω S gk tanh (kh) 2 (24)
where

= + - S C γ γ 1 (1 ) M ( 25
)
where S is a dimensionless reactance and C M is an added mass coefficient, C M = 1 for a single cylinder. This coefficient is chosen as a adjustable coefficient in [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF]. It is there fixed to C M = 0.29, in order to fit the experimental oscillating behaviour of the reflected coefficient versus frequency. Wave dissipation is taken into account through a linear or quadratic decay of the wave amplitude across the porous medium. For a linear decay, the complex wave number k p within the porous medium is of the form = + k k n (1 ik/ ) p w . n w corresponds to an attenuation rate per wavelength unit, it does not depend on frequency.

The dispersion relation ( 24) is a simplified form of

= ω Z igk tanh(kh) 2 (26)
proposed by Yu and Chwang [START_REF] Yu | Wave motion through porous structures[END_REF], where Z = f R + iS is the dimensionless impedance of the porous medium, f R being a dimensionless resistance. In this latter formulation, wave decay is directly connected to the wave wavenumber through the complex dispersion relation. These different approaches are discussed in [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF].

The array of cylinders considered in [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] consists of rows of cylinders with an angle of π/4 with the incoming wave direction. We here consider the case of cylinders of diameter D = 0.05 m, with 20 cylinders along x-axis. The distance between two successive cylinders along xaxis is Δx = 0.01 m, distance along y-axis is Δy = 0.064 m. This corresponds to a porous length L p = 1.20 m, its porosity γ = 0.7. The water depth h = 0.23 m. Taking advantage of the symmetry axes, the resolution with the present model can be done between two lateral boundaries, as shown in Fig. 14. Results in terms of reflection coefficient versus frequency are presented in Fig. 15a, considering either 4 or 5 modes (P = 3 or P = 4). Let us note that in the present configuration, only the first mode is propagative on the whole frequency range considered. Cylinders are discretized into n steps = 5. If we already observed a good general trend for P = 2 (not shown), convergence is assumed to be reached for P = 4, as shown in Fig. 15b, through analysis of the relative error for the reflection coefficient versus the number of modes at either f = 0.5 Hz (shallow water depth, kh = 0.50), f = 1.0 Hz (finite water depth, kh = 0.10), f = 1.5 Hz (finite water depth, kh = 2.14), f = 2.0 Hz (deep water, kh = 3.70). Both experimental results and the numerical results from Arnaud et al. [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] either taking into account dissipation with = n 25 w , as presented in Fig. 4, bottom, of their paper, or neglecting dissipation ( ⟶∞ n w ) are also presented in Fig. 15a. The semi-empirical method of [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] based on a fit of the added mass C M = 0.29 was found to fit correctly with the experiments on the whole frequency range. In addition, we observe that neglecting dissipation does not affect the location of the minima of reflection for this weakly dissipative medium case. The location of the minima of reflection for the present model sensibly fits those of [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] on the whole frequency range.

Neglecting dissipation, various formulations for the dispersion relation have been proposed (see [START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF] and references cited there). For vertical cylinders, the porous medium is considered as an anisotropic medium, derivations in [START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF] lead to the following dispersion relation:

⎜ ⎟ = ⎛ ⎝ ⎞ ⎠ ω S S gk tanh kh 2 (27)
which tends to [START_REF] Guyon | Physical Hydrodynamics[END_REF] in shallow water conditions. This formulation is found to provide a good agreement with sloshing test which allow determination of the frequency at resonance. In order to discuss the wave equivalent wavelength through the dense array of cylinders, reflection coefficient versus the dimensionless frequency kL p is presented in Fig. 16. Results of the present method are compared against predictions based on the use of dispersion relation [START_REF] Arnaud | Houle à la côte: Propagation, impacts et ouvrages innovants[END_REF] with respectively C M = 0. For both the methods of [START_REF] Yu | Wave motion through porous structures[END_REF][START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF] using C M = 1, which corresponds to the added mass of a single cylinder, we can observe that the minima of reflection are significantly shifted toward lower frequencies. This shift, which is the higher for the method of [START_REF] Yu | Wave motion through porous structures[END_REF], significates implicitly that the wave wavelenth (or its celerity) within the porous media is not accurately modelled, especially for the higher frequencies.

Other calculations with the present method where carried out with rows of cylinders along the x-axis, but with both the same porosity and the same specific surface (same number and size of cylinders). Results are presented in Fig. 17 for either aligned or staggered cylinders. Calculations are done with either one mode (P = 0) or five modes (P = 4). While considering only the first mode (P = 0), reflection coefficient does not depend on the array arrangement (aligned or staggered cylinders). Physically, higher modes are to be considered. We can consider here that convergence is reached for P = 4 as shown in Fig. 15.

We can then observe a difference in the location of the successive minima or maxima of reflection coefficient, which means that the wave equivalent wavelength within the porous medium depends on the arrangement of the cylinders. For staggered cylinders, minima of reflection are observed at lower frequency (higher wavelength) with respect to those for aligned cylinders. Arnaud et al. [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] had pointed out the influence of the porosity on the wave reflection, but had found that the specific surface did not play a significant role in the interference process, its role being essentially important on the wave damping. The cylinder array arrangement was the same in that work for the three cylinder diameters considered. We can deduce from the present results that the arrangement of the cylinders has a significant impact on the equivalent wave wavelength within the porous media, it plays in fact a role which can be compared to the tortuosity defined classically in mechanics for porous flows to characterize the permeability in addition to the porosity and the specific surface [START_REF] Guyon | Physical Hydrodynamics[END_REF]. The tortuosity τ = L′/L p where L′ corresponds to the tortuous trajectory of the fluid within the porous medium. In the present case, the tortuosity is then higher for staggered cylinders than for aligned cylinders.

Discussion and conclusion

Integral matching method is used for the numerical study of water wave scattering by arrays of surface-piercing structures. The method, already used for many years for rectangular breakwaters (see [START_REF] Dalrymple | Wave diffraction through offshore breakwaters[END_REF]), is here applied to various cases including sparse or dense surface-piercing structures with vertical walls. When applied to wave scattering by periodic arrays of cylinders, from long wave conditions (or narrow channels) to short waves (compared to the structure size or spacing), it is shown than far field effects as the reflection coefficient generally rapidly converge with both the number of the wave modes and the number of steps of discretization in the case of vertical structures with smooth contours. Convergence may however be checked carefully, especially when the matrix system becomes near singular as discussed at the end of Section 3.2.1 for resonance conditions for the third propagating mode. For long waves across cylinders arrays, we observe huge reflections at Bragg resonance as already discussed 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 incoming wave direction but also for other directions, according to resonance conditions for the oblique higher wave modes within the cylinders area. Resonance conditions also appear normally to the incoming wave direction when the channel width or the periodicity of the structures in this direction is a multiple of half the wavelength. For such corresponding frequencies, huge oscillations may be observed in a frequency band covering both parts of this resonance condition. Under the critical frequency, this oscillation is due to the evanescent studied by Dalrymple and Martin [START_REF] Dalrymple | Wave diffraction through offshore breakwaters[END_REF] for the first mode. When considering dense cylinder arrays, the structure can be considered as a porous medium (see [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF]). In this case, the obtained results fit with experiments of [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] in terms of the oscillatory behaviour of the reflection coefficient versus frequency. In the present method, the potential is calculated at the scale of the cylinders, contrary to the methods used in [START_REF] Yu | Wave motion through porous structures[END_REF][START_REF] Molin | On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders[END_REF] or [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] for which the porous medium is considered as an homogeneous medium at the scale of the wave, for which a wave wavenumber can be defined by a dispersion relation. Calculation of the reflection coefficient for another cylinder arrangement gives slightly different results, with a shift of minima and maxima of reflection versus frequency. In the experiments of [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF], three different porous media of identical porosity and different specific surface were considered. The cylinders arrangement was the same for the three media. Since the dispersion relation with the same fitted added mass coefficient C M = 0.29 for the whole cases was found to allow a good fit between model and experiments, Arnaud et al. [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] concluded that this coefficient did not depend on the specific surface for a given porosity and shape of the cylinders array. We observed here that this added mass however depends on the cylinders arrangement. This shows that a third permeability parameter, the tortuosity, has to be defined to fully define the porous properties in terms of permeability. The right equivalent wavelength within the porous medium has to take into account the cylinders arrangement. The present study focused on particular cases as described in the literature but is also applicable to structures with vertical walls of any shape in the horizontal plane. More complex structure geometries can then be considered, either for coastal or offshore purposes.

As mentioned in the introduction, dissipation effects may become significant for rather dense arrays of cylinders. In the present study, the theoretical model assumed perfect fluid and focused on wave scattering, neglecting damping effects. If one considers that dissipation remains small enough, dissipation effects can be taken into account through an exponentially decrease of the wave amplitude within the dense array of emerging vertical cylinders. Such a approach was used in [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays: interference process and specific surface effects on damping[END_REF] for the model used in their study, it is in progress for the model used in the present study.

  flux conditions are required at y = d m and y = d M .

  2) labelled 1 and 2, of respective widths d 1 = d M1d m1 (d m1 < y < d M1 ) and d 2 = d M2d m2 (d m2 < y < d M2 ), velocity and pressure matching conditions at their interface x = x 0 are the following:

Fig. 1 .

 1 Fig. 1. Sketch of the studied domains.

  In the case d m1 ⩽ d m2 and d M2 ⩽ d M1 depicted in Fig. 2, conditions are given by

1 .

 1 Construction of the matricial system Let us consider M domains, separated by N = M -1 discontinuities along the x-axis, x j , j = 0, …, N -1. Domain D 0 is defined for x < 0 and ⩽ , domain D M is for x ⩾ 0 and ⩽ superscript (j) refers to domain j. The discontinuities delimit (N -1) segments S j of length L j = x j+1x j . For each segment j, domains labelled D j are defined, bounded along y-axis by = y d m j ( ) and = y d M j For domain D 0 , the absence of incident wave of oblique incidence or the absence of divergency of the evanescent modes when x → -∞ leads to the following general expression for the free surface potential: is the coefficient for the incoming wave, -A n (0) , n = 0, …, (n prop -1) the coefficients for the scattered waves (n = 0 corresponds to the reflected wave, opposite to the incoming wave), and -A n (0) , n > (n prop -1) the coefficients for the evanescent waves. For domain D M , the absence of reflected wave of any incidence or the absence of divergency of the evanescent modes when x → ∞ leads to the following general expression n = 0, …, (n prop -1) are the coefficients for the scattered waves (n = 0 corresponds to the transmitted wave along the incoming wave direction), and -A n M ( ) , n > (n prop -1) the coefficients for the evanescent waves. (P + 1) complex coefficients are then to be calculated for each of these semi-infinite domains. General expressions for the free surface potentials are the following for domains D j , j = 1, …, N -1:

Fig. 2 .

 2 Fig. 2. Two successive domains of different width.

  for four wave periods around kb = π, respectively (a) kb -π =-0.1, (b) kb -π =-10 -8 , (c) kb -π = 10 -8 , (d) kb -π = 0.1. Only the first mode is propagative for cases (a) and (b) while the second mode is also propagative for cases (c) and (d). Both are represented in Fig.

Fig. 3 .

 3 Fig. 3. Sketch of half of the breakwater.

Fig. 4 .

 4 Fig. 4. Reflection coefficient versus kb: (-) first mode (n = 0); ( ) second mode (n = 1); ( ) third mode (n = 2); ( ) fourth mode (n = 3); relative reflected energy flux versus kb: (-.-) first mode (n = 0); (--) total energy. Vertical dashed lines correspond to the locations of the frequency cases presented in Fig. 5.

Fig. 5 .

 5 Fig. 5. Relative amplitude field relative to the incoming wave amplitude, (a) kb -π =-0.1, (b) kb -π =-10 -8 , (c) kb -π = 10 -8 , (d) kb -π = 0.1.

Fig. 6 .

 6 Fig. 6. Relative amplitude with respect to the incoming wave amplitude for ( ) kb -π =-0.1, ( ) kb -π =-10 -8 , ( ) kb -π = 10 -8 , (-) kb -π = 0.1 along the section, (a) y = 0; (b) y = 0.99.

Fig. 7 .

 7 Fig. 7. Sketch of a half-cylinder and its discretization.

Fig. 8 .

 8 Fig. 8. Relative surface wave amplitude field (a) n steps = 5, n modes = 5; (b) n steps = 5, n modes = 10; (c) n steps = 10, n modes = 5; (d) n steps = 10, n modes = 10.

  29 and C M = 1, and the dispersion relation (27) with C M = 1. The first five modes (P = 4) are considered for the present calculation, n step = 5 is used for the discretization of the cylinders.

Fig. 9 .

 9 Fig. 9. Relative wave amplitude along a x-transec with respect to the incoming wave amplitude, influence of the number of modes and of the number of steps of discretization of the cylinder, ( ) n steps = 5, n modes = 5; ( ) n steps = 10, n modes = 5; ( ) n steps = 5, n modes = 10; (-) n steps = 10, n modes = 10; (a) y = 0, (b) y = 0.6; vertical lines indicate the cylinder position.

Fig. 10 .

 10 Fig. 10. Relative error (a) of the amplitudes, (b) of the reflection coefficients versus the number of steps n steps . (a) (-*) x =-0.0135, y = 0; ( ) x = 0.5085, y = 0; ( ) x =-0.0135, y = 0.6; ( ) x = 0.5085, y = 0.6. (b) (-*) first mode; ( ) second mode; ( ) third mode.

Fig. 12 .

 12 Fig. 12. Reflection coefficient in the case of sparse array of 11 cylinders (a/ L = 0.10, L/W = 1) versus kL/π: (--) first mode, P = 2 (3 modes); (-) first mode, P = 3 (4 modes); ( ) second mode, P = 2 (3 modes); ( ) second mode, P = 3 (4 modes); ( ) third mode, P = 2 (3 modes); ( ) third mode, P = 3 (4 modes).

Fig. 13 .

 13 Fig. 13. (a) Reflection versus k x1 L/π: ( ) reflection coefficient for the second mode; ( ) total reflected energy; (-) reflected energy for the first mode; ( ) reflected energy for the second mode; ( ) reflected energy for the third mode. (b) Reflection versus k x2 L/π: ( ) reflection coefficient for the third mode; ( ) total reflected energy; (-) reflected energy for the first mode; ( ) reflected energy for the second mode; ( ) reflected energy for the third mode.

Fig. 11 .

 11 Fig. 11. Sketch of the aligned circular cylinders.

Fig. 14 .

 14 Fig. 14. Sketch of the staggered circular cylinders.

Fig. 16 .

 16 Fig. 16. Reflection coefficient versus dimensionless frequency, ( ) Molin et al. [12]; ( ) Yu and Chwang [25], C M = 1, without dissipation; (-) Arnaud et al. [7], C M = 0.29, without dissipation; ( ) present method, P = 4 (5 modes).

Fig. 17 .

 17 Fig. 17. Reflection coefficient versus dimensionless frequency, (-) staggered, one mode (P = 0); ( ) aligned, one mode (P = 0); ( ) staggered, three modes (P = 4); ( ) aligned, three modes (P = 4).

Fig. 15 .

 15 Fig. 15. (a) Reflection coefficient versus frequency, (-) Arnaud et al. [7], with linear dissipation ( = n 25 w ); ( ) Arnaud et al. [7], without dissipation ( ⟶∞ n w ); ( ) present method, P = 3 (4 modes); ( ) present method, P = 4 (5 modes); (*) experimental data from Arnaud et al. [7]. (b) Relative error versus the number of modes: (*) f = 0.5Hz; ( ) f = 1Hz; ( ) f = 1.5Hz; ( ) f = 2Hz.

Acknowledgements

Financial support from the University of Toulon for the visit of K. Belibassakis as invited Professor is acknowledged.