Walid Trabelsi

Laurent Michel
email: laurent.michel@isae.fr

Renaud Othomene

Effects of stitching on delamination of satin weave carbon-epoxy laminates under mode I, mode II and mixed-mode I/II loadings

Keywords: Delamination, Polymer matrix composites, Stitching, mixed-mode fracture, Woven fabrics

The objective of the present study is to characterize the effect of modified chain stitching on the delamination growth under mixed-mode I/II loading conditions. Delamination toughness under mode I is experimentally determined, for unstitched and stitched laminates, by using untabbed and tabbed double cantilever beam (TDCB) tests. The effect of the reinforcing tabs on mode I toughness is investigated. Stitching improves the energy release rate (ERR) up to 4 times in mode I. Mode II delamination toughness is evaluated in end-notched flexure (ENF) tests. Different geometries of stitched specimens are tested. Crack propagation occurs without any failure of stitching yarns. The final crack length attains the mid-span or it stops before and the specimen breaks in bending. The

ERR is initially low and gradually increases with crack length to very high values.

The mixed-mode delamination behaviour is investigated using a mixed-mode bending (MMB) test. For unstitched specimens, a simple mixed-mode criterion is identified. For stitched specimens, stitching yarns do not break during 25% of mode I ratio tests and the ERR increase is relatively small compared to unstitched values. For 70% and 50% of mode I ratios, failures of yarns are observed during crack propagation and tests are able to capture correctly the effect of the stitching: it clearly improves the ERR for these two mixed modes, as much as threefold.

Introduction

Laminated composites based on carbon fibres in a polymeric matrix are very attractive for applications in structures where high stiffness-to-weight and strength-to-weight ratios are important. However, a limitation of many composite systems is their poor resistance to impact. This is a major deficiency for laminated composites.

Various techniques have therefore been developed to enhance damage tolerance by increasing the resistance to delamination growth. These techniques involve either material improvement or fibre reinforcement. Material improvement normally requires increasing the fracture toughness of the matrix resin, as delamination initiation and propagation occur mainly in the matrix material.

However, with increasing demand for large allowable design strains and reliability of aerospace structural components, material improvement alone is not sufficient.

Use of through-the-thickness reinforcement can lead to significant improvements in interlaminar strengths and delamination resistance. This can be achieved by transverse stitching [START_REF] Mouritz | A review of the effects of stitching on the inplane mechanical properties of fibre-reinforced polymer composites[END_REF] and z-pinning [START_REF] Cartie | Delamination behaviour of z-pinned laminates[END_REF][START_REF] Grassi | Prediction of stiffness and stresses in z-fibre reinforced composite laminates[END_REF] of laminated preforms or by using textile process technologies like weaving, knitting and braiding [START_REF] Mouritz | A mechanistic approach to the properties of stitched laminates[END_REF].

Stitching is well adapted to thin composite laminates while z-pinning is more efficient for thicker laminates. Even if the stitching yarn makes up only a few percent of the total fibre content of the composite, this technique greatly improves the interlaminar fracture toughness in mode I [START_REF] Sharma | Effect of through-the-thickness stitching on impact and interlaminar fracture properties of textile graphite-epoxy laminates[END_REF][START_REF] Jain | On the effect of stitching on mode I delamination toughness of laminated composites[END_REF][START_REF] Chen | A new mode I fracture test for composites with translaminar reinforcements[END_REF][START_REF] Jain | Effect of reinforcing tabs on the mode I delamination toughness of stitched CFRPs[END_REF][START_REF] Sun | Effect of stitch distribution on mode I delamination toughness of laminated DCB specimens[END_REF][START_REF] Wood | The effect of stitch distribution on mode I delamination toughness of stitched laminated composites -Experimental results and FEA simulation[END_REF] and mode II [START_REF] Jain | On the effects of stitching in CFRPs -II. mode II delamination toughness[END_REF][START_REF] Sankar | Mode II delamination toughness of stitched graphite-epoxy textile composites[END_REF][START_REF] Sharma | Effect of through-the-thickness stitching on impact and interlaminar fracture properties of textile graphite-epoxy laminates[END_REF], as well as impact damage tolerance [START_REF] Liu | Delamination resistance in stitched and unstitched composite plates subjected to impact loading[END_REF][START_REF] Sharma | Effect of through-the-thickness stitching on impact and interlaminar fracture properties of textile graphite-epoxy laminates[END_REF]. The improvement in the delamination toughness in mode I and mode II was found to depend upon the thread type, stitch type and density, stitching distribution and density. Wood et al. [START_REF] Wood | The effect of stitch distribution on mode I delamination toughness of stitched laminated composites -Experimental results and FEA simulation[END_REF] reported that stitch distribution plays an important role in determining the steady-state strain-energy release rate. Experimental results show that, for similar stitch densities, stitches improve the critical strain-energy release rate in mode I, G IC . Jain et al. [START_REF] Jain | On the effect of stitching on mode I delamination toughness of laminated composites[END_REF] have found an up to 10-time increase in mode I fracture resistance of composite specimens made using a dry carbon fabric and resin transfer moulding (RTM) technique.

Regarding mode II fracture toughness, the authors [START_REF] Jain | On the effects of stitching in CFRPs -II. mode II delamination toughness[END_REF] reported that stitches improve the toughness up to 4 times. Sharma and Sankar [START_REF] Sharma | Effect of through-the-thickness stitching on impact and interlaminar fracture properties of textile graphite-epoxy laminates[END_REF][START_REF] Sankar | Mode II delamination toughness of stitched graphite-epoxy textile composites[END_REF] reported that the use of Kevlar-2790 as the stitching yarn improved the fracture toughness by about 15 times, Glass-1250 improved it by about 30 times, and Glass-750 by about 21 times.

For the development of an entirely composite aerospace structural component, a modified chain stitching technique has been chosen to reinforce some structural details such as the skin-stiffener cap joints. The general idea is to avoid any problems with delamination in these very critical areas. Therefore, the improvement in interlaminar toughness due to this specific reinforcement has to be measured. This has to be done not only in pure modes I and II but also in mixed mode I/II. To our knowledge, the stitching effect has not yet been addressed in the literature for mixed-mode conditions. Double cantilever beam (DCB) testing is an effective way to ascertain the mode I energy release rate (ERR) for unstitched laminated composites. It also can be used to evaluate the improved delamination toughness of stitched laminates. However, experimental studies observed that the stitched specimens fail in bending (high compressive bending stresses) using DCB testing [START_REF] Jain | On the effect of stitching on mode I delamination toughness of laminated composites[END_REF][START_REF] Chen | A new mode I fracture test for composites with translaminar reinforcements[END_REF]. Therefore, a new DCB specimen configuration based on the design developed by Guenon et al. [START_REF] Guenon | Toughness properties of a three-dimensional carbon-epoxy composite[END_REF] was used by Jain et al. [START_REF] Jain | Effect of reinforcing tabs on the mode I delamination toughness of stitched CFRPs[END_REF] and Wood et al. [START_REF] Sun | Effect of stitch distribution on mode I delamination toughness of laminated DCB specimens[END_REF][START_REF] Wood | The effect of stitch distribution on mode I delamination toughness of stitched laminated composites -Experimental results and FEA simulation[END_REF]. To prevent premature failure of the cantilever substrate beams, specimens were bonded with a pair of aluminium tabs to either side along with loading blocks in order to apply a higher bending load in DCB testing. These studies show that tabbed double cantilever beam (TDCB) specimens can effectively determine interlaminar toughness of stitched laminates.

Several test configurations have been proposed for mode II delamination toughness of carbon-epoxy composite laminates. The most commonly used of these is the end-notched flexure (ENF) configuration. General observations show that with unstitched specimens the cracking is unstable and instantaneously goes to the central loading pin, whereas the crack propagation in stitched beams is steady [START_REF] Jain | On the effects of stitching in CFRPs -II. mode II delamination toughness[END_REF][START_REF] Sankar | Mode II delamination toughness of stitched graphite-epoxy textile composites[END_REF]. The crack traction forces provided by stitches are sufficient to induce stable propagation and improve the toughness.

Different tests have been proposed for mixed mode [START_REF] Benzeggagh | A mixed-mode specimen for interlaminar fracture testing[END_REF][START_REF] Rugg | Mixed-mode delamination of polymer composite laminates reinforced through the thickness by z-fibres[END_REF][START_REF] Davies | Standard test methods for delamination resistance of composite materials: Current status[END_REF], but some of them have limitations and do not permit testing of material for a broad range of mixedmode ratios. The Standard ASTM Mixed-Mode Bending (MMB) test established for unidirectional carbon fibre tape laminates appears to be the best choice. The main advantage of this test method is the ability to obtain different mixed-mode ratios, ranging from pure mode I to pure mode II [START_REF] Reeder | Non-linear analysis and redesign of the mixed-mode bending delamination test[END_REF][START_REF] Reeder | A mixed-mode bending apparatus for delamination testing[END_REF]. This method can be used effectively to characterize mixed-mode interlaminar fracture of woven laminate specimens.

Consequently, the purpose of this study is, first, to examine the effect of stitching on mode I toughness. Unstitched carbon-epoxy laminate samples were tested under mode I using DCB and TDCB specimens to evaluate effects of the reinforcing tabs on the determination of mode I toughness. Then the interlaminar fracture toughness is investigated and evaluated for TDCB stitched specimens.

Second, the mode II fracture toughness is measured by subjecting unstitched and stitched specimens to end-notched flexure (ENF) tests. Third, the mixed-mode delamination behaviour of through-thickness reinforcement laminates is investigated using a mixed-mode bending (MMB) specimen. The fracture mechanism is analyzed and the energy required for the initiation and growth of an artificial crack is calculated.

To complement visual observations, acoustic emission is used here to monitor damage accumulation before and during crack propagation and helps to define crack initiation. As reported by Sharma and Sankar [START_REF] Sankar | Mode II delamination toughness of stitched graphite-epoxy textile composites[END_REF], the crack surfaces do not open during the ENF test, and it is difficult to control the crack length during the experiment. X-ray imaging is used to estimate crack length after tests.

Experimental results are then gathered for unstitched and stitched specimens. An empirical criterion is identified for the unstitched composite fracture behaviour.

Effects of stitching on these laminates are obtained and limits of this experimental study are given.

Materials and experiments

Material systems

In this study, the composite specimens were manufactured by Latécoère from carbon woven-fabric 5H satin (Hexforce ® G0926) and epoxy resin (Hexflow ® RTM 6) using a resin transfer moulding (RTM) technique. A 13-µm highperformance release Teflon film was inserted at the mid-plane of the dry fabric laminates to simulate the crack defect. Two different kinds of quasi-isotropic laminates were manufactured: one 3 mm thick (8 plies) and the other 6 mm thick (16 plies).

The dry 2D laminates were stitched, before moulding, with Kevlar ® stitching yarns. The type of stitch chosen for this study is the modified chain stitch ® (Figure 1). This process is easy to implement requiring only one side access to achieve stitching (one-side stitching). The stitching density chosen for this study is 2.5 st/cm². After consolidation, the plates were subjected to an ultrasonic C-scan to check quality and location of Teflon inserts. Following this process, the specimens were further cut and measured to the necessary dimensions for each type of test.

Specimens and testing procedures

Table I summarizes dimensions and specific features of all specimens tested in the study. For stitched specimens, the width was increased (from 25 mm to 35 mm) in order to get a line of 4 stitch points across the width and thus obtain the stitching density required for this study. The thickness values in this table define specimens 3 mm thick (8 plies) and 6 mm thick (16 plies).

Mode I delamination

ASTM D5528 [START_REF] Astm D | Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites[END_REF] describes the experimental determination of the opening mode I interlaminar fracture toughness, G IC . This method is used to ascertain the interlaminar fracture toughness of unstitched carbon-epoxy laminates. The crack length (a 0 ) is determined from the loading point to the crack tip (Figure 2). On this side, loading blocks are bonded to either face of the specimen.

For the TDCB tests, aluminium tabs, 10 mm thick, are bonded to either side of both stitched and unstitched specimens using adhesive film (Figure 2). Loading blocks are directly machined on the tabs.

Mode II delamination

Mode II tests were performed using a three-point bending test. The main advantage of the end-notched flexure (ENF) fixture is the ability to adjust the total span length (2L) of specimens (Figure 3).

The unstitched specimens were positioned in the fixture with a total span of 100 mm, and an initial crack length ranging from 15 mm to 35 mm.

For the stitched laminated composite specimens, the total span length varies between 130 mm and 190 mm for the two thicknesses (3 mm and 6 mm) studied here. The initial crack length depends on the span length and varies between 30 mm and 55 mm.

Mixed-mode delamination

The mixed-mode bending (MMB) test method is adopted for this study. It is a simple combination of the DCB (mode I) and ENF (mode II) specimen (Figure 4).

This standard ASTM test [START_REF] Astm D | Standard Test Method for Mixed Mode I-Mode II Interlaminar Fracture Toughness of Unidirectional Fiber Reinforced Polymer Matrix Composites[END_REF] is used to determine the strain-energy release rate (SERR) at various ratios of mode I and mode II loading. Three modal ratios were chosen (G I /G=25%, G I /G=50% and G I /G=70%). The initial crack length (a 0 ) ranged from 15 mm to 35 mm. Crack propagation is observed with the aid of an Olympus DP20 digital camera.

The system includes an Olympus SZX10 stereo microscope (×1 to ×5 magnification) able to capture a high-digital-resolution image in real time during the test.

During these tests, the crack length is measured by visual observation with the aid of white marks along the specimen edge for unstitched and stitched laminates.

Failure of stitch yarns is easily observed during mode I tests. In contrast, for mode II and mixed mode it is very difficult to conclude from visual observations. So, post-test X-ray imaging was used to determine stitch failures. Indeed, failed yarns absorb X-ray contrast agent and clearly appear on X-ray images (see Figure 16).

Analysis methods

The DCB test is used to determine the interlaminar fracture energy toughness, G IC , for unstitched laminates. A load P is applied to each arm of a specimen with a delamination length a. The deflection of the specimen is measured by the displacement of the crosshead, δ. The energy rate for mode I, G I , is calculated using the Modified Beam Theory (Eq. 1).

( )

∆ + = a b P I G 2 3 δ Eq. 1
where b is the specimen width and ∆ is the effective delamination extension to correct for rotation of DCB arms at the delamination front. It is determined experimentally by generating a least-squares plot of the cube root of compliance,

C 1/3
, as a function of delamination length (Figure 5). The compliance, C, is the ratio of the load point displacement to the applied load, δ/P. For the TDCB specimens, the energy release rate is calculated using the Dransfield et al. formulation [START_REF] Jain | Effect of reinforcing tabs on the mode I delamination toughness of stitched CFRPs[END_REF]. Assuming that each cracked substrate beam behaves as a cantilever of length a (the crack length), the deflection, δ (crackopening displacement), may be estimated by summing a bending deflection term with a shear deflection term as follows: The energy rate for TDCB test, G I , is defined by Griffith's energetic fracture criterion (Irwin-kies equation) [START_REF] Irwin | Analysis of stresses and strains near the end of crack traversing a plate[END_REF], it is calculated as follows:

( ) ( ) ( ) ( ) b h h G h h a P EI h h a P
( ) ( )       + + + + = b h h G EI h h a b P G a a I ) ( 5 6 12 2 2 χ Eq. 3
The expression for EI may be written as follows:

( )         - + = ∑ ∑ 2 3 12 i i i i i c y y A h b E EI Eq. 4                       +       +         + = = ∑ ∑ a C a c a C a i i i i i i i h E E b bh h h E E h b h h b E A y E A y 2 2
Eq. 5

Where y is the distance from the designated origin to the neutral axis (Eq. 5), y is the distance of the mid-point of each component from the same origin and A is the cross-sectional area of each component. Thus the expression for flexural rigidity may be simplified as :

              - +         +         +       - + = 2 3 2 3 2 12 2 12 y h h h E E b h E E b h y h b bh E EI a a C a a C a C
Eq. 6

The interlaminar shear modulus of the system, G 12 , is given by the inverse rule of mixtures (Eq. 7).

( )

        + + = a a c a G h G h h h G 12 12
Eq. 7

where G a is the shear modulus of the aluminium tabs.

For unstitched and stitched ENF specimens, the beam theory formula is used to calculate the critical strain-energy release rate in mode II.

( )

3 3 2 2 3 2 2 9 a L b a CP G II + = Eq. 8
For MMB specimens, the existing literature [START_REF] Reeder | Non-linear analysis and redesign of the mixed-mode bending delamination test[END_REF][START_REF] Kinloch | The mixed-mode delamination of fiber composite materials[END_REF] uses different methods to calculate the total mixed-mode strain-energy release rate, G. The MMB standard test method [START_REF] Kinloch | The mixed-mode delamination of fiber composite materials[END_REF] describes the calculation of the G I and G II with the Modified Beam Theory method. However these methods are designed for unidirectional continuous fibre tape laminates and their application in our case posed some difficulties, particularly in identifying a consistent bending modulus.

Therefore, for our stitched and unstitched multidirectional fabric laminates, the experimental compliance method already used for fabric laminates [START_REF] Benzeggagh | Measurement of Mixed-Mode Delamination Fracture Toughness of Unidirectional Glass/Epoxy Composites with Mixed-Mode Bending Apparatus[END_REF] is chosen to compute the energy release rate, G. This method considers that the mode ratio is independent of the delamination length and the crack propagation is similar to mode II.

=

Eq. 9

where m is the slope of the plot compliance (C=δ/P) versus a 3 . It depends on the mode ratio given by Eq 9.

where c is the lever arm length. G I and G II can be calculated using the following expressions:

1 + = Υ Υ G G I
Eq. 11

Υ + = 1 G G II Eq. 12 Υ =       + - ⋅ = 2 3 3 4 L c L c G G II I For L c ≥ 3
Eq. 10

Testing results and discussion

Unstitched and stitched mode I results

The results of DCB and TDCB test for unstitched specimens are shown in Figure 6. It depicts the load versus displacement for the two thicknesses (3 mm and 6 mm) studied here. Acoustic emission monitoring reveals intense acoustic activity for all specimens before the load decrease, and starting from the point at which the loaddisplacement curves become non-linear. The starting point of crack propagation is generally observed close to the top of (F/δ) curves. For DCB specimens, the load tends to decrease during crack extension but the propagation is unstable. Sudden load drops due to crack progress are observed and followed by slow recoveries of the load. It has been visually observed that such sudden crack progress is associated with bridging effects due to the fabric pattern. For TDCB specimens, the aluminium tabs reduce the effects of the fibre bridging and the load decreases more steadily and slowly.

Concerning the thickness effect (from 3 mm to 6 mm), the peak load increases with the thickness of DCB unstitched specimens. The reason for this is the increased stiffness, while for TDCB specimens increasing the composite thickness does not affect the load-displacement curve. The tabs increase the bending stiffness of the specimen, thus hiding the stiffness increase due to a thicker laminate.

The curves of mode I delamination toughness (G I ) versus crack length (a) for the unstitched specimens (3 mm and 6 mm thick) are shown in Figure 7. The composite toughness curves reveal a linear increase in G I with the crack growth up to a point where the evolution flattens. A certain crack extension, called the process zone by Jain et al. [START_REF] Jain | Effect of reinforcing tabs on the mode I delamination toughness of stitched CFRPs[END_REF], is needed to obtain stable propagation.

The value of G I at the plateau is referred to as the critical propagation toughness, G IC . The starting point of the plateau is shown by a vertical line (Figure 7). It can be clearly observed from the curves that the process zone of the unstitched 3-mmthick specimens increased from 40 mm to 90 mm and from 60 mm to 70 mm for the 6-mm-thick laminates.

Experimental results for unstitched specimens, both tabbed and untabbed, exhibit a good agreement for the G I values. For 3-mm-thick DCB and TDCB specimens, G IC is 803 and 853 J/m². The value at the plateau is 1017 J/m² (DCB) and 1024 J/m² (TDCB) for 6-mm-thick laminates.

These experimental data show that adding aluminium tabs causes an increase of the crack length needed to attain a steady-state behaviour and the critical propagation toughness, G IC .

For the stitched specimens, there is an intense acoustic emission as soon as the load-displacement curve becomes non-linear (Figure 8). At this point the crack propagates quickly before reaching the first stitching yarns. Then its propagation slows down but goes on steadily. More stitching points become active in the crack wake until the first one fails at a local opening displacement of around 1.5 mm (Figure 9). This failure is accompanied by a substantial load drop (Figure 8). As the crack front reaches the first row of stitches, the R-curve slope rises rapidly with the increase of bridging entities and G I increases up to a critical value of length propagation. G IC is achieved after 81 mm of crack propagation with a value of 3500 J/m² (Figure 8). At this length of propagation, 8 rows of two stitches are active in the crack wake; this is confirmed by visual observations (Figure 9). The scatter of data points after this point can be attributed to the fracture behaviour of the stitch fibres. In comparison to unstitched specimens, it is immediately obvious that the addition of the through-the-thickness reinforcement by stitching improves the interlaminar toughness in mode I as much as fourfold, from 853 to 3500 J/m².

Unstitched and stitched mode II results

For the unstitched specimens, as expected the crack propagation is unstable and leads to a sharp drop of the load. Close to this point, acoustic emission monitoring reveals intense activity (Figure 10).

For stitched specimens, depending upon the specimen geometry (L, h, a 0 ), either the crack propagates up to the mid-span of the specimen (Figure 10: c-d The maximum values of ERR, determined when the crack reaches the mid-span length, are higher than experimental data obtained with unstitched specimens.

Therefore, the stitches increase the ERR of woven laminates but much less than reported in the literature [START_REF] Jain | On the effects of stitching in CFRPs -II. mode II delamination toughness[END_REF][START_REF] Sankar | Mode II delamination toughness of stitched graphite-epoxy textile composites[END_REF], where the increase reached up to 15 times the fracture energy of unstitched specimens. It is important to note that there is no plateau in the G II R-curve where the G II is stable as is the case for mode I.

Even if the required energy for crack propagation is larger for stitched specimens, this additional energy is used to deform the stitches and is not high enough to break them. Indeed, X-ray images (Figure 13) obtained after test completion show that the stitches used in this study did not break. After testing, the cracked part of some specimens was cut. And it was impossible to separate by hand the two arms of these samples, showing that the stitching points remain safe. This also confirms visual observations made during the tests. To solve this problem and to obtain failure of stitches during experimental tests, changes were made in the geometry of specimens. The thickness (from 3 mm to 6 mm) and the mid-span length (up to L=95mm) were increased (see Figure 12, ENF3_6mm a 0 =45, L=95, b=35). However experimental tests showed that stitching yarns still did not break for this specimen. The maximum relative displacement of the two lips of the crack was of 1.05 mm (Figure 14) at the location of the first stitching row. That means this stitching is able to sustain a large extension before failure as it was observed for mode I. So a few more tests were performed with 3mm thick tabs bonded on each side of 3mm-thick specimen and here again crack propagations without any stitching failures were observed. This confirmed that under certain test configuration (i.e large span between supports) it is possible to obtain for the material under study steady crack propagation without any failure of stitching points.

Comparing to the standard results observed in the literature this behaviour has to be addressed. This might be due to a larger elongation capability of the stitching yarn or to the specific type of stitching: one side modified chain stitching or to a lack of tension of the yarn during the stitching process.

Unstitched and stitched mixed-mode I/II results

The load-versus-displacement curves for three MMB unstitched specimens are shown in Figure 15. The experiments presented in this figure correspond to the three modal ratios studied (G I /G=25%, G I /G=50% and G I /G=70%). For all the curves, the load increases in a single stage before progressively decreasing.

Observations and X-ray images show that the crack reaches the mid-span of the beam for all the unstitched specimens tested.

At G I /G=70%, the slope of the decrease flattens. The propagation is unsteady as it was observed for the pure mode I (Figure 6). When increasing mode II ratio the load decrease becomes more and more pronounced after damage initiation and the crack propagation is more and more unsteady.

Observations and X-ray images show that the crack reaches the mid-span of the beam for all the unstitched specimens tested. For stitched MMB specimens, the load-versus-displacement curves are shown in Figure 16 for the three modal ratios studied (25%, 50% and 70% of mode I). For G I /G=25% case, the load-displacement curve is almost linear. Then the load remains almost constant while the crack propagates quickly up to the mid-span of the beam. Observations and X-ray images show that stitching yarns did not break during this crack propagation.

In the case of the G I /G=50% MMB specimens test, the load-displacement curves show that the load decreases after damage initiation. At the end of the test, the average value of the crack length is 42 mm. The specimen fails in bending before the crack reaches the mid-span length. X-ray images show eight failure points of the stitching yarn.

For the G I /G=70% case, the final crack length presents an average value of 33 mm and the specimen fails in bending. X-ray images show that six stitching yarns are broken. For the cases of G I /G=70% and G I /G=50%, results show that the stitching improves the fracture delamination toughness. The increase of G is three times higher than with the unstitched specimens. For G I /G=25%, the value of G determined for stitched MMB specimen is only slightly higher than the fracture energy of unstitched specimens. The effect of stitching did not appear on experimental values for this last modal ratio. The reason is that the stitching yarns did not break. to define the fracture. The most widely used criterion to predict delamination propagation under mixed-mode loading, the power-law criterion, is established in terms of an interaction between the energy-release rates [START_REF] Whitcomb | Analysis of instability-related growth of a through-width delamination[END_REF]. It contains the G I /G IC and G II /G IIC terms (Eq. 13).

The power-law criterion obtained with α=1 and β=1 was found suitable to predict failure of composite laminates. This criterion simply normalizes each component of fracture toughness by its pure mode values [START_REF] Wu | Crack extension in fiberglass-reinforced plastics[END_REF].

The energy-release rate of 3-mm-thick unstitched specimens from experimental data is plotted on a GI-versus-GII chart as shown in Figure 17. For stitched specimens, experimental results determined on TDCB, ENF and MMB tests are plotted on the G I -versus-G II chart (Figure 18). For pure mode I, G I /G=70% and G I /G=50% ratios, good values of fracture energy are found.

Failures of yarns are observed during the propagation of the crack. Stitching clearly improves the fracture delamination toughness, as can be seen in Figure 18.

On the other hand, for the stitched ENF and 25% MMB specimens, the experimental results do not show any stitching yarn failures. As a consequence,

1 =         +         β α IIC II IC I G G G G Eq. 13 1 = + IIC II IC I G G G G
Eq. 14 the effect of stitching remains limited and improvement of ERR is rather small compared to the other ratios. So, at this point of the study it is not possible to provide a consistent mixed-mode fracture criterion for these stitched laminates including the fully role of the 3D reinforcement.

However what must be highlighted is that for mode II dominated propagation there exist the possibility to get large crack propagation without any failure of the stitching. This might be a restraint to the application of this technique to aeronautical structures where it is important to assess none propagation of the delamination at limit loads. For aeronautical structures, this stitching will play the role it is expected to have when mode I is critical. When mode II becomes predominant, even if the stitching doesn't improve largely the energy necessary to propagate the crack, it will likely lead to a stable propagation.

Conclusions

  strain-energy release rate in mixed mode G I : Strain-energy release rate in mode I G II : Strain-energy release rate in mode II G IC : Interlaminar fracture toughness in mode I G IIC : Interlaminar fracture toughness in mode II C: Compliance P: Applied load δ : Load point displacement b: Specimen width a 0 : Initial crack length a: Crack length h: Thickness of specimen arm L: Specimen mid-span length ∆: Effective delamination extension to correct for rotation of DCB arms at delamination front E: Flexure modulus EI: Flexure stiffness G 12 : Transverse shear modulus G 12c : Transverse shear modulus of composite laminates G a : Transverse shear modulus of the aluminium tab χ: Correction factor for TDCB specimens m: Slope of the plot compliance (C=δ/P) versus a 3 c: Lever arm length in mixed-mode test Y: Mixed-mode partition (G I /G II ) α, β : Material parameters used in power-law criterion

Figure 1 .

 1 Figure 1. Modified chain stitch.

Figure 2 .

 2 Figure 2. Schematic representation of (a) DCB specimen and (b) TDCB specimen for mode I tests.

Figure 3 .

 3 Figure 3. Schematic representation of ENF specimen for mode II test.

Figure 4 .

 4 Figure 4. Schematic representation of MMB specimen for mixed-mode test.

3. 3 .

 3 Test conditions All tests were performed in a testing machine equipped with two load cells (1 kN for DCB and MMB tests and 25 kN for TDCB and ENF tests). The specimen is pin-loaded at constant displacement rates of 2 mm/min and then 0.5 mm/min to ensure slow propagation. The load and displacement are measured by reading signals from a load cell and a linear variable displacement transformer (LVDT) into a digital data logger used to record acoustic emission data during the test.

Figure 5 .

 5 Figure 5. Correction factor for the modified beam theory (ASTM [20]).

  is the thickness of the aluminium tab, EI is the flexure stiffness; G 12 is the transverse shear modulus of the aluminium-composite beam and χ is a correction factor which takes into account the additional deflection due to the crack tip rotation of the beams. It can be determined from DCB and TDCB experimentally results of unstitched specimens.

Figure 6 .

 6 Figure 6. Load (P) and amplitude (acoustic emission) versus displacement (δ ) for unstitched specimens (a) 3-mm thick DCB (b) 3-mm thick TDCB (c) 6-mm thick DCB (d) 6-mm thick TDCB.

Figure 7 .

 7 Figure 7. G I versus crack length (a) for unstitched specimens (a) 3-mm thick DCB (b) 3-mm thick TDCB (c) 6-mm thick DCB (d) 6-mm thick TDCB.

Figure 8 .

 8 Figure 8. (a) Load (P) and amplitude (acoustic emission) versus displacement (δ ) for TDCB3_3mm stitched specimen (b) G I versus crack length (a) for 3-mm-thick TDCB stitched specimens.

Figure 9 .

 9 Figure 9. Photos displaying crack wake for 3-mm-thick stitched TDCB specimens.

Figure 10 .

 10 Figure 10. Load (P) and amplitude (acoustic emission) versus displacement (δ ) for ENF specimens (a) 3-mm-thick (a 0 =17, L=50, b=25) unstitched specimens (b) 6-mm-thick (a 0 =15, L=50, b=25) unstitched specimens (c) 3-mm-thick (ENF4_3mm a 0 =50, L=70, b=35) stitched specimens (d) 6-mm-thick (ENF2_6mm a 0 =35, L=75, b=35) stitched specimens.

Figure 11

 11 Figure11shows the G II R-curves obtained for the unstitched specimens. As observed for mode I, G II increases slightly by increasing the thickness (from 3 mm

Figure 11 .

 11 Figure 11. G II versus initial crack length (a 0 ) for ENF unstitched specimens.

Figure 12 .

 12 Figure 12. G II versus crack length (a) for ENF stitched specimens (3 mm and 6 mm thick).

Figure 13 .

 13 Figure 13. X-ray images and observations of ENF stitched specimens (a) 3 mm thick, (b) 6 mm thick.

Figure 14 .

 14 Figure 14. Relative displacement of the two lips of the crack at the location of the first stitching row for 6-mm-thick stitched ENF specimens.

Figure 15 .

 15 Figure 15. Load-versus-displacement curves for G I /G=70%, 50% and 25% (a 0 =25) for 3-mm-thick MMB unstitched specimens.

Figure 17

 17 Figure 17 presents X-ray images and test observations of MMB stitched specimens.

Figure 16 .

 16 Figure 16. Load-versus-displacement curves for G I /G=70%, 50% and 25% (a 0 =25) for 3-mm-thick MMB stitched specimens.

Figure 17 .

 17 Figure 17. X-ray images and observations of 3-mm-thick MMB stitched specimens (a) G I /G=70% (b) G I /G=50% (c) G I /G=25%.

Figure 18 .Figure 18 .

 1818 Figure 18. Mode I energy contribution, G I (Average values), versus mode I energy component G II (Average values), for 3-mm-thick unstitched and stitched specimens.

  For mode I, the experimental results for DCB and TDCB unstitched specimens are in a good agreement for the G IC values. Adding tabs increases the crack length needed to attain steady-state behaviour enabling measurement of the critical propagation toughness. For stitched specimens, the crack propagates quickly before reaching the first stitching yarns, and then after that point slowly and steadily. Yarn failure is associated with a substantial load drop on the load-versusdisplacement curve. Experimental results in mode I show that stitching generates an up-to-fourfold improvement of the energy release rate, from 853 to 3500 J/m². G IC is achieved after 81 mm of crack propagation. At this point, eight rows of two stitches are active in the process zone, with an opening relative displacement up to 1.5 mm.As expected, the crack propagation of ENF unstitched specimens is unstable, i.e. associated with a sharp load drop. The ERR is stable in mode II. For stitched specimens, either the final crack length attains the mid-span or it stops before and the specimen breaks in bending. Observations and X-ray images show that the stitches do not break even if the crack propagates up to the mid-span length and the relative shear displacement of the 1 st stitching row climbs up to 1mm. This specific behaviour may have several sources including type of yarn, of stitching point and even a problem of process and has to be addressed in the future. TheERR is initially low and gradually increases with crack length to very high values.There are no plateaus in the G II R-curve and the maximum values are higher than experimental data obtained with unstitched specimens.Concerning mixed-mode propagation of unstitched laminates a simple linear criterion has been identified. For stitched laminates, the first stage of the loaddisplacement curve of the G I /G=25% mode I ratio is almost linear. Then the load remains almost constant while the crack propagates quickly up to the mid-span of the beam without any failure of stitching yarns. The improvement in ERR is very slight compared to the values of unstitched specimens. For G I /G=70% and G I /G=50% mode I ratios, the load-displacement curves show that the load decreases after damage initiation. Observations and X-ray images reveal that the crack propagates and the stitching yarns are broken before the specimen fails in bending. For these two mixed-mode ratios the stitching improves the ERR up to a factor of three. 25% MMB and ENF tests reveal no failure of stitches, which disables a fullyeffect of stitching on ERR values to be obtained. This clearly reveals the fact that crack propagation is not stopped by stitching in mode II dominated delamination behaviour. On the contrary for TDCB, 70% and 50% MMB where the mode I ratio is higher crack propagation occurs with stitching failures. The reinforcement role of stitching is then fully achieved as measured by the large improvement of ERR values, up to threefold.

Table I .

 I Test specimens.

		Thickness	Length	Width	Initial crack
		(mm)	L (mm)	b (mm)	a 0 (mm)
	Unstitched DCB specimens	3 and 6	150	25	30 to 50
	Unstitched
	TDCB	3 and 6	150	25	30 to 50
	specimens
	Stitched TDCB specimens	3	150	35	30 to 40
	Unstitched ENF specimens	3 and 6	100	25	15 to 35
	Stitched ENF specimens	3 and 6	130 to 190	35	30 to 55
	Unstitched MMB specimens	3	100	25	15 to 35
	Stitched MMB specimens	3	100	35	15 to 35

  Table II lists the experimental results for MMB stitched and unstitched specimens.

Table II .

 II Total strain-energy release rate in mixed mode, G, for unstitched and stitched MMB specimens.

	Unstitched
	specimens
	G I /G	G	G I /G	G	G I /G	G
	(%)	(J/m²)	(%)	(J/m²)	(%)	(J/m²)
	70%	1340	50%	1600	25%	1510
	70%	1000	50%	1020	25%	1060
	70%	880	50%	900	25%	1290
	70%	1000	50%	810	25%	1230
	70%	950	50%	1070	25%	1450
	70%	780	50%	1230	25%	1890
	Average	990±190	Average	1100±280	Average	1410±290
	Stitched
	specimens
	G I /G	G	G I /G	G	G I /G	G
	(%)	(J/m²)	(%)	(J/m²)	(%)	(J/m²)
	70%	4530	50%	3190	25%	2200
	70%	3020	50%	3580	25%	1720
	70%	4100	50%	4590	25%	1950
	70%	1360	50%	4040	25%	1970
	70%	3940	50%	4310	25%	3520
	70%	2580			25%	1720
	Average	3640±810	Average	4020±540	Average	2120±680
	2.4

. Mixed-mode delamination criteria For

  unstitched and stitched composite laminates, it is interesting to know how the mode I energy contribution varies versus the mode II energy component. Many delamination fracture criteria were developed before obtaining consistent sets of mixed-mode experimental data. Early representations of the mixed-mode response were generally made by plotting G I and G II , data on a Cartesian coordinate system

Acknowledgment

This work was supported by Latécoère and Région Midi-Pyrénées.


Pierre Selva
email: pselva@mit.edu

Joseph Morlier

Yves Gourinat

Toward a three-dimensional finite-element model of the human inner ear angular accelerometers sensors

Keywords: Finite-element model, fluid structure interaction, inner ear, semicircular canals

The vestibular section of the inner ear is partly comprised of three semicircular canals that detect head angular acceleration. This system is filled with a Newtonian fluid and its functioning pertains on fluid-structure interactions. The aim of the present paper is first to present a two-dimensional model of a single semicircular canal using fluid-structural finiteelements simulations. Second, this two-dimensional model is extended to 3-D space and to the case where the entire set of canals is considered. To achieve this goal, we first develop a two dimensional finite-elements model of a single canal. Using a strong coupling between the fluid flow and the structural displacements and also an Arbitrary Lagrangian Eulerian (ALE) approach for the moving mesh, we analyze displacements of the cupulae and fluid velocity during head rotation. Second, this 2-D model is extended to a three-dimensional case by considering the entire set of canals. Preliminary results showing cupula deformation as well as fluid flow highlight a good correlation with the 2-D model.

Introduction

Daily human activity includes complex orientation, postural control, and movement coordination. All these tasks depend upon his perception of motion. The non-auditory section of the human inner ear, the vestibular system, is recognized as the prime motion sensing center. It represents an inertial measuring device which allows us to sense, in the absence of external sensory cues (vision, etc) self motion with respect to the six degree of freedom in space. The vestibular system, which constitutes a very small fluid-filled system the size of a pea, is comprised of two primary sense organs: three semicircular canals which detect head angular acceleration and two otolith organs which respond to linear accelerations of the head and to gravity. In particular, when our head experience a movement of rotation, the fluid in the canal lags behind due to its inertia and produces a force across a thick diaphragm called the cupula. Therefore this diaphragm, which completely occludes the canal, experiences a deflection in the opposite direction of head movement. The cupula covers a transverse ridge of sensory epithelium called the crista that contains thousand of sensory hair cells. When the cupula is deflected, hair cell sensory cilia bend as well. This initiates a complex transduction process in hair cells and vestibular afferent neurons, and a nervous signal is finally transmitted to the brain so as to provide a sensation of motion.

In the present paper, we are concerned with the modeling of the entire set of semicircular canals. Since the 30's, numerous models of the semicircular canal macromechanics have been suggested using different approaches. W. [START_REF] Steinhausen | Uber die Beobachtung der Cupula in den bogengansampullen des Lebenden Hechts[END_REF] formulated a classical torsion pendulum model for the dynamic behavior of a single SCC. This model, which has been the benchmark for subsequent works, consists of a single-degree of freedom overdamped springmass-damper system subject to mass-proportional inertia forcing. Several notable extensions have then been made to enhance this original model by relating the geometry and structure of the SCC to mass, stiffness, and damping parameters appearing in the model (e.g. [START_REF] Van Egmond | The mechanics of the semicircular canal[END_REF][START_REF] Groen | The function of the vestibular organ[END_REF][START_REF] Van Buskirk | The fluid mechanics of the semicircular canals[END_REF][START_REF] Oman | The influence of semicircular canal morphology on endolymph flow dynamics[END_REF][START_REF] Rabbitt | Biomechanics of the semicircular canals and otolith organs[END_REF]. Other models were based on the resolution of the fluid flow equation within the canal [START_REF] Van Buskirk | The effect of the Utricle on Fluid Flow in the Semicircular Canal[END_REF][START_REF] Van Buskirk | The biomechanics of the semicircular canals[END_REF][START_REF] Steer | The influence of angular and linear acceleration and thermal stimulation on the human semicircular canal[END_REF][START_REF] Oman | The influence of semicircular canal morphology on endolymph flow dynamics[END_REF][START_REF] Damiano | A singular perturbation model of fluid dynamics in the vestibular semicircular canal and ampulla[END_REF]Rabbitt et al., 1999). The three-dimensional model of [START_REF] Oman | The influence of semicircular canal morphology on endolymph flow dynamics[END_REF], in which the non-uniform geometry of a single canal was considered, probably constitutes the most compatible biophysical single-degree of freedom model of the SCC. Among all these models, some are formulated in one dimension, while few of them consider a three-dimensional geometry. Moreover, all of these models consider a single canal and most of them do not take into account the fluid-structure interaction but rather consider the influence of the cupula by a punctual elasticity.

The aim of this paper is to model the entire set of semircircular canals by taking into account fluid-structure interactions in order to investigate fluid flow and cupula motion during head rotation. To achieve this goal, we use the finite-element Comsol Multiphysics software as it permits to deal with different physics, and more particularly with FSI problems. The modeling strategy is as follows: first, we model a two-dimensional cross-section of the lateral semicircular canal using geometry and dimensions extracted from measured human data by [START_REF] Curthoys | Dimensions of the horizontal semicircular duct, ampulla, and utricle in the human[END_REF]; second, we extend this 2D model to a three-dimensional model of a single semicircular canal; and third, we develop a fully three-dimensional model by considering the three SCCs.

2D model

Geometry of the 2D model

We here consider a two-dimensional cross-section of the lateral semicircular canal. The geometry and all the associated dimensions (Fig. 2.12) are extracted from measured human data by [START_REF] Curthoys | Dimensions of the horizontal semicircular duct, ampulla, and utricle in the human[END_REF]. The canal consists of three main regions: the semicircular canal, the ampulla, and the utricle. The canal is filled of a water-like fluid, known as endolymph. The model also considers the cupula (solid) located in the ampulla which completely seals the canal.

Governing equations

We now present the governing equations describing the behavior of the cupula and of the endolymph. Typically, these equations are solved for displacement and for velocity and pressure respectively. We model the endolymph as an incompressible Newtonian fluid. The typical governing equations of fluid flow are described in terms of the following two dimensional Navier-Stokes equations:

( )

f f f f f f U U U F t ρ ρ σ ∂ + ⋅∇ = + ∇ ⋅ ∂ (1) 0 f U ∇ ⋅ = (2)
where f ρ is the fluid density, F is the volume force affecting the fluid, f σ is the Cauchy stress tensor, ( , )

f f f U u v =
is the velocity field, and p is the pressure. The first equation is the momentum transport equation, and the second is the equation of continuity for incompressible fluids. These equations describe how the velocity, pressure, and density of a moving fluid are related.

The structural deformations of the cupula are solved using an elastic formulation and a nonlinear geometry formulation to allow large deformations that may occur due to its very low stiffness. Neglecting body forces, the Navier equation of motion for the cupula can be written in terms of displacement vector ( , )

s s s U u v = as: 2 ² s s s U t ρ σ ∂ = ∇ ⋅ ∂ (3)
We consider an elastic isotropic material for the cupula so that the generalized Hooke's law is obtained:

2 ( ) s tr I σ µε λ ε = + (4)
where µ and λ are Lame's coefficients, which are related to Young's modulus of elasticity and Poisson's ratio ν , by the following equations:

2(1 )

E µ ν = + (5) (1 )(1 2 ) E ν λ ν ν = + - (6)
Regarding the question of exchanging information, typically in the form of boundary conditions at the interface, we adopt a two-way (or fully) coupled FSI. In this case, the response of the solid is strongly affected by the response of the fluid, and vice versa. In other words, fluid flow causes deformation of the structure. This deformation, in turn, changes the boundary conditions for the fluid flow. In addition, due to this coupling, conditions are needed to ensure that the fluid and structural domains will not detach or overlap during the motion (Fig. 2.13). For a viscous fluid, the coupling between fluid and structure requires that velocities coincide along the interface. In particular, the time derivatives of the structural displacements define the fluid's velocity so that:

s f dU U dt = (7)
In addition, the force exerted by the fluid on the solid boundary must be considered and is the negative reaction force on the fluid given by:

( ) ( ) T f f f n pI U U η   = -⋅ -+ ∇ + ∇     ( 8
)
where n is the outward normal vector to the boundary. This load represents a sum of pressure and viscous forces.

Boundary conditions

The aim of the model is to investigate fluid flow dynamics as well as cupula motion during head rotations. Therefore we apply a step in angular velocity to the semicircular canal. However, in order to avoid problems with discontinuity and to improve numerical reliability as well as convergence, we use smoothed Heaviside functions with a continuous second derivative that emulates the step of angular velocity.

Regarding the boundary conditions of the solid domain (cupula), we prescribe a rotational displacement at the top and bottom of the cupula and the fluid load is applied on the wetted surfaces of the cupula. As regard the fluid domain, we impose an angular velocity at the wall of the canal. Consequently the fluid in contact with the wall of the duct rotates at the same angular velocity of the canal, which is locally equivalent to a no-slip condition. We also consider an open boundary at the common crux outlet so the nodal velocities are left free to accommodate inflow and outflow of endolymph between the horizontal semicircular canal and the other canals. Finally, at the interface endolymph/cupula we apply the structural velocity of the cupula as well as the fluid load.

Moving mesh

In order to model effectively the FSI problem we consider an Arbitrary Lagrangian-Eulerian (ALE) method which is a common application in engineering used to solve problems pertaining to structure and fluid mechanics analysis [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF]. The ALE method employs the use of reference frames to represent the classical Lagrangian and Eulerian systems. The Lagrangian reference frame is used to study the structure problem while the Eulerian reference is used to study the fluid problem. Hence, the model combines the fluid flow with structural mechanics by using a moving mesh to make sure the fluid flow is deformed along with the cupula.

The model is divided into different subdomains so as to specify precisely how the mesh displacement is computed (Fig. 2). The imposed conditions of mesh displacement are as follows:

• Solid domain (cupula): the displacements ( , )

s s u v provided by the computation of the structure mechanics equations are imposed. In other words, a Lagrangian method is used where the mesh movement follows the material motion.

• Fluid domain near the cupula: as cupula deformation may affect fluid flow in its vicinity we define two subdomains on both sides of the cupula where the mesh is free to move. This means that the mesh is constrained only by the boundary conditions on the surrounding boundaries.

• Rest of the fluid domain: the displacement of rotation imposed to the whole structure is also applied to this subdomain. In other words, the mesh is not deformed and follows the rotation of the canal.

Simulations

Numerical solutions of the governing system of coupled nonlinear system partial differential equations (PDEs) are generated using finite-element analysis software (Comsol Multiphysics 3.5a) on a 2-D mesh with 2663 quadratics triangular elements that represent 16985 degrees of freedom. We consider a rotation of the structure at a constant angular velocity of / 2 π rad/s around a vertical axis passing through the center of the canal (Fig. 3). Note that this consideration does not exactly mimic a head rotation as each inner ear is located approximately 3 cm away from head vertical axis. However, as stated by many researchers [], we consider that the SCCs are only stimulated by head angular accelerations and thus we do not take into account tangential as well as normal acceleration components that appear during any rotation of the head. The imposed rotation is maintained during 15 s and the computation is performed till 30 t s = in order to investigate fluid dynamics and cupula motion after the deceleration phase. We consider a Young's modulus for the cupula of 5 Pa and a Poisson's ratio of 0.48 [START_REF] Selva | Mechanical properties and motion of the cupula of the human semicircular canal[END_REF].

3.1.Results at the onset of the imposed rotational motion

The displacement of the cupula is shown in figure 4. It can be observed that displacement begins at the base of the cupula close to the crista, which involves that initial movements of the cupula produce a shear type deformation right above the sensory epithelium. Thus, sensory hair cells are presumably stimulated as soon as head motion starts which highlights the hypersensitivity behavior of the sensor. After a certain limit is reached, about 0.1 s, maximal displacement spreads toward the center of the cupula. This behavior is consistent with previous studies of [START_REF] Mclaren | The configuration of movement of the semicircular canal cupula[END_REF] who measured the positions of oil droplets, which were injected in the cupula of the bullfrog, following the compression of the canal wall.

Velocity profiles of the fluid flow in the slender part of the canal are plotted in figure 5. In about 0.04 s -0.05 s the fluid flow is analogous to a Poiseuille flow. Indeed, the velocity profile tends to a parabola, with the fluid in the center of the canal having the greatest speed. This result is consistent with previous analytical studies. For instance, [START_REF] Groen | The function of the vestibular organ[END_REF] assumed a fully developed Poiseuille flow in a straight tube to investigate the dynamics of semicircular canal flow and cupula motion, while [START_REF] Van Buskirk | The fluid mechanics of the semicircular canals[END_REF] shown that endolymph volume displacement resulting from a step change in angular velocity under the non-steady state flow assumption can be approximate by the Poiseuille steady-state flow relation.

3.2.Dynamic behavior of the structure during the imposed rotational motion

We now study the response of the sensor during both the acceleration and deceleration phase of the imposed rotational motion. Figure 6 presents the fluid flow in the canal and cupula motion at the onset of the simulation. At the beginning of the rotation of the canal, the fluid lags behind due to its inertia. Therefore, fluid flow relative to the wall of the canal is oriented in the opposite direction of the imposed rotational motion during about 0.25 s (Fig. 6). Meanwhile, this flow, represented by the arrows in figure 6, exerts a pressure across the cupula, and thus deflects it in the opposite direction of rotation as well. The deflection of the cupula reaches a maximum value close to 15 µm for the set of elastic properties retained. At time t=0.3 s, even though the canal still experiences a rotational motion, the cupula starts to return to its rest position due to its elastic properties. In addition, because of the small diameter of the duct and the viscosity of the fluid, the latter tends to catch up with the rotation of the canal, eliminating little by little the relative movement between the fluid and the canal. One can notice that the maximum fluid velocity is decreased by a factor 20 between time instants 0.3 s and 0.6 s.

Figure 7 shows the dynamic behavior of the sensor when its rotation is suddenly stopped. The canal experiences a constant angular velocity until t=15 s, and then the movement of rotation is stopped in 0.3 s. At time t=15.1 s, the fluid is still in motion within the canal due to its inertia. As a consequence, the cupula, which was returned to its rest position, is deflected in the opposite direction than previously (Fig. 7). Once again, we can note that cupula deflection starts near the sensory epithelium and then spread toward its center. Finally, the cupula returns to its initial position which provokes a slight counter clockwise fluid flow.

Figure 8 shows the time-dependent displacement of a point located at the center of the cupula. It can be seen that the cupula experiences two deflections in opposite directions that are due to the acceleration and deceleration phase of the motion. This behavior is consistent with previous analytical models such as the classic torsion pendulum model [START_REF] Steinhausen | Uber die Beobachtung der Cupula in den bogengansampullen des Lebenden Hechts[END_REF][START_REF] Groen | The function of the vestibular organ[END_REF][START_REF] Mayne | A systems concept of the vestibular organs[END_REF]. Indeed, cupula displacement is traditionally described by two times constant: a short time constant which governs the fast deflection of the cupula at the onset of the rotational motion, and a long time constant which governs the slow return of the cupula to its rest position.

Preliminary results of the entire 3-D model

The final goal of the work presented in this paper is to obtain a complete fluid-structural finite-element model of the entire set of semicircular canals + cupulae + utricle. Thus we extend the previous 2-D model of a single canal to three-dimensional space by taking into account all the angular sensors.

The geometry of this 3-D model has been developed with the CAD software CATIA V5. Then, using the STEP file exchange format, the model has been imported into Comsol Multiphysics in which meshing, simulations, and data processing have been performed (Fig. 9). Note that dimensions of the vertical canals have been taken identical as the horizontal canal, and that all the canals are assumed orthogonal.

Figure 9b shows the mesh of this model which is comprised of 45,408 quadratic tetrahedral elements that represent 498,112 degrees of freedom. The large number of elements is primarily due to the very small diameter of the canal cross-section (≈0.32 mm) which is one order of magnitude lower than the diameter of the circular path formed by each canal (≈3.2 mm). In addition the connections between each canals and the utricle has to be meshed finely to improve elements quality and to ensure convergence of the simulation. In the present case, simulations of the model are run on a computer having two dual-core processors, 8 Go RAM, and a 64-bit linux operating system.

Results of the simulation in response to head angular velocity are presented in figure 10. In particular, field displacement of the cupulae and fluid velocity at time t=0.1 s are plotted. At this instant, the fluid flow is in the opposite direction of the imposed clockwise rotational motion. The cupula of the lateral is thus deflected as predicted by the previous 2-D model. Note that cupula of the anterior canal also experiences a slight deflection due to the gravity field. Cupula of the posterior canal remains at its rest position as there is no fluid flow within this canal and as this cupula is almost aligned with the gravity vector.

These results are very promising as they are in good correlation with the previous 2-D model. This first 3-D model will be used to run further simulations in order to study various cases of vestibular disease -such as the Ménière's disease which is due to an overpressure in the labyrinth -and to perform parametric analysis so as to study sensitivity effect, for instance of the mutual orientation of the canal plane.

Conclusion

The present paper deals with fluid-structural finite-elements model of the angular sensors of the human inner ear. These models were developed using the finite-element analysis software Comsol Multiphysics 3.5a. The endolymph has been modeled as an incompressible Newtonian fluid while the cupula was considered as an elastic solid. The fluid-structure interaction problem has been treated with a strong coupling between the fluid flow and structural displacements.

First, we have considered a two-dimensional model of a single canal and have investigated fluid flow as well as cupula motion during head rotation. Results of this model were in good agreement with previous analytical studies. We have shown that in about 0.04 s fluid flow in the slender part of the semicircular canal can be approximated by a Poiseuille flow. We have also demonstrated that for a step change in angular velocity cupula displacement begins near the sensory epithelium and then spreads towards its center. Finally, displacement of the cupula through time was in good agreement with the classic torsion-pendulum model, i.e. described by two exponential terms defined by a short and long time constant respectively.

Second, a three-dimensional finite-element model of the entire set of semicircular canals and cupulae has been proposed. The associated geometry has been constructed using a CAD software (CATIA V5) which has then been imported into the finite-element Comsol Multiphysics software. For a constant head rotation, promising results -in good correlation with the 2-D model previously developed -have been obtained as regard fluid flow and cupulae displacement. However further simulations for different profile of head motion need to be performed to validate this 3-D model. The main prospect of this work is to use the final 3-D model to simulate different kind of vestibular disorders such as the Ménière's disease, which is caused by an increase in pressure and volume of the endolymph in the whole vestibular system.

List of figures

Figure 1 .

 1 Figure 1. (a) Visualization of the inner ear and close-up sketch of the semicircular canals. (b) 2D model of the cross-section of the horizontal semicircular canal. Note that the cupula is modeled by a 400 µm thick section.

Figure 2 .

 2 Figure 2. Visualization of the model subdomains that have different conditions for mesh displacement. A prescribed displacement is applied to the mesh of the fluid domain. The computed displacement sU is applied to the solid domain while the mesh of the fluid subdomains on both sides of the cupula are free to move.

Figure 3 .Figure 4 .

 34 Figure 3. Rotational motion applied to the semicircular canal. (a) Moving mesh. (b) Profile of angular velocity applied to the structure.

Figure 5 .Figure 6 .

 56 Figure 5. Evolution of the velocity of the fluid in the slender part of the duct at the beginning of the rotational motion. The fluid flow becomes analogous to a Poiseuille flow in about 0.05 s.

Figure 7 .

 7 Figure 7. Fluid velocity (left) and cupula displacement (right) when the rotational motion is stopped.Visualization in the ALE reference frame (moving mesh). Due to its inertia the fluid is still in motion within the canal when the latter is stopped. Thus, the cupula is deflected in the opposite direction first deflection experienced at the onset of the rotational motion.

Figure 8 .Figure 9 .

 89 Figure 8. Displacement of the center of the cupula during a constant angular rotation which ends at time t=15 s.Because of the angular acceleration the cupula is deflected according to a fast time constant (≈6 ms) and then returns to its rest position according to a long time constant (≈6 s). The inverse effect is observed when the rotational motion is stopped, i.e. when the canal experiences a deceleration.Time (s)

Figure 10 .

 10 Figure 10. Results provided by the simulation of the final 3D model of the semicircular canals. (a) Field of fluid velocity at time t=0.1s. (b) Field of cupulae displacement at time t=0.1 s. The cupula of the lateral canal is deflected because of the imposed motion of rotation of the canals. The cupula of the anterior canal is slightly deflected at the beginning of the simulation because of the gravitational field. Note that cupula of the posterior canal does not experience any deflection as there is no fluid flow within this canal and as this cupula is almost aligned with the gravity vector.

Acknowledgements

The authors thank Nicolas Huc of COMSOL for tremendous support while performing the simulation study.


Keywords: Damping, power ratios, robustness, SHM

In composites structures, cracking, delamination will cause changes in the measured dynamic response of structure and so on experimentally modal parameters. Estimation of damping in structural control often poses a difficult problem especially using broadband experiments. If these estimations are faulty, it is difficult to propose a robust SHM algorithm. Recently H.P. Yin introduced the optimal power ratios damping estimator. A new theoretical basis of the bandwidth method for the damping estimation from frequency response functions (in case of a single degree of freedom system) has been proposed. The main goal of this paper is to study the robustness of this enhanced damping estimator (AIPR) on simulated signal (sampling frequency, Signal to Noise Ratio and damping level/density), and also compare its performance with industrial improved estimator like "Polymax" on experimental FRFs. The pole shifts would be studied as a change in the frequencydamping plane function of level and density of damage.

Introduction

Since the dynamical behavior of structures is influenced by damage, one is able to detect the occurrence and evolution of damage by monitoring the modal parameters. Besides the detection of structural damage, modal analysis also enables the localization and quantification of the damage since also spatial information is present in the model (mode shapes). Compared to nondestructive testing (ultrasound techniques, analysis of magnetic fields, radiology, thermal methods) vibration-based analysis allows a mixed global/local analysis of the structure with the potential of being applied in situ (laser vibrometer, optical sensor, etc) [START_REF] Doebling | Damage Identification and Health Monitoring of Structural and Mechanical Systems From Changes in Their Vibration Characteristics: A Literature Review[END_REF][START_REF] Sohn | A Review of Structural Health Monitoring Literature form 1996-2001[END_REF]. The investigation focuses in the estimation of damping ratio because among the modal parameters, it is the most difficult to model/estimate. We first analyze the reliability of damping estimation using 2 methods AIPR (Average Inverse Power Ratio) and RFP (Rational Fraction Polynomial) on analytical (supervised) data then we will compare AIPR Paper 57

On the Robustness of the Average Power Ratios in Damping Estimation: Application in the Structural Health Monitoring of Composites Beams results with classical POLYMAX results on benchmark studies of the effect of impact damages on modal parameters on composites and especially thin laminates beams [START_REF] Niemann | Damage localization using experimental modal parameters and topology optimization[END_REF][START_REF] Shahdin | Correlating low energy impact damage with changes in modal parameters: A preliminary study on composite beams[END_REF][START_REF] Shahdin | Damage monitoring in sandwich beams by modal parameter shifts: a comparative study of burst random and sine dwell vibration testing[END_REF][START_REF] Montalvão | A review of vibration-based structural health monitoring with special emphasis on composite materials[END_REF].In a previous work we study the influence of several parameters on damping estimation [START_REF] Morlier | Original statistical approach for the reliability in modal parameters estimation[END_REF]: sampling frequency and numbers of spectral lines (which induce the notion of frequency resolution), Signal to Noise Ratio and damping level and density. We compare the robustness and the accuracy of the damping ratio estimation for 3 different methods: RFP, Continuous Wavelet Transform (CWT) and Least Square Complex Exponential (LSCE), RFP robustness even on noisy data have been proved. In experimental data the influence of close modes and/or damping ratio level on the estimation have been highlighted for classical methods that why we wanted to compare our POLYMAX results with new SDOF method like AIPR.

Theoretical background

Tracking of Poles for system identification on FRF

Excellent survey of Experimental Modal Analysis can be found in [START_REF] Ewins | Modal Testing: Theory, Practice and Applications[END_REF][START_REF] Avitabile | Experimental Modal Analysis -A Simple Non-Mathematical Overview[END_REF] concerning method of excitation (shaker, hammer etc), types of excitations (burst random, chirp, stepped sine…) and choice of the physical quantity to be measured (displacement, velocity, acceleration). Identification of the frequency-response function (FRF) by means of spectral analysis is usually considered well-investigated. However, currently the choice of the identification parameters (record length, correlation window size, number of averages, etc.) is very subjective and weakly formulated. The choices made are largely empirical and require the experience of a skilled control engineer [START_REF] Arkova | Optimal spectral resolution in system identification[END_REF]. The equations of motion for a vibrating structure are commonly derived by applying Newton's second law :

) t ( f ) t ( x K ) t ( x C ) t ( x M = + + & & & (1)
The excitation forces f(t) and responses are functions of time (t), and M, C, K are the mass, damping and stiffness constants respectively. The equivalent frequency domain form of the dynamic model can be represented with Fourier transforms, as :

) K ) j ( C ) j ( M ( 1 ) j ( F ) j ( X ) j ( H 2 + + = = ω ω ω ω ω (2)
Modal parameter estimation is a special case of system identification where the a priori model of the system is known to be in the form of modal parameters. The identification process consists of estimating the modal parameters from Frequency Response Function (FRF) measurements. The current approach in modal identification involves using numerical techniques to separate the contributions of individual modes of vibration in measurements such as frequency response functions.

Each term of the FRF matrix can be represented in terms of pole location and a mode shape:

) ( ) ( ) ( * * k k k k j r j r j H λ ω λ ω ω - + - = (3)
Where, k r is the (n by n) residue matrix for the k th mode (* designates complex conjugates)

k λ (Pole value for mode k) = 2 1 k k k k j ζ ω ω ζ - + - (4)
with k ω the undamped natural frequency and k ζ the damping ratio of mode k.

It has been observed that presence of damage does not affect all the resonant modes equally. A basis for identifying those modes that are affected is therefore necessary for further analysis to pinpoint the damage location. This paper suggests the use of pole shifts as a means for identifying modes affected by damage [START_REF] Mitchell | Modal Test Methods-Quality, Quantity and Unobtainable[END_REF][START_REF] Richardson | Correlating minute structural faults with changes in modal parameters[END_REF] as impacts create friction zone between plies in the composites.

Rational Fraction Polynomial

The rational fraction polynomial is a SISO frequency domain method [START_REF] Richardson | Parameter estimation from frequency response measurements using rational fraction polynomials[END_REF] and is among the most used nowadays. It is based on the viscously damped model and the FRF is expressed as the ratio of two polynomials:

∑ ∑ = - = = N j j j N j j j j b j a j H 2 0 1 2 0 ) ( ) ( ) ( ω ω ω (5)
A linearized error function is defined between the measured FRF values and the model. Minimizing such an error function, a linear system of equations is obtained, from which the coefficients a k and b k are evaluated (Invfreqs function in matlab which identify discrete-time filter parameters from frequency response data). The resulting system is usually ill-conditioned. This situation is circumvented by reformulating the problem in terms of orthogonal polynomials. Knowing the resulting coefficients, the modal parameters are retrieved. A global version of this method also exists, i.e. its SIMO version, named global rational fraction polynomial (GRFP), where various FRFs are taken into account simultaneously [START_REF] Richardson | Global curve-fitting of frequency response measurements using the rational fraction polynomial method[END_REF]. A MDOF system is decomposed in the summation of the individual single degree of freedom contributions "SDOF" by residue decomposition (residue function in matlab). Then global parameters are extracted and averaged from every FRF using approximation of equation [START_REF] Shahdin | Correlating low energy impact damage with changes in modal parameters: A preliminary study on composite beams[END_REF]. In this work both methods (MDOF on analytical data, and SDOF on noisy data and experimental data) have been used.

Power Bandwidth Method and Inverse Power Ratio Method

The classical power bandwidth method [START_REF] Yin | A new theoretical basis for the bandwidth method and optimal power ratios for the damping estimation[END_REF] for the damping estimation from the amplitude of a FRF is based on the following formula

1 2 max - - = ρ ω ω ω ζ a b (6)
Where max ω is the maximum amplitude frequency, a ω and b ω are two frequencies at which the FRF has the same amplitude, ρ is the inverse power ratio defined as

2 2 max 2 2 max ) ( ) ( ) ( ) ( ) ( ) ( b a b a H H H H ω ω ω ω ω ρ ω ρ ρ = = = = (7)
Recently, an average inverse power ratio (AIPR) method has been proposed in [START_REF] Yin | An average inverse power ratio method for the damping estimation from a frequency response function[END_REF] for the damping estimation by using the following formula

1 2 max - - = - + ρ ω ω ω ζ (8)
Where ρ is an average inverse power ratio defined at two frequencies situated at the same distance δ from max ω

        + + - = 2 max 2 max 2 max 2 max ) ( ) ( ) ( ) ( 2 1 δ ω ω δ ω ω ρ H H H H (9) a b
In practice, the AIPR method is easier to use than the power bandwidth method which requires data interpolations for the determination of ω andω .

3 Supervised numerical simulation

Supervised damping estimation

This study is carried out on simulated FRF with the aim of identifying the factors which have the most significant effect on damping estimation algorithms [START_REF] Morlier | Original statistical approach for the reliability in modal parameters estimation[END_REF]. Here we compare AIRP and RFP. The analysis was applied to simulated analytical data set with known properties (Table 1). We can see that the 2 first modes are closed (coupled) whereas the third one is well separated. The equation ( 3) is used to generate analytical FRF using the variables of Table 1 as a sum of decoupled SDOF system.

Our numerical supervised experiments principally focus on the evaluation of AIPR sensitivity against frequency resolution, SNR (White Gaussian Noise). For comparison we added a well established algorithm RFP.

Damping estimation Vs frequency resolution

In order to show the Frequency Resolution dependence of the AIPR, 2 examples are shown on figure 1. It reveals as previously implemented (Yin) that the results is better with high frequency resolution and becomes interesting for low damping level (Figure 1). We can see on figure 1 that as damping diminish (peak sharpens) the AIPR method gives better estimation (at better frequency resolution in yellow dotted lines) even than RFP.

Damping estimation Vs SNR

In these experiments we use a constant Frequency Resolution of 0.25Hz, and for AIPR interpolation a resolution of 0.1 Hz. The basic principle of our method is to isolate each peak resonance around a bandwidth of 20 Hz for AIPR or RFP. Power ratio Damping Ratio (%) Figure 3: Damping estimation (2%) under different realization of a Gaussian noise process (SNR=30dB). RFP is able to fit the analytical model very well and so on limits the error in damping estimation whereas AIPR is underestimating the damping ratio.

According to 2% damping level tests (Figure 2 and 3), AIPR seems to be more sensitive to noise than RFP SDOF but gives good results (less than 20% errors) for a power ratio zone of 0.85-0.9. At very low SNR (20dB), several trial tests demonstrate (figure 4) that theoretical value is just between an overestimation (RFP) and an underestimation (AIPR). Same kind of results is observed for others modes. For very low damping level (0.1%), estimators indicate different behaviour. First they become to overestimate the damping ratio (Figure 5). RFP seems to to be closer (0.11%) whereas AIPR gives for low power ratio overestimate results around 0.12%. IRF is generated with a SNR of 60 dB (very low noise) (a). Damping estimation for low values 0.1% around the first mode (318 Hz) for noisy data (SNR of 20 dB) (b)

Finally for medium damping level (0.4%) that we are able to observe experimentally on composites materials, our numerical experiments demonstrate (Figure 6) that for very low SNR (noise corresponding to experimental results) AIPR is underestimating the theoretical results.

Experimental test rig

An interesting review of SHM in composites materials could be found in [START_REF] Zou | Vibration-based model-dependent damage (delamination) identification and health monitoring for composite structures. A review[END_REF].

Vibrations tests for measuring the modal parameters of the composite test specimens and impact tests in order to create damage in the specimens. The experimental set-up is based on Oberst method [START_REF] Wojtowicki | New approach for the measurements of damping properties of materials using Oberst beam[END_REF] and states that a free-free beam excited at its centre has the same dynamical behaviour as that of a half length cantilever beam. The test specimen is placed at its centre on a B&K force sensor 8200 which is then assembled on a shaker supplied by Prodera having a maximum force of 100N (Figure 7a). High quality frequency response functions are measured with the help of a Laser Vibrometer OFV-505 provided by Polytec. The centre of the test specimens is excited at Point 17 as shown in figure 7b anda The impact test system used to damage the composite beams is drop weight system.

The impactor tip has a hemispherical head with a diameter of 12.7 mm. The size of the impact window is 80 x 40 mm 2 which allows all the impact points to have the same boundary conditions and all the four ends are clamped.

The five composite specimens [START_REF] Shahdin | Significance of low energy impact damage on modal parameters of composite beams by design of experiments[END_REF] are impacted around the barely visible impact damage limit (BVID). BVID corresponds to the formation of an indentation on the surface of the structure, which can be detected by detailed visual inspection and can lead to high damage. In the aeronautical domain, BVID corresponds to an indentation of 0.3 mm after relaxation, aging etc (according to Airbus certifications). Impact energy of 10 J giving an initial indentation depth of 0.55 mm, shall be considered as the BVID limit. We focus on the results of beam 4 on the Table 2. In fact damage zones are not consistent with the energy of impact even with controlled system we used. Thus it appears local modes in the FRF different for each tested beam, but it is possible to update a finite element with equivalent damage with high correlation using only the first 4 bending modes according to Table 3.

Pole Shift results

Frequency and damping changes are studied with the help of the first four bending modes as they have the largest amplitudes for the type of test configuration presented in this article. As discussed previously, delamination induced damage in composites leads to an increase in damping and a decrease in natural frequency. This effect is more significant in the high frequency range. This fact is verified by our experimental results which show that the difference in natural frequencies between the damaged states (4 and 8) and the undamaged state for the first mode is very small. But this difference in frequencies increases for the higher modes. For the 2 nd and 3 rd bending modes, the variation of natural frequency as a function of the undamaged (0) and the two damage states (4 and 8) is presented in figure 8. show that the decrease in natural frequencies with the increase in damage is more significant in case of the higher impact energies. The damping ratios estimated by POLYMAX often tend to increase with impact number.

0 0,1 0,2 0,3 0,4 0,5 0,6 0
It can be seen from the results in figure 8 that the damping ratio increases with increase in damage in the beam 4. Others tested beams highlight the same effect. Furthermore, the change in damping ratios between the three states (0,4,8) for the 2 nd and 3 rd modes is very small whereas frequency shift is clear. However unlike natural frequencies, the increase in damping ratio between the damaged and the undamaged states is not always coherent with the impact energy level, due to the complex nature of damping and the difficulties in its estimation. But in case of Mode 4 damping ratio exhibits quasi linear dependence on the energy of impact.

Damping Estimation AIPR versus POLYMAX

In order to understand the slight shift in damping ratio for mode 2 and 3 estimate by POLYMAX [START_REF] Peeters | The PolyMAX frequency-domain method: a new standard for modal parameter estimation[END_REF] we analyse the experimental results with AIPR algorithm (isolate each peak resonance around a bandwidth of 20 Hz ) taking 10 random FRFs for mode 1 and 3 random FRFs for mode 2 (very noisy as you can see on Figure 9). Frequency (Hz) FRF (dB)

Figure 9 Typical FRF with local modes (small peaks) between the 4 first bending modes (higher amplitudes)

We can illustrate the results using the figure 10: as for medium damping level supervised data, the RFP method overestimates the damping (reference POLYMAX 0.28%) and AIPR method underestimate the reference damping estimation. We can also see that the variance in damping estimation is lower for low power ratio (<0.4). Finally all the results are resumed in figure 11 and globally AIPR (averaged over all the power ratio) results correlate well POLYMAX results and thus these 2 modes are not sensitive to damage comparing to mode 4 where damping ratio exhibits quasi linear dependence on the energy of impact and frequency shift is more important.

Figure 11: Damping estimation (Damping Ratio in abscissa in %) using AIPR and POLYMAX, for different number of impacts.

For engineer one advantage of AIPR is the simplicity: every one knows the classical bandwidth method (Power ratio of 0.5). But here instead of determining the bandwidth for a given value of power ratio which requires numerical interpolations, the method involves an average inverse power ratio calculated at two frequencies symmetrically located from a peak. In this way, the damping estimation from a FRF becomes straightforward for simple modal test cases and would be sufficiently accurate provided that the coupling effects are not strong and the frequency resolution is high enough to determine accurately the peak amplitude frequencies.

Conclusion

This paper focuses on the estimation of damping ratio often difficult to estimate accurately. We first analyze the robustness of AIPR method versus classical RFP method using analytical (supervised) data. Our results demonstrate different behaviour function of the damping ratio level. For low damping, low power ratio should be preferred (around 0.3), notably dealing with high noise. The experimental modal tests aim at showing the effect on modal parameters (pole shift) of multisite localized impacts on composites laminates. We focus especially in 2 modes where the pole shift was not clear using POLYMAX. The simplicity of AIPR method becomes here an advantage for rapid analysis. But the limitation is that a classification of "interesting FRFs" (on the 33 experiments) should be made for a good AIPR estimation (sufficient peak amplitude). We also choose to study the averaged value of AIPR, but it is interesting to see that for low power ratio the variance in the damping estimation is lower. Finally AIPR method well correlates POLYMAX method for noisy experimental data.

Figure 1 :

 1 Figure 1: Identified damping ratio on theoretical damping ratio 2% (a) and 0.5% (b) for several value of power ratio for first mode isolated at 318.6 Hz (a) and for third mode (b) at 508.5 Hz. . Thin red line is FR of 1.5 Hz, orange thick line is 0.625 Hz, yellow dotted line is 1 Hz. RFP is represented by blue line very close to theoretical value. AIPR is above 5% of error but always cross the real value at high power ratio.

Figure 2 :

 2 Figure 2: One realization of damping estimation for normal values 2% around the first mode (318 Hz). RFP (thin blue line) gives very good results whereas AIPR (black thick line) gives for low power ratio results around 2.05%. FRF is generated with a SNR of 60 dB (very low noise).

Figure 4 :

 4 Figure 4: Damping estimation (2%) under different realization of a Gaussian noise process (SNR=20dB). Theoretical value is just in the middle of the 2 estimators

  0

Figure 5 :

 5 Figure 5: Damping estimation for low values 0.1% around the first mode (318 Hz).RFP seems to to be closer (0.11%) whereas AIPR gives for low power ratio overestimate results around 0.12%. IRF is generated with a SNR of 60 dB (very low noise) (a). Damping estimation for low values 0.1% around the first mode (318 Hz) for noisy data (SNR of 20 dB) (b)

Figure 6 :

 6 Figure 6: Damping estimation (0.4%) under different realization of a Gaussian noise process (SNR=20dB). Theoretical value is just in the middle of the 2 estimators but the noise has less effect as damping becomes smaller (0.4%)

  high frequency resolution of (∆f = 0.25Hz) for precise modal parameter estimation. Burst random excitation is used which is a broadband type excitation signal [0-200 Hz]. The signal is averaged 10 times for each measurement point. Hanning windows are used for both the output and the input signals. Response is measured at 33 points that are symmetrically spaced in three rows along the length of the beam (figure7b). The modal parameters are extracted by a frequency domain estimation method (POLYMAX) based on an automatic extraction using stability diagram.

Figure 7 :

 7 Figure 7: Experimental set-up for vibration testing (a) and Composite beam specimen with location of damage, excitation and measurement points (b)

Figure 8 (

 8 Figure8(a) and (b) show that the decrease in natural frequencies with the increase in damage is more significant in case of the higher impact energies. The damping ratios estimated by POLYMAX often tend to increase with impact number.

Figure 10 :

 10 Figure10: Damping estimation (undamaged case) using AIPR and RFP for 3 random points on 33 measurements points around 3 rd bending mode (800 Hz). The red line is the average value for RFP (0.358%) and the green line is the average value for AIPR (0.183%) very close to LMS results.

Table 1 :

 1 Pertinent parameters used with analytical data

	Natural Frequency (Hz)	318.6 341.1 508.5
	Frequency Resolution (Hz)	0.25	0.625 1.5
	Damping Level (%)	0.1	0.5	2
	Signal to Noise Ratio (dB)	20	30	60

Table 2

 2

		Impact test parameters
	Beam No Energy of Impact (J)	Height	Velocity of impact
			(mm)	(m/s)
				measured
	3	10 (BVID)	552.9	3.24
	4	12	663.5	3.52
	5	14	774.1	3.84

Table 3

 3 Frequency estimation for the first 4 modes versus number of impacts

	Bending	Undamaged	4	8
	Modes		impacts	impacts
	1	41.5	39.8	39.1
	2	319.0	308.2	305.1
	3	802.1	773.1	760.2
	4	1543.7	1477.7	1465.9

University of Toulouse, INSA, UPS, Mines Albi, ISAE, Institute Clément Ader, France

University Paris-Est, UR Navier, Ecole des Ponts ParisTech, Marne-La-Vallée, France