Multi-layered multi-material finite element for crashworthiness studies

The materials used in the field of ground transportation are becoming increasingly complex and it is usual to find numerous composite parts. The potential for utilization is generally verified through the numerical simulations, especially due to the increase in severe safety criteria. This paper deals with a finite element intended to describe numerically the behaviour of multi-layered multi-materials (combination of metallic and composite plies). A complete description of this element is presented with the taking into account of the resin ductility and the asymmetrical tension/compression behaviour for composite plies. The latter are performed by the introduction of an isotropic hardening law and a compressive multiplier depending on strain. Some numerical validations will be produced to verify the developments.


Introduction
Composite materials have benefited from a remarkable growth in the field of transportation, such as aeronautic, aerospace, automobile, etc. Indeed, their advantages in low density and high specific properties give this type of material an increasingly important role. Therefore, new materials, particularly multi-layered multi-material, have emerged. In most cases, they are a combination of metallic and composite plies but it is also possible to incorporate them with foam materials. Many experimental studies have shown the importance of composite materials applied to the reinforcement of metallic components [1,2], especially in regard to their energy absorption capabilities and weight saving.
Today due to severe safety criteria, standard conformity, component behaviour prediction, etc. numerical simulations are being intensively developed, principally in the domains of transportation. So to foresee the behaviour of such metallic/composite materials, it seems paramount to have a reliable tool for rapid dynamic numerical simulations.
A new finite element has been developed within Pam Crash software [3], in collaboration with Engineering Systems International and Pam System International firms.
It allows the possibility of taking into account elastic plastic damaging behaviour for metallic plies and elastic damaging behaviour for composite plies. For these plies, two models are available: firstly, the composite ply is considered as a whole (homogeneous model [4,5]) and secondly, the behaviours of fibres and matrix are distinguished separately (bi-phase model [6][7][8]).
Most of the conducted numerical validations permit the necessary confirmation of such a finite element. Therefore, other types of ply (foam, fabric composite, etc.) have been integrated and this element currently offers a more variable panel for describing materials used in the fields of transportation.
However, the research program is restricted to the composite ply description and particularly to the homogeneous model optimization. Indeed, the complexity and diversity of these materials or the loading on structures implies that an elastic damaging model is sometimes not refined enough to account for experimental phenomena. For example, the matrix ductility, the strain rates or the delamination could have an effect on composite ply behaviour. With a view to integration of the strain rates dependence on composite behaviour, a plasticity algorithm and an asymmetric behaviour in tension/compression have been added to the element 131 for homogeneous model composite ply.
This paper starts with a presentation of the multi-layered multi-material finite element. Then it continues with a theoretical part concerning plastic phenomena of resin and asymmetrical compression/tension behaviour. After that, some examples of numerical validations are given. They are based on a tensile test on a ^45 2S laminate, a threepoint dynamic bending on a 0 2 90 2 S composite plate and an axial compressive test on steel/composite materials. These specific tests highlight the accuracy of the numerical achievements according to the described phenomena (matrix ductility for the tensile test, compression for the bending test and whole behaviour for the axial compressive test). Finally, the conclusions will guide future studies through a planned methodology in order to introduce strain rates effects for the composite ply.

Presentation of the multi-layered multi-material element
This element [3] is a shell of the Mindlin/Reissner type. It includes in-plane motion and features bi-linear interpolation functions with four nodes, which possess six degrees of freedom. Three different types of ply have been generated. The first two plies (metallic and "bi-phase" composite) derive from the descriptions of explicit code available. In addition, the last one ("homogeneous" composite) has been integrated in order to propose another approach for the composite material.
Each of these plies possesses individual mechanical characteristics and a single integration point located in the middle of the shell. Thus, the stacking and the nature of component materials are reproduced inside the shell element.

Metallic ply modelling
Metallic plies [9,10] display an elastic plastic behaviour, initially isotropic. The elastic behaviour is defined via Young's modulus, Poisson's ratio, the shear modulus, and the thickness of the ply, whereas the plastic behaviour is given by a stress-strain curve data. This curve is defined via seven points, for each of them the plastic tangent modulus and plastic stress is entered (Fig. 1). This type of plasticity definition permits the use of an analytical strain rate law that is here a Cowper-Symonds method: where s y is the yield stress, _ e the strain rate and p and D are parameters.

Composite ply modelling
Two possibilities are available to describe composite ply behaviour. The first method (bi-phase model) allows fibres and matrix behaviour to be distinguished; composite behaviour is then considered as the superposition of effects from the two phases. The second method (homogeneous model) proposes an overall description of the ply.

Bi-phase model
The bi-phase model [9,10] is a model of heterogeneous materials adapted to unidirectional composite materials reinforced by continuous fibres. This model [6,7] was introduced due to a crash behaviour study of car made entirely in composite [8]. Essentially, the key lies in distinguishing the two elementary layer components: fibres (uni-dimensional phase) and matrix (orthotropic phase) (Fig. 2).
Consequently, each phase features its specific rheological characteristics: fibres have a unidirectional behaviour that is brittle elastic with damage whereas the matrix has behaviour which is brittle orthotropic elastic or elastic damaging. A description of each law is given in tensile and compressive responses considering the composite asymmetric behaviour.
The effects' superposition of each phase allows the formulation of the material rigidity. However, stresses are calculated separately for each of these phases and damage  The description of behaviour of phases requires the knowledge of intrinsic composite elastic properties in the three directions of space (Young's modulus, shear modulus, Poisson's ratio). These characteristics are obtained through tensile, compressive, and bending tests on specific laminates: • tension and compression tests on an unidirectional composite with fibres at 0, 45 and 90Њ; • tension tests on laminated 0 2 90 6 S materials; and • tension tests on unidirectional and laminated notched materials, etc.
Elastic data of each phase [9,10] are obtained from the intrinsic elastic properties of the constituents and fibre volume fraction a f (Eq. (2)). The data obtained may be reduced with the hypothesis of isotropic transverse behaviour. After the elastic phase, damage inside the material is related to a drop of fibres and/or matrix rigidity. The two phases possess new moduli (C d ) which can be found linearly during loading from elastic moduli (C 0 ) and scalar damage parameter (d) according to C d C 0 1 Ϫ d (Fig. 3). The damage parameter (d) can be expressed as the sum of volumetric damage (d v ) due to a volumetric equivalent strain measure and shear damage (d s ) due to a shear equivalent strain measure. The evolution of damages d v and d s is linear between the threshold values of strain (e i , e l and e u ) and asymptotically reaches a value of one.

Homogeneous model
As opposed to the previous model, this modelling [4,5] proposes some homogenization of the elementary layer. In addition, it considers that composite behaviour in direction fibres is the same in tension and compression.
Damage is taken into account by two scalar variables (d and d H ). The latter expresses displayed experimental phenomena. The parameter d quantifies the damage which comes from the dissociation of fibres and matrix whereas the parameter d H is related to the damage due to the microcracking of matrix parallel to the direction of fibres (Fig.  4). That is the reason the parameters d and d H are applied, respectively, to shear and transverse modulus.
However, tensile and compressive transverse responses must be distinguished. Indeed, micro-cracks grow when the composite ply is under tensile transverse loading, but they close up in the other case. Moreover, modelling also includes the rupture of fibres. Nevertheless, this rupture rate is determined by considering the longitudinal strain limit because rupture of representative layers is brittle in most cases.
So the assigned moduli are calculated from the initial shear modulus G 0 12 and initial transverse modulus E 0 22 according to: Expression of the elementary ply behaviour is given by: The non-damaging thresholds Y and Y H for the elementary layer are given by Eq. (5). The multiplication factor b from Eq. (5) allows the representation of coupling between shear and transverse effects. Yt sup Moreover, laws of the damage evolution are chosen as being linear: Y H S and Y R have been introduced to limit the evolution of non-damaging thresholds Y and Y H ; in time. Indeed, ruptures can be reached even before d and d H reach "one" value.
Elastic characteristics (Young's modulus, shear modulus, Poisson's ratio) and parameters and b result from four experimental characterization tests (Table 1).

Validations of multi-layered multi-material finite element
In order to validate the developments achieved concerning this finite element, an experimental campaign has been followed [3,11]. It concerns the dynamic bending of aluminium/composite plates [12] and axial compressive tests on steel/composite tubes [13]. In this part, the summary of results demonstrates the necessity of such an element.

Three-point dynamic bending
The tested plates 300 × 200 mm 2 are a combination of aluminium (Al 2017TA) and glass E/epoxy composite (30% of fibres). The characteristics of the materials issued from experimental tests are presented in Table 2.
The stacking of plies appears under three shapes: • A mass of 12 kg (trapezoidal shape) impacts these plates under different velocities (4, 6.3 and 6.7 m/s) inside a vertical testing bench (Fig. 5). Experimental results include absorbed energy, applied loads, and maximum deflection determined from the displacement and accelerator sensors.
The behaviour of this multi-layered multi-material is numerically reproduced using element 131 (including biphase model and homogeneous model) and the superposition of two shell elements (the previous solution before the creation of element 131). The adopted meshing, boundary conditions, etc. are summarized in Fig. 6. Fig. 7 shows two major points. First, the numerical results derived from new element 131 come closer to experimental results because the real multi-layered multi-material behaviour is more suitably taken into account. Indeed, the discrepancy between results of the new finite element and superposed shells shows that the last one does not provide a sufficient solution in order to describe the right material rigidity. Secondly, some experimental facts such as delamination could also not be described adequately. Indeed, results for plates marked NS0/90 have suffered from ply detachment because of the crash. Some studies are in progress to predict this type of degradation.

Axial compression of square tubes
The tested tubes 250 × 60 × 60 mm 3 are a combination of one layer of steel (E24) reinforced by two, three, or four layer of glass E/epoxy composite (fibres are at 0Њ). The section and thickness are constant from the beginning to the end of tubes. No defaults are introduced in the creation of these tubes in order to initiate any effects. Table 3 gives the characteristics for both material types.
A mass (shape of a parallelepiped) of 264.4 kg crashes against tubes from a height of 2 m providing a velocity of 6.26 m/s. The testing bench and the sensors are the same as previously noted. Therefore, the experimental results of interest are maximum crushing displacements, maximum forces, time elapsed until maximum crushing, etc. Fig. 8 summarizes the numerical modelling approach. As for the three-point dynamic bending test, a quarter of tube is modelled due to the symmetry conditions.
The results from the numerical simulations (Table 4) allow the same conclusions as for the three-point dynamic bending test. The values of maximum displacement or of the applied mean force generate an error of no more than twelve percent compared with experimental values.

Conclusions on presentation of element 131
Experimental studies and numerical validations permit the conclusion that the development of such a finite element is necessary to describe rapid dynamic multi-layered multimaterial behaviour. Indeed, results generally show errors of less than twelve percent, and this forms an accurate validation. However, some new problems have appeared during

New elastic plastic damaging behaviour for homogeneous ply
In some cases of stacking, loading, etc. the ductility of resin can have a strong effect on composite behaviour. For example, if experimental conditions allow a weak loading of fibres, the behaviour of the resin "will be prevalent". Moreover, some studies on strain rate dependency are carried out in parallel. The first conclusions point to the high sensitivity of strain rate in previous cases where ductility of resin is dominant. It motivates the search for a better description of composite ply behaviour before the introduction of strain rate sensitivity.
It has also been decided to introduce the asymmetrical behaviour frequently noticed for tensile and compressive loading in the direction of fibres. Indeed, this could provide a difference in longitudinal Young's modulus because of fibre micro-buckling, micro-defaults at matrix/fibres interface, and fibres misalignments.
The integration of plasticity is proposed in the following parts as well as an account of asymmetrical tension/ compression behaviour.

Plasticity integration for homogeneous model of element 131
Studies concerning the plasticity of the elementary ply still draw inspiration from the modelling developed by Ladevèze [4,5]. The elastic damage model is identical to the model previously described. However, here the total strain is composed as a sum of elastic strain and plastic strain: e total e elastic ϩ e plastic 7 Elastic domain is described owing to a criteria function f chosen under a general form of Von Mises criteria, which is related to anisotropic materials. This function depends a priori on principal stresses. However, high rigidity of fibres locks plastic flow in their direction. Therefore, criteria function depends only on transverse and shear stresses. Effective stresses permit an agreement to plasticity/damage coupling: Plastic strain evolution is described in accordance with an isotropic hardening law (Eq. (9)). This choice is arbitrary but it allows the evaluation of resin ductility at the first approximation for dynamic studies. Moreover, it permits the simplifying of experimental characterization for newly introduced parameters.
In consequence, criteria function is: b, m, a 2 , and R 0 parameters are characteristics depending on materials and evaluated experimentally. The choice of isotropic hardening law implies that it is not necessary to add extra experimental tests. Even R 0 , b, and m are found from the existing test on laminate [^45] 2S , which is a "quasi-pure" shear case. So, from each load/unload cycles, R and p (by integration) can be found simply (Eq. (11)). Then linear interpolation of the curve log R versus log p gives the necessary coefficients of the hardening law.
As for weighting between transverse and shear stresses (a 2 ), the different load and unload values from the cyclic tensile test on the laminate [ ϩ 45] 8 allow its calculation:

Asymmetrical tension/compression behaviour for homogeneous model of element 131
During compressive loading, laminates can suffer from non-linear rigidity losses. This is generally a consequence of micro-buckling of fibres. Indeed, the misalignment of fibres or the micro-defaults at fibres/matrix interface can strongly influence the emergence of micro-buckling. Therefore, compressive behaviour must be distinguished from tensile behaviour. Drop in Young's modulus-or secant modulus because any damage appears in this direction-is a linear function of compressive stress [14]. The linear coefficient is called g. A similar approach has been introduced using the longitudinal strain: E compression 11 Ϸ E sec ant E initial compression 11 1 ϩ gE sec ant e 11 13 Therefore, Young's modulus is given by Eq. (14) if during calculation a compressive load is detected. Finally, the fact that the strain of rupture could be different in both cases of loading is also taken into account. Consequently, a new parameter is added to measure this compressive rupture strain rate. Note: Table 5 gives the materials properties derived from the new composite experimental characterization. The latter are used for all following numerical simulations.

Tensile test on laminate [^45] 2S
The tested composite material is a glass E/epoxy with 55% of fibres. This sample test piece has been chosen because plastic flow is predominant due to this type of resin. Fig. 9 effectively shows the necessity of taking into account plastic phenomena. An error appears in more than fifty percent of cases if there is no partition of total strain. Of course, this is a "limit" case because in practice, structural parts that use this laminate under the same loading conditions are rarely the same as this case. However, this suggests that plastic phenomena cannot be set aside and that the elastic/plastic damaging algorithm is advised in order to avoid problems of this type.

Dynamic bending on laminate [90 2 0 2 ] S
This experimental test shows the advantage of taking into account a different behaviour in tension and compression. Indeed, the plastic phenomena are not prevalent for this plate configuration. The composite used possesses a Young's modulus and a limit rupture strain that differs by about 40% for the tensile and compressive values. Therefore, the previous modelling overestimates the rigidity of the plies that are in compression. Damage grows prematurely and implies that the composite plate is ruined. The rupture deflection is less important than the case of new modelling (Fig. 10). Nevertheless, the rupture is not correct in both cases. Indeed, this type of test permits the appearance of delamination between lower plies. Micro-cracking appears parallel to fibres inside lower layers. It causes a change of the composite rigidity. Actually, the finite element could not describe this effect but studies are in progress.

Axial compressive test on steel/composite tube
For this last example, the tube is composed of one layer of steel (E24) and two layers of composite at 0Њ (glass E/epoxy with 55% of fibres). The numerical simulation conditions are identical to before: the dimensions of the tube are 200 × 60 × 60 mm 3 : The same experimental testing bench (mass 264.4 kg, impact velocity 6.26 m/s) and the same measurement system are still used. Both the models-elastic damaging and elastic plastic damaging-have been used.
For this type of test, both integrated phenomena appear. The loading allows the resin ductility to express itself, and   part of the specimen is in tension and part of it is in compression along the direction of fibres, during loading.
The crushing displacement (Fig. 11) becomes higher than the one predicted by the previous model. Indeed, kinetic energy is dissipated less quickly because of total strain partition (Fig. 12). The damage is calculated here not with total strains but with "real" elastic strains. As the contour of effective plastic strain shows (Fig. 13), plastic flow becomes important since the composite layers have practically plasticized p Ϸ 6% at the end of calculation).

Conclusions and future prospects
The field of transportation is showing an increasing trend in the usage of numerical simulations. Indeed, the prediction of structures behaviour represents one of the major concerns in order to satisfy the increasingly severe safety criteria. Therefore, the use of new materials, such as multi-layered multi-material, needs a particular numerical description.
The finite element displayed in this paper agrees with the requirements related to the study and the utilization of this type of material. The stacking of varied constituents is taken into account owing to the ply concept. These plies possess an accurate behaviour law in accordance with the nature of material used: for the metallic plies, an elastic plastic law that can integrate the strain rates effects is available. In the other case, i.e. for the unidirectional composite plies, it is possible to obtain an elastic damage law. Most of the numerical simulations show that there is a close connection between numerical and experimental results. This fact confirms the necessity of the development of such a finite element.
This paper also shows the improvements that have been made to the behaviour law of homogeneous composite ply. These improvements are composed, on the one hand, of the integration of the plastic phenomena essentially due to the resin; and on the other hand, of the taking into account of the asymmetrical behaviour in the fibre direction for the tensile and compressive responses. Thus, this element includes a large panel of behaviours in agreement with the nature of composite material: more or less brittle and ductile behaviours and fibres action which are more or less important. The conclusions emanating from experimental studies and numerical validations point out that such an algorithm is a precious precaution for better accuracy, control, etc. of the numerical simulations.
The future prospects issued from the studies purposed in this paper actually concern the strain rates effects and the delamination for unidirectional composite materials. These two phenomena seem to play a very important role in the behaviour of materials during dynamic impacts. Consequently, studies are in progress on these two subjects.