P Rahmé

C Bouvet

S Rivallant

V Fascio

G Valembois

Experimental investigation of impact on composite laminates with protective layers

Keywords: B, Impact behaviour C, Damage tolerance B, Delamination

This paper presents an experimental study of low energy impacts on composite plates covered with a protective layer. In service, composite materials are subjected to low energy impacts. Such impacts can generate damage in the material that results in significant reduction in material strength. In order to reduce the damage severity, one solution is to add a mechanical protection on composite structures. The protection layer is made up of a low density energy absorbent material (hollow spheres) of a certain thickness and a thin layer of composite laminate (Kevlar). Energy absorption ability of these protective layers can be deduced from the load/displacement impact curves. First, two configurations of protection are tested on an aluminium plate in order to identify their performance against impact, then the same are tested on composite plates. Test results from force-displacement curves and C-scan control are compared and discussed and finally a comparison of impact on composite plates with and without protection is made for different configurations.

Introduction

Composite materials are widely used nowadays in aeronautics. This growing interest is due to the relatively high strength/mass ratio of these materials compared to those of metals. These materials are sometimes subjected to low energy impacts during fabrication and in service as well which can have an influence on the residual mechanical properties of the structure [START_REF] Abrate | Impact on composite structures[END_REF]. In the field of aeronautics, depending on the impact detectability, damage tolerance is often the key factor in structure design [START_REF] Coutellier | Presentation of a methodology for delamination detection within laminated structures[END_REF][START_REF] Alderliesten | Damage tolerance of bonded aircraft structures[END_REF]. It is a known fact that damage due to impact reduces the residual compressive strength of a structure to less than 50% of the initial strength. The notion of impact detectability was introduced for aircraft certification needs. The minimum indentation damage that can be detected by visual evaluation is called ''Barely Visible Impact Damage'' (BVID) [START_REF] Moody | Width effects on the compression strength of composite sandwich test specimens after barely visible impact damage[END_REF], and in aeronautical standards, this threshold of detectability after few days of rest and humidity ageing is 0.3 mm depth. In the case of a non-detectable damage (indentation less than BVID), the structure must sustain the ultimate load. Consequently, damage tolerance in composites has been a subject of investigation for many years. Several authors, [START_REF] Mitrevski | The influence of impactor shape on damage to composite laminates[END_REF][START_REF] Li | Low-velocity impact-induced damage of continuous fiber-reinforced composite laminates Part II. Verification and numerical investigation[END_REF][START_REF] Davies | Impact damage prediction in carbon composite structures[END_REF][START_REF] Ballère | Damage tolerance of impacted curved panels[END_REF][START_REF] Abi Abdallah | Experimental analysis of damage creation and permanent indentation on highly oriented plates[END_REF][START_REF] Bouvet | Low velocity impact modelling in laminate composite panels with discrete interface elements[END_REF] have studied the resistance of composite structures against low velocity impacts.

The use of mechanical protection against impact is one other way of solving the problem of impact on composites. Although introduction of a mechanical protection against impact induces a non-negligible increase in weight of the structure, it nonetheless avoids or considerably reduces damage in the structure. Furthermore, use of mechanical protection improves the detectability of impact, as the protection is generally more deformable which causes a visible post-impact indentation mark.

Core material covered with skin seems to be a good candidate for protective layers. Core materials (honeycomb, foam, hollow spheres ...) are of significant scientific interest due to their good specific resistance and enhanced energy absorption capabilities and are already used to protect structures against impact. For example, some aircraft cockpits are protected against bird strikes by a layer of aluminium honeycomb covered by an aluminium skin to improve the capacity of energy absorption. Wang [START_REF] Wang | Impact behavior and energy absorption of paper honeycomb sandwich panels[END_REF] showed that the energy absorption tends to increase linearly with the increase in relative density of honeycomb cores. Therefore, increasing the relative density of honeycomb cores can efficiently improve the dynamic cushioning properties of the sandwich panels. Li et al. [START_REF] Li | Energy-absorption performance of porous materials in sandwich composites under hypervelocity impact loading[END_REF] used different materials in their study on high velocity impacts in which they determined and compared the energy absorption capability of these materials. They concluded that porous silicon carbide materials are best suited for absorbing impact energy. On the other hand, Apetre et al. [START_REF] Apetre | Low-velocity impact response of sandwich beams with functionally graded core[END_REF] studied the impact at low velocity on honeycomb sandwich beams with variable density and rigidity in the thickness. They showed that core with graded properties reduces maximum deformation corresponding to maximum impact load. Furthermore, Shin et al. [START_REF] Shin | An experimental study of low-velocity impact responses of sandwich panels for Korean low floor bus[END_REF] presented impact tests at low velocity on several types of sandwiches. They concluded that glass skin sandwich plates have a better resistance against impact compared to aluminium skins. In the same context, Petit et al. [START_REF] Petit | Impact and compression after impact experimental study of a composite laminate with a cork thermal shield[END_REF] showed that a layer of thermal protection (cork) is also a good mechanical protection against impact. The use of such protections increases the residual strength of compression after impact tests. In the same context, Rahmé et al. [START_REF] Rahmé | Mechanical protection for composite structures submitted to low energy impact[END_REF] compared the resistance against impact of three type of cores (foam, honeycomb and hollow spheres) having the same thickness or the same surface density. They concluded that hollow spheres give best efficiency against impact. Moreover, the concept of protective layers, consisting of epoxy filled glass microballoons covered with one or more layers of aluminium gauze was presented by Hart and Ubels [START_REF] Hart 't Wgj | Impact energy absorbing surface layers for protection of composite aircraft structures[END_REF]. Their work also showed that the resistance to impact damage of a composite structure can be improved by application of a protecting surface layer.

Anyways, the main objective is to obtain a good compromise between the residual strength and the mass of the structure, for a given impact energy (or a given BVID).

In this paper, two configurations of protective layers made of hollow spheres core (HS) and Kevlar skin are tested. Impact tests at low velocity (<10 m/s) and low energy (<100 J) are performed using a 20 mm diameter impactor in a drop weight setup. In the first step, composite plates are impacted without protection, then both configurations of layers are tested on aluminium and composite plates. A comparative study is then carried out based on the results of different configurations. Finally, the results of impact on composite plates with and without protection are compared and discussed. Furthermore, the observations of C-scan controls are also discussed.

Experimental study

Drop weight device

There are many testing procedures to simulate an impact on a structure; however, the drop weight remains the most used device [START_REF] Abrate | Impact on composite structures[END_REF]. Such a device has been used in this study to perform impact tests according to AITM 1-0010 standards [START_REF]Airbus Industrie Test Method: AITM 1-0010: Determination of compression strength after impact[END_REF]. The principle of this drop weight device is to drop a mass guided in a tube on a composite plate. Fig. 1a andb shows impacts on aluminium and composite plates respectively having protective layers. This device is dedicated to impact tests at low velocity (<10 m/ s) and low energy (<90 J). A 4.02 kg impactor with a 20 mm diameter hemispherical head is used. The impactor is instrumented with a load sensor, installed between the impactor head and its body. This KISTLER piezoelectric sensor having a maximum capacity of 120 kN is calibrated to measure the impact load. An optical sensor (Laser diode, Fig. 1) measures the initial velocity of the impactor before impact. The impactor displacement and velocity histories are obtained by integrating the impactor acceleration resulting from the measured load taking into account the initial velocity [START_REF] Petit | Impact and compression after impact experimental study of a composite laminate with a cork thermal shield[END_REF]. When impacting on composite plate, an additional optical sensor (OptoNCDT50) is used to measure the displacement in the impact direction of the non-impacted face of the composite plate. This displacement due to the bending of the plate helps in determination of indentation of the layers and then the corresponding energy absorption. This allows comparison of the different impact curves when supporting on aluminium and composite plates.

Protective layers and plates

Two configurations of protective layers have been tested. Configuration L1 is designed for 50 J energy impacts (cf. Fig. 2). It is composed of a 1.4 mm thick 0°/90°Kevlar woven fabric skin and of a polymer hollow spheres core made by ATECA Company. Sphere diameter is between 5.4 and 6 mm, and spheres are glued together and also with the skin. In the context of practical application on real world structures, for a simpler installation and an easier maintenance after impact, protective layers will not be glued to the structure. Consequently, an additionally thin 0°/90°ply of Kevlar woven fabric has been added between the hollow spheres and the plates. The patent ''Peau amortissante de protection de pièces composites'' no. 2 930 478 has been filed about these protective layers of hollow spheres [START_REF]Peau amortissante de protection de pièces composites[END_REF]. The mechanical properties of the skin and the hollow spheres cores used are provided in Tables 1 and2 respectively. For hollow spheres, the properties given in the table are homogenised characteristics, calculated from spheres layers compression tests. The Young's Modulus is the initial slope of the stress-strain curve, whereas the failure stress is the stress reached at the crushing plateau.

Configuration L2 is developed in order to withstand 90 J energy impacts. It consists of a 2.1 mm thick 0°/90°Kevlar woven fabric skin and two layers of hollow spheres core with an additionally thin 0°/90°ply as shown in Fig. 2. The design of this configuration is based on gradual stop of the impactor.

Table 3 gives the global surface density for each configuration of the tested protective layer along with material and geometrical information. In the case of tests on aluminium plates, the plate is 2 ] s according to AITM 1-0010 standards [START_REF]Airbus Industrie Test Method: AITM 1-0010: Determination of compression strength after impact[END_REF]. The thickness of one ply is 0.25 mm which amounts to a total plate thickness of 4 mm. The in-plane dimensions of the composite plate are 150 Â 100 mm 2 . During tests, composite plates are clamped using the device shown in Fig. 1b. A vertical load is applied on a clamping frame to prevent rebounds of the plate.

Test results

The specimens were impacted at several impact energies. Configuration 1 specimens were tested at 15, 30 and 50 J impact energies and composite plates without protection were also impacted at the same energy (15, 30 and 50 J). Impact energies of 15, 50 and 90 J were used for the thick protection layers (configuration L2). Displacement versus time curve is given in Fig. 3, for A-L1-30 J impact test (aluminium support and protective layer L1). Furthermore, load/displacement impact curves are presented in Fig. 4, for 30 J energy impact on configuration L1 specimen with aluminium support and composite plate. Tests on the other protective layers present similar curves. Both unfiltered and filtered impact forces are presented in Fig. 4a (15 kHz low pass filter), but onwards only filtered impact forces will be presented in his article. The initially measured impact velocity in this test is 3.9 m/s. The maximum displacement of the impactor is then about 6.2 mm in the case of aluminium support. This displacement is lower as compared to the 7.5 mm layer thickness represented on the curve by the vertical line. Hence, the impactor will not touch the aluminium plate in this case and the layer should be able to absorb the total impact energy.

The displacement of the non-impacted face of the composite plate is due to its bending. Using the LASER sensor, this displacement was measured. The indentation of the protective layer is then deduced by subtracting the plate bending displacement from the impactor displacement (Fig. 3), with the assumption that the thickness of the composite plate does not change (no local indentation of the plate). In the case of aluminium support, the layer indentation and the impactor displacement are the same. The maximum indentation of the impactor is around 5.4 mm, which appears in the protective layer when supported on composite plate.

The various steps of the impactor displacement through the protective layer can be identified during impact (Fig. 4a). First, there seems a very short elastic behaviour of the protective layer which is then followed by a progressive damage of both the skin and the hollow spheres. A major failure of the skin is then reached, characterised by a significant drop in load. After this break, the skin damage continues and the hollow spheres densification is held up to the maximum load. The indentation curve in Fig. 4b shows the same behaviour of the layer on composite plate. The maximum indentation is smaller than in previous case, due to the transformation of initial energy in damage and in elastic bending in the composite plate.

Fig. 5 shows the non-impacted faces of composite plates impacted at 50 J, with protection L1 (a), and without protection (b). No damage is visible on the first one (even if C-Scan shows internal delaminations: see further), whereas obvious longitudinal failure of the last ply is visible for the non-protected plate.

On the other hand, Fig. 6 shows the different protective layers after impact at maximum impact energy (50 J for L1 layer, 90 J Impact Force (kN) . Displacement (mm) .

Comp or Al plate

Hollow spheres

Hollow spheres

Comp or

Time (s)

Impactor LASER Indentation Impact Force Fig. 3. Force and displacements versus time curves for A-L1-30 J impact test. for L2). Impacts lead to hemispherical indentation marks with cross shaped failure of the skin. As shown on impact curves, the indentation is smaller in the case of impact on composite. The results observed on the different layers are compared and discussed in the next section.

Comparison and discussion

First, the different impact curves corresponding to protected and non-protected composite plates are presented in Fig. 7. Impact force is drawn as a function of the impactor displacement. For each protection configuration, the superposition of 15 and 50 J impact energy curves indicates good repeatability of the impact tests.

By carefully examining these curves, comparison of the behaviour of composite plates with or without protection can be carried out. The presence of protection reduces the impact force during the first 7 mm of displacement by mechanical fuse effect. Later, when the core layer is totally crushed, the load increases. Paradoxically, the load increases more in layer L1 than in the non-protected plate. The reason being that damage in non-protected plate reduces highly the rigidity under the impactor which then reduces the force. In case of protected plates, force spreading due to the protective layer leads to less critical damage due to which there remains a margin in the impact force which can still increase.

Concerning layer L2, at 50 J, the maximum force is not as high as in layer L1 test because the core is not entirely crushed under the impactor.

Fig. 8 compares the impact curves on protected aluminium and composite plates, for different impact energies. Fig. 8a shows the results for configuration L1 at 15 J and 50 J. The impact force is drawn as a function of indentation, that is crushing displacement in the core. In the first phase of the impact, the different curves show almost identical behaviour. Despite bending of the composite plate, the impact force is ruled by the indentation of the protective layer.

However, after 4 mm indentation, the behaviour of the layers becomes different on aluminium and on composite plates during the hollow spheres densification phase. Due to appearance of damage in composites, the force is lower in composite plate tests. For a given energy, the maximum indentation is also lower in the case of composite plate, due to absorption of energy due to damage in composites and release of bending energy. Similarly, Fig. 8b shows the curves of configuration 2 i.e., impact on aluminium and composite plate at different energies (15, 50 and 90 J). The curves are also similar except there is a difference in slopes due to densification of the hollow spheres. Note that the total impact energy is almost absorbed by the protective layer in case of the aluminium support (there is almost no impactor rebound), whereas, in composite plates a part of that energy is transformed into damage.

Table 4 presents the distribution of energies during a test. Impact energy which is the initial energy of the impactor is made up of the following:

-Energy dissipated in the protective layer: calculated from impact force and indentation.

-Elastic recovery energy (impactor rebound energy): calculated from force release. -Energy absorbed in the composite plate by damage: calculated from the subtraction of the two previous energies from impact energy. -Kinetic energy of the plate and energy dissipated in the test set-up can be neglected in this study.

Up to 50 J for L1 configuration, and 90 J for L2, the energy absorbed by the protective layer reaches about 50% or more, which proves the relevance of the use of such a protection.

To identify the importance of the protective layer C-scan observations are made for all tested composite plates. These C-scan observations help in determining the shape of the different delaminations and their locations in the plates. Fig. 9 presents the different C-scan snapshots at different energies (impacted face). It is clearly visible from these images that while impacting non-protected composite plate, damage is already present at 15 J impact energy. This damage increases with the rise in impact energy. For protected plates with the layers of configuration L1, no defect is seen for 15 and 30 J impacts. The defects appear at 50 J impact energy while nothing was visible with the naked eye on the external surfaces. Similarly, there are no defects on configuration L2 plates until 50 J. Fig. 10 presents the projected delaminated area of the impacted composite plates as a function of impact energy (a) and the energy absorbed by the composite plate (b). The first curve shows clearly the efficiency of the use of protective layers. The second one also shows that taking into account the only absorbed energy in the plate, behaviour is not the same between protected or non-protected plate.

A detail of delamination surface for composite plate with (L2) and without protection subjected to 30 and 90 J impacts respectively is presented in Fig. 11. The delaminated area in this case is respectively 1972 and 2031 mm 2 and the locations of different delaminations are also indicated in this figure. The differences between the two C-scan shots are clearly visible. In the case of non-protected plates, the damage is symmetrical about the impact point, and the classical delamination form of double helix is obtained. Moreover, every interface delamination is oriented by the fibres direction of the ply located below.

In the case of protected plate, a very asymmetric damage shape is observed on all tested specimen, and the location of delamination in the thickness of the laminate is limited to the central interfaces [START_REF] Alderliesten | Damage tolerance of bonded aircraft structures[END_REF][START_REF] Moody | Width effects on the compression strength of composite sandwich test specimens after barely visible impact damage[END_REF][START_REF] Mitrevski | The influence of impactor shape on damage to composite laminates[END_REF], like illustrated in Fig. 12. The first two interfaces of non-impacted face of the protected plates are largely protected, whereas usually it is the most damaged interface. Such observations have already been made in previous works of Petit et al. [START_REF] Petit | Impact and compression after impact experimental study of a composite laminate with a cork thermal shield[END_REF]. The physical explanation of the reason behind this specific damage shape is difficult. Some numerical investigations are in progress in order to understand the chronology of damage that leads to this phenomenon. About the impacted side, by taking into account only indentation in the plate not in the protective layer, it is observed that the permanent indentation is much smaller in the case of the protected plate.

As a result of this study, the proposed protective layers L1 can protect completely the composite structure without any delamination for impact energy lower than 30 J, for a 90% increase in mass. The layer L2 protects it against impacts at energy lower than 50 J, for a 150% increase in mass. The critical protective energies can be determined experimentally. It is between 30 and 50 J for protection L1 and between 50 and 90 J for protection L2. On the other hand, no damage was visible during visual inspection of the plate for 50 and 90 J corresponding respectively to protective layers L1 and L2. However, delamination at mid-thickness was observed during C-scan control. As observed in [START_REF] Petit | Impact and compression after impact experimental study of a composite laminate with a cork thermal shield[END_REF], mid-thickness delaminations should have less effect on the residual strength. However, in order to identify the exact effect of these defects on the residual resistance of the plates, tests of compression after impact (CAI) will be made in a future study.

Finally, this study shows that layers made of hollow spheres core offer a good candidate for protection of composite structures against impact. It is important to underline that for the purpose of In a subsequent study, tests of compression after impact will be made on impacted composite plates. These tests allow identification of the resistance/mass gain obtained when adding the protective layers. On the other hand, another study on modelling of impact on protective layers will also be carried out in order to optimise layers in function of given specifications.

Conclusions

In this paper, an experimental study of protective layers on composite structures submitted to low energy impacts is made. These layers are made up of hollow spheres core and Kevlar woven fabric skin. Two configurations of protective layers designed for 50 and 90 J have been impacted at several impact energies on aluminium support and on composite plates. The different impact curves show that the impactor does not touch the composite plates at the different impact energies used. Furthermore, no defects were visible during C-scan observations on the composite plates protected by: (i) configuration L1 and subjected to impact energy lower than 30 J, (ii) configuration L2 and impacted by energy lower than 50 J. However, delamination was observed when impacting plates protected by configuration L1 and L2, subjected to 50 and 90 J impact energies respectively.

Although the use of protective layer involves a significant increase in weight, this study shows the interest of this concept of protective layers made of hollow spheres core and woven fabric Kevlar skin for the protection of composite parts. Industrial applications for such a protective layer concept are aeronautical composite structural parts expensive to repair or replace, or composite structures for which the weight is not a major problem. It is also important to notice that in the field of aeronautics, this kind of protection is a very good indicator of impact, as the visible impact damage is in the layer, and not in the composite laminate.

Moreover, an advantage in the use of hollow spheres is their easier installation on complex shapes such as non-flat surfaces as well as circular ones. These protective layers can also be changed in case of impact, without changing or repairing the protected composite structure below. This shows the interest in using hollow spheres in the protection of structures against impact.

Of course, the authors are aware that currently the additional mass due to the protection is relatively important, and an optimisation phase of both materials and layer concept is required to obtain a better ratio between protection and weight.

  only resting on a flat rigid foundation. Concerning composite plates: material used is a carbon-epoxy UD prepreg of T700/M21. The stacking sequence for these quasi-isotropic plates is [

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Load/displacement curves for L1 protection configuration, at 30 J energy impact on aluminium support (a) and composite plate (b).

Fig. 7 .

 7 Fig. 7. Comparison of force-displacement curves for composite plates impacted with and without protection.

Fig. 8 .

 8 Fig. 8. Comparison of protective layers impact behaviour between aluminium (A) and composite (C) supports: configuration L1 (a) and configuration L2 (b).

Fig. 9 .

 9 Fig. 9. C-scan control for different tested composite plates (impacted face).

Fig. 10 .

 10 Fig. 10. Delaminated area in composite plate versus impact energy (a) and energy absorbed in the plate (b).

  energy absorption hollow spheres have an advantage over other materials (honeycomb, foam. . .) due to their easy installation on non-flat structures (e.g. cylindrical surfaces).

Fig. 11 .

 11 Fig.11. Focus on delaminated area for impacted composite plates: with protection L2 (90 J) -without protection (30 J).

Fig. 12 .

 12 Fig. 12. Location of delamination in the thickness for impacted composite plates: with protection L2 (90 J) -without protection (30 J).

Al plate 1.4 mm Kevlar thin Kevlar ply 2.1 mm Kevlar

  Fig. 2. The two configurations of protective layers.

	7.5 mm			13.8 mm
	Configuration 1	Configuration 2
	Table 1
	Skin properties.
	Material	Failure stress in compression (MPa)	Failure stress in tension (MPa)	Failure stress in shear (MPa)	Young's modulus (MPa)	Density (kg/m 3 )
	Kevlar (skin)	170		500	150	22,000	1330
	Table 2
	Core properties.
	Material		Failure stress in compression (MPa)	Young's modulus in compression (MPa)	Density (kg/m 3 )
	Hollow spheres (HS)	0.35		30		166

Table 3

 3 Characteristics of protective layers.

	Configuration	Core	Skins	Thickness (mm)	Surface density (g/dm 2 )
	L1	1 Â HS	1.4 mm Kevlar	7.5	55
	L2	2 Â HS	2.1 mm Kevlar	13.8	92

Table 4

 4 Energy repartition in impact tests (A for aluminium and C for composite).

	Name	Plate	Protection	Impact energy (J)	Layer absorbed energy (J)	Impactor rebound energy (J)	Plate absorbed energy (J)
	A-L1-15 J	A	L1	14.8	13.8	1	-
	A-L1-30 J	A	L1	29.9	28.6	1.3	-
	A-L1-50 J	A	L1	50	48.7	1.3	-
	A-L2-15 J	A	L2	14.2	11.6	2.6	-
	A-L2-50 J	A	L2	49.8	48.9	0.9	-
	A-L2-90 J	A	L2	88.6	87.5	1.1	-
	C-15 J	C	-	15.2	-	9.6	5.6
	C-30 J	C	-	30.3	-	10.5	19.8
	C-50 J	C	-	50.2	-	11.8	38.4
	C-L1-15 J	C	L1	15.5	11.3	2.2	2
	C-L1-30 J	C	L1	29.9	18	7.4	4.5
	C-L1-50 J	C	L1	49.5	23.3	13.2	13
	C-L1-70 J	C	L1	69.8	27.1	21.4	21.3
	C-L2-15 J	C	L2	14.6	7.4	5.4	1.8
	C-L2-50 J	C	L2	49.7	39.6	5.6	4.5
	C-L2-90 J	C	L2	89.7	60.2	13.4	16.1
	C-L2-107 J	C	L2	107.1	60.8	19	27.3

Acknowledgements

The authors would like to acknowledge Région Midi-Pyrénées and DRIRE for the funding of this study through EPICEA 2008 program, in project ''bielle HCM''.


Olivier Cherrier

V Pommier-Budinger

SILENCER WITH TUNABLE RESONATORS FOR ACTIVE NOISE CONTROL

Keywords: Resonators, acoustic impedance, holes profiles, active noise control normal incidence acoustic impedance (Z=R+jX) R: normalized acoustic resistance X: normalized acoustic reactance angular frequency pulsation c: sound speed k: wave number equal to / c 

The article deals with active noise control and presents a silencer made of tunable acoustic resonators to reduce variable tonal noises. The silencer is composed of baffles with acoustic resonators made of two superposed and identically perforated plates associated with cavities. One of the plates is mobile and its displacement is controlled by an electromechanical system, allowing changing the internal shapes of the holes of the perforated layers. Consequently the impedance of the resonators and their resonance frequencies can be controlled. These tunable resonators can be used for noise control of variable tonal noises such as those generated by fans for example. The advantages of the proposed system are the simplicity of the actuation and the compactness of the system. NOMENCLATURE l 1 : thickness of plate 1 l 2 : thickness of plate 2 L: cavity depth 2r: perforate hole diameter

INTRODUCTION

The article focus on silencers with acoustic resonators made of perforated sheets bonded onto honeycomb cavities for noise reduction. The advantage of these resonators is that they can be used in hard conditions. The main drawback is that they are selective in frequency. They are thus efficient for tonal noise reduction. In case of variable tonal noise, the efficiency is limited if the resonators are not tunable. One common solution is to control the cavities depth so as to tune the resonant frequency to a desired frequency. Another solution is to change the neck geometry of the resonators. Indeed, several works have shown the influence of the shapes of the holes in the perforated layers: Birnbach et al. [START_REF] Birnbach | Acoustic impedance and absorption coefficient of Helmholtz resonators modified by geometrical configuration variations of inlet neck of two parallel perforated layers[END_REF][START_REF] Birnbach | Increasing of sound absorption proprieties of acoustical panels based on Helmholtz resonators of double plates with orifice inlet neck by geometrical configuration variations[END_REF] studied a resonator with two perforated plates and an inlet air gap, Chanaud [START_REF] Chanaud | Effects of geometry on the resonance frequency of Helmholtz resonators[END_REF] studied the radiation impedance for geometries of non-circular orifices and Tang [START_REF] Tang | On Helmholtz resonators with tapered neck[END_REF] studied a resonator with a tapered neck. This article proposes the study of a baffle silencer with resonators tunable by the change of the neck shape. The resonators are made of two superposed and identically perforated plates associated with cavities [START_REF] Powell | NASA's subsonic jet transport noise reduction research[END_REF][START_REF] Cherrier | Active Absorber with variable resonance frequency for acoustic control[END_REF]. One plate is fixed and bonded to the cavities and the other plate is mobile. The present concept enables to change the internal shapes of the holes of the perforated layers and the resonance frequency of the system. The silencer can thus be used for noise control of variable tonal noises. The advantages of the proposed system are the simplicity of the actuation and the compactness of the system. The first part of the article concerns the design of tunable acoustic resonators. It will be shown that their resonance frequencies shifts to lower frequencies as the displacement of the mobile plate increases and that the resonators can be tuned to attenuate tonal noises that vary on a range of 1000 Hz. In the second part of the article, a silencer composed of baffles with tunable acoustic resonators is presented and experimental results show the efficiency of the system. The last part of the article describes a control law strategy of the active silencer to reduce variable tonal noises.

TUNABLE RESONATORS

Principle

The acoustic resonators made of perforated sheets bonded onto honeycomb cavities are selective in frequency. The acoustic properties are directly linked to the size, shape and porosity of holes and to the thickness of air gaps. In order to make the resonance frequency vary, the holes are generated by two perforated plates, side by side, one fixed and one mobile (Figure 1). The fixed plate (plate 1) is taped on cavity and the mobile plate (plate 2), on top, can move by translation in one direction. The perforations of the two sheets are identical: the same orifice shape (circular in this case) and the same porosity. The translation of the mobile plate generates neck geometry with an elliptic profile as shown in figure 2. The translation values vary from zero to the perforate hole diameter. The minimal distance between two holes is more than one diameter to avoid the case of two overlaps for the same perforation.

Simplified model

In order to describe the acoustic behavior of the variable frequency resonator, a simplified model is presented here, based on the normalized acoustic impedance at normal incidence [START_REF] Jouhaneau | Propagation des ondes acoustiques, Techniques de l'ingénieur (E.T.I. Sciences et Techniques[END_REF][START_REF] Pierce | Acoustics, Acoust. Soc. Am[END_REF][START_REF] Morse | Theoretical Acoustics[END_REF][START_REF] Ingard | Notes on Sound Absorption Technology[END_REF][START_REF] Rayleigh | The Theory of sound[END_REF]. For a resonator, it can be expressed as follows: The normalized impedance of an opened cavity is given by: cot( )

neck cavity Z ZZ  , (1)
cavity Z j kL  , (2)
The normalized impedance of the neck can be expressed by: -the inner surface resistance formed by the acoustic boundary layer in the perforations (plates 1 and 2) during the wave propagation, -viscous interactions near the orifice lips of the circular holes (on every side) and near the elliptic constriction lips between the two plates, -the mobile plate displacement. 2) and ( 5) can be linearized by considering that: >> and >>L.

neck neck neck Z R jX  (3)
(8) Resonance frequency of the resonator resonance f corresponds to the frequency for which the reactance of Z is null. So:

12 2 (l +l ) resonance ie c f L        (9)
When the mobile plate is in movement, previous works [START_REF] Cherrier | Active Absorber with variable resonance frequency for acoustic control[END_REF] have shown that the frequency resonance decreases and shifts towards the low frequencies.

Experimental results

In this paper, the following configuration is tested:

-thicknesses of the plates: 1mm -holes diameter: 1.5mm -cavity depth: 7mm -porosity: 6% For such a configuration, the estimated resonance frequency given by relation ( 9) is equal to 2822 Hz.

Tests are realized in the frequency band 400-3000Hz in an impedance tube. The open end of the tube enables to test samples with a sliding plate (Figure 3). The tube is equipped with two microphones Brüel & Kjaer 4187 for pressure measurements. At the opposite side of the tube, the loudspeaker generates a broadband random noise that propagates as plane waves at 120 dB (SPL in linear domain). The two microphones standard measurement method is used to compute transfer function H between the two microphones and consequently reflection coefficient R expressed by: [START_REF] Ingard | Notes on Sound Absorption Technology[END_REF] where X A and X B are the distances between the microphones and the tested sample.

    H e e H e R B A B A A X X jk X X jk jkX       2
With R, one can compute the normalized acoustic impedance:

R R Z    1 1 ( 11
)
and the absorption coefficient in normal incidence:

2 1 R A   (12)
Figure 4 shows the measured sound absorption coefficient versus the plate displacement. Results are only given for three plate positions to get the figure clear. However, the tests have been achieved for more displacements (every 0.1mm) and the results are similar whatever the displacement. When the mobile plate displacement is null, the measured resonance frequency is 2780 Hz, which is close to the computed value. A part of the difference can be explained by the accuracy in the positioning of the mobile plate displacement. When the mobile plate is moving and the orifice restriction becomes smaller, the resonance frequency of tunable resonators shifts to lower frequencies. The frequency shift is of 1000 Hz for the tested configuration. One can also notice that the absorption band widens with the orifice closing because of the increase of viscous interactions.

BAFFLE SILENCER WITH TUNABLE RESONATORS

Description of the silencer

The silencer is made of three baffles with tunable resonators on each side (Figures 5 and6). There are so six absorbent surfaces. The mobile parts of the three baffles are actuated by a single electromechanical linear actuator (Part B). The link between the actuator and the mobile parts of baffles is rigid, so the displacement of the actuator is the same as the ones of the baffles. A loudspeaker (Part A) is connected to one extremity of the silencer to generate noise. The characteristics of the baffles are as following:

-baffle sizes: 222mm*148mm -gap between 2 baffles: 22mm -variable resonance resonators of same sizes as in section 2

FIGURE 5: SCHEME OF THE SILENCER FIGURE 6: PICTURE OF THE SILENCER Experimental results

The measure of the acoustic attenuation consists in measuring the acoustic sound power of the source with and without the silencer. The acoustic attenuation is the difference between these two acoustic sound powers. These measurements are made on the frequency range 1600-3000 Hz using a sound level meter of type Bruel & Kjaer model 2250. The acoustic attenuation is measured for different values of the mobile plate displacement, from zero to the maximal value (1.4mm) (silencer in active mode). The acoustic attenuation is also measured for no displacement of the mobile plate (silencer in passive mode). Results are synthesized in Figure 7 and show the efficiency of the silencer to better attenuate the lower frequencies than the silencer with no active system. The maximal gain is at 2200Hz: the attenuation is increased by 17 dB with the silencer in active mode compared with the silencer in passive mode.

FIGURE 7: ATTENUATION OF THE SILENCER

CONTROL STRATEGY

The plant to control is composed of the electromechanical actuator linked to the moving part of the silencer. The objective of the control is to minimize the noise at the output of the silencer.

Control in position of the electromechanical actuator

The actuator must be controlled In order to ensure accuracy, the actuator will be controlled by a proportional-integral controller C(s):

s s C 1 10 ) (   (14)
The final scheme for the actuator control is given in Figure 8.

Silencer in active mode

Static model of the silencer Acoustic measurements are performed to establish a static model of the silencer. For frequencies between 2000 and 3000 Hz with a step of 50Hz, the attenuation produced by the silencer is measured for different displacements of the actuator and the displacement for which the attenuation is the highest is identified. A curve (displacement versus frequency) is extracted from these tests (Figure 9) and a polynomial expression of the acoustic behaviour of the silencer is established:

0.042786 5 - 8.2061e - 8 - 6.0139e .... 11 - 1.9119e - 15 - 2.2028e ) ( 2 3 4     f f f f f D (15)
This expression gives the displacement of the actuator that best attenuates an input tonal noise. It is a static model and it is assumed that the mechanical dynamic behaviour of the silencer is taken into account in the actuator model since this model is established from tests of the actuator connected to the silencer. The dynamic of the acoustic phenomena is negligible compared to the mechanical dynamic.

+ Measure ---Polynomial

FIGURE 9: ACOUSTIC BEHAVIOUR OF THE SILENCER Active control of the silencer

The proposed control strategy is as following (Figure 10): 1 -Measurement of the input noise with a microphone 2 -Computation of the tonal noise to be attenuated by a FFT analysis of the measured noise 3 -Computation of the displacement of the actuator that best attenuates this tonal noise by using the polynomial of relation (15)

In order to consider changes in the model of the silencer, the polynomial expression of the silencer behaviour can be up-dated regularly. This strategy has been tested and is efficient with tonal noises that vary between 2000 and 2800Hz. In passive mode, as the displacement of the mobile plate is null, the silencer is configured to attenuate mostly the frequency around 3000Hz. So the output noise is important when the input tonal noise is at 2000Hz and decreases as the frequency of the input noise increases and becomes closer to 2800Hz.

In active mode, the displacement of the mobile plate starts from zero (initial condition imposed by the system) and has to increase to become maximal as the frequency of the tonal noise is 2000Hz (see Figure 9). So, at the beginning of the measure, there is a transition phase for which the level of the noise is the same than in passive mode because the optimal mobile plate displacement is not yet reached. Then Figure 11(c) shows that the noise is much more attenuated than in passive mode as the displacement of the mobile plate is regulated with success (Figure 11(d)). For example, at 5s and for a frequency of 2300Hz, the input noise level is 0.5725V with the silencer in passive mode and 0.115V with the silencer in active mode, so there is so a gain of 14 dB. At the end of the measure, the noise measured in passive mode and in active mode is the same since the mobile plate is in the optimal position, i.e. null, in both cases.

FIGURE 11: EXPERIMENTAL RESULTS OF THE SILENCER WITH TUNABLE RESONATORS

CONCLUSION

The article has shown that resonators made of two superposed and identically perforated plates associated with cavities can be tuned by moving one of the two perforated plates. The resonance frequency then shifts to lower frequencies when the profile of the resonators' neck is modified to be restricted. A baffle silencer made with this kind of resonators can be efficient to attenuate tonal noise that varies on a frequency range of 1000 Hz. A control strategy based on the position control of the silencer and on the identification of the acoustic behavior of the silence can be implemented to achieve active noise control. Improvements could be made to the system by optimizing the actuator. The proposed active noise control system is easy to implement and is different of conventional active solutions based on the interference with secondary noise sources. It can be used for machines with variable tonal noises such as fans of aircraft engines and rotating machines controlled at variable speed.

FIGURE 1 :

 1 FIGURE 1: PRINCIPLE OF THE TUNABLE RESONATOR WITH TWO PERFORATED PLATES

  normalized acoustic resistance and reactance of the resonator neck.

  neck Xdepends on the length of the cavity, on the plate thickness and on the displacement of the mobile plate..

  by using results given[START_REF] Jouhaneau | Propagation des ondes acoustiques, Techniques de l'ingénieur (E.T.I. Sciences et Techniques[END_REF][START_REF] Pierce | Acoustics, Acoust. Soc. Am[END_REF][START_REF] Morse | Theoretical Acoustics[END_REF][START_REF] Ingard | Notes on Sound Absorption Technology[END_REF][START_REF] Rayleigh | The Theory of sound[END_REF][START_REF]Acoustics of ducts and mufflers with application to exhaust and ventilation system design[END_REF][START_REF] Xu | Dual Helmholtz resonator[END_REF]. When the holes of the two perforated plates are aligned (null displacement of the mobile plate), we get:

FIGURE 3 :

 3 FIGURE 3: SCHEMATIC OF EXPERIMENTAL SET UP

FIGURE 4 :

 4 FIGURE 4: MEASURED ABSORPTION COEFFICIENT FOR THREE MOBILE PLATE DISPLACEMENTS P1: mobile plate displacement = 0mm P2: mobile plate displacement = 0.8mm P3: mobile plate displacement = 1.4mm

  in position in order to ensure a good positioning of the mobile plates of the resonators. The model of the actuator with its load (the moving part of the silencer) is identified by tests in openloop. It can be described by transfer function H A (s): D the displacement and U the input voltage of the actuator.

  FIGURE 8: CONTROL OF THE ELECTROMECHANICAL LINEAR ACTUATOR

FIGURE 10 :

 10 FIGURE 10: CONTROL STRATEGY OF THE SILENCER

Figure 11 (

 11 c) shows the same voltage when the silencer is in active mode and when the mobile plate translates to track the tonal noise.

  Microphone at the output of the silencer in passive mode (Microphone at the output of the silencer in active mode (V) Mobile plate displacement in active mode (mm)(m)


F Zawawi

J Morlier

G Grondin

J M Moschetta

Design of Optimum Torsionally Flexible PropRotors for Tilt-Body MAVs

Keywords: tilt-body aircraft, micro-air vehicles, strip theory, minimum energy loss, low Reynolds number, optimum propeller/rotor blades, passively adaptive blades, long flight endurance

This paper presents a methodology to design the optimum proprotor for tilt-body microair-vehicles (TB-MAV) with efficient global propulsion system and long flight endurance in both cruise and hover modes. The TB-MAV developed at ISAE, which is called MAVion, was used as a baseline in the design process. To acquire maximum performance of TB-MAV's global propulsion system, an efficient optimization process of the proprotor propulsion system was carried out. The optimization process consists of two-step inverse design methods. The first step determines the optimal operating conditions in terms of power and rotational speed of proprotor and the second step designs the optimal blade geometry in terms of twist angle distribution. The optimal blade twist distribution along the blade was computed using the Glauert's strip theory for minimum energy loss condition. Meanwhile, the optimal operating conditions were determined by the motor outputs corresponding to high motor efficiency. A comparison of performance in terms of total efficiency and flight endurance between the optimized flexible proprotor, the optimized rigid proprotor, optimized propeller and optimized rotor is presented.

Introduction

Tilt-body micro-air-vehicles (TB-MAV) have been regarded over the last decade as promising candidates for convertible multi-tasking configuration as they intrinsically combine both vertical and horizontal flight capabilities. Such convertible MAVs can provide better observation abilities during hovering while offering faster forward movement during cruising. However, in these two flight modes the rotor operates under considerably different conditions where using a conventional rotor will require high power consumption or even inability to provide sufficient required thrust in the off-design cases. The substantial variations of blade pitch and twist that are required to offer highest possible efficiency at both flight modes require an active control system to change the shape according to the flight configurations. Nonetheless, the active blade shape modification concept is inapplicable for small-sized rotor (blade ~3mm of thickness) due to its complexity and the unavailability of control devices in the small size range, and requires non-mechanical control system to be used. Structural topology optimization where the mass density distributions along the blade is optimized according to the desired aerodynamic loadings [START_REF] Wang-Yu | Adaptive Bend-Torsional Coupling Wind Turbine Blade Design Imitating the Topology Structure of Natural Plant Leaves[END_REF] and off-axis carbon fiber composites for blade twisting moment introduction, are among the potential concepts to structurally design the passively adaptive blades. Thus, the idea behind the passive blade shape modification concept is to compute the optimal twist distribution for TB-MAV hingeless proprotor for both cruise and hover mode. The difference of optimized blade twist between these two flight modes and its corresponding resulting aerodynamic loadings to enable twist deformation, can be used as design inputs to structurally design the flexible blade by using inverse methods of blade structural stiffness distribution identification [START_REF] Li | Identification of Stiffness Distribution and Damage in Euler-Bernoulli Beams Using Static Response[END_REF].

Blade Twist Optimization

A numerical iterative method for computing the optimal blade twist angle distribution along the blade span in horizontal and vertical flight, with identical chord laws, was developed adapting the inverse design method of proprotor using Blade Element Momentum Theory (BEMT) based on Minimum Induce Loss (MIL) that satisfies the Betz's condition. The systems of equations used were proposed by Adkins and Liebeck [START_REF] Adkins | Design of Optimum Propellers[END_REF] and the design procedure is shown as a block diagram in Fig. 2 (lower part). This iterative design procedure determines the optimal blade twist angle along the radius (β(r)), the characteristics of the designed proprotor (C T , C P ) and its performance ( proprotor , FM). The radial and axial interference factors, Re, and relative Mach number depending on the operating condition are also iteratively computed in the design process. The details are described in [START_REF] Adkins | Design of Optimum Propellers[END_REF].

In the classical approach of propeller design and analysis, different assumptions are made, i.e. lift polar is a linear function, in order to obtain a closed loop solution or even to quickly converge the iterative process. Hence, in order to consider non-linear airfoil characteristics prevalent in low Reynolds number regime (Re<70,000), XFOIL, an airfoil design and analysis code developed by Drela [START_REF] Drela | XFOIL: An Analysis and Design System for Low Reynolds Number Airfoil, Low Reynolds Number Aerodynamics[END_REF] was integrated into the iterative process of design to eliminate this assumption as well as to precisely compute the aerodynamic coefficients. The airfoil sensitivity to the turbulence level and Reynolds number variation was eliminated by using cambered plate [START_REF] Hein | Hover Performance of a Micro Air Vehicle: Rotor at Low Reynolds number[END_REF] as a blade profile. The cambered plate profile also promises better performance than other airfoils in MAVs Reynolds number regime [START_REF] Laitone | Wind tunnel tests of wings at Reynolds numbers below 70,000[END_REF]. In this respect, the n crit parameter, used to define the turbulence level in XFOIL, was set at 0.1 for both flight modes. The cambered plate with thickness t/c of 2.5%, camber f/c of 3% and camber location xf/c of 30% was used in this study. For optimization, the input parameters are number of blades B, proprotor radius R, operational speed V, rotational speed Ω and mechanical power into the proprotor P mech . Constant non-dimensional chord c/R of 0.15 was used to ensure the activity factor AF of the blade is in the range of practical propeller, 90≤AF≤200 [START_REF] Raymer | Aircraft Design: A Conceptual Approach, 3rd[END_REF]. The proprotor diameter D used for MAVion was set at 18cm to ensure the momentum injected into the flow by the propwash covers almost the entire span of the MAVion's body. This will avoid the flow separation on the MAV's body while performing a transition between the two flight modes. To evaluate the performance of proprotor, the computed quantities are the proprotor in forward flight and the figure of merit FM in hover as defined in Eq. 1. where J is the advance ratio, C T is the thrust coefficient and C P is the power coefficient, which are respectively defined in Eq. 2 as:

J= 60. V Ω.D , C T = T ρ. Ω 60 2 .D 4 , C P = P mech ρ. Ω 60 3 .D 5 (2)
By knowing the performance of proprotor, the total in cruise and hover can be respectively computed by using relationships in Eq. 3.

η total(cruise) =η motor .η proprotor , η total(hover) =η motor .FM

Operating Condition Optimization

In order to control the global propulsion system to be in high performance in both flight modes, the design range of operating conditions was first specified. The operating conditions in the design process of optimal blade twist are the motor outputs in terms of P mech and Ω corresponding to the specified motor . The motor voltage V m must be carefully chosen to satisfy the design requirements. In this study, the initial V m was set up based on the performance data of MAVion which is the motor outputs at 10V with 60%≤ motor ≤ motor(max) were set, as can be seen in Fig. 4. However, the designed proprotor blade at given operating conditions from the motor outputs is subject to the design requirements. Since two proprotors were used in the MAVion in order to control the yaw direction, the T req at cruise flight was constrained by the half of MAVion's drag at the target cruising velocity V=16m/s. Whereas for the hover flight, T req was constrained by the half of the MAVion's weight of 350 grams. The design requirement for MAVion proprotor at each flight configuration is shown in Table 1. The design constraint considered in this study is the twist deformation between cruise and hover, ∆β= |∆β (cruise) -∆β (hover) |. The value needs to be as small as possible in order to be able to manufacture the passive control blade, also the desired twist deformation to be possibly adjusted by the aerodynamic loadings. For the design of this study, the ∆β was set as ≤8°. According to the theory of propeller design based on MIL, by increasing the proprotor Ω in the cruise can decrease the optimal blade twist angle. In contrast in the hover flight mode, the optimal blade twist increases as the Ω decreases. This will be advantageous in the design of optimal blades for proprotor where the ∆β can be minimized. From the first order motor model law, higher the voltage can provide maximum motor at higher Ω as illustrated in Fig. 3. Hence, by freezing the motor outputs at around maximum motor can promise long flight endurance due to low P total [8]. A direct current motor of AXI2204/54EVP with motor speed constant =1400 rev/min/Volt, zero load current I =0.35A and internal resistance R m =0.32 ohm, was used in this study. The motor output variables at constant V m were computed using relationships in Eq. 4 [8] and are plotted in Fig. 3 and Fig. 4.

P total = V m - Ω K v R m .v, P mech = V m - Ω R m , η motor = P total P mech (4)
Fig. 3 AXi2808/04 motor efficiency at different voltages Fig. 4 Motor output variables vs motor speed at 10V of applied voltage

From the known P total supplied into proprotor, the endurance of the flight is computed by the relationship in Eq. 5.

Endurance time, t end = E P total = C.K.V b P total [START_REF] Hein | Hover Performance of a Micro Air Vehicle: Rotor at Low Reynolds number[END_REF] where E is the battery energy stored, C is the battery capacity, K is the number of battery cells and V b is the nominal battery voltage. For MAVion, three cells of Polymer Lithium Ion battery (V b =3.7Volts/cell, C=1000mAh) were used in MAVion for its small size and light weight as well as high energy density for long endurance flight.

Results and Analysis

The developed iterative process of proprotor analysis using BEMT method was first validated with the results given by QPROP and corresponding experiment measurement. A 2-bladed rotor with a rectangular planform (R=0.1m, c=0.02m), flat plate blade profile with a thickness t/c=2.5% and blade pitch of 25° was used in the validation process. From the validation result in Fig. 5, it is found that the BEMT method with the integration of XFOIL shows good agreement with the experimental measurement. QPROP, a propeller analysis tool with simplified and limited lift polar, showed unacceptable correlations due to its incapability to predict aerodynamic coefficients at poststall region. The design of propulsion system in hover and cruise was done separately. The optimal blade twist within the specified design range of operating conditions was obtained. The results for hover case are shown by Fig. 6 and Fig. 7, and the results for cruise case are shown by Fig. 8 and Fig. 9. To generate the T req in a hover, the rotational speed of Ω=12,870 and the designed blade corresponding to this Ω is shown by a dashed line in Fig. 6. Whereas to attain T req in cruise at V c =16m/s, the proprotor needs to operate at Ω=13,500 can be observed in Fig. 9. The design at V c >16m/s was done to find the maximum V c that passive control blade can achieve. And it is found that the MAVion is not able to cruise beyond 16m/s using passively flexed proprotor. Fig. 10 represents the optimal blade twist distribution in hover and cruise modes that satisfy the design constraint (∆β≤8°) for flexible proprotor. The performance of the rigid proprotor was also computed. As for a rigid proprotor, the optimized blade twist of the flexible proprotor in cruise was chosen as a blade twist of the rigid proprotor, as shown in Fig. 10. The optimized blade geometry of single-point propeller/rotor design in terms of chord and twist angle distribution is shown in Fig. 11. And the performances of the optimized propeller and optimized rotor in cruise and hover respectively were used as a baseline of comparison with the designed proprotor as shown in Table 2. The optimized propeller able to hover the MAVion but consumes very high power which reduces the endurance, meanwhile, the optimized rotor stalls the MAVion in forward flight. The optimized rigid proprotor reduces vey slight its performance as compared to the optimized flexible proprotor even though it is needed to be operated at Ω=7,850 to obtain sufficient thrust in hover. The optimized flexible and rigid proprotor offers ~30% increment of total in hover as compared to the optimized propeller. It also promises flight endurance at both flight modes as long as the optimized rotor and optimized propeller in hover and forward flight respectively.

Conclusions

A practical and reliable methodology to design the proprotor for the application in tilt-body micro-air-vehicles (TB-MAV) with efficient global propulsion system and long flight endurance was presented. To acquire maximum total , the proprotor blade was optimized in terms of twist angle distribution along the blade using MIL condition to obtain the maximum proprotor while the operating condition was optimized to obtain the maximum motor. The presented efficient design methodology which combines both proprotor and motor design could be used as a practical design tool for the passive control proprotor development.

Fig. 1

 1 Fig.1 MAVion (Dimensions: 400mmx220mm) in hover configuration

Fig

  Fig. 2 Proprotor design methodology

Fig. 6 Fig. 7

 67 Fig. 5 Comparison of BEMT simulation and the corresponding measurements

Fig. 10

 10 Fig. 10 Blade geometry of flexible proprotor Fig. 11 Optimal blade shapes at single design point in hover and in cruise

Table 1 Design requirements for proprotor Flight configuration Velocity, V (m/s) Required thrust, Treq (N)

 1

	Hover	0	1.7
	Cruise	16	0.31

Table 2

 2 Results comparison of propulsion performance between propeller, rotor and proprotor

	Proprotor type	Flight mode	total (%)	P total (W) Endurance, t end (min)
	Single design point	Cruise	47.62	10.42	63.89
	(Optimized propeller)	Hover	0.058	127.44	5.23
	Single design point	Hover	61.74	20.33	32.75
	(Optimized rotor)	Cruise (stall)	-	-	-
	Flexible proprotor	Cruise	44.57	11.16	59.67
		Hover	31.78	25.22	26.40
	Rigid proprotor	Cruise	44.57	11.16	59.67
		Hover (Ω=7,850)	30.19	26.00	25.61


Ahmed Jhinaoui
email: ahmed.jhinaoui@inria.fr

Laurent Mevel
email: laurent.mevel@inria.fr

Joseph Morlier
email: joseph.morlier@isae.fr

Subspace Identification for Linear Periodically Time-varying Systems

Keywords: Helicopter dynamics, Time-varying systems, Periodic motions, Subspace methods, Stability analysis

In this paper, an extension of the output-only subspace identification, to the class of linear periodically time-varying (LPTV) systems, is proposed. The goal is to identify a useful information about the system's stability using the Floquet theory which gives a necessary and sufficient condition for stability analysis. This information is retrieved from a matrix called the monodromy matrix, which is extracted by some simultaneous singular value decompositions (SVD) and from a resolution of a least squares criterion. The method is, finally, illustrated by a simulation of a model of a helicopter with a hinged-blades rotor.

INTRODUCTION

Over the last forty decades, subspace identification methods have enjoyed some popularity and numerous applications of these methods have emerged in civil engineering, aeronautics and many other fields. The stochastic subspace identification (SSI) consists in obtaining the modal parameters (natural frequencies, modal damping ratios and mode shapes) of a system subject to ambient excitations, by some geometric manipulations and projections of data given by sensors measurements or input measurements. There exist, then, two types of identification algorithms: output-only and input-output algorithms. A comprehensive overview of different approaches of SSI can be found in [START_REF] Van Overschee | Subspace Identification for Linear Systems, Theory, Implementation[END_REF].

Unfortunately, most of research interest on subspace methods has been given to linear time-invariant (LTI) systems. In contrast, the literature on linear time-varying (LTV) case is not abundant. However, most physical phenomena exhibit time varying behaviors, mainly due to internal (fatigue...) or external (disturbances...) operating conditions. One particular subclass of LTV systems is the linear periodically time varying (LPTV) systems which are widely common in communications, circuit modeling and rotating machines such as wind turbines and helicopters' rotors.

In time-domain, the few works that have been carried out, in order to extend the subspace-based methods to these LPTV cases, can be categorized into two main approaches. The first approach consists in identifying the considered system recursively using adaptive algorithms which is appropriate only for slowly-varying dynamics and when a priori information about the variation behavior is avail-able, whereas, the second approach suggests to find a set of output subsequences that have time-invariant behaviors. This makes it possible to derive an identification algorithm for these subsequences, which is close to the classical timeinvariant algorithm. Among the few attempts on the subject, one can cite [START_REF] Liu | Identification of linear time-varying systems[END_REF] where a discrete time domain state space model, transition matrices and pseudo-modal parameters are used to describe and identify recursively periodic systems, and [START_REF] Verhaegen | A class of subspace model identification algorithms to identify periodically and arbitrarily time-varying systems[END_REF] in which subspace model identification algorithms that allow the identification of an LPTV state space model from a set of input-output measurements are presented.

In this paper, the focus is to develop a new identification algorithm for LPTV systems, with two main features: first, the algorithm can be applied when the system is subject to unknown excitations and only output data are available. Second, it aims at identifying useful parameters for stability analysis, by the Floquet theory, with no need to identify the system over a whole period as it is usually done in periodic cases.

The paper is organized as follows: in Section 2, a typical SSI algorithm is presented. Section 3 gives the essential elements of the Floquet theory. Then, Section 4 is devoted to the design of the new LPTV extended method. Finally, in Section 5 the efficiency of the method is tested on a numerical simulation of a model of a helicopter with a hinged-blades rotor.

TYPICAL OUTPUT-ONLY SUBSPACE IDENTIFICATION

In this section, a typical output-only SSI algorithm, based on covariance-driven data, from [START_REF] Benveniste | Single sample modal identification of a non-stationary stochastic process[END_REF] is recalled. Let consider the LTI discrete-time state space model of a given system:

z k+1 = F z k + v k y k = H x k + w k (1)
where x ∈ R n is the state vector, y ∈ R r the output vector or the observation, F ∈ R n×n the state transition matrix and H ∈ R r×n the observation matrix. The vectors v and w are noises assumed to be white Gaussian. For chosen parameters p and q such that min{pr, qr} ≥ n, the covariance-driven Hankel matrix below is built:

H p,q =     R 1 R 2 • • • R q R 2 R 3 • • • R q+1 . . . . . . . . . . . . R p+1 R p+2 • • • R p+q     (2)
The covariances of the output data write

R i = E(y k y T k-i )
, where E is the expectation operator. The R i 's can be estimated by Ri = 1 N N k=i+1 y k y T k-i , where N is the number of output measurements that are available (N 1). Then, let the estimated Hankel matrix be:

Ĥp,q = 1 N N -p k=q Y k + Y k -T (3)
where:

Y k + = y T k • • • y T k+p T , Y k -= y T k-1 • • • y T k-q T (4) Let G = E(z k y T k
) be the correlation between the state and the observation,

O p = H T , (HF ) T , • • • , (HF p ) T T and C q = F G, F 2 G, • • • , F q G the
p-th order observability matrix and the q-th order (shifted) controllability matrix. The computation of the R i 's leads to the decomposition (see [START_REF] Benveniste | Single sample modal identification of a non-stationary stochastic process[END_REF]):

H p,q = O p C q
(5) Therefore, the observability matrix O p can be obtained with a thin singular value decomposition (SVD) of the Hankel matrix H p,q and its truncation at the desired model order n [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]):

H p,q = U ∆V T = [ U 1 U 2 ] ∆ 1 0 0 ∆ 2 V 1 T V 2 T (6) O p = U 1 ∆ 1 1 2
(7) where ∆ 1 contains the first n singular values and U 1 the n first columns of U . The observation matrix H is extracted from the first r rows of the observability matrix O p . The state transition matrix F is obtained from a least squares resolution of [START_REF] Benveniste | Single sample modal identification of a non-stationary stochastic process[END_REF]):

O ↑ p F = O ↓ p ( 8
) where:

O ↑ p =     H HF . . . HF p-1     , O ↓ p =     HF HF 2 . . . HF p     (9)
Finally, the eigenstructure of the system (1) is retrieved from:

det(F -Iλ) = 0, F φ λ = λφ λ (10)
(λ, φ λ ) denote the eigenvalues and the eigenvectors (also, called mode shapes) of the system.

FLOQUET THEORY

The Floquet theory is a mathematical theory of ordinary differential equations (ODEs). Introduced by [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF], it is the first complete theory for the class of periodically time-varying systems. In this section, some of its essential elements, that are related to the study hereafter, are briefly reviewed. More details can be found in [START_REF] Dacunha | A unified floquet theory for discrete, continuous, and hybrid periodic linear systems[END_REF].

Let consider the periodic differential system:

ẋ(t) = A(t)x(t) (11)
where x ∈ R n is the state vector. The state transition matrix A ∈ R n×n is continuous in time (or at least, piecewise continuous) and periodic, of period T > 0. If an initial condition x(t 0 ) = x 0 is fixed, a solution of ( 11) is guaranteed to exist. Let Φ(t) be the matrix whose n columns are n linearly independent solutions of (11), Φ(t) is known as the Fundamental Transition Matrix (FTM). It has the properties:

Φ(t) = A(t)Φ(t), Φ(t + T ) = Φ(t)Φ(T ), ∀t

Stability Analysis

Let Q be the value of the fundamental matrix at t = T :

Q = Φ(T ), R = 1 T log(Q) (12)
Q is called the Monodromy Matrix . According to Floquet's theory, the dynamical system ( 11) is stable if and only if the eigenvalues of R are negative or, similarly, if the norms of the eigenvalues of Q are less than one.

Floquet Transformation

The Floquet transformation, also called the Lyapunov-Floquet transformation, gives an underlying autonomous system (a system with a constant state transition matrix with respect to time) that is equivalent to the initial periodic system i.e: the transformation is invertible.

If the change of variable x(t) = Φ(t)e -Rt z(t) is made, the theory insures that:

ż(t) = Rz(t) (13)
Assume that the equation of observation of the considered system is: y(t) = Cx(t), where C ∈ R r×n . The equivalent equation for the new variable z is:

ẏ(t) = Cx(t) = CΦ(t)e -Rt z(t) (14)
It is easy to demonstrate that Φ(t + T )e -R(t+T ) = Φ(t)e -Rt using the fact that R = 1 T log(Φ(T )). Therefore, any periodic system can be transformed into an equivalent autonomous system with an equivalent periodic observation matrix.

SUBSPACE IDENTIFICATION FOR LPTV

In this section, an extension of SSI to the linear periodically time-varying case is suggested. Let consider the periodic state-space system, of a period T :

ẋ(t) = A(t)x(t) + v(t), A(t + T ) = A(t), ∀t ∈ R y(t) = Cx(t) + w(t) (15)
where v and w are two white noises. As outlined in Section 3, an equivalent representation can be given by the Floquet's transformation:

ż(t) = Rz(t) + Φ(t)e -Rt -1 v(t) y(t) = C (t)z(t) + w(t), C (t + T ) = C (t), ∀t (16)
where C (t) = CΦ(t)e -Rt is periodic of period T . Sampling at rate 1 τ (τ > 0 and T assumed to be a multiple of τ ), yields the discrete time model below [START_REF] Kim | Robustness analysis of linear periodic time-varying systems subject to structured uncertainty[END_REF]):

z k+1 = F z k + Γ k w k y k = H k z k + v k (17)
The discretized monodromy matrix F = e Rτ , the discretized matrix of observation

H k = C (kτ ) and Γ k = (k+1)τ kτ e Rγ dγ are periodic. The discrete period is T d = T τ . If µ is an eigenvalue of the continuous system R, the corresponding eigenvalue λ of F is such that: λ = e τ µ .
According to the theory of Floquet, the continuous system is stable when all the eigenvalues µ of R are negative. The discrete system is, then, stable when the norms of the eigenvalues λ are inferior to one. The goal herein is to identify these eigenvalues and, thus, to draw conclusions the stability. p,q is defined as:

H (j) p,q = 1 N N -1 i=0 Y + j+iT d Y -T j+iT d (18)
Fact 1. Notice that T is assumed to be a multiple of τ . Otherwise, the time-invariant subsequences are the ensembles (y

H (j) p,q =     H j H j+1 F . . . H j+p F p     1 N F N -1 i=0 z j-1+iT d y T j-1+iT d • • • 1 N F q N -1 i=0 z j-q+iT d y T j-q+iT d + N ( 19
)
where N is a function converging to zero when N goes to infinity.

Proof. See the appendix. Corollary 3. As in the time-invariant case, the Hankel matrix H (j)

p,q converges to the decomposition below:

H (j) p,q = O (j) p C (j) q (20)
The j-th observability and (shifted) controllability matrices are defined as:

O (j) p = H T j , (H j+1 F ) T , • • • , (H j+p F p ) T T (21) C (j) q = F G (j-1) , F 2 G (j-2) , • • • , F q G (j-q) (22)
where G (j) is the state-output cross correlation of the jth invariant subsequence. It can be estimated by

G (j) = 1 N N -1 i=0 z j+iT d y T j+iT d . The observability matrix O (j)
p can be obtained as in Section 2 via an SVD of H (j) p,q and its truncation at the desired model order n:

H (j) p,q = U (j) ∆ (j) (V (j) ) T = U (j) 1 U (j) 2 ∆ (j) 1 0 0 ∆ (j) 2 (V (j) ) T (23)
where

U (j) 1
is the n first columns of U (j) , and ∆ (j) 1 the n upper singular values of ∆ (j) . The size of the Hankel matrix is (p + 1)r × qr. The complexity of the SVD does not depend on N and therefore, the computation is not cumbersome. The observation matrix H j is obtained from the first r rows of O (j) p . Notice that the extraction of the transition matrix from one observability matrix is no longer possible. Indeed, in order to get F , two successive matrices of Hankel should be computed. For the (j + 1)-th time-invariant subsequence, H (j+1) p,q writes:

H (j+1) p,q =     H j+1 H j+2 F . . . H j+p+1 F p     F G (j) • • • • • • F q G (j+1-q)
Differently from Section 2, the transition matrix F is a least squares solution of an equation involving two different observability matrices:

O ↑ p (j+1) F = O ↓ p (j)
(24) where:

O ↓ p (j) =     H j+1 F H j+2 F 2 . . . H j+p F p     , O ↑ p (k+1) =     H j+1 H j+2 F . . . H j+p+1 F p-1     Both O (j) p
and O (j+1) p should be expressed in the same basis of the state space system. In fact, the left part of the singular value decomposition gives the observability matrix up to some invertible (change of basis) matrix. For instance, in the time-invariant case, the estimated observability matrix Õp = U 1 ∆ 1

1 2 = O p T ,
where T is a non-singular transformation. In the case of timevarying systems, two successive SVD are needed. The corresponding invertible matrices T (j) and T (j+1) will be different:

Õ(j) p = U (j) 1 (∆ (j) 1 ) 1/2 = O (j) p T (j) (25) Õ(j+1) p = U (j+1) 1 (∆ (j+1) 1 ) 1/2 = O (j+1) p T (j+1) (26)
where none of the matrices T (j) or T (j+1) is known. The least squares resolution will results in:

F = ( Õ↑(j+1) p ) † Õ↓(j) p = (T (j+1) ) -1 (O ↑(j+1) p ) † O ↓(j) p T (j)
(27) where the symbol † denotes the Moore-Penrose pseudoinverse. The left and the right bases of the identified transition matrix are not the same. A solution for this problem is given in [START_REF] Liu | Identification of linear time-varying systems[END_REF]:

Let L (j) and L (j+1) the first l (l ≥ n r + 1) block rows of Õ(j)

p and Õ(j+1) p . For l = 2, and without loss of generality:

L (j) = H j H j+1 F T (j) (28) L (j+1) = H j+1 H j+2 F T (j+1) (29)
The 3 successive observations matrices can be considered as close: H j+2 H j+1 H j . Therefore:

(L (j) ) † L (j+1) = (T (j) ) -1 H j H j+1 F † H j+1 H j+2 F T (j+1) (T (j) ) -1 T (j+1) (30)
The number of block rows l is chosen such that L (j) is full rank. Finally, an approximate of the state transition matrix, with the same left and right transformation, writes:

F = (L (j) ) † L (j+1) F (T (j) ) -1 (O ↑(j+1) p ) † O ↓(j)
p T (j) (31) Once F and H j are computed, the modal structure of the system (17) are obtained from the resolution of det(F -Iλ) = 0, F φ λ = λφ λ and ϕ λ = H k φ λ .

Algorithm: To sum up, here are the steps of the suggested LPTV identification method:

• (N T d + p + 1) output-measurements are available.

The data are set in ensembles of time-invariant series (y j+iT d ) • compute, for the j-th and the (j +1)-th subsequences, the two successive Hankel matrices H (j)

p,q and H (j+1) p,q using the formula (18) • compute the SVD's of H (j) p,q and H (j+1) p,q . Then, retrieve Õ(j)

p and Õ(j+1) p as in ( 25) and (26) • compute the discrete transition (monodromy) matrix F , using ( 31) and (30) • given F and H j , compute the eigenstructure of the system (eigenvalues, observed eigenvectors)

ILLUSTRATIVE EXAMPLE

In this section, two illustrative examples are studied in order to test the suggested identification algorithm.

Simulation Model

The studied system herein is a model of a helicopter with a hinged-blades rotor. The system is modeled as in [START_REF] Coleman | Theory of self-excited mechanical oscillations of helicopter rotors with hinged blades[END_REF] and [START_REF] Sanches | Instability zones for isotropic and anisotropic multibladed rotor configurations[END_REF]. The helicopter's fuselage is considered to be a rigid body with mass M , attached to a flexible LG (landing gear) which is modeled by two springs K bx and K by , and two viscous dampers C x and C y as illustrated in Fig. 2. The

. Let Z = [ x y φ 1 • • • φ N b ]
T . A linear model for such rotating systems under free vibrations is defined by (see [START_REF] Christensen | Design of active controlled rotor-blade systems based on time-variant modal analysis[END_REF]):

M(t) Z(t) + C(t) Ż(t) + K(t)Z(t) = 0 (32)
where M(t), C(t) and K(t) are the system's mass, damping and stiffness matrices. The system can be written as in (15), with x = Z T ŻT T and

A(t) = 0 I -M(t) -1 K(t) -M(t) -1 D(t) A(t) = A(t + 2π Ω ). It is periodic of period T = 2π Ω .
For the study herein, the number of blades is fixed to N b = 4 and the numerical values used for the application are reported in Table 1 The helicopter model is simulated for Ω = 10 rad.s -1 . To get precise results, the subspace identification should be processed on very large dataset, in order to overcome the problem of bias due to the noise. Data points are generated over a samples number of 2000.T d , with the discrete period T d = 100 (the sampling period is τ = T /T d ). The order of the system is known and is equal to n = 12. The suggested identification algorithm is applied to the data as explained in Section 4 (with chosen q = p+1 = 11). The summary of the identified eigenvalues, and the true values (the true values are computed from a numerical resolution of the state equation of the system and the corresponding monodromy matrix, using the Matlab function for ordinary differential equation's resolution ode45), are given in Table 2.

Table 2. Identified eigenvalues vs. true ones

Identified eigenvalues

True eigenvalues 0.9985 ± 0.0225i 0.9987 ± 0.0218i 0.9988 ± 0.0196i 0.9988 ± 0.0198i 0.9988 ± 0.0138i 0.9995 ± 0.0132i 0.9998 ± 0.0151i 0.9996 ± 0.0150i 0.9992 ± 0.0149i 0.9996 ± 0.0150i 0.9995 ± 0.0142i 0.9995 ± 0.0141i

We get 12 eigenvalues that are conjugate. The differences between the obtained real parts and the true ones are less than 0.1 %. For the imaginary parts, the errors are greater.

The identification goes more accurate when the number of data samples is larger.

The SSI algorithm presented in Section 2 (from equation (3) to ( 10)), is also applied in order to compare it to the new suggested method. A total of 12 non conjugate modes are identified (see Table 3). Some of the identified eigenvalues have negative real parts. For a discrete system these real parts must be positive, because the discrete state transition matrix is an exponential of the continuous one. Therefore, applying the time-invariant SSI algorithm gives irrelevant results. This was predictable because a LTI algorithm is applied to a set of subsequences that have different dynamical behaviors. Finally, the suggested algorithm in Section ( 4) is applied to the case where Ω = 0.1 rad.s -1 in order to compare the identified eigenvalues to the true eigenvalues of the system when the rotor is turned off (Ω = 0). In fact, when Ω = 0, the state matrix is no longer time-dependent The identified eigenvalues for Ω = 0.1 rad.s -1 are, indeed, close to the true ones when the rotor is turned off. When Ω converges to 0, the results of the identification in the periodic case converges to the true time-invariant system.

A(t) = A.
Another feature of interest is to show that the identified state transition matrix is, indeed, invariant with respect to time. To verify it, the system is identified at different instants of time. The variations of the real parts and the imaginary parts of the identified eigenvalues with respect to time are reported in Fig. 3 The problem of identification for linear periodically timevarying is addressed. It has been shown that, in order to get an information about the stability of such systems, a so-called monodromy matrix should be identified. This matrix is derived from the Floquet transformation which gives an equivalent description of periodic systems. Based on this description, an extension of the covariance-driven SSI algorithm is suggested. The algorithm was successfully tested for simulation data.

Future works encompasses the estimation of the uncertainties about the identified parameters for the suggested identification algorithm, given herein.

Following

  Fig. 1. The j-th time-invariant subsequence time-invariant case, because each subsequence describes a different dynamical behavior from the others. The Hankel matrix should be redefined using each one of these timeinvariant subsequences, separately. For any j ≥ q and if (N T d + p) data samples are available, the j-th Hankel matrix H (j)

Fig. 2 .

 2 Fig. 2. Mechanical model of a helicopter with 4 blades rotor spinning with a velocity Ω, is articulated and the offset between the MR (main rotor) and each articulation is noted a. The blades are modeled by a concentrated mass m at a distance b of the articulation point. Torque stiffness and a viscous damping K φ and C φ are present into each articulation. The moment of inertia around the articulation point is I z . The degrees of freedom are the lateral displacements of the fuselage x and y, and the out-of-phase angles (φ k ) k=1•••N b , with N b the number of blades. Let Z = [ x y φ 1 • • • φ N b ]T . A linear model for such rotating systems under free vibrations is defined by (see[START_REF] Christensen | Design of active controlled rotor-blade systems based on time-variant modal analysis[END_REF]):

  and Fig. 4. These values are, as shown in the figure, constant over time. Therefore, the identified state transition matrix is constant as stated in (17).

Fig. 3 .

 3 Fig. 3. Eigenvalues variations: real parts of the eigenvalues vs. sample number

  jτ +iT τ ) i and the data Y + j+iT d and Y - j+iT d are replaced by Y + From the periodicity of the observation matrix H, i.e. H j+iT d = H j for all non negative integer i, one deduce the following property for the Hankel matrix:

	jτ +iT	and Y -jτ +iT	(where	is the floor
	τ	τ
	operator), in order to get the sum for approximately the
	same subsequence.
	Proposition 2.

Table 3 .

 3 Real parts of the identified eigenvalues using typical SSI

	0.8644 0.8644 -0.1108 -0.1108 -0.1948 0.0383
	0.0383 0.2527	0.2527	-0.0698	0.1752	0.0981

Table 4 .

 4 Therefore, the Fundamental Transition Matrix is equal to Φ(t) = e At and R = A. If the rotation velocity Ω = 0.1 rad.s -1 , the identified eigenvalues should then be close to the ones of e Aτ corresponding to the discrete state transition matrix for Ω = 0. The helicopter model is simulated for Ω = 0.1 rad.s -1 . Data points are generated over a samples number of 2000.T d , with the discrete period T d = 200. The results of the identification are reported in Table4. Identified eigenvalues at Ω = 0.1 rad.s -1 vs. True eigenvalues at Ω = 0

	Identified eigenvalues	True eigenvalues
	for Ω = 0.1	for Ω = 0
	0.5374 ± i0.7769	0.552 ± i0.7663
	0.5428 ± i0.7654	0.552 ± i0.7663
	0.8798 ± i0.4441	0.879 ± i0.4468
	0.8766 ± i0.4469	0.879 ± i0.4468
	0.8770 ± i0.4507	0.8786 ± i0.4474
	0.8759 ± i0.4489	0.8786 ± i0.4474

This work was supported by the European project FP7-NMP CP-IP 213968-2 IRIS. The authors acknowledge Leonardo Sanches for insightful discussions and comments.

Appendix A. PROOF OF PROPOSITION 2

The elements of the Hankel matrix in (18) are products between lagged outputs y k+j and y k-l . These products write:

The observation and the control matrices are periodic.

Then

Dividing by N , leads to: -l) when N goes large, where G (k-l) is the correlation between the (k -l)-th time-invariant subsequences of the state vector and the output vector

k-l+iT d converges to 0 when N goes large, because w is a white noise which is uncorrelated of any time-invariant subsequence of the output, idem for

k-l+iT d converges to 0 too (unless when j = l = 0, which is not the case herein), for the same reason of uncorrelation