Rémi Delille

Thierry Garcon

Grégory Haugou

Bruno Laurent

Denis Lesueur

Hervé Morvan

Frédéric Robache

Catherine Foucart

FSDT First-order Shear Deformation Theory MSE Mean Square Error NCF Non Crimp Fabric ODM_MS Onera Damage Model MicroStructure RTM Resin Transfer Moulding RVE Representative Volume Element

Keywords: EHSDT Enhanced Higher-order Shear Deformation Theory Composite Materials, Modelling, Damage Mechanics, Viscoelasticity, Delamination

First of all, this present research would not have been possible without the collaboration between Toyota Motor Europe (TME), the Région Nord-Pas-de-Calais and the University of Valenciennes and Hainaut-Cambrésis (UVHC).

I especially would like to acknowledge my industrial supervisors, Tsukatada Matsumoto and Dr. Ernesto Mottola, for their

Global overview of cracks following an in-plane shear loading on various

fabric reinforced polymers by observation of the cross-section. . . . . . . . [START_REF] Gilat | Strain rate sensitivity of epoxy resin in tensile and shear loading[END_REF] 3.5 Visualisation of the stored strain effect on the continuity of the stress response. [START_REF] Greenhalgh | The assessment of novel materials and processes for the impact tolerant design of stiffened composite aerospace structures[END_REF] 3.6 Visualisation of the stored strain effect on the continuity of the stress response. 44

Introduction

In these last decades, the use of Carbon Fabric Reinforced Polymers (CFRP) in the transport market increased very significantly. The high specific stiffness and strength as well as the great energy absorption of these materials widely encourage their diffusion. Substituting metallic components by CFRP, the structure becomes lighter and leads to a significant performance rise.

Previously limited to small runs (aeronautic, racing), last advances in manufacturing cycle, such as reduction of the resin viscosity or curing time improvements, lead to the use of CFRP for high volume automotive production. The CFRP components can be classified as structural or non-structural. This distinction is coherent because of the inherent differences in the technical uses. A structural component is dedicated to carry loads and/or dissipate energy during a crash and have a thick layered section (about several centimetres). Rather, the non-structural parts, such as the engine bonnet or the roof, are dedicated to cover important surfaces and are most subjected to low-velocity impacts (tools, pedestrian, etc.), which provides the context for this study. By their nature, their section is thin (of the order of the millimetre) and made up of just few plies.

Layered fabric composite materials have two major advantages of being used in the automotive industry. The first one refers to manufacturing. The fabric preform being used as raw material (provided by a supplier) potentially shapable, it reduces the number of manufacturing operations. The second one is the enhancement of the damage tolerance under impact by stacking textile preforms as outer layer.

INTRODUCTION

But, the behaviour understanding and modelling becomes essential for their implementation into the design loop, needed for the deployment on mass-produced vehicles. In order to ensure the protection of pedestrians and drivers/passengers in case of collision with a CFRP panel, a material model for layered fabric composites dedicated to the finite element analysis (FEA) of impacts is needed. By the multi-scale physics which occur in a composite material, three approaches exist to describe the macroscopic behaviour.

The homogenisation methods use the knowledge of the micro-structure and the constituent properties to approach the macroscopic behaviour thanks to analytical and semi-analytical methods.

An other method consists of finite element based virtual testing, from the microscopic scale, to predict the behaviour of the upper scale, until the macroscopic one. Both these physical approaches have a main disadvantage for the use in the industrial sector. The knowledge of the physical phenomena which occur at the micro-scale warrants the precision of the model, but they are extremely hard to identify experimentally.

The last approach is more attractive for industrial applications because based on a phenomenological approach. The material is seen as an homogeneous medium and the micro-structural discontinuities like the constituents or the cracks are considered as diffuse. Thus, the framework of continuum mechanics is applicable and the parameters may arise from classical experiments.

Objectives and outlines

This work deals with the behaviour modelling of a layered carbon fabric reinforced polymer made up various fabric preforms under isothermal condition for impact-like loading. Experimental tests at several loading conditions and strain rates, set up by decomposing the loading conditions of impacts and by using standards, are carried out on each preform independently and on the layered material investigated. According to the experimental observations a single constitutive model able to describe independently the mechanical behaviour of each preform is developed. This is achieved by highlighting common phenomena and by adapting the evolution laws at the preform considered. This constitutive layer model is established within the framework of continuum mechanics and the material parameters are identified by using the initial experimental campaign and only few additional tests. But the layer stacking leads to specific phenomena due to interlaminar interactions. The element formulations actually implemented in the commercial finite element codes, by considering the material as homogeneous, does not provide a realistic strain distribution through the thickness of the material. And it largely conditions the bending stiffness and the damage prediction of the laminate. A realistic strain field in the laminate is thus computed in the material model by using the classical strain provided by the element formulation and the current stiffness of the layers. It is implemented in a commercial explicit finite element code used in the automotive industry. The accuracy of the present model is then evaluated by comparison with the responses predicted by finite element simulations and the responses obtained by experimental tests using a high-speed jack facility. Thus, this work may be considered as the combination of three main axis defined below:

Identification -To set up experimental test campaign able to provide the needed parameters and to define the identification procedure.

Modelling -To develop and implement a material model in the framework of the automotive industry which satisfies the requirement of low-speed impact simulations on the carbon fabric reinforced polymer.

Validation -To validate the material model by comparing finite element simulations and experimental data.

This thesis work is structured around the key phenomena which affect the behaviour of the investigated material. Each phenomenon is then mentioned by following the three main axis.

Chapter 2 introduces the material from the micro-structure to the stacking sequence. The definition of an impact and the description of the induced loading is given. From the literature review, the expected behaviour is presented by using the characteristics of every scale. In a second part, a brief state of the art on the modelling techniques, used in the industry to simulate impacts on layered carbon fibre reinforced plastic, is presented. This chapter draws the motivations of the modelling choices presented in this work.

Chapter 3 presents the basis of the material model, namely the matrix damage. Firstly a review of the different continuum damage models is made. Then the formulation of an existing model which is extended to better correlate with the material investigated is presented. The third section is dedicated to the parameter identification. Finally the present matrix damage model is validated by comparison between the simulation results and the experimental data.

Chapter 4 is devoted to taking the strain rate dependency of the matrix into account. The different existing methods are presented and the adopted solution is described. The next section provides the method and the experimental tests needed the identification. Then the original coupling between the previous matrix damage model and this newly presented viscoelastic model is described and validated thanks to high-speed tests on an hydraulic jack.

INTRODUCTION

Chapter 5 focuses on the final fracture of the material. It can be divided in two sections: the failure criterion followed by a massive degradation of the properties. A state of the art of the different techniques classically used is given. A smeared crack formulation is retained and developed to match the strain-space formulation. Following the parameter identification, numerical simulations are carried out to check the validity of the model.

Chapter 6 breaks the framework of the continuum mechanics used in the previous chapters. It describes the method developed to take into account the interlaminar damage and to enhance the bending stiffness prediction of the layered material. This part of the model can be seen as a sub-formulation of an element but with a strong dependence to the material model. It acts like an interface able to bound classical shell formulation and more realistic laminate theory by using current internal variables of the material model. Firstly a review of the different laminate theory is given. The theory of the present formulation is explained, followed by the description of the implementation. The identification procedure is explained and finally experimental investigations are carried out in order to validate the model.

Finally Chapter 7 gives a conclusion of this work. The manuscript ends with perspectives for further study and some possibilities of improvements.

Introduction of materials

The Carbon Fabric Reinforced Polymers (CFRP) belong to the range of the composite materials. These heterogeneous materials are made of at least two phases. The first, the matrix, is a continuous phase which is intended for maintaining in place the second phase, the reinforcement, to create a new heterogeneous material. In this study, the reinforcement enhances the mechanical properties of the matrix, whereas the matrix leads to a better distribution of stresses in the reinforcement. As a result, the mechanical properties of the composite material are different from its constituents independently examined.

The fibrous nature of the reinforcement makes it possible to adapt the process from the textile industry to produce initial preforms (see section 2.1.2). Moreover, the fabric plies can be superimposed to shape components to the desired thickness and to provide the required mechanical properties. It seems that the material has a multi-scale physiognomy defined as follows [START_REF] Melro | Analytical and numerical modelling of damage and fracture of advanced composites[END_REF]:

Micro-scale -This is the heterogeneous scale of the constituents (matrix and reinforcement) and of their interactions.

Meso-scale -This is the scale of a single ply of the laminate. Due to its architecture, a fabric ply have orthotropic mechanical properties in an homogeneous way. They can be determined thanks to experimental investigations or through numerical simulations on a Representative Volume Element. But this last method requires good knowledge of the microscopic phenomena.

Macro-scale -This is the scale of the laminate. The material is assumed homogeneous and the material properties are consequently an average of the in-situ physical phenomena. This represents the behaviour seen by the external environment.

Constituents

In the case of the investigated CFRP, the reinforcement consists in continuous carbon fibres surrounded by an epoxy resin.

Carbon fibre

The polyacrylonitrile (PAN) based carbon fibres (5-15 µm diameter) are classified following their mechanical tensile properties. The classification methodology is presented in Table 2.1 and summarised in Figure 2.1. Their various mechanical properties are due to the different organisations of graphene sheets after carbonisation and eventually graphitisation [START_REF] Dupupet | Fibres de carbone. Techniques de l'ingénieur Matériaux composites : présentation et renforts, base documentaire[END_REF]. Only fibres of the standard elastic modulus type (HT) are used in the materials which are taken into account in this work. Their mechanical properties provided by the suppliers are given in Table 2.2.

Due to the various oxydation rates accross the radius of carbon fibres, the latter are transversely isotropic. Moreover Kimura and Kubomura [START_REF] Kimura | Mechanical properties and applications of pitch-based carbon fiber reinforced plastics (CFRP)[END_REF] have explained the non-linear elasticity by the rotation of graphite crystallites toward the fibre direction. The carbon fibre is brittle and its tensile strength is around twice as high as its compressive strength [START_REF] Peebles | Carbon fibres: Structure and mechanical properties[END_REF].

According to the observations of Zhou et al. [START_REF] Zhou | Tensile behavior of carbon fiber bundles at different strain rates[END_REF] and Wang et al. [START_REF] Wang | Statistical analysis on high strain rate tensile strength of t700 carbon fiber[END_REF], the behaviour of carbon fibres is insensitive to strain rate variations.

Epoxy resin

In the case of CFRP, the matrix is selected from two main families of organic polymers: the thermoplastics or the thermosetting plastics. Once heated, the thermoplastics becomes mouldable through changes of intermolecular interactions, which are restored upon cooling. This allows the remoulding and the recycling of thermoplatics. Instead, the hardening of thermosetting plastics is irremediable. The cure process definitely set the internal structure to a solid state. However, the thermosetting resins are still used primarily for the manufacture of CFRP because of their higher mechanical and thermomechanical properties. The polyepoxides, commonly known as epoxy, resins belong to the range of thermosetting polymers and are used as the matrix component of the studied material.

The epoxy resin is obtained through the catalytic polymerisation of a monomer by mixing it with a curative followed by a heating phase. In this work, the epoxy comes from the polymerisation of the EPIKOTE ™ 05475 by using the curative EPIKURE ™ 05500 from the Hexion ™ (formerly known as Momentive Specialty Chemical) company. The limited amount of time needed for the curing of this epoxy resin is beneficial to the automotive industry and its mass production.

The tensile behaviour of epoxy resins is generally defined as brittle [START_REF] Schieffer | Modélisation multiéchelle du comportement thermo-mécanique des CMO et prise en compte des effets du vieillissement thermique[END_REF] even if a light ductilty appears for thin coupons [START_REF] Fiedler | Failure behavior of an epoxy matrix under different kinds of static loading[END_REF] as seen in Figure 2.2. In case of compression tests of neat resin cylindrical specimens, the stress-strain curves highlight a strain-softening after yield followed by a re-hardening (Figure 2.3). This strain softening can be explained by local shearing and frictional forces inside the specimen (which becomes barrel shaped) and do not represent the uni-axial compressive properties. Fiedler et al. [START_REF] Fiedler | Failure behavior of an epoxy matrix under different kinds of static loading[END_REF] and Gilat et al. [START_REF] Gilat | Strain rate sensitivity of epoxy resin in tensile and shear loading[END_REF] have observed for torsion tests a totally plastic deformation after reaching a maximum shear strength (Figure 2.4). However, Fiedler et al. have noticed that the fracture occurs along the plane perpendicular to the direction of the maximal tensile stress and is typical of the brittle-type fracture.

Furthermore, the epoxy resins show a strong dependence on rate of deformation [START_REF] Körber | Mechanical response of advanced composites under high strain rates[END_REF][START_REF] Gilat | Strain rate sensitivity of epoxy resin in tensile and shear loading[END_REF][START_REF] Buckley | Deformation of thermosetting resins at impact rates of strain. part i: Experimental study[END_REF][START_REF] Gerlach | Experimental characterisation and constitutive modelling of RTM-6 resin under impact loading[END_REF]. For every type of loading, an increasing strain rate leads to a reduction of the plastic flow and an increase of the stiffness. On a qualitative level, this involves an increase of the yield, flow stress, the maximum shear stress, the failure strength and the elastic modulus (tensile, compression and shear). Regarding the failure strain, the latter decreases for increasing tensile strain rates and increases with the compression state. These effects can be observed in Figures 2.2, 2.3 and 2.4.

However, identification of matrix parameters from neat resin coupons and intended for characterisation of Fabric Reinforced Polymer (FRP) must be handled with great care. According to Schieffer [START_REF] Schieffer | Modélisation multiéchelle du comportement thermo-mécanique des CMO et prise en compte des effets du vieillissement thermique[END_REF], the physico-chemical state of epoxy varies from neat resin curing to resin as matrix of FRP. Thus, a move back and forth between the material scales is impossible without a verification that the properties are similar.

Fabric preforms

In case of high-performance FRP, many processes have been developed to arrange and organise the fibres inside the plies and/or through the thickness of the laminate. Because of nature of the reinforcements, most processes are coming from the textile industry. The followings provide a few specific examples: filament winding, weaving, braiding or knitting. Each one has advantages and drawbacks and the choice is strongly guided by their specific mechanical properties, the production range and the costs.

To design the investigated non-structural component, two dimensional non-crimp fabric and woven textile preforms are used. Both, as simple and common fabric preforms, are cost attractive, easy to handle. In addition, each one presents benefits in their mechanical behaviours, as described in the following paragraphs.

Non-crimp fabric preform

Non-crimp fabric (NCF) preforms are the fabrics which most closely match a layered unidirectional architecture. Various layers of unidirectional web of continuous fibres are held together thanks to a stitching or a knitting through their thickness.

The non-crimp preform under investigation in this thesis is a biaxial carbon fabric [±45°] with glass fibre knitting (see Figure 2.5) and is manufactured by Saertex ® by using Pyrofil ™ TR30S carbon fibres. Its weight is 304 g m -2 , including 300 g m -2 of carbon fibres. By the nature of its manufacturing process, the straightness of the fibres in a NCF is generally maintained. As a result, its stiffness can approach that of two unidirectional plies in transverse directions [START_REF] Lomov | Carbon composites based on multiaxial multiply stitched preforms. Part 1. Geometry of the preform[END_REF][START_REF] Wang | Mechanical Properties of Stitched Multiaxial Fabric Reinforced Composites From Mannual Layup Process[END_REF]. The knitting slightly enhances the damage tolerance as well [START_REF] Dexter | Mechanical properties and damage tolerance of multiaxial warp-knit composites[END_REF].

Woven textile preforms

Compared to a non-crimp fabric, the integrity of the woven architecture is maintained without secondary material by the mechanical interweaving of yarns. A yarn is a clumping of fibres called warp in the longitudinal direction (stretched during the manufacturing process) and weft in the transverse direction (woven between the warps). Lots of patterns exist due to the endless possibilities of interweaving (including three dimensional) but in this thesis the woven preforms use the same one: plain-weave. This is the simplest and most common weaves: it is made by interlacing alternately the warp and the weft and hence has the higher number of intersections.

However, both plain-weave preforms used in this work differ from each other by the amount of fibres per yarn. The first one uses 3000 fibres per yarn for both warp and weft, and will be called 3K woven (Figure 2.6), while the second one uses 12000 fibres per yarn and will be logically called 12K woven (Figure 2.7). The weight of the 3K woven preform is 212 g m -2 but its fibre weight is only 200 g m -2 . The difference is due to the spraying of thin drops of epoxy resin at the surface to maintain the yarn position during handling. The drops have a diameter of 500 µm and are visible in Figure 2.6a. This preform is manufactured by Engineered Cramer Composites, a division of C. Cramer & Co, by using Aksaca ® A-42 carbon fibres.

The weight of the 12K woven preform is 193 g m -2 and this preform is manufactured by Hexcel ® by using Torayca ® T700S carbon fibres.

The crimp (or waviness) of the yarns declines the load-carrying capability and favours stress concentrations. As a result, the stiffness and the strength of composite layers made up with woven preforms are reduced compared to those of unidirectional layers [START_REF] Naik | Numerical modelling of woven fabric composite materials[END_REF]. However, the non-straight interfaces and the crossing of yarns contain the crack propagation. The shear damage tolerance and the delamination toughness, essential for the preform properties under impact, are hence enhanced [START_REF] Bibo | The role of reinforcement architecture on impact damage mechanisms and post-impact compression behaviour[END_REF][START_REF] Greenhalgh | The assessment of novel materials and processes for the impact tolerant design of stiffened composite aerospace structures[END_REF]. Moreover, the compliance provided by the waving lead to greater elastic deformation before damage initiation. Therefore the woven layers dissipate more energy than unidirectional layers when impacted [START_REF] Bibo | The role of reinforcement architecture on impact damage mechanisms and post-impact compression behaviour[END_REF][START_REF] Greenhalgh | The assessment of novel materials and processes for the impact tolerant design of stiffened composite aerospace structures[END_REF].

Laminated

The advantage of the layered fibre reinforced polymers is the ability to assemble and tailor layers to create a material with the desired engineering properties.

The considered laminated structure is made up with 12K woven fabric for the lower ply, non-crimp fabric for the two internal plies and 3K woven fabric for the upper ply (Figure 2.8). This choice was taken by Toyota Motor Europe (TME) and is explained by the mechanical properties of the preforms given section 2.1.2: woven outside to limit the impact damages and NCF inside to get the higher stiffness. The 3K woven is used for the visible layer of the component to improve the visual aspect and has the highest crimp, ensuring the best impact tolerance. The layer located on the inner side of the car should not be directly subject to impact and is not apparent. The 12K woven preform, cheaper and less crimped (so less impact tolerant but more stiff), is suitable for this last layer.

Due to the anisotropic properties of the FRP preforms it is essential to define the orientation of each one to approach the expected mechanical behaviour. A lot of nomenclatures exist to denote the stacking sequence but in this thesis the method is based on the notation of Nettles and Marshall [START_REF] Nettles | Basic mechanics of laminated composite plates[END_REF]. Firstly the longitudinal axis of the laminate, or simply called laminate direction, is defined as a line passing through the centre of the pressure sensor in the middle of the plate and the centre of the injection point (please refer to Figure 2.9). Then the layer orientations are defined between box brackets by providing the angle between the longitudinal axis of the preform, or simply called material direction, and the longitudinal axis of the laminate, from the outermost ply until the bottom ply, by using commas to separate them. Note that clockwise rotations are positive angles, and counterclockwise rotations are negative angles. The longitudinal axis of the NCF is set as the carbon fibre direction of the front side (Figure 2.5a, where the knitting shows a single yarn), while the longitudinal axes of the 3K and 12K woven preforms are given by the warp direction. As a result the stacking sequence of the structure investigated is (0°) 12K , (±45°) NCF , (0°) 3K .

It is also interesting to write the stacking sequence according to the directions of reinforcement and not only the direction of the preforms. In this way the stacking sequence nomenclature provides direct information about the expected mechanical properties. Moreover it is more easily comparable to the laminate made up with unidirectional layers. But a problem is the interlacing of reinforcement for the woven preforms. Therefore, a slash "/" will indicate two directions of reinforcement which are interlaced, whereas a comma "," will again indicate that the next ply is deeper than the previous one. By following this original notation the stacking sequence of the plate structure investigated is (0°/90°) 12K , (+45°, -45°) NCF , (-45°, +45°) NCF , (0°/90°) 3K . From that the basic configuration of a quasi-isotropic laminate structure is recognisable. It explains the choice of this stacking sequence and the result can be seen

Single preforms

For the needs of this work, the behaviour of the various preforms (NCF, 3K woven or 12K woven) have to be independently investigated. Therefore plates made up with single preform have been manufactured in order to identify their own behaviour.

The usual stacking sequence (using the longitudinal axis of the preforms) of these plates made up with single preform is given in Table 2.3, while Table 2.4 gives the stacking sequence according to the directions of reinforcement. In these tables, the index "S" denotes a "symmetric" stacking sequence, in which case only the half of the ply directions are specified (the others are obtained by symmetry). If the midplane of the laminate lies at the midplane of a ply (odd number of plies) the last specified angle is overlined. Moreover a subscripted number denotes a repeating unit in the laminate.

Plate type

Usual stacking sequence The number of plies is determined to obtained a plate thickness around 2 mm.

Structure (0°) 12K , (±45°) NCF , (0°) 3K NCF [(+45°) 3 ] S 3K woven (0°) 3 , 0° S 12K woven [(0°) 4 ] S

Curing process

The high production rates in the automotive industry require quick and automated manufacturing process. Based on these considerations, the resin transfer moulding (RTM) is the preferred process to produce wide plates.

The RTM process consists of injecting a low-viscosity thermosetting resin in a closed heated mould, in which the reinforcements were initially laid and compressed by the mould closure. Once the mould is completely filled and that the required pressure is reached, the curing process will continue by prolonged heating. The stabilised plate is then demoulded and a new cycle can start.

The plates were injected at the Transportation Research & Application Centre (TRAC) of Hexion ™ in Duisburg (Germany). The project being at a investigational stage, there is currently no automated process available to make the preforms. So all preforms were cut, stacked and laid in the mould by hand.

Last advances in RTM process at TRAC make possible a total cycle time lower than 180 s in which 120 s are dedicated to the curing itself. Details of the different cycles are given Table 2.5. Regarding the plates thereby obtained, a particular attention should be paid to the fibre volume ratio.

The intrinsic properties of the different plates are provided in Table 2.6. The geometric description of the plates is shown in Figure 2.9 but needs some explanations. Because of lacks of preform in the peripheral area of the plates, test samples cannot be cut no less than 35 mm of the edge of the plates. Additionally, due to the injection point and to the pressure sensor, circular excess thicknesses of resin are present in the centre of the plate. Test samples cannot be taken off in this area too. Consequently, among the 1681 cm 2 of a plate only 1300 cm 2 are still usable.

Impacts

An impact, a collision between a moving projectile and a fixed object, is classified according to the mass and the velocity of the projectile [START_REF] Bland | Observations on the impact behaviour of carbonfibre reinforced polymers for the qualitative validation of models[END_REF]. This distinction is carried out because of a transition in the mechanical behaviour from low-velocity to high-velocity impacts. Its physical explanation relies on the wave propagation in the laminate and the power provided by the projectile. When the laminate is subjected to long or short impact times (respectively Figure 2.10c and 2.10b), the plate absorbs the projectile energy by shell deformation. In case of long impact times the structure remains at the mechanical equilibrium, whereas for short impact times bending waves arise and lead to locally increase the bending. In these loading cases, the transverse shear strain induced by the bending may lead to matrix cracks inclined at approximately 45°for the inner plies and transverse matrix cracks for the external plies under tension [START_REF] Richardson | Review of low-velocity impact properties of composite materials[END_REF]. It has been proved that these intralaminar matrix cracks initiate the delamination at the ply interfaces [START_REF] Choi | A New Approach toward Understanding Damage Mechanisms and Mechanics of Laminated Composites Due to Low-Velocity Impact: Part II-Analysis[END_REF][START_REF] Kashtalyan | The effect of delaminations induced by transverse cracks and splits on stiffness properties of composite laminates[END_REF]. However, for very short impact times the delaminated area is wider and cannot be explained only in terms of structural bending. Bland and Dear [START_REF] Bland | Observations on the impact behaviour of carbonfibre reinforced polymers for the qualitative validation of models[END_REF] have shown that the dynamic effects can no longer be neglected. A compressive wave is transmitted through the thickness of the plate. But at each ply interface and due to an impedance loss this wave is separated into transmitted and reflected waves. It results a critical tensile loading of the resin at the interface level that could lead to a delamination and a mode I opening [START_REF] Bouvet | Low velocity impact modeling in composite laminates capturing permanent indentation[END_REF][START_REF] Chang | Study on impact damage in laminated composites[END_REF].

In this study, only the behaviour of fabric reinforced polymers subject to low-velocity impacts is investigated. Excepted in case of massive car crash, low-velocity impacts mainly corresponds to the type of dynamic loading to which a non-structural automotive component is subject. Note that the background of this work is the use of the material model for simulations of pedestrian impacts.

The impact properties of layered composite materials are well reviewed by Richardson [START_REF] Richardson | Review of low-velocity impact properties of composite materials[END_REF] and Cantwell [START_REF] Cantwell | The impact resistance of composite materialsa review[END_REF]. The nonlinear material behaviour which leads to differences in the impact response of composites is attributed to fibre failure, intra-and interlaminar matrix cracking, fibre-matrix debonding and strain rate sensitivity of the matrix (Figure 2.11). Additionally, the textile composite materials are capable of large shearing prior to the ultimate failure due to the sliding and reorientation of yarns. The modelling of all these phenomena is therefore essential to describe the impact behaviour of layered fabric composites.

Modelling

In structural analysis, the material models are used to predict the macroscopic behaviour of a component by approaching the mechanical behaviour of the material. At the sight of the multi-scale physiognomy of the fibre reinforced polymers, their behaviour models can be classified according to the smallest scale where a parameter is needed. Therefore, a model is set as microscopic in case of parameters describing the constituents behaviour, mesoscopic in case of need of data on the plies, or macroscopic if the parameters describe the laminated structure.

Microscopic models

The microscopic models are based on the knowledge of the constituents and of their interactions. Theoretically the properties of the structural component can thus be predicted, even in case of redesign (stacking sequence, fabric preform, number of fibres per yarn). However, the parameter identification is still a hot topic. Firstly it is not easy to determine microscopic properties from standard facilities. It requires a massive investment in cutting-edge technology. Moreover, the constituents on their own can have different properties than in a composite material. For example Schieffer [START_REF] Schieffer | Modélisation multiéchelle du comportement thermo-mécanique des CMO et prise en compte des effets du vieillissement thermique[END_REF] has shown that the degree of cure of the resin differs from bulk specimen or taken in a composite, and the property identification of the fibre/matrix interface or interphase has not been resolved. Finally, the computing time of these models are huge and not adapted to structural analysis.

Mesoscopic models

The mesoscopic models are based on the hypothesis that a ply in a laminated material is homogeneous and can be described thanks to the continuum mechanics. The ply properties can be determined through classical experimental tests (traction, compression, bending, shear) on layered material made up with the same preform and with a simple and unidirectional stacking sequence ([±θ ] p ). For fixed ply geometry (same thickness, same fibre volume ratio, same preform) it becomes possible to predict the behaviour of a laminate with a given stacking sequence. Moreover, at this scale, most of the physical phenomena can be introduced in the model formulation (matrix damage, fibre failure, delamination). Thus, the mesoscopic scale is able to approach the mechanical behaviour of a laminate, while proposing affordable computation times.

Macroscopic models

The macroscopic models are based on an hypothesis of a homogeneous material too, but at the laminate level. These models are purely phenomenological and are very simple to identify and do not require massive computation times. However, these models are absolutely not predictive and the slightest change in the design of the material or on the structural component leads to launch a new experimental campaign. Consequently, it becomes difficult to use these models for design applications.

Choice of the modelling scale

The choice of the modelling scale fell naturally on the mesoscopic scale. The present model is dedicated to be used in an industrial framework. Thus, the parameter identification and the numerical simulations should not require unaffordable and inaccessible facilities. Microscale analysis have thus been excluded. Finally, one of the aims of this project is to approach as close as possible physics behind the material behaviour. Since the mesoscopic scale is consistent in the consideration of the intra-and interlaminar cracks and the strain-rate sensitivity, most of the physical phenomena can be taken into account. In this work, these phenomena will be considered as independent as possible in order to split the development in different modules. A bridge between the different bricks of the material model is done as a second step.

Finite strain framework

Given the strong anisotropic behaviour induced by the reinforcement of the fabric composite plies and the potential large rotation of the yarns in case of shear loading, it is essential that the formulation may be able to follow the directions of anisotropy for large strains. As a result, the base (or undeformed) and the current (deformed) configurations of the material cannot be assumed identical. This hypothesis is the basis of the small displacement and small deformation theory, making it not applicable. Consequently, it is essential to extend the model in finite strain.

With few exceptions, Lagrangian (or Material) coordinates are used to describe the deformations of a solid. This approach facilitates the formulation of material constitutive models since the position and the physical properties of solid particles are described according to a reference position of these material particles and time.

As a result, it is possible to distinguish three major configurations: base C 0 , reference C R and current C (Figure 2.12).

The choice of the reference configuration to express the strain and the stress quantities leads to strongly different formulations for the model.

Expression of tensors

In order to formulate the constitutive model in a finite strain framework, it is appropriate to recall the definition of further used strain and stress tensors. Let ⃗ dX be a position vector which describes material points in the undeformed configuration. The material points in the deformed configurations are now described by ⃗ dx (Figure 2.13). The change of the material points is defined by the function:

dS 0 ⃗ N ⃗ M ⃗ d f 0 F dS ⃗ n ⃗ m ⃗ d f ⃗ N ⃗ dX = dl 0 • ⃗ N ⃗ M ⃗ δ X = δ l 0 • ⃗ M F ⃗ n ⃗ dx = dl •⃗ n ⃗ m ⃗ δ x = δ l • ⃗ m
x i = m i (X I ,t).
(2.1)

The differential is:

dx i = ∂ x i ∂ X J • dX J (2.2)
which leads to the definition of the deformation gradient tensor F:

⃗ dx = F • ⃗ dX with F iJ = ∂ x i ∂ X J . (2.3)
Let dl 0 be the length of a position vector in the ⃗ N direction and dl this length after deformation such as

⃗ dX = dl 0 • ⃗ N and ⃗ dx = dl •⃗ n (2.4)
with ⃗ N and ⃗ n unit vectors. From this the stretch of this position vector is defined as:

λ ( ⃗ N) = dl dl 0 (2.5)
and the elongation as:

ε( ⃗ N) = dl -dl 0 dl 0 . (2.6)
Lastly, the shear angle is determined by:

γ( ⃗ M, ⃗ N) = ( ⃗ M, ⃗ N) -( ⃗ m,⃗ n). (2.7)
From now, a distinction should be made between the tensors expressed according to the base configuration or the current configuration.

First of all, suppose that the reference is the base configuration. The scalar product of the material vectors in the current configuration, evaluated according to these material vectors in the base configuration, is given by:

⃗ dx. ⃗ δ x = dx i .δ x i = F iJ • dX J • F iK • δ X K = dX J • F T Ji • F iK • δ X K = ⃗ dX • F T • F • ⃗ δ X = ⃗ dX • C • ⃗ δ X (2.8)
where C = F T F is the right Cauchy-Green deformation tensor. Similarly, the variation of the scalar product according to the material vectors in the base configuration is defined with:

⃗ dx. ⃗ δ x -⃗ dX. ⃗ δ X = dX J • (C JK -δ JK ) • δ X K = 2 • ⃗ dX • E • ⃗ δ X (2.9)
where E = 1 2 (C -I) is called Green-Lagrange strain tensor.

Now that the different strain tensors are defined, the stresses have to be evaluated according to both configurations too. ⃗ d f is a force which acts on the current body and is the only one measurable from the experiments. By expressing the force vector according to the base or the current configuration, this leads to:

⃗ d f = T •⃗ n • dS ⃗ d f = Π Π Π • ⃗ N • dS 0 (2.10)
with T the Cauchy stress tensor and Π Π Π the first Piola-Kirchhoff stress tensor (PK1). Note that whereas T is symmetric, Π Π Π because of its mixed nature is not.

Let ⃗ d f 0 be a virtual force, seen as the equivalent of ⃗ d f which may act on the reference configuration. ⃗ d f 0 has no physical existence and is the transposition of ⃗ d f in the base configuration:

⃗ d f 0 = F -1 ⃗ d f . (2.11)
Hence, a new stress tensor integrally based on the reference configuration can be defined through the relation:

⃗ d f 0 = Σ Σ Σ • ⃗ N • dS 0 (2.12)
where

Σ Σ Σ = F -1 •Π Π Π (2.13)
with Σ Σ Σ the second Piola-Kirchhoff stress tensor (PK2) and Σ Σ Σ is symmetric.

Furthermore, the time derivative of the scalar product of the material vectors according to both base and current configurations are required.

By definition of the position vector in the current configuration, its time derivative is:

⃗ dx = Ḟ • ⃗ dX = ḞF -1 • ⃗ dx = L • ⃗ dx (2.14)
where L = ḞF -1 is called the spatial velocity gradient. Consequently the time derivative of the scalar product of the material vectors in the current configuration, evaluated according to these material vectors in the base configuration is given by: d dt

⃗ dx. ⃗ δ x = ⃗ dx • ⃗ δ x + ⃗ dx • ⃗ δ x = L • ⃗ dx • ⃗ δ x + ⃗ dx • L • ⃗ δ x = ⃗ dx • L + L T • ⃗ δ x = 2 • ⃗ dx • D • ⃗ δ x (2.15)
where D = 1 2 L + L T is called the rate of deformation tensor and is the symmetric part of the spatial velocity gradient.

Finally we can determine the time derivative of the scalar product of the material vectors in the base configuration, evaluated according to these material vectors by using the definition of the Green-Lagrange tensor:

d dt ⃗ dx. ⃗ δ x = d dt 2 • ⃗ dX • E • ⃗ δ X + ⃗ dX. ⃗ δ X = 2 • ⃗ dX • Ė • ⃗ δ X.
(2.16)

Iterative formulations

By means of the Finite Element Method (FEM), the final solution is obtained from the base configuration C 0 by stepping computations in an incremental solution process. At step n, a current configuration C n is computed based on the reference configuration C R . The choice of this reference configuration leads to different Lagrangian formulations.

Updated Lagrangian formulation

In the updated Lagrangian formulation, the reference is assumed to be the previous current configuration C n-1 and is usually updated after each incremental step. Therefore, the new current configuration C n is directly getting from the previous current configuration C n-1 (Figure 2.14). The reference configuration being at a deformed state, the behaviour law is formulated in terms of the conjugate stress-strain rate pair (T, D).

C 0 Base Configuration C n-1 Reference Configuration C n Current Configuration
But this scheme leads to additional computations to achieve the incremental objectivity. In addition, objective stress rate as the commonly employed Green-Naghdi or Jaumann, do not correctly follow the directions of anisotropy of the composite material [START_REF] Boisse | Meso-macro approach for composites forming simulation[END_REF][START_REF] Ten Thije | Large deformation simulation of anisotropic material using an updated lagrangian finite element method[END_REF]. Approaches using specific hypoelastic models where the fibre direction is tracked during the deformation were proposed [START_REF] Badel | Rate constitutive equations for computational analyses of textile composite reinforcement mechanical behaviour during forming[END_REF][START_REF] Badel | Large deformation analysis of fibrous materials using rate constitutive equations[END_REF][START_REF] Spencer | Theory of fabric-reinforced viscous fluids[END_REF][START_REF] Khan | Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. application to the double dome benchmark[END_REF]. However these methods by addition of specific manipulation on frames may proved to be long to implement and complex to identify.

Total Lagrangian formulation

Otherwise, in the total Lagrangian formulation, the reference is assumed to be the base configuration at each step. Therefore, the current configuration C n is always getting from the base configuration C 0 , regardless the current configuration computed at the previous step ( The reference configuration being at the undeformed state, the behaviour law is formulated in terms of the conjugate stress-strain rate pair (Σ Σ Σ, Ė). Therefore, the incremental objectivity is ensured by the use of total Lagrangian tensors. Moreover, by referring all calculations to the base configuration, the directions of anisotropy are well-assessed throughout the shearing.

However, a problem potentially arises for large deformation by using a total Lagrangian formulation link to the scheme of the small deformation theory:

σ σ σ = C : ε ε ε ⇔ Σ Σ Σ = C : E.
(2.17)

In case of constant components for the stiffness tensor, a softening behaviour occurs for moderate compressive loading. To overcome this limitation, extended hyperelastic models based on invariants of the strain tensors E or C are developed [START_REF] Bonet | A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations[END_REF][START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF][START_REF] Charmetant | Hyperelastic model for large deformation analyses of 3d interlock composite preforms[END_REF][START_REF] Peng | A simple anisotropic hyperelastic constitutive model for textile fabrics with application to forming simulation[END_REF]. But the present model is based on a model which is initially formulated by following the scheme of the small deformation. For obvious reasons of development time, the main body of the model given by the literature had to be kept. Accordingly the model validity should be checked for the operating environment.

A couple reasons ensure the possibility to simply extend the law which were formulated in small strains/small displacements to a total Lagrangian formulation. Firstly, the fabric reinforced plastics sustain moderate shear strains before failure. The range of the shear strain remains in the monotonic part of the force/displacement curve. Obviously, the strain in the fibre directions raises no concern, knowing that the longitudinal strains remain small. Secondly, the problem arises for shear strain but the model is intended to be used for impact loading. De facto such levels of shear will not, in practice, be reached as a result of multiaxial loading and of earlier failure.

Coupling

The formulation used by the finite element software, especially in the case of explicit analysis, differs from the total Lagrangian scheme chosen for the material model. Hence a link between the total and the updated Lagrangian frameworks need to be established. This is done by converting the second Piola-Kirchhoff stress tensor to the Cauchy stress tensor at the completion of the material model execution by using the classical push forward operation:

T = 1 J • F •Σ Σ Σ • F T (2.18)
with J = det (F) the jacobian determinant of F. Hence, the total Lagrangian formulation can be used inside the material model without influence on the iterative scheme used in the finite element code. For relatively large shear strains, the total Lagrangian framework is well adapted to the modelling of composite materials. Moreover the small strains along the fibre directions do not influence the formulation. Finally, because of the iterative scheme based on the constant initial frame, the objectivity and the consideration of the true material directions are ensured.

CHAPTER 3

Constitutive modelling of the intralaminar matrix damage

The Onera Damage Model Microstructure is used as first brick of the model for laminated fabric reinforced polymers. This constitutive model, formulated in the strain-space which ensures its efficiency for explicit finite element analysis, belongs to the wide family of continuum damage models. By considering privileged direction of intralaminar damage model, the model is able to consider the effect of the damage on the stiffness. In this thesis, the stored strains are considered as representative of the position of the crack lips. Thus, after calculation of the stresses applied to the crack lips, a friction law has been implemented in order to represents the hysteresis loops during cyclic loading. The present model is applied to various fabric preforms (woven or non-crimp) and, because of similarities of the physical phenomena which occur in each investigated materials, is able to represent a realistic behaviour. This model is also expressed in the total Lagrangian framework to ensure the coherence of the model with the material directions. Moreover, the possibility of a shear locking, very common among the textile simulations, is introduced with its effect on the out-of-plane diffuse damage.

Introduction

The present model is established within the framework of the Continuum Damage Mechanics (CDM). It was first introduced by Kachanov [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF] and Rabotnov [START_REF] Rabotnov | Creep rupture[END_REF] by considering the damage as a distributed defects through defining thermodynamic state variables. These variables are categorised as observable (or measurable) state variables -such as strains, stresses or temperature -or internal state variables (not directly measurable) -such as damage.

Then, the material law is based on two potentials:

• a thermodynamic potential, function of the state variables yielding the state laws;

• a dissipation potential, representative of the energy dissipation (from the damage, the friction).

Thereafter, Lemaitre [START_REF] Lemaitre | Evaluation of dissipation and damage in metals submitted to dynamic loading[END_REF] introduced the concept of equivalence principle which gave a physical interpretation of damage variables. This idea is based on the definition of an effective stress tensor. It may be interpreted as the stress leading to the same amount of deformation by replacing the damaged material by a hypothetical virgin one.

Originally set for isotropic materials, such as metals, the use of CDM for anisotropic and composite materials was introduced by Chaboche [START_REF] Chaboche | Anisotropic creep damage in the framework of continuum damage mechanics[END_REF] and Ladevèze [START_REF] Ladeveze | Damage modelling of the elementary ply for laminated composites[END_REF]. The damage variables are given in null-, second-or fourth-order tensor forms. In case of privileged damage directions, scalar variables are sufficient to well-describe the crack influences on the material behaviour. On the other hand, when the damage direction depends on the loading direction second-or fourth-order tensor forms are used.

Besides, it is important to take into account the unilateral character of damage. The closure of the crack, given by the stress state applied to it, leads to the recovery (partial or total) of the initial stiffness of the material. This unilateral character may be easily introduced thanks to the use of the Macaulay brackets for the stress normal to the crack orientation [START_REF] Ladeveze | Damage modelling of the elementary ply for laminated composites[END_REF][START_REF] Abisset | On the validation of a damage mesomodel for laminated composites by means of open-hole tensile tests on quasi-isotropic laminates[END_REF]. However, it leads to an incorrect behaviour in case of multi-axial loading [START_REF] Chaboche | Damage induced anisotropy: On the difficulties associated with the active/passive unilateral condition[END_REF]. Two approaches were then proposed by Chaboche:

• The first one consists of closing the diagonal terms of the stiffness tensor [START_REF] Chaboche | Development of Continuum Damage Mechanics for Elastic Solids Sustaining Anisotropic and Unilateral Damage[END_REF]. As a result, the initial shear stiffness is not recovered. Physically this may be seen as crack closure with a perfect slippage of the lips.

• The second approach [START_REF] Chaboche | A new micromechanics based CDM model and its application to CMC's[END_REF] leads to the complete recovery of the initial stiffness. By comparison to the previous one, it may be seen as an infinite friction between crack lips. In order to do so, an additional internal variable, the stored strain, is added to the model. It may be seen as representative of the position of the lips at closure. By definition it ensures the continuity at closure, but leads to discontinuities at re-opening when closure and opening do not occur at same loading configuration.

After a given loading/unloading then followed by a relaxation, permanent strains are observed and are imputed to the presence of damage. It can be explained by a release of residual stresses inherent to the manufacturing process and because of different thermal dilatations between the constituents [START_REF] Schieffer | A coupled analysis of mechanical behaviour and ageing for polymer-matrix composites[END_REF][START_REF] Schieffer | Modélisation multiéchelle du comportement thermo-mécanique des CMO et prise en compte des effets du vieillissement thermique[END_REF], but also by friction effects and microscopic plasticity of the matrix in the vicinity of the cracks [START_REF] Ladeveze | Damage modelling of the elementary ply for laminated composites[END_REF][START_REF] Abisset | On the validation of a damage mesomodel for laminated composites by means of open-hole tensile tests on quasi-isotropic laminates[END_REF][START_REF] Maimí | Matrix cracking and delamination in laminated composites. Part I: Ply constitutive law, first ply failure and onset of delamination[END_REF]. Schieffer et al. [START_REF] Schieffer | A coupled analysis of mechanical behaviour and ageing for polymer-matrix composites[END_REF] models these permanent strains by means of residual strains evolving linearly with damage, whereas Ladeveze and LeDantec [START_REF] Ladeveze | Damage modelling of the elementary ply for laminated composites[END_REF], Abisset et al. [START_REF] Abisset | On the validation of a damage mesomodel for laminated composites by means of open-hole tensile tests on quasi-isotropic laminates[END_REF], Maimí et al. [START_REF] Maimí | Matrix cracking and delamination in laminated composites. Part I: Ply constitutive law, first ply failure and onset of delamination[END_REF] use a plastic formulation.

Further works took interest in micromechanical considerations such as the diffuse matrix damage or the fibre/matrix decohesion. A crack density, based on the stress and an energy balance, is introduced to take into account the ply thickness on the kinetics of damage [START_REF] Ladevèze | On a damage mesomodel for laminates: micro-meso relationships, possibilities and limits[END_REF][START_REF] Ladevèze | Towards a bridge between the micro-and mesomechanics of delamination for laminated composites[END_REF]. Huchette [START_REF] Huchette | Sur la complémentarité des approches expérimentales et numériques pour la modélisation des mécanismes d'endommagement des composites stratifiés[END_REF] integrated this microscopic crack density in a viscoelastic formulation to provide the latter additional non-linearity.

Marcin [START_REF] Marcin | Modélisation du comportement, de l'endommagement et de la rupture de matériaux composites à renforts tissés pour le dimensionnement robuste de structures[END_REF] adapted and extended the micromechanics based CDM model proposed by Chaboche and Maire [START_REF] Chaboche | A new micromechanics based CDM model and its application to CMC's[END_REF] to the fabric reinforced materials. While remaining within the mesoscopic scale, he introduced a coupling between in-plane loading and damage normal to the thickness direction. Moreover, and in a similar way to Schieffer et al. [START_REF] Schieffer | A coupled analysis of mechanical behaviour and ageing for polymer-matrix composites[END_REF] and Huchette [START_REF] Huchette | Sur la complémentarité des approches expérimentales et numériques pour la modélisation des mécanismes d'endommagement des composites stratifiés[END_REF], a coupling between the matrix damage and a spectral viscoelastic model is performed.

For the present work, the constitutive relation given by Marcin is recalled Section 3.2. In addition, an adaptation of the stored strains to take into account the friction mechanisms inside the composite materials is formulated. The possibility of a shear locking for textile preforms is then added to the model. This shear locking leads to a damage mechanism which is added to the previous formulation. The implementation into a finite element code is described in the next section. Section 3.5 concerns the identification of the model parameters. Finally the model, and particularly the friction mechanisms, is validated through simulations of quasi-static cyclic in-plane shear tests.

Continuum matrix damage model

This work relies on a closed version of the Onera Damage Model MicroStructure (ODM_MS) proposed by Marcin [START_REF] Marcin | Modélisation du comportement, de l'endommagement et de la rupture de matériaux composites à renforts tissés pour le dimensionnement robuste de structures[END_REF]. Although the model has to describe the physical mechanisms, it is intended to be used in an industrial environment. Hence, the ODM_MS being based at the mesoscopic scale and formulated in the strain-space (ideal for an implementation for finite element analysis), it was well suited for basis of the complete material model.

The matrix damage modelling is based on the assumption that the preferred damage directions at the mesoscopic scale correspond to the directions of reinforcement. This Moreover, this model is able to depict the effects of the crack closure on the recovery of the initial stiffness. For this purpose the internal variable so-called stored strain, representative of the position of crack lips at closing, is used. However, as depicted by Chaboche and Maire [START_REF] Chaboche | A new micromechanics based CDM model and its application to CMC's[END_REF], the re-opening may lead to stress discontinuities that may be critical for the stability of finite element analysis. A further development was to introduce friction effects but was never followed up.

Yet, these friction mechanisms are fundamental for structural simulations with complex loading, notably in case of positive/negative shear switchover. Consequently, an efficient and Coulomb-based friction formulation is introduced in the present model and is described in the section 3.2.7.

Another important improvement is the consideration of the shear locking which is specific to textile preforms. The effect of the shear locking on the damage kinetics are also considered and is given Section 3.2.8.

Another important improvement in case of structural simulations is the use of the finite strain framework all along the modelling work. Thus, the total Lagrangian formulation is used to ensure the objectivity and the directions of anisotropy are well-assessed.

Constitutive relation

The model is formulated in strain-space to ensure a good efficiency for FEA by using the Helmotz free energy as thermodynamic potential:

ψ = 1 2ρ ψ m + ψ 0 -ψ r -ψ s (3.1) with          ψ m = (E -E 0 ) : C : (E -E 0 ) ψ 0 = (E -E 0 ) : C 0 : E 0 + E 0 : C 0 : (E -E 0 ) ψ r = (E -E 0 ) : C 0 : E r + E r : C 0 : (E -E 0 ) ψ s = (E -E 0 ) : C 0 : E s + E s : C 0 : (E -E 0 ) (3.2)
with the thermodynamic potentials ψ m related to the matrix damage, ψ 0 related to the initial state of the material, ψ r related to the residual strains induced by microscopic plasticity of the matrix and ψ s related to the stored strains representative of the position of the crack lips.

Then the constitutive relation is obtained by derivation of this thermodynamic potential and is given by:

Σ Σ Σ = ρ ∂ ψ ∂ E = C : (E -E 0 ) -C 0 : (E r + E s -E 0 ). (3.3)
Σ Σ Σ and E are respectively the second Piola-Kirchhoff stress tensor and the Green-Lagrange strain tensor (see section 2.4.1). C and C 0 are both fourth-order tensors which characterise the stiffness of the material. While C 0 denotes the elastic stiffness tensor of the material, C

represents the effective (damaged) stiffness tensor, evolving with the damage. E 0 represents the strain state where the cracks close off. It is due to the difference between the coefficients of thermal expansion of the matrix and of the reinforcement which creates residual stresses during the manufacturing. The residual strain E r and the stored strain E s are defined in section 3.2.6, 3.2.7 and 3.2.5.

Thereafter, the Voigt notations are used to denote the second Piola-Kirchhoff stress tensor and the Green-Lagrange strain tensor. Both tensors are indeed symmetric and can be reduced to first-order tensors in order to simplify the establishment of the constitutive relation. Thus the previously mentioned second Piola-Kirchhoff stress tensor:

Σ Σ Σ =    Σ 11 Σ 12 Σ 13 Σ 21 Σ 22 Σ 23 Σ 31 Σ 32 Σ 33    (3.4)
can be reduced to the second Piola-Kirchhoff stress vector:

⃗ s = s 1 s 2 s 3 s 4 s 5 s 6 T = Σ 11 Σ 22 Σ 33 Σ 23 Σ 13 Σ 12 T . (3.5)
In the same way, the Green-Lagrange strain tensor:

E =    E 11 E 12 E 13 E 21 E 22 E 23 E 31 E 32 E 33    (3.6)
can be reduced to the Green-Lagrange strain vector: Please note that ⃗ e r , ⃗ e s and ⃗ e 0 are the respective Voigt notations of the strain Green-Lagrange tensors E r , E s and E 0 , and are determined in the same manner as ⃗ e (Equation 3.7). The shear components of ⃗ e 0 are considered null and the Green-Lagrange strain at crack closure is thus defined with Voigt notations by: ⃗ e 0 = e 0 1 e 0 2 e 0 3 0 0 0

⃗ e =
T (3.8)
where e 0 i are parameters of the model. Hence and by using the newly introduced stress and strain vectors, the constitutive relation can now be defined as follows:

⃗ s = C • ⃗ ē -C 0 • (⃗ e r +⃗ e s -⃗ e 0 ) (3.9) with ⃗ ē =⃗ e -⃗ e 0 , (3.10)
C 0 the elastic second-order stiffness tensor and C the effective second-order stiffness tensor. C 0 and C are respectively defined in section 3.2.2 and section 3.2.3.

Elastic stiffness tensor

As a first step, the purely elastic behaviour of a fabric layer is described thanks to the elastic stiffness and compliance tensors, respectively C 0 and S 0 = (C 0 ) -1 . In view of the architecture of the preforms with two orthogonal directions of reinforcement, the layers have orthotropic behaviours.

As a result in the material coordinate system of each preform previously defined, the compliance tensor takes the form:

S 0 =                   1 E 1 - ν 21 E 2 - ν 31 E 3 0 0 0 - ν 12 E 1 1 E 2 - ν 32 E 3 0 0 0 - ν 13 E 1 - ν 23 E 2 1 E 3 0 0 0 0 0 0 1 G 23 0 0 0 0 0 0 1 G 13 0 0 0 0 0 0 1 G 12                   , (3.11)
which leads to the expression of the elastic stiffness tensor:

C 0 =               1 -ν 23 ν 32 E 2 E 3 ∆ ν 21 + ν 31 ν 23 E 2 E 3 ∆ ν 31 + ν 21 ν 32 E 2 E 3 ∆ 0 0 0 ν 12 + ν 13 ν 32 E 3 E 1 ∆ 1 -ν 31 ν 13 E 3 E 1 ∆ ν 32 + ν 31 ν 12 E 3 E 1 ∆ 0 0 0 ν 13 + ν 12 ν 23 E 1 E 2 ∆ ν 23 + ν 13 ν 21 E 1 E 2 ∆ 1 -ν 12 ν 21 E 1 E 2 ∆ 0 0 0 0 0 0 G 23 0 0 0 0 0 0 G 13 0 0 0 0 0 0 G 12               (3.12)
where

∆ = 1 -ν 12 ν 23 -ν 23 ν 32 -ν 31 ν 13 -2ν 12 ν 23 ν 31 E 1 E 2 E 3 . (3.13)
Among the twelve coefficients introduced in these last tensors, three groups can be distinguished:

• E i , the elastic modulus along the axis i,

• G i j , the shear modulus of a plane normal to the axis i in the direction j,

• and ν i j , Poisson's ratios reflecting the deformation in direction j under imposed displacement in direction i.

Whereas the elastic shear modulus G i j are considered constant, the elastic behaviour in the directions of reinforcement is non-linear. It is imputed to both yarn flattening and intrinsic non-linear elasticity of carbon fibres (see Section 2.1.1). The elastic modulus E 1 and E 2 are hence function of the strain by following the relations

E i = E 0 i + α i • e i , i ∈ {1, 2}. (3.14)
where E 0 i and α i are material parameters. Symmetry of stiffness and compliance tensors is used to reduce the number of parameters. The elastic behaviour can be fully described by nine independent elastic coefficients and three Poisson's ratios are arbitrarily taken to be dependent of the others. They are obtained from the relation:

ν i j = ν ji • E i E j
with i ∈ {2, 3}, j ∈ {1, 2} and i ̸ = j.

(3.15)

Effective stiffness tensor

The damage is introduced by adding additional compliance to the elastic compliance tensor S 0 , and as first step in the establishment of the model, only the matrix damage is considered.

The effective compliance tensor is thus defined as C = ( S) -1 with S = S 0 + ∆S m . The additional compliance tensor due to the matrix damage is given by:

∆S m = ∑ i η m i d m i H m i (3.16)
where H m i is the compliance tensor associated with the damage variable d m i . This model uses three damage variables, two corresponding to the damage along the directions of reinforcement while the last one is used to describe the out-of-plane damage. η m i represents the crack closure index which varies from 0 (closed crack) to 1 (opened crack). It is defined by the following relation:

η m i =        1 if ∆e i ≤ ēi 1 2 1 -cos π 2 ēi +∆e i ∆e i if -∆e i ≤ ēi ≤ ∆e i 0 if ēi ≤ -∆e i . (3.17)
∆e i depicts a strain tolerance between a state where all cracks are closed and a state where all cracks are open. This dispersion being closely tied to the number of cracks, ∆e i is set dependent to the crack ratio by the relation:

∆e i = (1 + a m i • d m i )∆e 0 i (3.18)
where a m i and ∆e 0 i are parameters of the model.

In order to describe separately the different fracture modes the compliance tensor H m i is split into both compliance tensors:

H m i = H nm i + H tm i . (3.19)
H nm i depicts the additional compliance due to normal loading to the crack (mode I) and H tm i the additional compliance due to tangent loading to the crack (mode II and III). They are defined by

H nm 1 =           h I m(1) S 0 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0           , (3.20)
H tm 1 =            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h III m(1) S 0 55 0 0 0 0 0 0 h II m(1) S 0 66            (3.21)
and H nm 2 , H nm 3 , H tm 2 and H tm 3 are obtained by index permutations. h I m(i) , h II m(i) and h III m(i) are parameters of the model.

Damage evolution

By following the same approach than the additional compliance tensor, two driving forces affect the evolution of each damage variable: the normal y nm i and the tangential y tm i to the damage directions and are defined by

y nm 1 = 1 2 e + 1 C 0 11 e + 1 y tm 1 = 1 2 e + 5 C 0 55 e + 5 + b 1 e + 6 C 0 66 e + 6 (3.22)
where y nm 2 , y nm 3 , y tm 2 and y tm 3 are obtained by index permutations and b 1 , b 2 and b 3 are material parameters. The thermodynamic forces so defined are dependent of the effective strain tensor E + , which corresponds to the positive part of the spectral decomposition:

E + = P •    ⟨E I ⟩ 0 0 0 ⟨E II ⟩ 0 0 0 ⟨E III ⟩    • P T (3.23)
where E I , E II and E III are the eigenvalues of E and P the transformation matrix formed by the eigenvectors. Lastly e + i are the components of ⃗ e + , namely the Voigt notation of E + . From these thermodynamic forces the damage criterion is defined as follows:

F m i = f nm i (y nm i ) + f tm i (y tm i ) -d m i ≤ 0 (3.24)
with f the cumulative distribution function of Weibull:

f x i = d x c(i) •      1 -exp      -     y x i -y x 0(i) y x c(i)     p x i           (3.25)
where d x c(i) , y x 0(i) , y x c(i) and p x i are material parameters with "x" taking the value "nm" or "tm".

Residual strain

Following the appearance of damage, the strains do not recover potentially their initial state when the stresses are relaxed. Some phenomena close to the cracks and which occur at a microscopic scale such as micro-plasticity or debris inside gaps may prevent a complete closure.

Following the description given by Marcin [START_REF] Marcin | Modélisation du comportement, de l'endommagement et de la rupture de matériaux composites à renforts tissés pour le dimensionnement robuste de structures[END_REF] for the ODM_MS, the residual strains are determined by following the relation:

⃗ e r = S 0 • ∑ i ζ i η m i ḋm i C • H m i • C • ⃗ ē (3.26)
where ζ i is a parameter of the model.

However, for further developments dedicated to extend the formulation of the stored strains (Section 3.2.7), the evolution of the residual strains has to be split according to each crack direction and into two components. The first one corresponds to the evolution due to opening following a mode I (normal loading) of the cracks normal to the direction i, mentioned as ⃗ e r,nm i , while the second one corresponds to the opening following a mode II or III (tangential loading) of the cracks normal to the direction i, mentioned as ⃗ e r,tm i . They are determined by the following relations:

⃗ e r,nm i = S 0 • ζ i η m i ḋm i C • H nm i • C • ⃗ ē, (3.27)
⃗ e r,tm i = S 0 • ζ i η m i ḋm i C • H tm i • C • ⃗ ē, (3.28)
and the total evolution of the residual strains are then recovered by:

⃗ e r = ∑ i ⃗ e r,nm i + ⃗ e r,tm i . (3.29)

Stored strain

By contrast with the residual strain formulation, the evolution of the stored strains do not depend on the evolution of the damage ( ḋm i ) but on the evolution of the crack closure index ( ηm i ). It is given by the relation:

⃗ e s = -S 0 • ∑ i ηm i d m i C • H m i • C • ⃗ ē. (3.30)
The stored strain was initially introduced by Chaboche and Maire [START_REF] Chaboche | A new micromechanics based CDM model and its application to CMC's[END_REF] to ensure the recovery of the initial elastic stiffness after the damage deactivation. Moreover, it avoids as well a discontinuity of the (⃗ s,⃗ e) response for complex loading cases. To illustrate the effect of the stored strains, an example of complex loading case, which consists of cyclic simple shear and traction/compression, is applied to a single shell element (Figure 3.2). By looking at the stress response of the model, it appears that the implementation of the stored strains ensures the continuity of the shear stress in case of opening or closure of the cracks (Figure 3.3).

Introduction of friction mechanisms

Physically, the stored strain ⃗ e s can be regarded as representative of the position of the crack lips after closure. For this above formulation, the hypothesis of infinite friction coefficient for closed crack (no evolution of the stored strain after the crack closure) makes the model unable to describe a realistic shear behaviour. To further improve the matrix damage model, an original and minimalist way to introduce friction effects at crack lips is presented in the rest of this section.

First, it is convenient to temporary rewrite the constitutive equation 3.9:

⃗ s =⃗ s e +⃗ s m (3.31)
where

⃗ s e = S 0 + ∑ i d m i H m i -1 • ⃗ ē (3.32)
is the stress applied on the healthy zone of a representative volume element, and

⃗ s m = C 0 •⃗ e 0 + ∑ i (1 -η m i ) d m i C • H m i • C 0 • ⃗ ē -C 0 • (⃗ e r i +⃗ e s i ) (3.33)
is the stress applied on an area which surrounds the closed cracks normal to the i-direction.

Note that the open cracks are not considered in this constitutive relation since an open crack is unable to carry loads. Then, to further apply a Coulomb friction law for the cracks, the normal and the tangential stresses, respectively ⃗ s nm i and ⃗ s tm i , which act on the closed cracks are isolated. They are then given by:

⃗ s nm i = (1 -η m i ) d m i C • H nm i • C 0 • ⃗ ē -C 0 • (⃗ e r,nm i +⃗ e s,nm i -⃗ e 0,nm i ) (3.34)
with e 0,nm i j = δ i j • e 0 j where δ i j is the Kronecker symbol, and

⃗ s tm i = (1 -η m i ) d m i C • H tm i • C 0 • ⃗ ē -C 0 • (⃗ e r,tm i +⃗ e s,tm i ). (3.35)
As a result, a criterion which states that the cracks are subject to static friction (≤ 0) or to dynamic friction (> 0) is introduced and is given by:

ss i = ⃗ s tm i -µ i ∥⃗ s nm i ∥ . (3.36)
where µ i is a coefficient of friction.

In case of a positive friction criterion, the crack lips start sliding to recover an equilibrium state. In the present work, the simplest modelling of the friction effect is considered, and is based on two hypotheses:

• the dynamic friction coefficient is equal to the static friction coefficient µ i , thus the equilibrium state is recovered when the friction criterion ss i is null again;

• and the displacement of the lips, to recover the equilibrium state, is instantaneous.

It leads to the following equation:

ss i + ∆ ss i = 0, (3.37)
with ∆ ss i the instantaneous evolution of the friction criterion to recover an equilibrium state. Since the normal stress applied to a crack does not evolve during sliding, the instantaneous evolution of the friction criterion is equal to the norm of the variation of the tangential stress applied to the crack:

∆ ss i = ∆⃗ s tm i . (3.38)
Thus it is now possible to determine ∆ ⃗ ⃗ s tm i , the vector notation of the variation of the tangential stress applied to the crack, by following the relation:

∆⃗ s tm i = -ss i • ⃗ s tm i ⃗ s tm i . (3.39)
Finally, as ∆⃗ s tm i is only dependent on the tangential displacement of the crack lips and by using the equation (3.35), the evolution of the tangential stored strains due to the friction mechanisms are determined by using the following equation:

∆⃗ e s,tm i = -S 0 • ∆⃗ s tm i . (3.40)
However from a microscopic aspect the friction coefficient is dependent on the roughness of the crack lips. And there is definitely no one unique condition of surface, which leads to introduce a dispersion on the friction. Instead of considering all closed cracks together, N families are introduced and the stress applied to the closed cracks is assumed to be homogeneously distributed on them. Thus, for a family j of cracks normal to the i-direction, the normal and the tangential stress, respectively ⃗ s nm i j and ⃗ s tm i j , are given by:

⃗ s nm i j = (1 -η m i ) d m i C • H nm i • C 0 • ⃗ ē -C 0 • (⃗ e r,nm i +⃗ e s,nm i j -⃗ e 0,nm i ) (3.41)
and

⃗ s tm i j = (1 -η m i ) d m i C • H tm i • C 0 • ⃗ ē -C 0 • (⃗ e r,tm i +⃗ e s,tm i j ). (3.42)
Then, the evolution of the tangential stored strains ∆⃗ e s,tm i j is determined by following the procedure given by the equations 3.36 to 3.40. The only difference is the use of a coefficient of friction µ i j proper to each group. Finally the total stored strain evolution due to friction mechanisms is obtained by averaging every movements of crack lips:

∆⃗ e s,tm i = 1 N N ∑ j=1 ∆⃗ e s,tm i j . (3.43)
About the range of the coefficients of friction, a quadratic and homogeneous distribution on a given interval is carried out. Thus they are computed by following the relation:

µ i j = µ max i N 2 • j 2 with j ∈ [0, N -1]. (3.44)
However, as for the crack closure index, the distribution of the friction coefficients is sensitive to the damage. As a result the upper bound of their range is given by:

µ max i = (1 + m m i • d m i ) µ 0 i (3.45)
where m m i and µ 0 i are parameters of the model.

Particularities of the textile preforms

In the particular case of plies made up with textile preforms -where the yarns are interlacedthe shear locking phenomenon occurs when the adjacent yarns come into contact. It results a significant and progressive rise in the in-plane shear stiffness. In this model, the shear locking is introduced by using the following relation:

G 12 = G 0 12 + g l 12 • e 6 -e l 6 (3.46)
where G 0 12 , g l 12 and e l 6 are parameters of the model. Note that the operator ⟨•⟩ is called the Macauley bracket and ⟨x⟩ = x if x > 0, and null in others cases. Also, it has been proved that the large shearing ability of the textile plies leads to massive out-of plane damages even for purely in-plane loading (Figure 3.4). Its emergence is due to the opposite directions of rotation of transverse yarns. As a result, the matrix between both transverse yarns in the crossing area is subjected to torsion and causes breaking.

Marcin [START_REF] Marcin | Modélisation du comportement, de l'endommagement et de la rupture de matériaux composites à renforts tissés pour le dimensionnement robuste de structures[END_REF] introduced an in-plane/out-of-plane coupling by using a new thermodynamic force characterising these torsional modes. Consequently, the out-of plane damage is controlled by the previously defined normal and tangential forces, and the newly introducing coupling force:

     y nm 3 = 1 2 e + 3 C 0 33 e + 3 y tm 3 = 1 2 e + 4 C 0 44 e + 4 + b 3 e + 5 C 0 55 e + 5 y p 3 = 1 2 e + 6 C 0 66 e + 6 . (3.47)
The associated damage criterion is then actualised and given by:

F m 3 = f nm 3 (y nm 3 ) + f tm 3 (y nm 3 ) + f p 3 (y p 3 ) -d m 3 ≤ 0. (3.48)
But in addition to Marcin's formulation, the in-plane shear stiffness is considered sensitive to the out-of-plane damage since the rotation between both transverse yarns is made easier. Thus, the additional compliance tensor due to tangential loading is modified to take into account this influence and is now given by:

H tm 3 =            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h II m(3) S 0 44 0 0 0 0 0 0 h III m(3) S 0 55 0 0 0 0 0 0 h p m(3) S 0 66            . (3.49)
Therefore, the matrix damage model is able to consider a loss of the elastic stiffness and to take into account an eventual local plasticity around the cracks, the friction mechanisms and the shear locking of the textile preforms. Regarding the friction, the coefficients are evolving with the damage value.

Verification of the thermodynamic consistency

As a result of the addition of the friction mechanisms inside the material model, the consistency with regards to the thermodynamic laws needs to be checked.

As the analytical audit of the second law of thermodynamics is very complex, only basic audits based on extreme cases have be done.

The first one consists of a comparison of the stress response between the newly formulated material model with an infinite friction coefficient and the Onera Damage Model MicroStructure (ODM_MS). Because of the hypothesis of infinite friction for the formulation of the ODM_MS, the stress responses of both models have to be equivalent. The displacement which is imposed to a shell element is the same as the one which was presented Figure 3.2. The results in case of infinite friction are compared Figure 3.5 and are entirely similar. Rather, the second check consists of a comparison of the stress response between the newly formulated material model with an infinite friction coefficient and with an absence of friction at crack lips. The continuity of the shear stress response is preserved. The only difference concerns the shear stress evolution when the cracks are closed. The transmitted shear now corresponds to the stress which is transmitted by an healthy area of the material. As a consequence, the shear stress is smaller than in the case of infinite friction (Figure 3 The last check concerns the dissipated energy during crack lips displacements. This dissipated energy is represented graphically by hysteresis loop. Due to the clockwise orientation of these loops during friction mechanisms, and due to the symmetric behaviour of the friction mechanisms for positive or negative shear loading (which leads to recover the same stress state after a unloading/loading cycle), the dissipated energy is always positive. The second law of thermodynamics is thus ensured.

Implementation

The aim of the implementation of a material model in a finite element software is to provide, at the current time step, the stress at an integration point according to the strain value given by the element formulation. It is carried out by the implementation of the developed constitutive model in an incremental loop.

In an explicit scheme, the computation time is conditioned by the efficiency of the material model. It is thus essential to limit as much as possible the internal loops. By using a strain-space formulation, the explicit scheme (stress obtained from the strain without iteration) is maintained. Thanks to the strain field provided by the finite element software in the material coordinate system, the thermodynamic forces are determined, leading to the damage evaluation. The crack closure indexes are evaluated and then the residual and the stored strains are computed by means of an incremental scheme. Note that the total Lagrangian framework helps to maintain the model efficiency by avoiding the use of complex algorithm to ensure the objectivity. Once the damage values, the crack closure indexes, the residual and the stored strains are determined, it is possible to calculate the stress thanks to the constitutive relation. The details are given in Algorithm 3.1.

The model is implemented for both solid or shell elements. For the shell elements, it is possible to use this model under plane stress assumption. It uses the three-dimensional formulation described in this work and a stress return algorithm to evaluate the thickness reduction. Another possibility is to use a two-dimensional formulation (without consideration of the out-of-plane strain and stress) but with the possibility to take into account the transverse shear behaviour of the material.

Identification

The parameter identification of the matrix damage model is done by means of two different experimental tests carried out independently on the different preforms. These tests are closed to the standardised monotonic tensile and in-plane shear tests, respectively the NF EN Algorithm 3.1 Stress update algorithm with the matrix damage evaluation.

Step 1 Calculation of the effective strain ⃗ e + (Equation 3.23).

Step 2 Calculation of the thermodynamic forces y x i (Equations 3.22 and 3.47).

Step 3 Calculation of the damage d m i (Equation 3.25).

Step 4 Calculation of the crack closure index η m i (Equation 3.17).

Step 5 Calculation of the effective stiffness tensor C (Equation 3.16).

Step 6 Calculation of the residual strain ⃗ e r by following the iterative relation:

⃗ e r t+∆t =⃗ e r t + S 0 • ∑ i ζ i η m i ∆d m i C • H m i • C • 1 2 ⃗ ēt+∆t + ⃗ ēt .
Step 7 Calculation of the stored strain ⃗ e s due to crack closure by following the iterative relation:

⃗ e s t+∆t =⃗ e s t -S 0 • ∑ i ∆η m i d m i C • H m i • C • 1 2 ⃗ ēt+∆t + ⃗ ēt .
Step 8 Calculation of the stored strain evolution due to friction mechanisms ∆⃗ e s,tm i (Equations in Section 3.2.7) and addition to the previously computed stored strain.

Step 9 Calculation of the stress ⃗ s (Equation 3.9). ISO 527-4 and NF EN ISO 14129 tests, in order to be easily reproducible in an industrial framework. The approach consists in eliminating gradually and systematically the unknown parameters.

These tests are carried out at room temperature on an electromagnetic device (Sintech 20D) which ensures 100 kN for the maximal load capacity and at a speed of 2 mm min -1 . The specimens are cut using the water-jet technique and their shape and dimensions, identical for both tests, are shown in Figure 3

Longitudinal and transverse tensile tests

As a first step, the parameters describing the longitudinal and the transverse behaviours are determined through tensile tests in the fibre directions (Figure 3.8). Some materials, like the ones taken as example, show an elastic brittle behaviour for longitudinal and transverse loading. In this case, the standardised NF EN ISO 527-4 test is used as such. Otherwise, in case of matrix cracks appearing before the final fracture of the coupon (due to the fibre failure), the tensile loading is increased on a cyclical basis. This last method allows to quantify the damage evolution and its effect on the mechanical behaviour.

Due to small strains at failure and none apparent damage, two gauges have been used to measure:

• ε xx = ε( ⃗ X) the elongation in the centre of the coupon according to the ⃗ X direction,

• ε yy = ε( ⃗ Y ) the elongation in the centre of the coupon according to the ⃗ Y direction.

The force applied to the coupon F = ⃗ F • ⃗ X is measured with a 100 kN load cell and then the engineering stress is obtained by: where S 0 is the initial cross-section of the coupon.

σ xx = F S 0 (3.50) ⃗ X = ⃗ 1 ⃗ Y = ⃗ 2 (a) Longitudinal tensile test. ⃗ X ⃗ Y = ⃗ 1 ⃗ 2 90°( b) Transverse tensile test.
The parameters E 0 1 , α 1 , E 0 2 and α 2 are then obtained through linear regressions of the curves σ xx ε xx = f (ε xx ) from the longitudinal and transverse tensile tests, whereas ν 0 12 is obtained through the curve -

ε yy ε xx = f (ε xx ) from the longitudinal tensile test.

In-plane shear tests

Once the properties in the longitudinal and transverse directions are determined, the in-plane shear properties are obtained from a close version of the standardised NF EN ISO 14129 in-plane shear test. It consists of a tensile test according to an angle of 45°against the fibre directions (Figure 3.9).

In this configuration, the matrix damage is important and leads to substantial non-linearity in the mechanical behaviour. Therefore, to quantify the level of damage and its effects on the behaviour, the tensile loading is increased on a cyclical basis.

⃗ X ⃗ Y ⃗ 1 ⃗ 2 45°F
igure 3.9: Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for in-plane shear tests.

Due to the large shear strain ability of the fabric preforms, extensometers are preferred to measure:

• ∆L xx the relative displacement of both clips, initially set L 0 xx apart, of the longitudinal extensometer,

• ∆L yy the relative displacement of both clips, initially set L 0 yy apart, of the transverse extensometer.

The in-plane shear angle is then given by:

γ 12 = ∆L xx L 0 xx - ∆L yy L 0 yy . (3.51)
The force applied to the coupon F = ⃗ F • ⃗ X is measured with a 40 kN load cell and then the engineering in-plane shear stress is obtained by:

σ 12 = F 2 × S 0 (3.52)
where S 0 is the initial cross-section of the coupon. The initial slope of the curve σ 12 = f (γ 12 ) defines the initial in-plane shear stiffness G 0 12 . Then the identification of the parameters associated to the in-plane matrix damage is much more difficult due to the various coefficients which define the cumulative distribution functions of Weibull (Equation 3.25), but also due to the interdependence between the damage and the friction mechanisms. Consequently, these parameters are obtained through an optimisation study to minimise the Mean Square Error (MSE) between the experimental and the numerical responses F = f (∆L xx ).

In the particular case of FRP made up with textile preforms two optimisation studies are needed due to the important out-of plane damage after the shear locking. The shear locking angle and the resulting Green-Lagrange shear strain at shear locking can be determined through the direct observation of the force-displacement curve. At shear locking, the material becomes significantly more rigid and a slight re-hardening appears. Once the shear angle determined, the same optimisation procedure as explained above is used on the responses before the shear locking in order to describe the in-plane damage evolution (d m 1 and d m 2 ). Then an additional optimisation study is done on the out-of-plane damage evolution (d m 3 ) and the parameter g l 12 describing the non-linearity of the in-plane shear stiffness due to shear locking.

The out-of-plane damage induced by matrix torsion between both transverse yarns is considered to occur after the shear locking. It can be understood by the study of the unit-cell deformation of a textile preform. By means of Digital Image Correlation carried out during in-plane shear tests, the rotation of the yarns have been observed (Figure 3.10). The yarns start to rotate in the inter-yarn area. As soon as adjacent yarns come into contact (shear locking), the yarns start to rotate between them. This second mechanism leads to the torsion of the matrix between both transverse yarns, and consequently to the out-of-plane damage. By proceeding in this manner, step by step on simple experimental tests, the global singularity of the solution is preserved and the parameters maintain a physical relevance. Moreover in case of balanced fabrics, the number of parameters is considerably lowered.

The results of the parameter identification with ten families of friction mechanisms are shown in Table 3.1.

Validation of the identification procedure

The present matrix damage model is implemented in FORTRAN 90 in a user material subroutine as described section 3.4 for the explicit finite element code LS-DYNA ® . In this validation procedure, shell elements are used with the two-dimensional formulation of the material model. As a compromise between accuracy and efficiency, ten families of friction criterion have been used in this study.

This section presents the results of the identification procedure which have been presented Section 3.5 on three various fabric preforms, namely a bi-axial non-crimp fabric (see Figure 3.1a), a 3K plain-weave woven (see Figure 3.1b) and a 12K plain-weave woven. The optimisation procedures were carried out thanks to the commercial software LS-OPT ® to 2100 0 e l 6 0 0.17

0 b 1 = b 2 0 0 0 d nm c(1) = d nm c(2) 0 0 0 y nm 0(1) = y nm 0(2) 0 0 0 y nm c(1) = y nm c(2) 0 0 0 p nm 1 = p nm 2 0 0 0 d tm c(1) = d tm c(2)
0.6324 0.6745 0.6992

y tm 0(1) = y tm 0(2)
0.0041 0.0041 0.0041

y tm c(1) = y tm c(2)
2.124 1.729 0.04371 As a first step of the validation, a longitudinal tensile test of both preforms is simulated and the results are compared with the experimental tests. These results are shown Figures 3.11, 3.12 and 3.13. Because of limited non-linearity for loading in the fibre directions, the mechanical behaviour is well described. The in-plane shear behaviour, observed thanks to tensile tests with an angle of 45°a ccording to the longitudinal fibre direction, is quite a bit more challenging to model. Although various fabric plies were tested, the matrix damage model is able to well described their large non-linear behaviour as shown Figures 3.14 With few additional parameters compared to the previous formulations of the Onera Damage Model, because of the dependence between the friction mechanisms and the matrix crack density, the dissipated energy due to the hysteresis cycles can be assessed. Moreover, in case of positive/negative cyclic shearing this formulation avoids the previous limitations due to the hypothesis of infinite friction or null friction for the lips of the matrix cracks.

p tm 1 = p tm 2 0.2248 0.097 0.281 d tm c(3) / 0 / y tm 0(3) / 0 / y tm c(3) / 0 / p tm 3 / 0 / h I m(1) = h I m(2) 0 0 0 h II m(1) = h II m(2) 1 1 1.0467 h III m(1) = h III m(2) 0 0 0 d p c(3) / 0.3253 / y p 0(3) / 17.9 / y p c(3) / 20 / p p 3 / 1.751 / h p m(3) / 1 / ζ 1 = ζ 2 0 0 0.011 ζ 3 / 0 / m 0 1 = m 0 2 0.00838 0.0498 0.00558 m m 1 = m m
Finally, the use of the finite strain framework as well as the modelling of the shear locking phenomenon let to maintain the objectivity and the consideration of the real material orientation.

Concluding remarks on the matrix damage

The matrix damage model for fabric reinforced polymers is a close version of the Onera Damage MicroStructure model. A new evolution law for the stored strains has been presented. The friction mechanisms induced by the matrix cracks are now taken into account in a simple manner. It allows a good description of the shear behaviour, including the approximation of the dissipated energy due to the subsequent hysteresis loops.

Because of the large rotation of the yarns for shear loading, the model is extended in finite strain. The total Lagrangian formulation is used in order to well track the fibre orientation and ensure the objectivity. The shear locking is also introduced for the textile composites.

The procedure for the parameter identification is also provided. It consists of eliminating gradually and systematically the unknown parameters. It is done either by direct measurement on the experimental results, or by optimisation of restrained sets of parameters in order to keep the singularity and the physical meaning of the results.

From the numerical point of view, the model is implemented in the commercial finite element software LS-DYNA ® . It is validated through standardised tensile and in-plane shear experimental tests. The simulation results show the good efficiency of the proposed model for fully different fabric preforms, such as non-crimp or woven.

However, in order to fully simulate the behaviour of layered fabric composites, additional physical phenomena have to be taken into account. Such is the case, for instance, viscoelas-ticity, fibre failure or intralaminar damage. The newly introduced friction mechanisms are also as simple as possible (based on Coulomb criterion) and could be further improved.

Introduction

To introduce a strain-rate sensitivity for dynamic loading, phenomenological models exist and they describe empirically the dependence of the elastic modulus (and possibly also the damage evolution, the failure criterion, etc.) on the strain rate by polynomial or logarithmic functions [START_REF] Daniel | Strain-rate-dependent failure criteria for composites[END_REF][START_REF] Raimondo | Modelling of strain rate effects on matrix dominated elastic and failure properties of unidirectional fibrereinforced polymer-matrix composites[END_REF]. However, these models may suffer of instabilities during finite element analysis due to the difficulty to obtain a realistic instantaneous strain rate due to numerical instabilities.

In the wide family of the viscoelastic models used to described the strain rate sensitivity of the elastic behaviour, the unidirectional rheological ones are the most simple ones. They are based on the combination of two basic components, a purely elastic spring (Hooke element) and a purely viscous damper (Newton element), connected in parallel (Kelvin-Voigt model) or in series (Maxwell model). But these models cannot describe simultaneously creep and stress relaxation. In order to ensure this, an additional spring is added in series of the Kelvin-Voigt model or in parallel of the Maxwell model. They are respectively called the Boltzmann (Figure 4.1a) and the Zener ( model since the stress applied to each sub-element is equal to the total stress applied to the material. Instead, the generalised Maxwell model is preferred for strain-space formulations.

Their relaxation functions are given by the use of Prony series given by:

f (t) = N ∑ i=1 A i • e - t τ i . (4.1)
These models are widely used in the commercial finite element analysis software but suffer two drawbacks. The first one is the use of high-order Prony series N to approach the viscoelastic behaviour and, therefore, a significant number of parameters to identify. Second, as defined these generalised rheological viscoelastic models are linear and cannot represent the non-linear viscoelastic behaviour of the fibre reinforced polymers.

Another group is the family of the spectral models. Compared to the previously mentioned generalised rheological models, the spectral models provide a continuous spectrum of the relaxation times. The number of parameters to identify is thus reduced without losing accuracy for describing the viscous phenomena. By using a Gaussian description of the spectrum, Maire [START_REF] Maire | Études théorique et expérimentale du comportement de matériaux composites en contraintes planes[END_REF] introduces a spectral model for the fibre reinforced polymers. Thereafter, Rémy-Petipas [START_REF] Rémy-Petipas | Analyse et prévision du comportement à long terme des composites fibres de carbone/matrice organique[END_REF], Schieffer et al. [START_REF] Schieffer | A coupled analysis of mechanical behaviour and ageing for polymer-matrix composites[END_REF], Huchette [START_REF] Huchette | Sur la complémentarité des approches expérimentales et numériques pour la modélisation des mécanismes d'endommagement des composites stratifiés[END_REF], Berthe et al. [START_REF] Berthe | Thermo-viscoelastic modelling of organic matrix composite behaviour -Application to T700gc/M21[END_REF] have gradually improved the model by introducing thermal effects and damage coupling. The spectral models expressed as such, however, are suffering substantial computational times in case of explicit simulation scheme.

Functional formulations are also used to describe time irreversibility problems. These formulations rely on the basis that the instantaneous response of a material depends on the loading history. Thus, the Boltzmann superposition principle can be applied to viscoelastic stress analysis problems. Initially using a linear creep compliance in the formulation, Lou and Schapery [START_REF] Lou | Viscoelastic behavior of a nonlinear fiber-reinforced plastic[END_REF] introduced a non-linear viscoelastic formulation for the fibre reinforced polymers. This model was successfully used for various unidirectional composites [START_REF] Heil | The nonlinear viscoelastic response of resin matrix composite laminates[END_REF][START_REF] Tuttle | Prediction of the long-term creep compliance of general composite laminates[END_REF][START_REF] Zaoutsos | On the non-linear viscoelastic behaviour of polymer-matrix composites[END_REF] but essentially for creep simulations.

Balieu [START_REF] Balieu | Modèle viscoélastique-viscoplastique couplé avec endommagement pour les matériaux polymères semi-cristallins[END_REF] uses a functional formulation of the generalised Maxwell model introduced by Simo and Hughes [START_REF] Simo | Viscoelasticity[END_REF] to model the behaviour of semi-crystalline polymers at high strain-rates. He introduces also a non-linear viscoelastic behaviour by coupling the damage and the viscoelasticity. This formulation is very attractive by its computational efficiency, notably in explicit finite element simulations, and its simple implementation. But by using a generalised Maxwell model as basis it always remains to identify a consequent number of parameters. Despite this, this last model is chosen as viscoelastic model in the present work by its attractiveness for impact simulations. Moreover, a coupling with the matrix damage model is possible to introduce the non-linear viscoelastic behaviour of the fibre-reinforced polymer.

In a first section, the linear generalised Maxwell model is described, with the application of the Boltzmann superposition principle to obtain the constitutive equation. Then, a coupling with the previously established matrix damage model is proposed. Moreover, in a third section, the formulation of this model in the finite strain framework is discussed. The implementation of this newly formulated model into an explicit finite element code is given afterwards. The next section presents the parameter identification through straightforward Dynamical Mechanical Analysis. Finally, the model is validated by an experimental test campaign carried out on a high-speed hydraulic jack facility. The details of these tests are provided in a last section, with a particular attention to the scale effect, critical in the damage analysis of the fabric reinforced polymers due to substantial side effects.

Linear generalised Maxwell model

The generalised Maxwell model relies on the combination of only two basic elements: a spring called a Hooke element, and a dumper called a Newton element. By arranging these elements in series and in parallel according to a scheme given by the generalised Maxwell model, the model is able to describe the increase of the stiffness of polymers at increasing strain rates.

In this section the small strain formalism is adopted. The formulation of the viscoelastic model in the total Lagrangian framework will be discussed Section 4.4.

Linear viscoelasticity theory

The basic Hooke and Newton elements are used to describe the linear constitutive relations between the stress, the strain, and the strain-rate.

The Hooke (elastic) element, schematised by a spring (Figure 4.3a), represents the linear relation between the elastic stress σ e and the elastic strain ε e :

σ e = E • ε e (4.2)
with E the elastic constant. Rather, the Newton (viscous) element, schematised by a dumper (Figure 4.3b), represents the linear relation between the viscous stress σ v and the viscous strain rate εv :

σ v = λ • εv (4.3)
with λ the coefficient of viscosity. Thus, the strain ε of a Maxwell element, schematised by a Hooke element and a Newton element in series (Figure 4.3c), is the sum of the strain of each element:

ε = ε e + ε v , (4.4)
while the stress on each element is the same as the imposed stress:

σ = σ e = σ v . (4.5)
By differentiating the total strain equation (Equation 4.4) of the Maxwell element, by using

σ e E ε e σ e
(a) Hooke element the constitutive relations given by the Hooke and Newton elements (Equations 4.2 and 4.3) and because the stress applied to each element is equal to the total stress (Equation 4.5), the total strain rate is thus given by:

σ v λ ε v σ v (b) Newton element σ E ε e λ ε v σ (c) Maxwell element
ε = εe + εv = σ E + σ λ . (4.6)
To solve this differential equation, a simple loading case is considered in a first time. It will provide the basic concepts for a more general resolution. In this way, a constant strain ε 0 introduced by a increment at t = 0 is imposed on the Maxwell element (Figure 4.4). In

ε t 0 ε 0 σ t 0 σ 0 σ (t)
σ + σ τ = 0, (4.7)
whose solution after time integration is:

σ (t) = σ 0 • exp - t τ = E(t) • ε 0 . (4.8)
E(t) is called the relaxation modulus and is given by:

E(t) = E • exp - t τ . (4.9)
For an infinite relaxation time, the stress of a standard Maxwell element relax completely to zero. This behaviour is not representative of a solid behaviour but corresponds to a fluid behaviour. To sustain a substantial stress after a large relaxation time, a Hooke element is introduced in parallel to the previous Maxwell model (

ε = ε e + ε v , (4.10)
σ = σ e + σ ∞ = σ v + σ ∞ . (4.11)
After resolution of the differential equation, the stress at time t is given by:

σ (t) = σ ∞ + σ 0 • exp - t τ = E(t) • ε 0 , (4.12)
which leads to the subsequent relaxation modulus:

E(t) = E ∞ + E • exp - t τ . (4.13)
With the present viscoelastic model, the stress does not relax completely to zero for a strain increment, but reaches an asymptotic value σ ∞ (Figure 4.6) related to the so-called long-term modulus E ∞ by:

σ ∞ = E ∞ • ε 0 . (4.14)
It can be also noted that the model is relaxing through a single relaxation time. However it is not sufficient to describe the dynamic mechanical spectrum of polymers, and consequently the fibre reinforced polymers. To approach the real spectrum, the number of relaxation time is increased by using N Maxwell elements in parallel of the Hooke element (Figure 4.7). This model is called the generalised Maxwell model. The strain applied to the generalised

ε t 0 ε 0 σ t 0 σ 0 σ (t) σ ∞
σ . . . E N ε e N λ N ε v N E 1 ε e 1 λ 1 ε v 1 E ∞ ε . . . σ Figure 4.7: Generalised Maxwell model.
Maxwell model is equal to the strain applied to each branch:

ε = ε j = ε e j + ε v j ∀ j ∈ [[1, N]], (4.15)
and the total stress is the sum of the stress applied to each branch:

σ = σ ∞ + N ∑ j=1 σ j , (4.16)
where the subscript j indicates the properties of the j th Maxwell element. The stress at time t is then determined through the superposition theorem and is given by:

σ (t) = σ ∞ + N ∑ j=1 σ 0 • exp - t τ j , (4.17)
with τ j the relaxation time of the j th Maxwell element and leads to the relaxation modulus which is defined by:

E(t) = E ∞ + N ∑ j=1 E j • exp - t τ j , (4.18)
with E j the elastic modulus of the j th Maxwell element. This expression of the relaxation modulus follows the form of Prony series.

Functional formulation of the generalised Maxwell model

To simplify the resolution of the differential equations, a single strain increment was considered as loading of the viscoelastic model and the relaxation stress response was described by introducing a relaxation modulus (Equation 4.18). But the response of the viscoelastic model has to be determined for random loading cases. The Boltzmann superposition principle suggests that the response of a material to a strain increment is independent of responses due to strain increments which have been previously initiated. Thus, let σ k (t) the stress at time t due to a strain increment ∆ε k applied at a time ξ k previous to t. By considering for example two strain increments, the total stress at time t can be obtained by superposition as follows:

σ (t) = σ 1 (t) + σ 2 (t) = E(t -ξ 1 ) • ∆ε 1 + E(t -ξ 2 ) • ∆ε 2 . (4.19)
In a more general case, the total stress at time t is obtained by summing the effects of an infinite number of perturbations and is given by:

σ (t) = t -∞ E(t -ξ ) • dε(ξ ) = t -∞ E(t -ξ ) • dε(ξ ) dξ dξ (4.20)
By replacing the relaxation modulus by the form of Prony series (Equation 4.18) obtained for a generalised Maxwell model, the total stress becomes:

σ (t) = t -∞ E(t -ξ ) • dε(ξ ) dξ dξ = t -∞ E ∞ • dε(ξ ) dξ dξ + t -∞ N ∑ j=1 E j • exp - t -ξ τ j • dε(ξ ) dξ dξ = E ∞ • ε(t) + N ∑ j=1 t -∞ E j • exp - t -ξ τ j • dε(ξ ) dξ dξ = E ∞ • ε(t) + N ∑ j=1 l j (t) (4.21)
with

l j (t) = t -∞ E j • exp - t -ξ τ j • dε(ξ ) dξ dξ . (4.22)
l j (t) represents the stress at time t on the j th Maxwell element of the generalised Maxwell model. It is coherent with the previous constitutive equation of the Maxwell element as the stress vanishes for an infinite time.

Viscoelasticity coupled with the matrix damage

The Boltzmann integral is able to be applied to problems in three-dimensions. Consequently, by analogy with the previous uni-axial approach, the scalar values can be replaced by their tensor forms. Hence the total stress at time t results in:

⃗ σ = C ∞ :⃗ ε(t) + N ∑ j=1 l j (t) (4.23)
where

l j (t) = t -∞ C j • exp - t -ξ τ j : d⃗ ε(ξ ) dξ dξ ∀ j ∈ [[1, N]]. (4.24)
C ∞ and C i are respectively the long-term stiffness tensor and the i th viscoelastic stiffness tensor.

The total Helmholtz free potential of the viscoelastic model is therefore postulated as the sum of a long term component and the viscous components as follows:

ψ = 1 2ρ ψ ∞ + N ∑ j=1 ψ v j (4.25) with    ψ ∞ = ⃗ ε : C ∞ :⃗ ε ψ v j = t -∞ t -∞ d⃗ ε(χ) dχ dχ : C j • exp - 2t -ξ -χ τ j : d⃗ ε(ξ ) dξ dξ . (4.26)
The non-linear behaviour of the fibre reinforced polymers is now introduced through the coupling of the present viscoelastic model with the previous matrix damage model. In the first instance, the damage is considered to be acting only on the long term response of the composite material. Hence, the long term Helmhotz free potential is replaced by the Helmhotz free energy of the matrix damage model, whose the expression is reminded below:

ψ ∞ = ψ m + ψ 0 -ψ r -ψ s (4.27)
with the various Helmholtz free energy representative of the effective stiffness energy (ψ m ), the initial internal energy (ψ 0 ), the energy dissipated by local plasticity mechanisms (ψ r ) and the energy dissipated by friction mechanisms (ψ s ) during the damage emergence.

Formulated such as, the energy due to the viscous part of the model is absolutely predominant after damage. However, the strain rate dependency of the fibre reinforced polymers is due to the viscoelastic property of the resin used as matrix. It seems obvious and it has been proved that the matrix damage affects the viscoelastic response and that the viscoelasticity in a fibre reinforced polymer is non-linear. Schieffer et al. [START_REF] Schieffer | A coupled analysis of mechanical behaviour and ageing for polymer-matrix composites[END_REF] or Huchette [START_REF] Huchette | Sur la complémentarité des approches expérimentales et numériques pour la modélisation des mécanismes d'endommagement des composites stratifiés[END_REF] uses an effective viscoelastic stiffness tensor directly dependent of the damage variables defined by the matrix damage equations. However, this damage effect on the viscoelasticity is well-adapted in case of open cracks, but for the adopted formalism relying on closed matrix cracks and friction effect, this damage coupling leads to completely inaccurate behaviours. Because of that, the viscoelastic stresses are assumed to only be issued of the healthy areas of the material. Around the cracks the stresses are not transmitted because of energy dissipation due to crack lips movements. Thus, it is convenient to define the strain ε h really applied to

EXTENSION OF THE CONSTITUTIVE VISCOELASTIC MODEL IN FINITE STRAIN 67

the healthy material. By the following constitutive equation:

⃗ σ = C : (⃗ ε -⃗ ε 0 ) -C 0 : (⃗ ε r +⃗ ε s -⃗ ε 0 ) = C : (⃗ ε h -⃗ ε 0 ) , (4.28)
the healthy strain is thus given by:

⃗ ε h =⃗ ε - C-1 C 0 : (⃗ ε r +⃗ ε s -⃗ ε 0 ). (4.29)
The viscous Helmhotz free energy is now defined by the following relation:

ψ v j = t -∞ t -∞ d⃗ ε h (χ) dχ dχ : C j • exp - 2t -ξ -χ τ j : d⃗ ε h (ξ ) dξ dξ . (4.30)
Then, the constitutive relation is obtained by derivation of this thermodynamic potential and is given by:

⃗ σ = ρ ∂ ψ ∂⃗ ε = C∞ : (⃗ ε -⃗ ε 0 ) -C 0 ∞ : (⃗ ε r +⃗ ε s -⃗ ε 0 ) + N ∑ j=1 l m j (t) (4.31)
where

l m j (t) = t -∞ C j • exp - t -ξ τ j : d⃗ ε h (ξ ) dξ dξ ∀ j ∈ [[1, N]] (4.32)
and represents the viscous stress due to the strains in an healthy area of the fabric reinforced polymers. This is a pragmatic modelling of the strain-rate dependency of the fibre reinforced polymers but efficient for explicit finite element analysis. Moreover, the coupling with the matrix damage model is done without additional parameter.

Extension of the constitutive viscoelastic model in finite strain

The formulation of a viscoelastic model in a total Lagrangian framework have been studied extensively [START_REF] Tallec | Three-dimensional incompressible viscoelasticity in large strains: Formulation and numerical approximation[END_REF][START_REF] Simo | Viscoelasticity[END_REF]. In previous work, Flory [START_REF] Flory | Thermodynamic relations for high elastic materials[END_REF] introduced the volumetric/deviatoric multiplicative split. It relies on the hypothesis that the viscoelastic volumetric and deviatoric responses are fully uncoupled. Kaliske [START_REF] Kaliske | A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains[END_REF] proposed an efficient formulation which separates the relaxation tensor, independent of the deformation, and the nonlinear elastic material tensor, expressed according to an hyperelastic formulation. It makes easier the parameter identification by fully uncoupling the elastic and the viscoelastic terms.

In this work, the damageable viscoelastic model is extended in finite strain by replacing the engineering strain and stress tensor by respectively the Green-Lagrange strain tensor and the second Piola-Kirchhoff tensor. This assumption was made for two reasons. Firstly, the maximal shear strain is less than 20% for the studied materials. Thus, the spurious softening which appears at constant stiffness modulus for compressive loading does not occur. Moreover, due to the small displacements in the fibre directions, the material is nearby incompressible during the problematic large shearing. The decomposition of the behaviour in volumetric and deviatoric components is therefore not essential.

Finally the constitutive relation in the finite strain framework and by using the Voigt notation is given by:

⃗ s =⃗ s ∞ (t) + N ∑ j=1 ⃗ l m j (t). (4.33)
with

⃗ l m j (t) = t -∞ C j • exp - t -ξ τ j • d⃗ e h (ξ ) dξ dξ (4.34)
and ⃗ e h =⃗ e -C-1 C 0 : (⃗ e r +⃗ e s -⃗ e 0 ). (4.35)

Implementation

The numerical implementation of the viscoelastic model coupled with the matrix damage model into a finite element code is done by the determination of the stress at time t n+1 with the knowledge of the variables and of the strain at time t n and the time increment

∆t = t n+1 -t n : ⃗ s =⃗ s ∞ (t n+1 ) + N ∑ j=1 ⃗ l m j (t n+1 ). (4.36)
The long-term stress ⃗ s ∞ (t n+1 ) is directly obtained by the explicit formulation of the matrix damage model. It therefore remains to determine the different viscous stress ⃗ l m j (t n+1 ). Using the previously established Boltzmann heredity integral and the additive properties of the integrals, the viscous stress at time t n+1 becomes:

⃗ l m j (t n+1 ) = t n+1 0 C j • exp - t n+1 -ξ τ j • d⃗ e h (ξ ) dξ dξ = t n 0 C j • exp - t n + ∆t -ξ τ j • d⃗ e h (ξ ) dξ dξ + t n+1 t n C j • exp - t n+1 -ξ τ j • d⃗ e h (ξ ) dξ dξ . (4.37)
Using the multiplicative of the exponential, the viscous stress at time t n+1 is expressed in function of the viscous stress at the precedent increment:

⃗ l m j (t n+1 ) = exp - ∆t τ j t n 0 C j • exp - t n -ξ τ j • d⃗ e h (ξ ) dξ dξ + t n+1 t n C j • exp - t n+1 -ξ τ j • d⃗ e h (ξ ) dξ dξ = exp - ∆t τ j • ⃗ l m j (t n ) + t n+1 t n C j • exp - t n+1 -ξ τ j • d⃗ e h (ξ ) dξ dξ . (4.38)
Since the problem discretisation relies on the transition between differential values and discrete times by the following relation:

d⃗ e h dt = lim ∆t→0 ⃗ e h (t n+1 ) -⃗ e h (t n ) ∆t , (4.39)
the previously established viscous stress at time t n+1 (Equation 4.38) becomes:

⃗ l m j (t n+1 ) = exp - ∆t τ j • ⃗ l m j (t n ) + t n+1 t n C j • exp - t n+1 -ξ τ j • ⃗ e h (t n+1 ) -⃗ e h (t n ) ∆t dξ . (4.40)
Finally by integrating the previous relation, the exact regressive formula for the current value of the viscous stress quantities is given by:

⃗ l m j (t n+1 ) = exp - ∆t τ j • ⃗ l m j (t n ) + C j • 1 -exp -∆t τ j ∆t τ j • (⃗ e h (t n+1 ) -⃗ e h (t n )). (4.41)
Therefore, this viscoelastic model is set by N families of viscoelastic stiffness tensor C j , related to a relaxation time τ j provided as parameter. The viscoelastic tensor is given by:

C j =             E ve j(1) 0 0 0 0 0 0 E ve j(2) 0 0 0 0 0 0 E ve j(3) 0 0 0 0 0 0 G ve j(23) 0 0 0 0 0 0 G ve j(13) 0 0 0 0 0 0 G ve j(12)             (4.42)
where E ve j(i) and G ve j(ik) are parameters of the model. The generalised Maxwell model to take into account the strain-rate dependency behaviour of the fabric reinforced polymers has quickly emerged. The strain-space formulation of this viscoelastic model is very efficient in an explicit scheme because of the limitation of internal equilibrium loop. By using the strain field provided in finite element code in the material coordinate system and the viscoelastic stress at the previous time step, the actual viscoelastic stresses are updated. The details of the implementation procedure are given in Algorithm 4.1.

Algorithm 4.1 Computation of the viscoelastic stresses.

Step 1 Calculation of the healthy strain ⃗ e h (Equation 4.29).

Step 2 Loading of the previous configuration ⃗ l m j (t n ) and ⃗ e h (t n ).

Step 3 Calculation of the actual viscoelastic stress ⃗ l m j (t n+1 ) (Equation 4.41).

Step 4 Storage of the actual configuration ⃗ l m j (t n+1 ) and ⃗ e h (t n+1 ).

Step 5 Addition of the actual viscoelastic stress to the infinite stress (Equation 4.36).

Identification

Experimental procedure

The identification of the viscoelastic parameters is done by means of Dynamic Mechanical Analysis (DMA). Since the viscoelastic constants are not dependent of the strains, they are determined through the study of small amplitude tensile experiments. The material remains in the elastic part and the identification of the normalised relaxation modulus can be determined without damage correction. These tests are carried out at room temperature on an electromagnetic device (Instron E3000) whith a 3 kN load cell. The specimens are cut using the water-jet technique and their shape and dimensions are shown in Figure 4.8. No heels are used since the deformation remains small and in order not to distort measurements. The distance between the clamps was set to 10 mm. The Dynamic Mechanical Analysis are carried out by means of small amplitude cyclic tensile experiments. The cyclic deformation is introduced by the application of an imposed sinusoidal displacement with an angular frequency ω, which is set as follows:

ε(t) = ε 0 • cos (ωt) = ε 0 • Re {exp (iωt)} (4.43)
where ε 0 is the amplitude, t the time and Re {•} the real part of a complex number. In a permanent regime, the stress response of the material is sinusoidal with the same angular frequency as the one imposed for the strains. However, the stress output is not in phase with the strain input (Figure 4.9) and is thus given by:

σ (t) = σ 0 • cos (ωt + δ ) = σ 0 • Re {exp (i(ωt + δ ))} . (4.44)
By the ratio of stress and strain the complex modulus E * , which is consequently frequency dependent, is obtained:

E * (iω) = σ 0 ε 0 exp (iδ ) . ( 4

.45)

As every complex number, the complex modulus E * can be split into real and an imaginary components. The storage modulus E ′ = Re {E * } corresponds to the real part of the complex modulus which is significant of the pure elastic response (E ′ is in phase with the strain input). Rather, the loss modulus E ′′ = Im {E * } corresponds to the imaginary part of the complex modulus and is significant of the pure viscous response (E ′′ is out of phase with the strain input). By considering the split of the complex modulus given by:

E * = E ′ + iE ′′ , (4.46)
the storage and the loss modulus are thus defined by:

E ′ = σ 0 ε 0 cos (δ ) , (4.47)
E ′′ = σ 0 ε 0 sin (δ ) . (4.48)
The loss angle δ , defined as the ratio between the loss modulus and the storage modulus, namely:

tan (δ ) = E ′′ E ′ (4.49)
is representative of the viscous ability of the material. When the loss angle is null the material is purely elastic, whereas when δ tends towards π 2 the material is purely viscous. The experiments were conducted at nine frequencies between 0.01 and 30 Hz. It was not possible to test the present material for higher frequencies with the actual facility. The carbon fibre reinforced polymers are particularly stiff and machine limitations were reached.

Identification of the viscoelastic constants

From Section 4.2.2, the constitutive relation between the infinitesimal strain and the infinitesimal stress for an uni-axial model is given by:

σ (t) = t -∞ E(t -ξ ) • dε(ξ ) dξ dξ . (4.50)
By using the definition of the input sinusoidal strain ε(t) (Equation 4.43) in a complex form and by substituting it in the previous equation, the stress response is given by:

σ (t) = iωε 0 t -∞ E(t -ξ ) • exp (iωξ )dξ . (4.51)
Let u = tξ , with u = ∞ when ξ = -∞ and u = 0 when ξ = t, changing variables in the previous equation yield to:

σ (u) = iω • ε 0 exp (iωt) ∞ 0 E(u) • exp (-iωu)du = iω • ε(t) ∞ 0 E(u) • exp (-iωu)du. (4.52)
After replacing the variable u by t the complex modulus can be extracted such as:

E * = ∞ 0 E(t) • exp (-iωt)dt. (4.53)
E(t) is the linear relaxation modulus coming from the generalised Maxwell model and is formulated in terms of the Prony series (Section 4.2.1):

E(t) = E ∞ + N ∑ j=1 E j • exp - t τ j , (4.54)
with E ∞ the long term modulus, E j the elastic stiffness and τ j the relaxation time of the j th viscous mechanism. Consequently, the complex modulus is given by:

E * = E ∞ + N ∑ j=1 E j (ωτ j ) 2 1 + (ωτ j ) 2 + i N ∑ j=1 E j ωτ j 1 + (ωτ j ) 2 = E ′ + iE ′′ . (4.55)
The storage and the loss modulus are therefore expressed as:

E ′ (ω) = E ∞ + N ∑ j=1 E j (ωτ j ) 2 1 + (ωτ j ) 2 , (4.56)
E ′′ (ω) = N ∑ j=1 E j ωτ j 1 + (ωτ j ) 2 .
(4.57)

N couples of viscoelastic stiffness E j and relaxation time τ j are identified such that the average square deviation between the predicted moduli (by using Equations 4.56 and 4.57) and the measured storage E ′ exp and loss E ′′ exp moduli at M frequencies ω k such as: min

E j ,τ j M ∑ k=1   E ′ (ω k ) E ′ exp -1 2 + E ′′ (ω k ) E ′′ exp -1 2   (4.58)
is minimum.

The results of the parameter identification of the viscoelastic constants of the various preforms, by considering seven relaxation times, are provided in Table 4.1.

Identification of the intralaminar matrix damage parameters

Because the matrix damage is only acting on the long-term stiffness tensor (Section 4.3) which does not correspond to the initial elastic stiffness of the material, the parameters which are given in Table 3.1 cannot be used for the present model.

As a consequence, after the identification of the viscoelastic parameters, the intralaminar matrix damage parameters need to be re-evaluated with the same procedure as described in Section 3.5. Nothing change, except the initial (visco)elastic parameters now provided by means of DMA. The re-evaluated matrix damage parameters are given in Table 4.2.

Validation

Tensile tests with standard specimens

In order to evaluate the ability of the model to represent the strain rate sensitivity of the model, dynamic loading tests have been carried out at various strain rates. A hydraulic high speed device (Instron 65/20) with a 30 kN load cell is used for these tests. In order to ensure as much as possible a constant strain rate, the device is fitted with a piston system the stroke of which is used during speed increasing. Usually, specimens with reduced dimensions are used for dynamic testings. The dynamic facilities are not able to sustain the same maximal amount of force as the static ones. For this purpose, the specimen width is reduced to limit the specimen area. Furthermore, at constant loading speed the strain rate is artificially increased by reducing the specimen length.

0 0 0 h II m(1) = h II m(2) 1 1 1 h III m(1) = h III m(2) 0 0 0 d p c(3) / 0.
Berthe [START_REF] Berthe | Comportement thermo-visco-élastique des composites CMO -De la statique à la dynamique grande vitesse[END_REF] has shown that in case of fabric reinforced materials, the specimen length had to be at least two times greater than the specimen width. However the fabric reinforced polymers force to maintain a minimal width in order to present enough representative volume elements. An other aspect of the fabric reinforced polymers is critical in the choice of the coupon width: these materials show important edge damage and reduce the effective cross-section. If the ratio between the effective-cross section and the initial cross-section of the specimen is not maintained to an acceptable level, the results become non-representative of the real material behaviour (Figure 4.10). Despite all the technical drawbacks, the dynamic specimens are designed to be identical to the coupons used for the standardised quasi-static tensile tests (Figure 4.11) and are cut using the water-jet technique. Because of the original specimen dimensions for dynamic tests, none suitable clamping device were existing. Thus, as part of that work, a new clamping system has been designed. Three requirements specifications were formulated:

• to clamp standardised specimens;

• to apply homogenised clamp pressure in order to avoid a premature failure in the heels; • to reduce the mass as much as possible in order to limit the inertial effects. The adopted design solution consists of a prismatic clamping, where the clamping is carried out by screwing of the external parts. The Computer-Aided-Design-based digital mock-up is shown Figure 4.12. Since the strain rate sensitivity occurs mainly for shear loading (according to the material coordinate system) and due to the limitations of the hydraulic high speed device, only tensile tests at 45°according to the fibre directions are carried out. Except the speed loading these dynamic tests are similar to the NF EN ISO 14129 standard. They are carried out at room temperature and at 3 speed loadings (1.7, 41 and 1000 mm s -1 ) which correspond to approximate equivalent strain rates of 1 10 -2 , 3 10 -1 and 6.5 s -1 . Three coupons per speed loading are used to evaluate the repeatability of the results. The hydraulic high speed device, the clamping system and the optical measurements (laser extensometer and high-speed camera) device are shown under real conditions in Figure 4.13. The results were conclusive as no heel fracture occurred and as the system vibrations were sensibly contained.

Comparison of the numerical model with the experiments

The present damageable viscoelastic model is implemented in FORTRAN 90 in a user material subroutine as described section 4.5 for the explicit finite element code LS-DYNA ® .

The parameter identification of the viscoelastic model is carried out by using seven Maxwell elements. Among the studied fabric preforms, only the bi-axial non-crimp fabric and the 3K woven fabric were tested through Dynamic Mechanical Analysis (DMA). The 12K woven fabric presents a unit-cell of 15 mm width which does not allow to carry out a representative DMA test. Compared to other tests where the specimen width includes or is enlarged to include at least one unit-cell of the textile pattern, the load capacity of the electromagnetic device is too limited to enlarge the DMA specimens. As a consequence, the viscoelastic constants of the 12K woven fabric preform are considered identical to the 3K woven fabric.

This section presents the results of the identification procedure which have been presented Sections 4.6.2 and 4.6.3 on two various fabric preforms, namely a bi-axial non-crimp fabric (see Figure 2.5) and a 3K plain-weave woven (see Figure 2.6). The optimisation procedures were carried out thanks to the commercial software Matlab ® to minimise the Mean Square Error between the numerical and the experimental discrete relaxation spectra. The comparison between the viscoelastic model responses and the experimental data are shown Figure 4.14 for the storage modulii and Figure 4.15 for the loss modulii. As the storage modulii are much higher, and consequently being more influential on the mechanical response, than the loss modulii, the choice was done to better fit the evolution of the storage modulii.

As validation procedure, the dynamic tensile tests (presented Section 4.7.1) are simulated. Before comparing the numerical results with the experimental tests, some audits of the numerical model were carried out. Because of the sudden acceleration of the test device in case dynamic loading, the eigenfrequencies of the system are excited which can lead to a parasite noise in the measured response. Thus, a modal analysis of the system (high speed jack piston, clamping device, specimen and load cell) was made. The results demonstrate that the low eigenfrequencies (3.1 and 6.7 kHz) which may distort the force measurements are due to the mass of the clamping devices, and the stiffness of the specimen and the piston (Figure 4.16). Thus, unlike the simulations of the quasi-static tests, the model simplification by applying directly the loading rate on the specimen need to be verified. A first simulation of dynamic in-plane shear loading at 1000 mm s -1 , which includes the piston, the clamping devices and a bi-axial non-crimp fabric specimen, has been carried out (Figure 4.17). The vibration predicted by the modal analysis occurs and leads to a slightly noisy response as expected. Another simulation have been carried out, but this time on a single specimen clamped on both sides (Figure 4.18). The force response of both simulations are then compared and the results are shown (Figure 4. [START_REF] Buckley | Deformation of thermosetting resins at impact rates of strain. part i: Experimental study[END_REF]. Note that the displacements are obtained by placing a virtual extensometer in the middle of the specimen. Thus, the results are well compared in an area which is homogeneously deformed. Excepted the noise on the force response of the full device, the results are significantly close. Therefore, the following simulations will be carried out on a single specimen.

The results of the dynamic in-plane shear tests carried out on three various fabric reinforced polymers are given fabric preform) and Figure 4.22 (12K woven fabric preform). The oscillations on the numerical responses at 1000 mm s -1 is due to the use of the real displacements given by the experimental tests as input of the simulation. These displacements being noisy, the resulting responses are noisy too. Also, it is essential to recall that the coupling between the viscoelastic model and the matrix damage model is done without additional parameter. Moreover, the parameters of both modules of the model are identified independently: quasi-static cyclic tensile tests for the matrix damage and Dynamic Mechanical Analysis for the viscoelasticity. The viscoelastic parameters lead to an approximation of the real behaviour of the fabric reinforced polymers at various strain rate and need to be adapted for a good correlation with the experiments. The Dynamic Mechanical Tests are suffering of several drawbacks which can lead to the identification of non-representative parameters of the real strain-rate behaviour. The load capacity of the device is limited and the investigated materials are very stiff. Moreover, an important self-heating of the specimen occurs during DMA testing. Because of the epoxy response is temperature dependent, the identification may be distorted. By keeping in mind this context, the correlation between the experiments and the numerical model is good.

Due to a mistake in the calibration of the load cell, no data for the bi-axial non-crimp fabric material at 41 mm s -1 are available. Even so, the simulation response at this speed loading is plotted Figure 4.20. For the 3K woven fabric, the identification was more challenging. The quasi-static cyclic in-plane shear tests used to identify the matrix damage were not carried out in the same time that the dynamic test campaign. Thus, as seen Figure 4.21, an aging or environmental effects influence the material behaviour. As a consequence, the responses are too stiff and the viscoelastic parameters need to be re-evaluated.

Concluding remarks on the strain-rate dependency

The strain-rate sensitivity of the fabric reinforced polymers is greatly induced by the viscoelastic behaviour of the resin used as matrix. Among the various viscoelastic model, the generalised Maxwell model, also called Standard Linear Solid model, is used in the material law. This model offers the advantages of a very simple implementation into a finite element code, and most important a proved efficiency for explicit finite element simulations. This model avoids the use of internal equilibrium loops inside the material model. However, the generalised Maxwell model is a linear viscoelastic model, while the fibre reinforced polymers have proven that they present a non-linear viscoelasticity. Therefore, since the strain-rate dependency is essentially due to the matrix, the viscoelastic model is coupled with the matrix damage model. This coupling is introduced by considering an healthy area, only potential location where the matrix is stressed, because other mechanisms (such as friction) are predominant in the damage area and lead to a full relaxation of the matrix around the cracks (the cracks lips movements are preferred to the matrix strains).

The identification of the generalised Maxwell model is done by widespread Dynamic Mechanical Analysis. None additional parameter is needed for the coupling with the matrix damage.

This model is then validated through dynamic in-plane shear tests by using an high-speed jack device. In order to limit the side effects, very important in case of fabric composite, the standardised dimensions of the specimens for quasi-static experiments are maintained. A new clamping system have been designed to allow dynamic tests on standardised coupons. Moreover, this clamping device ensures an homogeneous pressure which avoids the premature failure in the heels. The simulations of the validation tests show good correlation on the various preforms.

The procedure for parameter identification of the viscoelastic model suffers of numerous drawbacks. The facility used for DMA have a limited load capacity and some materials cannot be tested (such as the 12K). Moreover, the frequency is limited to 30 Hz which restrained the identification of the viscoelastic spectrum to relatively low strain rates (less than 10 s -1 ). By using a time/temperature dependence, it could be possible to identify other parts of the spectrum by using a climatic chamber. Also, this temperature dependence could be implemented directly in the viscoelastic model. However, the vibrations during DMA tests lead to important self-heating which distorts the results. The identification of the viscoelastic parameters of the smaller relaxation times are therefore re-calibrated by using the dynamic in-plane shear tensile tests. In the future, another procedure to identify the viscoelastic parameters, such as creep tests or the SE Ė method which consists in plotting a three-dimensional surface given by the stress response and depending on the strain and the strain-rate, has to be investigated [START_REF] Lauro | The SEE method for determination of behaviour laws for strain rate dependent material: Application to polymer material[END_REF]. maximal strain criterion. The failure occurs when the component of the stress tensor or the component of the strain tensor reaches a threshold value. However, these failure criteria do not take into account the possible interactions between the failure modes. Hence in case of multi-axial loading the material strength is overestimated. However, due to their simplicity, these criteria are widely used for preliminary studies.

The energetic criteria such as Hill [START_REF] Hill | A Theory of the Yielding and Plastic Flow of Anisotropic Metals[END_REF] or the more general tensor polynomial failure criteria such as Tsai and Wu [START_REF] Tsai | A General Theory of Strength for Anisotropic Materials[END_REF] are based on the quadratic coupling between the different values of the stress tensors. They are expressed in the form:

F i j • σ i σ j + F i • σ i ≤ 1 (5.1)
where F i j and F i are the failure parameters of the material. However the quadratic formulation leads to some constraints on the parameters and can prevent an acceptable identification. The parameters should ensure the convexity of the failure envelope.

The last family of failure criterion relies on the uncoupling of the different failure mechanisms which can occur in the material. Yeh and Kim [START_REF] Yeh | The Yeh-Stratton Criterion for Composite Materials[END_REF] provides a first multiphysic criterion for the unidirectional composites. In last decades, the prediction of the failure of unidirectional plies has been greatly enhanced. Complex failure mechanisms, such as the fibre kinking, and of their interactions are nowadays well considered and predicted. It leads to efficient criterion for the unidirectional plies such as Puck [START_REF] Puck | Failure analysis of FRP laminates by means of physically based phenomenological models[END_REF] or LaRC04 [START_REF] Pinho | Failure models and criteria for frp under in-plane or three-dimensional stress states including shear non-linearity[END_REF]. Laurin et al. [START_REF] Laurin | A multiscale progressive failure approach for composite laminates based on thermodynamical viscoelastic and damage models[END_REF] introduces the influence of the matrix viscoelasticity on the final fracture and a failure criterion for the laminates. In the same vein, the NU criterion established by Daniel et al. [START_REF] Daniel | Strain-rate-dependent failure criteria for composites[END_REF] is given. Regarding the fabric reinforced polymers the state of the art is much more limited. Aiello [START_REF] Aiello | Utilisation des composites à matrice céramique SiC/SiC comme matériau de structure de composants internes du tore d'un réacteur à fusion[END_REF], provides a multi-envelope failure criterion for the fabric reinforced ceramics. However, the failure mechanisms are not the same between a resin matrix and a ceramic matrix. Other criteria, more adapted to the fabric reinforced polymers, were proposed by Key et al. [START_REF] Key | Progressive failure modeling of woven fabric composite materials using multicontinuum theory[END_REF] or Mallikarachchi and Pellegrino [START_REF] Mallikarachchi | Failure criterion for two-ply plain-weave CFRP laminates[END_REF].

However, these last multi-physics criteria are hard to precisely identify. They require specific experimental tests which are not easy to integrate in an industrial framework. Moreover, the maximum stress (or its dual form the maximal strain criterion), although this is a simplistic failure criterion, proves to be efficient and is in fact relatively accurate in case of the bi-axial fabric reinforced polymers. Consequently, this criterion is widely used in the literature [START_REF] Marcin | Modélisation du comportement, de l'endommagement et de la rupture de matériaux composites à renforts tissés pour le dimensionnement robuste de structures[END_REF][START_REF] Cousigné | Numerical modeling of nonlinearity, plasticity and damage in CFRP-woven composites for crash simulations[END_REF].

Once the fibre fracture initiation is predicted, the propagation needs to be taken into account. The failure, although sudden, is not immediate. Due to the amount of stiffness loss, the stress reaches locally a maximal value at the time of the failure initiation of the reinforcement. However, in the case of an implementation of the material law for finite element analyses, this material softening leads to a spurious mesh dependence phenomenon with a strain localisation. The loss of hyperbolicity of the governing equilibrium equations in case of softening creates this mesh dependency. Since the fibre failure leads to the degradation of the stiffness until the loss of the structural integrity, various numerical methods are able to limit the mesh dependency by managing this final stiffness erosion. Among these methods, three families of post-failure modelling techniques can be mentioned:

• non-local formulations,

• delay effect on the damage variables,

• smeared formulations.

An integral type non-local model [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Strömberg | FE-formulation of a nonlocal plasticity theory[END_REF][START_REF] Jirásek | Comparison of integral-type nonlocal plasticity models for strain-softening materials[END_REF][START_REF] Balieu | Modèle viscoélastique-viscoplastique couplé avec endommagement pour les matériaux polymères semi-cristallins[END_REF] consists in defining a damage variable by a spatial averaging. Thus the hypothesis of a local behaviour, which states that the thermodynamic potential in a point is only dependent on internal and external thermodynamic variables at the same point, is discarded. A diffusion model equation, relying on a reference length, is used. However in finite element analysis the non-local formulations require variables provided by several elements, or the development of adapted elements is not provided in commercial finite element codes.

The second approach consists in the limitation of the damage rate [START_REF] Ladeveze | Damage modelling of the elementary ply for laminated composites[END_REF][START_REF] Ladeveze | A damage computational approach for composites: Basic aspects and micromechanical relations[END_REF][START_REF] Allix | A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues[END_REF]. It can be physically understood by the fact that a crack tip cannot spread at infinite speed. However, although this approach limits the mesh dependency, it introduced a strain-rate dependency which increases the dissipated energy during failure. This viscous behaviour may be unrepresentative of the true material behaviour. Marcin [START_REF] Marcin | Modélisation du comportement, de l'endommagement et de la rupture de matériaux composites à renforts tissés pour le dimensionnement robuste de structures[END_REF] coupled both methods, non-local and delay effect, to approach in an accurate way the real fracture mechanisms, while ensuring the mesh independence.

The last method, the smeared crack formulation [START_REF] Pinho | Physically-based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking: Part I: Development[END_REF][START_REF] Cervera | Smeared crack approach: back to the original track[END_REF][START_REF] Camanho | Modeling the inelastic deformation and fracture of polymer composites -Part II: Smeared crack model[END_REF][START_REF] Cousigné | Numerical modeling of nonlinearity, plasticity and damage in CFRP-woven composites for crash simulations[END_REF], consists of the distribution of the fracture energy over the volume element. By knowing the fracture toughness and the element size it becomes possible to regularise the dissipated energy and to limit the effects of the mesh size on the failure behaviour. By its simplicity, this numerical method of damage regularisation is efficient in case of finite element simulations to limit the mesh dependence.

In this work, a modified maximal strain criterion is used to predict the failure initiation. This criterion, although simple, is rather accurate in case of the fabric reinforced polymers. Its formulation and the parameter identification procedure are given Section 5.2. Regarding the evolution law of the fibre damage, a smeared crack formulation is chosen. Its efficiency and its ease of implementation into the commercial finite element codes are great assets.

Moreover, this stiffness erosion is able to limit the mesh dependence due to the material softening which occurs at failure. The smeared crack approach and the original coupling with the global framework used in this work is described in Section 5.3. The implementation and some example of simulations with the present fibre failure model are then provided.

Maximal strain failure criterion 5.2.1 Formulation and implementation in the current framework

Because the final fracture of the bi-axial fabric reinforced polymer is due to the fibre failure, even for shear loading, the choice was made to use a simple failure criterion which relies on the strains in the reinforcement directions. The total Lagrangian framework (Section 2.4) used in this work is particularly suitable to the fibre failure prediction since the coordinate system axis are coincident with the material axis, namely the fibre directions, at any time. Hence, the longitudinal and the transverse components of the Green-Lagrange strain tensor E, respectively E 11 and E 22 are considered equivalent to the strains in the fibre directions. Note that in the interest of writing simplification, the Voigt notation which reduces the number of components of a symmetric tensor, are used for the Green-Lagrange strain tensor: (5.3)

E =    E 11 E 12 E 13 E 21 E 22 E 23 E 31 E 32 E 33    ( 5
Thus the longitudinal and transverse strains in the fibre directions are now respectively given by the components e 1 and e 2 .

By differentiating the tensile and the compressive fibre failures, four parameters are needed to establish the failure envelope: the longitudinal tensile failure strain e ft 1 , the longitudinal compressive failure strain e fc 1 , the transverse tensile failure e ft 2 and the transverse compressive failure strain e fc 2 . Please remind that unlike the maximal strain criterion applied to the unidirectional plies, no maximal shear strain criterion is used in the case of the bi-axial fabric reinforced composites. As a result, the maximal strain criterion is given by: Longitudinal fibres e fc 1 < e 1 < e ft 1

Transverse fibres e fc 2 < e 2 < e ft 2 (5.4) But this criterion lacks the consideration of shear strain, especially for mixed compressive/shear loading. The compressive fibre failure is generally induced by kinking, which can be seen as a local buckling of the yarns, and which is eased by shearing. The influence of shear loading is introduced by using the spectral decomposition of the strain tensor. In case of tensile loading, the fracture is predicted by using the positive part of the spectral decomposition, given by:

E + = P •    ⟨E I ⟩ 0 0 0 ⟨E II ⟩ 0 0 0 ⟨E III ⟩    • P T (5.5)
where E I , E II and E III are the eigenvalues of E and P the transformation matrix formed by the eigen vectors. Lastly e + i are the components of ⃗ e + , namely the Voigt notation of E + . Otherwise, the compressive fracture is predicted by using the negative part of the spectral decomposition, given by:

E -= P •    ⟨-E I ⟩ 0 0 0 ⟨-E II ⟩ 0 0 0 ⟨-E III ⟩    • P T . (5.6)
Here, e - i are the components of ⃗ e -, namely the Voigt notation of E -. Moreover, in order to ensure consistency between the matrix damage and the fibre damage formulations, the criterion will be expressed according to the thermodynamic forces associated to a fibre crack. They are defined by the following relations: Both fibre directions are tested (Figure 5.2) in order to obtain the longitudinal and the transverse tensile Green-Lagrange strains at failure, respectively e ft 1 and e ft 2 . However, the ultimate strain in the fibre directions is less than 2% and the small strain formalism (engineering strain ε and engineering stress σ ) can be used without compromise on the accuracy. In both cases and due to small strains, a gauge has been used to measure the longitudinal strain in the specimen coordinate system ε x = ε( ⃗ X), with ε( ⃗ X) the elongation in the centre of the coupon according to the ⃗ X direction. The force applied to the coupon F = ⃗ F • ⃗ X is measured with a 100 kN load cell and then the stress is obtained by:

Longitudinal fibres            Tension y ft 1 = 1 2 C 0 11 e 2 1 + c ft 1 • e + 1 2 c ft 1 + 1 Compression y fc 1 = 1 2 C 0 11 e 2 1 + c fc 1 • e - 1 2 c fc 1 + 1 Transverse fibres            Tension y ft 2 = 1 2 C 0 22 e 2 2 + c ft 2 • e + 2 2 c ft 2 + 1 Compression y fc 2 = 1 2 C 0 22 e 2 2 + c fc 2 • e -
⃗ X = ⃗ 1 ⃗ Y = ⃗ 2 (a) Longitudinal tensile test. ⃗ X ⃗ Y = ⃗ 1 ⃗ 2 90°( b) Transverse tensile test.
σ x = F S 0 (5.10)
where S 0 is the initial cross-section of the coupon. The initial longitudinal stiffness E 0 x is then defined as the initial slop of the stress-strain curve (σ x = f (ε x )).

The strain at failure ε ft

x is then defined as the strain when the maximal force is reached, which can be expressed as follows:

ε ft x = ε x (t ft ) with F(t ft ) = max t F(t).
(5.11)

The thermodynamic forces associated to the tensile fibre failure y ft 1 and y ft 2 are then obtained by following the Equation 5.8 whith:

Longitudinal tensile test e ft 1 ≈ ε ft x C 0 11 = E 0 x Transverse tensile test e ft 2 ≈ ε ft x C 0 22 = E 0 x .
(5.12)

In-plane shear tests

The test configuration of the in-plane shear tests are the same as the tensile tests previously described. Only the direction of loading according to the fibre direction is changing (Figure 5.3).

⃗ X ⃗ Y ⃗ 1 ⃗ 2 45°F
igure 5.3: Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for in-plane shear tests used for failure identification.

Because of large strain for in-plane shear tests of fabric reinforced polymers, extensometers have been used to measure the longitudinal strain in the specimen coordinate system ε x = ε( ⃗ X), with ε( ⃗ X) the elongation in the centre of the coupon according to the ⃗ X direction, and the transverse strain in the specimen coordinate system ε y = ε( ⃗ Y ), with ε( ⃗ Y ) the elongation in the centre of the coupon according to the ⃗ Y direction.

Since the model is expressed in the total Lagrangian framework the Green-Lagrange strain tensor have to be defined from the measurements:

F = 1 + ε X 0 0 1 + ε Y . (5.13)
After a change of basis to be placed in the material frame, the tensor becomes:

F =     1 + ε X + ε Y 2 ε X -ε Y 2 ε X -ε Y 2 1 + ε X + ε Y 2     .
(5.14)

By using Equation 2.9, the total Lagrangian tensor is determined in the material coordinate system:

E m = 1 2     ε X + ε Y + ε 2 X + ε 2 Y 2 ε X -ε Y + ε 2 X -ε 2 Y 2 ε X -ε Y + ε 2 X -ε 2 Y 2 ε X + ε Y + ε 2 X + ε 2 Y 2     .
(5.15)

Thus, it becomes possible to compute E 11 and E + 11 by spectral decomposition (or e 1 and e + 1 in Voigt notations). By using Equations 5.7 and 5.8 the parameter c ft i is finally given by:

c ft i = e 2 i -e ft i 2 e ft i 2 -e + i 2 ∀i ∈ {1, 2}. (5.16)
The parameters c f c i can be identified by following the same procedure but by using in-plane compressive shear tests. Theses tests have not been carried out in this work.

Compression tests

These tests are carried out at room temperature on an electromagnetic device (Sintech 20D) which ensures 100 kN for the maximal load capacity and at a speed of 1 mm min -1 . The specimens are cut using the water-jet technique and their shape and dimensions are shown in As in the case of tensile tests, both fibre directions are tested (Figure 5.6) and the ultimate strains remaining at value less than 2%, the small strain formalism is used. In case of a balanced fabric, one direction is sufficient to identify the failure criterion. In order to measure the strain and to control the eventual buckling of the coupon, a gauge is glued on both sides of the coupons. A buckling criterion, which states when a test is valid, is thus defined by:

⃗ X = ⃗ 1 ⃗ Y = ⃗ 2 (a) Longitudinal compression test. ⃗ X ⃗ Y = ⃗ 1 ⃗ 2 90°( b) Transverse compression test.
r buck = max ε b x -ε a x ε b x + ε a x ≤ 0.2 (5.17) with • ε a x = ε a ( ⃗ X)
the elongation in the centre of the coupon according to the ⃗ X direction given by the gauge on the front side,

• ε b x = ε b ( ⃗ X)
the elongation in the centre of the coupon according to the ⃗ X direction given by the gauge on the back side.

The problem of buckling for compression tests can be assessed Figure 5.7. The strain of the specimen is then given by the mean of both strains provided by the gauges:

ε x = ε a x + ε b x 2 .
(5.18)

Finally, the same procedure as that used for the tensile failure (Equations 5.10 to 5.12) is used to obtain the thermodynamic forces associated to the compressive fibre failure, namely y fc 1 and y fc 2 .

5.3 Smeared crack approach for the final stiffness erosion

Basic concepts of the smeared crack formulation

During finite element analysis, the softening due to the fibre failure induced a strain localisation on the element where the fibre damage appears. Hence, the energy dissipation by the fibre failure is dependent on the element size. The method of regularisation, usually named the smeared crack method proposed by Bažant and Oh [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF] and used in a first attempt for the fibre reinforced materials by Pinho et al. [START_REF] Pinho | Physically-based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking: Part I: Development[END_REF], consists of ensure a constant energy released per unit area, independently from the element dimensions.

Let U the energy dissipated by one element during the fibre failure. By considering an element with a volume V , the total energy dissipated per unit volume is defined as follows:

U vol = Γ ∞ + Γ dam = U V . ( 5

.19)

As represented Figure 5.8, Γ ∞ is the dissipated energy in case of instantaneous fibre failure and corresponds to the energy stored by the elastic deformation in the fibre direction, and Γ dam is the additional energy which is dissipated by a progressive fibre failure. By considering the area of the fracture plane A, transverse to the fibre direction, a dissipated energy per unit area Γ, namely the fracture toughness which is a material constant, is defined by the relation:

U = AΓ, (5.20)
which leads to the following definition:

Γ = U vol V A = (Γ ∞ + Γ dam )l (5.21)
where l is the element characteristic length. Since l represents the length of the element in the crack direction, and as the crack may feature a random propagation, an averaging of the element characteristic length is used to remain as general as possible. Therefore, Bažant and Oh [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF] defined the characteristic element length of a four nodes square element by: l = 1.12 √ A IP (5.22) with A IP the area associated with each integration point. Also Pinho et al. [START_REF] Pinho | Physically-based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking: Part I: Development[END_REF] defined the characteristic element length for triangular elements:

l ≈ 1.52 √ A IP .
(5.23)

Evolution of the fibre damage and formulation in the current framework

In order to simplify the subsequent equations and enhance the understanding, the adopted formalism for the establishment of the fibre damage model is based on the matrix damage model. The viscoelastic coupling will be briefly given at the end of this section.

Consequently, the constitutive equation after introduction of the fibre damage is given in the total Lagrangian framework and by using the Voigt notation by:

⃗ s = C • (⃗ e -⃗ e 0 ) -C f • (⃗ e r +⃗ e s -⃗ e 0 ) (5.24)
where the effective stiffness tensor C takes into consideration both matrix and fibre damages and is defined as follows:

C = S 0 + ∆S m + ∆S f -1 (5.25)
with S 0 the elastic stiffness tensor, ∆S m the additional compliance due to the matrix damage and S f the additional compliance due to the fibre damage. C f corresponds to the effective stiffness tensor due to the fibre failure such as:

C f = S 0 + ∆S f -1 . ( 5.26)
As a recall, ⃗ e 0 corresponds to the strain vector such as the residual stress induced by the difference between the thermal dilatation properties of the constituents, ⃗ e r denotes the eventual plastic deformations around the matrix cracks and ⃗ e s depends of the position of the matrix crack lips.

The compliance induced by the fibre damage is expressed according to fibre damage variables d f i and compliance tensor H f i associated to the damage variable d f i , such as:

∆S f = ∑ i d f i H f i . (5.27)
In the case of bi-axial fabric preforms, two damage variables are used to describe two fibre directions. Therefore, the compliance tensors associated to the damage variables are defined by:

H f 1 =            h I f (1) S 0 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h III f (1) S 0 55 0 0 0 0 0 0 h II f (1) S 0 66            , (5.28)
and

H f 2 =            0 0 0 0 0 0 0 h I f (2) S 0 11 0 0 0 0 0 0 0 0 0 0 0 0 0 h III f (2) S 0 44 0 0 0 0 0 0 0 0 0 0 0 0 0 h II f (2) S 0 66            , (5.29)
where

h I f (i) , h II f (i) and h III f (i)
are parameters of the model. The damage evolution is defined by independent criteria, which ensure to verify the second law of thermodynamics, on every damage variables, such as:

Longitudinal fibres Tension F ft 1 (y ft 1 , d f 1 ) = f ft 1 (y ft 1 ) -d f 1 ≤ 0 Compression F fc 1 (y fc 1 , d f 1 ) = f fc 1 (y fc 1 ) -d f 1 ≤ 0 Transverse fibres Tension F ft 2 (y ft 2 , d f 2 ) = f ft 2 (y ft 2 ) -d f 2 ≤ 0 Compression F fc 2 (y fc 2 , d f 2 ) = f fc 2 (y fc 2 ) -d f 2 ≤ 0 . (5.30)
The damage evolution law f fx i , with i ∈ {1, 2} and x being replaced by the symbols "f" for the tensile damage and "c" for the compressive damage, is set by the function:

f fx i (y fx i ) = κ fx i < y fx i -y fx 0(i) > (5.31)
with κ fx i a constant defining the damage rate, y fx i the thermodynamic force (Equation 5.7) and y fx 0(i) the thermodynamic force at failure initiation based on a close version of the maximal strain criterion (Equation 5.8).

The coupling between the current fibre damage formulation and the smeared crack approach leads to find the appropriate value of the constant κ fx i . To simplify the equations, a uni-axial case is used to lead to the definition of the damage rate and the matrix damage is not considered. Therefore, the constitutive equation between the stress and the strain becomes:

σ = C • (ε -ε 0 ) -C f • ε 0 11 = C f • ε. (5.32)
By introducing the damage variable d and the damage effect h affecting the elastic stiffness C 0 , it yields to:

σ = 1 (1 + d • h) C 0 • ε. (5.33)
Then the additional energy which is dissipated by a progressive fibre failure Γ dam (Figure 5.8) is obtained by integrating stress on the strain interval [ε f0 , ε fF ], with ε f0 the strain at failure and ε fF the strain at integration point deletion or computation termination:

Γ dam = ε fF ε f0 σ • dε = ε fF ε f0 C 0 1 + d • h • ε • dε. (5.34)
The definitions of the damage variable and of the thermodynamic force yield to:

Γ dam = ε fF ε f0 C 0 1 + κ f • h • ( 1 2 C 0 • ε 2 -y f 0 ) • ε • dε. = ε fF ε f0 C 0 1 + κ f • h • ( 1 2 C 0 • ε 2 -y f 0 ) • κ f • h •C 0 κ f • h •C 0 • ε • dε. = 1 κ f • h ln(1 + κ f • h • ( 1 2 C 0 • ε 2 -y f 0 ) ε fF ε f0 = 1 κ f • h • ln(1 + κ f • h • ( 1 2 C 0 • ε fF 2 -y f 0 ) (5.35)
In order to regularise the strain at element deletion ε fF , a dependence to the element size is introduced by the relation:

ε fF = l ref l • ε f,ref F (5.36)
where ε fF,ref is a parameter of the model and l ref is a reference size of element for which the strain at element deletion is equal to the parameter ε fF,ref . By using Equations 5.21 and 5.36 in Equation 5.35, a relation between the damage rate κ f , the element size l and the parameters of the model is given by:

k(κ f ) = ln 1 + κ f • h • 1 2 C 0 • l ref 2 l 2 • ε f,ref F 2 -y f 0 -κ f • h • Γ l -Γ ∞ = 0. ( 5.37)
By taking back the case of a bi-axial fabric preforms and by analogy with the previous uni-directional case, the various fibre damage rate κ ft 1 , κ fc 1 , κ ft 2 and κ fc 2 are given by the following relation:

k fx i (κ fx i ) = ln 1 + κ fx i • h I f (i) • 1 2 C 0 ii • l ref 2 l 2 • ε fx,ref F(i) 2 -y fx 0(i) -κ fx i • h I f (i) • Γ f x(i) l -Γ ∞ f x(i) = 0. (5.38)
With i ∈ {1, 2} and x being replaced by the symbols "f" for the tensile damage and "c" for the compressive damage, the thermodynamic force at damage initiation y fx 0(i) , the strain at element deletion ε fx,ref F(i) , the fracture toughness Γ f x(i) and the reference element size l ref are parameters of the fibre damage model. The dissipated energy in case of instantaneous failure

Γ ∞ f x(i)
is computed during calculation or can be considered as parameters which is determined by a preliminary calculation.

Then κ fx

i is determined thanks to the Newton-Raphson method by the iterative algorithm such as:

κ fx i (n + 1) = κ fx i (n) - k fx i (κ fx i (n)) k fx i ′ (κ fx i (n)) (5.39)
which is repeated until a sufficiently accurate value is reached. Note that k fx i ′ (κ fx i (n)) is given by:

k fx i (κ fx i ) = h I f (i) • 1 2 C 0 ii • l ref 2 l 2 • ε fx,ref F(i) 2 -y fx 0(i) 1 + κ fx i • h I f (i) • 1 2 C 0 ii • l ref 2 l 2 • ε fx,ref F(i) 2 -y fx 0(i) -h I f (i) • Γ f x(i) l -Γ ∞ f x(i) = 0.
(5.40)

Identification of the post-failure parameters

The accurate identification of the post-failure parameters, namely the critical damages d fx c(i)

and the fracture toughnesses Γ f x(i) , with i ∈ {1, 2} and x being replaced by the symbols "t" for the tensile damage and "c" for the compressive damage, is quite challenging. Whereas the critical damage d fx c(i) influences only the strain at integration point deletion, the identification of the fracture toughnesses Γ f x(i) requires non standardised tests. Three different methodologies can be used: four-points bending with pre-crack, tensile or compression on specimens with open holes or the use of Compact Tension (CT) or Compact Compression (CC) specimens.

Pinho et al. [START_REF] Pinho | Developing a four point bend specimen to measure the mode I intralaminar fracture toughness of unidirectional laminated composites[END_REF] proposed a four point bending test with a pre-crack on the side in tension and at the position of the maximal deflection (Figure 5.9). From the total applied load at failure and from the nominal dimensions of the specimens the critical energy release rate can be obtained. [START_REF] Soutis | Failure Prediction Technique for Compression Loaded Carbon Fibre-Epoxy Laminate with Open Holes[END_REF]. The damage localisation leads to a progressive failure of the specimen around the hole. As the failure is not immediate, the dissipated energy can be determined. The last solution, which is an optimisation of open-hole testing as presented by Pinho et al. [START_REF] Pinho | Fracture toughness of the tensile and compressive fibre failure modes in laminated composites[END_REF], used specific thick specimens to precisely identify the fracture toughness (Figure 5.10). However, the need of thick specimens is not possible for the investigated materials. Only thin plates are available.

Another solution is the use of tension or compression tests on open-holes specimens

In the case of this work, none of these tests have been carried out. As a consequence, the identification of the fracture toughness energy is not carried out.

Implementation

In the present model, the fibre damage is evaluated at each integration point. If a fibre damage variable d f i reaches a critical value set by the smeared crack model, then the integration point is flagged as failed. The managing of the element deletion according to the number of failed integration point is set directly in the finite element code. In addition to the strain tensor, the characteristic element length is needed in order to compute the dissipated energy per volume unit. In the explicit finite element code LS-DYNA ® the element size is provided as input parameter of the user material subroutine. Then the damage calculation follows the procedure given by Algorithm 5.1. Algorithm 5.1 Evaluation of the fibre damage.

Step 1 Spectral decomposition of the strain tensor (Equation 5.5 and 5.6).

Step 2 Loading of the variable at the previous time step.

Step 3 Evaluation of the thermodynamic forces (Equation 5.7).

Step 4 Calculation of the fibre damage rate (Equation 5.39).

Step 5 Evaluation of the fibre damage (Equations 5.30 and 5.31).

Step 6 Failure of the integration point if

d f i > κ ft i • h I f (i) • 1 2 C 0 ii • l ref 2 l 2 • ε ft,ref F(i) 2 -y ft 0(i)
or

d f i > κ fc i • h I f (i) • 1 2 C 0 ii • l ref 2 l 2 • ε fc,ref F(i)
2y fc 0(i) .

Step 7 Storage of the variables.

Step 8 Determination of the additional compliance due to the fibre damage (Equation 5.27).

Verification of the implementation by benchmark simulations 5.5.1 Mesh independence

In order to verify the independence to the mesh size of the energy dissipated by fibre failure on a given material volume, numerical simulations of tensile tests on 25 × 25 mm 2 initial area were carried out (Figure 5.11).

The responses given by the simulation of tensile tests along the fibre direction with a damage rate independent from the element size (Figure 5.12a) or with fibre damage regu- larisation (Figure 5.12b) show the difference between both fibre damage models. Whereas the dissipated energy with a constant damage rate tends quickly to the Γ ∞ , the dissipated energy in case of instantaneous fibre failure, the smeared crack model stabilises the energy dissipated during fibre failure It is important to clarify that the strain localisation is not avoided by the smeared crack method. The final fracture continue to be carried out by one element line transverse to the loading direction (Figure 5.13). Moreover, this localisation leads to spurious oscillations at failure because one element focuses all strains at failure.

Shear influence on the final fracture

After identification of the parameters c fx i to adapt the failure to shear loading, the procedure of validation adopted in the previous section 5. [START_REF] Allix | A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues[END_REF].1 is applied to in-plane shear loading. The resultants of the simulations for different element sizes are shown Figure 5.14.

In the case of shearing, the failure is more progressive because of sharing of the forces on the reinforcements. This behaviour has been experimentally observed and the simulation correlates with the expected behaviour.

However, the energy which is dissipated at failure seems random and very sensitive to the mesh size. In fact, the results are converging and this phenomenon can be explained by the observation of the crack obtained numerically (Figure 5.15). The failure is spread across the mesh and the opening is less clear than for the tensile tests. It leads to an important mesh dependency, besides the strain localisation.

Concluding remarks on the final fracture

In case of bi-axial fabric composites, the final fracture of the material is induced by the fibre fracture. The study of the failure of the fibres is based on two axis: the failure initiation and the post-failure behaviour. The failure initiation is predicted thanks to a modified maximal strain criterion. A spectral decomposition of the strain tensor is carried out in order to determine its positive part and its negative part. The components of the positive part in the fibre directions are used to predict the tensile fibre failure, whereas the components of the negative part in the fibre directions are used for the compressive failure. Hence, shear strain of the fabric ply influences the fibre failure.

The post-failure behaviour is managed through introduction of fibre damage variables. Due to the important loss of stiffness, a softening behaviour occurs. In case of implementation in a finite element code, the softening leads to a spurious mesh dependency because of a strain localisation on an element. The smeared crack approach is used to address this issue. It consists of managing the energy dissipation by using the constant fracture toughness, property of the material, weighted by the element size. Therefore for a constant volume of an element patch, the dissipated energy becomes constant.

In the future, the experimental evaluation of the fracture toughness will require specific attention. Other methods to control the fibre failure and to avoid a mesh dependence, such as a non-local formulation, can possibly be implemented. Other failure criterion can also be studied. Additional tests with loading according to a random angle with respect to the fibre directions, and bi-axial experiments are needed to further enhance the failure prediction.

CHAPTER 6

Laminate theory and modelling of the interlaminar damage

The mechanisms of deformation in a laminated fabric reinforced polymer are complex due to the heterogeneity of the stiffness through the thickness of the material. Hence, the usual Firstorder Shear Deformation Theory (FSDT), so-called Mindlin-Reissner, used as formulation of most shell elements in commercial finite element codes is not sufficient to describe the bending behaviour of the laminated structure. Moreover, due to the uniform transverse shear stress, the delamination, critical for the integrity of a composite structure, cannot be predicted. An Enhanced Higher-order Shear Deformation Theory (EHSDT) which is able to consider eventual interface imperfections is used as new displacement theory. This theory requires as many degrees of freedom as the FSDT. The determination of this realistic displacement field is introduced directly at the level of the material model. The bridge between the FSDT variables provided by the element formulation and the present theory is achieved by a strain energy equilibrium. By the knowledge of the true strain field, the stress at the integration points and the interface damage can then be evaluated. This theory, in association with the the intralaminar model for fabric reinforced polymers previously introduced in the previous chapters of the manuscript, is finally used to predict the behaviour of an industrial structure submitted to impact-type loading.

Introduction

Inside a laminated fibre reinforced polymer, a failure of the matrix between two adjacent plies leads to a massive degradation of the mechanical properties. This interface stiffness degradation, so-called delamination, is essential for the prediction of the mechanical behaviour of laminated structure, especially in case of impact loading.

The modelling of interface defects in case of laminated materials is nowadays still challenging for efficient numerical simulations in an industrial framework. Because the delamination is localised on the interfaces and due to the crack area widespread compared to the characteristic length of the material, the Continuum Damage Mechanics is not appropriate for the interlaminar damage modelling.

The delamination prediction is achieved by determining the stresses applied to the ply interfaces and by applying a failure criterion. The most common numerical method consists of the use of the mesoscopic scale and of a fine meshing of the out-of-plane dimension of the laminate. Once the stresses at ply interfaces, a cohesive law, with or without thickness, is used to describe the interlaminar behaviour. Therefore, it requires several elements -solid, shell or thick shell -through-the-thickness of the laminate using a ply constitutive law, coupled by inserting cohesive elements between two adjacent plies. The cohesive laws for fibre reinforced polymers have been, and still are, widely studied [6,[START_REF] Johnson | Computational methods for predicting impact damage in composite structures[END_REF][START_REF] Greve | Delamination testing and modelling for composite crash simulation[END_REF][START_REF] Pinho | Formulation and implementation of decohesion elements in an explicit finite element code[END_REF][START_REF] Bouvet | Low velocity impact modeling in composite laminates capturing permanent indentation[END_REF][START_REF] May | Numerical evaluation of cohesive zone models for modeling impact induced delamination in composite materials[END_REF][START_REF] Long | Delamination prediction in composite laminates under low-velocity impact[END_REF]. It has been proven that the delaminated area and the material responses are well predicted by using this method. However, the cohesive-based models suffer of two main drawbacks. As proven by Choi et al. [START_REF] Choi | A New Approach toward Understanding Damage Mechanisms and Mechanics of Laminated Composites Due to Low-Velocity Impact: Part II-Analysis[END_REF], the intralaminar and interlaminar matrix damage are coupled and is explained by the propagation of a transverse ply crack until the adjacent ply. The use of different elements dedicated to model either the ply behaviour or the interface behaviour makes difficult to couple both formulations. Moreover, for industrial applications, the design of the finite element model leaves the standards. A thin plate (such as a roof, a side door or a car hood in the automotive industry) is modelled by mean of shell elements and to design specific models with several elements through-the-thickness and cohesive elements may be tedious. In the interests of simplification and efficiency, discrete interface (springer-like) models can be used [START_REF] Allix | Mesomodeling of delamination: towards industrial applications[END_REF][START_REF] Ladevèze | A Computational Damage Micromodel of Laminated Composites[END_REF][START_REF] Bouvet | Low velocity impact modelling in laminate composite panels with discrete interface elements[END_REF].

Another solution is the formulation of specific finite elements dedicated to the simulation of laminated materials. An extended review of the different displacement-based formulations have been carried out by Khandan et al. [START_REF] Khandan | The development of laminated composite plate theories: a review[END_REF]. The classical laminate theory is based on the Kirchhoff-Love theory which assumes an homogeneous in-plane deformation and ignores the transverse shear strain. However composite plates are subject to transverse shear and is even a critical factor in the development of interlaminar damage. Mindlin-Reissner theory [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF][START_REF] Mindlin | Influence of rotary inertia and shear on flexural motions of isotropic elastic plates[END_REF] includes a constant transverse shear through the thickness which leads to constant cross-sectional rotations and a linear in-plane displacements variation. This theory is referred as the First-order Shear Deformation Theory (FSDT). It is nevertheless unable to give a physically acceptable distribution of transverse shear strain since the traction-free conditions on the top and bottom surfaces is not satisfied. For that, Reddy [START_REF] Reddy | A Simple Higher-Order Theory for Laminated Composite Plates[END_REF], Lo et al. [START_REF] Lo | A High-Order Theory of Plate Deformation-Part 2: Laminated Plates[END_REF] propose a third-order (so-called parabolic) shear deformation theory. But the stress continuity at the layer interfaces is not respected due to the lack of consideration of the strain discontinuities induced by the stiffness heterogeneity of laminated composite plates. Layerwise theories formulated by Reddy [START_REF] Reddy | A generalization of two-dimensional theories of laminated composite plates[END_REF] are able to represent the stress continuity at layer interfaces but require as much families of degrees of freedoms as the number of plies. These theories suffer of huge computational costs. Sciuva [START_REF] Sciuva | Bending, vibration and buckling of simply supported thick multilayered orthotropic plates: An evaluation of a new displacement model[END_REF], Sciuva and Icardi [START_REF] Sciuva | Analysis of thick multilayered anisotropic plates by a higher order plate element[END_REF] reduces the number of degree of freedoms for a zigzag theory as many as used by FSDT but loses boundary traction-free conditions. The composite plate theory given by Cho and Parmerter [START_REF] Cho | An efficient higher-order plate theory for laminated composites[END_REF][START_REF] Cho | Efficient higher order composite plate theory for general lamination configurations[END_REF] and called Enhanced Higher-order Shear Deformation Theory (EHSDT) solves the drawbacks of the previously mentioned theories. With the same number of degrees of freedom as the FSDT, the boundary traction-free conditions and the stress continuity at the layer interfaces are preserved. However, even if the in-plane stresses are very close to the exact solutions, the constitutive equations cannot predict an accurate transverse shear stress distribution. In order to accurately predict the interlaminar stresses, equilibrium equations have to be considered, either by using non-local method, or by using quadratic shell elements.

After further researches, Cho and Kim [START_REF] Cho | Higher-Order Zig-Zag Theory for Laminated Composites With Multiple Delaminations[END_REF] extend the EHSDT to take into account delamination by introducing interface imperfections. Despite the increasing complexity of the laminate theory, the number of degrees of freedom is preserved. Kim and Cho [START_REF] Kim | Efficient Higher-Order Shell Theory for Laminated Composites with Multiple Delaminations[END_REF], Oh et al. [START_REF] Oh | Dynamic analysis of composite plate with multiple delaminations based on higher-order zigzag theory[END_REF] used this displacement field to formulate new shell elements for finite element analysis of delaminated composite materials. Although efficient, these shell formulations are considering delamination as preliminary state of the material. The interlaminar damage is defined during the pre-processing phase and are not able to evolve during the simulation. Other formulations have been recently developed to introduce a real-time monitoring of delaminated area [102,[START_REF] Eijo | Delamination in laminated plates using the 4-noded quadrilateral QLRZ plate element based on the refined zigzag theory[END_REF]. However, a non-standard connection between the material model and the shell element formulation has to be established in order to transfer damage variables. The implementation in a commercial finite element software is complex and is outside a standard framework.

An efficient method proposed by Kim et al. [START_REF] Kim | Reconstruction of First-Order Shear Deformation Shell Theory for Laminated and Sandwich Shells[END_REF][START_REF] Kim | Enhanced modeling of laminated and sandwich plates via strain energy transformation[END_REF][START_REF] Kim | Efficient analysis of laminated composite and sandwich plates with interfacial imperfections[END_REF] provides a bridge between the First-Order Shear Deformation Theory (so-called Mindlin-Reissner) and the previously mentioned Enhanced Higher-order Shear Deformation Theory. Since the shell elements based on the Mindlin-Reissner are widely implemented in commercial finite element codes, this method takes advantages of both robustness of time-tested and enhanced shell formulations and accuracy of the EHSDT. The method proposed by Kim et al. consists of the calculation during a pre-processing phase of the shear correction factor representative of the bending behaviour of the investigated laminate. This factor is determined through transverse shear energy equivalence and least square approximation of strains and stresses between the FSDT and the EHSDT. By assuming that the FSDT is an averaging of a three-dimensional displacement field, and after numerical simulations with the appropriate shear correction factor for the Mindlin-Reissner shell elements, a realistic displacement field, and consequently the associated stress field, can be recovered during a post-processing phase. However the real-time prediction of delamination needed for impact simulations is not achieved with this modelling framework.

Thus, the present work concerns the implementation of this strain energy equivalence between FSDT and EHSDT directly in the user material subroutine. By using widespread Mindlin-Reissner shell elements of commercial finite element codes, it becomes possible to integrate delamination predictions directly at the level of the material model. Moreover, a complete coupling between the intralaminar and the interlaminar models become available. Therefore, the proposed formulation for the modelling of composite laminated structures is halfway between a material model and a finite element formulation, but is implemented only at the level of the user material subroutine which is easily transferable to the industrial sector.

The three-dimensional displacement and strain fields of an healthy laminated material are presented in a first section. Then, the interface imperfections are introduced in the previously described formulation by using damageable spring-layer model. As prevision of the coupling with a First-Order Shear Deformation Theory used in most of the shell elements in the commercial finite element codes, the average displacements of the Enhanced Higher-Order Shear Deformation Theory are provided and the strain energy equilibrium is established. The implementation scheme into the commercial finite element code LS-DYNA ® is described. In a next section, the identification procedure of the out-of plane and delamination parameters is given. Finally, the present model is validated through comparison between the predicted behaviour and the experimental response of an industrial layered structure made up with fabric reinforced polymer plies.

Internal displacement and strain fields of laminated structures

A thin composite plate is considered, consisting of a finite number P of stacked plies and with a thickness h. The out-of-plane position of the upper interface of the k-th ply is denoted by x 3(k) , whereas the bottom interface is denoted by x 3(k-1) (Figure 6.1).

1

(

1 -Ω) h 2 = x 3(P) Ref. plane 0 -1 -(1 + Ω) h 2 = x 3(0) Ω x 3 h x 3(1)
x 3(2)

x 3(P-1)

x 3(k) . . . . . . In order to formulate the most general displacement and strain field, the position of the reference plane Ω is introduced. By setting Ω coincident with the reference plane of the shell element section in the configuration of the finite element analysis, it becomes possible to define some change of sections in the material (Figure 6.2). To simplify the writing of the equations, Greek values take values of 1 or 2 and Latin indices take values of 1, 2 or 3. Moreover, the summation convention on repeated indices is used.

Ω

Ref. plane

Displacement field

The present theory relies on a third-order shear deformation theory [START_REF] Reddy | A generalization of two-dimensional theories of laminated composite plates[END_REF] as shown Figure 6.3 whose displacements are given by:

u α (x i ) = u 0 α (x β ) + ψ α (x β ) • x 3 + ξ α (x β ) • x 2 3 + φ α (x β ) • x 3 3 , (6.1)
u 3 (x i ) = u 0 3 (x β ) (6.2)
with u 0 i (x α ) the displacements of the reference plane at the position given by x α . The stiffness heterogeneity through the thickness of the layered material induces discontinuities in the displacement field at the ply interfaces. These discontinuities are inserted in the present formulation by a zigzag approach [START_REF] Cho | An efficient higher-order plate theory for laminated composites[END_REF][START_REF] Cho | Efficient higher order composite plate theory for general lamination configurations[END_REF] using slop angle S (k) α at the k-th interface level as follows:

u α (x i ) = u 0 α (x β ) + ψ α (x β ) • x 3 + ξ α (x β ) • x 2 3 + φ α (x β ) • x 3 3 + P-1 ∑ k=1 S (k) α (x β ) • x 3 -x 3(k) (6.3)
where < • > are the Macaulay brackets. Note that the out-of plane displacements remain the same. A representation of these slop angles are given Figure 6.4.

Traction-free conditions on the top and bottom surfaces

By applying the traction-free conditions on both top and bottom surfaces, it is possible to simplify the previous equations and to introduce the reference plane. For that purpose, the transverse shear angle need to be defined. In the present theory, the hypothesis of small transverse shear deformations is made. However, that does not prevent the use of the large strain framework for the in-plane strains. So, the transverse shear angle is given by:

γ α3 (x i ) = u 3,α + u α,3 = u 0 3,α (x β ) + ψ α (x β ) + 2 • ξ α (x β ) • x 3 + 3 • φ α (x β ) • x 2 3 + P-1 ∑ k=1 S (k) α (x β ) • H x 3 -x 3(k) (6.4)
with H(X) the Heaviside function equal to one if X > 0 and else null. Then, the traction-free conditions yields to:

x 3(P) = 0 ⇒ γ α3 (x i ) = u 0 3,α (x β ) + ψ α (x β ) + 2 • ξ α (x β ) • (1 -Ω) • h 2 +3 • φ α (x β ) • (1 -Ω) 2 • h 2 4 + P-1 ∑ k=1 S (k) α (x β ) = 0, , (6.5)
x 3(0) = 0 ⇒ γ α3 (x i ) = u 0 3,α (x β ) + ψ α (x β ) -2 • ξ α (x β ) • (1 + Ω) • h 2 +3 • φ α (x β ) • (1 + Ω) 2 • h 2 4 = 0. . ( 6.6)
By subtracting Equations 6.5 and 6.6, the parameter ξ α (x β ) is defined as follows:

ξ α (x β ) = 1 2 • h 3 • φ α (x β ) • h 2 • Ω - P-1 ∑ k=1 S (k) α (x β ) , (6.7)
and after replacing it in Equation 6.6 the following condition on the displacement of the reference plane is obtained:

u 0 3,α (x β ) + ψ α (x β ) = 3 • h 2 4 • Ω 2 -1 • φ α (x β ) - 1 + Ω 2 • P-1 ∑ k=1 S (k) α (x β ). (6.8)
Finally by replacing the terms ξ α (x β ) and u 0 3,α (x β ) + ψ α (x β ) in Equation 6.4, the transverse shear angle become:

γ α3 (x i ) = 3 • h 2 4 • Ω 2 -1 + h • Ω • x 3 + x 2 3 • φ α (x β ) + P-1 ∑ k=1 S (k) α (x β ) • - 1 + Ω 2 - x 3 h + H x 3 -x 3(k) .
(6.9)

Continuity conditions of the transverse shear stresses

By applying the conditions of continuity of the transverse shear stresses at the layer interfaces, the slop angles will be expressed according to material constants and the parameter φ α (x β ).

For the needs of writing, the stiffness matrix of the plies expressed in the global coordinate system are split into in-plane and out-of-plane parts. Thus, the in-plane stiffness tensor Q (m) and the out-of plane stiffness tensor Q s(m) of the m-th ply in the global coordinate system are given by:

Q (m) =    Q (m) 11 Q (m) 12 0 Q (m) 12 Q (m) 22 0 0 0 Q (m) 66    , Q s(m) = Q (m) 44
Q (m) 45
Q (m) 45
Q (m) 55
. (6.10)

At the m-th interface located at the out-of-plane coordinate x 3(m) , the transverse shear stresses infinitely close to the interface and located either in the upper ply or in the lower ply are admitted to be equal (continuity conditions of the transverse shear stresses):

σ 13 (x α , x 3 = x + 3(m) ) = σ 13 (x α , x 3 = x - 3(m) ) σ 23 (x α , x 3 = x + 3(m) ) = σ 23 (x α , x 3 = x - 3(m) ) (6.11)
By using the stress/strain constitutive equation, it yields to:

Q (m+1) 44 • γ 23 (x α , x 3 = x + 3(m) ) + Q (m+1) 45 • γ 13 (x α , x 3 = x + 3(m) ) = Q (m) 44 • γ 23 (x α , x 3 = x - 3(m) ) + Q (m) 45 • γ 13 (x α , x 3 = x - 3(m) )
, (6.12)

Q (m+1) 45 • γ 23 (x α , x 3 = x + 3(m) ) + Q (m+1) 55 • γ 13 (x α , x 3 = x + 3(m) ) = Q (m) 45 • γ 23 (x α , x 3 = x - 3(m) ) + Q (m) 55 • γ 13 (x α , x 3 = x - 3(m) ) . (6.13)
It is convenient for the next step to express literally the transverse strains close to the interfaces. Thus, they are given by:

         γ α3 (x β , x 3 = x + 3(m) ) = 3 • h 2 4 • Ω 2 -1 + h • Ω • x 3(m) + x 2 3(m) • φ α (x β ) + P-1 ∑ k=1 S (k) α (x β ) • - 1 + Ω 2 - x 3(m) h + m ∑ k=1 S (k) α (x β ) (6.14)          γ α3 (x β , x 3 = x - 3(m) ) = 3 • h 2 4 • Ω 2 -1 + h • Ω • x 3(m) + x 2 3(m) • φ α (x β ) + P-1 ∑ k=1 S (k) α (x β ) • - 1 + Ω 2 - x 3(m) h + m-1 ∑ k=1 S (k) α (x β ) (6.15)
Then by introducing the transverse shear stiffness jump at the m-th interface ∆Q (m)

i j = Q (m+1) i j -Q (m)
i j and by using Equation 6.9 to express the transverse shear strains, these relations become:

                                                             m-1 ∑ k=1 S (k) 2 (x α ) • ∆Q (m) 44 • - 1 + Ω 2 - x 3(m) h + 1 +S (m) 2 (x α ) • Q (m+1) 44 + ∆Q (m) 44 • - 1 + Ω 2 - x 3(m) h + P-1 ∑ k=m+1 S (k) 2 (x α ) • ∆Q (m) 44 • - 1 + Ω 2 - x 3(m) h + m-1 ∑ k=1 S (k) 1 (x α ) • ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h + 1 +S (m) 1 (x α ) • Q (m+1) 45 + ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h + P-1 ∑ k=m+1 S (k) 1 (x α ) • ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h = -3 • ∆Q (m) 44 • h 2 4 • Ω 2 -1 + h • Ω • x 3(m) + x 2 3(m) • φ 2 (x α ) -3 • ∆Q (m) 45 • h 2 4 • Ω 2 -1 + h • Ω • x 3(m) + x 2 3(m) • φ 1 (x α ) , (6.16)
                                                             m-1 ∑ k=1 S (k) 2 (x α ) • ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h + 1 +S (m) 2 (x α ) • Q (m+1) 45 + ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h + P-1 ∑ k=m+1 S (k) 2 (x α ) • ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h + m-1 ∑ k=1 S (k) 1 (x α ) • ∆Q (m) 55 • - 1 + Ω 2 - x 3(m) h + 1 +S (m) 1 (x α ) • Q (m+1) 55 + ∆Q (m) 55 • - 1 + Ω 2 - x 3(m) h + P-1 ∑ k=m+1 S (k) 1 (x α ) • ∆Q (m) 55 • - 1 + Ω 2 - x 3(m) h = -3 • ∆Q (m) 45 • h 2 4 • Ω 2 -1 + h • Ω • x 3(m) + x 2 3(m) • φ 2 (x α ) -3 • ∆Q (m) 55 • h 2 4 • Ω 2 -1 + h • Ω • x 3(m) + x 2 3(m) • φ 1 (x α )
. (6.17)

For further considerations these relations are presented in a matrix form as follows:

A • S(x α ) = B 1 • φ 1 (x α ) + B 2 • φ 2 (x α ) (6.18) with A = (A i, j ) ∈ R (2P-2)×(2P-2) (6.19)
where the even lines of A are given by:

                                       A 2m-1,1→m-1 = ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h + 1 A 2m-1,m = Q (m+1) 45 + ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h A 2m-1,m+1→P-1 = ∆Q (m) 45 • - 1 + Ω 2 - x 3(m) h A 2m-1,P→P+m-2 = ∆Q (m) 44 • - 1 + Ω 2 - x 3(m) h + 1 A 2m-1,P+m-1 = Q (m+1) 45 + ∆Q (m) 44 • - 1 + Ω 2 - x 3(m) h A 2m-1,P+m→2P-2 = ∆Q (m) 44 • - 1 + Ω 2 - x 3(m) h , (6.20)
u (d) 3 (x γ ) = R (d) 33 (x γ ) • σ (d) 33 (x γ , x 3(d) ), (6.34)
with R In this work, mode I opening are not taken into account for the delamination. As a consequence, R (d) 33 (x γ ) = 0. Under this condition and by using Equation 6.32, the transverse shear angle is given by:

γ α3 (x i ) = u 3,α (x i ) + u α,3 (x i ) = Φ αβ ,3 (x 3 ) • φ β (x γ ) (6.35)
where

Φ αβ ,3 (x 3 ) = 3 • h 2 4 • (Ω 2 -1) + h • Ω • x 3 + x 2 3 • δ αβ + P-1 ∑ k=1 a (k) αβ • - 1 + Ω 2 - x 3 h + H x 3 -x 3(k) .
(6.36)

Since the transverse shear stress is continuous across the thickness and is given in the k-th ply by:

σ (k) α3 (x i ) = Q s(k) αβ • γ β 3 (x i ) (6.37)
and by substituting Equations 6.35 and 6.33 into Equation 6.32, the in-plane displacement according to the out-of-plane position is now given by:

u α (x i ) = u 0 α (x λ ) -u 0 3,α (x λ ) • x 3 + Φ αβ (x 3 ) • φ β (x λ ) + P-1 ∑ d=1 R (d) αγ • Q s(d) γω • Φ ωβ ,3 (x 3(d) ) • φ β (x λ ) • H x 3 -x 3(d) .
(6.38)

Finally and following the factoring of φ β (x γ ), the displacement field is given in a simple form by:

u α (x i ) = u 0 α (x γ ) -u 0 3,α (x γ ) • x 3 + Ψ αβ (x 3 ) • φ β (x γ ) (6.39)
with

Ψ αβ (x 3 ) = Φ αβ (x 3 ) + T αβ (x 3 ) (6.40)
where T αβ (x 3 ) is function of the interface damages and is given by:

T αβ (x 3 ) = P-1 ∑ d=1 R (d) αγ • Q s(d) γω • Φ ωβ ,3 (x 3(d) ) • H x 3 -x 3(d) . (6.41)
However, by this formulation and because of the zigzag function and the eventual delamination, the variable u 0 α does not necessary characterise the displacement of the reference plane. Therefore, an integral form C w αβ is introduced to convert the displacement component of the reference plane u 0 α into a variable, such as:

Ψ αβ (x 3 ) = Φ αβ (x 3 ) + T αβ (x 3 ) -C w αβ (6.42)
with

C w αβ = 1 h • x 3(P)
x 3(0)

Φ αβ (x 3 ) + T αβ (x 3 ) • dx 3 . (6.43)
By means of this equation, the laminate displacement field with possibility of delamination is fully described and written in a simple manner. It is interesting to note that Equation 6.39 is similar to the Mindlin-Reissner theory. The only change is Ψ αβ which depends on the out-of-plane position.

6.4 Strain energy equilibrium 6.4.1 Averaged displacements and strain fields in least square sense

In order to use the present theory in combination with finite elements using a formulation based on the Mindlin-Reissner theory, both have to present the same amount of strain energy. Kim [START_REF] Kim | Reconstruction of First-Order Shear Deformation Shell Theory for Laminated and Sandwich Shells[END_REF], Kim and Cho [START_REF] Kim | Enhanced modeling of laminated and sandwich plates via strain energy transformation[END_REF], Kim et al. [START_REF] Kim | Efficient analysis of laminated composite and sandwich plates with interfacial imperfections[END_REF] proposed a methodology to achieve this balance by considering the First-order Shear Deformation Theory (FSDT) as the average displacement in the least square sense of the Enhanced Higher-order Shear Deformation Theory (EHSDT). Therefore, the average displacements are given by the Mindlin-Reissner theory, such as:

v α (x i ) = v 0 α (x β ) + x 3 • θ α (x β ), (6.44)
v 3 (x i ) = v 0 3 (x β ). (6.45)
By applying the least square approximation, the displacement variables of the FSDT can be expressed according to the displacement variables of the EHSDT. Thus, the relations between the variables are given by: min

v 0 α x 3(P)
x 3(0)

∥u α -v α ∥ 2 • dx 3 → v 0 α = u 0 α + x 3(P)
x 3(0)

Ψ αβ • dx 3 • φ β , (6.46)
min θ α

x 3(P)

x 3(0)

∥u α -v α ∥ 2 • dx 3 → θ α = -u 0 3,α + 12 h 3 • x 3(P)
x 3(0)

x 3 • Ψ αβ • dx 3 • φ β , (6.47)
min v 0 3 x 3(P)
x 3(0)

∥u 3 -v 3 ∥ 2 • dx 3 → v 0 3 = u 0 3 . (6.48)
The relation between the engineering strains of both theories can also be obtained and are given in a simple form by Kim et al. [START_REF] Kim | Efficient analysis of laminated composite and sandwich plates with interfacial imperfections[END_REF]. Note that when the balance between both strain energies is achieved for the engineering strains, it is also the case for strains in the finite strain framework (and notably Green-Lagrange strains). Thus, the engineering strains given by the FSDT, are expressed in a vector form by:

⃗ γ e = ⃗ γ 0 + x 3 • ⃗ λ , ⃗ γ s = ⃗ γ 0 s (6.49)
where ⃗ γ e represents the in-plane strain and ⃗ γ s represents the transverse shear strain. The in-plane shear of the reference plane ⃗ γ 0 , the curvature ⃗ λ and the transverse shear strain ⃗ γ 0 s vectors are respectively given by:

⃗ γ 0 = v 0 1,1 v 0 2,2 (v 0 1,2 + v 0 2,1 ) T , (6.50)
⃗ λ = θ 0 1,1 θ 0 2,2 (θ 0 1,2 + θ 0 2,1 ) T , (6.51)
⃗ γ 0 s = θ 1 + v 0 3,1 θ 2 + v 0 3,2 T . (6.52)
In the same manner, the engineering strains given by the EHSDT are expressed in a tensor form by:

⃗ ε e =⃗ ε 0 + x 3 •⃗ κ + Ψ Ψ Ψe •⃗ κ h , ⃗ ε s = Ψ Ψ Ψs •⃗ ε 0 s (6.53)
where

⃗ ε 0 = u 0 1,1 u 0 2,2 (u 0 1,2 + u 0 2,1 ) T , (6.54)
⃗ κ = -u 0 3,11 u 0 3,22 2 • u 0 3,12 T , (6.55)
⃗ κ h = φ 1,1 φ 2,2 φ 1,2 φ 2,1 T , (6.56)
Ψ Ψ Ψe =    Ψ 11 0 0 Ψ 12 0 Ψ 22 Ψ 21 0 Ψ 21 Ψ 12 Ψ 11 Ψ 22    , (6.57)
⃗ ε 0 s = φ 1 φ 2 T , (6.58)
Ψ Ψ Ψs = Ψ 11,3 Ψ 12,3 Ψ 21,3 Ψ 22,3 . (6.59)
Finally the strain relationships between both theories are given by:

⃗ ε 0 = ⃗ γ 0 -Ξ Ξ Ξ •⃗ κ h , (6.60)
⃗ κ 0 = ⃗ λ -Γ Γ Γ •⃗ κ h , (6.61)
⃗ ε 0 s = Γ Γ Γ-1 s •⃗ γ 0 s (6.62) with Ξ Ξ Ξ = x 3(P) x 3(0) Ψ Ψ Ψe • dx 3 , (6.63)
Γ Γ Γ = 12 h 3 •

x 3(P)

x 3(0)

x 3 • Ψ Ψ Ψe • dx 3 , (6.64)
Γ Γ Γs = 12 h 3 •Γ Γ Γ (6.65)
where the components of Γ Γ Γ are defined as follows:

Γ αβ = x 3(P)
x 3(0)

x 3 • Ψ αβ • dx 3 . (6.66)

Strain energy transformation

The strain energy which results of a displacement theory is deduced by means of the laminate constitutive equation. Therefore, in case of the FSDT, this laminate constitutive equation is expressed in a tensor form by:

⃗ N ⃗ M = A B B T D ⃗ γ 0 ⃗ λ , ⃗ V = Ĝ •⃗ γ 0 s . (6.67)
The tensors representative of the tension stiffness A, the tension-bending coupling stiffness B, the bending stiffness D and the transverse shear stiffness Ĝ are obtained by integrating the stiffness on the thickness of the laminate, such as:

A = x 3(P)
x 3(0)

Q • dx 3 , B = x 3(P)
x 3(0)

x 3 • Q • dx 3 , D = x 3(P)
x 3(0)

x 2 3 • Q • dx 3 , (6.68)
and Ĝ =

x 3(P)

x 3(0)

Q s • dx 3 . (6.69)
It yields to the following expression for the strain energy of the FSDT:

2 • Û = ⃗ γ 0T • A •⃗ γ 0 + 2 •⃗ γ 0T • B • ⃗ λ + ⃗ λ T • D • ⃗ λ +⃗ γ 0 s T • Ĝ •⃗ γ 0 s . (6.70)
Regarding the EHSDT, the laminate constitutive equation becomes:

     ⃗ N ⃗ M ⃗ R      =    A B E B T D F E T F T H         ⃗ ε 0 ⃗ κ ⃗ κ h      , ⃗ V = G •⃗ ε 0 s , (6.71)
with E =

x 3(P)

x 3(0)

Q • Ψ Ψ Ψe • dx 3 , F = x 3(P)
x 3(0)

x 3 • Q • Ψ Ψ Ψe • dx 3 , H = x 3(P) x 3(0) Ψ Ψ ΨT e • Q • Ψ Ψ Ψe • dx 3 , (6.72)
and G =

x 3(P)

x 3(0) Ψ Ψ ΨT s • Q s • Ψ Ψ Ψs • dx 3 . (6.73)
The stiffness tensors E, F and H are representative of the higher-order zigzag effect through the thickness of the laminate. The tensors A, B and D remain the same as used in the FSDT. By contrast, the transverse shear stiffness varies from one theory to the other. Usually, the transverse shear energy of the FSDT is corrected by a shear correction factor. In the case of the EHSDT, the transverse strain energy is realistic of the real distribution and do not need any correction. Consequently, the strain energy of the EHSDT is given by:

2 •U = ⃗ ε 0T • A •⃗ ε 0 + 2 •⃗ ε 0T • B •⃗ κ +⃗ κ T • D •⃗ κ +⃗ ε 0 s T • Ĝ •⃗ ε 0 s +2 •⃗ ε 0T • E •⃗ κ h + 2 •⃗ κ T • F •⃗ κ h + •⃗ κ T h • H •⃗ κ h . (6.74)
Then by using Equations 6.60, 6.61 and 6.62 and by replacing the deformation variables of the EHSDT by deformation variables of the FSDT, the strain energy of the EHSDT becomes:

2 •U = 2 • Û + 2 • Ũ (6.75)
with Û the strain energy of the FSDT (Equation 6.70) and where Ũ is given by:

2 • Ũ = 2 •⃗ ε 0T • -A • Ξ Ξ Ξ -B • Γ Γ Γ + E •⃗ κ h + 2 •⃗ κ T • -B T • Ξ Ξ Ξ -D • Γ Γ Γ + F •⃗ κ h +2 •⃗ κ T h • Ξ Ξ ΞT • A • Ξ Ξ Ξ + 2 • Ξ Ξ ΞT • B • Γ Γ Γ -2 • Ξ Ξ ΞT • E + Γ Γ ΓT • D • Γ Γ Γ -2 • Γ Γ ΓT • F + H •⃗ κ h . (6.76)
Therefore, by determining appropriate tensors Ξ Ξ Ξ and Γ Γ Γ in order to make Ũ as close to zero as possible, the strain energy of the EHSDT will be equivalent to the strain energy of the FSDT, provided through the Mindlin-Reissner variables by the finite element. By this energy balance, it is now possible to realise a coupling between the element formulation and the appropriate strain field for the laminate inside the material model. The algorithm used to determine both tensors Ξ Ξ Ξ and Γ Γ Γ will be provided Section 6.6. Now that Ξ Ξ Ξ and Γ Γ Γ are determined to satisfy the strain energy equilibrium, the displacement field can be recovered by the relation:

u α (x i ) = v 0 α (x γ ) -v 0 3,α (x γ ) • x 3 + (Ψ αβ (x 3 ) -Ξ αβ ) • φ β (x γ ) (6.77)
where

φ α (x γ ) = Γ -1 αβ • (θ β + v 0 3,β ) = Γ -1 αβ • γ 0 sβ . (6.78)
and by considering:

Ξ Ξ Ξ =    Ξ 11 0 0 Ξ 12 0 Ξ 22 Ξ 21 0 Ξ 21 Ξ 12 Ξ 11 Ξ 22    , (6.79) and Γ
Γ Γ =    Γ 11 0 0 Γ 12 0 Γ 22 Γ 21 0 Γ 21 Γ 12 Γ 11 Γ 22    . (6.80)
And finally the deformation gradient, used for a finite strain formulation, is obtained by spatial derivation of the displacement field and is given by:

u α,β (x i ) = v 0 α,β (x µ ) -v 0 3,αβ (x µ ) • x 3 + (Ψ αγ (x 3 ) -Ξ αγ ) • Γ -1 γω • γ 0 sω,β (x µ ), (6.81)
u α,3 (x i ) = -v 0 3,α (x µ ) + Ψ αγ,3 (x 3 ) • Γ -1 γω • γ 0 sω (x µ ), (6.82)
u 3,α (x i ) = v 0 3,α (x µ ). (6.83)
It is essential to note that the transverse shear strain distribution obtained from the constitutive relation is not exactly representative of the exact transverse shear strains through the thickness of the laminate. Unfortunately, an equilibrium equation is needed to recover a distribution very close to the exact solution. It is given in the three-dimensional case by the equation:

σ i j, j = 0. (6.84)
It can be achieved only by the knowledge of the spatial derivatives of the in-plane stress field, not available in the present model. For that, a non-local method or quadratic shell elements could be used in further studies.

Delamination criterion

In the case of this work, the modelling of the interface layer have not been widely studied.

Hence, a simple delamination criterion associated to the damage evolution have been used.

The delamination criterion is defined by a quadratic form on the maximal interlaminar shear stresses. Thus, it is defined at the k-th interface in the global coordinate system by: , the ultimate transverse shear stresses of the k-th interface, are parameters of the model.

µ (k) = σ (k) zx σ max(k) zx 2 + σ (k) zy σ max(k) zy 2 -1 (6.
Then the compliance of the string layer interface is dependant of the delamination criterion and is given by:

R (k) = R c(k) • 1 -exp - < µ (k) > τ (k) (6.86)
with R c(k) a model parameter which defines the critical compliance and τ (k) a model parameter which influences the rate of interlaminar damage.

The compliance of the interface layer is considered independent of the loading direction.

Implementation

This part of the material model acts as an interface between the element formulation and the constitutive equation of the material model. Because the strain equilibrium between both element formulation and EHSDT is ensures by this theory, the stresses at the integration points calculated from the laminate strain theory can be keep as such and can be given as output without any adaptation. As a consequence, the post-processing of the finite element analysis is easier. So, this module acts directly on the input parameters of the user material subroutine to provide the laminate strain field to the stress update function.

From the element formulation, the deformation gradient and the spatial derivative of the deformation gradient on a point of the reference plane are obtained. This step depends on the finite element code used to implement the material model.

The procedure used to compute the strain field at the increment n and according to variables expressed at the previous increment n -1 is given in Algorithm 6.1. For the first increment, the stiffness are initialised thanks to the elastic stiffness provided in the material card. For the integration on the thickness of the laminate, a Simpson's rule is used.

In order to achieve Step 6 of Algorithm 6.1 the difference between both strain energies provided by the FSDT and the EHSDT and given Equation 6.76 should reach a negligible value. Kim [START_REF] Kim | Reconstruction of First-Order Shear Deformation Shell Theory for Laminated and Sandwich Shells[END_REF], Kim and Cho [START_REF] Kim | Enhanced modeling of laminated and sandwich plates via strain energy transformation[END_REF] proposed a solving procedure by using the two following equations:

-A • Ξ Ξ Ξ -B • Γ Γ Γ + E = 0, -B T • Ξ Ξ Ξ -D • Γ Γ Γ + F = 0 (6.87)
which are not directly solvable. To minimise the errors for the in-plane strains, Ξ Ξ Ξ has to be as close to null tensor as possible. Consequently, the iterative tensor forms, where (i) indicates the i-th increment, of the relationship tensors Ξ Ξ Ξ and Γ Γ Γ are given by:

Γ Γ Γ = D -1 • D (i-1) -B • Ξ Ξ Ξ(i-1) , (6.88)
Ξ Ξ Ξ(i) = A -1 • E (i-1) -B • Γ Γ Γ , (6.89)
and in-plane warping function is given by:

Ψ (i) αβ (x 3 ) = Φ αβ (x 3 ) + T αβ (x 3 ) -C w αβ + Ξ (i) αβ . (6.90)
In order to minimise the in-plane normal strain rather the transverse shear strain, Ξ

αβ is defined as follows:

Ξ (i) 11 = Ξ(i) 11 , Ξ (i) 12 = Ξ(i) 14 , Ξ (i) 21 = Ξ(i) 23 , Ξ (i) 22 = Ξ(i) 22 . (6.91)
Then the iterative procedure is given in Algorithm 6.2.

Algorithm 6.1 Determination of the laminate strain field from Mindlin-Reissner variables.

Step 1 Determination of the transverse shear stiffness tensors Q s(m) in the element coordinate system by using those of the previous increment given by the stress update function.

Step 2 Determination of the in-plane stiffness tensors Q (m) in the element coordinate system by using those of the previous increment given by the stress update function.

Step 3 Calculation of the A,B and D tensors (Equation 6.68).

Step 4 Calculation of the slope angle coefficients (Equation 6.25).

Step 5 Calculation of C w (Equation 6.43) by using compliance of the interface layer computed at the previous increment.

Step 6 Calculation of Ξ Ξ Ξ and Γ Γ Γ (Algorithm 6.2) by using compliance of the interface layer computed at the previous increment.

Step 7 Calculation of the deformation gradient at layer interfaces (Equations 6.81 to 6.83).

Step 8 Calculation of the strains at layer interfaces.

Step 9 Calculation of the stresses at layer interfaces.

Step 10 Evaluation of the interlaminar damage (Equations 6.85 and 6.86).

Step 11 Storage of the variables.

Step 12 Calculation of the deformation gradient at current integration point.

Step 13 Calculation of the strains at current integration point.

Step 14 Calculation of the strains at current integration point in the material coordinate system.

Algorithm 6.2 Iterative procedure to achieve a strain energy balance between both FSDT and EHSDT theories.

Step 1 Initialisation of Ξ Ξ Ξ(0) = 3 × 4null matrix.

Step 2 Calculation of Ψ Ψ Ψ(i-1) e (x 3 ) (Equations 6.90 and 6.57).

Step 3 Calculation of the E (i-1) and F (i-1) tensors (Equation 6.72).

Step 4 Calculation of Γ Γ Γ and Ξ Ξ Ξ (Equations 6.88 and 6.89).

Step 5 Update Ξ Ξ Ξ(i) using Ξ Ξ Ξ(i+1) = Ξ Ξ Ξ(i) + Ξ Ξ Ξ(i-1) .

Step 6 If Ξ Ξ Ξ(i) -Ξ Ξ Ξ(i-1) Ξ Ξ Ξ(i) < tolerance then continue, else go to step 2.

Step 7 If Γ11 Γ22 ≥ Γ33 Γ34 , then

Γ 11 = Γ11 , Γ 22 = Γ22 , Γ 12 = Γ14 , Γ 21 = Γ23 else Γ 11 = Γ33 , Γ 22 = Γ34 , Γ 12 = Γ14 , Γ 21 = Γ23 . Step 8 Calculation of Ξ Ξ Ξ = A -1 • E -B • Γ Γ Γ .

Identification of the out-of-plane shear properties

The parameter identification of the out-of-plane properties and the interlaminar damage model is done by means of a bending experimental test carried out independently on the different preforms. These tests are in conformity with the standardised interlaminar shear tests NF EN ISO 14130 in order to be easily reproducible in an industrial framework. It consists of a three-point bending where the supports are spaced of only 10 mm to favour the transverse shear strain inside the laminate specimen (Figure 6.6). These tests are carried out at room temperature on an electromagnetic device (Sintech 20D) which ensures 100 kN for the maximal load capacity and at a speed of 1 mm min -1 . The specimens are cut using the water-jet technique and their shape and dimensions, are shown in Figures 6.7 By considering the in plane behaviour fully identified, the first step consists of the determination of the transverse shear stiffness, namely G 13 and G 23 , through interlaminar shear tests in the fibre directions (Figure 6.9). During elastic deformation, the transverse shear stiffness is the only parameter which have an influence on the bending behaviour 20 mm 30 mm Figure 6.8: Specimen geometry used for interlaminar shear tests of 12K woven preforms.

⃗ X = ⃗ 1 ⃗ Y = ⃗ 2 (a) Longitudinal interlaminar shear test. ⃗ X ⃗ Y = ⃗ 1 ⃗ 2 90°( b) Transverse interlaminar shear test. ⃗ X ⃗ Y ⃗ 1 ⃗ 2 45°( c)
In-plane/interlaminar shear test. Figure 6.9: Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for interlaminar shear tests. and which have not been identified. Thus, G 13 and G 23 are obtained through an inverse engineering procedure.

According to the standard NF EN ISO 14130, a first approximation of the delamination criteria σ (k) zα is obtained, regardless the loading direction, as follows:

σ max zα = 3 4 • F max b • h (6.92)
where F max is the maximal force applied to the central pinch, b is the width of the specimen (10 mm) and h is the thickness of the composite.

The bending direction influences a lot the damage of the preforms. In case of bending at 45°according to the fibre direction, which leads to a coupling between in-plane/transverse shear strain of the preform, none intralaminar failure occurs. The only damage visible is interlaminar failure as shown Figures 6.10b and 6.11b. Moreover, the delamination is progressive and the evolution of pure interlaminar damage (Equation 6.86) can be well identified. Hence, this test is ideal for the identification of the delamination criterion and the evolution of the compliance of the interface layer. The results of this identification procedure are provided in Table 6.1.

Paramaters

NCF 3K woven 12K woven

σ max(k) zx = σ max(k) zy 75 58 45 R c(k) 0.4 0.1 0.4 τ (k) 5 5 5
Otherwise, in case of bending in the fibre direction, a massive transverse crack occurs through the thickness of the plies (Figures 6.10a and 6.11a). This transverse ply failure, due to the transverse shear (easily recognisable by the failure angle) leads to a sudden and important loss of bending stiffness. In addition, it has been proved that an intralaminar damage will initiate a delamination [START_REF] Choi | A New Approach toward Understanding Damage Mechanisms and Mechanics of Laminated Composites Due to Low-Velocity Impact: Part II-Analysis[END_REF][START_REF] Kashtalyan | The effect of delaminations induced by transverse cracks and splits on stiffness properties of composite laminates[END_REF]. From that, this test is not adapted to the characterisation of the interface layer. However, it provides useful information about the intralaminar damage and notably the effect of the transverse shear strains, but also about the intra/interlaminar damage coupling.

In this work, the intra/interlaminar damage coupling has not been studied even if the user material subroutine is ready to implement it thanks to damage variable transfer between intraand interlaminar model. Moreover, the intralaminar fracture direction due to transverse shear stresses is not considered as privileged direction of damage and the influence of this failure can only be approximated by the three damage variables already implemented.

The simulations of the different interlaminar shear tests were done by means of the present material model implemented in FORTRAN 90 in a user material subroutine as described section 4.5 for the commercial explicit finite element code LS-DYNA ® . As shown Figures 6.12 and 6.14, the numerical responses of the interlaminar shear tests in the fibre direction are close to the experimental responses. However, since only the delamination is taken into account, the loss of stiffness at delamination is not sufficient. The transverse matrix cracking and the coupling with the delamination have to be considered to well describe the bending behaviour. The numerical simulations of the interlaminar shear tests at 45°shows a better correlation with the experimental responses (Figures 6.13 and 6.17

Evaluation of the prediction ability on an industrial structure

The bi-axial non-crimp, the 3K woven and the 12K woven fabrics constitute the various preforms of a plate used as non-structural component in the automotive industry (Section 2.1.3).

The stacking sequence of this laminate is given by (0°) 3K , (±45°) NCF , (0°) 12K as shown Figure 6. [START_REF] Bouvet | Low velocity impact modeling in composite laminates capturing permanent indentation[END_REF] defined as the 0°direction of the roof structure. As a consequence, the 45°direction is defined by the following stacking sequence: (45°) 3K , (0°) NCF , (90°) NCF , (45°) 12K .

In order to validate the material model, two types of experimental tests have been carried out on the investigated structure: the standard interlaminar shear test and a dynamic bulge test in order to approach impact-type loading. Then, and by using exclusively the parameters provided by the identification procedure on the specimens made up with single preforms, numerical simulations have been used to predict the material behaviour. The evaluation of the model is carried out through comparison of the experimental and numerical responses, as well as comparison of the damaged area. Experimentally, the damage visualisation is carried out by means of C-Scan with a 3-GHz sensor, enabling a resolution of 100 µm. These ultrasonic inspections have been carried out by the IEMN-DOAE (UMR CNRS 8520) laboratory at the University of Valenciennes and Hainaut-Cambrésis.

Regarding the transfer of the parameter identified for the single preforms to the structure material, two remarks have to be done. Firstly, the fibre volume ratios between the single preforms and the structure are not strictly the same (please refer to Table 2.6). Consequently, the choice was done to retain the thickness of the plies, such as in single preforms. Thus, the thickness of the structure for the simulations is given by 0.25 + 2 × 0.4 + 0.3 = 1.35 mm against 1.40 mm for the real structure plate. Also, the properties of the interface layer are not strictly the same since the coincident preforms differ. In the simulations of the structure plate, the properties of the interfaces are taken as the mean values of the properties provided by the adjacent preforms.

Standardised interlaminar shear tests

As first validation step, interlaminar shear tests in both 0°and 45°according to the fibre directions are carried out and simulated without parameter adaptation. The experimental procedure of these tests have been provided Section 6.7 and remains the same. The coupon size is the same as for the 3K woven and NCF specimens (Figure 6.7).

The finite element model is presented Figure 6.19. The shell element formulation which has been used is the Belytschko-Tsay with reduced integration. The shell elements defining the specimen part have a size of 1 mm used in the industry. In regards to the results, the bending stiffness and the emergence of interlaminar damage are relatively well predicted. However the simulations lack to predict the post-failure behaviour, since the transverse shear matrix fracture is not modelled. Moreover, in case of the interlaminar shear test at 0°, the bending induces a premature fibre failure of the 12K ply. The thinness of this woven ply and the small gap between the yarns may lead to an excessive stress concentration and the fibre failure. This fibre failure can be seen in yellow on the ultrasonic inspections of the specimens after testing in Figure 6.29. Also, the influence of the transverse shearing on the fibre failure is not taken into account in the present fibre failure criterion.

The variety of preforms used in the structure induced a stiffness heterogeneity inside the laminate which influences the strain distribution, and consequently the stress through the thickness of the structure. Figures 6.22 and 6.23 provide an example of stress distribution before and after delamination for 0°structure specimen exposed to bending. On the other hand, Figures 6.24 and 6.25 provide an example of stress distribution before and after delamination for 45°structure specimen. These figures are showing the stress values at integration points, when two integration points per ply are used for the section of the shell. The emergence of sliding at the interface located in the middle of the specimen is clearly spotted on the figures. Compared to the damaged area predicted by finite element analysis (Figures 6. [START_REF] Chaboche | A new micromechanics based CDM model and its application to CMC's[END_REF] ures 6.29 and 6.30) is located on one interface layer and on one side of the specimen. This is due to the combination of two phenomena. The first one is the approximation of the transverse shear strain field by the actual theory. The second one is the use of a deterministic approach for the numerical simulations. The presence of eventual defects is ignored and leads to an "ideal" behaviour. The introduction of fuzzy variables could lead to a more realistic, but less predictable, behaviour. Otherwise, the progressive and diffuse delamination in the 0°c oupon due to plies in the middle oriented at 45°according to the bending and the sudden failures in the 45°coupon, located at different interfaces, are well reported. Otherwise, and as given by the ultrasonic inspections, no damage is present before failure. Interlaminar damage Figure 6.30: Ultrasonic inspections for damage visualisation on a structure specimen after interlaminar shear tests according to the 45°structure direction.

Healthy interfaces

Controlled impact tests on layered structure

Dynamic bulge tests have been used to approach the low-speed impact loading. Compared to the usual impact tests carried out on fall tower, the bulge tests are carried out on a hydraulic high speed device (Instron 65/20) in order to ensure as much as possible a constant strain rate. Moreover, for a speed loading less than 300 mm s -1 , it is possible to impose a cyclic loading which allows damage controls at given displacements. The bulge montage used for these tests is shown Figure 6.31. The circular area which is inside the clamping system has a diameter of 110 mm, whereas the pinch has a diameter of 40 mm. Three speed loading have been tested: 1.7 mm s -1 , 22 mm s -1 and 300 mm s -1 . After a first test which leads to the complete failure of the plate, the displacement value when the first fibre failure occurs is imposed as maximal displacement of the pinch. Thus, different plates have been inspected by means of C-scans to evaluate the eventual damage before the first fibre failure.

According to the ultrasonic inspections, the fibre failure of the upper ply is the first damage mechanism which occurs in the plate for the bulge tests. Despite of the transverse shear, no delaminations are observed.

Regarding the numerical simulations, the finite element model is presented Figure 6.32. The shell element formulation which has been used is the Belytschko-Tsay with reduced integration. The shell elements defining the specimen part have a size of 1 mm. The correlation between the numerical and the experimental responses are provided Figure 6.33. As recall, the results are predictions of the real behaviour of the plate based on a independent parameter identification of the ply behaviours. The results are relatively accurate but an unexpected re-hardening occurs numerically for a displacement of 2.25 mm. The cause has at this time not been found. Regarding the fibre failure, the same problem as for the interlaminar shear tests is present. The transverse shearing affects the ultimate failure strain and this phenomenon is not modelled. Thus, the fibre failure has not occurred in the present simulation.

Concluding remarks on the simulation of laminated structures

In a laminated composite, and because of stiffness heterogeneity through the thickness, the strain field is complex. The First-order Shear Deformation Theory (FSDT), which is widely used for shell element formulations in the finite element codes, leads to an incorrect bending stiffness and the impossibility to predict delamination due to transverse shear stresses.

In order to preserve the industrial framework and to avoid the use of several elements through-the thickness of the laminate material, a coupling between the usual Mindlin-Reissner (FSDT) element formulation and a higher-order zigzag deformation theory is implemented at the level of the material model. Moreover and without addition of new degrees of freedom, interface imperfections can be taken into account in the strain field. The first validation tests show promising results.

Thanks to an iterative scheme, a strain energy balance between the element formulation and the material model is achieved. Thus, by using traditional shell elements of commercial finite element code, it becomes possible to obtain a more realistic displacement field as input of the stress update algorithm (defined by the constitutive modelling of the behaviour of the fabric reinforced polymers). Moreover, the delamination is also predicted by calculation of the stresses at the interfaces. However, this method still suffers of a lack of efficiency. The numerous integration carried out by a Simpson's rule uses a lot of CPU time. Although numerous improvements have been done, such as the computation of the strain field, only one time per element instead of one time per integration point, the expression of the analytical solution of the integration is essential for a use in an industrial framework.

In regard to the delamination prediction, the transverse shear strain distribution obtained by the constitutive equation of the Enhanced Higher-order Shear Deformation Theory (EHSDT) are not describing precisely the exact solution. For a more accurate prediction, a strain equilibrium equation has to be used. However it requires the use of the spatial derivative of the stress field, not directly available in the finite element codes. A non-local formulation or the use of a quadratic shell element has to be used in parallel with the actual theory to accurately predict the delamination evolution.

Actually, the model is limited to a two-dimensional modelling, with consideration of the transverse shear stresses. An important way of improvement will be the development of this model in a three-dimensional framework, which will permit in addition the consideration of the mode I interlaminar crack opening. Most general simulations, such as high-speed impacts, crush, on most general structures with presence of curvatures will be possible.

Finally, a coupling between the intralaminar and the interlaminar damage can be easily implemented. Because of the presence of the strain theory and of the internal variables of the plies, such as matrix damage, inside the material model, it becomes possible to establish a bridge between both intra-and interlaminar damages in a finite element analysis. For that, a study on the laws which govern this coupling and on their identification have to be done.

CHAPTER 7

Conclusions and future perspectives

In order to reduce the energy consumption in the transportation industry (e.g. in automotive sector), the fabric reinforced polymers are more and more used because of their high specific stiffness and specific strength. Moreover, these materials present many damage mechanisms which lead to an important energy dissipation in case of crash and an improved impact tolerance compared to unidirectional composites. In order to ensure the protection of pedestrians and drivers/passengers in case of collision with non-structural components, such as a car hood, a roof or other car body elements, the behaviour of a layered composite structure made up with various fabric preforms subject to low-speed impact loading has been investigated in this thesis.

By decomposition of a low-speed impact loading in simple mechanical solicitations, the different physical mechanisms which occur during impacts have been highlighted and identified. Through the experimental observations and the literature review, a laminated fabric reinforced polymer presents intralaminar matrix damages, is sensitive to the strain-rate, suffers of delamination (or interlaminar damage) at the layer interfaces and finally fails because of the the fracture of the reinforcement, namely carbon fibres in the investigated material. Regardless the preform used as reinforcement in a ply, these mechanisms are strictly the same. All these mechanisms have been studied, modelled, identified as independently as was possible. Finally, the models of the different physical phenomena have been assembled and implemented into a user material subroutine of the commercial finite element code LS-DYNA ® . The important anisotropy as well as the eventual large rotation of the yarns during shearing yields to use a total Lagrangian framework. Hence, the material directions are well described all along the deformations and the objectivity is ensured.

The intralaminar damage model relies on the Onera Damage Model MicroStructure (ODM_MS). This strain-space formulated constitutive model is based on the assumption of privileged matrix crack directions. These directions corresponds to the fibre directions and a third damage variables is representative of the out-of plane intralaminar matrix cracking. This model has been extended in this thesis to take into account the friction effects at crack lips. Thus, the in-plane behaviour of the preforms, but also hysteresis loop observed in case of cyclic loading and evolution of the inelastic strains, both assumed to be due to friction mechanisms, have been successfully modelled by this updated matrix damage model.

The strain-rate dependency is introduced by a rheological viscoelastic model. To keep the strain-space formulation to ensure the efficiency of the model for explicit finite element analysis, the generalised Maxwell model is used. The coupling with the matrix damage model is introduced by considering that the strains which are applied to the viscoelastic model are the strains applied on an healthy area of the material. The model were validated by means of dynamic in-plane shear tests on an hydraulic high-speed device. Because of the specific specimen dimensions, keep identical to the quasi-static ones to limit scale effects, a new clamping system adapted to the composite materials have been designed. The numerical simulations of these tests show a good correlation with the experiments.

The fibre failure is implemented by considering two aspects: the failure criterion, and the fibre damage evolution. The fibre criterion derives from the maximal strain criterion. The maximal strain in the fibre direction is considered for every load cases as the ultimate strain leading to the final fracture of the ply. To introduce the dependence of shearing on the fibre failure criterion, the spectral decomposition of the strain tensor in a positive and in a negative part is used. Since the fibre failure leads to an important loss of stiffness, a softening occurs and causes a numerical strain localisation. To achieve an apparent mesh independence of the global amount of dissipated energy, without prevents the strain localisation nevertheless, a smeared crack approach is used. The energy dissipated is regularised by means of the introduction of a length representative of the element size in the formulation.

Finally, a laminate presents a specific distribution of strain through the thickness due to the stiffness heterogeneity. In the industrial framework where the use of single shell elements through the thickness of a plate is favoured, the bending behaviour is difficult to predict. In order to combine both usual element formulation and realistic strain distribution as input of the constitutive intralaminar model, a Enhanced higher-order Shear Deformation Theory has been implemented at the level of the user material subroutine. By the balancing of the strain energy between both theories, it becomes possible to use simultaneously both formulations: First-order Shear Deformation Theory (so-called Mindlin-Reissner theory) at the element level and Enhanced higher-order Shear Deformation Theory at the material level. Moreover, this newly introduced displacement field is able to consider the effect of interface degradation. Thus, the delamination is also taken into account. In addition, and because the strain distribution is computed in the user material subroutine, it becomes very simple to eventually couple the intralaminar damage with the interlaminar damage. This model is identified thanks to standardised interlaminar shear tests on single components. Finally the model is validated by means of interlaminar shear tests and dynamic bulge tests on industrial structure plates. The bulge test is used in order to test the layered composite plate under impact-type loading, but by controlling some input variables such as the strain-rate in order to be less dispersive.

Future perspectives

Future research can be focused on each physical phenomena which have been studied in this work.

In case of the intralaminar matrix damage, friction mechanisms have been implemented through the simplest friction law: Coulomb's law. These mechanisms could be improved by further differentiating static and dynamic friction. Influence of the strain-rate on the friction mechanisms may be studied to improve the modelling for dynamic loading cases. Experimentally these mechanisms can be observed by disrupted dynamic tests (by means of crossbow tests or by adding a fuse on the hydraulic high speed facility). The behaviour during unloading is thus due to viscoelasticity and friction mechanisms. Otherwise, the distinction between the stored strains, due to friction mechanisms, and the residual strains, due to local plasticity around the cracks, is complex. A possibility to isolate the mechanisms is to change the hydrostatic pressure inside the fabric reinforced materials to enforce the crack opening. Thus, bi-axial tests could be envisaged to understand more precisely the influence of the residual and the stored strains independently. Also, in case of bending in the fibre direction, a massive transverse crack occurs through the thickness of the plies and with an approximate angle of 45°. This crack is actually not well considered in the present model. A new privileged direction of matrix cracking could be added. It is also possible to consider the matrix damage under a tensor form in order to take into account the loading direction on the crack angle as described by Marcin [START_REF] Marcin | Modélisation du comportement, de l'endommagement et de la rupture de matériaux composites à renforts tissés pour le dimensionnement robuste de structures[END_REF]. It could finally be appropriate to reconsider the friction mechanisms as due to microscopic matrix damage. Hence, a distinction could be done between massive transverse matrix crack of the ply and the friction mechanisms due to local cracking.

The viscoelastic model parameters are identified by Dynamic Mechanical Analysis tests. However the large stiffness of the fibre reinforced plastics leads to limitation in the range of frequency because of device limitations. Moreover, the composite specimens because of internal friction are subject to important self-heating. Because of temperature dependency of the polymer matrix, the identification becomes incorrect. A special attention has to be paid on an enhanced parameter identification of the viscoelastic model. It can be done either by a modification of Dynamic Mechanical Analysis protocol, or by identification through other experimental procedures. For dynamic properties, dynamic tensile tests could be used to identify a surface depending on the stress, the strain and the strain rate from tests at different strain rates and by using the strain field heterogeneity (SEE method [START_REF] Lauro | The SEE method for determination of behaviour laws for strain rate dependent material: Application to polymer material[END_REF]).

The fibre failure criterion suffers of lack in case of transverse shear loading. The bending induces a strain and damage localisation leading to a premature fibre failure. This effect is not taken into account in the actual failure criterion. For fabric reinforced polymers, the literature does not provide numerous criteria but could be studied and implemented in the future. Otherwise, the simple criterion used in this model, based on the unique consideration of the fibre tension and some degradation of the ultimate strain by shear loading, could be adapted to differentiate in-plane and transverse shearing. It could be carried out by introducing a tensorial form as defined by Aiello [START_REF] Aiello | Utilisation des composites à matrice céramique SiC/SiC comme matériau de structure de composants internes du tore d'un réacteur à fusion[END_REF]. The dissipated energy by the fibre failure has not been identified precisely yet. Numerous experimental tests exist and need to be carried out.

Finally, the strain field reconstruction for laminated composite materials do not describe precisely the transverse shear strain distribution. To achieve a precise description of the strain field through this theory, an equilibrium equation is needed. It needs data on the spatial derivative of the stress fields which can be obtained by means of a non-local approach or a quadratic finite element formulation. The delamination is introduced by means of a spring-layer interface. In this thesis, the delamination criterion and the interlaminar damage evolution is introduced in a simple manner. The identification and the implementation of a physical interface law for the spring-layer interface is an important way of improvement. Moreover, the model is actually only formulated by neglecting transverse tensile and compressive strains. That makes the mode I opening of the interface impossible. The formulation and the implementation of this model in a full three-dimensional framework could permit a more general use of the present model. As last perspective, the integral forms used for the strain energy equilibrium have to be expressed analytically. Actually, they are computed through the CPU expensive Simpson's method which needs to be removed in order to reach an industrial efficiency.

Résumé

Les composites à matrice organique et renforcés par des préformes textiles (CMORT) sont en passe d'être déployés sur les véhicules de grandes séries pour réduire leur poids. Lorsqu'ils sont soumis à des impacts basse vitesse ces matériaux présentent des comportements complexes qui doivent être précisément modélisés et prédis au moyen de simulations par éléments finis. Dans ce but, un modèle matériau a été développé et implémenté dans un code éléments finis commercial. Soumis à un impact basse vitesse, un CMORT présente quatre mécanismes physiques majeurs qui altèrent la rigidité initiale du matériau : fissuration matricielle intralaminaire, rupture des fibres, délaminage et sensibilité à la vitesse de déformation. L'endommagement matriciel est modélisé grâce à un modèle constitutif reposant sur la mécanique de l'endommagement des milieux continus. Basé sur l'Onera Damage Model, il prend en compte les mécanismes de friction aux abords des fissures. La sensibilité à la vitesse de déformation est introduite au moyen d'un modèle de Maxwell généralisé. Ensuite, un critère de rupture est utilisé pour prédire l'initiation de la rupture des fibres et l'endommagement des fibres qui en découle est régularisé par l'utilisation d'un modèle de rupture progressive. Finalement, afin de prédire précisément le comportement hors-plan d'un stratifié, le calcul d'une distribution de déformation réaliste à travers l'épaisseur est réalisé au niveau du modèle matériau. Cette modélisation est capable de prendre en compte les effets du délaminage en utilisant seulement un élément coque. De plus, l'intégralité du modèle est formulé suivant la description Lagrangienne totale afin d'assurer l'objectivité et la cohérence matérielle durant la simulation. La procédure d'identification, ainsi que les tests de validation et les corrélations essais/simulations sont décrits pour chaque mécanisme physique. Enfin, le modèle est évalué au travers de la prédiction du comportement d'une structure automobile industrielle.

Mots-clés : Composites ; Modélisation ; Endommagement ; Viscoélasticité ; Délaminage

3. 7 47 3. 8 48 3. 9 58 4. 2 58 4. 3 61 4. 4 61 4. 5 62 4. 6 63 4. 7 100 5. 8 101 5. 9 117 6. 2 117 6. 3 118 6. 4 119 6. 5 2 . 2 2 . 3 FT 0 il 6 Green-iim( 3 )i 0 ii

 747848958258361461562663710081019117211731184119522230630 Tensile and in-plane shear specimen geometry. . . . . . . . . . . . . . . . . Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for tensile tests. . . . . . . . . . . . . . . . . . . . . . . . . Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for in-plane shear tests. . . . . . . . . . . . . . . . . . . . . 48 3.10 Local rotations of the preforms obtained through Digital Image Correlation observation during in-plane shear tests. . . . . . . . . . . . . . . . . . . . . 50 3.11 Reaction force comparisons between the numerical model and the experimental data for tensile tests in the fibre direction of the bi-axial non-crimp fabric. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.12 Reaction force comparisons between the numerical model and the experimental data for tensile tests in the fibre direction of the 3K plain weave woven fabric. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.13 Reaction force comparisons between the numerical model and the experimental data for tensile tests in the fibre direction of the 12K plain weave woven fabric. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.14 Reaction force comparisons between the numerical model and the experimental data for in-plane shear tests of the bi-axial non-crimp fabric. . . . . 54 3.15 Reaction force comparisons between the numerical model and the experimental data for in-plane shear tests of the 3K plain weave woven fabric. . . 54 3.16 Reaction force comparisons between the numerical model and the experimental data for in-plane shear tests of the 12K plain weave woven fabric. . 554.1 Standard rheological models for the viscoelasticity of solids. . . . . . . . . Generalised rheological models for the viscoelasticity of solids. . . . . . . . Basic rheological models for the viscoelasticity of solids. . . . . . . . . . . Strain and stress histories of a Maxwell element in a relaxation test. . . . . Zener model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Strain and stress histories of the Zener model in a relaxation test. . . . . . . Generalised Maxwell model. . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.8 DMA specimen geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.9 Sinusoidal strain input and stress response of a visoelastic material. . . . . . 72 4.10 Highlighting of the scale effect by the evolution of the specimen responses on in-plane shear tests. Two different specimen widths have been tested: 10 mm and 25 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 In-plane shear specimen geometry for dynamic tests. . . . . . . . . . . . . 4.12 Prismatic/screwing clamping system for dynamic tensile tests. . . . . . . . 4.13 High speed jack facility for dynamic tensile tests on fabric reinforced polymers. 4.14 Measured storage modulus of various preforms and the Prony series fit. . .4.15 Measured loss modulus of various preforms and the Prony series fit. . . . . 4.16 Modal analysis of the dynamic device by taking into consideration the piston, the clamping devices, the load cell and the specimen. . . . . . . . . . . . . 4.17 Visualisation of the non-homogeneous velocity field during simulation of a dynamic in-plane shear test on the fully modelled device. . . . . . . . . . . 4.18 Visualisation of the in-plane shear stress during simulation of a dynamic in-plane shear test on a single specimen. . . . . . . . . . . . . . . . . . . . 4.19 Comparison of the force responses of a dynamic in-plane shear test given by a simulation where the full device is modelled, or by a simulation where only the coupon is considered. . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Comparison between the experimental tests and the numerical model on in-plane shear tests at various dynamic strain rates for the non-crimp fabric composite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 Comparison between the experimental tests and the numerical model on in-plane shear tests at various dynamic strain rates for the 3K woven fabric composite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Comparison between the experimental tests and the numerical model on in-plane shear tests at various dynamic strain rates for the 12K woven fabric composite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Tensile specimen geometry used for failure identification. . . . . . . . . . . 5.2 Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for tensile tests used for failure identification. . . . . . . . . 5.3 Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for in-plane shear tests used for failure identification. . . . . 5.4 Compression specimen geometry for failure identification. . . . . . . . . . 5.5 Overview of the Wyoming modified ITRII compression test fixture. . . . . 5.6 Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for compression tests used for failure identification. . . . . . 5.7 Highlighting of the buckling phenomenon which may occur during compressive tests by plotting stress-strain responses of biaxial non-crimp fabric specimens. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Overview of the share of the dissipated energy in case of tensile fibre failure. Four point bending test and loading for the identification of the fracture toughness. Data reproduced from [95]. . . . . . . . . . . . . . . . . . . . . 106 5.10 Compression Tensile (CT) and Compression Compresson (CT) tests for the identification of the fracture toughness. Data reproduced from [94]. . . . . 106 5.11 Overview of the meshes used for tests of fibre failure regularisation thanks to the smeared crack model. . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.12 Responses of tensile tests along the fibre direction of various patch of different element sizes with or without damage regularisation. . . . . . . . . . . 109 5.13 Visualisation of the strain localisation even in the presence of fibre damage regularisation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.14 Responses of in-plane shear tests of various patch of different element sizes with damage regularisation. . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.15 Visualisation of the strain localisation during in-plane shear test even in the presence of fibre damage regularisation. . . . . . . . . . . . . . . . . . . . 111 6.1 Geometry and configuration of a laminated plate. . . . . . . . . . . . . . . Influence of the position of the reference plane Ω on the design of the laminate.Representation of the third-order displacement field through the thickness of the laminate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Representation of the enhanced higher-order displacement field through the thickness of the laminate. . . . . . . . . . . . . . . . . . . . . . . . . . . . Representation of the enhanced higher-order displacement field through the thickness of the laminate with delamination. . . . . . . . . . . . . . . . . . 125 6.6 3-points bending montage for interlaminar shear tests on fabric reinforced polymers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 6.7 Specimen geometry used for interlaminar shear tests of non-crimp fabric and 3K woven preforms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 6.8 Specimen geometry used for interlaminar shear tests of 12K woven preforms. 138 6.9 Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for interlaminar shear tests. . . . . . . . . . . . . . . . . . . 138 6.10 Influence of the bending direction on 3K woven fabric specimens. . . . . . 140 6.11 Influence of the bending direction on non-crimp fabric specimens. . . . . . 140 List of tables 2.1 Types of carbon fibres. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Properties of the carbon fibres. . . . . . . . . . . . . . . . . . . . . . . . . 7 Stacking sequences of the plate types. . . . . . . . . . . . . . . . . . . . . 15 2.4 Stacking sequences according to the directions of reinforcement of the plate type. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.5 RTM cycle informations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.6 Intrinsic plate informations. . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1 Parameters for the intralaminar matrix damage model for the various fabric preforms. Expressed in SI (N, mm and s). . . . . . . . . . . . . . . . . . . 51 4.1 Viscoelastic parameters of the damageable viscoelastic model for the various fabric preforms. Expressed in SI (N, mm and s). . . . . . . . . . . . . . . . 75 4.2 Intralaminar matrix damage parameters of the damageable viscoelastic model for the various fabric preforms. Expressed in SI (N, mm and s). . . . . . . . 76 6.1 Interlaminar matrix damage parameters for the various fabric preforms. Expressed in SI (N, mm and s). . . . . . . . . . . . . . . . . . . . . . . . . Component of a second-order tensor x x Fourth-order tensor ∂ (•) ∂ a , (•) ,a Derivative of (•) with respect of a ( •) Material time derivative of (•) (•) T Transpose of (•) : Double contraction of tensors |x| Absolute value of x ∥⃗ x∥, ∥x∥ Euclidean norm of vectors and tensors < x > Macauley bracket, i.e. < x >= |x|+x 2 δ i j Kronecker symbol, i.e. δ i j = 1 for i = j, else null ⃗ e Green-Lagrange strain vector (Voigt notation) Roman Symbols dl Length of a position vector dl 0 after deformation dl 0 Length of a position vector in the ⃗ N direction in the base configuration ⃗ dX Position vector of material points in undeformed configuration ⃗ dx Position vector of material points in deformed configuration Deformation gradient tensor E Green-Lagrange strain tensor ⃗ s Second Piola-Kirchhoff stress vector (Voigt notation) curing Curing temperature t curing Curing time t handling Handling time t injection Injection duration t total Process time V f ,theo Theoretical fibre volume ratio NOMENCLATURE XIX W f ,theo Theoretical fibre weight ratio Intralaminar matrix damage parameters e Green-Lagrange strain at crack closure in the iof the in-plane shear stiffness due to shear locking e Lagrange shear strain at shear locking (Voigt notation) G 13 Out-of-plane longitudinal shear modulus G 23 Out-of-plane transverse shear modulus ∆e 0 Initial Green-Lagrange strain interval of crack closure in the i-direction a m Effect of the damage on the strain interval of crack closure in the i-direction h I m(i) Effect of mode I loading on cracks normal to the i-direction h II m(i) Effect of mode II loading on cracks normal to the i-direction h III m(i) Effect of mode III loading on cracks normal to the i-direction h p Effect of the out-of-plane damage on the in-plane shear stiffness b i Weight of the transverse shear strain on the thermodynamic force in the i-direction d x c(i) Critical damage value for cracks normal to the i-direction p x Severity of damage emergence for cracks normal to the i-direction y x c(i) Thermodynamic force regulation parameter for cracks normal to the i-direction y x 0(i) Thermodynamic force at damage activation for cracks normal to the i-direction ζ i Residual strain rate for cracks normal to the i-direction µ Initial friction coefficient for cracks normal to the i-direction m m Damage effect on the friction coefficient for cracks normal to the i-direction Viscoelastic parameters E ve j(i) Viscoelastic stiffness in direction i related to the j-th viscoelastic mechanism G ve j(ik) Viscoelastic shear stiffness in directions i and k related to the j-th viscoelastic mechanism τ j Relaxation time related to the j-th viscoelastic mechanism Fibre failure parameters c fx i Influence of shearing on fibre failure in the i-direction y fx 0(i) Thermodynamic force at damage activation for fibre failure in the i-direction Interlaminar matrix damage parameters R c(k) Critical compliance of the k-th interface σ max(k) zα Ultimate transverse shear stress of the k-th interface τ (k) Rate of interlaminar damage of the k-th interface Acronyms / Abbreviations CDM Continuum Damage Mechanics CFRP Carbon Fabric Reinforced Polymer DMA Dynamic Mechanical Analysis FEA Finite Element Analysis FRP Fabric Reinforced Polymer CHAPTER 1

Figure 2 . 1 :

 21 Figure 2.1: Fibre grades according to the tensile modulus and the tensile strength. Data reproduced from [typ].

Figure 2 . 2 :

 22 Figure 2.2: Tensile stress-strain response of epoxy at different strain rates. Data reproduced from [42].

Figure 2 . 3 :

 23 Figure 2.3: Compressive stress-strain response of epoxy at different strain rates. Data reproduced from [41].

Figure 2 . 4 :

 24 Figure 2.4: Torsion stress-strain response of epoxy at different strain rates. Data reproduced from [43].

Figure 2 . 5 :

 25 Figure 2.5: Overview of the non-crimp fabric preform.

Figure 2 . 6 :

 26 Figure 2.6: Overview of the 3K woven textile preform.

Figure 2 . 7 :

 27 Figure 2.7: Overview of the 12K woven textile preform.

  Figure 2.8.

Figure 2 . 8 :

 28 Figure 2.8: Overview of the preform sequence inside the considered laminated structure: (0°/90°) 12K , (+45°, -45°) NCF , (-45°, +45°) NCF , (0°/90°) 3K (given from the lowest ply to the uppermost ply).

( a )Figure 2 . 10 :

 a210 Figure 2.10: Response types during impact on plates. Data reproduced from [87].

Figure 2 . 11 :

 211 Figure 2.11: Experimental observations of intralaminar matrix damage (blue), fibre failure (green) and delamination (red) which occur during bending of a laminate composite.

Figure 2 . 12 :

 212 Figure 2.12: The major configurations used in finite strain theory.

Figure 2 . 13 :

 213 Figure 2.13: Deformation of a body with representation of base and current configurations.

Figure 2 . 14 :

 214 Figure 2.14: The major configurations used in Updated Lagrangian framework.

Figure 2 Figure 2 . 15 :

 2215 Figure 2.15: The major configurations used in Total Lagrangian framework.

  assumption is confirmed by basic experimental observations on two different fabric preforms after an in-plane shear test (Figure 3.1). (a) Initial view of biaxial non-crimp fabric (left) and once damaged (right). (b) Initial view of 3K plain-weave woven fabric (left) and once damaged (right).

Figure 3 . 1 :

 31 Figure 3.1: Global overview of cracks following an in-plane shear loading on various fabric reinforced polymers by observation of the external surface.

e 1 e 2 e 3 e 4 e 5 e 6 T= E 11 E 22 E 33

 6112233

4 Figure 3 . 2 :

 432 Figure 3.2: Imposed complex displacements for visualisation of the stored strain effect.

Figure 3 . 3 :

 33 Figure 3.3: Visualisation of the stored strain effect on the continuity of the stress response.

  Figure 3.4: Global overview of cracks following an in-plane shear loading on various fabric reinforced polymers by observation of the cross-section.

Figure 3 . 5 :

 35 Figure 3.5: Visualisation of the stored strain effect on the continuity of the stress response.

Figure 3 . 6 :

 36 Figure 3.6: Visualisation of the stored strain effect on the continuity of the stress response.

. 7 .Figure 3 . 7 :

 737 Figure 3.7: Tensile and in-plane shear specimen geometry.

Figure 3 . 8 :

 38 Figure 3.8: Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for tensile tests.

Figure 3 . 10 :

 310 Figure 3.10: Local rotations of the preforms obtained through Digital Image Correlation observation during in-plane shear tests.

  Square Error between the numerical and the experimental responses. The sampling on the parameter range is done by using a D-Optimal method and the optimisation algorithm which has been used is the Adaptive Simulated Annealing.

Figure 3 . 11 :

 311 Figure 3.11: Reaction force comparisons between the numerical model and the experimental data for tensile tests in the fibre direction of the bi-axial non-crimp fabric.

  , 3.15 and 3.16.

Figure 3 . 12 :Figure 3 . 13 :Figure 3 . 14 :Figure 3 . 15 :Figure 3 . 16 :

 312313314315316 Figure 3.12: Reaction force comparisons between the numerical model and the experimental data for tensile tests in the fibre direction of the 3K plain weave woven fabric.

Figure 4 .

 4 1b) models. However these models σ σ (a) Boltzamnn model σ σ (b) Zener model

Figure 4 . 1 :

 41 Figure 4.1: Standard rheological models for the viscoelasticity of solids.

Figure 4 . 2 :

 42 Figure 4.2: Generalised rheological models for the viscoelasticity of solids.

Figure 4 . 3 :

 43 Figure 4.3: Basic rheological models for the viscoelasticity of solids.

Figure 4 . 4 :

 44 Figure 4.4: Strain and stress histories of a Maxwell element in a relaxation test.

Figure 4 . 5 )Figure 4 . 5 :

 4545 Figure 4.5: Zener model.

Figure 4 . 6 :

 46 Figure 4.6: Strain and stress histories of the Zener model in a relaxation test.

Figure 4 . 8 :

 48 Figure 4.8: DMA specimen geometry.

Figure 4 . 9 :

 49 Figure 4.9: Sinusoidal strain input and stress response of a visoelastic material.

Figure 4 . 10 :

 410 Figure 4.10: Highlighting of the scale effect by the evolution of the specimen responses on in-plane shear tests. Two different specimen widths have been tested: 10 mm and 25 mm.

Figure 4 . 11 :

 411 Figure 4.11: In-plane shear specimen geometry for dynamic tests.

Figure 4 . 12 :

 412 Figure 4.12: Prismatic/screwing clamping system for dynamic tensile tests.

( a )

 a Overview of the dynamic tensile tests with high speed camera and optical extensometer for displacement field measurements. (b) Details on the specific clamping system used for dynamic composite tensile tests.

Figure 4 . 13 :

 413 Figure 4.13: High speed jack facility for dynamic tensile tests on fabric reinforced polymers.

Figure 4 . 14 :

 414 Figure 4.14: Measured storage modulus of various preforms and the Prony series fit.

Figure 4 . 15 :

 415 Figure 4.15: Measured loss modulus of various preforms and the Prony series fit.

Figure 4 . 16 :

 416 Figure 4.16: Modal analysis of the dynamic device by taking into consideration the piston, the clamping devices, the load cell and the specimen.

Figure 4 . 17 :

 417 Figure 4.17: Visualisation of the non-homogeneous velocity field during simulation of a dynamic in-plane shear test on the fully modelled device.

Figure 4 . 18 :

 418 Figure 4.18: Visualisation of the in-plane shear stress during simulation of a dynamic in-plane shear test on a single specimen.

Figure 4 .Figure 4 . 19 :

 4419 Figure 4.19: Comparison of the force responses of a dynamic in-plane shear test given by a simulation where the full device is modelled, or by a simulation where only the coupon is considered.

Figure 4 . 20 :

 420 Figure 4.20: Comparison between the experimental tests and the numerical model on in-plane shear tests at various dynamic strain rates for the non-crimp fabric composite.

Figure 4 . 21 :

 421 Figure 4.21: Comparison between the experimental tests and the numerical model on in-plane shear tests at various dynamic strain rates for the 3K woven fabric composite.

Figure 4 . 22 :

 422 Figure 4.22: Comparison between the experimental tests and the numerical model on in-plane shear tests at various dynamic strain rates for the 12K woven fabric composite.

. 2 )⃗ e = e 1 e 2 e 3 e 4 e 5 e 6 T= E 11 E

 2611 which becomes the Green-Lagrange strain vector ⃗ e 22 E 33 2E 23 2E 13 2E 12 T .

Figure 5 . 1 :

 51 Figure 5.1: Tensile specimen geometry used for failure identification.

Figure 5 . 2 :

 52 Figure 5.2: Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for tensile tests used for failure identification.

Figure 5 Figure 5 . 4 :

 554 Figure 5.4: Compression specimen geometry for failure identification.

  The coupon is clamped on both sides on 50 mm length. The test setup used for the present work is a Wyoming modified ITRII compression test fixture (Figure5.5) in accordance with the NF EN ISO 14126 standard.

Figure 5 . 5 :

 55 Figure 5.5: Overview of the Wyoming modified ITRII compression test fixture.

Figure 5 . 6 :

 56 Figure 5.6: Material coordinate system ( ⃗ 1, ⃗ 2) according to the specimen coordinate system ( ⃗ X, ⃗ Y ) for compression tests used for failure identification.

Figure 5 . 7 :Figure 5 . 8 :

 5758 Figure 5.7: Highlighting of the buckling phenomenon which may occur during compressive tests by plotting stress-strain responses of biaxial non-crimp fabric specimens.

Figure 5 . 9 :

 59 Figure 5.9: Four point bending test and loading for the identification of the fracture toughness. Data reproduced from [95].

Figure 5 . 10 :

 510 Figure 5.10: Compression Tensile (CT) and Compression Compresson (CT) tests for the identification of the fracture toughness. Data reproduced from [94].

Figure 5 . 11 :

 511 Figure 5.11: Overview of the meshes used for tests of fibre failure regularisation thanks to the smeared crack model.

  No regularisation of the fibre damage (κ ft i = cte). Regularisation of the fibre damage (κ ft i = f (l)).

Figure 5 . 12 :

 512 Figure 5.12: Responses of tensile tests along the fibre direction of various patch of different element sizes with or without damage regularisation.

Figure 5 . 13 :

 513 Figure 5.13: Visualisation of the strain localisation even in the presence of fibre damage regularisation.

Figure 5 . 14 :

 514 Figure 5.14: Responses of in-plane shear tests of various patch of different element sizes with damage regularisation.

Figure 5 . 15 :

 515 Figure 5.15: Visualisation of the strain localisation during in-plane shear test even in the presence of fibre damage regularisation.

Figure 6 . 1 :

 61 Figure 6.1: Geometry and configuration of a laminated plate.

Section 1 Section 2 (

 2 a) Ω = 0 for both sections. Ω = -1 for both sections.

Figure 6 . 2 :

 62 Figure 6.2: Influence of the position of the reference plane Ω on the design of the laminate.

x 3 u α u 0 αFigure 6 . 3 :

 063 Figure 6.3: Representation of the third-order displacement field through the thickness of the laminate.

Figure 6 . 4 :

 64 Figure 6.4: Representation of the enhanced higher-order displacement field through the thickness of the laminate.

  γ ) the compliance of the d-th interface layer. Cheng et al.[START_REF] Cheng | Modelling of weakly bonded laminated composite plates at large deflections[END_REF], Sciuva et al.[START_REF] Sciuva | Effects of Interfacial Damage on the Global and Local Static Response of Cross-Ply Laminates[END_REF] established a correlation between the compliance of the interface layer and the physical displacements at interfaces. if R (d) αβ (x γ ) and R (d) 33 (x γ ) → 0, the coupling between two adjacent plies is perfect. When delamination occurs R (d) αβ (x γ ) and R (d)33 (x γ ) rise to a given value. In case of infinite value for the compliance of the spring-layer interface, the plies are moving independently and without friction.

  zα the transverse shear stresses at the out-of-plane position x

Figure 6 . 6 :

 66 Figure 6.6: 3-points bending montage for interlaminar shear tests on fabric reinforced polymers.

Figure 6 . 7 :

 67 Figure 6.7: Specimen geometry used for interlaminar shear tests of non-crimp fabric and 3K woven preforms.

  (a) Transverse shear fracture inside 3K woven fabric specimen for interlaminar shear tests in the fibre direction. (b) Delamination inside 3K woven fabric specimen for interlaminar shear tests with a angle of 45°according to the fibre direction.

Figure 6 . 10 :

 610 Figure 6.10: Influence of the bending direction on 3K woven fabric specimens.

( a )

 a Transverse shear fracture inside noncrimp fabric specimen for interlaminar shear tests in the fibre direction. (b) Delamination inside non-crimp fabric specimen for interlaminar shear tests with a angle of 45°according to the fibre direction.

Figure 6 . 11 :

 611 Figure 6.11: Influence of the bending direction on non-crimp fabric specimens.

Figure 6 . 12 :

 612 Figure 6.12: Comparison between the experimental tests and the numerical model of transverse shear tests in the fibre direction for the non-crimp fabric composite.

Figure 6 . 13 :Figure 6 . 14 :Figure 6 . 15 :Figure 6 . 16 :Figure 6 . 17 :

 613614615616617 Figure 6.13: Comparison between the experimental tests and the numerical model of transverse shear tests at 45°for the non-crimp fabric composite.

Figure 6 . 18 :

 618 Figure 6.18: Split view of the stacking sequence of the investigated structure.

Figure 6 . 19 :Figure 6 . 20 :Figure 6 . 21 :

 619620621 Figure 6.19: Overview of the model for simulations of interlaminar shear tests.

  -plane position (mm) Longitudinal in-plane stress (MPa)

Figure 6 . 22 :

 622 Figure 6.22: Longitudinal in-plane stress through the thickness of a 0°specimen before the first delamination (displacement of 0.5 mm).

  and 6.27), the real interlaminar damage visualised by mean of ultrasonic measurements (Fig--plane position (mm) Longitudinal in-plane stress (MPa)

Figure 6 . 23 :

 623 Figure 6.23: Longitudinal in-plane stress through the thickness of a 0°specimen after delamination (displacement of 0.6 mm).

  -plane position (mm) Longitudinal in-plane stress (MPa)

Figure 6 . 24 :

 624 Figure 6.24: Longitudinal in-plane stress through the thickness of a 45°specimen before the first delamination (displacement of 0.6 mm).

  -plane position (mm) Longitudinal in-plane stress (MPa)

Figure 6 . 25 :

 625 Figure 6.25: Longitudinal in-plane stress through the thickness of a 45°specimen after delamination (displacement of 0.7 mm).

Figure 6 . 26 :

 626 Figure 6.26: Visualisation of the interlaminar damage predicted by simulation of interlaminar shear tests according to the 0°direction. The three interfaces are given from the lower interface to the upper one (from left to right).

Figure 6 . 27 :

 627 Figure 6.27: Visualisation of the interlaminar damage predicted by simulation of interlaminar shear tests according to the 45°direction. The three interfaces are given from the lower interface to the upper one (from left to right).

Figure 6 . 28 :

 628 Figure 6.28: Ultrasonic inspections on a structure specimen before interlaminar shear tests where no damage is apparent.

Figure 6 . 29 :

 629 Figure 6.29: Ultrasonic inspections for damage visualisation on a structure specimen after interlaminar shear tests according to the 0°structure direction.

Figure 6 . 31 :

 631 Figure 6.31: Overview of the bulge montage for multi-axial tests on the structure specimen.

Figure 6 . 32 :

 632 Figure 6.32: Overview of the model for simulations of bulge tests.

Figure 6 . 33 :

 633 Figure 6.33: Comparison between the experimental and the predicted numerical responses of bulge test at 300 mm s -1 for the investigated structure.

Table 2

 2

	Fibre type	Notation	Tensile Modulus Tensile Strength
			GPa	MPa
	Ultra high elastic modulus High elastic modulus Intermediate elastic modulus Standard elastic modulus Low elastic modulus	UHM HM IM HT LM	≥ 600 ≥ 350 & ≤ 600 ≥ 280 & ≤ 350 ≥ 200 & ≤ 280 ≤ 200	≥ 2500 ≥ 2500 ≥ 3500 ≥ 2500 ≤ 3500

.1: Types of carbon fibres.

Table 2 . 2 :

 22 Properties of the carbon fibres.

	Fibre name	Tensile Modulus Tensile Strength Strain Density Diameter GPa MPa % g cm -3 µm
	Aksaca ® A-42	240	4200	1.8	1.76	6.9
	Pyrofil ™ TR30S	235	4410	1.7	1.79	7
	Torayca ® T700S	230	4900	2.1	1.79	7

Table 2 . 3 :

 23 Stacking sequences of the plate types.

Table 2 . 4 :

 24 Stacking sequences according to the directions of reinforcement of the plate type.

	Plate type	Reinforcement stacking sequence
	Structure NCF	(0°/90°) 12K , (+45°, -45°) NCF , (-45°, +45°) NCF , (0°/90°) 3K [(±45°) 3 ] S
	3K woven 12K woven	(0°/90°) 3 , 0°/90° [(0°/90°) 4 ] S	S

Table 2 . 5 :

 25 RTM cycle informations.

	Plate type	t injection t curing t handling t total T curing s s s s • C
	Structure NCF 3K woven 12K woven ≈ 26 ≈ 19 ≈ 27 ≈ 26	120 120 120 120	≈ 10 ≈ 149 120 ≈ 10 ≈ 157 120 ≈ 10 ≈ 156 120 ≈ 10 ≈ 156 120

Table 2 . 6 :

 26 Intrinsic plate informations.

	Plate type	Thickness m f	m r V f ,theo W f ,theo
		mm	g	g	%	%
	Structure	1.40	168 185 38.5	48.8
	NCF	2.30	306 270 45.4	55.5
	3K woven	2.10	250 260 38.3	48.6
	12K woven	2.00	245 250 41.7	47.9

  .6).

	300		s 1 whithout friction s 6 whithout friction			s 1 with null friction s 6 with null friction
	250
	200
	-250 -200 -150 -100 -50 0 50 100 150 Stress (MPa)	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
							Time (s)

Table 3 . 1 :

 31 Parameters for the intralaminar matrix damage model for the various fabric preforms. Expressed in SI (N, mm and s).

	Paramaters	NCF	3K woven 12K woven
	E 0 1 = E 0 2 α 1 = α 2	55562 277722	46589 147722	49212 286678
	ν 21	0.13	0.10	0.05
	G 0 12 g l 12	3000 0	4910	3236

Table 4 .

 4 1: Viscoelastic parameters of the damageable viscoelastic model for the various fabric preforms. Expressed in SI (N, mm and s).

	Paramaters	NCF	3K woven 12K woven
	τ 1	100	100	100
	G ve 1(12)	351.19	325.46	164.55
	τ 2	17.78	17.78	17.78
	G ve 2(12)	137.3	23.59	87.5
	τ 3	3.16	3.16	3.16
	G ve 3(12)	174.22	103.3	155
	τ 4	0.562	0.562	0.562
	G ve 4(12)	207.2	73.68	126.70
	τ 5	0.1	0.1	0.1
	G ve 5(12)	122.6	84.74	109.6
	τ 6	0.018	0.018	0.018
	G ve 6(12)	120.58	76.046	159.72
	τ 7	0.00316 0.00316	0.00316
	G ve 7(12)	81.421	39.93	247.36

Table 4 . 2 :

 42 Intralaminar matrix damage parameters of the damageable viscoelastic model for the various fabric preforms. Expressed in SI (N, mm and s).

	Paramaters	NCF	3K woven 12K woven
	E 0 1 = E 0 2 α 1 = α 2	61062 277722	46589 147722	49212 286678
	ν 21	0.13	0.10	0.05
	G 0 12 g l	3000	3000	3000

Table 6 . 1 :

 61 Interlaminar matrix damage parameters for the various fabric preforms. Expressed in SI (N, mm and s).

  and the thickness of the plate is around 1.40 mm. This stacking sequence is

	(N)
	Force	800
		600
		400
		0 200	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7 Experiment Simulation	0.8	0.9
							Displacement (mm)

Acknowledgements

CHAPTER 5

Constitutive modelling of the final fracture

The final fracture of the bi-axial fabric reinforced polymers is due to the failure of fibres in one direction. To well-describe the fibre failure and the dissipated energy, two mechanisms have to be taken into account: the fracture initiation and the post-failure behaviour. The failure criterion is based on the simple maximal strain criterion, applied in the fibre directions in tension and compression. A double spectral decomposition leads to the introduction of a shear dependence on the ultimate strain of the fibres. Then, this criterion is placed in the adopted formalism (Onera Damage Model). The strain-softening due to the important loss of stiffness induced a mesh dependency phenomenon in the finite element analysis. This mesh dependence is managed by using a energy dissipated per unit volume. This method, so-called smeared crack, uses the element size given by the finite element code, to regularise the fibre damage rate.

Introduction

In the general case of the layered continuous fibre reinforced polymers the final fracture is induced by the fibre failure. The modelling of the failure lies in understanding two mechanisms: the failure initiation and the failure propagation.

The damage initiation is evaluated by a failure criterion. Various criteria have been proposed in the literature. The most classical are the maximal stress criterion, or its dual where y ft 1 and y ft 2 are the thermodynamic forces associated to the tensile fibre failures y fc 1 and y fc 2 are the thermodynamic forces associated to the compressive fibre failures and C 0 is the elastic second-order stiffness tensor. The parameters c ft 1 , c fc 1 , c ft 2 and c fc 2 adjust the influence of shearing on the fibre failure. From that, the thermodynamic forces at failure initiation are given by: Longitudinal fibres Tension y ft 0(1) = .

(5.9)

Identification

The number of parameters to identify in order to describe the failure envelope is restricted to 8 parameters (Equations 5.8 and 5.7): y ft 0(1) , y fc 0(1) , y ft 0(2) , y fc 0(2) , c ft 1 , c fc 1 , c ft 2 and c fc 2 , namely the thermodynamic forces and the parameters describing the influence of shearing related to longitudinal and transverse tensile and compressive fibre failure. In case of a balanced fabric preform, the number of parameters to identify is reduced to four: the thermodynamic forces and the shearing effect parameters related to the tensile and compressive fibre failure.

The parameter identification of the fibre failure criterion is done by means of three different experimental tests carried out independently on the different preforms. These tests are the standardised monotonic tensile, in-plane shear and compression tests, respectively the NF EN ISO 527-4, the NF EN ISO 14129 and NF EN ISO 14126 tests, in order to be easily reproducible in an industrial framework. Each test provides the value of one parameter which makes identification immediate.

Tensile tests

These tests are carried out at room temperature on an electromagnetic device (Sintech 20D) which ensures 100 kN for the maximal load capacity and at a speed of 2 mm min -1 . The and the odd lines by:

and where S, B 1 and B 2 are defined as follows:

. . .

. .

. . .

and

. . .

As a result, the slop angles can be determined by means of only the material constants (transverse shear stiffness), the geometric description of the laminate (thickness and out-ofplane position of the plies) and the variable φ α (x β ) as follows:

which leads to:

where the a

αβ parameters are function of the material parameters around the k-th ply interface.

As a consequence, the displacement field can now be expressed in a simple form and is given by:

where Φ αβ (x 3 ) is a function of the out-of-plane coordinate and represents the in-plane displacements through the thickness of the laminate due to transverse shearing, and is given by:

with δ αβ the Kronecker symbol defined by: δ αβ = {1 for i = j, 0 for i ̸ = j . (6.29)

Influence of the interlaminar damage

Cho and Kim [START_REF] Cho | Higher-Order Zig-Zag Theory for Laminated Composites With Multiple Delaminations[END_REF] presented an extension of the Enhanced Higher-order Shear Deformation Theory to take into account the delamination. Possible jumps in the slipping or in the opening at the interfaces are introduced through the displacement variable u

i (x α ) as shown Figure 6.5. As a consequence the displacement field is now given by:

.5: Representation of the enhanced higher-order displacement field through the thickness of the laminate with delamination.

Since the slop angles are not affected by delamination, the simplified in-plane displacement field is now expressed such as:

(6.32)

In order to eliminate the newly introduced degrees of freedom, a springer-layer model is employed and the interlaminar displacement variables at interfaces are expressed as follows:

Abstract

Carbon Fabric Reinforced Polymers (CFRP) will soon used in high volume automotive production in order to reduce the vehicle weight. For safety and design reasons, their complex behaviours under low-speed impacts, such as pedestrian impacts, need to be accurately modelled and predicted by finite element simulations. For this purpose, a material model dedicated to explicit finite element simulations has been developed and implemented in a commercial finite element code. Subject to low-speed impacts, the CFRP shows four different physical mechanisms which alter the initial stiffness of the material: intralaminar matrix cracks, fibre failure, delamination and strain-rate sensitivity. The intralaminar damage is modelled through constitutive equations based on the continuum damage theory. It is based on the Onera Damage Model, but with the consideration of friction mechanisms between crack lips in order to represent the hysteresis loops in case of cyclic loading. The strain-rate sensitivity is introduced by means of the rheological generalised Maxwell viscoelastic model. Regarding the fibre damage, a failure criterion based on the strain of the fibre direction is introduced. The energy release due to the fibre failure is also regularised thanks to a smeared crack approach. Finally, in order to welldescribed the out-of-plane behaviour, such as bending, of a laminated CFRP material, a recomputation of a realistic strain field through-the thickness of the laminate is introduced at level of the material model. Based on strain energy equilibrium between usual shell element theory and higher-order zigzag theory, this formulation is able to consider delamination at ply interfaces by using only one shell element through-the-thickness of a laminate. In addition, the model is placed in a total Lagrangian framework to ensure both objectivity and material coherence. The identification procedure, with the needed experimental tests, as well as validation tests and experimental/numerical correlations are given for all physical mechanisms previously described. Finally, this model is evaluated through the behaviour prediction of an industrial structure.