Gwendoline Arnaud
email: gwendoline.arnaud@univ-tln.fr

Vincent Rey
email: vincent.rey@univ-tln.fr

Julien Touboul
email: julien.touboul@univ-tln.fr

Damien Sous
email: damien.sous@mio.osupytheas.fr

Wave propagation through dense vertical cylinder arrays: 3D experimental study

Keywords: Wave, Porous medium, Dissipation, Diffraction, Reflection, Integral matching numerical method 1 Wave propagation

The purpose of this research work is to study the diffraction of surface gravity waves propagating through rectangular porous medium in three dimensions. The considered porous structure consists in dense arrays of surface piercing vertical cylinders. Experiments for different regular wave conditions have been carried out, especially for three wave frequencies. The experimental data of wave refraction-diffraction and reflection have been compared to computed results from potential linear theory solved with an integral matching method. Comparison with a previous 2D study about wave propagation through porous medium in a 10 m long wave flume is also discussed in order to highlight the refraction-diffraction effect due to the porous structure.

Introduction

Coastal countries are generally facing the over-urbanization of their coastlines and have to consider the increasing risks of damages in the nearshore due to climate ered. Arrays of vertically emerging bodies have also 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 threedimensional bodies in water waves: an exact algebraic method[END_REF]; [START_REF] Ohl | Water wave diffraction by a cylinder array. Part 1. Regular waves[END_REF] and references cited herein). However, these studies, which are often based on wave field expansion in a circular coordinate system, assume sparse arrays of cylinders. Application to either sparse or dense emerging vertical cylinder arrays was recently proposed by [START_REF] Rey | Water wave scattering by dense or sparse arrays of surface-piercing bodies by integral matching method[END_REF] for uneven periodic cylinders normally to the incoming wave direction. For dense arrays of vertical emerging cylinders, both cylinder diameter and spacing are small compared to the wavelength. With this approximation, the cylinder array can be considered as a porous medium, as done by [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF] in a one-dimensional case.

The purpose of the present work is to improve the understanding of wave propagation through a porous structure involving a cylinder array homogeneous at the wave scale, in terms of coastal protection. To do so, the general expressions of velocity potentials are expanded in a cartesian frame of reference for each finite or semi-finite medium composing the wave propagation media. The present analytical model relies on the complex dispersion relation presented by [START_REF] Yu | Wave motion through porous structures[END_REF] for the wavelength determination within the porous medium. Since the porous media is considered as homogeneous at the wave scale a linear damping approach is used to take into account wave dissipation inside the porous media. Wave damping parametrization takes advantage of the careful experiments including various wave steepnesses and periods carried out in the 2D case by [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF] for similar structure.. The problem solution is obtained by integral matching method resolving continuity equations of velocity and pressure at vertical boundaries between domains. The model is confronted to a corresponding physical model, and the performance of the system is interpreted by intercomparing this model with experimental results.

Experimental setup

Wave basin

Experiments have been carried out in a basin 10 m long and 2.60 width, of the Engineering School of the University of Toulon (SeaTech). At one end, a "corner" type wave-maker is used to generate regular waves in the range 0.7 -2.5 Hz. At the other end, a parabolic porous beach is installed to maximize incoming waves dissipation and avoid spurious reflection. The typical water depth is one meter but in the present study an additional bottom was used in order to fulfill finite water depth conditions. The bottom is then made of a sloping ramp from the wavemaker and the horizontal part supporting the porous medium and imposing a constant water depth h = 0.23 m.

Porous medium

The model porous medium is made of a vertical cylinder array. Cylinders are uniformly staggered along two perpendicular axes. These axes form a 45 the longitudinal axis of the structure. A detailed study of such structure is presented in [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF] in the 2D case. In the present study, the cylinder diameter is D = 0.032 m. The length of the structure is L = 0.30 m and the width is 2d p = 1.20 m. The distances between cylinders along horizontal axis are dx = 0.0384 m and dy = 0.0365 m, respectively (see Fig. 2). The porosity γ and the specific surface s are defined in Eqs. 1 and 2, respectively.

γ = 2d p L -N π(D/2) 2 2d p L (1) s = πDN 2d p L (2)
where N is the number of cylinders of porous medium, here N = 120.

In the present study γ = 0.7 and s = 33 m -1 .

Sensors and experimental conditions

The free surface has been measured with thirteen acoustic wave probes deployed in the near vicinity of the porous structure (see Fig. 1). The sampling frequency is 40 Hz. Measurements of the wave amplitude are carried out around the structure according to a grid of 0.10 × 0.10 m 2 by moving the racks of gauges. For this experiment, wave frequency and wave height are chosen in intermediate to almost deep water conditions. Frequencies and amplitudes corresponding to the dimensionless wavelength inside and around the porous medium are given in Table 1. Note that the conditions studied concern non breaking waves. λ and λ p are the

Frequency

Amplitudes (m) h/λ h/λ p f = 1 Hz 0.010 0.015 0.020 0.18 0.19 f = 1.5 Hz 0.010 0.015 0.020 0.34 0.38 f = 1.8 Hz 0.010 0.015 0.48 0.53 wave wavelengths outside and inside the porous structure, respectively.

Theoretical Models

The case of an incident wave in the Ox direction is considered. The wave propagates towards positive x direction in the longitudinal axis of the basin. Oy is orthogonal to the incident wave direction with y = 0 corresponding to the median axis of the basin. Oz is vertical upwards with z = 0 corresponding to the position of the still water level. The coordinate system (O, x, y, z) is a direct orthonormal coordinate system. Using symmetry properties with respect to y = 0, the problem is solved for y ≥ 0. The problem solution is found writing a general expression of the potential in four domains:

• Domain 1: upstream, 0 < y < d B , x < 0, • Domain 2: porous medium, 0 < y < d p , 0 < x < L p , • Domain 3: sidelong, d p < y < d B , 0 < x < L p , • Domain 4: downstream, 0 < y < d B , x > L p ,
and solving the continuity equations between adjacent domains (Fig. 3).

General equations of velocity potentials

For each rectangular domain i of width d i , d m < y < d M , and of length L i , the general expression of the velocity potentials satisfying the impermeability conditions at y = d im and y = d iM are of the form:

Φ i (x, y, z, t) = cosh[k i (z + h)]φ i (x, y)e iωt = cosh[k i (z + h)] ∞ n=0 A - i,n e -ik ixn x + A + i,n e +ik ixn x ψ i,n (y)e iωt (3) i = 1, ..., 4, with ω 2 S r = gk i tanh k i h (4)
with k ixn = k i cos θ n and k iyn = k i sin θ n , n is a priori infinite (but practically truncated to an order n = P ). For domains 1, 3 and 4, S r = 1 in contrast to domain 2 where S r = 1 + C m 1-γ γ with C m = 0.29 (see [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF]). The wave number is k i in domains i = 1, ...4. In order to match the impermeability condition which assumes a zero normal velocity, at y = d im ,

ψ i,n (y) = cos [k iyn (y -d im )]
(5) For an impermeable boundary at y = d iM , k iyn must satisfy:

k iyn = nπ d iM -d im (6)
The surface deformation η i (x, y, t) is given by:

∂η i ∂t = ∂Φ i ∂z z=0 (7)
Wave amplitude is related to the potential amplitude through the following relation:

A ± j = iω k j sinh k j h a ± j (8) Since ∇ 2 Φ i = 0, k ixn = k 2 i -k 2 iyn 1 2 (9)
According to the k xn expression (eq. 9), from a given number of order n p > 0, k xn becomes pure imaginary. The direction of propagation θ n for each mode n, n ≤ n p , with respect to the Ox direction of incident wave is given by:

θ in = arctan k iyn k ixn
For n > n p , the modes do not correspond anymore to propagative waves but to evanescent waves, and Eq. ( 9) can be written in the following form:

k ixn = i k 2 i -k 2 iyn 1 2 (11)
For an impermeable boundary y = d im and a permeable boundary y = d iM , the function ψ n can still be written:

ψ n (y) = cos [k iyn (y -d im )] (12)
For a permeable boundary in y = d im and an impermeable boundary y = d iM , the function is written:

ψ n (y) = cos [k iyn (y -d iM )] ( 13
)
Expressions for domain 1 (x < 0, 0 < y < d B ):

The incident wave is known. Writing that only the incident wave (θ = 0) propagates toward x > 0, and that evanescent modes do not diverge at -∞, the general expression of the reduced potential φ 1 (x, y) is the following:

φ 1 (x, y) = A - 10 e -ik 1x0 x ψ 10 (y) + ∞ n=0 A + 1n e +ik 1xn x ψ 1n (y) (14)
with

ψ 1n (y) = cos [k 1yn (y)] and k 1yn = nπ d B ( 15
)
Expressions for domain 3 (0

< x < L p , d p < y < d B ):
The general expression of reduced potential φ 3 (x, y) is given by:

φ 3 (x, y) = ∞ n=1 A - 3n e -ik 3xn x + A + 3n e +ik 3xn (x-Lp) ψ 3n (y) (16)
with

ψ 3n (y) = cos [k 3yn (y -d B )] ( 17
)
k 3yn is chosen of the form k 3yn = nπ d B , which corresponds to the decomposition mode in the absence of porous medium.

Expressions for domain 4 (x > L p , 0 < y < d B ): Since there is no wave reflection from +∞ whatever θ angle, and evanescent modes do not diverge at +∞, the general expression of the reduced potential φ 4 (x, y) is of the form:

φ 4 (x, y) = ∞ n=0 A - 4n e -ik 4xn (x-Lp) ψ 4n (y) (18)
with

ψ 4n (y) = cos [k 4yn (y)] and k y4n = nπ d B ( 19
)
Expressions for domain 2 (0 < x < L p , 0 < y < d B ):

The general expression of the reduced potential φ 2 (x, y) writes:

φ 2 (x, y) = ∞ n=0 A - 2n e -ik 2xn x + A + 2n e +ik 2xn (x-L) ψ 2n (y) (20)
with

ψ 2n (y) = cos [k 2yn (y)] (21)
As for domain 3, k 2yn is chosen of the form k y2n = nπ d B . In [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF] both linear and quadratic expression for wave damping were discussed. We consider here the damping is linear because the porous structure length L p is considered small compared to the wave's wavelength. A constant attenuation rate of wave amplitude is assumed in the direction of propagation within the porous medium dissipation. It is neglected elsewhere (domains 1, 3 and 4). For a partially stationary wave of the form η 2 = a - 2 e i(ωt-k 2 x) + a + 2 e i(ωt+k 2 x) , the expressions of a - 2 and a + 2 for a porous medium of length L p such as 0 x L p write:

a - 2 (x) = a - 2 (0)e -k 2,d x a + 2 (x) = a + 2 (L p )e -k 2,d (Lp-x) (22)
where k 2,d = k 2 /n w is the damping rate per unit length, n w is chosen constant. It is obtained empirically from comparison to experimental data (see [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF]).

Numerical method of resolution

We present here not only the case of porous structure but also the case of a porous structure with impervious lateral boundaries (at y = ±d p ) since both are computed for comparisons with experimental data.

Porous structure

The continuity conditions for pressure p = -ρ ∂Φ ∂t and velocity at vertical boundaries are applied to the vertical extremities of the porous medium: For x = 0:

Φ 1 = S r Φ 2 and ∂Φ 1 ∂x = γ ∂Φ 2 ∂x for 0 y d p Φ 1 = Φ 3 and ∂Φ 1 ∂x = ∂Φ 3 ∂x for d p y d B (23) For x = L p : S r Φ 2 = Φ 4 and γ ∂Φ 2 ∂x = ∂Φ 4 ∂x for 0 y d p Φ 3 = Φ 4 and ∂Φ 3 ∂x = ∂Φ 4 ∂x for d p y d B (24) For y = d p S r Φ 2 = Φ 3 for 0 x L p γ ∂Φ 2 ∂y = ∂Φ 3 ∂y for 0 x L p (25)
Continuity equations between domains:

The weak formulation of momentum and pressure continuity conditions is used by double integration over both y and z directions, for n = 0, ..., P based on eigenfunction expansion matching method. Integration along the z-axis on the water depth gives the following constants K and K p .

K = 0 -h cosh[k(z + h)]dz and K p = 0 -h cosh[k p (z + h)]dz
where k = k i for i = 1, 3, 4 and k p = k 2 . At boundaries x = 0 and x = L p , classical integral matching method is used, with the choice of orthogonal eigenfunctions for the integral formulation. For the boundary conditions at x = 0 and x = L p , basic functions ψ j,n , n = 0, ..., P are used (j = 1, 4 for x = 0 and x = L p respectively). For the boundary conditions at y = d p , we consider the functions Ch n (x) = cosh[l n x] with l n = nπ Lp , which form an orthogonal set since,

Lp 0 Ch n .Ch m dx = 0 if n = m. ( 27
)
The continuity equation can then be written, for x = 0: After a development of the potential expressions ( 14), ( 16), ( 18) and (20), a set of 6(P + 1) linear equations, with 6(P + 1) unknowns A ± are then to be solved.

K d B 0 φ 1 ψ 1n (y)dy = S r K p dp 0 φ 2 ψ 1n (y)dy + K d B dp φ 3 ψ 1n (y)dy K d B 0 ∂φ 1 ∂x ψ 1n
Without dissipation by the porous medium (n w → ∞), the wave energy flux along the Ox direction, over the basin width d B , is conserved:

|A - 10 | 2 = |A + 10 | 2 + |A - 40 | 2 + 1 2 np n=1 |A + 1n | 2 + |A - 4n | 2 (31)
3.2.2. Porous structure with impervious walls at y = d p

In order to fulfill the zero normal velocity at y = d p ,

k 2yn = nπ d B -d p and k 2yn = nπ d p ( 32
)
The number of unknown coefficients is 6(P + 1). For x = 0:

     K dp 0 φ 1 ψ 1n (y)dy = S r K p dp 0 φ 2 ψ 1n (y)dy

Results and discussions

Wave amplitude fields around the porous medium have been measured for the wave conditions given in Table 1. Our first aim was to analyse the wave scattering by the porous structure including reflection, refraction-diffraction and dissipation phenomena. Let us note that the side walls of the basin may affect the wave behaviour near the structure, since d B is not large compared to d p . Such a configuration can also describe the case of a periodic detached breakwater system [START_REF] Rey | Water wave scattering by dense or sparse arrays of surface-piercing bodies by integral matching method[END_REF]. Intermediate water depth conditions are considered for f = 1 and 1.5 Hz, while f = 1.8 Hz nearly corresponds to deep water conditions since h/λ = 0.48 and h/λ p = 0.53 where λ p = 2π/k p . For the three wave conditions, the wave's wavelength remains much higher than the cylinders size and spacing. The length of the porous structure remains lower than the wavelength for the three frequency cases, while the width 2d p = 1.20 m is of the order of the wavelength for f = 1 Hz, twice the wavelength for f = 1.5 Hz and five times for f = 1.8 Hz.

Maps of wave amplitude

Figures 4 to 6 depict 2D maps of relative wave amplitudes for each wave frequency, for incoming wave amplitude 0.010 m for f=1Hz and 0.015 m for f = 1.5 Hz and f = 1.8 Hz. Results concerning the other amplitudes (see table 1) are shown and discussed in sections 4.2 and 4.3 for longitudinal and transversal axes, respectively. Relative wave amplitude corresponds to the ratio between the wave amplitude measured around the porous medium and the wave amplitude measured at the position of the reference wave gauge upwave the structure (see Fig. 1). For each case, plots of wave amplitude from the experimental results are compared to numerical results for the porous structure. The experimental wave fields are obtained after using a cubic spline interpolation technique.

For the first studied frequency f = 1 Hz and a = 0.010 m, the wave amplitude response upstream the porous structure along the central axis of the basin presents successive maxima and minima, representing a wave reflection phenomenon (Fig. 4, top). Downstream the porous structure the wave amplitude is rather homogeneous; a part of the wave being transmitted through the porous structure. Wave refractiondiffraction is observed on both edges of the structure.

Figure 4 (bottom) represents the theoretical model with a porous structure. The wave reflection is well represented regarding the 3D wave amplitude pattern.

The wave amplitude field downstream of the porous structure appears to be smooth do the transmission of a significant part of the wave energy through the structure. We can notice however that the amplitude variations are more pronounced for the experimental case.

For f = 1.8 Hz, the experimental result (Fig. 6, top) shows wave reflection as for f = 1 Hz but the extrema of wave amplitude upstream the porous structure are more numerous due to a higher wavelength. Downstream the structure the wave rays are curved particularly near to structure sides. The transmission is also lower than for the first studied frequency. The theoretical model (Fig. 6, bottom) fits correctly the experimental results. The oscillations of wave amplitude upstream the structure are well represented as the wave rotation downstream the porous medium.

For the wave frequency f = 1.5 Hz, the experimental data give another pattern compared to other cases. The wave amplitude oscillations is still observed along the wave propagation direction but another amplitude oscillation is observed along the transverse direction. Downstream of the structure, such a wave behaviour is also observed even if it is less pronounced due to wave dissipation and backscattering. Such transverse oscillations are well predicted qualitatively by both theoretical models for the porous medium (Fig. 5, bottom).

To summarize these results, the first overall observation is the presence of significant wave reflection and refraction-diffraction by the porous medium together with a partial transmission which smooths the diffraction pattern downwave the structure. Upstream the structure (y < 0), the series of successive minima and maxima can be seen as the result of the partial reflection by the porous medium even if diffraction effects arising from the structure edges may modify the wave patterns. These spatial modulations of wave amplitude upstream the porous block are observed for the three wave conditions, with an increasing number of oscillations for increasing wave frequency. Downstream of the porous medium, refraction-diffraction effects can be observed. The wave direction varies along the basin width, propagating faster on both sides of the structure and forming a cylindrical-type wave. The wave energy downstream the porous medium results from the combination of transmission through the porous medium and diffraction at its edges, producing a more complex wave structure than for an impervious structure. However, some discrepancies are observed between experimental data and theoretical results. They are certainly due to wave scattering from the beach and from the basin walls. Indeed, the basin and the beach are designed for wave propagation along the x-axis, scattered waves by the porous structure may be reflected by the beach, not designed for oblique incidence.

The theoretical results for the porous structure are found to provide a correct representation of the surface wave field obtained for the three studied frequencies. Since a given dissipation rate n w [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF] per unit wavelength is considered, higher dissipation is observed for the higher frequency case, the porous length L p being kept constant in the experiments. For f = 1.5 Hz, we can also observe max- ima and minima of amplitude along the cross section ascribed to resonance effects since the gap between the side walls and the porous medium (0.7 m) corresponds nearly to the wave wavelength (λ = 0.68 m) and the porous width (2d p = 1.20 m) to the wave wavelength inside the porous medium (λ p = 0.61 m).

Wave amplitude along a longitudinal transect

Measurements along the section y = 0.10 m, i.e. near the longitudinal axis of the basin, are compared to the 3D model results for the porous structure. The results are also compared to the 2D theoretical model neglecting structure edges lateral effects (d p → ∞) [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF].

Overall, Figure 7 shows that the theoretical results for the 3D case for a porous structure are in good agreement with experimental data both upstream and down- Upstream the structure, the amplitude of wave modulation is greater when increasing frequency. This is induced by the frequency dependency of the reflected wave amplitude due to interference processes shown in [START_REF] Arnaud | Wave propagation through dense vertical cylinder arrays : Interference process and specific surface effects on damping[END_REF]. Downstream the structure, the sheltering effect of the porous medium progressively decreases and the wave amplitude increases slowly with some weak oscillations.

Frequency effects are highlighted when comparing porous 2D and 3D model. For frequencies f = 1 Hz (Fig. 7, top) and f = 1.8 Hz (Fig. 7, bottom) both models upstream the structure are very close. This indicates that side effects are weak along the median axis. By constrast, the case f = 1.5 Hz (Fig. 7, middle) highlights the influence of the finite width of the structure with significant differences between 3D and 2D cases. This discrepancy can be explained by the presence of transverse resonance in the 3D case. This resonant mode may then be due to either or both the basin width or the gap between the basin and side walls of the porous structure.

In addition, one note also that downstream the structure for f = 1 Hz large amplitude oscillations are measured but not provided by the models. They may be due to resonance behaviour between the structure and the absorbing beach used in the basin which is not perfect at low frequencies. The theoretical result with a wave reflection of the beach of 20% highlights the presence of an amplitude modulation downstream according to the observed modulation. However the predicted oscillation remains much lower even if beach reflection is taken into account in the model (Fig. 7, top). A possible explanation for this for this discrepancy might be related to a ydependent reflection of the beach, resulting from wave focusing behind the structure.

Wave amplitude along a transversal transect

The comparison between measurements and numerical results along the transversal direction of the basin around the porous structure is presented in Fig. 8.

These plots highlight the 3D effects of the porous structure on the wave propagation. The position of the transects are chosen to corresponds to minimum and maximum upwave the cylinder array, respectively. The general tendency of the wave amplitude is well represented by the model, the minima and maxima are well positioned along the basin width although the experimental results show higher amplitude for f = 1 Hz and smaller amplitude for f = 1.8 Hz.

Discussion on the wave resonance

In order to better understand the observed transverse resonance occurring for f = 1.5 Hz, we present here results from a 3D model configuration neglecting the fluxes through the side walls of the porous medium. A comparison of measured wave amplitude and both theoretical cases (porous and impervious lateral boundaries) for the three studied frequencies is depicted in Fig. 9.

For frequency f = 1 Hz (Fig. 9, top),the wave amplitude modulations are similar for both theoretical models. This indicates that reflection and transmission are not influenced by the exchanges through the side walls of the porous structure. For f = 1.5 Hz (Fig. 9, middle), significant discrepancies are observed between model configurations with strong differences of wave modulations along the incoming wave direction, upstream and downstream the porous structure. As already mentioned in the previous section, the model which considers porous side walls suits better with experimental data, particularly upstream of the structure with a good representation of the extrema. This indicates that fluxes and/or pressure drops across the porous side walls have a significant influence at this frequency due to transverse resonance, while they are much less significant for f = 1 Hz and f = 1.8 Hz.

For the higher frequency, f = 1.8 Hz (Fig. 9, bottom) both models give same results upstream of porous medium and are in good agreement with experimental data. This means that upstream the center part of porous structure there is no 3D effect as for f = 1 Hz. Downwave the structure, the modulation is found to be higher when the structure sides are considered as porous walls. Both theoretical trends fit with experimental data even if oscillations observed experimentally may be due to beach scattering as already mentioned in the previous subsection (subsec. 4.2). The diffraction effects are then significant downstream also in its central part because the wave wavelength is of the same order as the porous length.

In the absence of the transverse resonance, wave behaviour upstream the struc- ture far from its edges can be well described by the 2D model. By contrast, in the presence of transverse resonance, the whole wave field is affected by 3D effects.

Conclusions

This combined experimental and theoretical study was dedicated to the analysis of the interactions between regular waves and porous structure of finite dimensions.

Measurements of the wave field in the presence of the porous structure allowed to evidence the effect of a porous structure on the wave scattering. Further insight was provided by the theoretical approaches, which compare 2D and 3D models. The 3D model provides a satisfactory representation of the wave field around the structure, confirming the validity of the approach. The discrepancies observed between model and measurements downstream the structure especially at the lower frequency are attributed to wave scattering by the absorbing beach which hardly dissipates wave energy at low frequency in the experiments.

Transverse resonances related to both structure width and the basin side walls are shown by the measurements and accurately modelled. Two configurations, with porous and impervious lateral boundaries, are used in the model to better understand such processes. The resonance effect by the porous structure is amplified when side-walls are considered as impervious. The difference between both hypothesis in theoretical models is significant only for the case f = 1.5 Hz due to the presence of transverse resonance, whereas the difference is rather negligible for both other studied wave conditions.

On the whole, model gives a good prediction of wave transformation processes through and around porous structure of finite dimensions. Further efforts will be engaged to extend the model to more complex configurations of porous medium, in particular including inhomogeneous porosity or specific surface within the structure.

4

  Arnaud G., Rey V., Touboul J. and Sous D.

Figure 1 .

 1 Figure 1. Experimental set-up

Figure 3 .

 3 Figure 3. Scheme of half-width of tank with different domains considered in the theoretical model (left); Scheme of domain definition in the model (right)

  for n = 0, ..., P . Similarly for x = L p :S r K p dp 0 φ 2 ψ 4n (y)dy + K d B dp φ 3 ψ 4n (y)dy = K d B 0 φ 4 ψ 4n (y)dy γK p dp 0 ∂φ 2 ∂x ψ 4n (y)dy + K d B dp ∂φ 3 ∂x ψ 4n (y)dy = Kd B 0 ∂φ 4 ∂x ψ 4n (y)dy (29) n = 0, ..., P and, for y = d p : S r K p Lp 0 φ 2 Ch n (x)dx = K Lp 0 φ 3 Ch n (x

  ψ 1n (y)dy = γK p dp 0 ∂φ 2 ∂x ψ 1n (y)dy + K d B dp ∂φ 3 ∂x ψ 1n (y)dy Wave propagation through dense vertical cylinder arrays: 3D experimental study 11 For x = L p : condition (31) is still verified.

  Figure 4. Wave amplitude for f = 1 Hz experiments (top), theory (down)

  14 Arnaud G., Rey V., Touboul J. and Sous D.

Figure 5 .

 5 Figure 5. Wave amplitude for f = 1.5 Hz experiments (top), theory (down)

  Figure 6. Wave amplitude for f = 1.8 Hz experiments (top), theory (down)

March

16

  Arnaud G., Rey V., Touboul J. and Sous D.

Figure 8 .

 8 Figure 8. Transversal cuts of the basin for f = 1 Hz (top), f = 1.5 Hz (middle) and f = 1.8 Hz (down)

  Wave propagation through dense vertical cylinder arrays: 3D experimental study 5

	March	18,	2019 14:39 WSPC/INSTRUCTION	FILE
	Arnaud-et-al˙IJOCE˙Review

• angle with

Table 1 .

 1 Wave conditions

  Wave propagation through dense vertical cylinder arrays: 3D experimental study 9

	March	18,	2019 14:39 WSPC/INSTRUCTION	FILE
	Arnaud-et-al˙IJOCE˙Review

Arnaud G., Rey V., Touboul J. and Sous D.

Arnaud G., Rey V., Touboul J. and Sous D.

Acknowledgements

The authors would like to thank the Region Provence Alpes Côte d'azur and the company ACRI-IN for their financial contributions of the PhD thesis scholarship of Gwendoline Arnaud that allowed this work. The authors are grateful Aimed Ajroud for his help to the physical models conception.


L'archive ouverte pluridisciplinaire

Il. Access Point association in WLAN-The two modes of operation defined in the 802.11 standard are described: the Infrastructure and the Ad-Hoc modes, where there is no access points. In this paper, we consider only the infrastructure mode. The procedures that a ST A should perfonn in order to associate with an AP such as scanning, authentication ( [START_REF] Walker | Unsafe at any key size: An analysis of the WEP encapsulation[END_REF]) and association itself are reminded. Two scanning procedures are defined: active and pa.<:sive scanning. A STA, after scaiming and mak:ing a list of al! available APs in the network, will choose an AP to associate with. Which critcrion must the STA use to select this AP? This criterion is not specifically slated in 802.11 standard. DifTerent possible method of selection are applied in manufacturcrs' dcviccs. In most prcscntly-proposcd 802.11 products, the STA chooses to associa.te with the AP that has the best pe1formanccs in the tcrm of Signal to Noise Ratio (SNR) at the ST A geographical position.

We highlight the inconvenients of the cmrently-used methods of selection of the access points by the wireless stations and try to find and to propose solutions that can rectify this situation. III. Model and results-We simula.te a simple WLAN model for differcnt critcrias of AP selection: received power lev el, receivcd si gn al-to-intereference ratio and AP-ST A distance. Channels are a'>sîgned to APs according to 802.11 standard. We remînd in the paper APs channel assignment mies.

To calculate the Signal to Noise Ratio (SNR) betwecn each station and cach access point, we consider that Noise = sum of the transmitted power of each AP using the same channel * link gain between this AP and the wireless station.

In the paper, simulation results are commented. For example, Figure 2 shows the AP cclls (callcd Basic Service Arca or BSA) for a non-uniform distribution of wireless stations, where the best SNR criterion is used. W e see that this type of AP selection leads to a cell that bas a high concentration of stations. In fact, cells shapes are indendent of the concentration of wireless stations. A WLAN cell with a high number of stations might have quickly bad perfonnances: low data rate and high del ay s. In the paper we propose some solutions such as transmitted power decrease for the AP that bas a high data load.

IV. Conclusion

In this paper, we highlight the fact that if the load distribution isn't con;;idered for AP selection in a WLAN, this may lead to inefficient use of the bandwidth.

Sorne solutions that can rectify this situation are proposed and discussed at the end of the article: power variations and proposing and implementing load distribution algorithms in the 802.11 standard are enhancements that can lead to a gn:ater capacity and better perfo1mances. Research papers proposing solutions are mentioned (e.g., [START_REF] Wang | Integration of SNR, Load and Time in Handoff. Initiation for Wireless LAN[END_REF])

Access

  Point Association in IEEE 802.11 WLAN Loutfi Nuaymi -Jinane El-Sayah ENST Bretagne 2, Rue de la Chataigneraie, B.P. 78, Cesson Sévigné 35512 Cedex, FRANCE e-mail: loutfi.nuaymi(Wenst-bretagne .fr ; j inanesayah@hotmail.com Keywords-WLAN, IEEE 802.11, selection, association, access points Subject-Here is an extended abstract of the work we would like to present at ICTT A'04 (April 04) I. Introduction-ln fEEE 802. l l [ 1] Wireless Local Area Networks (WLAN), each wireless station (STA) is associated to one and only one AP (Access Point). The selection of the AP for every wireless stations is an important procedure. A proper choice incrcases the total capacity of 802.11 WLAN and then leads to efficient use of the available frequency bandwidth. A bad cell selection can cause wireless networks blockings that could be avoided ([3]) as we show it in this paper.

Fig. 2 -

 2 Fig. 2-For a non-uniform distribution of wireless stations, shape ofWLAN cells for an AP selection based on the best received SNR. The central cell bas a high number of stations. This leads to bad performances.


Alexandra Lacombe

Yann Landon

Manuel Paredes

Clément Chirol

Audrey Benaben

Numerical evaluation of residual stress induced by reaming of aluminium 2024-T351

Keywords: Reaming, numerical simulation, residual stress, aluminum

  L'archive ouverte pluridisciplinaire

Introduction

Fastening holes in mechanical structures are critical areas for the initiation of fatigue damage, especially for metal parts. The experimental study of fatigue strength and failure modes of bolted joints made of aluminum alloy parts is the subject of a lot of experimental works [START_REF] Krishnakumar | Fatigue behaviour of aluminium alloy 7075 bolted joints treated with oily film corrosion compounds[END_REF][START_REF] Croccolo | Failure analysis of bolted joints: Effect of friction coefficients in torque-preloading relationship[END_REF][START_REF] Hoang | Rupture and damage mechanism analysis of a bolted assembly using coupling techniques between A.E. and D.I[END_REF][START_REF] Chakherlou | An experimental investigation of the bolt clamping force and friction effect on the fatigue behavior of aluminum alloy 2024-T3 double shear lap joint[END_REF]. Different failure modes are observed. Some failure occurs in the bearing mode, with cracks initiating at the edge of the hole and others occurs in the fretting mode, with cracks initiating away from the hole boundary.

Compressive residual stress at the vicinity of a hole resulted in a significant increase in the fatigue lives of holed parts [START_REF] Lai | Fatigue properties of holes with residual stresses[END_REF]. Indeed, cracks find it much more difficult to propagate through a compressive field. Thus, in the aeronautical industry, processes have largely been developed in order to introduce a compressive residual stress field at the vicinity of the holes. Interference fit fastener [START_REF] Chakherlou | An investigation about interference fit effect on improving fatigue life of a holed single plate in joints[END_REF][START_REF] Lanciotti | The effect of interference-fit fasteners on the fatigue life of central hole specimens[END_REF] or cold expansion process [START_REF] Ozdemir | Effect of expansion technique and plate thickness on near-hole residual stresses and fatigue life of cold expanded holes[END_REF][START_REF] Lacarac | Fatigue crack growth from plain and cold expanded holes in aluminium alloys[END_REF][START_REF] Chakherlou | The effect of cold expansion on improving the fatigue life of fastener holes[END_REF], which consist in introducing a conical mandrel with a high level of interference, allow to improve greatly fatigue life of aluminum alloy joints. Compressive residual stress fields introduced by these two processes are well known to be extensive and to reach high values. However, other processes involved in the machining of the fastener hole, such as drilling and reaming, also leave a region of disturbed material along the edges of the hole that can cause a residual stress field. Even if these induced residual stress fields are much lower than those left by a process such as cold expansion (Figure 1), they affect the fatigue behavior of the holed part [START_REF]Assessment of residual stresses and hole quality on the fatigue behavior of aircraft structural joints 1[END_REF]. Therefore, they need to be better apprehended in order to be accounted for the life-prediction methodology. This paper focuses on the modification of the stress state of an aluminum 2024-T351 part during the reaming process which has been weakly studied. This alloy was selected because of its use in the majority of the current fleet of commercial aircraft. Only the numerical part of the study is presented in this paper.

Generation of residual stress during reaming is a complex phenomenon to simulate numerically. It involves both mechanical and thermal phenomena. On the one hand, simulation of metal cutting, which is associated with large strains, requires Eulerian methods. On the other hand, the prediction of residual stress, which needs knowledge of the loading history of each material particle, requires Lagrangian methods. Thus, several authors choose to develop « decoupled » numerical models, composed of a first model simulating the cutting, allowing them to obtain an equivalent thermomechanical loading, which they apply in a second model directly on the finished part. This is for example the case of Girinon [START_REF] Girinon | Développement d'un modèle numérique sur l'étude de l'intégrité des surfaces en perçage[END_REF], who, in his work on numerical simulation of drilling in hard metals, develops a rigid ALE model of material removal, then extracts the thermal loading on the part which is then applied in a second model pure-ly Lagrangian. Other authors, such as Valiorgue [START_REF] Valiorgue | A new approach for the modelling of residual stresses induced by turning of 316L[END_REF] in his work on the simulation of the turning of hard metals, free themselves from the modeling of material cutting by experimentally identifying the thermomechanical equivalent loading to apply in their Lagrangian model. Still others adopt even simpler reaming modeling strategies. Achard [START_REF] Achard | Numerical modelling of the cold expansion process in mechanical stacked assemblies[END_REF], to simulate the reaming process after cold expansion on hard metals, chooses to simply deactivate the elements located in the reaming zone of the part in its finite element model and define a calculation step allowing the part to reach a new static equilibrium.

Most authors studying the numerical simulation of machining processes consider a material devoid of mechanical and thermal stress. However, in the aeronautical industry, parts do not present a zero stress state before reaming since they undergo at least a drilling operation, and sometimes a cold expansion operation. Thus, it is interesting to carry out a study on the impact of the reaming process taking into account the history of the loads applied to the piece. For this, a numerical model considering a prestressed part is developed. The initial stress state of the part is introduced through the simulation of the cold expansion process.

Cold expansion simulation

The Finite Elements Model (FEM) representing the cold expansion process developed with Abaqus by Achard [START_REF] Achard | Numerical modelling of the cold expansion process in mechanical stacked assemblies[END_REF] is taken over and adapted to the expansion in an aluminum alloy part in order to generate a stress state in the part before reaming. The model, presented in Figure 2, represents the workpiece and the different parts of the tooling : the conical mandrel, the sleeve and the jaw. All the parts respect the dimensions imposed by the industrial specifications to obtain a final hole diameter of 6.35 mm and the thickness of the part is 6 mm. The simulated expansion rate is 6% in order to generate important residual stress fields in the part to better observe the influence of the reaming process in the rest of the study. The model is 2D axisymmetric, the part geometries, the loadings and the boundary conditions being compatible with this behavior hypothesis.

The steel mandrel can be considered and modeled as a rigid body in relation to the expanded aluminium part. The elastic behavior of the other three parts is considered as isotropic. The Young's modulus (Table 1) are calculated on the basis of results from quasi-static characterizations and data used by Achard [START_REF] Achard | Numerical modelling of the cold expansion process in mechanical stacked assemblies[END_REF] in his work. The steel jaw which undergoes only weak plastic strain is modelled with a simple elastic constitutive law. The mechanical behavior of the steel sleeve and the aluminum workpiece which undergo severe plastic strain is described by elastoplastic constitutive laws. The plastic behavior of these parts is defined on the basis of results from uniaxial tensile strength tests. The material behavior of the workpiece is modelled by a kinematic hardening and the material behavior of the sleeve is modelled by an isotropic hardening. In order to simulate the kinematics of the expansion, a displacement of -40 mm along the model axis of revolution (z) is imposed on the mandrel and the lower face of the jaw is fixed. Four contact areas are considered with a "surface-to-surface" discretization method (mandrel/sleeve, sleeve/part, mandrel/jaw and jaw/sleeve).

The mesh, presented in Figure 3, is refined and composed of biquadratic elements with eight nodes (CAX8) in the areas of high plasticity (sleeve and hole edge area of the part). These elements which imply stress and strain calculations at 8 integration points per element allow providing accurate results that take the plastic phenomena into consideration well. The rest of the workpiece and the jaw are meshed with bilinear elements with four nodes using a reduced integration (CAX4R). These elements which use a reduced integration imply a stress and strain calculation at only one integration point by element. This allows saving computational time without losing accuracy, the areas being slightly plastically deformed. The calculation step corresponding to the simulation of the cold expansion process with the axial displacement of the mandrel is followed by a step simulating the withdrawal of the sleeve which is called "recovery". During this step, the contacts between the elements are deactivated and no loading is applied. This allows the part to achieve a new stress state.

The whole calculation is resolved according to an implicit time integration method.

Reaming simulation

The residual stress state obtained from the expanded part after the recovery step constitutes the initial stress state for the simulation of the reaming process.

Reaming is simulated by deactivating an elements band along the height of the part at the hole edge with a width of 0.35 mm. Two strategies for simulating the reaming process are compared. The first one is similar to the one implemented in Achard's work [START_REF] Achard | Numerical modelling of the cold expansion process in mechanical stacked assemblies[END_REF]: all the elements of the reaming area are deactivated in the same calculation step. This modeling strategy allows to obtain the final residual stress state after reaming, it will be called "global reaming" in this study. The second strategy consists in deactivating the elements by several steps along the height of the part in order to simulate the progressive feed of the tool during reaming. 19 calculation steps, each corresponding to the deactivation of around ten elements in the thickness are carried out (Figure 4). This allows obtaining the history of the stress state in the part during reaming. This will be called "progressive reaming".

Thus, the two modelling strategies do not take into account the thermal effects occurring during reaming. The aim is to study only the impact of material removal on the stress state of the part which cannot be studied separately experimentally. Both modelling strategies should provide similar results. In the case of progressive reaming, the machining kinetics is exploited to define the calculation steps (height of deactivated elements and step duration). This time effect will fully make sense when another model involving progressive thermal loading will be implemented.

Results and discussion

During cold expansion simulation, residual stresses are generated in the part along the three directions: axial, radial and circumferential. Axial residual stresses are largely lower than radial and circumferential stresses, so this study focuses on the latter two.

Figure 5 shows the radial and circumferential stress states in the part after cold expansion and global and progressive reaming simulations. These results show a non-homogeneous distribution of the residual stresses in the part after the cold expansion and reaming processes. The method for deactivating the elements in order to simulate the reaming process does not seem to affect the final stress state of the part.

The low influence of the reaming simulation on the residual stress state of the part can be noted. The radial residual stress field undergoes a progressive offset in the radial direction during the reaming simulation. The two compression points on the entry and exit faces are slightly relaxed and the radial residual stress tends to homogenize in the thickness of the part. This relaxation could be related to the reduction of the axial strain at the hole edge on the entry and exit faces which is generated during the cold expansion process and which is sometimes called "volcano effect". Concerning the circumferential residual stress field, this one is simply truncated in the reaming area. Thus, a higher compressive stress state favoring a longer fatigue life is obtained at the hole edge.

The radial and circumferential residual stresses on the entry face and middle section of the part after cold expansion, global and progressive reaming simulations are shown in Figure 6. Radial stresses are compressive and have a characteristic profile. They cancel at the hole edge and at the free outer edge of the part. Circumferential stresses are of mixed nature with a compressive peak at the hole edge and then a tensile zone far from the hole. A relaxation of the radial residual stress is observed during reaming simulation with a 7 % reduction in the compressive peak on the entry face and 15 % on the middle section. The modification of the circumferential residual stress is significant only on the entry face with a 6 % increase and a small offset of the compressive peak. Even though the level of residual stress changes during the reaming simulation, the heterogeneity of the stress state in the part remains. Indeed, the residual stress curves obtained after the cold expansion and reaming simulations present similar shapes. The changes in residual stress states are located in the hole edge area.

The results obtained with the two simulation strategies for reaming are identical. This confirms the equivalence of the two methods to obtain the final residual stress state of the part. The progressive reaming method seems nevertheless more interesting because it allows accessing to the intermediate stress states during the process. These useful results could be considered to gradually apply an equivalent thermo-mechanical load in order to envisage a more realistic simulation of the reaming process.

The evolution of the residual stress fields during progressive reaming is studied in order to better understand the impact of the process on the material characteristics of the part. For this purpose, radial and circumferential stresses are measured on 9 elements distributed according to the height and depth of the part. The location of these elements is described in Figure 7. The residual stress evolution for each element is presented in Figure 8. On these graphs, a dotted line indicates the position of the elements where the stress results are measured in relation to the part depth.

The evolutions of the residual stress during the reaming process present different shapes depending on the elements considered. A stress balancing takes place gradually when the elements bands are deactivated for the elements located at the same level in the part thickness. This balancing is mainly carried out in the radial direction and intensifies by approaching the hole edge. Indeed, the elements located exactly at the hole edge (E1, E4 and E7) show a sudden stress evolution. The elements a little farther away from the hole edge (E2, E5 and E8) present a smoother evolution. Finaly, the elements very far from the hole edge (E3, E6 and E9) show almost linear evolutions. The comparison of the numerical results obtained with experimental results is not yet possible. Indeed, it is impossible to implement experimentally the mechanism of material removal without involving thermal phenomena. Therefore, the results obtained are compared with previous numerical work as a first validation step. Kang [START_REF] Kang | Three-dimensional finite element analysis of the cold expansion of fastener holes in two aluminum alloys[END_REF] in his work develop a 3D FEM simulating the entire cold expansion process including hole expansion, recovery and finish reaming. A quarter of 2024-T351 aluminum part is modeled with a thickness of 6.35 mm and a final diameter hole of 7.26 mm. The reaming is simulated through the deactivation of mesh elements in the same way as in this study. Kang provides results of circumferential residual stress along the entry face and the middle section. The result curves present a general shape similar to the one obtained in this study with a compressive peak at the hole edge and a tensile zone far from the hole. However, the shape of the compressive peak varies. This can be explained by the number of points considered for the plotting of the result graphs which depends on the mesh size. The residual stress levels obtained by Kang are of the same order of magnitude as those presented in this study. The differences may be due to the different dimensions used.

Thus, the model developed in this study seems relevant. Besides, it seems more interesting to gain in computing time compared to the 3D model.

Conclusion

An enhancement in the reaming simulation is proposed compared to Achard's simulation strategy [START_REF] Achard | Numerical modelling of the cold expansion process in mechanical stacked assemblies[END_REF]. This new simulation which consists of progressively deactivating mesh elements allows to access to the history of the residual stress state of the part.

The results obtained show only a small change in the residual stress state of the part during the reaming simulation. It means the material removal has a weak influence on the residual stress state of the part. Thus, if a significant change in the residual stress state occurs during the reaming process, it should mainly be related to the thermo-mechanical phenomena. Therefore, a new strategy of reaming simulation according to the principle of applying a thermo-mechanical load as in Valiorgue's works [START_REF] Valiorgue | A new approach for the modelling of residual stresses induced by turning of 316L[END_REF] is envisaged in order to take into account the mechanical cutting and thermal heating phenomena.

Nevertheless, the interpretation of the results emphasize a certain change of the residual stress state in the part during the reaming simulation especially at the hole edge area. It means to envisage a more realistic reaming simulation by applying an equivalent thermo-mechanical load on the hole edge, this one should not be dissociated from the material removal. Thereby, in the envisaged reaming simulation, the loading and the material removal will be applied synchronously and gradually in order to achieve an accurate final stress state in the part.

A reaming simulation on a pre-stressed part obtained after drilling simulation is also envisaged in order to study the evolution of the residual stress state in the hole edge zone during the whole machining process of the hole.

Finally, reaming tests with residual stress measurements in the part are envisaged in order to correlate the experimental measurements with numerical results. Residual stress measurements will be performed using the X-ray diffraction method. This technique has the advantage of not being destructive but it is difficult to implement in the case of a textured material. This is the case for the studied specimens which have been rolled and which present elongated grains in the rolling direction. For that, a second technique which is destructive will also be used for residual stress measurements: the incremental hole drilling method.

Fig. 1 .

 1 Fig. 1. Residual stress induced by machining (a) and cold expansion process (b) in a 2024-T3 aluminum alloy part [11]

Fig. 2 .

 2 Fig. 2. Finite element model of the cold expansion process.

Fig. 3 .

 3 Fig. 3. Mesh of the numerical model

Fig. 4 .

 4 Fig. 4. Representation of the reaming area and the deactivated elements in the first calculation step for the progressive reaming (Dimensions in mm)

Fig. 5 .

 5 Fig. 5. Radial (a) and circumferential (b) stress states obtained in the part at the end of the cold expansion (1 st column), global (2 nd column) and progressive (3 rd column) reaming simulations

Fig. 6 .

 6 Fig. 6. Evolution of radial (a) and circumferential (b) residual stresses along the entry face (1 st column) and middle section (2 nd column) of the part after cold expansion and global and progressive reaming simulations

Fig. 7 .

 7 Fig. 7. Selected elements locating (center) in mm on the part after cold expansion simulation for constraint results reporting

Fig. 8 .

 8 Fig. 8. Evolution of radial (a,b,c) and circumferential (d,e,f) residual stresses during progressive reaming simulation

Table 1 .

 1 Material characteristics of the numerical model parts

Workpiece (AA 2024 T-351) Steel sleeve (15-5 PH) Steel jaw (AISI 301)

	Young's modulus (MPa)	67 800	194 501	205 000
	Poisson's ratio	0.33	0.272	0.29
	Yield stress (MPa)	340	1024	/


Michel Barbeau

Joaquin Garcia-Alfaro

Evangelos Kranakis

Geocaching-inspired Resilient Path Planning for Drone Swarms

Keywords: Autonomous aerial vehicle, drone formation control, drone swarm, goal location, quadcopter, information sharing, localization, location, path planning

A path planning algorithm for drone swarms is presented. From the outset, none of the drones knows the path and final destination. Together, they collectively determine and unravel step-by-step the waypoints and final destination, resolving a localization problem at each step. It is a sharedinformation path planning algorithm. The algorithm is faulttolerant and resilient to drones falling victim of attacks to their positioning system. It is shown that correctly functioning drones navigate the path provided that the number of faulty drones is less than n-d 2 , where n is the total number of drones and d, equal to two or three, is the dimension of the space navigated by the drones. We validate the algorithm with appropriate simulations, implemented over OMNeT++ and GNSSim, which allow building network simulations including GPS attacks (e.g., jamming and spoofing attacks). The OMNeT++ models and GNSSim functions are linked together.

I. INTRODUCTION

We present a path planning algorithm for drone swarms. Together they collectively determine and uncover step-by-step the path and final destination, resolving a localization problem at each step. Ignorance about the path and final destination is a desired feature to collectively handle fault-tolerance and resilience to attacks.

Our work finds inspiration in geocaching, a well-known outdoor recreational activity, in which participants use a GPS receiver as well as a variety of scene cues (e.g., clues and references to landmarks) and other navigational techniques in order to collectively hide and seek objects while at the same time navigating a waypoint trajectory.

The algorithm presented allows for shared-information path planning through waypoints and is fault-tolerant and resilient to drone attacks. More specifically, non-faulty drones can correctly navigate the path provided that the number of faulty drones is less than n-d 2 , where n is the total number of drones and d is the dimension of the space navigated (d = 2, 3).

Section II surveys related work. Section III presents our path planning algorithm. Fault tolerance and resilience to malicious attacks are discussed in Section IV. Section V validates our results with simulations. Section VI concludes the paper.

II. RELATED WORK

Path planning for drones involves making them move from one point to another, while avoiding obstacles, following walls or moving consistently as a team, according to a pattern. Radmanesh et al. have published a survey on path planning with obstacle avoidance for drones [START_REF] Radmanesh | Overview of path-planning and obstacle avoidance algorithms for UAVs: A comparative study[END_REF]. Path planning with a coverage goal is discussed in a survey authored by Otto et al. [START_REF] Otto | Optimization approaches for civil applications of unmanned aerial vehicles (UAVs) or aerial drones: A survey[END_REF]. Zhao et al. survey path planning involving computational intelligence [START_REF] Zhao | Survey on computational-intelligencebased UAV path planning[END_REF]. General path planning approaches find inspiration in ideas developed for classical robotics, such as the ones using artificial potential function [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF], random trees [START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planning[END_REF] or Voronoi diagrams [START_REF] Lavalle | Planning Algorithms[END_REF]. Path planning may be considered together with team work and formation control [START_REF] Turpin | Capt: Concurrent assignment and planning of trajectories for multiple robots[END_REF]. Some approaches have been adapted to quadcopters [START_REF] Rizqi | Path planning and formation control via potential function for uav quadrotor[END_REF], as well. Similarly to other path planning algorithms, the drones have a common goal, i.e., a location or coverage of an area. In contrast to the others, in our algorithm the drones do not know what the exact goal is, i.e., the final location, until the very end.

In the geocaching game analogy, a participant navigates to a location and tries to find a hidden container. It is given geographical coordinates and details about the location, clues and references to landmarks. In our algorithm, the geocaching game analogy is present when a drone participates to the resolution of a waypoint. For an individual drone, it means finding a new point, which is used with other points found by other drones. It is given a relative starting point and a translation, which, as in geocaching, may also be clues and references to landmarks. For that aspect, scene analysis, landmark search and navigation methods for drones can be used in combination with our algorithm [START_REF] Fuentes | Entropy-based search combined with a dual feedforward-feedback controller for landmark search and detection for the navigation of a UAV using visual topological maps[END_REF]- [START_REF] Tariq | Non-GPS positioning systems: A survey[END_REF].

Our research is also related to drone formation control. The authors of [START_REF] Park | Drone formation algorithm on 3d space for a drone-based network infrastructure[END_REF] assume that the signal propagation model of the drone has the shape of a sphere and analyze network capacity allocation in drone-based network infrastructure; they propose a drone formation algorithm that determines the 3D geographic location of each drone. In [START_REF] Shrit | A new approach to realize drone swarm using ad-hoc network[END_REF], the authors show how to operate a swarm by human piloting a drone (the leader) while the remaining followers are autonomous; they propose a solution in order to synchronize and orchestrate a swarm of drones, based only on ad hoc communications to position drones. In [START_REF] Duan | ?hybrid particle swarm optimization and genetic algorithm for multi-uav formation reconfiguration[END_REF], the authors formulate the multi-UAV formation reconfiguration problem as an optimal control problem with dynamical and algebraic constraints and provide a hybrid particle swarm optimization and genetic algorithm.

III. SHARED-INFORMATION PATH PLANNING

A path is represented as a number of waypoints. At every waypoint, the drones collectively determine the next step by solving a localization problem. In the sequel, we discuss the localization problem, waypoint representation and path modelling.

A. Localization Problem

Let n be a positive integer. We define the localization problem as the collective solution by n drones of the location of an unknown point Q. In the two-dimensional Euclidean space, the drones find Q on the perimeter of a circle S, see Figure 1. Q is the point on the intersection of the perimeter of circle S and straight line denoted by L( c, v). Vector c specifies the centre of the circle S. Unit vector v is a direction vector originating in the centre of the circle. The location of Q is where the supporting line of v intersects with S. It is represented by a pair of coordinates (q 1 , q 2 ). Initially, the circle S is unknown to the drones. Let i denote a drone index, with i ∈ {1, 2, . . . , n}. Each drone i is given the direction vector v and information that determines a point P i on the perimeter of the circle S. The information about the point P i is owned only by drone i and not shared with other drones. For brevity, in the sequel, we abbreviate the localization problem as L(P 1 , . . . , P n ; v).

c v O Q = (q 1 , q 2 ) S L( c, v) r
Recall from Euclidean geometry that a circle is uniquely determined by three pairwise different points on its perimeter. Therefore if n ≥ 3, then there is a unique circle S determined by the n points. Upon determination of the points P 1 , . . . , P n , the drones broadcast their positions to each other. They can all determine Q as the point at which the straight line L( c, v) intersects with the perimeter of circle S. The localization problem is thus solved.

B. Representing Waypoints

Building upon the localization problem L(P 1 , . . . , P n ; v), we show how to represent any desired movement of a drone swarm. To this end, we show for any two different points Q, Q how to map Q to Q . Suppose that the drones have

C. Representing Paths

Using waypoint representation, we determine a sequence of consecutive steps forming a path

Q 0 , Q 1 , . . . , Q k-1 = Q.
The drones start from an initial configuration and inductively reach a specific destination as a group. The drones start from point Q 0 determined by the solution of the initial instance of the localization problem L(P 1,0 , . . . , P n0,0 ; v 0 ), see Figure 3.

Given that the solution of the i-th intermediate localization problem L(P 1,i , . . . , P ni,i ; v i ) is the i-th intermediate point Q i , they solve the i + 1st instance of the localization problem L(P 1,i+1 , . . . , P ni+1,i+1 ; v i+1 ) and obtain the i

+ 1st interme- diate point Q i+1 . The final desired destination Q = Q k-1 is reached when i = k -1.
Observe that in each iteration of the localization procedure some of the drones may malfunction, due to faults or malicious attacks of their positioning system. They may well withdraw from the localization problem. Hence, the number of drones participating in consecutive rounds may change, e.g., from the i-th to the i + 1-st the number of participating mobiles may change from n i to n i+1 , respectively. In case of an attack, the number of victims depends on the disruption power that the adversary has. Fault tolerance and resilience to attacks is discussed further in the upcoming section.

IV. FAULT TOLERANCE AND RESILIENCE TO ATTACKS

We analyze the tolerance to faults and resilience to malicious attacks of the shared-information path planning algorithm. Hardware failures may occur. Besides, drone positioning systems may be attacked, e.g., GPS spoofing attacks can be perpetrated or visual references can be perturbed. Another plausible attack is the capture and reverse engineering of drones by an adversary. In Euclidean space, information from three different non-faulty drones is sufficient to determine the path. The risk of success for this attack is low assuming drone capture and reverse engineering take much longer time than the mission of the drones.

The shared-information path planning algorithm, described in Section III, can be generalized to a d-dimensional space, where d is two or three. From Euclidean geometry, it is known that a sphere in d-dimensional space is uniquely determined by d + 1 pairwise different points on its outline (perimeter or surface, respectively). Therefore if n ≥ d + 1, then there is a unique geometric object S determined by the n points. When d = 2 the object is a circle, and when d = 3 it is a sphere. It follows that if the drones broadcast their own positions to each other then they can all determine Q as the point at which the line denoted by L( c, v) intersects with the outline of object S, where c is its centre.

Suppose that at most f of the n drones are faulty in the sense that they do not provide their correct location to the rest of the swarm either because they are malfunctioning or their positioning system is the victim of an attack. The remaining n -f drones are reliable and broadcast their correct location to the peers.

Is it possible for the reliable drone to correctly solve the localization problem and point Q despite the presence of the f faults? Observe that communication may be reliable but data may not be because some drones may have their local system hacked. The main point here is that although the drones know from the broadcast each other's coordinates they do not know in advance which among them are reliable and which are not. In particular, although the drones can broadcast their coordinates to the other drones correctly and compute the correct set of coordinates of all the drones, they have no way of verifying their authenticity. Our motivation is inspired from the work of Kleinberg [START_REF] Kleinberg | This Will Make You Smarter: New Scientific Concepts to Improve Your Thinking[END_REF] which proposed the "circle" paradigm and looked at the case of non-faulty drones (f = 0).

We show that the localization problem is correctly solved assuming that the majority of the drones are reliable. We make this precise by proving the following theorem. (Part 1) First, let us consider the case where n ≤ 2f + 2 (see Figure 4 for an arrangement of n = 8 drones with f = 3 faulty drones; black dots represent reliable drones and black squares unreliable drones). If n -f ≤ 2 then there is more than one circle passing through the reliable drones and so the point Q cannot be uniquely determined. So, without loss of generality we may assume that n -f ≥ 3. In this case, the reliable drones determine a unique circle, say S. The remaining f faulty drones may therefore use any two among the n -f drones (a total of n-f 2 possibilities) to form a different circle S that also specifies a different point Q. This gives rise to two circles. The first one S determined by the f +2 drones and the second S by the n -f drones. Moreover, since n -f ≤ f + 2 the drones cannot determine the circle formed by the reliable drones. (Part 2) Clearly, the n -f reliable drones lie on the outline of a d-dimensional object (here we consider the case of d = 2). We would like to show that if n -f > f + 2, then this circle is unique and the faulty drones cannot determine a different circle. If n > 2f + 2 (see Figure 5 for an arrangement of n = 11 drones with f = 3 faulty; black dots represent reliable drones and black squares unreliable drones) then the reliable n-f drones determine a unique circle that involves all the nonfaulty drones. This is because n -f > 2 and from Euclidean geometry we know that a circle is uniquely determined by three pairwise different points on its perimeter.

The f faulty drones may try to produce a different circle S so as to mislead the reliable drones. To this end the unreliable drones can only co-operate with at most two reliable drones (without their knowledge) in order to form a different circle that uses f + 2 drones. Clearly, the faulty drones cannot cooperate with three or more reliable drones, because we know from Euclidean geometry that three points in the plane determine a unique circle, which is also the circle determined by the n-f reliable drones. However, by assumption n > 2f +2. Therefore the correct circle is the one determined by the n-f reliable drones which by our assumption are in the majority since n -f > f + 2.

V. SIMULATIONS

A. Simulation tools

In order to validate our work, the path planning algorithm presented in Section III has been integrated into the OMNeT++ (Objective Modular Network Testbed in C++) framework [START_REF] Varga | An overview of the OMNeT++ simulation environment[END_REF], [START_REF]The OMNeT++ Network Simulation Framework[END_REF]. The implementation leverages a series of shared libraries over INET [START_REF]Last Access[END_REF], an OMNeT++ model suite for the simulation of wired, wireless and mobile networking protocols, including UDP, TCP, SCTP, IP, IPv6, Ethernet, PPP, 802.11, MPLS, OSPF, and many others.

Together, OMNeT++ and INET allow the development of discrete event simulations, in which models are described in the C++ language [START_REF] Heijmans | OMNeT++ Discrete Event Simulation System -User Manual[END_REF] and network properties are specified in the NED (NEtwork Description) language [START_REF]NED language overview[END_REF]. In a nutshell, INET simulations consist of networking modules that communicate with each other through message passing (either via OMNeT++ gates, or through direct messages). The use of add-ons, e.g., via shared libraries, makes possible complementing OMNeT++ code with functionality from other simulation engines. This add-on architecture can be used to extend the INET framework for the development of satellitebased communication simulations. As a result, we can use frameworks like OS3 (Open Source Satellite Simulator) [START_REF] Niehoefer | The CNI open source satellite simulator based on OMNeT++[END_REF] and GNSSim (Global Navigation Satellite Simulation System) [START_REF] Jahan | Gnssim: An open source gnss/gps framework for unmanned aerial vehicular network simulation[END_REF], [START_REF] Javaid | Analysis of global positioning system-based attacks and a novel global positioning system spoofing detection/mitigation algorithm for unmanned aerial vehicle simulation[END_REF]. Together, they allow the development of OMNeT++ simulations combining satellite and navigation communication protocols, and attacks like GPS jamming and GPS spoofing.

OS3 includes functionality to import real satellite tracks and weather data. This functionality is then used by GNSSim to implement a series of navigation services, such as GPS, GLONASS, Galileo, and BeiDou [START_REF] Hofmann-Wellenhof | GNSS-global navigation satellite systems: GPS, GLONASS, Galileo, and more[END_REF]. The implementation of two GPS related attacks, GPS spoofing and GPS jamming, are also included in the last version of GNSSim [START_REF] Javaid | Analysis of global positioning system-based attacks and a novel global positioning system spoofing detection/mitigation algorithm for unmanned aerial vehicle simulation[END_REF]. These two attacks are used in our simulation experiments to validate the feasibility of our path planning algorithm to correct the intentional disruption of a series of malicious drones in a swarm, under a linear path mobility experiment. The threat model consists of zombie drones, under the control of a remote adversary conducting a cyber-physical attack [START_REF] Rubio-Hernan | Revisiting a Watermark-Based Detection Scheme to Handle Cyber-Physical Attacks[END_REF]- [START_REF]Adaptive control-theoretic detection of integrity attacks against cyber-physical industrial systems[END_REF]. The zombie drones perpetrate GPS attacks on their neighboring drones of the swarm in order to disrupt their mission. The victims of the attacks are captured by the zombies and taken away from the swarm.

B. Simulation scenario and early results

We report next our simulation scenario and some numeric simulation results. Firstly, we summarize below the main entities and details underlying our simulation scenario:

• Zombie drones are drones from the initial swarm which get hacked by a remote adversary (e.g., a remote attacker conducting a cyber-physical attack [START_REF] Rubio-Hernan | Revisiting a Watermark-Based Detection Scheme to Handle Cyber-Physical Attacks[END_REF]- [START_REF]Adaptive control-theoretic detection of integrity attacks against cyber-physical industrial systems[END_REF]). Zombie drones become perpetrators of GPS attacks (i.e., jamming and spoofing attacks) with the objective of disrupting the remaining drones of the swarm. • There is a variable number of zombie drones under the control of a remote adversary. Each zombie drone captures a variable number of neighboring drones (i.e., non-zombie drones that are within the coverage range of the GPS attacks of the zombie drones). • The victims of the zombie drones are pulled out from their normal trajectory, i.e., from the swarm. They are referred to as captured drones. • Each attack is perpetrated only at a given waypoint, i.e., at each waypoint, new drones are compromised, and may become either zombie drones (i.e., victims of the

A Q 0 Q 1 Q 2 ... Q k B
(a) Simulation scenario baseline (without attacks) remote adversary) or captured drones (i.e., victims of the zombie drones). The remaining drones which successfully escape from the remote adversary and zombie drones, are denoted as survivors.

A Q 0 Q 1 Q 2 Q k B attack
• Each zombie drone leaves the swarm in order to conduct the GPS attacks (both GPS jamming and spoofing), as indicated in [START_REF] Javaid | Analysis of global positioning system-based attacks and a novel global positioning system spoofing detection/mitigation algorithm for unmanned aerial vehicle simulation[END_REF]. In other words, the zombie drones never recover from their attack. Hence, the attack perpetrated by the zombie drones is self-destructive.

Figure 6 depicts a more elaborated representation of our simulation scenario. In Figure 6 (a), we assume a swarm of n drones. All the drones start the mission at point A. They must cooperate to reach point B as final destination, after visiting k intermediate waypoints (i.e., Q 0 , Q 1 , Q 2 , . . . , Q k ). They execute the path planning algorithm presented in Section III in order to identify the waypoints and complete the mission. Under the absence of the adversary, all the n drones reach point B, in the end.

Parameter Value

Mobility type of satellites SarSGP4Mobility

Mobility type of drones PathPlanningMobility . Sample visualization capture of our ongoing simulation testbed using OMNeT++ [START_REF] Varga | An overview of the OMNeT++ simulation environment[END_REF], OS3 [START_REF] Niehoefer | The CNI open source satellite simulator based on OMNeT++[END_REF] and GNSSim [START_REF] Jahan | Gnssim: An open source gnss/gps framework for unmanned aerial vehicular network simulation[END_REF], [START_REF] Javaid | Analysis of global positioning system-based attacks and a novel global positioning system spoofing detection/mitigation algorithm for unmanned aerial vehicle simulation[END_REF]. Some additional videocaptures of our ongoing simulation testbed are available at http://j.mp/gnssimuav.

In Figure 6 (b), the previous scenario changes due to the presence of a remote adversary which is capable of infecting some of the drones in the swarm, with the objective of disrupting their mission. Such drones, i.e., the zombie drones, are now under the control of the remote adversary. Commanded by the adversary, the zombie drones perpetrate the GPS attacks (both jamming and spoofing attacks) in order to disrupt their neighboring drones (i.e., drones under the coverage zone of the zombies). The victims of the GPS attacks, i.e., the captured drones, get disrupted and also fail at reaching the following waypoint. In the end, they fail their mission. The remaining survivor drones (i.e., those not affected by the adversary or the zombie drones) succeed at unraveling the intermediate waypoints and successfully reach the final destination (point B). The defense strategy depicted in Figure 6 (b) also assumes that the victims of the adversary, i.e., the zombie and captured drones, get lost forever.

Figures 7 and8 show some other important parameters used in our OMNeT++ simulations, as well as some captures of the OMNeT++ GUI using OS3 [START_REF] Niehoefer | The CNI open source satellite simulator based on OMNeT++[END_REF] and GNSSim [START_REF] Jahan | Gnssim: An open source gnss/gps framework for unmanned aerial vehicular network simulation[END_REF], [START_REF] Javaid | Analysis of global positioning system-based attacks and a novel global positioning system spoofing detection/mitigation algorithm for unmanned aerial vehicle simulation[END_REF]. Some simulation videocaptures are avail- able at http://j.mp/gnssimuav. Figure 9 plots the numeric results. At every waypoint, the number of zombies is a random variable that follows a Poisson distribution with parameter λ 1 . The number of victims per zombie is also a random variable that follows a Poisson distribution with parameter λ 2 . The following λ 1 , λ 2 pairs have been simulated: (1, 1), [START_REF] Radmanesh | Overview of path-planning and obstacle avoidance algorithms for UAVs: A comparative study[END_REF][START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planning[END_REF], [START_REF] Otto | Optimization approaches for civil applications of unmanned aerial vehicles (UAVs) or aerial drones: A survey[END_REF][START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planning[END_REF] and [START_REF] Otto | Optimization approaches for civil applications of unmanned aerial vehicles (UAVs) or aerial drones: A survey[END_REF][START_REF] Jayatilleke | Landmark-based localization for unmanned aerial vehicles[END_REF]. For a mission to succeed, the number of drones at the end (at point B) must be greater than zero, since no waypoints are assumed to be computed in point B, as depicted in Figure 6. For given values of parameters λ 1 and λ 2 , the simulation confirms that the success rate grows consistently with the number of drones; while greater values for the parameters λ 1 and λ 2 translate in higher impact and less chances of mission success.

VI. CONCLUSION

We have presented an information sharing path planning algorithm for drone swarms. The drones in the swarm do not have any knowledge about the entire path, nor the final destination. They must determine and unravel a series of intermediate steps by conducting a collective process, resolving a localization problem at each step. Likewise, they must work collectively. At each step, the algorithm executed by the drones is geocaching inspired. Geocaching is a well-known outdoor recreational activity in which participants use GPS receivers (or any other alternative navigational technique) as well as a variety of scene cues (e.g., visual clues or references to landmarks) to collectively seek objects while navigating a waypoint trajectory. In our work, the series of waypoints allows the drones to build a shared-information path planning process. The goal is to have a fault-tolerant process, resilient to traditional in drone scenarios, such as GPS jamming and spoofing. Our algorithm have been formalized, analyzed and validated via numeric simulations.

Fig. 1 .

 1 Fig. 1. In Euclidean space with origin O, the point Q is on the intersection of the line of action of vector v, i.e., L( c, v), and perimeter of the circle S.

Fig. 2 .

 2 Fig. 2. Given points Q, Q a unique circle can be determined. It is formed by the new positions of the drones (depicted as squares) in such a way that the point Q lies on its perimeter.found the location Q in a previous round, i.e., solution of L(P 1 , . . . , P n ; v). The drones resolve a new instance of the localization problem L(P 1 , . . . , P n ; v), depicted in Figure2. The hollow squares represent the new positions that the drones have to find, i.e., P 1 , . . . , P n . In reference to point Q, point P i results from the transformation c + re jθi , such that θ i are n pairwise different angles in the interval [0, 2π], j = √ -1 and r is the radius of the circle. They locate a new point Q . Each drone needs to know the transformation c+re jθi before going on a mission. Assuming drones have scene analysis capability, this transformation may take the form of visual cues, such as landmarks or buildings.

Fig. 3 .

 3 Fig. 3. A path consisting of four hops, as traversed by the drones. The drones start from point Q 0 . In each instance, they use a direction vector v to compute an intermediate destination point Q i on the perimeter of a circle. They determine their new positions and again compute the next intermediate destination using the next destination vector. This is repeated until the final destination point Q is reached.

  1 , . . . , P n ; v) for n drones at most f of which are faulty. Let v be a direction vector in d-dimensional space known to all the drones. 1) If n ≤ 2f + d, then there is a d-dimensional object S with center c on whose outline at least f + d (which is ≥ n -f ) of the drones can be located in such way that the intersection of the straight line L( c, v) and outline of S is not the correct destination point Q. 2) If n > 2f + d, then there is a unique object S with center determined by c on whose outline at least n -f (which is > f + d) of the drones are located and the correct destination point Q can be uniquely determined as the intersection of the straight line L( c, v) with the outline of S, despite the presence of the f faulty drones. Proof: We give the proof only for points in the plane (d = 2). However, the general case in d-dimensional space is merely a reformulation of this special case.

Fig. 4 .

 4 Fig. 4. An arrangement of n = 8 drones with f = 3 faulty. Black dots represent reliable drones and black squares faulty drones.

Fig. 5 .

 5 Fig. 5. An arrangement of n = 11 drones with f = 3 faulty. Black dots represent reliable drones and black squares unreliable drones.

Fig. 6 .

 6 Fig. 6. Simulation scenario. (a) depicts a swarm of n drones, starting at point A and cooperating (via the path planning algorithm presented in Section III) to reach point B, after visiting k intermediate waypoints (i.e., Q 0 , Q 1 , Q 2 , . . . , Q k ). (b) depicts a series of zombie drones (under the control of the remote adversary) and captured drones (disrupted by GPS jamming and spoofing attacks perpetrated by the zombie drones). Both victim types in (b) fail at reaching the waypoints of the path, and get lost forever. Only a few survivor drones from the original swarm succeed at reaching the final destination.

Fig. 7 .

 7 Fig. 7. GNSSim parameters used in our simulations.

Fig. 8

 8 Fig.8. Sample visualization capture of our ongoing simulation testbed using OMNeT++[START_REF] Varga | An overview of the OMNeT++ simulation environment[END_REF], OS3[START_REF] Niehoefer | The CNI open source satellite simulator based on OMNeT++[END_REF] and GNSSim[START_REF] Jahan | Gnssim: An open source gnss/gps framework for unmanned aerial vehicular network simulation[END_REF],[START_REF] Javaid | Analysis of global positioning system-based attacks and a novel global positioning system spoofing detection/mitigation algorithm for unmanned aerial vehicle simulation[END_REF]. Some additional videocaptures of our ongoing simulation testbed are available at http://j.mp/gnssimuav.

Fig. 9 .

 9 Fig. 9. Simulation assuming the number of zombies per attack follow a Poisson distribution with parameter λ 1 and number of victims per zombie follows a Poisson distribution with parameter λ 2 .

Acknowledgements -The authors gratefully acknowledge financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC), the Cyber CNI Chair of the Institut Mines-Télécom (cf. http://chaire-cyber-cni.fr/), and the European Commission (H2020 SPARTA project, under grant agreement 830892).