# A new modelling of plasticity coupled with the damage and identification for carbon fibre composite laminates

Ahmed Boutaous, Bernard Peseux, Laurent Gornet, Abdelkader Bélaidi

► To cite this version:

Ahmed Boutaous, Bernard Peseux, Laurent Gornet, Abdelkader Bélaidi. A new modelling of plasticity coupled with the damage and identification for carbon fibre composite laminates. Composite Structures, 2006, 74 (1), pp.1-9. 10.1016/j.compstruct.2005.11.004. hal-01004914

HAL Id: hal-01004914

https://hal.science/hal-01004914

Submitted on 16 Mar 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.

Public Domain
A new modeling of plasticity coupled with the damage
and identification for carbon fibre composite laminates

Ahmed Boutaous<sup>a</sup>, Bernard Peseux<sup>b</sup>, Laurent Gornet<sup>b</sup>, Abdelkader Bélaidi<sup>c</sup>

<sup>a</sup> Faculty of Science, Department of Physics, U.S.T.O BP 1505, Oran El Mnaouar, Algeria

<sup>b</sup> Institute of Research in Civil and Mechanical Engineering, UMR CNRS, 6183 Nantes, France

<sup>c</sup> Department of Physics, E.N.S.E.T BP 1425, Oran, Algeria

The research described herein was aimed at gaining a better understanding of the behavior up to the fracture of laminate composite structures. An elastoplastic damage model was developed to describe the behavior of these structures in the greatest detail possible. The damage is modeled on each elementary constituent. Then a complete model at meso (intermediate) scale is obtained by applying a homogenization method.

This model takes into account the transverse tensile and shear damage of matrix and fibre-matrix interface, as well as the nonlinear behavior of the fibre under axial compression. Plasticity develops only in the matrix and its flow is blocked in the direction of the fibres. It is modeled by an anisotropic yield criterion (of the HILL type) which takes not only the transverse and shear stress into account but which also includes isotropic and nonlinear kinematic hardening.

The model developed was included in the nonlinear finite element “Castem 2000” to allow comparisons between numerical and experimental results. There is a good agreement between the simulated and experimental data. However, the model could advantageously be extended to take the phenomena such as viscoplasticity into account.

Keywords: Composite materials; Laminates; Elementary layer; Fibre-matrix; Damage; Plastic

# 1. Introduction

To have a better reliability of the structures made out of composite materials, research was undertaken to simulate their behavior. One essential characteristic of this behavior is the appearance of various degradations well before the total rupture of the composite.

The authors of Refs. [1,2] developed models in mechanics accounting for degradation within these materials. This degradation was taken into account in the form of damage. D. Gilletta modeled the elementary layer as a whole [1], it takes into account of the damage in transverse tensile and shearing. This model does not seem sufficient with the sight of experimental results, Le Dantec [2] developed a second

model based on the first one. In this model, the damage of the matrix and the interface matrix fibre are taken into account, as well as plasticity with isotropic hardening. Moreover these models are written in an implicit way, i.e. in strain and ignore a certain number of other phenomena existing within this elementary layer like plasticity with nonlinear kinematic hardening and the loss of rigidity of fibres.

It was proved to be necessary to better determine the behavior of the laminated structures, to establish in a computer code using finite element “Castem 2000” developed by the C.E.A, a more refined model of degradation. Thus, not only the damage in shearing and in transverse traction will be taken into account, but also behavior in longitudinal compression. To supplement this model, plasticity existing within the matrix was coupled with the damage which is based on anisotropic criteria of Hill, with isotropic
and kinematic nonlinear hardening. This modeling is the subject of this paper.

## 2. General framework

The model presented is developed within the framework of the irreversible phenomena of thermodynamics, where the thermal phenomena are neglected, under the assumption of the generalized plane strain ($\sigma_{33} = 0$); it is formulated in stress.

The laminated composite materials made up of one-way folds being transversely isotropic compared to the direction of fibres [3], the volumic energy of elastic strain which is taken as potential thermodynamics is given by:

$$2W_e = \frac{\sigma_{11}^2}{E_{11}} + \frac{\sigma_{22}^2}{E_{22}} - \frac{2\nu_{12}\sigma_{11}\sigma_{22}}{E_{11}} + \frac{\sigma_{12}^2}{G_{12}} + \frac{\sigma_{13}^2}{G_{13}} + \frac{\sigma_{23}^2}{G_{23}} \quad (1)$$

with

$$\frac{\nu_{12}}{E_{11}} = \frac{\nu_{21}}{E_{22}}$$

where the elastic law of behavior is the expression of the stress vs strain.

$$\epsilon_{ij}^e = - \frac{\partial W_e}{\partial \sigma_{ij}} \left\{ \begin{array}{l} \epsilon_{11}^e = \frac{1}{E_{11}} (\sigma_{11} - \nu_{12}\sigma_{22}) \\ \epsilon_{22}^e = \frac{1}{E_{22}} (\sigma_{22} - \nu_{21}\sigma_{11}) \\ \epsilon_{12}^e = \frac{\sigma_{12}}{2G_{12}}, \quad \epsilon_{13}^e = \frac{\sigma_{13}}{2G_{13}}, \quad \epsilon_{23}^e = \frac{\sigma_{23}}{2G_{23}} \end{array} \right.$$

By inversion of this law, the strain is expressed according to the stress as follows:

$$\begin{aligned} \sigma_{11} &= C_{11}(\varepsilon_{11} + \nu_{21}\varepsilon_{22}) \\ \sigma_{22} &= C_{22}(\varepsilon_{22} + \nu_{12}\varepsilon_{11}) \\ \sigma_{12} &= G_{12}\gamma_{12}, \quad \sigma_{13} = G_{13}\gamma_{13}, \quad \sigma_{23} = G_{23}\gamma_{23} \end{aligned}$$

with

$$C_{11} = \frac{E_{11}}{1 - \nu_{12}\nu_{21}}; \quad C_{22} = \frac{E_{22}}{1 - \nu_{12}\nu_{21}} \\ \gamma_{12} = 2\varepsilon_{12}; \quad \gamma_{13} = 2\varepsilon_{13}; \quad \gamma_{23} = 2\varepsilon_{23}$$

And the volumic energy of elastic strain (1) as a function of the deformations is written as:

$$2W_e = \sigma_{11}\varepsilon_{11} + \sigma_{22}\varepsilon_{22} + \sigma_{12}\gamma_{12} + \sigma_{13}\gamma_{13} + \sigma_{23}\gamma_{23} \quad (2)$$

where

$$2W_e = C_{11}\varepsilon_{11}^2 + 2\nu_{12}\varepsilon_{11}\varepsilon_{22} + C_{22}\varepsilon_{22}^2 + G_{12}\gamma_{12}^2 + G_{13}\gamma_{13}^2 + G_{23}\gamma_{23}^2$$

By taking into account the relations of symmetry coming from the existence of a thermodynamic potential we, necessarily, have:

$$\frac{C_{11}}{C_{22}} = \frac{E_{11}}{E_{22}} = \frac{\nu_{12}}{\nu_{21}}$$

It follows: $\nu_{21}C_{11} = \nu_{12}C_{22}$.

And according to the definition of $C_{11}$ we notice that:

$$C_{11}(1 - \nu_{12}\nu_{21}) = E_{11} \Rightarrow C_{11} = E_{11} + \nu_{12}\nu_{21} = E_{11} + \nu_{12}^2 C_{22}$$

In this form, the value of the strain $\sigma_{22} = C_{22}(\varepsilon_{22} + \nu_{12}\varepsilon_{11})$ clearly appears.

The elastic deformation energy (1) can thus be now written as:

$$\begin{aligned} 2W_e &= E_{11}\varepsilon_{11}^2 + \nu_{12}^2\varepsilon_{11}^2 + 2\nu_{12}C_{22}\varepsilon_{11}\varepsilon_{22} + C_{22}\varepsilon_{22}^2 \\ &\quad + G_{12}\gamma_{12}^2 + G_{13}\gamma_{13}^2 + G_{23}\gamma_{23}^2 \\ 2W_e &= E_{11}\varepsilon_{11}^2 + C_{22}(\varepsilon_{22} + \nu_{12}\varepsilon_{11})^2 \\ &\quad + G_{12}\gamma_{12}^2 + G_{13}\gamma_{13}^2 + G_{23}\gamma_{23}^2 \end{aligned}$$

Concerning the elastoplastic law of behavior of the elementary layer [4], it will be broken up into elastic strain and plastic stress according to the traditional form:

$$\left\{ \begin{aligned} \varepsilon_{11} &= \frac{\sigma_{11}}{E_{11}} - \frac{\nu_{12}\sigma_{22}}{E_{11}} + \varepsilon_{11}^p \\ \varepsilon_{22} &= \frac{\sigma_{22}}{E_{22}} - \frac{\nu_{12}\sigma_{11}}{E_{22}} + \varepsilon_{22}^p \\ \varepsilon_{12} &= \frac{\sigma_{12}}{2G_{12}} + \varepsilon_{12}^p \end{aligned} \right.$$

The worked out model, describing the undamaged behavior of an elementary layer will now be presented.

## 3. Modeling of the damage of the elementary layer

One seeks to build a model on the level meso based on the characteristics of degradations on the microlevel.

According to Ref. [5], it was observed that the damage appears in the form of cracks located in the matrix and laid out parallel to fibres.

The damage is supposed to occur only in shearing and transverse tensile. Indeed in transverse compression cracks tend to be closed again, and do not create an additional damage.

Moreover in the direction of fibres there is a fall of the longitudinal module in localized compression because of buckling of the fibres.

On the level of the elementary layer, the modules $G_{12}$ and $C_{22}$ are affected by two variables of damage $d_{12}$ and $d_{22}$, and $E_{11}$ is affected by a variable $\xi_{11}$ representing the loss of rigidity of fibres into longitudinal compression and brittle fracture in tensile in the direction of the fibres. Thus, the elastic modulus of the materials can be written as:

$$\left\{ \begin{aligned} E_{11} &= E_{11}^0 (1 - \xi_{11}) \\ C_{22} &= C_{22}^0 (1 - d_{22}) \quad \text{if } \varepsilon_{22} + \nu_{12}\varepsilon_{11} > 0 \\ C_{22} &= C_{22}^0 \quad \text{if } \varepsilon_{22} + \nu_{12}\varepsilon_{11} \le 0 \\ G_{12} &= G_{12}^0 (1 - d_{12}) \end{aligned} \right. \quad (3)$$

where the index “0” indicates the undamaged state of the material, and for lack of experimental information on
other modules. The latter are assumed not prone to damage.

$$v_{12} = v_{12}^{0}; \quad G_{13} = G_{13}^{0}; \quad G_{23} = G_{23}^{0}$$

In such a material, initially transverse isotropic $G_{12}^0 = G_{13}^0$ do not necessarily remain with transverse isotropy $G_{12}^0 \neq G_{13}^0$.

Also the module $E_{22}$ and the Poisson's ratio $v_{21}$ are affected by the damage. However, let $d''$ be the variable affecting $E_{22}$ i.e.

$$E_{22} = E_{22}^{0}(1 - d'') \quad (4)$$

We get:

$$v_{21} = v_{12} \frac{E_{22}}{E_{11}} = v_{12}^{0} \frac{E_{22}^{0}(1 - d'')}{E_{11}^{0}} = v_{21}^{0}(1 - d'') \quad (5)$$

$d_{22}$ and $d''$ are related by the relation existing between $C_{22}$ and $E_{22}$ that allows to clarify $d_{22}$ according to $d''$.

$$\Rightarrow d_{22} = \frac{d''}{1 - v_{12}^{0} v_{21}^{0} (1 - d'')} \quad (6)$$

the relative error obtained by replacing $d''$ by $d_{22}$ is then

$$\frac{d_{22} - d''}{d''} = \frac{v_{12}^{0} v_{21}^{0} (1 - d'')}{1 - v_{12}^{0} v_{21}^{0} (1 - d'')} \quad (7)$$

It is maximum on $[0, 1]$ in $d'' = 0$ and decreases to “0” in $d'' = 1$ the maximum value is equal to:

$$\frac{v_{12}^{0} v_{21}^{0}}{1 - v_{12}^{0} v_{21}^{0}} \quad (8)$$

which remains weak as long as:

$$v_{12}^{0} v_{21}^{0} \ll 1 \quad (9)$$

which is checked for the studied material.

The variables $d''$ and $d_{22}$ are thus appreciably identical.

### 3.1. Homogenization of the elementary layer

The calculation of homogenization is carried out on the modules $E_{22}$ and $G_{12}$. The method of the asymptotic developments for the periodic mediums is used [6,7], indicated by $d_{12}$ and $d_{22}$ the loss of rigidity of the layer [2]:

$$\left\{ \begin{array}{l} \frac{1}{G_{12}(1 - d_{12})} = \frac{v_m}{G_{12}^m(1 - D_{12})} + \frac{2}{ek_{66}(1 - \delta_{12})} \\ \frac{1}{E_{22}(1 - d_{22})} = \frac{v_m}{E_{22}^m(1 - D_{22})} + \frac{1}{ek_{22}(1 - \delta_{22})} \end{array} \right. \quad (10)$$

and

$$v_m + v_f = 1 \quad (11)$$

where

$v_m$: volumic fraction of the matrix,

$v_f$: volumic fraction of fibre,

$e$: thickness of the fold,

$D_{ii}$: damage of the matrix,

$\delta_{ii}$: damage of the interface.

By carrying out a limited development in the neighborhood of “0”, the damage of the layer is then given by:

$$\left\{ \begin{aligned} d_{12} &= \frac{v_m G_{12} D_{12}}{G_{12}^m} + \frac{2 G_{12} \delta_{12}}{ek_{66}} \\ d_{22} &= \frac{v_m E_{22} D_{22}}{E_{22}^m} + \frac{E_{22} \delta_{22}}{ek_{22}} \end{aligned} \right. \quad (12)$$

Since some parameters were difficult to determine, it was decided to obtain a relatively simple model and conform to that developed by [2], to take the damage of the elementary layer as a sum of the damage due to the matrix and the fibre-matrix interface, respectively.

$$d_{12} = D_{12} + \delta_{12} \quad (13)$$

$$d_{22} = D_{22} + \delta_{22} \quad (14)$$

Finally, in the model on the meso scale there are only three internal variables: two variables of damage $d_{12}$ and $d_{22}$, and one variable $\xi_{11}$ translating the loss of rigidity fibres in compression.

### 3.2. Model of damage of the elementary layer

$$\left\{ \begin{aligned} \sigma_{11} &= E_{11}^{0}(1 - \xi_{11})\varepsilon_{11}^{e} + (v_{12}^{0})^{2}C_{22}^{0}(1 - d_{22})\varepsilon_{11}^{e} \\ &\quad + v_{12}^{0}C_{22}^{0}(1 - d_{22})\varepsilon_{22}^{e} \\ \sigma_{22} &= v_{12}^{0}C_{22}^{0}(1 - d_{22})\varepsilon_{11}^{e} + C_{22}^{0}(1 - d_{22})\varepsilon_{22}^{e} \\ \sigma_{12} &= G_{12}^{0}(1 - d_{12})\gamma_{12}^{e} \end{aligned} \right. \quad (15)$$

The volumic energy of elastic strain damaged layer is:

$$2W_e = E_{11}^{0}(1 - \xi_{11})\varepsilon_{11}^{e} + C_{22}^{0}(\varepsilon_{22} + v_{12}^{0}\varepsilon_{11})_{-}^{2} + C_{22}^{0}(1 - d_{22}) \\ \times \langle \varepsilon_{22} + v_{12}^{0}\varepsilon_{11} \rangle_{+}^{2} + G_{12}^{0}(1 - d_{12})\gamma_{12}^{e} + G_{13}^{0}\gamma_{13}^{e} + G_{23}^{0}\gamma_{23}^{e}$$

Thermodynamic variables associated with the internal variables are thus defined by:

$$\left\{ \begin{aligned} Y_{12} &= -\frac{\partial W_e}{\partial d_{12}} = \frac{1}{2} G_{12}^0 \gamma_{12}^2 \\ Y_{22} &= -\frac{\partial W_e}{\partial d_{22}} = \frac{1}{2} C_{22}^0 \langle \varepsilon_{22} + v_{12}^0 \varepsilon_{11} \rangle_+^2 \end{aligned} \right. \quad (16)$$

$\langle X \rangle_+$ is the positive part of $X$:

$$\langle X \rangle_{+} = \begin{cases} X & \text{if } X > 0 \\ 0 & \text{if } X < 0 \end{cases}$$

and the threshold variables defining the fields of no damage are given by:

$$Y_m = \max_{\tau \le t} \sqrt{Y_{12} + bY_{22}}$$

→ Damage of the matrix

$$Y_i = \max_{\tau \le t} \sqrt{Y_{12}}$$

→ Damage of the interface in shearing

$$Y'_{i} = \max_{\tau \le t} \sqrt{Y_{22}}$$
where “max” defines the maximum value of the function considered.

→ Damage of the interface in transverse tensile

$$ \mu = \frac{\langle \sigma_{22} \rangle_{-}}{\sigma_{22}^{c}} - \left( \frac{\sigma_{22}}{2\sigma_{12}^{r}} \right)^{2} - 1) + \left( \frac{\langle \sigma_{22} \rangle_{-}}{2\sigma_{12}^{r}} \right)^{2} + \left( \frac{\sigma_{12}}{\sigma_{12}^{r}} \right)^{2} $$

→ Hashin criterion of matrix in compression [8].

The laws of evolution, have the form:

$$ \left\{ \begin{aligned} \dot{d}_{12} &= \dot{\lambda} \frac{\partial g}{\partial Y_{12}} \\ \dot{d}_{22} &= \dot{\lambda} \frac{\partial g}{\partial Y_{22}} \quad \text{avec } \dot{\lambda} > 0 \end{aligned} \right. \qquad (17) $$

### 3.2.1. The model

The model integrates the different rupture criteria. As the interface has a fragile behavior in traction transverse, it is natural to introduce an indicator of rupture such as:

$\chi_i = 0$ as long as there is no rupture, i.e. as long as $Y'_i < Y'_r$.

In the same way for the rupture of the matrix in transverse compression, if the Hashin criterion is verified the matrix breaks down, which results in defining $\chi_\mu$, indicator of rupture of the matrix, as follows:

$\chi_\mu = 0$ as long as there is no rupture, i.e. as long as $\mu < 1$.

If $\mu \ge 1$ then $\chi_\mu = 1$. Finally for the rupture of fibres, the selected criterion is a criterion of different maximum deformation in tensile and compression, that is to say: $\chi_f = 0$ as long as there is no rupture, i.e. as long as $\varepsilon_{11}^c < \varepsilon_{11} < \varepsilon_{11}^t$.

If $\varepsilon_{11} \ge \varepsilon_{11}^t$ or $\varepsilon_{11} \ge \varepsilon_{11}^c$ then $\chi_f = 1$.

Where $\varepsilon_{11}^t$ and $\varepsilon_{11}^c$, respectively, indicate the longitudinal deflections with rupture in tensile and compression, respectively.

### 3.3. Behavior of each component

The interface is elastic and prone to damage in shearing and elastic fragile in transverse tensile.

The matrix is prone to damage in tensile and fragile in compression (Hashin criterion).

The fibres are elastic linear fragile in tensile and fragile elastic nonlinear in compression.

The experiments showed that the variables of damage evolve in a linear way compared to the thresholds variables. These thresholds are introduced so as to ensure the irreversibility of the phenomena, through terms “max”. The laws of evolution of the damages and the fall of rigidity in longitudinal compression are selected as follows:

• In the longitudinal direction:

$$ \xi_{11} = \begin{cases} 0 & \text{if } \varepsilon_{11} > 0 \text{ and } \chi_f = 0 \\ -\gamma\varepsilon_{11} & \text{if } \varepsilon_{11} < 0 \text{ and } \chi_f = 0 \\ 1 & \text{if } \chi_f = 0 \end{cases} $$

• In shearing:

$$ d_{12} = \begin{cases} \frac{\langle Y_m - Y_0 \rangle_+}{bY_c} + \frac{\langle Y_i - Y'_0 \rangle_+}{Y'_c} & \text{if } \chi_i = 0 \\ 1 & \text{if not} \end{cases} $$

• In the transverse direction: If $\sigma_{22} \ge 0$

$$ d_{22} = \begin{cases} \frac{\langle Y_m - Y_0 \rangle_+}{Y_c} + \chi_i = 0 & \text{if } \chi_i = 0 \\ 1 & \text{if not} \end{cases} $$

If $\sigma_{22} \le 0$, $d_{22} = \chi_\mu$

with

$Y_0$: initial threshold of the matrix damage,

$Y'_0$: initial threshold of damage of the interface,

$Y_c$: resistance of the matrix,

$Y'_c$: resistance of the interface in tensile transverse,

$Y'_r$: resistance of the interface to shearing,

$\gamma$: loss of longitudinal rigidity in compression,

$b$: coupling coefficient which quantifies the relative influence of $\sigma_{12}$ and $\sigma_{22}$ in the degradation of material.

All these parameters were identified by Le Dantec [2].

The model of damage being elaborate, it is necessary to define the criterion of rupture.

### 3.4. Criterion of rupture and unstable condition

To have the criterion of rupture, we should calculate the instability condition given by:

$$ \det_{\tilde{\sigma} \to \infty} (K_{\text{tg}}) = 0 $$

After having developed the calculations, the conditions to satisfy could be written as follows.

$$ \left\| \begin{array}{l} P \frac{\langle \tilde{\sigma}_{22} \rangle_+^2}{C_{22}} + R \frac{\tilde{\sigma}_{12}^2}{G_{12}} + (Q - RP) \frac{\langle \tilde{\sigma}_{22} \rangle_+^2 + \tilde{\sigma}_{12}^2}{G_{12}C_{22}} = 1 \\ \text{or } E_{11} + \tilde{\sigma}_{11} \gamma H(-\varepsilon_{11}) H(\varepsilon_{11} - \varepsilon_{11}^c) = 0 \end{array} \right. $$

where the coefficients $P$, $Q$ and $R$ depend on the characteristics of the material and the expression of $\tilde{\sigma}$.

The introduction of the term $H(\varepsilon_{11})$ enables to avoid the singularities when the damage reaches “1” [9].

### 3.5. Analysis and results

The model of damage of the elementary layer worked out was established in the computer code by finite elements “Castem 2000”. A rectangular and isotropic transverse plate in T300/914 was composed of eight layers according to the sequence of stacking [+45°/−45°]$_2$s, of dimension (0.001 * 0.012 * 0.024) meters. One quarter of the plate (rectangular grid four quadratic elements with 16 nodes) was tested in uniform uniaxial tensile of density 12 N/m$^2$ (Fig. 1).
<figure><img src="image_1.png" /><figcaption>Fig. 1. Charge and discharges experiment on [+45°/−45°]₂s in tensile.</figcaption></figure>

According to this curve the structure presents a no linear behavior, moreover, during the discharges, some residual stress as well as variations of slopes translating the rigidity losses of the structures appear. A calculation based on the model of damage and simulating this experiment was carried out. The results obtained are presented in Fig. 2.

The analysis and the comparison of the numerical and experimental results confirm the presence of damage within

the structure laminated composites. Indeed, the introduction of the model of damage accounts well for the loss of rigidity of materials appearing in the experiment through variations of slope and the nonlinear character of the behavior of the structure.

However, the fact that in simulations carried out, the behavior is more rigid than at the time during the tests, and that the residual stress are nonexistent, lead us to believe that the damage developed within this material cannot be studied within the framework of elasticity. It would be thus interesting in order to better account for the real behavior of the laminated composites, and in particular of the existence of the residual deformations to model plasticity in the composites.

### 4. Modeling of plasticity with isotropic hardening and nonlinear kinematic assumptions

The plasticity is due to the matrix and the interface matrix fibre, the flow is blocked in the direction of fibres $\dot{e}_{11}^p = 0$ [10]. The model uses only the quantities of shearing and transverse, the damage is taken into account by the introduction of $\frac{\sigma}{1-d}$ in the criterion of plasticity. Moreover, in this model nonlinear kinematic hardening is introduced. The whole is placed within the framework of small transformations, plane strain and nonassociated plasticity [11,12].

The criterion of plasticity is written as:

$$f(\sigma, X, R) = \sqrt{H(1-d_{12}) \left[ \frac{\sigma_{12} - X_{12}}{1-d_{12}} \right]^2 + A^2 H(1-d_{22}) \left[ \frac{\sigma_{22} - X_{22}}{1-d_{22}} \right]^2} - R - R_0$$

<figure><img src="image_2.png" /><figcaption>Fig. 2. Charge and discharges with damage model on [+45°/−45°]₂s in tensile.</figcaption></figure>

The introduction of the terms $H(1-d)$ allows to avoid the singularities when the damage reaches value “1”. Because in this case there is rupture.

$$H(1 - d_{12}) = 0 \text{ if } d_{12} = 1 \text{ else } H(1 - d_{12}) = 1.$$

$$H(1 - d_{22}) = 0 \text{ if } d_{22} = 1 \text{ else } H(1 - d_{22}) = 1.$$

The laws of normality are given by:

$$\dot{e}_{22}^p = \dot{\lambda} \frac{\partial f}{\partial \sigma_{22}} = \frac{\dot{\lambda}}{f + R + R_0} A^2 \frac{\sigma_{22} - X_{22}}{(1 - d_{22})^2}$$

$$\dot{e}_{12}^p = \dot{\lambda} \frac{\partial f}{\partial \sigma_{12}} = \frac{\dot{\lambda}}{f + R + R_0} \frac{\sigma_{12} - X_{12}}{(1 - d_{12})^2}$$

$$\dot{p} = -\dot{\lambda} \frac{\partial f}{\partial R} = \dot{\lambda}$$

$$\dot{\alpha}_{22} = -\dot{\lambda} \frac{\partial f}{\partial X_{22}} = \dot{\lambda} \frac{\partial f}{\partial \sigma_{22}} = \dot{e}_{22}^p$$

$$\dot{\alpha}_{12} = -\dot{\lambda} \frac{\partial f}{\partial X_{12}} = \dot{\lambda} \frac{\partial f}{\partial \sigma_{22}} = \dot{e}_{12}^p$$

and the complementary laws by:

$$\dot{X}_{22} = c \dot{\alpha}_{22} - a X_{22} \dot{p} = c \dot{e}_{22}^p - a X_{22} \dot{p}$$

$$\dot{X}_{12} = c \dot{\alpha}_{12} - a X_{12} \dot{p} = c \dot{e}_{12}^p - a X_{12} \dot{p}$$
To obtain “X” terms of nonlinear kinematic hardening, we have to apply a method of integration, “Euler’s explicit method” and R = βp<sup>α</sup>.

R₀: initial threshold of plasticity.

“c”, “a” (parameters of nonlinear kinematic hardening) and “α”, “β” (parameters of isotropic hardening), characteristics identified on a laminate T300/914 with [+45°, -45°]<sub>2s</sub>.

Finally to determine “λ” the plasticity multiplier, it is necessary to solve f = ŷ = 0.

$$
\begin{align*}
f + \frac{\partial f}{\partial \sigma} d\sigma + \frac{\partial f}{\partial R} dR + \frac{\partial f}{\partial X} dX &= 0 \\
\Rightarrow f + \frac{\partial f}{\partial \sigma} \Lambda d\epsilon^c + \frac{\partial f}{\partial R} \frac{\partial R}{\partial p} dp + \frac{\partial f}{\partial X} \frac{\partial X}{\partial \alpha} d\alpha &= 0 \\
\Rightarrow f - \frac{\partial f}{\partial \sigma} \Lambda \frac{\partial f}{\partial \sigma} d\lambda - \frac{\partial R}{\partial p} d\lambda + \frac{\partial f}{\partial X} \left(c \frac{\partial f}{\partial \sigma} - aX\right) d\lambda &= 0
\end{align*}
$$

Finally

$$
d\lambda = \frac{f}{\frac{\partial f}{\partial \sigma} \Lambda \frac{\partial f}{\partial \sigma} - \frac{\partial R}{\partial p} + \frac{\partial f}{\partial X} \left( c \frac{\partial f}{\partial \sigma} - aX \right)}
$$

4.1. Method of identification of the parameters of the model of plasticity

4.1.1. Model with nonlinear kinematic hardening

To identify the constants “c” and “a” of the model with nonlinear kinematic hardening, a tensile test on the laminate 300/914 with [+45°/-45°]<sub>2s</sub> was studied. This test is interesting, because it was noticed that it was equivalent to a case of quasi-pure shearing, moreover it is rather easy to obtain an analytical solution of this problem.

Indeed, in the case of a quasi-pure shearing σ₁₁ = σ₂₂ = 0, the total strain and stress are related according to the formulas of rotations on the strain and the local stress by:

$$
\left\{
\begin{array}{l}
\sigma_{12} = \frac{\sigma_{xx}}{2} \\
\varepsilon_{12} = \frac{\varepsilon_{xx} - \varepsilon_{yy}}{2}
\end{array}
\right.
\qquad (18)
$$

The elastoplastic problem in the case of a quasi-pure
shearing is written as:

$$
\gamma = \gamma_e + \gamma_p
$$
$$
\sigma_{12} = G(1 - d_{12})\gamma_e
$$

with

$d_{12} = g\gamma_e$

and

$$
g = \sqrt{\frac{G}{2} \left( \frac{1}{bY_c} + \frac{1}{Y'_c} \right)}
$$

$$
\left.
\begin{aligned}
f &= \frac{\sigma_{12} - X_{12}}{1 - d_{12}} - R_0 \\
\dot{\gamma}_p &= \frac{\dot{\lambda}}{1 - d_{12}} \\
\dot{X}_{12} &= c \dot{\gamma}_p - a X_{12} \dot{\lambda} \\
\Rightarrow \dot{X}_{12} &= [c - a(1 - d_{12})X_{12}] \dot{\gamma}_p
\end{aligned}
\right\}
$$

It remains to obtain a curve making it possible to deduce
the constants “a” and “c”.

Carrying out the variable change $u = (1 - d_{12})X_{12}$, $\frac{\dot{X}_{12}}{\dot{\gamma}_p} = c - au$ is a linear function, where “c” represents the ordinate at the origin and “a” the slope of the curve representative of this linear function.

It now remains to apply this change of variable to all the
equations to calculate:

$$
\frac{\dot{X}_{12}}{\dot{\gamma}_p}
$$

where

$$
\begin{align}
t &\rightarrow \sigma \rightarrow \gamma_e^{(19)} \rightarrow u^{(20)} \notag \\
\sigma &= G(1 - g\gamma_e)\gamma_e \tag{19}
\end{align}
$$

$$
u = (1 - d_{12})X_{12} = (1 - g\gamma_e)^2(G\gamma_e - R_0) \quad (20)
$$

However: X₁₂ = σ - R₀(1 - d₁₂) in the elastic yield and:
d₁₂ = gγₑ, where

$$
\frac{d}{dt} = \frac{d}{du} \frac{du}{d\gamma_e} \frac{d\gamma_e}{d\sigma} \frac{d\sigma}{dt} \quad (21)
$$

Thus

$$
\frac{d\gamma_p}{du} = \frac{dX_{12}}{du} \frac{1}{c - au}
$$

$$
\Rightarrow c - au = \frac{du}{\frac{dX_{12}}{du}} \qquad (23)
$$

However

$$
\frac{dX_{12}}{du} = \frac{dX_{12}}{d\gamma_e} \frac{d\gamma_e}{du} \quad \text{and} \quad \frac{d\gamma_p}{du} = \frac{d\gamma_p}{d\gamma_e} \frac{d\gamma_e}{du} \qquad (24)
$$

Thus

$$
\frac{dX_{12}}{du} \frac{dX_{12}}{d\gamma_e} \frac{d\gamma_e}{du} = \frac{dX_{12}}{d\gamma_p} \frac{d\gamma_p}{du} \quad (25)
$$

$$
\Rightarrow \frac{dX_{12}}{du} = \frac{dX_{12}}{d\gamma_e} \frac{d\gamma_e}{du} \qquad (26)
$$

It thus remains to deduce $\frac{d\gamma_p}{d\gamma_e}$ from experimental results so we use an approximation by cubic splines of the curve $\gamma_p = fc(\gamma_e)$, deduced from the experimental results as explained below.
### 4.1.2. Methodology of identification

From the test (Fig. 1), a series of points ($\sigma_{xx}$, $\epsilon_{xx}$, $\epsilon_{yy}$) are recorded. It is possible to calculate $\sigma_{12}$ and $\gamma_{exp}$ according to (18) (Fig. 3).

Knowing $\sigma$, $\gamma_e$ is calculated by (19).

Once $\gamma_e$ obtained, for each one of these deformations it is necessary to calculate:

* $\gamma_p$ by: $\gamma_p = \gamma_{exp} - \gamma_e$.

* $u$ by: $ru = (1 - d_{12})X_{12} + (1 - g\gamma_e)^2(G\gamma_e - R_0)$.

* $\frac{dX_{12}}{d\gamma_e}$ by: $\frac{dX_{12}}{d\gamma_e} = (1 - g\gamma_e)G - g(G - R_0)$.

* $\frac{d\gamma_p}{d\gamma_e}$ by the splines method according to the experimental curve $\gamma_p = fc(\gamma_e)$ (Fig. 4).

<figure><img src="image_3.png" /><figcaption>Fig. 3. Experimental curve of one T300 [+45°/-45°]$_2$s in tensile.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4. Plastic strain function of the elastic strain.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. Identification of the parameters of the model of plasticity to nonlinear kinematic hardening.</figcaption></figure>

Finally, a curve ($\frac{\dot{X}_{12}}{\dot{\gamma}_p}$, $u$) is obtained (Fig. 5). Once this curve is plotted, we take the line passing through these points and deduce “c” as the ordinate at the origin and “a” the slope of this curve. We consider only the linear part of the curve where the stress is sufficiently significant so that hardening is important. Thus we take the slope line presented in Fig. 5.

The results obtained are: “c = 1.054 × 10¹⁰ Pa” and “a = 275.834”.

### 4.1.3. Model with isotropic hardening

In this case, the isotropic hardening is given by the law $R = \beta p^\alpha$ with “$\beta$” and “$\alpha$” parameters to be identified. So, the tensile test on the laminate [+45°/-45°]$_2$s in T300, studied before is used. Eqs. (18) remain always valid. However, it is now assumed that only isotropic hardening occurs.

The equations of the problems are thus:

$$ 
\begin{aligned}
\gamma &= \gamma_e + \gamma_p \\
\sigma_{12} &= G(1 - d_{12})\gamma_e \\
d_{12} &= g\gamma_e \quad \text{and} \quad g = \sqrt{\frac{G}{2}\left(\frac{1}{bY_c} + \frac{1}{Y'_c}\right)} \\
f &= \frac{\sigma_{12}}{1 - d_{12}} - R - R_0 \quad \text{with } R = \beta p^\alpha \\
\dot{\gamma}_p &= \frac{\dot{\lambda}}{1 - d_{12}} \\
\dot{p} &= \dot{\lambda}
\end{aligned} 
\quad (27) 
$$

It thus remains to find a simple curve between $R$ and $p$ to identify “$\beta$” and “$\alpha$”, however $\log R = \log \beta - \alpha \log p$. Therefore if it is possible to show a curve $\log R - \log p$, “$\alpha$” and “$\beta$” can be identified. Indeed $\log \beta$ represents the ordinate at the origin and $\alpha$ the slope of the curve. We have to obtain the curve $\log R - \log p$.
<figure><img src="image_6.png" /><figcaption>Fig. 6. Identification of the parameters of the model of plasticity to isotropic hardening.</figcaption></figure>

### 4.1.4. Methodology of identification

It is possible to easily obtain $\log R = fc(\gamma_e)$ indeed, by (19) it comes:

$$ 
\begin{aligned}
R &= G\gamma_e - R_0 \\
\Rightarrow \log R &= \log(G\gamma_e - R_0)
\end{aligned}
$$

It is thus necessary to try to have

$\log p = fc(\gamma_e)$

$$ \text{However } \frac{dp}{d\gamma_e} = \frac{d\gamma_p}{d\gamma_e} (1 - g\gamma_e). $$

<figure><img src="image_7.png" /><figcaption>Fig. 7. Charge and discharge on T300 [+45°/-45°]₂s in tensile. Models damage crossed to plasticity with nonlinear kinematic hardening.</figcaption></figure>

And by integrating by parts:

$$ p = [(1 - g\gamma_e)\gamma_p] + g \int \gamma_p d\gamma_e $$

Having the experimental curve $\sigma_{12} = fc(\gamma_{exp})$ (Fig. 3), it is easy to deduce the curve from it $\gamma_p = fc(\gamma_e)$ (Fig. 4) and to know $\int \gamma_p d\gamma_e$ by using trapezoidal integration method. Thus it is easy to calculate $p = fc(\gamma_e)$ and $\log p = fc(\gamma_e)$.

Having $\log R$ and $\log p$ function of “$\gamma_e$”, the curve log $R - \log p$ (Fig. 6) is obtained, and it only remains to identify “$\log \beta$” as the ordinate at the origin and “$\alpha$” as the slope passing through these points. As for nonlinear kinematic hardening, we consider the field where the stress is significant, so that hardening is important, and corresponds to the linear part of the curve (Fig. 6). The values obtained are: “$\alpha$ = 0.481” and “$\beta$ = 636 × 10⁶ Pa”.

## 5. Conclusion

The analysis and the comparison of the numerical and experimental results confirm the presence of damage within the laminated composites structures. Indeed the introduction of the damage model accounts well for the loss of rigidity of materials appearing in the experiment through variations of slope and nonlinear character of the behavior of the structure. Then a calculation of charge and discharge with the nonlinear kinematic hardening model was carried out. The results obtained are presented in Fig. 7. By comparing these results with the experiment, it is interesting to note that the residual stress is in good agreement with the experimental stress. Moreover, hysteresis loops which were inexistent before start to appear. However these loops do not have the same amplitude as those of the experiment.

With the sight of these results, it is possible to believe that the identified parameters account relatively well for the real behavior of the laminated composites structures studied and that other phenomena are to be taken into account in the modeling of the composites structures, such as: aging by updating the coefficients characteristic and the viscoplasticity.

Parameters of T300/914:

<table><tr><td>ν</td><td>0.33</td></tr><tr><td>E<sub>1</sub></td><td>0.148 × 10<sup>12</sup> (Pa)</td></tr><tr><td>E<sub>2</sub></td><td>0.108 × 10<sup>11</sup> (Pa)</td></tr><tr><td>G<sub>12</sub></td><td>0.58 × 10<sup>10</sup> (Pa)</td></tr><tr><td>G<sub>23</sub></td><td>0.36 × 10<sup>10</sup> (Pa)</td></tr><tr><td>Y<sub>0</sub></td><td>300 (Pa<sup>1/2</sup>)</td></tr><tr><td>Y<sub>c</sub></td><td>3600 (Pa<sup>1/2</sup>)</td></tr><tr><td>Y'<sub>0</sub></td><td>0 (Pa<sup>1/2</sup>)</td></tr><tr><td>Y'<sub>c</sub></td><td>7500 (Pa<sup>1/2</sup>)</td></tr><tr><td>Y'<sub>r</sub></td><td>450 (Pa<sup>1/2</sup>)</td></tr><tr><td>σ<sup>t</sup><sub>11</sub></td><td>0.15 × 10<sup>10</sup> (Pa)</td></tr><tr><td>σ<sup>t</sup><sub>11</sub></td><td>0.15 × 10<sup>10</sup> (Pa)</td></tr><tr><td>σ<sup>t</sup><sub>22</sub></td><td>0.55 × 10<sup>8</sup> (Pa)</td></tr><tr><td>σ<sup>c</sup><sub>22</sub></td><td>0.205 × 10<sup>9</sup> (Pa)</td></tr></table>
<table><tr><td>σ<sub>12</sub><sup>r</sup></td><td>0.105 × 10<sup>9</sup> (Pa)</td></tr><tr><td>α</td><td>1.0</td></tr><tr><td>γ</td><td>22.4</td></tr><tr><td>ε<sub>11</sub><sup>t</sup></td><td>0.011</td></tr><tr><td>ε<sub>11</sub><sup>c</sup></td><td>0.011</td></tr><tr><td>A<sup>2</sup></td><td>0.38</td></tr><tr><td>R<sub>0</sub></td><td>40 × 10<sup>6</sup> (Pa)</td></tr><tr><td>α<sub>plast</sub></td><td>0.481</td></tr><tr><td>β</td><td>636 × 10<sup>6</sup> (Pa)</td></tr><tr><td>c</td><td>1.054 × 10<sup>10</sup> (Pa)</td></tr><tr><td>a</td><td>275.834</td></tr></table>

References

[1] Giletta De St. Joseph D. Modélisation mécanique et identification de la couche élémentaire. Thèse du L.M.T, 1985.

[2] Le Kantec E. Contribution à la modélisation du comportement des composites stratifiés. Thèse Univ. P.M. Curie, 1989.

[3] Ladeveze P, Allix O, Deü J-F, Lévèque D. A mesomodel for localisation an damage computation in laminates. Comput Methods Appl Mech Eng 2000;183:105–22.

[4] Reese S. Meso-macro modeling of fibre-reinforced rubber-like composites exhibiting large elastoplastic deformation. Compos Sci Technol 2001;61:2223–30.

[5] Herakovitch Carl T, Crobert D, Gasser A, Guitard L. Damage evolution in laminates with fiber rotation. Compos Sci Technol 2000;60:2781–9.

[6] Allen David H. Homogenization principles and their application to continuum damage mechanics. Compos Sci Technol 2001;61: 2223–30.

[7] Devries F, Dumontet H, Duvaut G, Léné F. Homogenization and damage for composite structures. Int J Numer Methods Eng 1989;27:285–98.

[8] Barbero EJ, Reddy JN. Nonlinear analysis of composite laminates using a generalized laminated plate theory. AIAAJ 1990;28(11):1987–94.

[9] Vechery G, Hamelin P, Massard T, Aivazzadeh S. Calcul des structures en matériaux composites, Paris, 1987.

[10] Murray Y, Schwer L. Implementation and verification of fiber composite damage models. AMD ASME 1990;167:21–30.

[11] Lemaitre J. Coupled elastoplasticity and damage constitutive equa- tions. Comput Methods Appl Mech Eng 1985;51:31–49.

[12] Lemaitre J, Chaboche JL. Mécanique des matériaux solides, Dunod, 1988.