# Multi-layered multi-material finite element for crashworthiness studies

Daniel Coutellier, Patrick Rozycki

► To cite this version:

Daniel Coutellier, Patrick Rozycki. Multi-layered multi-material finite element for crashworthiness studies. Composites Part A: Applied Science and Manufacturing, 2000, 31 (8), pp.841-851. 10.1016/S1359-835X(00)00022-1. hal-01007092

HAL Id: hal-01007092

https://hal.science/hal-01007092

Submitted on 1 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Multi-layered multi-material finite element
for crashworthiness studies

D. Coutellier<sup>a</sup>, P. Rozycki<sup>b</sup>

<sup>a</sup>Lamih-Genie-Mecanique UMR CNRS 8530, Université de Valenciennes et du Hainaut Cambrésis, 59313

Valenciennes Cedex 9, France

<sup>b</sup>Pam System International, 20, rue Saarinen—Silic 303, 94588 Rungis Cedex, France

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.

*Keywords:* : Layered structures; Multi-material; Impact behaviour; Finite element analysis

# 1. 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 confor-
mity, component behaviour prediction, etc. numerical
simulations are being intensively developed, principally
in the domains of transportation. So to foresee the beha-
viour of such metallic/composite materials, it seems para-
mount 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 beha-
viours of fibres and matrix are distinguished separately
(bi-phase model [6–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 vari-
able 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 homo-
geneous 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 beha-
viour. With a view to integration of the strain rates depen-
dence on composite behaviour, a plasticity algorithm and an
asymmetric behaviour in tension/compression have been
added to the element 131 for homogeneous model com-
posite ply.

This paper starts with a presentation of the multi-layered
multi-material finite element. Then it continues with a theo-
retical part concerning plastic phenomena of resin and
<figure><img src="image_2.png" /><figcaption>Fig. 1. Multi-linear stress-strain law for metallic plies.</figcaption></figure>

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 three-point 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.

## 2. 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.

### 2.1. 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:

$$ \sigma \quad \sigma_y = 1 + \left( \frac{\dot{\epsilon}}{D} \right)^{1/p} \qquad (1) $$

where $\sigma_y$ is the yield stress, $\dot{\epsilon}$ the strain rate and $p$ and $D$ are parameters.

### 2.2. 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.

#### 2.2.1. 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

<figure><img src="image_3.png" /><figcaption>Fig. 2. Bi-phase composite model.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Determination of bi-phase model damage.</figcaption></figure>

can be propagated independently according to the selected criteria.

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, 90, 6, 15] 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 $\alpha^f$ (Eq. (2)). The data obtained may be reduced with the hypothesis of isotropic transverse behaviour.

$$E_{11}^f \quad \alpha^f \times E_{\text{true}}^f$$

(where $E_{\text{true}}^f$ is “real” Young's modulus of fibres)

$$E_{11}^m \quad E_{11}^{\text{UD}} - E_{11}^f; \quad E_{22}^m \quad \frac{E_{22}^{\text{UD}}}{1 + (\nu_{12}^{\text{UD}})^2 \frac{E_{22}^{\text{UD}} E_{11}^{\text{UD}}}{E_{11}^{\text{f}}(E_{11}^{\text{UD}} - E_{11}^{\text{f}})}}$$

$$\nu_{12}^m \quad \nu_{12}^{\text{UD}}; \quad \nu_{21}^m \quad \frac{\nu_{21}^{\text{UD}}}{1 - E_{11}^f \frac{1 - \nu_{12}^{\text{UD}} \nu_{21}^{\text{UD}}}{E_{11}^{\text{UD}}}} \quad G_{12}^m \quad G_{12}^{\text{UD}} \quad (2)$$

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 ($\epsilon_i, \epsilon_1$ and $\epsilon_u$) and asymptotically reaches a value of one.

### 2.2.2. 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'$). 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'$ is related to the damage due to the micro-cracking of matrix parallel to the direction of fibres (Fig. 4). That is the reason the parameters $d$ and $d'$ 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
<figure><img src="image_5.png" /><figcaption>Fig. 4. Damage modes inside elementary ply for the homogeneous model.</figcaption></figure>

shear modulus $G_{12}^0$ and initial transverse modulus $E_{22}^0$ according to:

$$ \begin{cases} E_{22} & E_{22}^0(1-d') & \text{if } \sigma_{22} > 0 \\ E_{22} & E_{22}^0 & \text{if } \sigma_{22} \le 0 \end{cases} \quad G_{12} \quad G_{12}^0(1-d) \quad (3) $$

Expression of the elementary ply behaviour is given by:

$$ \left\{ \begin{matrix} \epsilon_{11} \\ \epsilon_{22} \\ 2\epsilon_{12} \end{matrix} \right\} = \begin{bmatrix} \frac{1}{E_{11}^0} & -\frac{\nu_{12}^0}{E_{11}^0} & 0 \\ -\frac{\nu_{12}^0}{E_{11}^0} & \frac{1}{E_{22}^0(1-d')}\frac{\langle\sigma_{22}\rangle_+}{|\sigma_{22}|} + \frac{1}{E_{22}^0}\frac{\langle\sigma_{22}\rangle_-}{|\sigma_{22}|} & 0 \\ 0 & 0 & \frac{1}{G_{12}^0(1-d)} \end{bmatrix} \times \left\{ \begin{matrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{matrix} \right\} \quad (4) $$

where

$$ \langle a \rangle_+ = \begin{cases} a & \text{if } a > 0 \\ 0 & \text{if } a < 0 \end{cases} \quad \text{and} \quad \langle a \rangle_- = \begin{cases} 0 & \text{if } a > 0 \\ a & \text{if } a < 0 \end{cases} $$

The non-damaging thresholds ($\bar{Y}$ and $\bar{Y}'$) 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.

$$ \begin{align} \bar{Y}(t) &= \sup_{\tau \le t} \sqrt{\frac{1}{2} \frac{\sigma_{12}^2}{G_{12}^0(1-d)^2} + b \frac{\langle\sigma_{22}^2\rangle_+}{E_{22}^0(1-d')^2}} \tag{5} \\ \bar{Y}'(t) &= \sup_{\tau \le t} \sqrt{\frac{1}{2} \frac{\langle\sigma_{22}^2\rangle_+}{E_{22}^0(1-d')^2}} \end{align} $$

Moreover, laws of the damage evolution are chosen as being linear:

$$ \begin{align} d \quad &\frac{\langle \bar{Y} - Y_O \rangle_+}{Y_C} \text{ if } d < 1 \text{ and } \bar{Y}' < Y'_S \text{ and } \bar{Y} < Y_R \notag \\ \text{else } d \quad &1 \tag{6} \\ d' \quad &\frac{\langle \bar{Y}' - Y'_O \rangle_+}{Y'_C} \text{ if } d' < 1 \text{ and } \bar{Y}' < Y'_S \text{ and } \bar{Y} < Y_R \notag \\ \text{else } d' \quad &1 \notag \end{align} $$

$Y'_S$ and $Y_R$ have been introduced to limit the evolution of non-damaging thresholds $\bar{Y}$ and $\bar{Y}'$, in time. Indeed, ruptures can be reached even before $d$ and $d'$ reach “one” value.

Elastic characteristics (Young's modulus, shear modulus, Poisson's ratio) and parameters $Y_O$, $Y_C$, $Y'_O$, $Y'_C$, $Y_R$, $Y'_S$ and $b$ result from four experimental characterization tests (Table 1).

## 2.3. Validations of multi-layered multi-material finite element

In order to validate the developments achieved concerning

<table><caption>Table 1
Data parameters corresponding to the appropriate experimental test</caption><thead><tr><th>Rupture tensile<br>test on [0]<sub>8</sub></th><th>Cyclic tensile<br>test on [±45]<sub>2S</sub></th><th>Cyclic tensile<br>test on [±67.5]<sub>2S</sub></th><th>Cyclic tensile<br>test on [+45]<sub>8</sub></th></tr></thead><tbody><tr><td>E<sub>11</sub></td><td>G<sub>12</sub></td><td>E<sub>22</sub></td><td>Check of G<sub>12</sub> and E<sub>22</sub></td></tr><tr><td>ν<sub>12</sub> fracture limit</td><td>Y<sub>C</sub></td><td>Y'<sub>C</sub></td><td>Check of Y<sub>C</sub> and Y'<sub>C</sub></td></tr><tr><td>ε<sub>11</sub></td><td>Y<sub>O</sub></td><td>Y'<sub>O</sub></td><td>Check of Y<sub>O</sub> and Y'<sub>O</sub></td></tr><tr><td></td><td>Y<sub>R</sub></td><td>Y'<sub>S</sub></td><td></td></tr><tr><td></td><td></td><td>b</td><td></td></tr></tbody></table>
<table><caption>Table 2
Data parameters required for the numerical simulations of the dynamic bending test</caption><thead><tr><th></th><th colspan="3">Glass E/epoxy (30% of fibres)</th></tr><tr><th></th><th>Aluminium</th><th>Bi-phase model</th><th>Homogeneous model</th></tr></thead><tbody><tr><td>Density</td><td>2790 kg/m<sup>3</sup></td><td>1620 kg/m<sup>3</sup></td><td>1620 kg/m<sup>3</sup></td></tr><tr><td>Young's modulus</td><td>E = 72748 MPa</td><td>E<sub>11</sub> = 553 MPa<br>E<sub>22</sub> = 1653 MPa</td><td>E<sub>11</sub> = 22753 MPa<br>E<sub>22</sub> = 1955 MPa</td></tr><tr><td>Shear modulus</td><td></td><td>G<sub>12</sub> = 6487 MPa</td><td>G<sub>12</sub> = 6487 MPa</td></tr><tr><td>Poisson's ratio</td><td>ν = 0.34</td><td>ν<sub>12</sub> = 0.23</td><td>ν<sub>12</sub> = 0.23</td></tr><tr><td>Yield stress</td><td>σ<sub>y</sub> = 227 MPa</td><td>—</td><td>—</td></tr><tr><td>Multi-linear law<br>for metallic ply or<br>damage<br>parameters for<br>composite plies</td><td></td><td>Fibres</td><td>Matrix</td></tr><tr><td></td><td>E<sub>1</sub> = 25090 MPa<br>σ<sub>1</sub> = 355 MPa<br>E<sub>2</sub> = 640 MPa<br>σ<sub>2</sub> = 355 MPa<br>E<sub>3</sub> = 181 MPa<br>σ<sub>3</sub> = 416 MPa<br>(No strain rate law input)</td><td>ε<sub>f</sub><sup>f</sup> = 1.90%<br>ε<sub>f</sub><sup>i</sup> = 1.95%<br>ε<sub>u</sub><sup>f</sup> = 2.00%<br>d<sub>f</sub><sup>i</sup> = 0.001<br>d<sub>u</sub><sup>f</sup> = 0.999</td><td>ε<sub>si</sub><sup>m</sup> = 1.349%<br>ε<sub>si</sub><sup>m</sup> = 1.354%<br>ε<sub>su</sub><sup>m</sup> = 1.420%<br>d<sub>si</sub><sup>m</sup> = 0.001<br>d<sub>su</sub><sup>m</sup> = 0.210</td><td>Y<sub>c</sub> = 1.924 MPa<sup>1/2</sup><br>Y<sub>O</sub> = 0.189 MPa<sup>1/2</sup><br>Y<sub>c</sub><sup>f</sup> = 1.935 MPa<sup>1/2</sup><br>Y<sub>c</sub><sup>f</sup> = 0.375 MPa<sup>1/2</sup><br>b = 1.672<br>Y<sub>s</sub><sup>f</sup> = 0.181 MPa<sup>1/2</sup><br>Y<sub>R</sub> = 0.618 MPa<sup>1/2</sup><br>ε<sub>11</sub><sup>R</sup> = 2.1%</td></tr></tbody></table>

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.

### 2.3.1. Three-point dynamic bending

The tested plates (300 × 200 mm²) 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:

* [Al 0 Al 0 Al]—mark S—;

* [Al 90 0 Al 90 0 Al]—mark NS90/0—; and

* [Al 0 90 Al 0 90 Al]—mark NS0/90—.

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 bi-phase 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

<figure><img src="image_6.png" /><figcaption>Fig. 5. Experimental testing bench.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Quarter meshed plate of the three-point dynamic bending test.</figcaption></figure>

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 experi-
mental 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.

2.3.2. Axial compression of square tubes

The tested tubes (250 × 60 × 60 mm³) 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.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Maximal deflection results for the different cases of study.</figcaption></figure>
<table><caption>Table 3
Date parameters required for the numerical simulations on the axial compressive test</caption><thead><tr><th></th><th colspan="3">Glass E/epoxy (30% of fibres)</th></tr><tr><th></th><th>Steel</th><th>Bi-phase model</th><th>Homogeneous model</th></tr></thead><tbody><tr><td>Density</td><td>7800 kg/m<sup>3</sup></td><td>1620 kg/m<sup>3</sup></td><td>1620 kg/m<sup>3</sup></td></tr><tr><td>Young's modulus</td><td>E = 200000 MPa</td><td>E<sub>11</sub> = 553 MPa<br>E<sub>22</sub> = 1653 MPa</td><td>E<sub>11</sub> = 22753 MPa<br>E<sub>22</sub> = 1955 MPa</td></tr><tr><td>Shear modulus</td><td></td><td>G<sub>12</sub> = 6487 MPa</td><td>G<sub>12</sub> = 6487 MPa</td></tr><tr><td>Poisson's ratio</td><td>ν = 0.3</td><td>ν<sub>12</sub> = 0.23</td><td>ν<sub>12</sub> = 0.23</td></tr><tr><td>Yield stress</td><td>σ<sub>y</sub> = 250 MPa</td><td>-</td><td>-</td></tr><tr><td>Multi-linear law<br>for metallic ply or<br>damage<br>parameters for<br>composite plies</td><td>E<sub>1</sub> = 49560 MPa<br>σ<sub>1</sub> = 290 MPa<br>E<sub>2</sub> = 12390 MPa<br>σ<sub>2</sub> = 310 MPa<br>E<sub>3</sub> = 4375 MPa<br>σ<sub>3</sub> = 345 MPa<br>E<sub>4</sub> = 2000 MPa<br>σ<sub>4</sub> = 365 MPa<br>E<sub>5</sub> = 437.5 MPa<br>σ<sub>5</sub> = 785 MPa<br>Strain rate law:<br>Cowper-<br>Symonds<br>D = 802 s<sup>-1</sup><br>p = 3.585</td><td>Fibres<br>ε<sub>f</sub><sup>l</sup> = 1.90%<br>ε<sub>f</sub><sup>l</sup> = 1.95%<br>ε<sub>f</sub><sup>u</sup> = 2.00%<br>d<sub>f</sub><sup>l</sup> = 0.001<br>d<sub>f</sub><sup>u</sup> = 0.999</td><td>Matrix<br>ε<sub>si</sub><sup>m</sup> = 1.349% ε<sub>vi</sub><sup>m</sup> = 0<br>ε<sub>si</sub><sup>m</sup> = 1.354% ε<sub>vi</sub><sup>m</sup> = 0<br>ε<sub>su</sub><sup>m</sup> = 1.420% ε<sub>vu</sub><sup>m</sup> = 0<br>d<sub>si</sub><sup>m</sup> = 0.001 d<sub>vi</sub><sup>m</sup> = 0<br>d<sub>su</sub><sup>m</sup> = 0.210 d<sub>vu</sub><sup>m</sup> = 0</td><td>Y<sub>C</sub> = 1.924 MPa<sup>1/2</sup><br>Y<sub>O</sub> = 0.189 MPa<sup>1/2</sup><br>Y'<sub>C</sub> = 1.935 MPa<sup>1/2</sup><br>Y'<sub>C</sub> = 0.375 MPa<sup>1/2</sup><br>b = 1.672<br>Y'<sub>S</sub> = 0.181 MPa<sup>1/2</sup><br>Y'<sub>R</sub> = 0.618 MPa<sup>1/2</sup><br>ε<sub>11</sub><sup>R</sup> = 2.1%</td></tr></tbody></table>

## 2.4. 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 multi-material behaviour. Indeed, results generally show errors of less than twelve percent, and this forms an accurate validation. However, some new problems have appeared during

<figure><img src="image_9.png" /><figcaption>Fig. 8. Quarter meshed box with boundary conditions.</figcaption></figure>
<table><caption>Table 4
Error percentage between experimental and numerical results for square tubes</caption><thead><tr><th rowspan="2">Number of shell elements for<br>numerical simulation: 600</th><th colspan="2">Homogeneous model</th><th colspan="2">Bi-phase model</th><th rowspan="2">Experiment</th></tr><tr><th>Numerical<br>values</th><th>Error with<br>respect to<br>experiment (%)</th><th>Numerical<br>values</th><th>Error with<br>respect to<br>experiment (%)</th></tr></thead><tbody><tr><td>Maximum displacement (mm)</td><td>157.23</td><td>6.98</td><td>157.11</td><td>7.05</td><td>169.03</td></tr><tr><td>Mean force (N)</td><td>3018.60</td><td>4.53</td><td>28362.53</td><td>10.29</td><td>31615.8</td></tr><tr><td>Maximum crushing time (ms)</td><td>50.00</td><td>0.40</td><td>45.70</td><td>8.96</td><td>50.2</td></tr></tbody></table>

these experimental tests. Matrix ductility, compressive loading non-integrated for homogeneous model, etc. need to be implemented in order to better control numerical simulations and their accuracy, especially for entire composite sample.

### 3. 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.

#### 3.1. 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:

$$ \epsilon^{\text{total}} = \epsilon^{\text{elastic}} + \epsilon^{\text{plastic}} \quad (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:

$$ \tilde{\sigma}_{11} \quad \sigma_{11} \quad \tilde{\sigma}_{22} \quad \frac{\langle \sigma_{22} \rangle_{+}}{1 - d'} + \langle \sigma_{22} \rangle_{-} \quad \tilde{\sigma}_{12} \quad \frac{\sigma_{12}}{1 - d} \quad (8) $$

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.

$$ R \quad \beta p^m + R_0 \quad \begin{cases} \beta, m: \text{ material constants} \\ p: \text{ effective plastic strain} \\ R_0: \text{ initial yield stress} \end{cases} \quad (9) $$

In consequence, criteria function is:

$$ f(\tilde{\sigma}, R) \quad \sqrt{\tilde{\sigma}_{12}^2 + a^2 \tilde{\sigma}_{22}^2 - R} \quad (10) $$

$\beta, 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$, $\beta$, 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.

$$ R \quad \frac{\sigma_{12}}{(1 - d)} - R_0 \quad \dot{p} \quad \dot{\epsilon}_{12}^p (1 - d) \quad (11) $$

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:

$$ a^2 \quad \frac{\dot{\epsilon}_{22}^p (1 - d')^2}{2 \dot{\epsilon}_{12}^p (1 - d)^2} \quad (12) $$

#### 3.2. Asymmetrical tension/compression behaviour for homogeneous model of element 131

During compressive loading, laminates can suffer from
<table><caption>Table 5
Data parameters required for the numerical simulations of the dynamic bending test</caption><thead><tr><th>Elastic parameters</th><th>Damage parameters</th><th>Plastic parameters</th><th>Compressive parameters</th></tr></thead><tbody><tr><td>ρ = 1970 kg/m³</td><td>Y<sub>C</sub> = 2.128 MPa<sup>1/2</sup></td><td>R<sub>0</sub> = 14.62 MPa</td><td>E<sub>11</sub> = 27260 MPa</td></tr><tr><td>E<sub>11</sub> = 43340 MPa</td><td>Y<sub>O</sub> = 0.169 MPa<sup>1/2</sup></td><td>β = 436</td><td>γ = 6.272 × 10<sup>-6</sup> MPa<sup>-1</sup></td></tr><tr><td>E<sub>22</sub> = 11180 MPa</td><td>Y'<sub>C</sub> = 2.922 MPa<sup>1/2</sup></td><td>m = 0.3247</td><td>ε<sup>R</sup><sub>11</sub> = 3.0%</td></tr><tr><td>G<sub>12</sub> = 3959 MPa</td><td>Y'<sub>O</sub> = 0.169 MPa<sup>1/2</sup></td><td>a<sup>2</sup> = 0.3633</td><td></td></tr><tr><td>ν<sub>12</sub> = 0.31</td><td>b = 1.715</td><td></td><td></td></tr><tr><td>ε<sup>R</sup><sub>11</sub> = 1.5%</td><td>Y'<sub>S</sub> = 0.380 MPa<sup>1/2</sup></td><td></td><td></td></tr><tr><td></td><td>Y<sub>R</sub> = 1.368 MPa<sup>1/2</sup></td><td></td><td></td></tr></tbody></table>

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 $γ$. A similar approach has been introduced using the longitudinal strain:

$$E_{11}^{\text{compression}} \approx E^{\text{sec ant}} \quad E_{11}^{\text{initial compression}}(1 + \gamma E^{\text{sec ant}} \epsilon_{11}) \quad (13)$$

Therefore, Young's modulus is given by Eq. (14) if during calculation a compressive load is detected.

$$E_{11}^{\text{compression}} \approx E^{\text{sec ant}} \frac{E_{11}^{\text{initial compression}}}{1 - \gamma E_{11}^{\text{initial compression}} \epsilon_{11}} \quad (14)$$

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.

### 3.3. Numerical validations of leading development concerning element 131

This part starts with some examples about the integration

of the newness added to the homogeneous model. They are related to:

* the tensile test on laminate [±45]₂s;

* the three-point dynamic bending test on a plate entirely in the composite;

* the previous axial compressive test on the steel/composite tube.

*Note:* Table 5 gives the materials properties derived from the new composite experimental characterization. The latter are used for all following numerical simulations.

#### 3.3.1. Tensile test on laminate [±45]₂s

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.

#### 3.3.2. Dynamic bending on laminate [90<sub>2</sub> 0<sub>2</sub>]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

<figure><img src="image_10.png" /><figcaption>Fig. 9. Tensile test on a [±45]₂s glass E/epoxy (60% of fibres).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Applied force versus deflection for dynamic bending plate.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Crushing displacement curves.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Kinetic, internal and total energies curves.</figcaption></figure>

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.

### 3.3.3. 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³. 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

<figure><img src="image_14.png" /><figcaption>Fig. 13. Contour of plastic effective strain.</figcaption></figure>
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 \approx 6\%$ at the end of calculation).

## 4. 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.

## Acknowledgements

The authors would like to thank Pam System International and Engineering Systems International companies for the trust that they have shown in us since the beginning of our collaboration, particularly M. Haug E. and M. Ni X. Moreover, they would like to express their deepest gratitude to the “Conseil Régional du Nord-Pas-de-Calais” as well as to the CNRS.

## References

[1] Hanefi E. Etude du comportement au crash et de l'absorption d'énergie dans les structures en matériaux légers. Doctoral dissertation, INSA, Lyon, 1996.

[2] Wang X. Comportement à l'écrasement de structures tubulaires en multi-matériaux. Doctoral dissertation, l'École Centrale, Lyon, 1991.

[3] Gauthier C. Contribution à la modélisation en crash des structures stratifiées métal/composite. Développement d'un élément de coque multicouches multi-matériaux. Application à des essais dynamiques. Doctoral dissertation, University of Valenciennes and Hainaut-Cambrésis, 1996.

[4] Ladeveze P, Le Dantec E. Damage modeling of the elementary ply for laminated composites. Composites Science and Technology 1992;43:257-67.

[5] Ladeveze P. Sur la mécanique de l'endommagement des composites. In: Pluralis, editor. J.N.C.5. Paris, 1986, p. 667-83.

[6] Haug E, De Rouvray A. Crash response of composite structures, Structural crashworthiness and failure. Amsterdam: Elsevier, 1992 (chap. 7; p. 237-94).

[7] Haug E, Fort O, Trameçon A, Watanabe M, Nakada I. Numerical crashworthiness simulation of automotive structures and components made of continuous fiber reinforced composite and sandwich assemblies. S.A.E. International Congress & Exhibition, Detroit, 1991, p. 245-58.

[8] De Rouvray A, Haug E. Failure of brittle and composite materials by numerical methods. In: Wierzbicki T, editor. Structural failure, New York: Wiley, 1989. p. 237-92 (chap. 7).

[9] Pam Crash™—ESI Group Software Product Company. Theory manual, 1992.

[10] Pam Crash™—ESI Group Software Product Company. Solver notes manual, 1998.

[11] Coutellier D, Gauthier C, Ravalard Y, Ni X, Haug E. Un élément multi-matériaux pour l'étude du comportement en crash des structures stratifiées: synthèse et perspectives. In Hermès, editor. Calcul des structures GIENS'97. Revue européenne des Eléments finis, 7 (1-3), 1998, p. 177-92.

[12] Coutellier D, Gauthier C, Ravalard Y, Ni X, Haug E. Modélisation en flexion des multicouches multi-matériaux en dynamique rapide. Journées Nationales sur les Composites Paris 1996;10:969-80.

[13] Coutellier D, Gauthier C, Ravalard Y, Ni X, Haug E. Simulation numérique du comportement des multicouches multi-matériaux sous impact. Revue des Composites et des Matériaux Avancés 1996;6(2):199-216.

[14] Le Dantec E. Contribution à la modélisation du comportement mécanique de composites stratifiés. Doctoral dissertation, University of Pierre et Marie Curie, Paris 6, 1989.