Jérémie Derré
email: jeremie.derre@onera.fr

Frank Simon
email: frank.simon@onera.fr

CONCEPT OF "FRACTAL" HELICOPTER TRIM PANEL

The reduction of structural vibrations and acoustic radiation has become increasingly challenging due to mass reduction and the use of lightweight materials, especially in the aerospace industry. This paper presents the vibroacoustic behaviour of sandwich constructions locally overloaded by a pre-fractal mass distribution. The structural mechanical model consists in a homogenised beam in simple bending, approximated by finite differences and solved numerically to form a modal basis (eigenvalues and eigenvectors). The self-similar mass distribution, inspired from the Cantor set, is presented and integrated using local overloadings within the structure, which leads to a strong localisation of the beam transverse displacement. An analysis of the integrated density of states points out localised modes and their types (central or lateral localisation). Finally, acoustic coupling is modelled through a radiation resistance matrix from a discretisation of the structure into rectangular elementary radiators. Simulations exhibit a smaller radiation efficiency and lower modal critical frequencies. These results are encouraging for reducing the acoustic radiation of lightweight structures.

NOMENCLATURE [I N-1 ] identity matrix of dimension (N -1) [M]

symmetric matrix [R] radiated resistance matrix [Z] complex impedance matrix α loc local overloading mass ratio α tot global overloading mass ratio Θ N-1 null vector of dimension (N -1) p e radiator element wall pressure U transverse displacement vector v e radiator element surface velocities ω angular frequency (rad.s -1 ) ρ homogenised material density (kg.m -3 ) ρ 0 air density at NTP (kg.m -3 ) σ acoustic radiation efficiency σ m mode m acoustic efficiency c 0 sound velocity in air at NTP (m.s -1 ) D approximate total bending stiffness (Pa.m 4 ) D sand equivalent bending stiffness (Pa.m [START_REF] Langley | Wave transmission through onedimensional near periodic structures: optimum and random disorder[END_REF]

) u(x) beam transverse displacement (m) u(x i ) point x i displacement (m) u z (x, t) time dependant transverse displacement (m) v 2
n (ω) normal averaged quadratic surface velocity x i segment delimiting point

INTRODUCTION

In the aerospace industry, sandwich honeycomb core panels have been widely used for decades within secondary aircraft parts (e.g., trim panels), because of their high static bending stiffness at a lesser mass cost. The type of sandwich considered in this paper consists of two thin but sufficiently stiff laminates bonded to each side of a relatively thick and lightweight core. Nowadays, they are still of great interest and developments are continuously made to reinforce them, in order to finally get them integrated as primary parts (NHIndustries NH90 for instance, fuselage or structural elements). The weight reduction is unfortunately responsible for a rise in the panel acoustic transmission over the critical frequency, [START_REF] Simon | Optimization of sandwich trim panels for reducing helicopter internal noise[END_REF] which has to be compensated.

In noise reduction, several strategies are based on the passive principle of structural vibration reduction to minimize the radiated noise. Damped viscoelastic material can be added onto the sandwich external skin, either covering all the panel surface or localised as small patches. Another technique used in the helicopter industry consists in filling out the honeycomb core with small elements such as empty spheres made of elastomer. Nonetheless, instead of reducing the global vibration amplitudes, the wave propagation in the vibrating material can be modified by introducing additional elements within the structure. The additional elements (irregularities, heterogeneities or others) can be distributed within the initial structure, following a periodic, [START_REF] Mead | Plates with regular stiffening in acoustic media: Vibration and radiation[END_REF][START_REF] Mead | Wave propagation in continuous periodic structures : Research contributions from Southampton, 1964-1995[END_REF] random, [START_REF] Langley | Wave transmission through onedimensional near periodic structures: optimum and random disorder[END_REF] as well as selfsimilar patterns. [START_REF] Alexander | Density of states on fractals : "fractons[END_REF][START_REF] Bertaud | A numerical study of 1D self-similar waveguides: Relationship between localization, integrated density of state and the distribution of scatterers[END_REF][START_REF] Bertaud Du Chazaud | Elastic waves in 2D self-similar structures : localisation, vibrational integrated density of state and distribution of the scatterers[END_REF] The structure itself can be composed of several sub-parts connected together to form a multi-span system. [START_REF] Bouzit | Localization of vibration in disordered multi-span beams with damping[END_REF][START_REF] Bouzit | An experimental investigation of vibration localization in disordered multi-span beams[END_REF] In both cases, material and geometrical discontinuities create diffraction, reflection, transmission and even dissipation phenomena, which lead to the apparition of stop-bands and pass-bands in the frequency content, as well as a concentration of the mechanical energy and thus the localisation of the structure mode shapes. Therefore, the structure can be considered as a mechanical filter for wave propagation.

Recent studies have shown encouraging mechanical results on self-similar configurations such as string [START_REF] Bertaud | A numerical study of 1D self-similar waveguides: Relationship between localization, integrated density of state and the distribution of scatterers[END_REF] and membrane. [START_REF] Bertaud Du Chazaud | Elastic waves in 2D self-similar structures : localisation, vibrational integrated density of state and distribution of the scatterers[END_REF] Simulations identify localisation phenomena of the transverse displacement, confirmed by some experimental data. Integrated density of states are analysed, as well as the heterogeneity distribution role and the corresponding modification on modal shapes. The potential of such configurations motivates its application to the helicopter industry. Trim panels are responsible for the most part of radiated noise within the aircraft cabin. These panel surface areas are about a square meter. The global aim is to reduce the generated sound of vibrating panels, in the frequency bands of interest (main gear box noise or aerodynamic excitation for instance).

In this paper, the key idea is to use the geometrical cell network of the honeycomb core, and to fill out some cells by following a self-similar scheme. Indeed, the transition between two media of different natures and densities (empty cell and cell fully filled of spheres) could be seen for a bending wave as a mechanical impedance discontinuity, leading to multiple scattering phenomena. Thus, the idea is to use the properties of the full cell network to create localised modes, resulting of wave pseudo-stationarity between heterogeneities.

There are two main particularities to the present study. The first one is the application of the well-known theory of pre-fractal distribution to a composite material, i.e. a sandwich beam. The second one is that the loading is realized within the material core, meaning that the structure load-bearing capability is not altered, which is of primordial importance from a mechanical point of view. Sample panels have already been realized and characterized in transmission with promising results. [START_REF] Onera/Ateca | Panneau insonorisant[END_REF] From an industrial point of view, the manufacturing of such panels has been carried out manually without any specific issues.

The paper is organized in five sections. Following the introduction, the developed mechanical model is first presented for a monodimensional homogenised structure in Section 2. In Section 3, the heterogeneity integration into the model is described. Then, the self-similar distribution inspired from the Cantor set is briefly introduced, followed by simulation results of a localised mode for several orders of the self-similar scheme. Section 4 analyses the integrated density of states, identifying localised and classical extended modes. The vibroacoustic behaviour is presented in Section 5, and modal acoustic radiation efficiency is simulated and discussed. Finally, the conclusion and potential extensions are summarized in Section 6.

MECHANICAL MODEL OF A HOMOGENISED BENDING BEAM

The first step is a presentation of the mechanical model of a non-overloaded sandwich structure governed by bending dynamics.

Composite homogenisation

A typical sandwich trim panel construction is a core material trapped between two laminate sheets of glassfiber. Herein, the core is made of an aramid fibre paper impregnated with a heat resistant phenolic resin and arranged in honeycomb hexagonal cells. The purpose of such a core is to maintain distance between the laminates and to some extent, to resist shear deformation. [START_REF] Backström | Modelling the vibration of sandwich beams using frequencydependent parameters[END_REF] In the modelling of acoustic radiation created by mechanical vibration (fluid-structure interaction), bending waves have a predominant role in the sound frequency domain. [START_REF] Fahy | Sound and structural vibration -Radiation, transmission and response[END_REF] Under Nilsson assumptions [START_REF] Anders | Wave propagation in and sound transmission through sandwich plates[END_REF] (homogeneity and isotropy of all layers and skins assumed thin), an equivalent bending stiffness of the whole structure can be derived using Huygens theorem. The resulting equivalent bending stiffness D sand is given by

D sand = E c I c + E s I s , = E c L y h 3 c 12 + 2 E s L y h 3 s 12 + L y h s h c +h s 2 2 , (1)
where Young's modulus and second moment of area (of the cross section) of the core (c) and of the skins (s) are identified respectively by E and I, as well as thicknesses of elements h and the beam width is L y .

Kim and Hwang [START_REF] Kim | Effect of debonding on natural frequencies and frequency response functions of honeycomb sandwich beams[END_REF] examined the validity conditions of the approximation used in such a model. Bending stiffness of the skins projected to the neutral line of the sandwich beam is predominant and thus, the approximate total bending stiffness could be written as

D = E s L y h s (h c + h s ) 2 2 , (2)
where the thickness (h c + h s ) confers the ensemble a large equivalent bending stiffness, called sandwich effect. In the following of the paper, the structure is considered as an homogeneous equivalent material, described with the previous Eq. ( 2).

Due to its complex constitution, the mechanical behaviour of a sandwich beam involves different parts of the structure depending on frequency. However, the global behaviour at high frequencies is governed by the bending dynamics of the laminates, [START_REF] Backström | Modelling the vibration of sandwich beams using frequencydependent parameters[END_REF] which gives consistency to our approach. Indeed, the strong localisation phenomena are expected for high mode orders. As long as the first model is a beam, the frequencies of interest would then be shifted towards to higher frequencies due to its small width.

Simple bending beam

Under the previous hypotheses and approximations, a fourth-order beam theory is used. Considering a straight slender Euler-Bernoulli beam, and timedependency taken as harmonic with angular frequency ω, the transverse displacement of the beam u z (x, t) = u(x) × exp ( jωt) is governed by equilibrium Eq. ( 3)

D ρS d 4 u(x) dx 4 -ω 2 • u(x) = 0 , (3)
with the cross section S = hL y (total beam thickness by its width) and ρ the density of the material.

Discretisation

The choice of derivative approximation and structure discretisation by finite difference (abbreviated FD) method is mainly justified by its convenience and the flexibility of its implementation regarding the different orders of self-similar distributions which are integrated within the structure.

Given a beam under the aforementioned hypothesis, of length L x , the structure is divided into N segments of constant length h x as L x = h x × N. A set of (N + 1) points (x i ) delimiting the segments is defined as x i = (i -1) × h x , with x 1 and x N+1 the beam extremities. Finally, the transverse displacement at position x i is written as u(x i ) = u i .

Following the classical FD technique, such as Subrahmanyam and Leissa, [START_REF] Subrahmanyam | An improved finite difference analysis of uncoupled vibrations of cantilevered beams[END_REF] Eq. ( 3) is approximated with a second order central scheme, which gives the equilibrium around position x i as

D ρS h 4 x u i-2 -4u i-1 + 6u i -4u i+1 + u i+2 -ω 2 • u i = 0 . (4)

Boundary conditions

In this study, the beam is considered as simply supported in its width and free in its length. To integrate boundary conditions, a method based on the creation of off-domain supplementary points is used. These points are linked by FD relationships and substituted, then integrated to the equations of the physical domain.

Matrix notation

Once the boundary conditions are integrated, the following system (5) is obtained, with (N -1) equations

                                                                           (i = 2) : D h 4 x [5u 2 -4u 3 + u 4 ] -ω 2 n ρS • u 2 = 0 , (i = 3) : D h 4 x [ -4u 2 + 6u 3 -4u 4 + u 5 ] -ω 2 n ρS • u 3 = 0 , ... (i) : D h 4 x [u i-2 -4u i-1 + 6u i -4u i+1 + u i+2 ] -ω 2 n ρS • u i = 0 , ... (i = N-1) : D h 4 x [u N-3 -4u N-2 + 6u N-1 -4u N ] -ω 2 n ρS • u N-1 = 0 , (i = N) : D h 4 x [u N-2 -4u N-1 + 5u N ] -ω 2 n ρS • u N = 0 . (5)
The ω n are discrete. Indeed, the system dimension is n and thus the number of degrees of freedom is fixed. The previous system can be rewritten as Eq. ( 6), where matrix [M] is a symmetric pentadiagonal matrix, [I N-1 ] is the identity matrix of dimension (N -1), U is the vector of transverse displacements and Θ N-1 the null vector of (N -1) components

D ρS .h 4 x M -ω 2 n I N-1 U = Θ N-1 . (6)
Eq. ( 6) is an eigenvalue problem, which is solved to form a modal basis of the structure (eigenvalues and eigenvectors for natural frequencies and mode shapes respectively).

MECHANICAL BEHAVIOUR OF A STRUCTURE LOADED BY A SELF-SIMILAR DISTRIBUTION

This part aims at describing how the heterogeneities are integrated within the structure and how they are modelled in the FD scheme. The self-similar distribution is then presented. Finally, the shape of a given localised mode is computed for different pre-fractal orders.

Integration of heterogeneities following a pre-fractal scheme

The patent [START_REF] Onera/Ateca | Panneau insonorisant[END_REF] presents the concept of elements distributed following a self-similar scheme in the core of composite sandwich panels. Preliminary experiments featured localisation phenomena and identified the added mass as the key element of overloaded structures.

Integration of a single heterogeneity

Given that the main mechanism identified is mass adding, it is modelled as the local overload of a discrete element, and can represent a heterogeneity within the honeycomb core. As mentioned above, transition between two different media creates an impedance breakdown. Therefore, the addition of a small mass M over on a discrete element k is modelled as in Eq. ( 7)

D h 4 x u k-2 -4u k-1 + 6u k -4u k+1 + u k+2 -ω 2 n ρS + M over h x • u k = 0 . (7)
The local overloading mass ratio is defined as Eq. ( 8)

α loc = ρS ρS + M over h x , (8)
leading to the following formulation

α loc D ρS h 4 x u k-2 -4u k-1 + 6u k -4u k+1 + u k+2 -ω 2 n • u k = 0 , (9)
which is the k-line of matrix [M] multiplied by α loc in Eq. ( 6).

Self-similar scheme

The studied structure is the homogenised beam as described previously, overloaded by a self-similar distribution inspired by the Cantor set. This iterative construction process has an homothetic ratio of three: at each iteration (order), the element is sub-divided into three parts and the central one is removed. The first three orders of this motif are presented on Figure 1, as well as order 0, which is the initial structure. The distribution used is strongly inspired from this set, except that only the extremities of the sub-sets are marked and correspond to local overloadings. Is called fractal the object obtained when a self-similar pattern construction is infinitely repeated, and this is why in the present case the motif is called pre-fractal with a finite iteration order n frac . Based on the construction structure, [START_REF] Bertaud | A numerical study of 1D self-similar waveguides: Relationship between localization, integrated density of state and the distribution of scatterers[END_REF] a self-similar beam of pre-fractal order n frac is overloaded by p masses with

p = n frac i=1 2 i = 2 × (2 n frac -1) . (10)

Integration of the pre-fractal distribution

For a given n frac pre-fractal order, p local masses are integrated on the extremities of the sub-sets defined previously. This leads to multiplications of the p corresponding lines of matrix [M] by α loc . It is noticeable that the discretisation in N segments plays a crucial role regarding the size of the discrete elements. Indeed, even though masses are taken as points, they stand for a segment of equivalent length h x . The discretisation is therefore of primordial importance and a key feature in the model for simulations. Figure 2 sketches simply-supported beams with the previously described Cantor-like set, from order 0 to 3. Overloaded masses are represented by black dots.

Localisation of transverse displacement and mechanical behaviour

Figure 3 presents a simulation of shapes of mode 26 in the non-overloaded case and for beams overloaded by self-similar distribution at order 2 and 3, computed with a finite difference method, with N = 243 segments. All quantities of the numerical simulations are summarised in Table 1 on the appendix. Mode identification and numbering corresponds to the number of transverse nodes, as for a classical homogeneous beam. In this paper, classical beam refers to a non-overloaded beam.

The global overloading mass ratio α tot of the beam is defined as the ratio of added overload masses on non- overloaded beam mass as

α tot = p × M over ρL x . ( 11
)
Total beam overload ratio α tot is about 32% at order 2 and about 75% at order 3.

The mode shapes of Figure 3 are normalized by Euclidean norm equals to unity. The displacements of the central part have a significantly larger amplitude than those on lateral side parts. For the second-order structure, a localisation ratio (ratio between maximal amplitudes of the central part and of the lateral parts) close to three can be observed; for the third-order structure, it is about four. It can also be noted that the global behaviour of the curve does not present any singular point. Indeed, bending waves are not completely stopped at the heterogeneities but partially reflected and transmitted.

The local overloading ratio α loc remains constant for the different orders presented and is about 0,071. In the case of a helicopter sandwich material, this corresponds to a core cell filled with a material of density about 3800 kg.m -3 , which lies within the ceramic ranges.

ROLE OF THE PRE-FRACTAL DISTRIBUTION ON INTEGRATED DENSITY OF STATES

The phenomenon of mode localisation can be clearly identified on integrated density of state (abbreviated DOS) curves, as illustrated on Figure 5. • In the low frequency region, the DOS curve of overloaded beam is quasi-identical to the classical beam, meaning a linear curve of constant slope (in logarithmic axis). The only impact of the overloading in these modes is a modification of the y-intercept. The fractal does not influence these modes, which are so-called extended modes. In the case of the order-three Cantor-like distribution, this corresponds to the first 20 modes.

• The first slope breakdown corresponds to the first modes with a localisation of the transverse displacement (mode 21 in the present study case), abbreviated localised modes. Indeed, these modes exhibit a mechanical energy concentration in the central part of the beam. For mode 22, central-to-lateral amplitude ratio is about three.

• The second slope breakdown appears at mode 23, for two modes localised on the non-overloaded lateral areas of the pre-fractal beam. These modes have a strong localisation ratio (more than 3.5) as well as very close eigenvalues (1% relative difference), which means an almost null slope on the DOS curve.

• The third slope breakdown appears for mode 25, linked to new central-part localised modes. The integrated DOS slope follows the same slope augmentation than for the first slope breakdown.

• Finally, mode 27 has a frequency almost equals to the non-overloaded beam mode 27 (about 1% of relative difference). Indeed, the overload distribution has almost no impact on mode shape because over masses are located on vibration nodes, which do not affect the dynamic behaviour.

If Figure 5 is then analysed, plotting integrated DOS of the 80 first modes, the behaviours previously identified duplicate with frequency growth. So the pattern of the 27 first modes of Figure 4 is valid for the 27 following modes (from 28 to 54), so on and so forth (zones are identified on the curve by grey vertical dash lines). Integrated DOS curves give information about the localised or extended nature of the mode, as well as the localisation type (central or lateral) of transverse displacement and about the overloading distribution.

Here are all the elements of a periodic network coupled with a multi-scale distribution : this is the pre-fractal structure signature.

VIBROACOUSTIC COUPLING AND ACOUSTIC RADIATION

This part studies the modification of the acoustic response of a baffled structure, from the mechanical behaviour of the overloaded beam. Modes with a strong localisation of the transverse displacement presents an acoustic radiation attenuated in some frequencies regions.

Acoustic radiation modelling

As demonstrated by Elliott and Johnson, [START_REF] Elliott | Radiation modes and the active control of sound power[END_REF] total acoustic radiated power can be calculated either by mode amplitudes of the vibrating structure or wall celerities from a network of elementary radiating elements.

Without full analytical formulation, the second way seems more relevant. Far-field acoustic radiated power can be evaluated from Radiation Resistance Matrix (abbreviated RRM). Assuming a time-harmonic motion, the vibrating surface is divided into N rad rectangular elementary radiators, considered as monopoles vibrating at normal surface velocities of complex amplitudes described in column vector v e . If the wall sound pressures acting on each grid element are also assembled into a column vector p e , the far-field total acoustic sound power P ac (ω) radiated by a baffled structure can be formulated as the summation of the radiated powers of each element, leading to Eq. ( 12)

P ac (ω) = S 2 N rad {v H e (ω).p e (ω)} , (12)
where superscript H refers to the Hermitian transpose. Elementary radiator sizes have to be smaller than both structural and acoustic wavelengths (typically about the quart of them) and the structure is considered as baffled, under the classical hypotheses. [START_REF] Fahy | Sound and structural vibration -Radiation, transmission and response[END_REF][START_REF] Elliott | Radiation modes and the active control of sound power[END_REF] The sizeto-frequency limitation might be approached as high order modes getting closer, but the acoustic efficiency tending to unity physically means that the modal critical frequency has been reached (which can be proved by an analytical study).

Sound pressure vector can be described by the complex impedance matrix relation as

p e (ω) = [Z(ω)] v e (ω) . (13)
Then, the radiated resistance matrix [R] = S 2N rad

[Z] can be written as Eq. ( 14)

[R(ω)] = ω 2 S 2 e ρ 0 4πc 0                        1 sin(kr 12 )
kr 12

• • • • • • 1                        (14)
to finally recast Eq. ( 12) into

P ac (ω) = v H e (ω) [R(ω)] v e (ω) . (15)

Acoustic modal efficiency simulation

Results are expressed as acoustic radiation efficiencies. Indeed, it is a more fundamental concept since it is independent from the absolute size of the radiating element and from the characteristic impedance of the propagation medium. [START_REF] Wallace | Radiation resistance of a baffled beam[END_REF] It is defined as

σ(ω) = P ac (ω) P vib (ω) = P ac (ω) ρ 0 c 0 S < v 2 n (ω) > , (16)
where ρ 0 c 0 is the characteristic impedance of the air and P vib the vibrational reference power generated by a piston uniformly vibrating with < v 2 n (ω) > the normal averaged quadratic surface velocity.

Figure 6 presents a simulation of acoustic radiation efficiency function of frequency of mode 26, for a non-overloaded beam and a pre-fractal beam of order 3, as well as mode 8 of a non-overloaded beam of length L x /3. Here the results are presented in terms of modal acoustic efficiency, written σ m , because it is taken for a given mode m.

Analysis and discussion

Acoustic efficiency of a pre-fractal structure mode 26 is close to the mode 8 of a non-overloaded beam of length L x /3, as plotted on Figure 6. Indeed, by comparing both mode shapes from Figure 3, the central part of the over-loaded beam can be identified as the principal vibrating zone and thus predominant in the acoustic radiation, which corresponds to the mode 8 of a non-overloaded beam of length L x /3. It also corresponds to the central part of the Cantor set, of length a third of L x . However, the efficiency of the self-similar beam is still slightly smaller, the lateral parts also contribute to reducing the global acoustic radiation.

Moreover, the frequencies of non-overloaded beam mode 8 and order 3 pre-fractal beam are similar (difference of about 5%) and so on for their modal critical frequencies. That reinforces the idea that these modes present some similarities in their acoustic radiation.

Finally, by comparing self-similar and nonoverloaded beams of same length, the overloaded beam acoustic efficiency is slightly smaller in low frequencies. Efficiency then reaches the cross-over region and becomes larger until the modal critical frequency. The waviness of the non-overloaded beam is identified as side radiation lobes, [START_REF] Wallace | Radiation resistance of a baffled beam[END_REF] and it is even more tortuous as the mode is high-order. Yet, it is noticeable that the overloaded beam only presents a few, as the mode 8 of the beam of length L x /3, which consolidates the idea that the acoustic radiation of the pre-fractal beam is driven by the central part of mode 8, where the transverse displacement has a larger amplitude (localised part).

CONCLUSION

This paper has presented the vibroacoustic model of a sandwich beam overloaded by a pre-fractal distribution of masses. The sandwich structure is modelled as a homogenised material under simple bending dynamics, then discretised and approximated by a finite difference method.

Heterogeneities are integrated within the material such as additional masses distributed following a Canto-like set. The eigenvalue problem is then solved numerically to obtain the in vacuo modal basis.

Modal shapes have been presented for different orders. The chosen mode presents the localisation behaviour of the transverse displacement in the central part of the beam. The impact of the fractal order at iso local-overloading is then investigated at order 2 and 3, and a large localisation increase is identified by comparing maximal amplitudes of central and lateral parts.

Integrated DOS enables to feature localised modes through slope breakdowns as well as the type of localisation (central or lateral). The modes for which the over-masses are situated on vibration nodes are also identified because frequencies of non-overloaded and pre-fractal beams match as their mode shapes do.

Afterwards, fluid-structure coupling has been modelled through elementary radiators and the acoustic efficiency extracted from radiated acoustic power is calculated using a radiation resistance matrix. The acoustic radiation efficiency of the overloaded beam is close to the one of a one-third length beam (global behaviour, waviness and modal critical frequency), while remaining slightly smaller. The pre-fractal beam and its localised modes appear to present a smaller acoustic radiation compared to a non-overloaded beam in a given frequency range.

The interest of using such self-similar construction process is intrinsically linked to its definition. Indeed, the heterogeneity disposition goes further than a simple periodic distribution, by adding a homothetic ratio linked to the scheme order n frac . This adds the repetition of slope breakdowns in the integrated DOS, as illustrated on Figure 5. Furthermore, the over-mass distribution presents a real physical sense, contrary to current reduction techniques, as viscoelastic patches for instance.

Finally, the masses are added within the honeycomb core, so within the structure, contrary to viscoelastic treatments that are glued to the external skins. Therefore, they have a direct action on the mechanical waves in the material, acting on the wave propagation with phenomenons induced by the heterogeneities.

Figure 1 :

 1 Figure 1: Original Cantor set, orders 0 to 3.

Figure 2 :

 2 Figure 2: Self-similar loaded beam sketch, orders 0 to 3.

Figure 3 :

 3 Figure 3: Transverse deflection (mode 26) : (-) classical beam (analytical), (-•-) pre-fractal beam of order 2 (FD) and (-o -) pre-fractal beam of order 3 (FD).

Figure 4

 4 Figure 4 is first analysed because plotted in logarithmic scale, which emphasises the phenomena representation. It presents the first set of extended and localised modes. An extended mode is a classical modal mode whose maxima stand at the same amplitude. Several areas of the DOS behaviours are identified, directly linked to the self-similar overloaded distribution:

Figure 4 :

 4 Figure 4: Natural frequencies of first 30 modes : (--) classical beam, (-•-) pre-fractal beam (order 3), and (-) homogeneous beam of equivalent mass of pre-fractal beam (log scales).

Figure 5 :

 5 Figure 5: Natural frequencies of the first 80 modes : (--) classical beam, (-•-) pre-fractal beam (order 3), in linear scale.

Figure 6 :

 6 Figure 6: Acoustic radiation efficiency : mode 26 (-) classical beam, mode 26 (•) pre-fractal beam (order 3) and mode 8 (--) of a non-overloaded beam of length L x /3.

ACKNOWLEDGEMENTS

This research has been supported jointly by the Midi-Pyrénées region and ONERA (The French Aerospace Lab).


Sara Bouraine
email: s_bouraine@yahoo.fr

Thierry Fraichard
email: thierry.fraichard@inria.fr

Ouahiba Azouaoui
email: oazouaoui@cdta.dz

Real-time Safe Path Planning for Robot Navigation in Unknown Dynamic Environments

Keywords: mobile robotics, motion planning, motion safety, dynamic environments

This paper solves a motion planning problem from a motion safety perspective, where a variant of the classical Rapidly exploring Random Tree (RRT) approach [1] called psafe RRT is proposed. The exploration of the search space is similar to RRT, however, the highlight of p-safe RRT is the integration of passive motion safety. The basic principle of this safety level is to guarantee that the system can brake down and stop before collision. P-safe RRT extends a tree through the state time space, where tree's nodes and primitives are checked for passive motion safety. The computed trajectory is passively safe and drives the robot from its initial state to the goal state. The developed algorithms have been tested in simulation scenarios; featuring both fixed and moving objects with unknown trajectories for a car-like robot with a limited field of view.

INTRODUCTION

Few years ago, Autonomous Ground Vehicles (AGVs) were a challenge, but, today it becomes a reality: many car dealers currently work on AGVs projects (e.g. BMW, Toyota, Renault, Volkswagen, Peugeot, General Motors, Mercedes-Benz, or Tesla). Several demonstrations have been carried out on real roads, i.e. in the presence of other cars and pedestrians (both with unknown behaviour), having only a partial knowledge about the environment. However, the risk of accident remains present (see [START_REF] Fletcher | The MIT-Cornell collision and why it happened[END_REF]). Even if these systems are autonomous, up to this day, cars are never unmanned [START_REF] Thrun | What we're driving[END_REF]; will the driver seat ever be empty? Therefore, guaranteeing motion safety remains an open problem in such situations.

There is a rich literature on motion planning, most of approaches compute a compete path from an initial position of the robot to the goal position based on the environment representation (which is usually a priori known or built during a first phase of exploration) [START_REF] Latombe | Robot Motion Planning[END_REF][START_REF] Choset | Principles of Robot Motion: Theory, Algorithms, and Implementation[END_REF][START_REF] Lavalle | Planning Algorithms[END_REF][START_REF] Laumond | Robot Motion Planning and Control[END_REF][START_REF] Siegwart | Introduction to Autonomous Mobile Robots[END_REF][START_REF] Ge | Autonomous Mobile Robots: Sensing, Control, Decision Making and Applications[END_REF]. Generally, the environment is considered as static. These methods can be adapted to dynamic environments but adding the temporal dimension increases the complexity of the problem. Other methods as PRM (Probabilistic Roadmap) [START_REF] Kavraki | Probabilistic roadmaps for path planning in high-dimensional configuration spaces[END_REF] and RRT [START_REF] Lavalle | Randomized kinodynamic planning[END_REF][START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planning[END_REF] have more performances in high dimensional configuration spaces. However, PRM requires a priori knowledge about the environment or an initial map constructed off-line. Otherwise, RRT is more suitable for unknown environments. Besides, the kinodynamic constraints of the system are explicit-ly taken into account and the incremental tree construction features better the strong changes of the environment compared to PRM. Many variants of RRT have been developed over years. First works [START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planning[END_REF][START_REF] Lavalle | Randomized kinodynamic planning[END_REF] were very expensive in computing time as the tree is expanded over the entire workspace. More recent works have combined RRT with other methods like PRM to solve the problem [START_REF] Bekris | Multiple query motion planning using single query primitives[END_REF][START_REF] Akinc | Probabilistic Roadmaps of Trees for Parallel Computation of Multiple Query Roadmaps[END_REF]. All these works are applied in static environments. Other extensions of RRT have been also proposed to drive the robot in dynamic environments; in [START_REF] Bruce | Real-time randomized path planning for robot navigation[END_REF], RRT has been extended to ERRT (Extended RRT), where environment's changes are featured by interleaving planning and execution phases. An anytime RRT version was proposed in [START_REF] Ferguson | Anytime RRTs[END_REF]; at each time instance, a new plan is computed based on the environment information update. Another variant of RRT is Multipartite RRT, which combines ERRT and anytime RRT [START_REF] Zucker | Multipartite RRTs for Rapid Replanning in Dynamic Environments[END_REF]. It is suitable for environments containing moving objects with a priori known behaviour. There are other extensions of RRT (ex. [START_REF] Macek | Motion Planning for Car-Like Vehicles In Dynamic Urban Scenarios[END_REF][START_REF] Heon-Cheol | Grafting: A Path Replanning Technique for Rapidly-Exploring Random Trees in Dynamic Environments[END_REF][START_REF] Ma | Efficient Sampling-Based Motion Planning for On-Road Autonomous Driving[END_REF]), but, in most cases, the future behaviour of moving objects is not considered at all or assumed a priori known.

Generally, guaranteeing motion safety requires a priori knowledge about the overall trajectory of each moving object to ensure that a collision will never occur: that is what is called absolute motion safety. However, when the future behaviour of objects is unknown, this form of safety is impossible to guarantee [START_REF] Fraichard | Will the driver seat ever be empty?[END_REF]. Therefore, following the opinion of some authors [START_REF] Macek | Towards safe vehicle navigation in dynamic urban scenarios[END_REF], it is better to settle for weaker levels of motion safety: it is better to guarantee less than to guarantee nothing. In [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF], the proposed solution is to guarantee that the robot will be at rest if a collision took place. This is what we called passive motion safety.

The contribution of this paper is a variant of the classical RRT method [START_REF] Lavalle | Randomized kinodynamic planning[END_REF]; a new approach p-safe RRT is proposed. It is based on building an incremental tree in the state time space. The generated nodes and trajectory primitives are checked for passive motion safety based on a verification algorithm of passively safe states that is proposed in [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF]. So, the planner computes a passively safe trajectory to the goal.

The paper is organized as follows. Section 2 explains the paradigm collision-free vs. passive motion safety. In Section 3, the classical RRT principle is first presented, then the proposed extension p-safe RRT is described in detail. Finally, the validation of the developed algorithms is illustrated in Section 4 thanks to a simulation implementation.

COLLISION-FREE VS. SAFETY OF MOTION

From a motion safety perspective, previous works (e.g. [START_REF] Pallottino | Decentralized cooperative policy for conflict resolution in multivehicle systems[END_REF][24][START_REF] Molinos | Dynamic obstacle avoidance based on curvature arcs[END_REF][START_REF] Lalish | Decentralized reactive collision avoidance for multivehicle systems[END_REF][START_REF] Bekris | Safe and distributed kinodynamic replanning for vehicular networks[END_REF][START_REF] Gopalakrishnan | Time scaled collision cone based trajectory optimization approach for reactive planning in dynamic environments[END_REF][START_REF] Blanco | Extending obstacle avoidance methods through multiple parameter-space transformations[END_REF]) proposed approaches that guarantee collision-free motion. However, it is not always sufficient; collision occurs regardless the undertaken action. For example, in the case of a vehicle traveling at high speed, confronted to a wall, it can be in a situation where collision is inevitable given its dynamic constraints (it is not possible to avoid the wall or to stop before collision). In other words, the vehicle is placed in an inevitable collision state, i.e. whatever the future trajectory of the vehicle, collision occurs. Therefore, to guarantee motion safety, inevitable collision states should be avoided.

Given the perceptual limitations of the robot and the unknown future behaviour of the moving objects present in the environment, it is not possible to guarantee absolute motion safety (collision will never occur given a priori known model of the future). This is why, we settled for weaker level of motion safety, but stronger given the harsh constraints related to the robot and the environment. This level of safety is called passive motion safety (p-safety). It is passive in a sense that the robot takes its own responsibility to not to be harmful with respect to its surrounding environment, regardless if the other obstacles will collide with him. Passive safety guarantees that if a collision takes place, the robot will be at rest. From this strategy, the braking inevitable collision states concept (braking ICS) was proposed [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF]. A braking ICS (denoted ICS b ) is a state for which whatever the future trajectory of the robot, a collision occurs before the robot is stopped. A state that is ICS b -free is a p-safe state.

To verify if a state is ICS b or not (by duality p-safe or not), a braking ICS checker (called ICS b -CHECK) was developed in [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF]. The developed planner presented in section 3 is based on ICS b -CHECK to check the p-safety of the planned trajectory.

THE DEVELOPED MOTION PLANNING TECHNIQUE

This paper proposes a new motion planning approach called p-safe RRT. It is inspired from RRT method [START_REF] Lavalle | Randomized kinodynamic planning[END_REF]. However, unlike RRT, p-safe RRT has to guarantee passive motion safety criterion and takes into account the field-of-view limits and occlusions besides to the future behaviour of moving objects.

The classical RRT

RRT is a diffusion technique based on probabilistic sampling of the search space [START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planning[END_REF][START_REF] Lavalle | Randomized kinodynamic planning[END_REF]. It has been used to solve motion planning problems in high dimensional configuration spaces. A tree is incrementally expended through the space (set of nodes related by primitives) for planning a path from an initial state to a goal state . The basic principle of this technique is explained in the algorithm presented in figure 1. First, the tree (noted TREE) is initialized by the initial state of the robot. It represents the root node . The function TREE_EXPANSION grows the tree from toward randomly selected nodes (using the function RANDOM_TARGET): when a random target is generated, the function NEAREST_NEIGHBOR selects the nearest node in the tree relatively to . is then extended toward by applying a control that is generally selected from a set of possible controls (thanks to NEAREST_CTRL). This selection is based on a minimized Euclidean distance between and . The resulting node is ; it is added to the tree with its corresponding new primitive . If the expansion of the tree is obstructed by an obstacle, is not generated. This process is repeated until the state is reached (i.e. when the distance between and is lower than a certain threshold).

Generally, RRT and its variants [START_REF] Lavalle | Randomized kinodynamic planning[END_REF][START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planning[END_REF][START_REF] Bekris | Multiple query motion planning using single query primitives[END_REF][START_REF] Akinc | Probabilistic Roadmaps of Trees for Parallel Computation of Multiple Query Roadmaps[END_REF] solve motion planning problem only in the configuration space or the state space, the time parameter is not considered at all. Furthermore, most of these works are applicable in static environments. Even if RRT has been extended for dynamic environments [START_REF] Bruce | Real-time randomized path planning for robot navigation[END_REF][START_REF] Ferguson | Anytime RRTs[END_REF][START_REF] Zucker | Multipartite RRTs for Rapid Replanning in Dynamic Environments[END_REF][START_REF] Macek | Motion Planning for Car-Like Vehicles In Dynamic Urban Scenarios[END_REF][START_REF] Heon-Cheol | Grafting: A Path Replanning Technique for Rapidly-Exploring Random Trees in Dynamic Environments[END_REF][START_REF] Ma | Efficient Sampling-Based Motion Planning for On-Road Autonomous Driving[END_REF], as the future behaviour of moving objects is not considered or a priori known, motion safety cannot be guaranteed.

To solve such issue, we propose a modified version of RRT; p-safe RRT. It extends a tree in the state time space and computes a p-safe trajectory that drives the robot to the goal in a safe manner.

P-safe RRT

The developed technique p-safe RRT constructs a tree thanks to an exhaustive exploration of the state time space, i.e. using a fixed set of feasible controls . The expanded tree is checked for passive motion safety using the algorithm ICS b -CHECK proposed in [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF] based on a model of the future (this point is detailed in [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF]). Besides ICS b -CHECK, an important property should be established concerning the psafety guarantee of a trajectory, it is expressed as follow: Property 1: A trajectory is p-safe (i.e. not a braking ICS) if P is collision-free and the final state of P is p-safe.

For more details concerning the proof of this property, or other p-safety guarantee properties and definitions, the reader is referred to [START_REF] Bouraine | Passively Safe Partial Motion Planning for Mobile Robots with Limited Field-of-Views in Unknown Dynamic Environment[END_REF].

The p-safe RRT principle is explained through the algorithm presented in figure 2. As for RRT, p-safe RRT is initialized by the root node . In addition, the algorithm has as inputs and . During the planning time, each node in the tree is expanded thanks to the function PSAFE_TREE_EXPANSION according to a set of controls . Each control leads to the generation of a new primitive and a new node (i.e. the final state of ) (via GENERATE_TRAJ_PRIM). is checked for p-safety based on ICS b -CHECK (BRAKING_ICS_CHECK) and property 1. In the case no p-safe primitive is found for a given node , a braking trajectory is generated to guarantee the motion safety (using GENERATE_BRAKING_TRAJ). At the end of the process, different p-safe trajectories are possible, but only one trajectory is selected. The function SELECT_BEST_TRAJ is responsible of this task using an optimization function which is based on the Euclidean distance (∆ ) between the trajectory final state and the goal state and the time cost of the trajectory (∆ ). Weighting factors and are respectively associated to ∆ and ∆

. is expressed as follows: This algorithm illustrates the general principle of p-safe RRT technique, however, for motion safety reasons, the planning process can be interleaved with the execution process like in [START_REF] Bruce | Real-time randomized path planning for robot navigation[END_REF] or instead a partial motion planning can be applied like in [START_REF] Bouraine | Passively Safe Partial Motion Planning for Mobile Robots with Limited Field-of-Views in Unknown Dynamic Environment[END_REF][START_REF] Bouraine | contrıbutıon a la planıfıcatıon de mouvements en envıronnements dynamıques pour des robots mobıles de type voıture : cas du robucar[END_REF], where a new plan is periodically computed until the robot reaches its goal.

= ∆ + ∆ (1) PSAFE_RRT_CONSTRUCTION

RESULTS

In order to illustrate p-safe RRT performances and demonstrate its motion safety guarantee, p-safe RRT algorithm has been implemented in simulation scenarios that feature fixed and moving objects with arbitrary trajectories. The application has been developed in C++ and it has been tested on a laptop equipped with an Intel Core i7 (1.6GHz, 4GB RAM, OS: Linux).

4.1

The system's dynamic model P-safe RRT has been tested for a car-like robot that operates in a state time space. The robot state is a 5 tuple = ( , , , , ). ( , ) are the Cartesian coordinates (related to the midpoint of the rear axle of the robot), is the robot's orientation, its linear velocity and its steering angle. The dynamic of the robot is governed by the following differential equations:

= / 0 0 + 0 0 0 1 0 + 0 0 0 0 1 (2)
Where denotes the robot's wheelbase, is the linear acceleration and is the steering angle velocity. These last two parameters represent the control inputs that should be applied to the system. The robot behaviour is limited by kinodynamic constraints that should be respected, namely:

| | ≤ , | | ≤ , | | ≤ , ≤
, where = 20 / , = 0.314 , = 7 / , = 0.314 / .

P-safe RRT

To illustrate how the developed technique works, it has been tested in a scenario called simple urban scenario (see figure 3.a); it contains fixed and moving objects with an arbitrary behaviour. In this case, no information concerning the future behaviour of the moving objects is available, only the maximum velocity is considered ( = 20 / ). The robot has a limited field-of-view with maximum radius of 80 and it moves in a 2D workspace with 180 × 180 dimensions. The radius of the disk objects is 2.5 . P-safe RRT has to drive the robot starting from initial position until a predefined goal while considering the passive motion safety guarantee. Both seen and unseen objects (perception limits and occlusions) are considered. For safety purpose, p-safe RRT computes at first ICS b set (see figure 3.b), the black region represents the forbidden part of the environment that should be avoided by the robot.

The tree is expanded in the space following a set of controls ( , ) selected as follows: = , i.e. a constant maximum linear acceleration is considered and ∈ [-, ], i.e. a constant steering angle velocity is selected in this interval of values. [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF] for more details). The robot's field-of-view is represented by the magenta segments (perception limits) and the green segments (occlusions).

Figure 4 illustrates the behaviour of the robot further to the implementation of psafe RRT in the simple urban scenario. The executed trajectory is represented by the thick trace behind the system while the trajectory in front of the robot is the planned part. Based on the ICS b set computed given an updated model of the future, the robot moves until it reaches the goal while avoiding seen and unseen objects (the limits of the field of view). For example, in the second snapshot, the robot escaped an ICS b region resulting from perception limits (unexpected objects). In the third snapshot, it escaped a fixed object and an unexpected objects region. The built p-safe RRT trees corresponding to these two snapshots are respectively illustrated in figure 5.a and figure 5.b. Each tree's primitive and node is checked for passive safety. Only p-safe nodes are expandable. For an ICS b node, a collision-free braking trajectory is generated. This experiment shows the performance of p-safe RRT to guarantee a passively safe behaviour for the robot. Even the avoidance is not possible, the robot brakes down to remain in a p-safe state. ry is generated). The selected trajectory (to be executed) is represented in red.

CONCLUSION

AGVs navigation in real road is still challenging up to now. In this paper, we proposed a solution by developing a passively safe motion planner; p-safe RRT. It is a variant of the classical RRT technique that is based on the guarantee of a passively safe motion. This navigation system takes into account the dynamics of both the robot and the environment, where an unknown model of the future is considered. P-safe RRT has been validated through simulation tests; the robot showed a passively safe behaviour by avoiding perceived objects and unexpected objects that are potentially dangerous from a motion safety point of view.

A logical future work to this one is the implementation of this planning system on an experimental platform. On the other hand, concerning the approach concept, the tree expansion strategy used in this paper is basic. The exploration of the space time could be improved by means of other existing expansion techniques or strategies.

Fig. 1 .

 1 Fig. 1. The classical RRT algorithm.

Fig. 2 .

 2 Fig. 2. P-safe RRT algorithm.

Fig. 3 .

 3 Fig. 3. (a) Simple urban scenario; polygons are the fixed objects and black discs are the moving objects with their corresponding trajectories represented in blue. The robot is the red disc at the center. (b) ICS b set for a given position of the robot; the forbidden states (ICS b ) are represented in black while the p-safe states are represented in white (see[START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF] for more details). The robot's field-of-view is represented by the magenta segments (perception limits) and the green segments (occlusions).

Fig. 4 .

 4 Fig. 4. Snapshots of p-safe RRT at work in the simple urban scenario (the goal to reach is represented by the red point) (see text).

Fig. 5 .

 5 Fig. 5. P-safe RRT trees corresponding to the second (a) and third snapshot (b) of figure 4. Tree's primitives represented in blue are p-safe, those in cyan are ICS b (not p-safe) and those in magenta characterize collision-free braking trajectories (for each ICS b node, a braking trajectory is generated). The selected trajectory (to be executed) is represented in red.

(a) (b)