# Modelling the hysteresis composite behavior using an elasto-plasto-damage model with fractional derivatives

Alina Krasnobrizha, Patrick Rozycki, Pascal Cosson, Laurent Gornet

► To cite this version:

Alina Krasnobrizha, Patrick Rozycki, Pascal Cosson, Laurent Gornet. Modelling the hysteresis composite behavior using an elasto-plasto-damage model with fractional derivatives. 20th International Conference on Composite Materials - ICCM20, Jul 2015, Copenhagen, Denmark. hal-01240728

HAL Id: hal-01240728

https://hal.science/hal-01240728v1

Submitted on 9 Dec 2015

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 CC BY 4.0 - Attribution - International License
MODELLING THE HYSTERESIS COMPOSITE BEHAVIOR USING AN
ELASTO-PLASTO-DAMAGE MODEL WITH FRACTIONAL
DERIVATIVES

Alina Krasnobrizha¹, Patrick Rozycki¹, Pascal Cosson¹, Laurent Gornet¹

¹Institut de recherche en Génie Civil et Mécanique (GeM) - UMR CNRS 6183, Ecole Centrale de
Nantes, 1 rue de la Noë, BP 92101, 44321 Nantes Cedex 3, France

Email: alina.krasnobrizha@ec-nantes.fr,

patrick.rozycki@ec-nantes.fr, pascal.cosson@ec-nantes.fr, laurent.gornet@ec-nantes.fr

**Keywords:** Composite material, Behavior law, Fractional derivative, Hysteresis, Characterisation

ABSTRACT

This paper deals with a new model of an orthotropic elementary ply, which consists of two sub-
models. The first one is dedicated to treat the behavior during loading. The elastic and inelastic strains
are computed as well as the in-ply damages. The strain-rate sensitivity can be also taken into account.
The second sub-model involves a fractional derivate approach to describe the hysteretic behavior
observed during unloading/reloading. Instead of using the viscoelastic strain, we express it by a
fractional derivative. The behavior law requires few parameters. They can be easily identified through
an optimization algorithm from the experimental data. The second sub-model calculates different
quantities of dissipated energy during the unloading. The model is validated for thermoset and
thermoplastic composite materials.

# 1 INTRODUCTION

The extensive use of composite materials in industrial applications requires a better understanding
of their mechanical behavior. Composites materials are anisotropic and heterogeneous materials. A
complex model is required to adequately describe their behaviors. There are a lot of works concerning
the behavior of unidirectional or woven composites with thermosetting matrix under quasi-static
loadings [1-12]. These models are able to describe damages and inelastic strains that are appearing in
shear and transverse directions. The strain rate dependent models were developed for unidirectional
composites [13-15], and for woven composites [16, 17].

The previous developed models are well describing different physical phenomena such as the
elastic and inelastic strains, the damages whatever the orthotropic directions of composite material.
However, they cannot represent the hysteretic behavior during unloading/reloading. The work is based
on fractional derivative approach which is a good issue for some natural structures and modern
heterogenic materials such as elastomers and polymers. There is a significant number of works about
the fractional viscoelastic constitutive equations for different materials under various types of loading.
The fractional Zeiner model was applied by Caputo to represent behavior of glass and some metals
[18]. Accordance between the results of an empirical model based on fractional derivative and the
experimental data for some polymers and elastomers was presented by Bagley and Torvik in [19].
Rouse found a relation between molecular theory of polymers and viscoelastic law based on fractional
derivatives [20]. Rabotanov presented a generalized rheological model to describe behavior of
hereditary medium [21]. The fractional derivative approach was successfully applied to represent the
rheological behavior of organic glass, elastomers, polyurethane, polyisobutylene, amorphous solid
polymers in a wide temperature range [22-27]. The hysteresis cycles were represented by Caputo for
some metals under fatigue loading [28]. The constitutive model based on the fractional derivative
approach for composite materials under cycling quasi-static shear loading was proposed by Mateos
[29, 30]. This model is also able to represent hysteresis composite behavior using a few material
parameters. The elastic and irreversible strains, damage and strain rate effects are taken into account
but not separately. The purpose of our work is to consider a similar approach which allows us to compute the amount of dissipated energy caused by each of these phenomena.

In this paper we combine an elastic-plastic damage behavior law [6] with strain-rate sensitivity [16] and a fractional derivative approach. A complete composite behavior model with accurate dissipations is needed to compute self-heating during fatigue test. A method based on self-heating tests has been developed [31] and provides a fast identification of the endurance limit of the composite material in comparison with conventional methods (S-N curves). It requires an analysis and an understanding of all the thermo-mechanical phenomena appearing in the material during the test, especially in the hysteresis loops. To simplify the analysis and to increase the potential of this method, numerical simulations of the fatigue tests could lead to a better estimation of the different amounts of dissipated energy. Therefore, the thermo-mechanical finite element simulations using the proposed constitutive model with accurate dissipations can be applied to fatigue problems.

## 2 THEORETICAL MODEL FOR COMPOSITE PLY

The constitutive model is based within the framework of thermodynamics isothermal irreversible processes. The orthotropic mechanical behavior of a ply of a composite laminate is described at the mesoscopic scale. We consider a plane-stress state. Subscripts 1 and 2 represent the fiber and the transverse directions, respectively. The continuum damage mechanics theory describes damages such as matrix micro-cracking and fiber/matrix debonding. An isotropic hardening model is used for plasticity. The viscoelastic effects are expressed by fractional derivatives.

The material parameters of the law are obtained from an experimental campaign based on cyclic tensile tests. The fiber orientation at 0° allows characterizing the longitudinal behavior of elementary ply and the ± 45°, the plane shear behavior. To characterize the transverse behavior, we can choose 90° orientation. However, the material in this direction has a brittle behavior and parameters cannot be determined precisely. Therefore, it is suitable to use a test with a ±67.5° laminate to provide a coupling between the transverse direction and shear.

Analyzing the experimental campaign results, the biggest non-linear effects appear during quasi-static shear test. So our work is focused on the simulation of [±45]s tensile test. The experimental data of a unidirectional carbon/epoxy composite (Fig. 1) and woven carbon/PA66 material were used to study thermoset and thermoplastic material behavior, respectively.

<figure><img src="image_2.png" /><figcaption>Figure 1: Shear stress-strain curve for carbon/epoxy composite.</figcaption></figure>

The developed constitutive model is composed of two sub-models. The first one deals with damage of the elementary ply with isotropic hardening. And the second one describes the hysteresis behavior based on a fractional derivative approach.

Within the framework of thermodynamic irreversible theory, we choose Helmholtz potential depending on the internal variables:

$$ \rho\psi = \rho\psi(\varepsilon_{12}^e, d_{12}, p) \quad (1) $$
where $\rho\psi$ is a volumetric density of free-energy, $\epsilon_{12}^e$ is the elastic shear strain, $d_{12}$ is a damage internal variable and $p$ is a cumulative plasticity.

To formulate constitutive equations in terms of strain, the Helmholtz potential has been chosen [8, 32]. On the contrary, the potential of Gibbs provides a problem formulation in stress.

## 2.1 Damage model

Under the quasi-static shear loading, the damage mainly appears in the form of fiber-matrix debonding. Following the second law of thermodynamic, the expression for elastic strain energy of the damaged material in shear case can be written as:

$$W_D = \frac{1}{2} \sigma : \varepsilon = \frac{1}{2} G_{12}^0 (1 - d_{12}) (2 \varepsilon_{12}^e)^2 \quad (2)$$

where $G_{12}^0$ is an initial shear modulus.

The damage variable $d_{12}$ represents a loss of material stiffness which can be determined by the diminution of the shear modulus during experiment:

$$d_{12} = 1 - \frac{G_{12}^i}{G_{12}^0} \quad (3)$$

where $G_{12}^i$ is the current shear modulus of the loop.

We introduce the following notation for effective stress:

$$\tilde{\sigma}_{12} = \frac{\sigma_{12}}{1 - d_{12}} \quad (4)$$

Therefore, the constitutive law can be deduced:

$$\epsilon_{12}^e = \frac{\sigma_{12}}{2G_{12}^0(1 - d_{12})} = \frac{\tilde{\sigma}_{12}}{2G_{12}^0} \quad (5)$$

The conjugate quantity of the thermodynamical force $Y_{12}$ is given by (6). As Helmholtz potential is used, the associated thermodynamical force depends on strain $\gamma_{12}^e = 2\epsilon_{12}^e$ and is independent of the damage variable $d_{12}$.

$$Y_{12} = \frac{\partial \rho \psi}{\partial d_{12}} = - \frac{\partial W_D}{\partial d_{12}} = \frac{1}{2} G_{12}^0 (2 \epsilon_{12}^e)^2 = \frac{1}{2} G_{12}^0 (\gamma_{12}^e)^2 \quad (6)$$

The threshold of undamaged zone is defined like the maximal thermodynamical force in shear for all previous time $\tau$ up to the current time $t$:

$$\bar{Y}_{12} = \max_{\tau \le t} \sqrt{Y_{12}} \quad (7)$$

Generally, the damage evolution is approximated by a linear function of the thermodynamical force. In order to describe more composite materials behaviors, the evolution damage law can be fitted by a polynomial approximation as (8) [33]:

$$d_{12} = \sum_{n=1}^{N} a_n (\bar{Y}_{12} - \bar{Y}_0)^n \quad (8)$$

if $d_{12} < 1$ and $\bar{Y}_{12} < \bar{Y}_{12}^{\text{failure}}$; otherwise $d_{12} = 1$

where $\bar{Y}_0$ is the initial threshold of damage and $\bar{Y}_{12}^{\text{failure}}$ is the failure-damage threshold. The constants in this law are material parameters.
## 2.2 Plasticity modelling and damage-plasticity coupling

In order to take into account inelastic strain caused by material damage, the effective stress (4) is introduced in a plasticity model. Assuming an isotropic strain hardening, the elastic domain is defined by the yield function:

$$f = \sqrt{\bar{\sigma}_{12}^2 - R(p) - R_0} \qquad (9)$$

where $R_0$ is the yield stress and the function $R(p)$ is a material characteristic function, which depends on the cumulative plastic strain defined by $p$. Generally, the hardening function $R(p)$ is approximated by a power law:

$$R = \beta p^k \quad \text{with} \quad p = \int_{0}^{\varepsilon_{12}^p} (1 - d_{12}) d \varepsilon_{12}^p \qquad (10)$$

where $\beta$ and $k$ are material parameters, identified from the experimental data. In order to suit to more composite material behaviors, a polynomial approximation of the hardening function $R(p)$ is used in following form [33]:

$$R = \sum_{m=1}^{M} b_m p^m \qquad (11)$$

The material parameter identification procedure is described in the next sections. The resulting stress-strain curve is obtained for unidirectional carbon/epoxy composite (Fig. 2). The model represents elasto-plastic damage behavior of composite. The simulation results are in good agreement with experimental data. However, hysteresis loops cannot be represented by this model.

<figure><img src="image_3.png" /><figcaption>Figure 2: Comparison of experimental and simulated behavior.</figcaption></figure>

## 2.3 Fractional derivative model

The hysteresis loops appearance is associated with energy dissipation in composite under cyclic quasi-static shear loading. Hysteresis is a hereditary phenomenon, i.e. the past loading history has to be taken into account. The idea is to introduce fractional derivatives in the constitutive equation to provide a material history dependence. The Riemann-Liouville definition of the fractional derivative is [34, 35]:

$$D^{\alpha} f(x) = \frac{1}{\Gamma(1-\alpha)} \frac{d}{dx} \int_{0}^{x} (x-t)^{-\alpha} f(t) dt \qquad (12)$$

where $D^\alpha$ is a fractional derivative of order $\alpha$, and $\Gamma$ is the Gamma-function defined by:
$$ \Gamma(z) = \int_{0}^{\infty} e^{-x} x^{z-1} dx, \quad z \in \mathbb{R}_{+} \qquad (13) $$

Assuming $\epsilon_{(t=0)} = 0$, we obtain the stress/strain relation in terms of fractional derivative in the context of linear viscoelasticity:

$$ \sigma(t) = K_0 D^{\alpha} \epsilon(t), \quad \text{for } 0 < \alpha < 1 \qquad (14) $$

where $K_0$ is a material parameter.

This expression represents a constitutive stress/strain relation for materials with memories effects within the framework of small strains, and complies with Volterras' principles of hereditary mechanics [36-38]. The viscoelastic state takes an intermediate form between Hook's law for $\alpha = 0$ and Newton's rheological law with $\alpha = 1$.

In order to represent hysteretic behavior observed in experimental test (Fig. 1), we propose to express the strain by a fractional derivative. We assume that both damage and inelastic strain stay constant during the unloading. Therefore, the strain corresponding to the hysteresis loop can be expressed such as:

$$ \gamma_{12}^{\text{test}} = \frac{\sigma_{12}^{\text{test}}}{G_{12}^{0}(1 - d_{12})} \qquad (15) $$

where $\sigma_{12}^{\text{test}}$ is the experimental value of stress.

The figure 3 represents the experimental elastic strain $\gamma_{12}^{\text{test}}$ and the numerical elastic strain $\gamma_{12}^{\text{elast}}$ calculated by the elasto-plastic damage model (Fig. 2). The difference between the curves corresponds to the magnitude of viscoelastic strain $\gamma_{12}^{\text{ve}}$ which provides the appearance of hysteresis loops. To represent this non-linearity, we assume that total strain within hysteresis loop is defined by:

$$ \gamma_{12}^{\text{test}}(t) = \gamma_{12}^{\text{elast}}(t) + \gamma_{12}^{\text{ve}}(t) = A + BD^{\alpha}\gamma_{12}^{\text{elast}}(t) \qquad (16) $$

where $A, G_{12}^1, \alpha$ are parameters of fractional derivative model, $D^\alpha$ is a fractional derivative in the Riemann-Liouville sense (12). The constitutive relation in terms of fractional derivatives defined by:

$$ \sigma_{12}(t) = G_{12}^{0}(1 - d_{12})(\gamma_{12}^{\text{elast}}(t) + \gamma_{12}^{\text{ve}}(t)) = G_{12}^{0}(1 - d_{12})A + G_{12}^{1}D^{\alpha}\gamma_{12}^{\text{elast}}(t) \qquad (17) $$

<figure><img src="image_4.png" /><figcaption>Figure 3: Elastic strain.</figcaption></figure>

## 2.3 Coupling of two sub-models

We combine the two sub-models. The switching between the models is performed automatically, depending on the sign of the yield function (9) and its derivative. If $f = 0$ and $\dot{f} = 0$, the elasto-plastic damage model is used. During unloading if $f < 0$ or $f = 0$ and $\dot{f} < 0$ and during reloading if $f < 0$
or $f = 0$ and $\dot{f} > 0$, we consider that damage and plastic strain stay constant, thus we switch to the fractional derivative model.

### 3 IDENTIFICATION OF PARAMETERS

In the following section, the identification procedure concerning in-ply damage, isotropic strain hardening and fractional derivative model is described. The suggested methodology of experimental identification is applied to unidirectional carbon/epoxy composite under shear quasi-static loading (Fig. 1).

#### 3.1 Damage and plastic strain evaluation

The Fig. 4 corresponds to the evolution of the shear damage evaluated by (3) versus the thermodynamical associated force issued from (7). We choose to interpolate this experimental data by a 4th order polynom. The values of material parameters are presented in table 1.

<figure><img src="image_5.png" /><figcaption>Figure 4: Damage master curve of elementary ply.</figcaption></figure>

<table><caption>Table 1: Damage model parameters.</caption><thead><tr><th>Engineering constant</th><th>Unit</th><th>Identification</th></tr></thead><tbody><tr><td>Y&#773;<sub>0</sub></td><td>[√MPa]</td><td>0.188</td></tr><tr><td>a<sub>1</sub></td><td>[(1/√MPa)<sup>4</sup>]</td><td>-0.0921</td></tr><tr><td>a<sub>2</sub></td><td>[(1/√MPa)<sup>3</sup>]</td><td>-0.0463</td></tr><tr><td>a<sub>3</sub></td><td>[(1/√MPa)<sup>2</sup>]</td><td>0.1096</td></tr><tr><td>a<sub>4</sub></td><td>[1/√MPa]</td><td>0.4005</td></tr></tbody></table>

The cumulative plastic strain $p$ can be obtained using identified damage parameters. The strain hardening function $R(p)$ is determined by fitting the measured values of the effective stress and the cumulative plastic strain $p$ in accordance with (18). The fitting of experimental points (Fig. 5) is performed by a 5th degree polynom. The corresponding coefficients are presented in table 2.

$$R(p) = \frac{\sigma_{12}}{(1 - d_{12})} - R_0 \approx \sum_{m=1}^{5} b_m p^m \quad (18)$$
<figure><img src="image_6.png" /><figcaption>Figure 5: Plasticity mater curve of elementary ply.</figcaption></figure>

<table><caption>Table 2: Model parameters associated with plastic strain.</caption><thead><tr><th>Engineering constant</th><th>Unit</th><th>Identification</th></tr></thead><tbody><tr><td>b<sub>1</sub></td><td>[MPa]</td><td>292.66</td></tr><tr><td>b<sub>2</sub></td><td>[MPa]</td><td>-52.777 10<sup>6</sup></td></tr><tr><td>b<sub>3</sub></td><td>[MPa]</td><td>4.65 10<sup>9</sup></td></tr><tr><td>b<sub>4</sub></td><td>[MPa]</td><td>-192.8 10<sup>9</sup></td></tr><tr><td>b<sub>5</sub></td><td>[MPa]</td><td>3029.45 10<sup>9</sup></td></tr></tbody></table>

## 3.2 Identification of fractional derivative model parameters

In order to determine fractional derivative model parameters, an optimization problem has been resolved. The fractional derivative is a non-local space/time operator and depends on the material loading history. The elastic strain fractional derivatives are calculated for each loop separately. The fractional model parameters are determined from the last point of each loop. The zero initial conditions are required to avoid a singularity in fractional calculus. To insure that, the fractional derivatives are calculated on a time-interval corresponded to “previous loop loading-unloading-reloading” phase. In our case the elastic strain $\gamma_{12}^{\text{elast}}$ has a zero initial value for each cycle. Once the fractional derivatives are calculated using previous loading history, the optimal solution can be found within a time-interval corresponded to “unloading-reloading” or to the hysteresis loop (Fig. 6). The error function $\delta$ is computed by adding all the individual point errors inside considering interval:

$$ \delta = \frac{\sqrt{\sum_{i=1}^{N} (\gamma_{12}^{\text{test}} - \gamma_{12}^{\text{model}})^2}}{N} \qquad (19) $$

where $\gamma_{12}^{\text{test}}$ is the experimental strain determined by equation (15), $\gamma_{12}^{\text{model}}$ is the shear strain calculated by the mathematical model (16) and $N$ is a number of time-points inside the referring interval. In the further calculus, the differentiable strain function $\gamma_{12}^{\text{elast}}$ is replaced by a piecewise approximation.
<figure><img src="image_7.png" /><figcaption>Figure 6: Elastic strain referring to a hysteresis loop.</figcaption></figure>

### 3.3 Numerical evaluation of fractional model

The fractional derivative can be calculated numerically by the different algorithms. In general, the L1-algorithm (20) is used to calculate Remann-Liouville fractional derivative [35]. It has the following form:

$$D^{\alpha} f(t)_{L1} = \frac{(\Delta t)^{-\alpha}}{\Gamma(2-\alpha)} \left[ \left( \frac{1-\alpha}{N^{\alpha}} \right) f_0 + \sum_{j=0}^{N-1} (f_{N-j} - f_{N-j-1}) ((j+1)^{1-\alpha} - j^{1-\alpha}) \right] \quad (20)$$

where $\Delta t$ is the time-step, N the number of time-points and $\Gamma$ is the Gamma-function defined by (13).

The alternative forms of Remann-Liouville fractional derivatives are

$$D^{\alpha} f(t)_{M1} = \frac{d}{dt} \left( \frac{t^{1-\alpha}}{\Gamma(\alpha)} \int_{0}^{1} f\left(t\left(1-\nu^{\frac{1}{1-\alpha}}\right)\right) dv \right) \quad (21)$$

$$D^{\alpha} f(t)_{M2} = \frac{d}{dt} \left( \frac{2t^{1-\alpha}}{\Gamma(\alpha)} \int_{0}^{1} v f\left(t\left(1-\nu^{\frac{2}{1-\alpha}}\right)\right) dv \right) \quad (22)$$

The expression (21) and (22) can be easily implemented in numerical code if an analytical form of the differentiable function is known.

The Remann-Liouville fractional derivative can be calculated analytically for certain functions [34, 39]. The fractional derivative of a polynomial function $f(t) = mt^{\lambda}$ of order $\alpha$ is:

$$D^{\alpha} f(t)_{analyt} = m \frac{\Gamma(\lambda + 1)}{\Gamma(\lambda + 1 - \alpha)} t^{\lambda - \alpha} \quad (23)$$

The numerical accuracy of the proposed method (Fig. 7) is estimated on the example of the fractional derivatives of linear function $f(t) = t$ of order $\alpha = 0.5$ within the interval [0,1] with a time-step $dt=0.05$. The integrals in the M1, M2 methods are calculated by Gaussian quadrature. The interval [0,1] is divided in ten sub-integrals and five gauss-points are used in each sub-interval.

The relative error function is:

$$\delta = \frac{\sqrt{\sum_{i=1}^{N} (D^{\alpha} f(t)_{analyt} - D^{\alpha} f(t)_{approx})^2}}{\frac{N}{D^{\alpha} f(t)_{analyt}}} 100 \quad (24)$$
where $D^\alpha f(t)_{\text{analyt}}$ is an analytical fractional derivative (23), $D^\alpha f(t)_{\text{approx}}$ corresponds to the $L_1$, $M1$ or $M2$ mathematical approximations, $N$ is a number of points inside considering interval. The relative error values are presented in table 3. The biggest error is observed near zero. The $M1$-method is the most accurate and is used in further calculus.

<table><caption>Table 3: Estimation of numerical error values</caption><thead><tr><td>Approximation</td><td>Error value (%)</td></tr></thead><tbody><tr><td>L1-algorytm</td><td>0.94</td></tr><tr><td>Method 1</td><td>0.29962</td></tr><tr><td>Method 2</td><td>0.29971</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7: Linear function $f(t) = t$ and its fractional derivative of order $\alpha = 0.5$.</figcaption></figure>

## 4 RESULTS

The initial shear modulus $G_{12}^0 = 3562$ MPa and yield stress $R_0 = 15.77$ MPa are identified from experimental data (Fig. 1). The damage parameters are obtained from damage master curve (Fig. 4). The approximation of 4th order polynomial for damage evaluation is used. The polynomial coefficients and material parameters are presented in table 1. The plasticity evaluation (Fig. 5) can be characterized by the 5th order polynomial. The corresponding coefficients are presented in table 2. The fractional derivative parameters $A$, $B$, and $\alpha$ are determined by resolving an optimization problem with constraints for each loop. The parameter $A$ has a quadratic dependence from damage and the fractional derivative order $\alpha$ depends linearly on the damage evaluation (Fig. 8). Parameter $B= 0.85$ is constant.

<figure><img src="image_9.png" /><figcaption>Figure 8: Fractional derivative parameters $A$, $\alpha$ evaluation with damage.</figcaption></figure>
Taking into account the previous assumptions and the identified material parameters, the simulation of stress-strain curves is performed. The simulation results are in good agreement with experimental data for thermoset composite material (Fig. 9). The simulation results for thermoplastic are presented in Fig. 10. The significant plastic effects in matrix and the modification of the fibers orientation during the loading provide the differences between experimental and simulation data. The asymmetric hysteresis loops cannot be treated by the proposed model. To improve the numerical results, additional variables should be introduced in constitutive relations.

<figure><img src="image_10.png" /><figcaption>Figure 9: Comparison between experimental and simulated behavior of thermoset material.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Comparison between experimental and simulated behavior of thermoplastic material.</figcaption></figure>

# 5 CONCLUSION

The proposed elementary-ply model under shear loading takes into account damage, plastic strain and hysteresis mechanisms. An isotropic strain hardening is used for plasticity. The hysteresis behavior is modeled by a fractional derivative approach. Few parameters are required to represent hysteresis loops. To determine these parameters, an optimization problem has been resolved. The fractional derivatives can be simply programmed in Matlab or FORTRAN. The model validation has been demonstrated for thermoset and thermoplastic carbon fibers composite materials.

This constitutive model is able to describe behavior of unidirectional or woven composite materials under cyclic quasi-static shear loading. One of the main advantages of the proposed model is that the elasto plastic damage model [6] and strain rate dependent model [16] are classic and are completed by the fractional derivative model. The numerical implementation is simple and the computational cost is low. This model could be easily used for a wide range of applications such as the thermo-mechanical finite element simulations of a fatigue test.
REFERENCES

[1] P. Ladevèze, Sur la mécanique de l'endommagement des composites, *Journées Nationales sur les Composites - JNC 5, Paris, France, Septembre 9-11*, 1986, pp.667-683.

[2] O. Allix, P. Ladevèze, E. Le Dantec, E. Vittecoq, Damage mechanics for composite laminates under complex loading, *IUTAM/ICM Symposium on Yielding. Damage and Failure of Anisotropic Solids, Grenoble, France, August 24-27, 1987*, pp.551-569.

[3] P. Ladevèze, Sur la modélisation et le calcul des structures composites, *Journées Nationales sur les Composites - JNC 6, Paris, France, October 11-13*, 1988, pp.679-712.

[4] E. Ledantec, Contribution à la modélisation du comportement mécanique des composites stratifiés, *PhD thesis*, Université de Pierre et Marie Curie, 1989.

[5] O. Allix, P. Ladevèze, E. Le Dantec and E. Vittecoq, Damage mechanics for composite material under complex loading, *Yielding, Damage and Failure of Anisotropic Solids, EGF5* (Edited by JP Boehler), 1990, pp. 551-596.

[6] P. Ladevèze, E. Le Dantec, Damage modelling of the elementary ply for laminated composites, *Composites Science and technology*, **43**, 1992, pp. 257-267.

[7] P. Ladevèze, Towards a fracture theory, *Proceedings of the 3rd International Conference on Computational Plasticity, Part II, Cambridge, U.K.*, 1992, 1369-1400.

[8] L. Gornet, Simulation des endommagements et de la rupture dans les composites stratifiés, *PhD thesis*, École Normale Supérieure de Cachan, 1996.

[9] P. Ladeveze, O. Allix, L. Gornet, D. Leveque, L. Perret, *Damage Mechanics in Engineering Materials, Studies in Applied Mechanics*, **46**, Elsevier Publisher, 1998.

[10] C. Hochard, P.A. Aubourg, J.P. Charles, Modelling of the mechanical behavior of woven-fabric CFRP laminates up to failure, *Composites Science and Technology*, **61**, 2001, pp. 221-230.

[11] A.F. Johnson, A.K. Pickett, P. Rozycki, Computational methods for predicting impact damage in composite structures, *Composites Science and Technology*, **61** (15), 2001, pp. 2183-2192.

[12] C. Hochard, J. Payan, Design and computation of laminated composite structures, *Composites Science and Technology*, **65**, 2005, pp. 467-474.

[13] P. Rozycki, Contribution au développement de lois de comportement pour matériaux composites soumis à l'impact, *PhD thesis*, Université de Valenciennes et du Hainaut-Cambrésis, 2000.

[14] P. Rozycki, D. Coutellier, Sensibilité à la vitesse de déformation des matériaux composites unidirectionnels: Applications au verre E/époxy à 60% de fibres, *Revue Européenne des éléments Finis*, **11**, 2002, pp. 263-276.

[15] D. Coutellier, P. Rozycki, Mise en évidence de la sensibilité à la vitesse de déformation des composites unidirectionnels verre époxy. Proposition d'un modèle numérique, *Revue des composites et des matériaux avancés*. **14**, 2004, pp. 31-50.

[16] S. Marguet, P. Rozycki, L. Gornet, A rate dependent constitutive model for carbon-fiber reinforced plastic woven fabrics, *Journal Mechanics of Advanced Materials and Structures*. **14**, 2007, pp. 619-631.

[17] S. Marguet, Contribution à la modélisation du comportement mécanique des structures sandwichs soumises à l'impact, *PhD thesis*, Ecole Centrale de Nantes, 2007.

[18] M. Caputo, F. Mainardi, A new dissipation model based on memory mechanism, *Pure and Applied Geophysics*, **91**(1), 1971, pp.134–147.

[19] R.L. Bagley and P.J. Torvik, On the Fractional Calculus Model of Viscoelastic Behavior, *Journal of Rheology*, **30**(1), 1986, pp. 133–135.
[20] P.E. Rouse, The theory of the linear viscoelastic properties of dilute solutions of coiling polymers, *Journal of Chemical Physics*, **21**, 1953, pp.1272-1280.

[21] Yu.N. Rabotanov, Ravnovesie uprugoj sredi s posledstviem, *Prikladnaya matematika i mechanika*, **12**, 1948, pp.53-62 (Russian).

[22] G.W. Scott Blair, J. E. Caffyn, An application of the theory of quasiproperties to the treatment of anomalous stress-strain relations, *Philosophical Magazine*, **40**, 1949, pp. 80–94.

[23] R.L. Bagley, P.J. Torvik, A Theoretical Basis for the Application of Fractional Calculus to Viscoelasticity, *Journal of Rheology*, **27(3)**, 1983, pp. 201-210.

[24] F. Gorenflo, F. Mainardi, *Fractional Calculus: Integral and Differential Equations of Fractional Order*, Fractals and Fractional Calculus in Continuum Mechanics, pp. 223-276, Springer, Wien and New York, 1997.

[25] Chr. Friedrich, Relaxation and retardation functions of the Maxwell model with fractional derivatives, *Rheologica Acta*, **30**, 1991, pp. 151-158.

[26] H. Schiessel, A. Metzler, A. Blumen, T.F. Nonnenmacher, Generalized viscoelastic models: Their fractional equations with solutions, *Journal of Physics A: Mathematical and General*, **28**, 1995, pp. 6567-6584.

[27] W.G. Glockle and T.F. Nonnenmacher, Fractional integral operators and Fox functions in the theory of viscoelasticity, *Macromolecules*, **24**, 1991, pp. 6426-6434.

[28] M. Caputo, J.M. Carcione, Hysteresis cycles and fatigue criteria using anelastic models based on fractional derivatives, *Rheologica Acta*, **50**, 2011, pp. 107-115.

[29] M. Mateos, L. Gornet, P. Rozycki, L. Aretxabaleta, Hysteretic shear behaviour of fibre-reinforced composite laminates, *16th European Conference on Composite Materials*, Sevilla, Spain, June 23-27, 2014.

[30] M. Mateos, L. Gornet, H. Zabala, L. Aretxabaleta, P. Rozycki, P. Cartraud, Hysteretic behavior of fiber-reinforced composites, *15th European Conference on Composite Materials*, Venice, Italy, June 24-28, 2012.

[31] L. Gornet, O. Westphal, C. Burtin, J.L. Bailleul, P. Rozycki, L. Stainier, Rapid Determination of the High Cycle Fatigue Limit Curve of Carbon Fiber Epoxy Matrix Composite Laminates by Thermography Methodology: Tests and Finite Element Simulations, *Procedia Engineering*, **2013**, *66*, 2013, pp. 697-704

[32] C. Linhone, Un modèle de comportement elastoplastique endommageable pour les structures en composites stratifies, *PhD thesis*, Université de Pierre et Marie Curie, Paris 6, 1995.

[33] P. Rozycki, Loi de comportement pour composite et outils numériques en dynamique rapide, *Habilitation à Diriger des Recherches*, Ecole Centrale de Nantes, soutenance prévue septembre 2015.

[34] S.G. Samko, A.A. Kilbas, O.I. Marichev, *Fractional integrals and derevatives. Theory and applications*, Gordon and breach science publishers, Amsterdam, 1993.

[35] K.B. Oldham, J. Spanier, *The fractional calculus*, Academic Press, New York and London, 1974.

[36] A.N. Gerasimov, Obobshenie linejnih zakonov deformacii i ih prilozhenie k zadacham vnutrennego treniya, *Prikladnaya matematika i mechanika*, **12**, 1948, 251-260 (Russian).

[37] V. Volterra, Theory of Functionals and of Integral and Integro-Differential Equations, Courier Corporation, 2005

[38] V.V. Uchaikin, *Metod drobnih proizvodnih*, Artishok, Ulianovsk, 2008, (Russian).

[39] P. Cosson, Contribution à la modélisation du comportement mécanique des solides viscoélastiques par des opérateurs différentiels d'ordre non entier, *PhD thesis*, Ecole Centrale de Nantes, 1995.