J Limido

C Espinosa

M Salaun

C Mabru

R Chieragatti
email: remy.chieragatti@isae.fr

J L Lacome

Metal cutting modelling SPH approach

Keywords: high speed machining, smoothed particle hydrodynamics, SPH, chip formation, tool wear, cutting forces, milling

The purpose of this work is to evaluate the use of the smoothed particle hydrodynamics (SPH) method within the framework of high speed cutting modelling. First, a 2D SPH based model is carried out using the LS-DYNA® software. The developed SPH model proves its ability to account for continuous and shear localised chip formation and also correctly estimates the cutting forces, as illustrated in some orthogonal cutting examples. Then, the SPH model is used in order to improve the general understanding of machining with worn tools. At last, a hybrid milling model allowing the calculation of the 3D cutting forces is presented. The interest of the suggested approach is to be freed from classically needed machining tests: Those are replaced by 2D numerical tests using the SPH model. The developed approach proved its ability to model the 3D cutting forces in ball end milling.

Biographical notes: Jérôme Limido has experience from research and advanced engineering within aerospace applications. His main work has focused on processes involving large deformations, both experimentally and numerically. He has special interests in advanced numerical methods (smoothed particle hydrodynamics) and fatigue of materials. He is R&D responsible at IMPETUS Afea France and teaches Advanced Computational Mechanics and Numerical Methods at ISAE. C. Espinosa, after an A-level diploma in Mathematics-Physics (1982), two licence degrees at the University of Bordeaux I in Theoretical Mathematics and Applied Mathematics (1986), he received his Master degrees at this University (M1 of today Master 1987Master , M2 1988) ) in Applied Mathematics. He obtained his PhD in Theoretical and Computational Mechanics (1991) for the development of a non-linear damage anisotropic behaviour for aeronautical composite materials subjected to low velocity impacts. He was hired by DYNALIS in 1992, in charge of LS-DYNA distribution in France. He is an ISAE Associate Professor, for teaching and research in Advanced Computational Mechanics and Numerical Methods.

Michel Salaun teaches mathematics and scientific computing at ISAE. His major field of research is computational mechanics. He developed a mixed finite element for thin elastic shells. This formulation was used for the development of an Eulerian approach for large displacements of thin shells. He works also on numerical modelling of incompressible flows (stream functionvorticity and vorticity-velocity-pressure formulations). During the last years, he was interested in the application of extended finite element methods (XFEM) for the modelling of cracks in thin structures, first, in the frame of bidimensionnal linear elasticity and, second, for the Kirchhoff-Love plate model.

Catherine Mabru is an Associate Professor at ISAE. She teaches Elasticity, Thermo-elasticity and Plasticity. Her research mainly deals with fatigue of metallic materials and focuses on the interaction between surface condition and mechanical behaviour. Coupling experimental and numerical approaches, she particularly worked on the effect of machining on fatigue strength and the influence of anodic treatment on thermal fatigue degradation mechanisms of anodic film, both concerning aluminium alloys. Recently, she was interested in studies that aim to isolate surface and subsurface behaviour by analysing thin samples strength under thermo-mechanical fatigue loading.

Remy Chieragatti is an Associate Professor at ISAE. He teaches Fatigue and Mechanics. His research mainly deals with fatigue of metallic materials and focuses on the interaction between surface condition and mechanical behaviour.

Introduction

Machining is the most used process in industrial components production, while involving very complex fully 3D physical phenomena. A considerable effort has been made by the scientific community to develop accurate theoretical and numerical machining models, each one aiming to represent the relevant phenomena of the physics at the considered scale. Among the objectives of these models reproducing this complexity and linking the quality of the surface to the machine control are often encountered. The focus of this work is to propose a numerical model able to reproduce the efficient cutting physics and quality of the generated surface and link local machining models to the machine driving parameters, in the final goal of predicting the effect of the surface quality on the component fatigue life. The application field is limited here to aeronautical representative samples. To be more precise, the purpose of this work is not to illustrate the use of a (n + 1)th numerical model, nor to represent the whole complexity of the physics in the tool-matter interaction. The originality is to propose a new methodology of machine control to 'optimise' in a way the quality of the machined surface regarding fatigue life duration, through the use of an appropriate numerical model, under constraints of being relevant and mechanically consistent with the physics of two mechanical scales. The invisible part of this work is the mathematics, computational mechanics (theory) and code developments (numeric) done prior to the numerical tests. The fatigue life prediction part of the work that uses the numerical tests and the machine driving 'optimisation' phase are not presented here.

Why the proposed methodology overcomes the key limitations of existing numerical works? Because numerical models are themselves limited compared to fine theoretical approaches, and due to additional computational and numerical limitations, numerical models are mostly limited to the orthogonal cutting framework, especially when it is aimed to compare numerical results and real complex and uncertain experimental measurements. Thus, 2D plain strain studies are commonly dealt. In these conditions, the full 3D and very complex interacting machine control cutting parameters can be resumed to two of them considered as independent, the cutting speed (Vc) and the feed (f). In that field, most of the numerical models are based on Lagrangian or arbitrary Lagrangian Eulerian finite element methods (ALE FEM). Lists of available publications are of course very large. Taking the risk of not being exhaustive, let us cite only two of them [START_REF] Bil | A comparison of orthogonal cutting data from experiments with three different finite element models[END_REF] and [START_REF] Marusich | Effects of friction and cutting speed on cutting force[END_REF]. These references papers in particular show two major difficulties closely related to FEM. First contact and friction models must be introduced thus leading to the necessity of accounting for the thermo-mechanical exchange complexity of machining through the interface description, furthermore a line in 2D. In most cases, the Coulomb model is used. It is very simple to implement and use but the contact complexity is, to us, poorly resituated. Moreover, the friction parameter is shown to be a very sensible coefficient often used to tune the FEM cutting forces on experimental results. Thus, FEM can not be predictive. The workpiece/chip material separation model is the second aspect of the modelling difficulties. Lagrangian based FE methods present the disadvantage of leading to large grid distortions since the matter and the grid are strictly connected in the absence of additional 'opening' technique. Use of ALE in the FEM remedy since the matter-grid connection is somewhat relaxed and it allows for the boundary of the matter to evolve while remaining Lagrangian. Remeshing techniques are associated to remapping methods. This is hugely time consuming, and in order to be not energy dissipative and accurate enough, this implies very fine mesh sizes to follow the matter fluxes that open to question the consistency between the solid continuum mechanics and numerical hypotheses.

Regarding the real 3D industrial machining cases (milling for example), numerical models existing at the beginning of this work are not yet rather powerful because of the huge computing time and still remaining numerical problems. Analytical models are often preferred [START_REF] Altintas | Manufacturing Automation-Metal Cutting Mechanics Machine Tool Vibrations and CNC Design[END_REF].

The developed approach is based on the smoothed particle hydrodynamics (SPH) method in the frame of the LS-DYNA® hydrodynamic software (Hallquist, 2006). It is shown in Section 2 that SPH is able to model cutting, avoids the need of contact algorithm or remap/remesh additional techniques. In Section 3, SPH cutting model applications are outlined and compared with other numerical or experimental data. The use of SPH as a numerical tool for a better understanding of the chip formation with worn tools is also presented. Finally, in Section 4, our hybrid (analytic-SPH) 3D milling model is presented. This orthogonal cutting based methodology gives as a result a predicted machined surface geometry without the need of 3D numerical cutting model or real test.

SPH cutting model

Basic principles of the SPH method

SPH method is a meshless Lagrangian technique. Material properties and state variables are approximated by their values on a set of discrete points, or SPH particles. This avoids the severe problems of mesh tangling and distortion which usually occur in Lagrangian analyses involving large deformation. For more details on the method used, the reader can refer to [START_REF] Halquist | LS-DYNA Theory Manual[END_REF] and [START_REF] Liu | Smoothed particle hydrodynamics: a meshfree particle method[END_REF].

Model description

Several assumptions were made in order to reduce the model size and the computation time, allowing the development of a useful tool.

The model is implemented in the orthogonal cutting framework, thus in 2D (a 3D model can easily be developed but we will show in Section 4 that a coupled analytical-SPH approach could be more efficient, see Figure 2). The tool is supposed rigid and has an imposed velocity. During the process, the computation time is reduced by using an imposed tool velocity ten times higher than the real velocity. This assumption is usually used in simulation of stamping processes. It is valid as long as the accelerated mass is low and the material behaviour is slightly influenced by the strain rate.

Accurate and reliable flow stress models are highly necessary to represent work material behaviour under high speed cutting conditions. In our model, the constitutive material law proposed by [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF] is used and all associated parameters result from the literature.

SPH cutting model capabilities compared to classical FEM approach

An overview of the main specificities of the SPH cutting model approach compared to classical Lagrangian FE models is summarised in Table 1. SPH method applied to machining modelling involves several advantages. First, high strains are easily handled. The particles move relatively to each other in a disordered way during the deformation. This can be considered as a particles rearranging without topological restriction, thus no remeshing is needed.

Another advantage induced by SPH is the 'natural' chip/workpiece separation. The relative motion of the particles creates the opening. The new free surfaces are given by the particles positions. Thus, the workpiece matter 'flows naturally' around the tool tip (Figure 1). In the same way, the SPH method presents an original aspect regarding contact handling. Indeed, when a workpiece particle 'sees', in its neighbouring, tool particles, then workpiece particle circumvents the tool. Thus, friction is modelled as particles interactions and friction parameter does not have to be defined. SPH friction modelling must be studied in-depth but it offers a very interesting alternative to traditional definitions.

Some specificities of the SPH model implementation

The implemented SPH model within the framework of LS-DYNA reveals some specific aspects.

Renormalised SPH formulation

Due to the lack of neighbours, classical SPH approximation is not correctly calculated on the boundaries. Indeed, a particle located on a free edge only interacts with existing particles located on one side of the free edge. LS-DYNA proposes an alternative formulation to particle SPH approximation which corrects this deficiency. This formulation is called renormalisation. It is based on the works of [START_REF] Randles | Smoothed particle hydrodynamics: some recent developments and applications[END_REF] and [START_REF] Vila | SPH renormalized hybrid methods for conservation laws: applications to free surface flows[END_REF]. The renormalised formulation improves the SPH method consistency (order 1). Two advantages are obtained: a better precision when the particles are disordered and a better approximation of the calculated quantities at the boundary.

In the machining modelling case, standard SPH underestimates the chip curve (Figure 3). The use of the renormalised formulation allows a better modelling of the chip curve and a better evaluation of the matter state on the surface (Figure 3).

Numerical instabilities

The SPH method is prone to two types of numerical instabilities: tensile instability and zero energy modes. These aspects are discussed in details in [START_REF] Belytschko | Stability analysis of particle methods with corrected derivatives[END_REF].

We faced numerical instability of the SPH method for the cutting model. Figure 4 shows an example of numerical fractures. Increasing the physical velocity of the tool decreases these instabilities. Currently, the SPH method in LS-DYNA does not provide any advanced solution to solve this problem but a work is in progress.

SPH cutting model applications

Two applications are presented here: the first one reproduces a continuous chip process and the second one a shear localised chip. They are studied and compared to experimental results and numerical FEM (AdvantEdge; [START_REF] Marusich | Effects of friction and cutting speed on cutting force[END_REF] results. AdvantEdge is a reference in machining simulations. It is an explicit dynamic, thermo-mechanically coupled FEM package specialised for metal cutting. Here, comparisons are carried out on the chip morphology, stress distribution and the specific cutting forces. Other comparisons can be found in [START_REF] Limido | SPH method applied to high speed cutting modelling[END_REF].

Continuous chip: Al6061-T6

We first present the result of a cutting problem using material Aluminium alloy Al6061-T6. The process parameters are defined as follows: speed 10 m/s, rake angle 5 , feed 250 µm and edge radius 25 µm. The material model parameters result from [START_REF] Lesuer | Experimental investigations of material models for Ti-6AL-4V titanium and 2024-T3 aluminum[END_REF]. All the AdvantEdge and experimental results presented in this part are based on [START_REF] Marusich | Effects of friction and cutting speed on cutting force[END_REF].

Chip morphology and stress distribution

LS-DYNA and AdvantEdge chip morphology results are presented in Figure 5. The considered Aluminium alloy produces continuous chip in the speed and feed range studied [START_REF] Marusich | Effects of friction and cutting speed on cutting force[END_REF]. LS-DYNA and AdvantEdge models results are in agreement with these experimental observations. The AdvantEdge model overestimates the chip thickness and the LS-DYNA SPH model underestimates it.

Shear localised chip: Ti6Al4V and worn tool

In this section, the Ti6Al4V alloy machining is used under dry conditions. When machining titanium alloys with conventional tools, the tool wear rate progresses rapidly, and it is generally difficult to achieve a cutting speed of over 60 m/min. The tool wear analysis method developed here is based on a chip formation analysis using experimental and numerical results. The main idea is to use the developed SPH model to study the cutting forces and the chip formation with new and worn tools and compare these results to available experimental data. This work was completed in collaboration with the Lamefip ENSAM Bordeaux France [START_REF] Limido | Toward a better understanding of tool wear effect through a comparison between experiments and numerical model[END_REF].

The studied cutting conditions are summarised in Table 2. The tool geometry before and after machining are shown in Figure 7.

Chip morphology

Under the cutting speed and feed range studied, Ti6Al4V alloy produces shear localised chips (see Figure 8). Shear localised chips are characterised by oscillatory profiles. They result from adiabatic shear band formation in the primary shear zone of the workpiece material. The chip type obtained for the SPH model is in accordance with experiments (Figure 8). The SPH model shows a shear band thickness increase and a decrease of frequency of chip segmentation with wear, which is also observed experimentally. It is possible to identify (Figure 9) a metal dead zone in the chip formation mechanisms with worn tools. This area is reduced to a triangular zone in front of the tool tip which moves with the tool. This physical mechanism was already observed in experiments [START_REF] Kountanya | A high-magnification experimental study of orthogonal cutting with edge-honed tools[END_REF]. The SPH approach is able to predict this stagnation zone because of its meshless nature.

Cutting forces

Experimental results, illustrated in Figure 10, show a great influence of wear on the feed force (about 90%) and a less important increase of the cutting force (about 20%). The SPH cutting model results show a mean feed force increase of about 95% with tool wear and of the cutting force of about 10% (Figure 10). So, the SPH model is able to correctly predict these variations, induced by the tool wear by only taking into account the tool tip shape changes. We think that the metal dead zone formed in front of the tool tip is the main cause of the large feed force increase [START_REF] Limido | Toward a better understanding of tool wear effect through a comparison between experiments and numerical model[END_REF].

3D Milling forces via 2D SPH

Principle of the 3D cutting forces determination

The milling force modelling is not a simple extrapolation of the orthogonal cutting case because locally oblique cutting occurs. Recently, 3D finite elements milling models were developed but the computational times necessary are too important for exploitation within an industrial framework. The alternative suggested here is mainly based on a coupling of [START_REF] Altintas | Manufacturing Automation-Metal Cutting Mechanics Machine Tool Vibrations and CNC Design[END_REF] approach and our 2D SPH model. Altintas approach principle is to consider the total cutting pressures as the sum of the local cutting pressures along the edge(s). In other words, the cutter is represented as a sum of elementary tools.

In the developed model, the local cutting forces are calculated starting from the local cutting conditions (thickness and depth of chip) and from the cutting coefficients by equation (1) [START_REF] Altintas | Manufacturing Automation-Metal Cutting Mechanics Machine Tool Vibrations and CNC Design[END_REF]. Let us note that tangential components (t index), axial (a index) and radial (r index) of the cutting coefficients are separate in components related to shearing (c index) and in components related to the edge action (e index). The total cutting forces are then determined by integration of the equation (2). In order to obtain the total cutting forces by this method, it is necessary to determine on the one hand the local cutting conditions and on the other hand the coefficients of cut.

k k dF k k dF dF - - - ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = - - ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ - ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ψ ψ ψ ψ ψ ψ (

Local cutting condition determination

We developed a model which represents the intersection tool/workpiece like the intersection between a B-Rep surface (discretised tool) and a Z-Map volume (discretised workpiece). We chose to use this model in order to determine the local cutting conditions. Indeed, the B-Rep representation of the tool already proposes a discretisation of the cutting edges. Our local 'tools' are thus the triangles pairs of discretisation of the total tool. The local cutting conditions are given from the 'erased' chip volume by the travel of the local or elementary tool for the considered time step.

) ( ) ( ) chipi ti tci tei i chipi ri rci rei i chipi ai aci aei i V F K K s R z t V F K K s R z t V F K K s R z t = + Δ = + Δ = + Δ ω ω ω (3) F ti , F ri , F ai
tangential, radial, and axial elementary (i) tool force The suggested method for determination of the local cutting conditions is very simple (Z-Map+B-Rep). It is thus possible to evaluate the local cutting conditions for various cutting configurations. Nevertheless, analytical solutions exist for 'simple' milling conditions like rectilinear slotting [START_REF] Altintas | Manufacturing Automation-Metal Cutting Mechanics Machine Tool Vibrations and CNC Design[END_REF]. The developed model results are compared (Figure 12) with those of the reference analytical model developed by [START_REF] Altintas | Manufacturing Automation-Metal Cutting Mechanics Machine Tool Vibrations and CNC Design[END_REF]. The studied case is related to a rectilinear ball end slotting at constant cutting parameters. We note that the approximation induced by the developed model represents a weak error compared to the reference model (Figure 12). Let us note that the cutting coefficients used for this comparison are those given by [START_REF] Altintas | Manufacturing Automation-Metal Cutting Mechanics Machine Tool Vibrations and CNC Design[END_REF] experiments: K tc = 2172.1 N/mm2, K rc = 848.90 N/mm2, K ac = -725.07 N/mm2, K te = 17.29 N/mm, K re = 7.79 N/mm, K ae = -6.63 N/mm.

K tci , K tei , K rci , K rei , K aci ,
x i i i i i r i y i i i i i t i i i i a i z F k k F F k k F k k F F - - - ⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎢ ⎥ ⎢ ⎥⎢ ⎥ = - - ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ - ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ∑ θ ψ ψ ψ θ ψ ψ ψ θ (4)

Cutting coefficients identification

The cutting coefficients [equation ( 1)] are classically determined from many milling tests for which the cutting forces are measured. This approach is known as mechanistic [START_REF] Altintas | Manufacturing Automation-Metal Cutting Mechanics Machine Tool Vibrations and CNC Design[END_REF]. It permits to obtain very precise results but the determined cutting coefficients are only valid for the considered triplet 'range, tool, matter'. A simple change like the cutter diameter requires new experimental tests. The predictive character of this approach is thus very limited. [START_REF] Lee | Prediction of ball-end milling forces from orthogonal cutting data[END_REF] proposed a very interesting alternative to the mechanistic approaches. It permits to determine the cutting coefficients for many cutting conditions only on the basis of orthogonal cutting tests. This approach can be divided into three principal steps summarised here (more precise details can be found in [START_REF] Lee | Prediction of ball-end milling forces from orthogonal cutting data[END_REF][START_REF] Altintas | Manufacturing Automation-Metal Cutting Mechanics Machine Tool Vibrations and CNC Design[END_REF]:

• Orthogonal cutting force measurement, cutting and edge force parting, cutting coefficients identification K te , K re and K ae [equation ( 5)]. This cutting coefficients definition is obtained on the basis of following simplifying assumption: the angle of the primary shear band in orthogonal cutting is equal to the angle of the primary shear band of oblique cutting; the normal rake angle is equal to the rake angle of orthogonal cutting; the chip flow angle is equal to the angle of obliquity; average shear stress in the primary shear band and the friction coefficient are identical in the orthogonal case and the oblique case.

i K i K i i i i K i - + = + - + - = + - + - - = + - + τ β α β φ φ β α β τ β α φ φ β α β τ β α β φ φ β α β (

From 2D SPH model to 3D oblique cutting

The identification method of the cutting coefficients which we have just presented is classically based on orthogonal cutting experiments results. We propose to replace these physical tests by numerical tests on the basis of model SPH 2D presented in Section 2. The principal interest is to avoid the numerous orthogonal cutting tests which are necessary when a large variety of tools is studied. Indeed, cutting forces and chip thickness vary, in particular, according to the geometry of the orthogonal cutting tool (edge radius, rake angle). However, we showed into Section 2 that 2D SPH model permits to correctly take into account the tool geometry effect on the cutting forces. Thus, it seems interesting to replace the experimental phase of the cutting coefficients identification by a numerical SPH phase.

This suggested approach is evaluated in the following via the cutting coefficient identification in the case of aluminium 7075-T6 alloy.

SPH Model was used in order to determine the cutting forces and the chip thickness for 3 different depths of cut (or feed) (Table 3). The JC material model parameters used result from work of [START_REF] Fang | Method for modeling material constitutive behavior[END_REF] (Table 4). Obtained SPH Results are gathered in Figure13. The equations of the linear regression curves represented in grey permits to obtain K te = 15 N/mm and K re = 29N/mm. K ae is often negligible and fixed at a zero value [START_REF] Lee | Prediction of ball-end milling forces from orthogonal cutting data[END_REF]. Equations ( 6), ( 7) and (8) permit to obtain the Table 5 values. Equation (9) gives the three cutting coefficients K tc , K rc , K ac . Let us note that these coefficients vary according to the helix angle (i) as illustrated in Figure 14. We validate the proposed approach by comparison to a reference case resulting from [START_REF] Lamikiz | Cutting force estimation in sculptured surface milling[END_REF] work. This case corresponds to a 7075-T6 aluminium alloy ball-end milling of a slope with a 15° angle (Figure 15). The cutting condition studied is defined by the Table 6. The 3D model results are based on the cutting coefficients determined via the 2D SPH model (Figure 14).

The numerical and experimental results obtained by [START_REF] Lamikiz | Cutting force estimation in sculptured surface milling[END_REF] are compared with those of the 3D model developed, see Figure 16. We can note that our model permits to obtain comparable results to the Lamikiz model. Maximum obtained error compared to the tests is approximately 25% for the studied case. This one corresponds to the F y component and we can notice that this error is weaker for the F x and F z components. Let us note that the input data necessary to the implementation of the 3D milling forces calculation are limited to the geometrical definition of the cutter, the machining parameters and the JC material model parameters for considered alloy. The developed approach thus permits to limit in a very important way the quantity of needed tests. It is possible to use this model in order to evaluate the tool geometry influence on the cutting forces. An example concerning the tool radius influence in the case of a plane milling is illustrated in Figure 17. We note that the normal force on the generated surface (F z ) is very influenced by the ball-end tool radius. It would be possible to study the influence of each geometrical parameter for a tool family, as well as the influence of the machining conditions. Thus, determining optimum machining conditions could be possible, but that exceeds the framework of this study.

Conclusions

The results of the LS-DYNA SPH model were compared with experimental and numerical data. This study shows the relevance of the selected numerical tool. The SPH model is able to predict continuous and shear localised chips. The model also correctly estimates the cutting forces without introducing an adjusting friction parameter. Thus, comparable results compared to machining dedicated codes are obtained.

The SPH 2D cutting model has also been implemented as a helpful tool for understanding chip formation. The meshless nature of the method makes possible to represent a dominating physical phenomenon in the chip formation with strongly worn tool: the metal dead zone.

We also validated a 3D hybrid milling force model that couple our 2D SPH model and a semi-analytical model. This hybrid model permits to study industrial cases without building a 3D large and time consuming numerical model. It could be interesting to include our approach in a machining strategy optimisation study.

Figure 1 Figure 2

 12 Figure 1 Chip separation (see online version for colours)

Figure 3

 3 Figure 3 Standard/renormalised SPH formulation in metal cutting (see online version for colours)

Figure 4

 4 Figure 4 Numerical instabilities (see online version for colours)

  Figure 5(b) shows the Von Mises stresses during the continuous chip formation. Primary and secondary shear zone can be easily identified.

Figure 5

 5 Figure 5 Al6061-T6, (a) chip thickness AdvantEdge/LS-DYNA comparisons (b) SPH model Von Mises stress (see online version for colours)

Figure 6

 6 Figure 6 Al6061-T6 predicted cutting forces comparisons AdvantEdge/LS-DYNA

Figure 7

 7 Figure 7 Scanning electron microscopy images of (a) a new (b) a worn tool

Figure 8 Figure 9

 89 Figure 8 Experimental/numerical chip formation using (a) a new (b) a worn tool geometry (see online version for colours)

Figure 10

 10 Figure 10 Tool wear influence on cutting forces: SPH model compared to experimental data Cutting force increase with tool Wear

Figure 11

 11 Figure 11 Matter removal by an elementary tool of the cutter (see online version for colours)

(

FFigure 12

 12 Figure 12 Developed model compared to Altintas' (2000) results on a ball-end milling case (see online version for colours)

  force (orthogonal cutting) K tc , K rc cutting force coefficient (tangential and radial) K te , K re edge force coefficient (tangential and radial) h uncut chip thickness (feed in orthogonal cutting, also noticed f) b width of cut (also notices w in orthogonal cutting)The cutting forces are supposed to be linearly dependent on the depth of cut (or feed in orthogonal cutting). Edge forces are identified by extrapolation of the measured data for a null depth of cut (f = 0).• Identification of the cutting parameters: primary shear band angle [equation (6)],average shearing in the primary shear band [equation (7force (feed force in orthogonal cutting) h uncut chip thickness (feed in orthogonal cutting) Transformation from orthogonal cutting (2D) to oblique cutting (3D): cutting coefficient identification (9).

  9) K tc , K rc , K ac cutting coefficients (tangential, radial and axial) φ c primary shear band angle τ s mean shear stress in the primary shear band β n normal friction angle (tan β n = tan β a cos i) α r rake angle i local helix angle or angle of obliquity

Figure 13 TableFigure 14

 1314 Figure 13 Cutting forces predicted by SPH 2D model, α = 0 r = 20 μm AA7075-T6

Figure 15

 15 Figure 15Machining strategy of the studied case

Figure 17

 17 Figure 17 Tool radius effect on milling forces (model results) (see online version for colours)

Table 1

 1 SPH and classical approach comparison

		Lagrangian FE models	SPH cutting model
	Large deformation process	Adaptative remeshing algorithm	SPH meshless nature
	New free surfaces creation	Continuous remeshing and	Particles separation
		fracture model
	Contact	Friction Coulomb approach	Particles interactions

Table 2

 2 Cutting conditions

	Tool material	Tungsten carbides Grade H13A WC-6Co K20
	Rake angle	0
	Clearance angle	11
	Cutting speed	60 m/min
	Feed	0.3 mm

Table 3

 3 Cutting parameters

	Case	Feed	Rake angle	Edge radius
	1	50 μm	0°	20 μm
	2	200 μm	0°	20 μm
	3	250 μm	0°	20 μm
	Table 4 JC model parameters
	A	B	n	C	m
	496 MPa	310 MPa	0.3	0.0	1.2
	Source: Fang and Fronk (2007)

Table 6

 6 Validation case definition

	Mill	Radius R 0	Flutes number N f	Helix angle i 0	Feed f	Depth of cut Ap
	Ball-end	4 mm	2	30	0.15 mm/dent	2mm

Figure 16 3D model validation (see online version for colours)


Pierre Selva
email: pierre.selva@ct-medysysae.com

Olivier Cherrier
email: olivier.cherrier@isae.fr

Valérie Budinger $
email: valerie.budinger@isae.fr

Frédéric Lachaud
email: frederic.lachaud@isae.fr

Joseph Morlier
email: joseph.morlier@isae.fr

A. J. M. Ferreira (Editor

Keywords: Composite structures, monitoring, EM impedance, neural networks

This paper presents a structural health monitoring (SHM) method for in-situ damage detection and localization in carbon fibre reinforced plates (CFRP). The detection is achieved using the electromechanical impedance (EMI) technique employing piezoelectric transducers as high-frequency modal sensors. Numerical simulations based on the finite element method are carried out so as to simulate more than a hundred damage scenarios. Damage metrics are then used to quantify and detect changes between the electromechanical impedance spectrum of a pristine and damaged structure. The localization process relies on artificial neural networks (ANN) whose inputs are derived from a principal component analysis of the damage metrics. It is shown that the resulting ANN can be used as a tool to predict the in-plane position of a single damage in a laminated composite plate.

INTRODUCTION

Structural health monitoring employing smart materials is increasingly gaining popularity among engineering communities. In particular, the electromechanical impedance (EMI) based SHM technique possesses distinct advantages such as the ability to detect incipient damage, use of non-intrusive piezoelectric transducers and potentially low-cost applications. That is, the advantage of using piezoelectric sensors for damage detection resides in their highfrequency capability, which exceeds by orders of magnitudes the frequency capability of conventional modal analysis sensors. Thus, they are able to detect changes in the highfrequency structural dynamics at local scale which are directly associated with the presence of incipient damage.

Damage indicators derived from the measured electromechanical impedance are commonly used to either detect the presence of damage in any structure or provide information about damage localization in a one-dimensional structure. However, as soon as a two-dimensional structure is considered such indicators fail to furnish enough information for damage localization. Therefore, in the present paper, artificial neural networks are utilized in order to predict the in-plane damage position in a laminated composite plate.

PROPOSED METHOD

The use of neural networks requires a large set of data in order to learn a process as well as to improve their capacity to generalize and predict the process they have been trained for. In the present case, ANNs are trained to predict the position of a single damage in a composite plate. In order to generate a significant dataset relative to different damage localization, a coupled-field finite-element model of the EMI technique is first developed (figure 1). The resulting electrical impedance at the terminal sensors is then utilized as a damage indicator. Finally, damage metrics derived from the impedance spectrum are used as inputs to train, validate and test the ANN.

FINITE ELEMENT MODELING

The proposed method lies on a harmonic analysis of the structure under study. This structure must be modeled with accuracy since it is well-known that a simple modelling of the structure with perfect boundary conditions and standard material properties gives results far from measurements. The material properties of the model such as damping ratios, Young modulus, dielectric properties of the sensors are thus adjusted so as to match the experimental impedance spectrum of a pristine structure. The FE model is composed of composite shell elements as well as 3D piezoelectric elements.

LOCALISATION BY ARTIFICIAL NEURAL NETWORK

The localization process relies on an artificial neural network (ANN) whose inputs are derived from a principal component analysis of the damage metrics. First, the optimal number

Q i I ω = Q ˆ I V Z / = Z ˆ Nodal electric charge extraction ) ( ˆω i Q
of neurons to assign to the hidden layer is computed. Second, neural networks are trained and validated using random input vectors taken from the dataset. Finally, the trained neural networks are tested using samples partitioned from the main dataset. The testing data is not used in training in any way and hence provides an "out-of-sample' dataset to test the network on. This gives a sense of how well the network is doing.

CONCLUSION

The damage detection and localization method based on EMI and ANNs is tested on a carbon fibre reinforced plate. The experimental results show that ANNs can be used as a tool to predict the in-plane position of a single damage in a laminated composite plate taking into account an important database of damages scenarios created from finite-element simulations.

Figure 1 :

 1 Figure 1: Modeling principle of the EMI technique.


Ahmed Jhinaoui

Laurent Mevel

Joseph Morlier

A CUSUM Test with Sliding Reference for Ground Resonance Monitoring

Ground resonance is potentially destructive oscillations that may develop on helicopters rotors when the aircraft is on or near the ground. Therefore, this unstable phenomenon has to be detected before it occurs in order to be avoided by the pilot. To predict the zones of instability, works have generally relayed on off-line modal analysis of the helicopter model. Unfortunately, this off-line analysis is not sufficiently reliable. The subspace-based cumulative sum CUSUM test, able of on-line monitoring, is a good alternative which permits -at once-to avoid the system identification for each flight point and to have more robust detection, with reduced costs. In this paper, we describe an alternative test-with a moving reference this time-in order to kill wrong alarms or premature responses that are observed for fixed-reference tests. Numerical results reported herein are driven from simulation data.

INTRODUCTION

Ground resonance is a recurrent phenomenon of instability for helicopters. It is due to the coupling between the lagging motions of the rotor blades and the fuselage in-plane oscillations.

Since the works of Coleman and Feingold [START_REF] Coleman | Theory of self-excited mechanical oscillations of helicopter rotors with hinged blades[END_REF], in which a description of this phenomenon and a mechanical analysis were first given, numerous methods -analytic and numeric-have emerged and many papers have been published on the subject [START_REF] Eckert | Analytical and a numerical ground resonance analysis of a conventionally articulated main rotor helicopter[END_REF]. The main goal of these contributions was to give a sufficiently accurate mechanical modeling and, then, to determine the instability zone-that is the values of the rotor angular velocity at which the resonance may occur (i.e. the system is unstable when one or more of the damping coefficients become negative). The disadvantage of this off-line and deterministic modal analysis is that it does not take into account the uncertainties for complex systems like helicopters. In fact, the dynamical behavior of this class of systems is function of much randomness (structural uncertainties, loads, erosion, fatigue…) so that the margins of stability could be affected [START_REF] Khodaparast | Propagation of structural uncertainty to linear aeroelastic stability[END_REF]. The structure has, therefore, to be identified continuously in-flight and often tested on ground, which is costly in money and in time.

The subspace-based algorithm of detection offers an interesting alternative to deal with this problem. The main idea behind this method is to compute a criterion of instability at a stable reference and then, by some distance formula, determine when this distance is significantly different from zero. A comprehensive study of fault detection method can be found in [START_REF] Basseville | Detection of Abrupt Changes-Theory and Applications, Information and System Sciences Series[END_REF][START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF].

A statistical cumulative sum (CUSUM) test [START_REF] Mevel | Fast in-flight detection of flutter onset : A statistical approach[END_REF] can be build for this subspace approach, in order to track eventual changes recursively i.e. in real-time. In previous works [START_REF] Jhinaoui | Ground resonance monitoring for hinged-blades helicopters: statistical CUSUM approach[END_REF], the authors have investigated the capacities of this method to detect ground resonance. Results have shown that the test responds close to the instability but there still be a slight premature response before. This was predictable; the reference is taken far from the resonance region, so any slight change in the stability criterion (which is the value of the damping ratios) engenders a slight response of the test, and these responses are then cumulated in time, so that they could trigger a wrong alarm.

We describe herein an adaptive CUSUM test, which updates the reference recursively (sliding reference) and thus, kill all wrong responses. The paper is organized as follows: first, the principle of CUSUM subspace-based test is explained. Then, we give the analogous test for the case of sliding reference. Finally, both of two algorithms are applied to a simulation data of a helicopter with hingedblades rotor.

SUBSPACE-BASED DETECTION WITH SLIDING REFERENCE

The subspace-based fault detection and isolation has enjoyed some popularity since its introduction in the seventies, and has found its application in many fields such as civil engineering and aeronautics. These methods are derived from the subspace-based identification which is reminded below, in its covariance-driven version. For more extended explanation, one can refer to [START_REF] Van Overschee | Subspace Identification for Linear Systems, Theory, Implementation, Applications[END_REF].

Subspace Identification

Let consider the linear discrete system:

   + ⋅ = + ⋅ = + + k k k k k k v x H y w x F x 1 1 (1)
Where n x ℜ ∈ is the state vector, The number of sensors r is chosen so that it is inferior The classical subspace identification method consists in building the Hankel matrix filled with the output covariances. Then, from a well-known factorization of this matrix into a product of the observability matrix and the controllability matrix, one can deduce the eigenstructure of the system in [START_REF] Coleman | Theory of self-excited mechanical oscillations of helicopter rotors with hinged blades[END_REF].

Let cov H be the covariance-driven Hankel matrix of dimension qr r p × + ) 1 (

, with p and q the tail length of output data.

            + + + + = q p p p q q R R R R R R R R R L M O M M L L 2 1 1 3 2 2 1 cov H (2)
Where

R i = E (y k y k-i T
) is the correlation of the output data and E is the expectation operator. For a large number of data N, R i could be estimated by

T i k N i k k i y y N R - + = ∑ = 1 1 .

One can easily demonstrate that cov

H posses the factorization property:

q p C O H ⋅ = +1 cov
, with O and C, respectively, the observability and the Controllability matrices.

And:

(3)

The subspace spanned by the left part of Hankel decomposition, namely the observability matrix, contains all the information about the eigenstructure of the system. This matrix can be obtained from a Singular Values Decomposition (SVD) of cov H :

[ ] T V U U       Σ Σ = 0 1 0 1 cov H , 2 / 1 1 1 Σ = U O (4)
Then, to extract the matrices of transition and observation F and H, a least square minimization is made:

O O = F , with             = -1 p HF HF H M O ,             = p HF HF HF M 2 O (5)
Once F and H found, the eigenvalues and the observed eigenvectors (λ, φ λ ) are computed by resolving the equations:

The couple (λ, φ λ ) is the eigenstructure of the system. It is stacked into the vector

vec θ Λ   =   Φ  
where ᴧ is the vector whose elements are the eigenvalues λi and Φ is the matrix whose columns are the mode shapes φ i .

Subspace-Based Fault Detection

The fault detection consists in monitoring the eigenstructure and determining if any change

            = + p p FH FH F M 1 O λ λ λ φ ϕ φ λ H I F I F = = - = - , 0 ) λ ( , 0 ) det(
has occurred on it. For that, a reference state θ 0 and some distance, from this reference, called residual are defined.

That distance is chosen as the product between a left kernel S of the matrix of observability (or of cov H ) at the reference, and this matrix at the current state θ:

) ( ) ( 1 0 θ θ + ⋅ p T O S or ) ( ) ( , 1 cov 0 θ θ q p T H S + ⋅
, S is taken so that S T .S=Id. This distance is null when the current eigenstructure θ is close to θ 0 , and different from zero, else. We build then the residual below [START_REF] Mevel | Fast in-flight detection of flutter onset : A statistical approach[END_REF]: [START_REF] Jhinaoui | Ground resonance monitoring for hinged-blades helicopters: statistical CUSUM approach[END_REF] Where

T T p k T k k y y Y ) ... ( + + = and T T q k T k k y y Y ) ... ( 1 + - -=
and N the number of output data we have. This residual is not useful for real-time detecting of instability on helicopters. First, because it is computed once one has all the N data. And second, because for an aircraft, the change on an eigenvalue does not mean that there is no more stability. In fact, the eigenvalues of a helicopter changes all the time with the rotor angular velocity; the helicopter is unstable when one of the damping ratios is negative. The residual we have should then be expressed as function of the damping coefficients ρ i , and recursively for online detection.

) ) ( ( ) ( ) ( ) ( 0 0 1 0 - + - ⋅ ⋅ Σ ⋅ = k k T T k Y Y S vec J Z θ θ ρ ρ
, J and Σ are the sensitivity and the covariance of the residual. Consistent estimates of these matrices ∧ J and 1 ^-Σ are given in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF].

The CUSUM test to decide whether a change has occurred or no is (see [START_REF] Mevel | Fast in-flight detection of flutter onset : A statistical approach[END_REF][START_REF] Jhinaoui | Ground resonance monitoring for hinged-blades helicopters: statistical CUSUM approach[END_REF] for more details):

The two hypotheses to test are      ≈ 0 : 0 :

) ( 1 ) ( 0 0 0 f ρ ρ n n g g E H E H (8)
In order to have a response at the instability, the reference eigenstructure has to be taken close to the instable state. This is not possible for systems like helicopters, because it may lead to the destruction of the apparel. The reference is then taken far from resonance. In this case, any slight change in the damping would lead to a slight response of the test; these parasite responses are then cumulated and could trigger a false alarm. The modal analysis of a hinged blades helicopter given later in this paper will show that the damping ratios changes are not unimportant. So, the parasite responses would be significant.

CUSUM Test with Sliding Reference

The idea of tests with sliding reference derives from the adaptive algorithms which were investigated in some works [START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF][START_REF] Goethals | Recursive Output Only Subspace Identification for In-flight Flutter Monitoring[END_REF][START_REF] Zouari | An Adaptive Statistical Approach To Flutter Detection[END_REF]. These algorithms consist in subspace tracking by updating it,

) ) ( ( 1 )) ( ) ( ( 0 cov 0 T k p N q k k T T N Y Y S vec N H S vec N - - = + ∑ = ⋅ = θ θ θ ζ ) ( ) ( ) ( ) ( max ) ( ) ( ) ( ) ( 0 0 0 0 0 0 0 2 1 0 ρ ρ ρ ρ ρ ρ θ ρ n n n k p n k q n k p n q k k n Sum T g Sum T Z Sum - = = ⋅ Σ = - ≤ ≤ - = - ∑
using new coming data from sensors. The subspace which has to be tracked, in our case, is the left kernel S(θ 0 ). This kernel is now computed for a sliding reference θ n . The recursive residual Z k writes this time:

)

) ( ( ) ( ) ( ) ( 1 - + - ⋅ ⋅ Σ ⋅ = k k k T k k T k k Y Y S vec J Z θ θ ρ ρ
The kernel S k is updated with the IV-PAST method that is investigated in [START_REF] Gustafsson | Instrumental variable subspace tracking using projection approximation[END_REF][START_REF] Goethals | Recursive Output Only Subspace Identification for In-flight Flutter Monitoring[END_REF] and is computed for the sample k-l with k the current sample and l some time lag fixed by the user. The algorithm can be described as follows:

• Compute an initial kernel S 0 for some data tail N, (y q …y N-p ) • Then with IV-PAST, if we have S n ,(computed for data y from sample q+n-l-L to sample q+n-lp, L is the length of the sliding window) S n+1 is computed • J n and Σ n -1 are estimated as in the fixed-reference case but using S n in calculus, for the current sample n • The test is applied to compare the sliding reference at n-l and the current state n

The two hypotheses to test are      ≈ 0 : 0 :

) ( 1 ) ( 0 f n n n n g g E H E H ρ ρ (10)
The utility of this algorithm is shown below on the application to a helicopter simulation data. It indeed permits to kill any premature response, and only responds when the system becomes unstable.

HINGED BLADES HELICOPTER MODEL

We give herein the equations of motion for helicopter's ground resonance and a modal analysis of stability. Further mechanical explanations could be found in [START_REF] Jhinaoui | Ground resonance monitoring for hinged-blades helicopters: statistical CUSUM approach[END_REF]. The class of helicopters considered herein is the one with in-plane hinged blades rotors. The model below is known to be a sufficiently precise description for ground resonance studies [START_REF] Coleman | Theory of self-excited mechanical oscillations of helicopter rotors with hinged blades[END_REF].

Fig. 1. Helicopter's mechanical model-3 blades

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 x and K y , and two viscous dampers C x and C y as

n n n k p n k q n k p n q k k n Sum T g Sum T Z Sum - = = ⋅ Σ = - ≤ ≤ - = - ∑ max 2 1
illustrated in Fig. 1. The 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=1…Nb, with N b the number of blades.

State Vector Model

Applying the theorem of Lagrange to this model, one can demonstrate that [START_REF] Jhinaoui | Ground resonance monitoring for hinged-blades helicopters: statistical CUSUM approach[END_REF]: [START_REF] Zouari | An Adaptive Statistical Approach To Flutter Detection[END_REF] Where z=x+iy and η the Coleman [START_REF] Coleman | Theory of self-excited mechanical oscillations of helicopter rotors with hinged blades[END_REF] coordinate, such that:

0 ) Im( ) ( Re ) Im( ) ( Re ) Im( ) ( Re . . . . .. .. .. .. =                         ⋅ +                         ⋅ +                         ⋅ η η η η η η l y x K l y x D l y x M r r r
∑ - = + ⋅ = 1 0 ) 2 ( b b N k N k t i k b e N ib π ω β η
. The matrices M r , D r and K r write:

                      + + + + = 2 2 0 0 0 0 0 0 0 0 b I m m b I m m m N M m m N M m M z z b b r                       + - + = 2 2 2 2 ) ( 2 0 0 ) ( 2 0 0 0 0 0 0 0 0 b C b I m b I m b C N C N C D z z b y b x r β β ω ω                       - - + - - - + = ) 1 ( 0 0 ) 1 ( 0 0 0 0 0 0 0 0 2 2 2 2
( . ) ( t X C t Y =
, where:

      ⋅ - ⋅ - = - - r r r r D M K M I A 1 1 0 ) (ω , we choose ] 0 [I C = .
Sampling the system at some rate τ yields the equations:

k k k k k k w X C Y v X F X + = + = + + . ). ( 1 1

ω

, v and w are Gaussian noises with zero mean, and F=e τA (12)

Modal Analysis

The criterion of instability here is the drop of one or more damping ratios ρ i from positive values to negative ones (ρ i <0). Given the matrix A, variation of these coefficients with the angular velocity can be plotted. They are computed using the eigenvalues λ i of A:

The plot is shown in Fig. 2 with the structural properties reported in Table 1. Fig. 2 illustrates that the damping coefficients become positive and they are varying smoothly with rotor's angular velocity. The second modal damping coefficient ρ 2 (mode 2) changes from positive values to negative ones, from ω=1.62 rad/s to ω =2 rad/s, which is a criterion for ground resonance. It is this damping coefficient which will be monitored by CUSUM test.

A question to ask would be: why using the CUSUM test if we have the interval of instability from the modal analysis of the mechanical model? The answer is that this analysis is based on a deterministic and off-line method and does not take that may affect the margins of stability. In instability are highly sensitive to small changes flutter; the use of a statistical approach w Fig. 2. Damping coefficients vs. blades angular velocity

Simulation Results for CUSUM Test

To test the performances of the sliding approach to the fixed model above is simulated at a rate τ state of reference) to ω=2 rad/s by a step of 0.01 rad/s. For are simulated. The fixed reference test responds for zoom shows that the response occurs before. This was predictable shown in Fig. 2) is continuously dropping; any change on its value engenders a response of the test.

The response of the sliding test proves that, indeed, this approach kills any parasite response and the alarm is only triggered on near the instability, and not far from it.

CONCLUSION

The problem of detecting the ground resonance is addressed. An adaptive a order to perform the response time and tested with simulation data. Future works encompasses a generalization of these methods to anisotropic blades that leads to a ACKNOWLEDGEMENTS Zoom method and does not take into consideration randomness of stability. In [START_REF] Khodaparast | Propagation of structural uncertainty to linear aeroelastic stability[END_REF] it is shown for the flutter phenomenon that airspeeds of highly sensitive to small changes in aircraft structure. Ground resonance is similar to approach will then make the detection robust.

Fig. 2. Damping coefficients vs. blades angular velocity

Simulation Results for CUSUM Test

To test the performances of the sliding approach to the fixed-reference one, the helicopter model above is simulated at a rate τ=0.02 s. The angular velocity varies from ω=1 rad/s (taken as the =2 rad/s by a step of 0.01 rad/s. For each value of velocity, 1000 output samples The fixed reference test responds for ω=1.6 rad/s, corresponding to the velocity of resonance, but a zoom shows that the response occurs before. This was predictable, the second damping coefficients, as shown in Fig. 2) is continuously dropping; any change on its value engenders a response of the test. response of the sliding test proves that, indeed, this approach kills any parasite response and the alarm is only triggered on near the instability, and not far from it.

The problem of detecting the ground resonance is addressed. An adaptive algorithm is proposed in order to perform the response time and tested with simulation data. Future works encompasses a generalization of these methods to a more complex model for helicopters rotors, which is the model of anisotropic blades that leads to a Linear Periodically Time-Varying system. consideration randomness and uncertainties phenomenon that airspeeds of Ground resonance is similar to reference one, the helicopter ω=1 rad/s (taken as the each value of velocity, 1000 output samples Fig. 4. Sliding reference CUSUM test response vs. sample corresponding to the velocity of resonance, but a , the second damping coefficients, as shown in Fig. 2) is continuously dropping; any change on its value engenders a response of the test.

response of the sliding test proves that, indeed, this approach kills any parasite response and the lgorithm is proposed in order to perform the response time and tested with simulation data. Future works encompasses a more complex model for helicopters rotors, which is the model of

  vector. The vectors w and v are two white Gaussian noises with zero means.

  be the state vector. The system can be written as a Linear Parameter-Varying (LPV): the equation of observation

)

Fig. 3 .

 3 Fig. 3. Fixed reference CUSUM test response vs. sample

Fig. 3 .Fig. 4 .

 34 Fig. 3. Fixed reference CUSUM test response vs. sample

Table 1 .

 1 Structural properties for hinged-blades helicopter with 4 blades

	Structure variable	Value(unit)
	Blade mass
	m	31.9 Kg
	Fuselage mass
	M	2902.9 Kg
	Blade stiffness
	Kβ	200 N/m
	Main LG stiffness
	Kx	3200 N/m
	Ky	3200 N/m
	Blade damping
	Cβ	15 N.s./m
	Main LG damping
	Cx	300 N.s/m
	Cy	300 N.s/m
	Lengths
	a	0.2 m
	b	2.5 m
	Moment of Inertia at articulation point
	Iz	259 Kg/m2

This work was supported by the European project FP7-NMP CPIP 213968-2 IRIS.


I Tawk

J Aubry

P Navarro

J.-F Ferrero

S Marguet

S Rivallant

S Lemaire

P Rauch

Study of impact on helicopter blade

Keywords:

This article presents a study of damage in structures that are similar to helicopter blade sections, subjected to an impact. These complex composite structures were impacted by a steel ball of 125 g at impact speed ranging from 30 to 130 m/s. This led to properly highlight the kinematics of the impact and to define the sequence of the damage's mechanisms. An explicit FE model is also presented. The damage modelling of the roving is performed through a scale change. It allows a good representation of observed experimental behaviour. As the mesh density is low, it can be used for the modelling of a real structure.

Introduction

To reduce certification and development costs, computational methods are required by the aeronautical industry that are able to predict the structural integrity of composite structures under impact from various bodies, such as birds, hailstones, small metal parts and other soft or hard objects.

The complexity of impact problems is increased when composite materials are involved, due to the distinct failure mode that may occur. The comprehensive review by Abrate [START_REF] Abrate | Impact on laminate composites[END_REF][START_REF] Abrate | Impact on laminate composites, recent advances[END_REF][START_REF] Abrate | Impact on composite structures[END_REF] discusses impact failure mechanisms and summarises impact modelling approaches, based mainly on analytical models.

Recent advances in understanding the damage mechanisms of laminated composite [START_REF] Talreja | Damage mechanics of composite materials[END_REF][START_REF] Jonhson | Computational methods for predicting impact damage in composite structure[END_REF][START_REF] Hu | A new cohesive model for simulating delamination propagation in composite laminates under transverse loads[END_REF][START_REF] Tawk | Composite delamination modelling using a multi layered solid element[END_REF][START_REF] Tay | Progressive failure analysis of composites[END_REF][START_REF] Fleming | Delamination modeling of composites for improved crash analysis[END_REF] coupled with development of the performance and capability of computers, offer the possibility of avoiding many experimental tests by using impact simulation. However, numerical results should be examined cautiously. They must, in addition, be validated by experimental tests. A key point in this research is the development of models that are usable in an industrial structure modelling.

The objective of this study is to propose a methodology of calculation that can be used by industry for the analysis of a rigid body impact on complex composite structures such as helicopter blades. Very few impact studies have been conducted on complex composite structures. They deal generally with soft body impact, such as gelatine, hail or ice [START_REF] Johnson | Modelling soft body impact on composites structures[END_REF][START_REF] Anghileri | A survey of numerical models for hail impact analysis using explicit finite element codes[END_REF][START_REF] Kim | Experimental investigation of high velocity ice impacts on woven carbon/epoxy composite panels[END_REF].

In this study, impact tests on structures similar to the blade sections are performed. They are carried out using a steel ball and for a wide range of speed. These tests allow the identification of the kinematic mechanisms of degradation. Following this, a model is proposed which represents the observed experimental phenomena. It relies on a multi scale approach. Some models of realistic damage that can be used on industrial structures are proposed. The numerical model is validated by comparison with experimental results.

Experimental test

As shown on the schema, the studied specimens are composed of a main spar with unidirectional glass-epoxy, a skin, a rib with hybrid composite glass-epoxy and carbon-epoxy, a polyurethane foam core, and a protective stainless steel that covers the front edge (see Fig. 1). The length of the specimen was 378.2 mm, and with a maximum thickness of 22 mm.

The tests were conducted using a gas gun. The projectile used was a steel ball with a diameter of 30 mm and a mass of 110 g. This steel ball was stabilized inside the gas gun by a polyurethane foam support. The impact speed was regulated by controlling the gas pressure. The specimen was positioned inside a safety box. A Fastcam Photron camera with high capacity was used to record the moment of impact. The recording speed used was 30,000 fps (frame per second). Fig. 2 summarises the impact characteristics and the energy balance established (impact energy, energy of the specimen and of the ball after impact, and energy absorbed by the specimen).

The energy absorbed by the specimen during the impact varied between 71% and 76% of the energy provided by the projectile. The time of contact between the projectile and the specimen varied between 0.15 ms and 2.35 ms. The analysis of the pictures obtained with the fast camera during a test can be used to define the kinematics of the impact. The photographs in Fig. 3 present the impact test for a specimen with 550 J. The time between two images was 6.66 Â 10 À5 s.

The projectile impacted the specimen and penetrated into the front edge. Firstly, the shockwave formed at the impact point and spread transversally. After contact with the rib, it then propagated along the length of both ends of the blade. This spread varied with impact energy and corresponded to the separation of the skins from the foam. The time that the spread took to arrive at the edges of the blade was 0.5 ms. The indentation caused a strong plastic deformation of the stainless steel protection which opened and lead to the separation of the skin from the foam. This phenomenon was due to a rupture of a small thickness of the foam bonded with the skin. The observed bending of the front edge decreased after 0.5 ms due to the elastic return of the structure. The experimental tests carried out allowed an understanding of the observable degradations on the test specimens after impact. The depth of the penetration, as well as the deformation of the front edge, was measured using a three-dimensional measurement machine (DEA-Hexagon Metrology). These measurements give the depth, as well as the length, of residual deformation (see Fig. 4). The behaviour is quasi linear for all tests. To quantify the zone of delamination, a tapping was carried out on the two faces of the test specimens. Fig. 5 shows the evolution of the delaminated area as the impact energy increased. In order to analyse the internal damage of the structure, the test specimens were stabilized by inclusion in a resin, and then cut every 10 mm in parallel to the symmetry plane. The observation of the cross sections located in the plane of impact allows the determination of the damage that appears after the degradation mechanism (see Fig. 6).

For an impact energy of 20 J, the damage generated in the structure is limited to a light staking under the impact point of stainless steel protection. When the speed of the projectile increases, it penetrates increasingly into the leading edge of the blade. This results in a plastic deformation of the stainless steel protection, as well as a greater degradation of the roving. Following this, the rupture of the skin/foam interface is initiated by the debonding of the interface between skin and the roving. It is necessary to reach an energy of 500 J to observe a separation between the stainless steel protection and the skin. This separation remains only in the zone impacted by the projectile and does not propagate longitudinally.

The skins are partly or completely torn behind the stainless steel protection throughout the length of the specimen. In this zone, to pass under the stainless steel protection, the skin undergoes a strong localised curvature that can have an influence on the behaviour observed.

For a relatively low energy level, the roving is fissured at the point of impact. For an impact energy of 200 J, the first cracks appear at 45°from the direction of impact. By increasing the impact velocity, the number of cracks in the resin increase while following the same angle. Packages of fibres will move according to this gliding plane. This kinematics of displacement is found for all the impacted blades. For an impact greater than 800 J, localised micro-cracking of the resin and displacement of fibre packages located on the front of the spar are observed. More cracks appear on the rear end of the roving, which causes the complete rupture of the spar. The foam in contact with the spar is also significantly damaged.

The damage observed for cross sections located slightly further from the plane of impact confirms these observations. For impacts greater than 300 J, the energy is sufficient to debond the skins of the roving up to the end of the test specimen. Micro-cracking of the resin and rupture of the fibre plies is limited to the zone located at the impact point. However, debonding can propagate longitudinally and can reach the ends of the specimens.

Numerical modelling

An F.E. model the impact on the blade's structure was developed. It must represent the physical phenomena observed in experimental tests. One difficulty results in the compromise between mesh sizing and convergence/calculation time. This can be very limiting for the modelling of the complete blade. The complexity of the phenomena that must be modelled and the scales distinct to which the phenomena of damage appear, force the use of very fine meshing to obtain an acceptable modelling, which involves prohibitive calculation times.

It is then important to set up effective strategies of modelling, with a macroscopic scale, which enable the modelling of these phenomena without the use of fine mesh. These developments were performed using an explicit code RADIOSS.

The main originality of the model is in the modelling of the roving. The model has to represent the damage of the matrix (micro-cracking and cracking) and the fragmentation in packages of fibres, whilst maintaining fibres with no damage.The idea was to model the matrix using 3D elements and the packages of fibres using 1D elements with a surface acceptable compared with the observed experimental results (/ = 0.9 mm). The fibre surface represents 65% of the spar surface, which is a typical percentage of fibre. During an impact, the matrix is damaged independently and the displacement of fibres is managed by slipping and contact between the 1D elements (see Fig. 7).

The impactor was modelled as a rigid body. The resin, as well as the stainless steel protection, was modelled by elasticplastic properties with isotropic damage. The composite skin was defined by a conventional property of orthotropic shells. The rupture was managed with strain. The modelling of the roving was carried out using a change scale. An equivalent Young modulus of the glass-fibre was obtained by applying the homogenisation laws. An elastic-plastic law, with damage, was used to model the resin. In order to accurately represent the behaviour of the fibres, being quasi linear until the rupture, a bilinear law with damage and with strain rupture was defined. The foam was modelled by a special law, which is able to differentiate between the compression and traction behaviour.

The numerical model (see Fig. 8) was made up of 104,810 nodes, consisting of 53,408 solid elements, 8819 shell elements and 20,509 bar elements. The length of the blade section with a grid was 280 mm. The definition of different interfaces allows the automatic definition of contacts.

Results and discussion

The analysis of results obtained from the model allows us to understand the kinematics of the impact. Fig. 9 shows the behaviour of the blade section during the impact for an impact energy of 550 J. The penetration of the front edge at the impact point is on the curve as a function of time.

The image at t = 0.06 ms shows the generation of the wave on the front edge caused by the contact with the ball. The amplitude of this wave increased and the wave propagated in the direction of the impact and towards the two ends of On the graph, which gives the penetration of the front edge as a function time, four zones can be observed:

Mass

À A first zone of load increasing during 0.35 ms, which corresponded to the local phase of damage. The penetration of the ball involved the rupture of the resin, the plastic deformation of the stainless steel protection, and then delamination between the skin and the roving which generated a swelling of the section. The maximum swelling was reached at approximately 0.38 ms. À A second zone where the penetration of the front edge stabilised for a duration of 0.25 ms. During this time, delamination propagated throughout the entire area between the front edge and the rib. À A third zone during the next 1 ms which represented the elastic return of the structure. À A fourth zone where the blade took its final form.

To validate the model, the results obtained numerically were compared with experimental measurements. Fig. 10 shows the speed of the blade after impact (a), the depth of penetration (b), the length of the residual deformation in the front edge (c), as well as the surface of the skin foam debonding as a function of impact energy (d). The results obtained are rather close. The maximum relative error for the blade speed after impact is 8.48%, for an impact energy of 854 J (V = 119 m/s). It is 5.68% for the penetration in the front edge, 4% for the length of residual deformation of the front edge and finally, 7.5% for the delamination area. These results show that the strategy used for the model is valid and that the model provides a good representation of the experimental ruptures for energy levels less than 1000 J. To test the robustness of the model, an impact of 520 J on a specimen with two plies of skin was modelled. The skin was made up of one ply of glass with 0°/90°orientation and one ply of carbon fabric with ±45°orientation. For the experimental test, this specimen was instrumented with strain gauges bonded to the structure in the axis of the impact, on the level of the spar, at the level of the rib and behind the rib.

Graphs (B, N, and D) shown in Fig. 11 presents the strains measured by the gauges and those calculated numerically as well as a zoom for the first 0.5 ms of the impact. This duration represents the time taken for the shock wave to reach the end of the specimen. For the three gauges, the numerical strain results are close to those measured in the experiments. The model reacted in the same way to that of the test specimen. At the beginning of the contact, the slope of the load is identical to that measured in the experiments.

The slope undergoes a break at t = 0.11 ms, which corresponds numerically to the beginning of the delamination at the skin-foam interface behind the spar. The elements of the skin-roving interface were broken and the resin was damaged locally under the impact point. The delamination of the skin-foam interface propagates by rupture of the interface elements, and reaches the rib. This phase corresponds to the unloading of the skin, which in the experiment, resulted in a major reduction in the strains measured by the gauges. This phenomenon is well represented numerically.

If one considers the strains measured behind the rib (D), the behaviour of the strains obtained numerically is very close to that measured in the experiment. When delamination reaches the rib, a decrease in the strains can be observed, which was also observed experimentally.

Upon impact, most of the energy is dissipated at the point of impact by plastic deformation of the stainless steel protection and by damage of the roving. Taking into account this damage is essential in the calculation result. The difficulty is to represent it in a finite element model of the complete structure. The developed strategy provides a solution to this problem. For the roving, the scale of the modelling is not the classical scale related to the material. It is related to the physical phenomena observed experimentally and used to represent the behaviour of the roving during an impact. This type of approach can be developed to represent different types of phenomena (buckling or damage) that appear on a smaller scale than that envisaged for the modelling of an industrial structure.

Impact energy (J) Numerical results

Experimental results

Impact energy (J)

Numerical results

Experimental results

Numerical results

Experimental results

Impact energy (J)

Conclusion

Many experimental tests were carried out to further understand the phenomena concerned during and after a frontal impact to a blade. From the observations, the mechanism of damage could be defined. Firstly, the front edge is the first element that undergoes damage, and then the rise of the impact energy level causes a plastic deformation of the stainless steel protection, which leads to the delamination skin-foam of the structure. Moreover, the penetration of the projectile damages the roving by cracking the resin, which causes the slipping of the packages of fibres on both sides of the impacted zone. A new strategy for the modelling of the impact on blades was proposed. This strategy consists of modelling the roving by the superposition of 3D and 1D elements. The 3D elements represent the epoxy's resin, while the 1D bars represent the packages of fibres. This approach was validated by the results obtained through the modelling of an impact on the blade specimen for various speeds. These results are satisfactory. The developed model allows a good representation of all the experimental ruptures for energy levels less than 1000 J.

The developed method precisely evaluates the energy dissipated by the damage to the roving. It allows this phenomenon to be taken into account in a scale model of the structure without greatly increasing the computing time.

This approach can be used industrially in the phase of pre-dimensioning to model the behaviour of the blade subjected to an impact. The developed model does not take into account the delamination of the two interfaces skin-rib and skin-roving. To manage these phenomena, a specific volume multi-layer element, dedicated to the modelling of the delamination, is currently in development within the team. [ms] 0 0,1 0,2 0,3 0,4 0,5 Fig. 11. Experimental and numerical strains in the axis of the impact versus time.

Fig. 1 .

 1 Fig.1. Section of the specimens and materials used.

Fig. 3 .

 3 Fig. 3. Successive photographs of an impact (67 ls).

  Residual penetration (b) Residual length deformationImpact energy (J) 0

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

 4567 Fig. 4. Observable degradation versus impact energy.

3D

Fig. 10 .

 10 Fig. 10. Numerical/experimental results.

  Fig. 8. Specimen mesh.

		T=0 ms				T 0 06 T=0.06 ms				T=0.16
		Stainless steel		i t f 3D interface elements l t		UD: Roving
		T=0.38 ms T 0.38 ms				Foam			Foam			T 0 68 ms T=0.68 ms
			UD: Roving g			Skin Skin	Rib Rib Penetration
		T=1 ms T=1 ms				T=1 4 ms T=1.4 ms					T=2 94 ms T=2.94 ms
	35 (m/s) 30						16	(mm)
	25							12
	20
	15							8
	5 10				Numerical results Experimental results	4
	0							0
	0	200	400		600	800	1000	0	200	400	600	800	1000	1200
		S i (a) Specimen speed after impact d ft i t			(b) Residual penetration p
	160	(mm)							(mm 2 )
	120
	80
	40
	0				Impact energy (J)
	0	200	400	600	800	1000	1200	0	200	400	600	800	1000	1200
		(c) Residual deformation				(d) Skin foam debonding

Fig. 9. Behaviour of the blade during the impact.