Christine Espinosa
email: christine.espinosa@isae.fr

Frédéric Lachaud

Jérome Limido

Jean-Luc Lacome

Antoine Bisson

Miguel Charlotte

J.-L Lacome

Coupling continuous damage and debris fragmentation for energy absorption prediction by cfrp structures during crushing

Keywords: Composite, Crushing, Damage mechanics, Computational modeling, Discrete particles, Energy absorption

Energy absorption during crushing is evaluated using a thermodynamic based continuum damage model inspired from the Matzenmiller-Lubliner-Taylors model. It was found that for crash-worthiness applications, it is necessary to couple the progressive ruin of the material to a representation of the matter openings and debris generation. Element kill technique (erosion) and/or cohesive elements are efficient but not predictive. A technique switching finite elements into discrete particles at rupture is used to create debris and accumulated mater during the crushing of the structure. Switching criteria are evaluated using the contribution of the different ruin modes in the damage evolution, energy absorption, and reaction force generation.

Introduction

New generations of aircraft structures are more and more using composite materials for primary structures, especially carbon fiber reinforced plastics (CFRP), thanks to their high specific energy absorption to low mass to strength ratio. For certification to crash without landing gears, it is necessary to design energy-absorbing structures that are linked to the fuselage both maximizing the absorbed energy and minimizing the total mass [START_REF] Wiggenraad | Finite element methodologies development to simulate the behaviour of composite fuselage structure and correlation with drop test[END_REF][START_REF] Mccarthy | Numerical investigation of a crash test of a composite helicopter subfloor structure[END_REF]. The quality of the optimum design of energy absorption structure in driven by two criteria: the structure should stay in a stable energy absorption process, and a maximum volume of material should be damaged during this process [START_REF] Farley | Crushing characteristics of continuous fiber-reinforced composite tubes[END_REF][START_REF] Hull | A unified approach to progressive crushing of fiberreinforced composite tubes[END_REF]. The design must take into account the complexity of loading conditions and the competition between structural long-dimensions behaviors and the local structural thickness-dimensions damages [START_REF] Qiao | Design optimization of fiber reinforced plastic composite shapes[END_REF][START_REF] Feraboli | Development of a corrugated test specimen for composite materials energy absorption[END_REF] (Fig. 1).

According to Farley and Jones [START_REF] Farley | Crushing characteristics of continuous fiber-reinforced composite tubes[END_REF], one can distinguish three structural modes of crushing for composite tubes reinforced with continuous fibers: transverse shearing, lamina bending, and local buckling. The energy absorbed by the structure during the process and its global resistance are measured by the specific energy absorption (SEA) and specific sustained crush stress (SSCS). The SEA corresponds to the energy required to crash 1 kg of material during the stabilized crushing phase. To link the stabilized crushing phase to the initiation of the crushing process, the Trigger ratio (TR) is another criterion that is used to qualify the structure design. It is computed as the ratio between the peak force (Fmax) and the mean force (Fmoy) during the stabilized crushing phase as presented in Eq. ( 1):

T R = F max F moy ( 1
)
This ratio is a macroscopic measure of the initiation quality. If we are to look at the initiation phase, we have to take into account the material degradation process and its effect on macroscopic forces redistributions. It is well known that composite materials suffer from five damage mode that are at the origin of a loss of mechanical strength and that are due to local loss of continuity of fibers (tensile rupture or buckling) or fiber/matrix links (debonding) inside the volume of the plies, Fig. 1 Left Composite fuselage analysis [START_REF] Wiggenraad | Finite element methodologies development to simulate the behaviour of composite fuselage structure and correlation with drop test[END_REF]. Right Load buckling of some fronds [START_REF] Feraboli | Development of a corrugated test specimen for composite materials energy absorption[END_REF] Fig. 2 Left Schematic representation of crash zone initiation [START_REF] Hull | A unified approach to progressive crushing of fiberreinforced composite tubes[END_REF]. Right Quasi-static observation of initiation with a chamfered sample [START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF] Fig. 3 Frond creation in the initiation phase for a chamfered sample plate [START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF] and from resin cracks through or between the plies (transverse cracking or delamination). Besides, it appears that for a same material (T700M21 CFRP for example), the crushing process produces different kinds of debris which can accumulate and create artificial and non-permanent inserts acting as oscillatory cutting blades during the process (Fig. 2). Indeed, initiation is the most critical phase [START_REF] Hull | A unified approach to progressive crushing of fiberreinforced composite tubes[END_REF][START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF] (Figs. 2,3), so that it is often necessary to both choose the adequate structure material couple, and also to design the initiator.

Nowadays crash-oriented design of vehicles or aeronautical structures or structural sub-components is realized through comparison of experimental and numerical testing on elementary geometries [START_REF] Mccarthy | Numerical investigation of a crash test of a composite helicopter subfloor structure[END_REF][START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF][START_REF] Davies | Impact on composite structures[END_REF][START_REF] Joosten | Experimental and numerical investigation of the crushing response of an open section composite energy absorbing element[END_REF][START_REF] Xiao | Modeling energy absorption with a damage mechanics based composite material model[END_REF][START_REF] Gregor | Simulation of progressive damage development in braided composite tubes under axial compression[END_REF]. Efforts have been made to represent material properties degradations, debris wedge generation, and residual shape of structures. To experimentally study the effect of the initiator for typical aeronautical components, Guillon [START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF] realized different tests and tried to link the observed phenomena with the registered load in time, and with macroscopic properties that we can extract of it.

Regarding numerical modelling, if it is desired to model both the material degradation and the structural efficiency in term of energy absorption, it is necessary to use methods that are able to represent both scales of behaviors and to calculate local as well as global quantities and fields.

As far as is known by the authors, numerical studies devoted to the analysis of the structural behavior and design of tubes or plates during crushing found in the literature are based upon the hypothesis of a 2D or multilayer shell behavior for the structure. Under this hypothesis, behavior laws for the plies material can be reduced to a plane stress or plane strain hypothesis and shell or 3D finite elements are used [START_REF] Joosten | Experimental and numerical investigation of the crushing response of an open section composite energy absorbing element[END_REF][START_REF] Gregor | Simulation of progressive damage development in braided composite tubes under axial compression[END_REF]. To numerically represent fractures creation and propagation, essentially two techniques are used: decohesion (loss of links between elements) and element kill (suppression of elements). Erosion by the element kill technique consists in eliminating from the computation the elements that suffers too large strains or stresses, simulating the local loss of mechanical strength [START_REF] Mccarthy | Numerical investigation of a crash test of a composite helicopter subfloor structure[END_REF][START_REF] Xiao | Modeling energy absorption with a damage mechanics based composite material model[END_REF]. The criterion is often based on a mean state of damage for composites which is directly linked to the mechanical strength. To represent the structure pealing [START_REF] Israr | Finite element simulation of 0/90 CFRP laminated plates subjected to crush-ing using a free-face-crushing concept[END_REF] between sub-lamina (often of different orientation only) or to represent wedge debris or fracture [START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF][START_REF] Joosten | Experimental and numerical investigation of the crushing response of an open section composite energy absorbing element[END_REF], a cohesive element is suppressed when the maximum stress is reached in the corresponding opening mode. Cohesive criteria can also be coupled with damage in sub-lamina [START_REF] Xiao | Modeling energy absorption with a damage mechanics based composite material model[END_REF][START_REF] Gregor | Simulation of progressive damage development in braided composite tubes under axial compression[END_REF][START_REF] Espinosa | A finite deformation continuum/discrete model for the description of fragmentation and damage in brittle materials[END_REF][START_REF] Guinard | A 3D damage analysis of low-velocity impacts on laminated composites[END_REF]. These studies are focused on the global properties of the structure and are able to create debris if the crack path is known a priori.

Other numerical methods exist to represent the generation of fragments, chips or debris from structures subjected to impulsive or high transient loadings, in particular the particlebased methods. The advantage of these methods is that they do not suffer from localized large deformations as do the mesh-based methods. The Lagrangian SPH particle method for example [START_REF] Liu | Smoothed Particle hydrodynamics: a meshfree particle method[END_REF] allows to use classical continuum mechanics material models, so that it is appropriate to self-generate a natural separation of matter and fragments as soon as the material reaches a certain state criteria defined in the frame of the classical continuum mechanics [START_REF] Limido | A new approach of high speed cutting modelling: SPH method[END_REF]. The SPH method is in particular well suited and robust enough to represent matter de-cohesion, holes, or matter accumulation due to shear and compression in ductile materials, since it is a continuous particle method (through the kernel approximation). But SPH becomes instable in tension and needs regularization if the ruptures/cracks that generate the fragments are due to a brittle (instantaneous stress relaxation after rupture) or semi brittle behavior (delayed stress relaxation process after process) [START_REF] Lacome | SPH formulation with Lagrangian Eulerian adaptive kernel[END_REF][START_REF] Michel | Hypervelocity impacts on thin brittle targets: experimental data and SPH simulations[END_REF]. Furthermore, SPH simulations require huge models and big computation efforts (time and space). To keep advantage of the continuous particle method approximation to represent large to very large and severe deformations up to rupture while limiting numerical instabilities, and allowing linking the particles to classical finite elements, Johnson and his wo-workers have developed the generalized particle algorithm (GPA) [START_REF] Johnson | A generalized particle algorithm for high velocity impact computations[END_REF][START_REF] Johnson | Numerical algorithms and material models for high-velocity impact computations[END_REF]. The method is able to deal with unequal spacing of particles and recent developments allow switching finite elements into SPH-like particles when a failure (or a non-admissible state) criterion is reached. Matter is then kept during fragmentation. The latter method is called CPEM for Combined Particle-element Method [START_REF] Johnson | A combined particleelement method for high-velocity impact computations[END_REF]. The method is indeed very efficient in modelling a perforation process and seems to capture very well the distribution of matter ejected from the impacted or from the bottom faces of a target impacted at a high velocity. As can be seen in [START_REF] Johnson | Another approach to a hybrid particle-finite element algorithm for high-velocity impact[END_REF], it is nevertheless still difficult to determine the residual size (mass) and velocity of the generated fragments, since fragments are generated by brittle or semi brittle (cohesive like) phenomena. To represent the size of fragments and distribution of matter, another particle based method seems to be useful, the material-point method (MPM). The method has been introduced by D. Sulsky et al. [START_REF] Sulsky | A particle method for historydependent materials[END_REF]. The MPM method is inspired from the particle-in-cell method used for fluid mechanics [START_REF] Sulsky | Application of a particlein-cell method to solid mechanics[END_REF][START_REF] Sulsky | MPM simulation of dynamic material failure with a decohesion constitutive model[END_REF][START_REF] Chen | An evaluation of the Material Point Method[END_REF]. It uses an Eulerian grid that computes the matter state variables (velocity, pressure, and internal energy) without any mesh distortion. The particles are Lagrangian material points which velocities and displacements are tracked through the grid and used to compute the particles state. Information exchanges on state variables between materialpoint particles are thus done in an indirect way using the grid as a support. This presents the advantage of stabilizing the computation (in particular the time step [START_REF] Ma | Comparison study of MPM and SPH in modeling hypervelocity impact problems[END_REF]) and allows the computation of non-cohesive material behaviors [START_REF] Bardenhagen | The materialpoint method for granular materials[END_REF] such as that of snow [START_REF] Stomakhin | A material point method for snow simulation[END_REF], seeds in a silo [START_REF] Wieckowski | The material point method in large strain engineering problems[END_REF], or soils [START_REF] Andersen | Material-Point Analysis of Large-Strain Problems: Modelling of Landsildes[END_REF][START_REF] Ikj | Formulation of a Dynamic Material Point Method (MPM) for Geomechanical Problems[END_REF]. Comparisons of simulations results obtained with a SPH model and a MPM model for a case of hyper velocity impact are presented by Zhang [START_REF] Ma | Comparison study of MPM and SPH in modeling hypervelocity impact problems[END_REF]. The method is more stable especially for large deformations due to high tensions, but needs a special treatment of cracks and contacts. As an example, Ambati et al. used the MPM method to simulate the orthogonal cutting process [START_REF] Ambati | Application of material point methods for cutting[END_REF]. The chip formation was well recovered. As in classical FE models, it was nevertheless necessary to set friction coefficients to get the right cutting force and the right shear bands formation, whereas the friction coefficient is an output in SPH simulations [START_REF] Limido | Metal cutting modelling SPH approach[END_REF]. The only need on SPH simulations is to set enough particles to catch strain localization, which again gives huge models and is not efficient in computation time because the time step becomes very small. For crack initiation and propagation, the MPM is attractive as it is possible to take into account the effect of a local failure on far-filed stress [START_REF] Schreyer | Modeling delamination a s a strong discontinuity with the material point method[END_REF] thus avoiding localization effects that can be observed for example when using cohesive elements in classical finite element simulations of impacts on composite structures. But the treatment of explicit cracks needs to set discontinuous velocities in the support grid [START_REF] Nairn | Material point method calculations with explicit cracks[END_REF][START_REF] Nairn | Analytical and numerical modeling of R curves for cracks with bridging zones[END_REF]. New developments have been done to compute stress concentrations in 3D structures containing a crack placed between the particles [START_REF] Guo | Three-dimensional dynamic fracture analysis using the material point method[END_REF], or multiple cracks supported by the particles for brittle materials [START_REF] Li | Modelling brittle impact failure of disc particles using material point method[END_REF]. Keeping in mind the idea that local stresses must be correctly computed to create and propagate cracks, and that the rupture can be supported by particles instead of being supported by surfaces between the particles, recent research activities have been devoted to the quality of numerical integration of meshess methods [START_REF] Babuska | Effect of numerical integration on meshless methods[END_REF]. These works show the necessity to link the computation of the large deformations due to plasticity or damage (up to rupture) in the continuum description of the material, and its representation through Lagrangian material points.

To model both the local material behavior and the structural fragmentation process a very interesting application similar to what can happen in the crushing process of our composite plates, is the case of fragmentation of structures under impulsive high energy loading, typically either hyper velocity impacts or explosions. Two strategies are retained here as the reference. The approach of Banerjee [START_REF] Banerjee | Material point method simulations of fragmenting cylinders[END_REF][START_REF] Banerjee | Simulation of impact and fragmentation with the material point method, computational physics[END_REF] consists in modelling the structure with discrete material-point. As a consequence, the structure is allowed to fragment following free paths from particles to particles that are set broken when their state variables reach continuum thermodynamically based criteria. The Eulerian grid that is used as a support for the material points is also used to compute the loading from the surrounding gas of explosion. On the opposite, Borvik and Wadley [START_REF] Borvik | A discrete particle approach to simulate the combined effect of blast and sand impact loading of steel plates[END_REF][START_REF] Wadley | Deformation and fracture of impulsively loaded sandwich panels[END_REF] use discrete particles for the gas and the sand clouds generated by the explosion and iso-geometric finite element with a nodal splitting technique to fragment the structure. The high order iso-geometric computation in the structure allows large deformations and the nodal splitting technique allows the structure to create fragments that have a minimum size under interest. This approach is more appropriate to our objectives and is easier to implement in the industrial finite element codes that are expected to be used.

Regarding the crushing process initiation, it is still a challenge to find the method or the methods that represent the fragment sizes, masses, velocities, and compressibility after it has been created because the debris itself becomes part of the crushing process. The computational strategy presented in this paper is a new methodology capable to simulate the different ruin and fragmentation modes during the crushing process of composite plates that self-initiates and propagates debris generation from local material degradations, and predict the global structural efficiency of the structure. A cor-Fig. 4 Experimental setup used by Guillon for experimental crush tests [START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF] puscular method that switches finite elements into discrete particles is used for small debris. The analysis presented in this paper is focused on the capabilities of the methods to reproduce the strong interactions between the physical ruin modes and to self both initiate and propagate fractures in a full 3D simulation. This is why the first part of this paper is devoted to the analysis of the experimental analysis of the crushing process it-self. The second part presents the particle method. The third one presents the strategy that has been adopted here to represent the continuous and discontinuous behaviors of a composite sample plate during crushing, and their ability to predict energy absorption and the fragmentation processes. The final part of this paper presents the simulation results and examines the added value of using discrete particles to adequately represent debris accumulation and their effect on the crushing process and energy consumption computation. It is concluded that discrete particle methods in conjunction with classical finite elements methods are the best suited combination to represent both the global structural behavior and the local fragmentation of composites while keeping a certain computational efficiency. efficiency.

Technical discussion

Experimental observations

Results under interest are experimental results available in [START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF] on T700GC/ M21e unidirectional laminates with numerical simulations of T800S/M21e laminated plates with the same layup [(0/45/90/-45/0)2]s. Samples are made of 268 g/m prepreg tapes (35 % of resin), so that the global density is 1535 kg/m 3 . The experimental drop test setup used by Guillon is presented on (Fig. 4). Samples are 160 mm (0 • fiber direction) × 60 mm (90 • perpendicular direction) chamfered rectangular plates so that initiation is conditioned by the wedge angle, here at 45 degree. Each ply is 0.27 mm thick and the laminate thickness is 5.4 mm.

The reference experimental results are presented in (Fig. 3). It can be seen that the initiation of crushing concerns the four first images. Initiation can be cut into several stages:

-0-contact: at the first beginning of the contact, the force grows very rapidly and the extremity of the plate is destroyed -1-first damage: the force has reached about half the maximum value and the outer 0 • ply delaminates; transverse shearing appears in lamina neighbors to the external 0 • ply and some fragments are created -2-peak: the force reaches the maximum value and about one half thickness of the wedge has been eroded, the plate is slightly flexed and touches the right horizontal guide; a long crack (delamination) appears at an interface situated near the middle thickness; the outer left lamina bend -3-full contact: the force is maintained when the external right ply touches the rigid base; the plate is slightly bend and touches the left horizontal guide; fragments of +45 and -45 lamina are created by transverse shearing; the middle 0 • lamina starts a local buckling; a lot of small debris are generated -4-two halves: the plate is separated into two half parts; the left part lamina bend and debris are generated in the confined region under the bending zones; the right half lamina starts bending on a long dimension rod between the rigid base and a damaged zone near the horizontal guide (fracture in the outer right 0 degree lamina); the force just dropped down from a value (point 3.5) and crossing the displacement measures on the pictures with the force curve, one can suggest that this is due to the fracture in the outer 0 degree lamina that has been created by the contact with the horizontal guide -5-stationary crush; one can see the big fracture in the outer right half plate sliding on the rigid base; fragments are generated in the two confined regions under the bending lamina with a bigger size than debris at stage 0 to 3.

It was concluded from the experiment observations of Guillon that the global behavior is the same in all cross sections of the sample [START_REF] Guillon | Etude des mécanismes dabsorption dénergie lors de lécrasement progressif de structures composites base de fibre de carbone[END_REF]. The whole structure and the heavy carriage fall together from a height of 1.5 m. No indication has been found on the effective carriage mass used in the so called UDFoG-DYNA-CH experimental tests.

A basis of five ruin modes

The risk of rupture of a composite structure is evaluated using structural criteria based on stresses (as the famous criteria proposed by Hashin [START_REF] Hashin | Failure criteria for unidirectional fiber composites[END_REF]), and rely on the hypothesis that the matter is continuous to compute the loading path up to the rupture. These criteria are discussed by different authors among whom the work of Davila et al. [START_REF] Davila | Failure criteria for FRP laminates[END_REF] or Brewer [START_REF] Brewer | Failure of graphite/epoxy induced by delamination[END_REF]. These criteria have in common the will to link the material elementary ruin modes and the structural rupture prediction. The question that is addressed in the present work is to find five basis ruin criteria of the structure that could drive coupling between inside plies continuous damage (in the matter volume) and discontinuous fragmentation of a CFRP structure during crushing. We propose here five independent ruin modes: three sublamina modes corresponding to damage in fibers, matrix and fiber matrix interface, one delamination mode (cohesive elements at the interface between lamina) and one interaction mode. These modes must be considered as an elementary basis from which more sophisticated modes can be derived.

As the nonlinear elastic anisotropic homogenized continuum damage mechanics material model has been developed based on the developments of Matzenmiller-Lubliner-Taylor [48], in the frame of the irreversible thermodynamics. Both models emphasize that the material remains elastic and keep its orthotropic directions after damage. Thus, as done classically, damage is modeled through its effect on the elastic rigidity loss in further loading or unloading until the damage reaches the value of 1 meaning rupture. The model has been generalized by Xiao et al. [START_REF] Xiao | Progressive damage and delamination in plain weave S-2 glass/SC-15 composites under quasi-static punch-shear loading[END_REF] for 3D woven composites under high energy impacts and by others such as Ilyas [START_REF] Ilyas | Damage modeling of carbon epoxy laminated composites submitted to impact loading[END_REF][START_REF] Ilyas | Modeling aeronautical composite laminates behavior under impact using a saturation damage and delamination continuous material model[END_REF]. The CODAM model developed by Mc Gregor [START_REF] Gregor | Simulation of progressive damage development in braided composite tubes under axial compression[END_REF] is inspired from the Matzenmiller CDM model and has been developed for braided composites in thin structures modeled with shell elements. The MAT162 material model available in LS-DYNA (http://www.lstc.com) uses the same behavior law formulation. The model proposed by Ilyas [START_REF] Ilyas | Damage modeling of carbon epoxy laminated composites submitted to impact loading[END_REF][START_REF] Ilyas | Modeling aeronautical composite laminates behavior under impact using a saturation damage and delamination continuous material model[END_REF] distinguishes six damage variables affecting the diagonal components of the homogenized compliance S of the ply. In Eq. ( 3), coefficients with a 0 superscript refer to undamaged material properties.

d i , i = 1.6 (2) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 11 22 33 12 13 23 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 E 0 11 (1-d 1 ) -ν 0 12 E 0 11 -ν 0 13 E 0 11 0 0 0 -ν 0 21 E 0 22 1 E 0 22 (1-d 2 ) -ν 0 23 E 0 22 0 0 0 -ν 0 31 E 0 33 -ν 0 32 E 0 33 1 E 0 33 (1-d 3 ) 0 0 0 0 0 0 1 2G 0 12 (1-d 4 ) 0 0 0 0 0 0 1 2G 0 13 (1-d 5 ) 0 0 0 0 0 0 1 2G 0 23 (1-d 6 ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ * ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ σ 11 σ 22 σ 33 σ 12 σ 13 σ 23 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (3)
In addition, the five failure modes are defined as follows:

f j , j = 1.5 ( 4 )
f 1 : tensile failure criterion:

< σ 11 > X T 2 + σ 2 12 + σ 2 13 S 2 f s -r 2 1 = 0 ( 5 )
f 2 : compressive failure criterion:

< -σ 11 + <-σ 22 -σ 33 > 2 > X C 2 -r 2 2 = 0 ( 6 )
f 3 : spherical crushing failure criterion:

< -(σ 11 + σ 22 + σ 33 ) > 3Z C 2 -r 2 3 = 0 ( 7 )
f 4 : transverse cracking criterion:

<σ 22 > Y T 2 + < -σ 22 > Y C 2 + σ 12 S 12 + < -σ 22 > tan(φ) 2 + σ 23 S 23 + < -σ 22 > tan(φ) 2 -r 2 4 = 0 ( 8 )
f 5 : delamination failure criterion:

< σ 33 > Y T 2 + σ 13 S 13 + < -σ 33 > tan(φ) 2 + σ 23 S 23 + < -σ 33 > tan(φ) 2 -r 2 5 = 0 ( 9 )
In these equations:

-<> defines the Macaulay brackets (positive part or ramp function) -X T and X C are respectively tensile/compressive failure stresses in fiber direction -Y T and Y C are respectively tensile/compressive failure stresses in transverse direction -Z T and Z C are respectively tensile/compressive failure stresses in out-of-plane direction -S f s is the debonding strength of fiber matrix interface -S 12 , S 23 and S 13 are respectively shear damage threshold stresses r j ∈ [1, +∞] is called the limit load ratio -The coefficient φ is given in order, first, to increase the failure shear stress when normal compressive stress occurs and secondly to represent friction

The limit load ratios r j are critical parameters since their evolutions will affect the d i damage variables in Eq. [START_REF] Gregor | Simulation of progressive damage development in braided composite tubes under axial compression[END_REF].They determine the volume of each space of admissible stresses but don't affect their shapes as we can see on (Fig. 5).

The effects of strain rate on failure/threshold stresses development are taken into account through Eq. [START_REF] Xiao | Modeling energy absorption with a damage mechanics based composite material model[END_REF], where C i and ˙ re f are material parameters, for strain rates higher than 400/s (that is the lowest strain rate than can be obtain with the experimental SHPB test equipment, see [START_REF] Ilyas | Damage modeling of carbon epoxy laminated composites submitted to impact loading[END_REF]):

σ strength dyn = σ strength stat 1 -C i ln ˙ ˙ re f ( 10
)
At the undamaged state, each limit load ratio r j is set to 1. Limit load ratios are increased to take into account the current state of damage using the first thermodynamics principle [48] (Fig. 5). Failure modes are supposed to be independent from each other, meaning that the initial situation where each limit load ratio is set to 1 is the only possibility to describe the undamaged material. Under this assumption, these failures modes must be seen as a basis derived from the one described in Sect. 2.1.When the current computed stress state goes beyond the admissible stress envelop of a failure mode, the corresponding r j is increased so that to scale back the stress to an admissible value (radial return onto the scalar function) [START_REF] Xiao | Progressive damage and delamination in plain weave S-2 glass/SC-15 composites under quasi-static punch-shear loading[END_REF] (Fig. 6).

The evolution functions φ j use the updated values of the limit ratios to compute the contribution of the corresponding failure modes to the different damage variables. Since the behavior is supposed to stay elastic, scaling down the stress on the admissible envelop and scaling up the damage induced by the evolution of r j must be coherent. It is easy then, using this condition, to derive several solutions for the evolution function (very similar to a generic form). In this study, we have chosen to use the following expression which is coherent with the assumptions of preserved elasticity and anisotropy after damage:

d i (σ, ˙ ) = min ⎛ ⎝ 5 j=1 q i j (m, r j ), d max i ⎞ ⎠ , φ j = 1 -e 1 m (1-r m j ) , d max i ≤ 1 (11)
The coupling matrix q i j is a part of the model and was originally set to:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d 1 d 2 d 3 d 4 d 5 d 6 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 1 1 0 0 0 0 1 1 0 0 0 1 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ * ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ φ 1 φ 2 φ 3 φ 4 φ 5 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ (12)
This is an improvement from what can be found in earlier models, in which the damage variable d 5 is not coupled with the limit load ratio r 5 . This modification has been done since the cracks inside unidirectional ply may not be normal to the plate surface, especially for low velocity impact loading. Moreover, it has been supposed that there is only one m parameter for all failure criteria, which allows an easier identification.

Effects of localized temperature are not taken into account in the behavior model during the process, neither on the evolution of material parameters such as the Young's modulus, nor on the stress-strain relations or ruin criteria. Indeed, the effect of temperature on the apparent Young's moduli of the unidirectional CFRP ply has been proven to be low (a few percent) because temperature affects essentially the resin behavior while the composite ply Young's moduli are essentially due to the carbon fibers resistance. Carbon fibers moduli are not affected by local temperature increase because the loading condition make them brake in compression or shear in a brittle way before the local temperature has reached a significant value for carbon [START_REF] Didierjean | tude du comportement des matriaux composites carbone/poxy en environnement hygrothermique[END_REF]. As a consequence, as for machining in aluminum work pieces, event though the local temperature could reach high values due to strain localization, it is dissipated in the fragmentation process and debris are evacuated rapidly so that the whole crushing process can be considered as adiabatic. In some cases of quasi-static loading such as a composite corner unfolding (L-Shape structures), it has been proven that it is necessary to take into account for the residual curing stresses in the identification of the rupture yield stresses Z T and Y T which appear to have different values [START_REF] Michel | Experimental and numerical investigation of delamination in curved-beam multidirectional laminated composite specimen[END_REF]. For the dynamic crushing process, these considerations were not taken into account.

Orthotropic damage driven fracture modeling

The key issue in fragmentation is to be able to transform a material behavior into a series of domains while keeping each domain a realistic matter resistance and kinematics state properties while releasing the tensile strengths at domains connections. A special attention is paid here at choosing appropriate numerical methods to disconnect domains and create smaller matter domains that correspond to physical debris or fragments. For computational efficiency, it is aimed to keep as long as possible in the simulation large 3D finite element domains to represent a continuous behavior for the composite structure.

As other commercial codes Impetus Afea Solver (http:// www.impetus-afea.com) offers different ways to create fractures. De-cohesion between domains is generated following different criteria: isotropic or anisotropic Cockcroft-Latham failure criterion, plastic strain failure criterion, Johnson-Cook failure criterion, geometric failure strain criterion. Elements can be eroded following the criteria, or they can be split and new nodes are created to allow matter opening. Among other features of Impetus Afea Solver, GPU computing allows much more rapid computation durations and makes it possible to use higher degree isogeometric interpolation finite elements. It has been chosen here to try different combinations of these methods with the previous described CDM material model to generate the different debris.

From the continuous volumes to particle-based small debris

To avoid sudden drops down and repeated shocks in the contact force between eroding parts in Impetus Afea Solver, ele-ments that should be killed can be switched into discrete rigid particles instead of being simply suppressed. When erosion is generated, each integration point of the element becomes a discrete particle. The rupture criterion must be reached at least at 25 % of the Gauss points for the element to be eroded which means two Gauss points among eight for order one fully integrated 3D elements. Finite element to particle switching is illustrated on (Fig. 7). The discrete particle-based method [START_REF] Olovsson | A particlebased approach to close-range blast loading[END_REF] has been developed to circumvent numerical difficulties arising from modelling large flows and contact of fluid or gas-like materials on constraining complex geometries of light and highly deformable envelopes (e.g. movable boundaries with wrinkles and close shapes). It is based on the kinematic molecular theory. Molecules are rigid sphere that follow the Newtons law of motion and collide, collisions being perfectly elastic. Particles represent several molecules, and are set as spheres whose radius can vary from a minimum radius R min to twice this value. The discrete particle method follows the method proposed by Cundall and Strack [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Their behavior follows a classical linear model as proposed by Salot [START_REF] Salot | Influence of relative density on granular materials behavior: DEM simulations of triaxial tests[END_REF]: an elastic and a damping interaction in the normal and tangent direction. Each particle created from an eroded element is given a rigidity k 0 , a mass m 0 and a radius R min , and a velocity. The chosen rigidity in our simulations was k 0 = 10 G Pa. The mass and velocity are automatically computed to conserve the kinetic energy of the former element. The mass is equally distributed over the particles and the velocity is computed using the iso-geometric interpolation functions and the barycentric coordinates of the Gauss points. For order one elements, the eight integration points become eight particles. Each particle is then given 1/8 of the mass of the element m e . The minimum radius is computed using the current density ρ e of the material and the computed mass of the discrete particles m 0 : Contact between particles follows the Coulomb law. The tangent force is limited to μ times the normal force with μ the coefficient of friction. Additionally all forces are equal to zero when the particles are not penetrating one another (Fig. 8). In our simulations, following the experimental observations of the reference case, μ was set to 10 %. ζ is the fraction of critical damping. It was computed from the previous values following the formula given on Fig. 8.

m 0 = m e 8 = 4 3 π R min ρ e ( 13
)

Cohesive failure

Tied contact with failure between adjacent originally discontinuous groups of elements is a feature available in Impetus Afea Solver as in other FE commercial codes. The failure option is equivalent to a tied contact that can break and is defined as a MERGE FAILURE in Impetus Afea Solver. We used the MERGE FAILURE COHESIVE option. Failure is generated using a classical quadratic failure criterion using mode I and mode II strength, and related G I C and G I I C energy release rates.

ζ σ σ n f ail 2 + ζ τ τ s f ail 2 ≥ 1 ( 1 4 )
where ζ is a scale factor accounting the inability to resolve stress concentrations at coarse element grids:

ζ = max 1, δ δ re f ( 15
)
δ is the local characteristic element size on the slave side of the merge interface. The stress unloading from failure is a linear function of the crack opening distance. It is defined such that the consumed energy per unit area of cracking G is:

G = σ σ n f ail G I 2 + τ τ s f ail G I I 2 (16)

Proposed strategy

The CDM equations have been implemented as user defined material models in LS-DYNA, SAMCEF, Abaqus, and Impetus Afea Solver. LS-DYNA was chosen at the beginning because the code is very stable and efficient for crash applications. Impacts and compression after impacts simulations have been done using our user defined material model with LS-DYNA [START_REF] Ilyas | Modeling aeronautical composite laminates behavior under impact using a saturation damage and delamination continuous material model[END_REF][START_REF] Chen | A numerical approach for analyzing post-impact behavior of composite laminate plate under in plane compression[END_REF], as well as L-shape stiffener quasi-static and dynamic unfolding [START_REF] Ch | Modélisation numérique de linitiation et de la propagation des dommages lors du dépliage de cornières composites aéronautiques Comparaison des prévisions théoriques[END_REF]. In these cases, the damage was completely modeled using the continuous through ply damage model. A first attempt to strongly couple through plies damage with interfaces openings has been done for an impact case of study, using a user version of SAMCEF code. The CDM model was used for damage in the plies and as the initiation criterion of neighbor cohesive elements using a user defined mixed mode cohesive model [START_REF] Lachaud | Impacts on fuselages: a trial in using numerical simulation to predict residual strength[END_REF]. The final state of damage after impact was very well predicted, as for the CDM model alone. The non-local coupling was proven to be efficient and the simulation was very fast (some minutes compared with hours for the CDM only). In the present work, we use the CDM model to damage the plies. We used the erosion technique that transforms finite elements into discrete rigid particles when a chosen cumulated continuous damage criterion is reached. Three different criteria have been tested in the presented work:

-Criterion I: the material in-plane shear modulus G 12 has dropped down to zero, which means that the damage variable d 4 has reached 1:

d 4 = 1 (17)
-Criterion II: the in-plane E 11 fiber or E 22 matrix moduli have dropped down to their minimal values (not zero in axial fiber direction to keep a residual compressive strength) :

d 1 = d max 1 or d 2 = 1 ( 18
)
-Criterion III: the limit load ratios for matrix cracking and delamination have reached their limit value (see [START_REF] Ilyas | Damage modeling of carbon epoxy laminated composites submitted to impact loading[END_REF] for details)

φ 4 = 1 and φ 5 = 1 ( 19
)
Three computational strategies are used to create fractures: -A progressive finite element to particle switching is controlled by the damage state computed with the 3D continuous damage material model only (CDM) -Delamination is created using a cohesive failure contact between the damageable plies while damage in plies will not create fracture -Both cohesive delamination and finite element damage driven erosion are simultaneous. Two kinds of soft coupling are evaluated, each using a different criterion for FE to particles switching.

The four numerical simulation cases presented in this paper are summarized in (Table 1).

Numerical models and input data

Wedged plate crush FE model

The numerical model presented here is set in the framework of plane strain hypothesis, defining symmetry constraints on the lateral boundaries of the 3D one element in depth finite element model. In the simulation presented here, a value of 72 kg was chosen upon the possible values. In order to reduce the model size, the length of the plate in the model was shorter than in the experiment that is 100 mm instead of 160 mm. The free height is fixed to 20 mm in the model.

The finite element model does not always contain specific cohesive elements. When used, the cohesive elements are classical ones that are available in the commercial code. We have limited our work to elements of order one, which allows Fig. 10 Orientation of fiber directions us to keep a mesh that is fine enough for the observation of both material and structural failure modes, with a relatively good compromise between precision and computation time. Observations and analyses presented here are made from one lateral face as presented on (Fig. 9).

Finite elements for the composite structure are eight-nodes hexagonal fully integrated elements, with a characteristic size of 0.27 mm. To generate the chamfer, lines of elements in the stacking plies are one element shorter from top to bottom along the z direction of the frame, creating a scale shape (Fig. 10).

The composite structure is the only deformable part of the model. A flexural modulus is introduced in place of the Young axial modulus. The reader is invited to read reference [START_REF] Ilyas | Damage modeling of carbon epoxy laminated composites submitted to impact loading[END_REF] for further details. The heavy carriage is perfectly tied at the right end of the composite structure. This system has an initial velocity of 5.4 m/s. The mass of the carriage has been scaled down to the appropriate value for a 0.27 mm depth model, and represent a physical mass of 72 kg. Contacts are defined between the composite structure and the rigid guides, and between the structure and the rigid base. A static friction coefficient of 10 % is introduced on the rigid base, and 8 % on the rigid lateral guides. Faces of elements perpendicular to the y direction are constrained with symmetry conditions to obtain a 2D plane strain model.

Inputs for the CDM material model and the cohesive failure

The identification procedure requires a three steps experimental campaign: = 0.87 corresponding to a maximal density of cracks observed in experiments. Note also that damage on the fiber direction cannot go higher than a limit value so that a residual compression will always be preserved.

Even if smaller values of m in Eq. ( 11) (m = 1 or m = 4) resulted in a better stress-strain curve approximation, m = 10 has been chosen to obtain the correct maximum strain at rupture in dynamic compression (see Fig. 11). Then the model is considered to have less strength than the real T800S/M21e. One can also notice on (Fig. 12) that the T700GC/M21e presents a quasi-static tensile behavior different than the T800S/M21e, similar to hardening of metallic materials.

Table 2 provides the input data for the CDM material model parameters. Indeed, the brittleness coefficients m used in Eq. ( 11) are very sensitive to small variations and the SHPB characterization tests were difficult to realize on composite samples. Critical strain energy release rate calculations are obtained with quasi-static tests. The initiation values of 450 J/m 2 and propagation values of 800 J/m 2 , for G I C , have been reported by Prombut [61] for the T700GC/M21e. These values are in close comparison other test campaigns. A value of 765 J/m 2 was chosen for GIC which gives the admissible values of 50-65 mm for the crack length a. Validation simulations have been done with both implicit and explicit Matlab models in the framework of plane strain hypothesis, and with commercial codes (SAMEF, LS-DYNA, and ABAQUS) for 3D composite laminates. In order to obtain a numerical model that is valid for transient explicit computations of impact or crash, a double cantilever beam (DCB) test case from [START_REF] Pinho | Formulation and implementation of decohesion elements in an explicit finite element code[END_REF] was investigated. Analysis of the sensitivity of results on several numerical parameters e.g. damping, hourglass control, applied displacement and applied velocity, have been investigated and are presented in [START_REF] Ilyas | Dynamic delamination using cohesive finite elements[END_REF] (Table 3).

0 22 = E 0 33 = 7.64 GPa X T = 2.2 GPa X C = 1.2 GPa E f =

Comparison of experimental and numerical failure modes

Results of model 1 (CDM erosion) and model 2 (cohesive failure and CDM without fracture) are presented here. Figure 13 shows the first three stages computed with model 1.

The model does not show any of the behavior detailed in the analysis of the reference experimental case. It is too brittle. Analyzing the contribution of each limit load ratio in the evolution of damage variables d i shows that compression contribution computed by Eqs. ( 6) and ( 7) initiate this too brittle behavior. Slight oscillations appear in the global plate movements, but they are too small to create another shape of eroded area in the contact zone between the plate and the rigid base. Mechanical loss of resistance can be a reason of element erosion. That is why criterion II or III were chosen to erode the elements in the further models 3 and 4. To avoid instabilities due to the high level of compression at the real beginning of the crush, the coupling matrix has been changed in model 3 and 4: compression ruin modes driven by equations do not affect G 12 . It was also decided to adapt the model 4 to disable damage corresponding to cohesive failure that is damage d 3 , d 5 and d 6 are set to zero (Fig. 14).

Model 2 shows global and local buckling from the first beginning of crushing on the rigid base. The global flexural behavior of the plate that was observed in the experiment is also reproduced by the cohesive failure model 2. The problem of this model is that fully damaged finite elements stay in the simulation, and create too much instability in the crushing process. In order to get both the stable fragmentation process of the pure CDM model 1, and the local and global flexural behaviors of the plate obtained with the pure cohesive model 2, the two coupling models 3 and 4 are tested. Figure 15 provides results of successive erosion and delamination with model 3, and (Fig. 16) provides same results with model 4.

It can be seen that model 3 is too brittle again since elements are eroded on the basis of the loss of mechanical elasticity in the fiber direction and in the perpendicular in-plane direction. Erosion is more stable and the simulation run further for criterion III (E 22 , limit load ratio for perpendicular matrix cracking and limit load ratio for delamination reaching limit values). The three types of ruptures observed in tests are present: transverse shearing, lamina bending and local buckling. Furthermore the global opening modes and debris generation process is very similar to what is observed in the tests of Guillon. In model 4 the simulation reached the stabilized crushing mode as presented on (Fig. 17) where particles have been blanked for clarity of the pictures.

Model 4 is considered to be the most representative of what was observed by Guillon. 0 degree lamina is subjected to lamina bending or local buckling while other directions lamina is essentially subjected to transverse shearing. Debris of 45/90/-45 lamina are half the plate thickness long most of the time, which is characteristic of this ruin mode, and delamination about 0 degree lamina are stopped by the circumferential directions lamina. The simulation however did not exhibit the bending of outer right half part lamina as was observed in the experiment. Apart from the form of the failure criteria or limit values of the material, many reasons could be invoked to explain this difference: free space between the guides and the plate could be different in the model than in the experiment, the mass of the heavy carriage could be different as well, the rigid base could be not completely rigid during the experiment, and all conditions that are not perfectly controlled during tests and that must be rigorously defined as inputs in a computational model. Damage d 1 is almost not visible meaning that the elements that reach the limit value are eroded. Damages d 2 and d 4 reach slightly the same values at the same locations, meaning the strong interaction between ruin modes in the plies and at the interfaces, and a predominance of φ 4 and φ 5 on the loss of mechanical strength potential. Damage zones are visible in the outer lamina where the plate impacted the guides, as in the experiment.

Energy absorption and crushing forces prediction

Energy accumulated in elastic deformation and damage, called the internal damage energy E int , is computed at each time step n by summing the elastic deformation energy over all the N non-eroded finite elements in the composite plate (see Eq. 20). The elastic deformation energy of one non-eroded finite element at one time step is classically com- puted by the Gaussian quadrature of integral form of the elastic energy over the eight integration points. It uses the iso-geometric interpolation functions (of order one in our case) and the barycentric coordinates of the Gauss integration points (x; y; z) ∈ {-1/ √ 3; 1/ √ 3} 3 to obtain stresses and strains at each integration point. The internal damage energy versus time is reported for each of the four models on (Fig. 18). The energy of model 4 is taken as a reference for the comparison.

E n int =

N e=1 e σ n n dx dy dz [START_REF] Johnson | Numerical algorithms and material models for high-velocity impact computations[END_REF] One can isolate on the energy curve of model 4 three time periods each composed of a small unstable part and then a quite stable part between stages 1 and 4.5, stages 4.5 to 6 and stages 6 to 8. Oscillations are quite always 0.08 ms long, while stable parts are about 0.12 ms long. It can be seen that model 1 reaches quite instantaneously a value of about 4J at stage 1 of the crushing process, and remains at this plateau all the crushing long. This is coherent with the localization of fragmentation at the tip of the sample in the vicinity of the rigid support. The depth of fragmentation and the debris are so small that no oscillations are visible on the energy curve. Model 2 follows almost the same dissipation tendency than model 4 in the first time period and becomes unstable at stage 4. It is concluded that stage 4 is the point where matter erosion is necessary in the simulation, and the point at which coupling between in plane damage due to out of plane delamination becomes the controlling phenomenon. Model 3 also follows model 4 in the first time period, remains stable in the second time period, but absorbs less energy staying at a plateau level characteristic of localized fragmentation at the sample tip. The energy curves decreases when the process becomes unstable at stage 6. It is concluded that in plane The vertical displacement has been registered on the heavy carriage center of mass u carriage . For each case, the crushing force has been derived from the total amount of accumulated deformation energy E int during the crush using Eq. ( 21). The factor 60/0.27 is used to scale up the elementary force obtained with the model that is 0.27 mm thick, because the real thickness of the crushing plate is 60 mm. Curves for the four models are shown on (Fig. 19). It can be noticed on (Fig. 18) and on (Fig. 19) the correspondence with the stages described in Sect. 5. For comparison purposes, all time axes have been reduced to the first stages of the crushing process for which the simulation did not fail.

F cr ushing = 2 E int u carriage 60 0.27 (21)
As for the energy, the same stages can be observed on the crushing force. Model 4 exhibits both the effects of out of plane to in plane damage, and on a reciprocal phenomenon from damage to out of plane fracture. The peak crushing force computed by the model is about 41 kN (Fig. 20), while it is about 32 kN in the experiment. The mean crushing force computed by model 4 is about 10 to 6 kN in the stabilized parts, while it is about 10 to 5 kN in the experiment. Even though both are slightly overestimated, the computed pick and mean forces are of the same order of magnitude than the experiment ones.

Discussion

CDM manages distributed defects in volumes of materials by decreasing the mechanical elastic moduli. If used alone, the proposed CDM has proved its to represent the global behavior of wedge plates during crush if it is expected to represent the three types of ruin that are involved in experiments: transverse shearing, lamina bending and local buckling. The same limitation was found by [START_REF] Israr | Finite element simulation of 0/90 CFRP laminated plates subjected to crush-ing using a free-face-crushing concept[END_REF]. The solution that was found by Israr was to artificially stop the lateral deformations of 90 • plies in order to allow an axial crushing up to 90 or 95 % enabling the initiation of delamination in the neighboring cohesive elements. In our case, no numerical supplementary yields for strains or stresses where introduced. Using only cohesive failure is not adequate either since the real process also creates debris that are partially ejected and partially confined under the bending lamina. In the case of crushing plates, it is necessary to take into account for the re-distribution of loads between the laminas coming from decohesions and debris generation and accumulation, in a closed loop interaction. The proposed strategy that switches finite elements into particle allows to keep the low cost finite elements in the model as long as possible and to create matter separations and continuous debris which lengths are coherent with experimental observations of Farley and Jones [START_REF] Farley | Crushing characteristics of continuous fiber-reinforced composite tubes[END_REF]. Regarding the debris dimensions and accumulation during crushing, the proposed strategy is considered to be more adequate than the CPEM method proposed by Johnson [START_REF] Johnson | A combined particleelement method for high-velocity impact computations[END_REF] that gives either large clouds of small debris or even single particles, or large unbroken shrapnel. Results are comparable in quality with those obtained with the MPM method by Banerjee regarding the number of fragments [START_REF] Banerjee | Material point method simulations of fragmenting cylinders[END_REF]. The difference seems essentially to state in the necessity for the MPM method to model the entire surrounding environment in order to catch the right velocity fields whereas in our model only the structure is modelled. Macroscopic discontinuities that are clearly visible on pictures of experiments are strongly interacting with micro cracking. Here, a soft coupling model has been proposed that uses in plane damage to release the thermodynamic potential of the matter of the plies, and a mixed-mode cohesive failure for the interfaces. The in-plane damage model takes into account the load path redistribution due to through ply cracking and delamination rates, through a coupling matrix between five bases ruin modes and intrinsic damage variables, but in a passive way. If we compare with the very interesting work of Li [START_REF] Li | Modelling brittle impact failure of disc particles using material point method[END_REF] who uses the MPM method for brittle materials, it can be noticed that our model does not need any distribution default neither to initiate nor to propagate cracks, even if the coupling is weak. Indeed the transient dynamics introduces enough (and probably too much) differences in the integration points states of stresses and strains.

It is not really surprising that the computed force signal is not perfectly comparable with the one of Guillons experiments. We did not try to fit both of them in this study. The first reason of the discrepancies is that the materials are different and it is known that both materials do not have the same damage resistance under dynamic loading [START_REF] Ilyas | Damage modeling of carbon epoxy laminated composites submitted to impact loading[END_REF]. In this context, it is possible to compare qualitatively the behaviors but not quantitatively until a complete characterization of the T700GC/M21e has been done. Furthermore, the finite element model itself will create lower contact forces because even though finite element are switched into particles after erosion, the bulk resistance of the particle is not really the one of the real matter and higher values for the contact stiffness slow down the time step size. The purpose of this work was to analyze the energy absorption phenomena and the effect of coupling between fragmentations and debonding of the plies on the absorption. It has been shown that the global damage shapes are well qualitatively represented. The degree one for the finite element interpolation could not be enough to capture the stresses and strains gradients. Local bending is not well computed, and it is not possible for example to capture kink band phenomena which seem to be replaced by FE to particle switching instead. The 2D plane strain hypothesis is also very restrictive. It is then not necessary at this level of study, to fit exact failure criteria, especially when they are only weakly coupled as it was the case here.

Summary and conclusions

In this work we have used a continuous damage model coupled in a weak form with a cohesive failure model available in Impetus Afea Solver to compute the damage induced fracturing of a composite plate during edge crushing. The CDM model has been implemented as a user defined material model in different commercial computation codes. LS-DYNA was used to identify the material and damage parameters and with success to predict damage after low velocity impact or compression after impact residual strength. Samcef was used to strongly couple damage in plies and interface cohesive failure or low velocity impacts. Impetus Afea Solver has been used in this study to determine the combination of criteria that could be used to reproduce fragmentation and energy absorption of a composite plate during crushing. The best combination that has been found is use continuous damage to decrease in-plane elastic moduli of plies and use the cohesive failure to open interfaces using criteria of model 4. In model 4, finite element to particles switching is achieved when in plane damage is maximum. The best found criterion for switching was driven by the loss of thermodynamic internal potential of the matter due to ruin rates of through ply cracking and delamination. Using model 4, the three kinds of ruin modes and the stable evolution of the crushing force(displacement) and energy(time) are qualitatively reproduced. Next work will be oriented on strong coupling, and localization of load path to quantify the effect of local hinges or bending on the stability of the crushing process.

Fig. 5 2 Fig. 6

 526 Fig. 5 Tensile failure criterion hull with 2 parameters (stresses are in GPa)-the cylinder axis corresponds to σ 11 and the other axis to σ 12 and σ 13 -Left r 1 = 1; Center r 1 = 1.5; Right r 1 = 2

Fig. 7

 7 Fig. 7 Element erosion in Impetus Afea Solver-Left the bottom element e1 meets an obstacle; Center four of its Gauss points (in red) reach the criterion; Right the bottom element is eroded and all its Gauss points are converted into free particles (in green). (Color figure online)

Fig. 8

 8 Fig. 8 Contact between discrete particles

Fig. 9 Fig. 11 Fig. 12 Left

 91112 Fig. 9 Numerical model description: Left Dimensions; Right 3D plane strain finite element model

4 = 0. 87 S 10 ρ = 1 , 7 Table 3 Fig. 13

 4871017313 Fig. 13 Successive views of erosion at 0.09, 0.19, 0.39, and 0.78 ms with model 1 (CDM erosion)

Fig. 14 Fig. 15 Fig. 16 Fig. 17

 14151617 Fig. 14 Successive views of erosion at 0.2, 0.4, and 0.65 ms with model 2 (cohesive failure and CDM with no failure)

Fig. 18

 18 Fig. 18 Simulated absorbed energy versus time

Fig. 19

 19 Fig. 19 Simulated crushing force versus time

Fig. 20

 20 Fig. 20 Simulated force versus displacement from model 4

Table 1

 1 Numerical models presented

	Model	Continuous	CDM erosion	Cohesive
	identification	damage	criterion	failure
	number	mechanics	EF to particles
		(CDM)	switching
	1	X	Criterion I
	2	X		X
	3	X	Criterion II	X
	4	X	Criterion III	X

Table 2

 2 T800S/M21e input data for the model E 0 11 = 165 GPa E

Acknowledgments This work has been partially funded by the Région Midi-Pyrénée which is thanked here for support to the EPICEA project MODCOMP.