# Hysteresis behaviour modelling of woven composite using a collaborative elastoplastic damage model with fractional derivatives

Alina Krasnobrizha, Patrick Rozycki, Laurent Gornet, Pascal Cosson

► To cite this version:

Alina Krasnobrizha, Patrick Rozycki, Laurent Gornet, Pascal Cosson. Hysteresis behaviour modelling of woven composite using a collaborative elastoplastic damage model with fractional derivatives. Composite Structures, 2016, 158, pp.101-111. 10.1016/j.compstruct.2016.09.016. hal-01580754

HAL Id: hal-01580754

https://hal.science/hal-01580754

Submitted on 2 Sep 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
# Hysteresis behaviour modelling of woven composite using a collaborative elastoplastic damage model with fractional derivatives

Alina Krasnobrizha, Patrick Rozycki, Laurent Gornet, Pascal Cosson

LUNAM Université, Ecole Centrale de Nantes, Institut de Recherche en Génie Civil et Mécanique (GeM), UMR CNRS 6183, BP 92101, 44321 Nantes Cedex 3, France

The collaborative model is a new tool to represent visco-elastoplastic damaged material behaviour including hysteresis loops for woven composites under cyclic loading. The model consists of two sub-models. The first one describes a material behaviour during loading path. The elastic and the in-elastic strains as well as the in-ply damage are computed. The second part deals with hysteresis modelling during unloading path using a fractional derivative approach. Just a few parameters are required to represent the hysteresis loops. The material strain-rate effect is also taken into account by two sub-models. The parameters identification methodology is presented and the fractional derivative approximations are detailed. The collaborative model is validated for the woven composite with thermoset and thermoplastic matrixes. The precise estimation of the dissipated energy is possible.

**Keywords:**

Composite material

Behaviour law

Hysteresis

Fractional derivative

Characterisation

## 1. Introduction

The extensive use of composite materials in the industrial applications requires a careful examination of their mechanical properties. The material behaviour must be validated experimentally for different types of loadings (static, dynamic, fatigue) depending on the material application. Then, a behaviour model should be developed to simulate structures' response as a supportive tool to the experimental method. As the composites are anisotropic heterogeneous materials, complex models are required to adequately describe their response. Some works referred by [1-6] is concerned with the composites behaviour modelling in the mesoscale. Another option is to use the micromechanics models [7,8]. Based on the in-ply damage scenarios at the microscale, a recent work [8] considers the coupling between the delamination and micro cracking. These models take into account only the damage propagation, the material hardening and the in-elastic strains' appearances. Furthermore, the polymer matrix composites have the viscoelastic behaviour due to the nature of matrix and this behaviour is evident in the material's strain-rate sensibility. The dependence of a material's strain rate sensitivity on the viscoelastic effects was first taken into account in the works [9-11] for the

unidirectional composites and in the works [12-14] for the woven composites.

However, the coupling between the cyclic loading damages and the viscoelastic effects has not been considered in the classical damage modelling [1-14]. It must be highlighted that the damage propagation is strongly linked with the material energy dissipation and, thus, with the hysteresis loops appearance. Thus, the precise modelling of hysteresis loops is needed. Therefore, in this paper we propose a collaborative model based on a fractional derivative approach.

A fractional derivative theory is a promising technique to describe the history-dependent phenomena like viscoelasticity. A significant number of fractional viscoelastic models have been developed and successfully applied to model a mechanical response of natural structures and modern heterogenic materials such as elastomers and polymers [15-17]. The classical rheological models are rewritten in terms of fractional derivative like the Kelvin-Voigt [18,19], Maxwell [20-22] and Zener [23,24] fractional models that allows to obtain creep and relaxation functions adapted for the polymers' response. Rabotnov has developed a hereditary theory where he had introduced a generalized fractional rheological model and a new class of hereditary functions in his works [25-27]. The Rabotnov's theory is widely used to describe the behaviour of the polymers, metals and concrete. Bagley and Torvik established the fractional law in the frequency domain to describe the behaviour of polymers and elastomers [28]. The link between molecular structure of polymers and fractional derivative law is given in the works [29-31] for the different types of polymers. These relations were
obtained by using the Rouse molecular theory [32]. The hysteresis loops' appearance is another viscoelastic effect which can be observed in the viscoelastic materials response and then, it can be treated using the fractional derivative approach. Caputo has described the hysteresis loop using the fractional derivatives in the frequency domain to determine a fatigue limit for a few metals [33]. Mateos [34] proposed a fractional derivative model with a few material parameters to describe the hysteresis composite behaviour in the temporal domain. This model also includes strain-rate dependence and can be applied for dynamic loading [35]. The elastic and the irreversible strains, the damage and the strain rate effects are taken into account into [34] and [35] but not the material hardening. Moreover, these phenomena that are cited are inseparable making it difficult to calculate their individual contribution in the material's energy dissipation. Hence this paper proposes a collaborative model to describe the complete composite behaviour during cyclic loading including the representation of the hysteresis phenomena during the loading-unloading paths. The proposed model is composed of the elastoplastic damage behaviour law [3] with possible strain-rate sensitivity [9] and including a newly proposed fractional derivative approach. This model allows us to distinguish the dissipation due to the material damage, plasticity and viscoelasticity.

This paper is organized as follows. In Section 2, the essential mathematical elements of the fractional calculus are introduced and the collaborative model for woven composite materials is proposed in the Section 3. The parameter identification procedure is developed in the Section 4 and the validation of the collaborative model is demonstrated in the Section 5. The Section 6 concludes this paper and notes the various engineering application of this model.

## 2. Introduction in fractional calculus

The advent of advanced materials requires new mathematical methods to describe their behaviours. One of the ways is to increase the number of parameters in the existing constitutive laws. This leads to sophisticated models with a significant number of parameters which sometimes require complex identification procedures. On the other hand, new methods may be introduced as the fractional derivative theory.

Fractional operators are an extension of the order of the integrals and derivatives to the fractional order. By their definition, the fractional operators take into account the past loading history and are suitable to describe the viscoelastic phenomena. The solution of fractional differential equations is expressed by the fractional functions [36,37] whereas the solution of ordinary differential equations is an exponential basis. Therefore, the fractional behaviour laws are widely used to describe behaviour of heterogeneous materials and fractal structures such as polymers, alloys and geological strata [27,38,39].

There are a few definitions of the fractional integrals and derivatives which are equivalent for most of the functions but they do not always lead to the same results. In this work the classical Riemann and Liouville definition (6) [40,41] is used. To deduce this definition, the formula of nth repeated integral of the continuous function $f(x)$ (Cauchy's formula) is used:

$$ \int_a^x dx \int_a^x dx \dots \int_a^x f(x) dx = \frac{1}{(n-1)!} \int_a^x (x-t)^{n-1} f(t) dt, n \in \mathbb{N} \quad (1) $$

To generalize the Eq. (1), the natural number $n$ is replaced by a real positive number $\alpha$. It is known that the extension of the factorial function to real and complex numbers is the gamma function $\Gamma(z)$ defined by:

$$ \Gamma(z) = \int_{0}^{+\infty} e^{-x} x^{z-1} dx, z \in \mathbb{R}_{+}^{*} \quad (2) $$

The gamma function $\Gamma(z)$ satisfies the following recurrence relation:

$$ \Gamma(z+1) = z\Gamma(z), \quad z \in \mathbb{R}_{+}^{*} \quad (3) $$

The fractional Riemann-Liouville integral of order $\alpha$ (4) is obtained from the generalized form of the Eq. (1) by using the property of the gamma function (3)

$$ (I_a^\alpha f)_{RL}(x) \stackrel{\text{def}}{=} \frac{1}{\Gamma(\alpha)} \int_a^x \frac{f(t)}{(x-t)^{1-\alpha}} dt, \quad \alpha > 0 \quad (4) $$

where $(I_a^\alpha f)_{RL}(x)$ is the fractional integral of order $\alpha$ of the function $f(x)$.

The fractional derivative is an inverse operator to the fractional integral. Let us consider $\alpha$ which is a positive real number and $n$ which is an integer that satisfies the inequality $n-1 < \alpha < n$. The fractional derivative of order $\alpha$ is defined as the fractional integral of order $n - \alpha$ derived $n$ times. For example, if $\alpha$ belongs to the open interval $(0, 1)$, the fractional derivative of order $\alpha$ of the function $f(x)$ is the fractional integral of order $1 - \alpha$ derived one time:

$$ (D_a^\alpha f)_{RL}(x) \stackrel{\text{def}}{=} \frac{d}{dx} (I_a^{1-\alpha} f)_{RL}(x) \quad (5) $$

The alternative form of the previous expression is:

$$ (D_a^\alpha f)_{RL}(x) \stackrel{\text{def}}{=} \frac{1}{\Gamma(1-\alpha)} \frac{d}{dx} \int_a^x \frac{f(t)}{(x-t)^\alpha} dt, \quad 0 < \alpha < 1 \quad (6) $$

where $(D_a^\alpha f)_{RL}(x)$ is the fractional derivative of order $\alpha$ in the sense of Riemann-Liouville.

### 2.1. Numerical evaluation of fractional derivatives

In order to implement the fractional derivatives in a numerical code, different algorithms exist. The L1-approximation [42] of the Riemann-Liouville fractional derivative (6) $(D_0^\alpha f(x))_{L1}$ has a following form:

$$ (D_0^\alpha f(x))_{L1} = \frac{1}{\Gamma(1-\alpha)} \frac{f(0)}{x^\alpha} + \frac{1}{\Gamma(2-\alpha)} \sum_{j=0}^{N-1} \frac{f(x-x_{j+1})-f(x-x_j)}{\Delta x} [x_{j+1}^{1-\alpha} - x_j^{1-\alpha}] \quad (7) $$

where $\Delta x$ is a step interval, $N$ is a number of interval points and $\Gamma$ is the gamma function defined by (2).

The Riemann-Liouville fractional integral (4) can be rewritten in the alternative forms [43]:

$$ (I_0^\alpha f(x))_{M1} = \frac{x^\alpha}{\Gamma(\alpha+1)} \int_0^1 f(x(1-v^{\frac{1}{\alpha}})) dv \quad (8) $$

$$ (I_0^\alpha f(x))_{M2} = \frac{2x^\alpha}{\Gamma(\alpha+1)} \int_0^1 f(x(1-v^\alpha)) dv \quad (9) $$

The formulas (8) and (9) can be easily implemented in the numerical code if an analytical expression of the function $f(x)$ is known.

According to the expression (6), the fractional derivative is the first derivative of the $M1$ (8) or $M2$ (9) integrals. The $M1$ and $M2$ fractional derivative approximations are calculated by a central difference scheme:

$$ (D_0^\alpha f(x))_{M1} = \frac{(I_0^{1-\alpha} f(x + \Delta x))_{M1} - (I_0^{1-\alpha} f(x - \Delta x))_{M1}}{2\Delta x} \quad (10) $$

$$ (D_0^\alpha f(x))_{M2} = \frac{(I_0^{1-\alpha} f(x + \Delta x))_{M2} - (I_0^{1-\alpha} f(x - \Delta x))_{M2}}{2\Delta x} \quad (11) $$

In order to estimate the numerical accuracy of the proposed approximations, the fractional derivatives of the linear function $f(x) = x$ are computed by the methods (7), (10) and (11) and compared with the analytical results (12) [36].
$$ (D_0^\alpha f)_{RL}(x) = \frac{1}{\Gamma(2-\alpha)} x^{1-\alpha} \qquad (12) $$

The order of fractional derivatives is set $\alpha = 0.3$. The computational interval $x \in (0, 1)$ is discretized by $N = 20$ points. The results are illustrated in Fig. 1.

Accuracy of the proposed methods is estimated by using the relative error $\delta_i^{relative}$ (13) and the total error $\delta$ (14) functions.

$$ \delta_i^{relative} = \frac{|D^\alpha f(x_i)_{analyt} - D^\alpha f(x_i)_{approx}|}{D^\alpha f(x_i)_{analyt}} \quad 100\% \qquad (13) $$

$$ \delta = \frac{1}{N} \sum_{i=1}^{N} \delta_{i}^{relative} \qquad (14) $$

where $D^\alpha f(x_i)_{analyt}$ is the fractional derivative computed analytically by the Eq. (12), $D^\alpha f(x_i)_{approx}$ corresponds to the $L1$, $M1$ or $M2$ approximations defined by the Eqs. (7), (10), (11) respectively.

The integrals (8) and (9) are computed using the Gaussian quadrature. The computed interval (0, 1) is divided in ten subintervals and five gauss-points are used in each sub-interval.

The results are presented in the Fig. 2 and in the Table 1. The biggest error is observed near the zero point. The total relative errors of $M1$ and $M2$ algorithms are less than 1%. The $M1$-method is the most accurate and is used in further calculations.

### 3. Theoretical model for composite ply

The material which is used in this work is a carbon woven fibre epoxy matrix composite. The traction cyclic tests have been carried out on an Instron 3369 machine with a strain rate of 5 mm/min in two different orthotropic directions: fibre and shear directions. The specimens' geometry and their dimensions for $[0]_8$ and $[45]_8$ woven laminates are shown in the Figs. 3 and 4 respectively. The nominal thickness of the samples is 1.8 mm.

The experimental result of the tensile cyclic test for the woven carbon/epoxy composite is illustrated in the Figs. 5 and 6. Subscripts 1 and 2 represent the warp and the weft directions, respectively. The composite is considered to be perfectly balanced and thus the longitudinal and transverse behaviours are considered to be equivalent. Typically, the experimental tests show an elastic brittle response in longitudinal (transverse) direction (Fig. 5). However, shear behaviour has a non-linear character. The irreversible strains, the in-ply damage (which is observed from a regular decrease of shear modulus) and the hysteresis loops are presented in the stress-strain curve (Fig. 6).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Different fractional derivative numerical approximations for a function $f(x) = x$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Relative errors comparison for the fractional derivatives of the linear function.</figcaption></figure>

<table><caption>Table 1
Total relative errors comparison for the linear function $f(x) = x$.</caption><tr><td>δ<sub>M1</sub> (%)</td><td>δ<sub>M2</sub> (%)</td><td>δ<sub>L1</sub> (%)</td></tr><tr><td>0.2611</td><td>0.2612</td><td>5.712</td></tr></table>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Specimen geometry for $[0]_8$ laminate.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Specimen geometry for $[45]_8$ laminate.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Longitudinal traction curve for the carbon/epoxy composite.</figcaption></figure>

The collaborative model is developed to represent the elasto-plastic damaged material behaviour as well as the hysteresis loops for woven composites under cyclic loading. The model is based on the collaboration of two sub-models. The first one deals with a composite behaviour during loading path. The elastic and the in-elastic strains are computed as well as the in-ply damage. The
<figure><img src="image_7.png" /><figcaption>Fig. 6. Shear stress-strain curve for the carbon/epoxy composite.</figcaption></figure>

second sub-model involves the fractional derivative approach to represent the hysteresis loops during unloading path. The mesoscale model is developed under a state of plane stress. The constitutive equations are deduced within the framework of irreversible thermodynamics processes using the local state method.

### 3.1. Thermodynamic aspects

Within the framework of irreversible thermodynamics, a Helmholtz potential $\rho\psi$ is chosen as a function of internal variables:

$$ \rho\psi = \rho\psi(\epsilon^e, d_i, p) \quad (15) $$

where $\epsilon^e$, $d_i$, $p$ are internal variables associated with elastic strain, damage in the orthotropic directions and cumulated plasticity respectively.

### 3.2. Damage modelling

Under the given loading, the damage appears in the elementary ply. It leads to the degradation of the mechanical properties and then to the material failure. The continuum damage mechanics theory is used to describe the degradation of elementary ply [44]. Using the effective stresses notation [45], the degraded stiffness matrix C is introduced as follows:

$$ C = \begin{pmatrix} (1 - d_{11})C_{11}^0 & v_{21}^0 C_{11}^0 & 0 \\ v_{12}^0 C_{22}^0 & (1 - d_{22})C_{22}^0 & 0 \\ 0 & 0 & (1 - d_{12})G_{12}^0 \end{pmatrix} \quad (16) $$ 
with $$ \begin{cases} C_{11}^0 = \frac{E_{11}^0}{1 - v_{12}^0 v_{21}^0} & \text{and } C_{22}^0 = \frac{E_{22}^0}{1 - v_{12}^0 v_{21}^0} \\ v_{21}^0 E_{11}^0 = v_{12}^0 E_{22}^0 \end{cases} $$

where $v_{12}^0$ and $v_{21}^0$ are the Poisson's ratios, $E_{11}^0$ and $E_{22}^0$ are longitudinal and transverse Young's moduli and $G_{12}^0$ is a shear modulus.

The damage internal variables ($d_{11}$, $d_{22}$, $d_{12}$) represent the material degradation in orthotropic directions. These parameters are set between 0 and 1:

* When the material is not degraded: $d_{ij} = 0$, $\{i, j = 1, 2\}$

* If there is a “complete” damage: $d_{ij} = 1$, $\{i, j = 1, 2\}$.

The constitutive equations are deduced using the second principal of thermodynamics. The elastic strain energy of the damaged material $W_e^d$ (equivalent to the Helmholtz potential $\rho\psi$) has a following form:

$$ W_e^d = \frac{1}{2} \{ C_{11}^0 (1 - d_{11}) (\epsilon_{11}^e)^2 + C_{22}^0 (1 - d_{22}) (\epsilon_{22}^e)^2 \\ + 2 v_{21}^0 C_{11}^0 \epsilon_{11}^e \epsilon_{22}^e + G_{12}^0 (1 - d_{12}) (2 \epsilon_{12}^e)^2 \} \quad (17) $$

where $\epsilon_{11}^e$, $\epsilon_{22}^e$ and $\epsilon_{12}^e$ are the component of the elastic strain vector. The stress-strain relation is:

$$ \sigma = \frac{\partial W_e^d}{\partial \epsilon^e} \Rightarrow \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sqrt{2}\sigma_{12} \end{pmatrix} = \begin{pmatrix} (1 - d_{11})C_{11}^0 & v_{21}^0 C_{11}^0 & 0 \\ v_{12}^0 C_{22}^0 & (1 - d_{22})C_{22}^0 & 0 \\ 0 & 0 & 2(1 - d_{12})G_{12}^0 \end{pmatrix} \begin{pmatrix} \epsilon_{11}^e \\ \epsilon_{22}^e \\ \sqrt{2}\epsilon_{12}^e \end{pmatrix} \quad (18) $$

where $\sigma$ is the stress vector.

The thermodynamic forces associated with internal variables $d_{ij}$ are defined as following:

$$ Y_{ij} = - \frac{\partial W_e^d}{\partial d_{ij}} \Rightarrow \begin{aligned} Y_{11} &= \frac{1}{2} C_{11}^0 (\epsilon_{11}^e)^2; & Y_{22} &= \frac{1}{2} C_{22}^0 (\epsilon_{22}^e)^2; & Y_{12} &= \frac{1}{2} G_{12}^0 (2\epsilon_{12}^e)^2 \end{aligned} \quad (19) $$

where $Y_{11}$, $Y_{22}$ and $Y_{12}$ are thermodynamic forces associated with damage internal variables $d_{11}$, $d_{22}$ and $d_{12}$.

The associated thermodynamic forces characterize the damage propagation. The state of damage can only grow [1,3,46] and therefore, the threshold of undamaged zone $\bar{Y}_{ij}$ is defined as a maximal thermodynamic force for all previous time ($\tau$) up to the current time ($t$) [12]:

$$ \bar{Y}_{ij} = \sup_{\tau \le t} (Y_{ij}(t)), \quad \{i, j = 1, 2\} \quad (20) $$

The coupling between the internal damage variables and the associated thermodynamic forces is determined by the experimental data fitting and the functions given below are frequently used for this:

Linear law

$$ d_{ij} = \frac{\sqrt{\bar{Y}_{ij}} - \sqrt{\bar{Y}_{ij}^0}}{\sqrt{\bar{Y}_c}} \quad \text{if } d_{ij} < 1 \text{ and } \bar{Y}_{ij} < \bar{Y}_{ij}^R; \text{ otherwise } d_{ij} = 1 \quad (21) $$

Heaviside function $H$ (This function is used to describe elastic brittle behaviour)

$$ d_{ij} = H \left( \sqrt{\bar{Y}_{ij}} - \sqrt{\bar{Y}_{ij}^R} \right) \quad \text{if } d_{ij} < 1 \text{ and } \bar{Y}_{ij} < \bar{Y}_{ij}^R; \text{ otherwise } d_{ij} = 1 \quad (22) $$

where $\bar{Y}_{ij}^0$ is initial damage thresholds and $\bar{Y}_{ij}^R$ is failure-damage thresholds in the orthotropic directions. The constants in these laws are the material parameters.

The Eqs. (21) and (22) are not always able to adequately represent the damage propagation especially in the composite material with thermoplastic matrix. Thus, the polynomial law is proposed [47] to increase the performance of the elastoplastic damage model:

$$ d_{ij} = \sum_{n=1}^{N} a_{ij}^{n} \left( \sqrt{\bar{Y}_{ij}} - \sqrt{\bar{Y}_{ij}^{0}} \right)^{n} \quad \text{if } d_{ij} < 1 \text{ and } \bar{Y}_{ij} < \bar{Y}_{ij}^{R}; \text{ otherwise } d_{ij} = 1 \quad (23) $$

### 3.3. Plasticity modelling and damage-plasticity coupling

The experimental data shows the irreversible strains appearing mainly in shear [6,9,12,13]. Thus, plastic flow is considered to be blocked in fibre directions:
$$ \epsilon_{11}^p = \epsilon_{22}^p = 0; \epsilon_{12}^p \neq 0 \qquad (24) $$

where $\epsilon_{11}^p, \epsilon_{22}^p$ and $\epsilon_{12}^p$ are the components of the plastic strain.

The damage and plasticity coupling is made using the effective stress notation $\tilde{\sigma}_{12}$ (25).

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

The isotropic strain hardening is assumed. The elastic domain is defined by the yield function $f$:

$$ f = |\tilde{\sigma}_{12}| - R(p) - R_0 \qquad (26) $$

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

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

where $\beta$ and $k$ are material parameters identified from the experimental data.

The material parameter identification is detailed in the next sections. The resulting stress-strain shear curve is obtained by the elastoplastic damage model [3] for the carbon/epoxy woven composite (Fig. 7). The simulation results are in a good agreement with the experimental data. However, the hysteresis loops cannot be reproduced by this model.

### 3.4. Fractional derivative model to represent hysteresis loops

The woven composites have a hysteresis behaviour for certain loading types such as a cyclic shear test. The reason of this is that the matrix is subjected a bigger part of the shear load and as a result, the viscoelastic response of the polymer matrix can be observed as well as the gradual development of microscopic damage in the matrix or along matrix-fibre interface. This leads to the appearance of hysteresis loops which is associated with material energy dissipation. Since we are dealing with a viscoelastic phenomenon, the loading history has to be taken into account. That is why the fractional derivatives are introduced in the behaviour law to provide the material history dependence.

To deduce a behaviour law during an unloading path let us consider one hysteresis loop. According to the elastoplastic damage model, the plastic strain stays constant and the elastic strain is a linear function of time. However, a non-linearity in strain is observed in the experimental curve during the unloading/reloading path (Fig. 8).

Commonly, an additional viscoelastic strain vector $\epsilon^{ve}$ is introduced to take into account the observed non-linearity. In this case,

the total strain $\epsilon_{12}^t$ is a sum of three strains: elastic $\epsilon_{12}^e$, plastic $\epsilon_{12}^p$ and viscoelastic $\epsilon_{12}^{ve}$:

$$ \epsilon_{12}^t = \epsilon_{12}^e + \epsilon_{12}^p + \epsilon_{12}^{ve} \qquad (28) $$

Generally, this approach leads to the complex models with a large number of parameters.

An alternative approach is proposed in this work. In order to model the viscoelastic effect, the fractional derivatives are introduced in the governing law. The total strain within the hysteresis loop $\epsilon_{12}^{loop}$ is defined as following:

$$ \epsilon_{12}^{loop} = A + 2BD^{\alpha} \epsilon_{12}^{e}(t) \qquad (29) $$

where $\epsilon_{12}^e$ is the elastic strain determined by the elastoplastic damage model, $D^\alpha$ is the Riemann-Liouville fractional derivative (6) and A, B and $\alpha$ are the fractional model parameters.

As the plastic flow stays constant, the stress is expressed by the elastic law:

$$ \sigma_{12}(t) = G_{12}^{0}(1 - d_{12})\epsilon_{12}^{loop}(t) \qquad (30) $$

By substitution the Eq. (29) in the Eq. (30), the constitutive law is:

$$ \sigma_{12}(t) = G_{12}^{0}(1 - d_{12})A + 2G_{12}^{1}D^{\alpha}\epsilon_{12}^{e}(t) \qquad (31) $$

with $G_{12}^1 = G_{12}^0(1 - d_{12})B$.

The Eq. (31) is a non-standard form of Kelvin-Voight fractional model. The parameters A, B and $\alpha$ is considered to be the functions of damage $d_{12}$ and consequently of the elastic strain $\epsilon_{12}^e$ taking into account that shear damage is a function of the shear elastic strain: $d_{12} = d_{12}(\epsilon_{12}^e)$. In order to signify meanings of the fractional parameters, the influence of each term in the Eq. (31) is considered.

* The fractional derivative of the elastic strain $D^\alpha \epsilon_{12}^e(t)$ provides a strain non-linearity. Hence, the hysteresis loop appears (Fig. 9). The fractional order $\alpha$ adjusts the loop size. The growth of $\alpha$ leads to the increasing of loop area and consequently the dissipation of viscoelastic energy increases.

* The factor $G_{12}^1$ governs the slope of the loop (Fig. 10). The damage level is adjusted in respect to the elastoplastic damage model.

* The term $G_{12}^0(1 - d_{12})A$ is a connecting stress that links two sub-models: the fractional and the elastoplastic damage models (Fig. 11). Thus the position of the loop is adjusted and the stress becomes a continuous function.

### 3.5. Coupling of two sub-models

In this part, the coupling between two sub-models is considered.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison of experimental and simulated behaviour by elastoplastic damage model in shear for the carbon/epoxy composite.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Total strain composition.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Graphical representation of the fractional derivative $\sigma_{12}(t) = 2D^\alpha c_{12}^\alpha(t)$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Graphical representation of the term $\sigma_{12}(t) = 2G_{12}^1 D^\alpha c_{12}^\alpha(t)$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Graphical representation of the fractional model $\sigma_{12}(t) = G_{12}^0 (t - d_{12})A + 2G_{12}^1 D^\alpha c_{12}^\alpha(t)$.</figcaption></figure>

Let us consider the typical response in shear of the carbon/ epoxy woven composite (Fig. 12). The collaboration between two sub-models is ensured by the yield function *f* (26). There are two possible cases:

1. If the plastic yield function *f* with its derivative *f'* are zero (*f* = *f'* = 0), the plastic flow and damage propagation are present in the material and the elastoplastic damage model is used. The model represents the envelope of the stress-strain curve which is marked by the blue triangles on the Fig. 12.

<figure><img src="image_13.png" /><figcaption>Fig. 12. Coupling between the elastoplastic damage and the fractional models. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

2. If the plastic yield function *f* is negative or if the yield function *f* is zero and its derivative *f'* is nonzero (*f* < 0 or *f* = 0 and *f'* < 0 or *f* = 0 and *f'* > 0), then a viscoelastic unloading is to be considered. The plasticity and damage levels remain a constant, thus we switch to the fractional derivative model. This model represents the hysteresis loops which are marked by the red circles in the Fig. 12.

### 3.6. Methodology of fractional derivatives implementation

Let us consider the “unloading-reloading” intervals on the stress-strain curve (Fig. 12) which are indicated by the red circles. One of the intervals is represented for the “time-strain” axes in the Fig. 13. The first time point of unloading path is referred to the time $t_M$ (the beginning of the hysteresis loop) and the last time point is $t_N$ (the end of the hysteresis loop). Within the time interval [$t_M, t_N$], the yield function *f* satisfies the conditions: *f* < 0 or *f* = 0 and *f'* < 0 or *f* = 0 and *f'* > 0 and the fractional model is used. Thus, following the Eq. (31), the fractional derivatives $D^\alpha c_{12}^\alpha(t)$ is calculated within the time interval [$t_M, t_N$]. Due to the mathematical particularity of fractional derivatives [43], the differentiable function $c_{12}^e(t)$ must be vanished at the first calculation point $t_M$. Otherwise, the Riemann-Liouville fractional derivative of this function will tend to infinity at $t_M$: $D^\alpha c_{12}^e(t_M) \to \infty$. In our case this condition is not satisfied: $c_{12}^e(t_M) \neq 0$. To avoid this singularity, the time interval is extended by adding the previous loading region. Hence, the new computational interval is referred to the [$t_L, t_N$] time interval (Fig. 13) where the condition $c_{12}^e(t_L) = 0$ is satisfied. This technique

<figure><img src="image_14.png" /><figcaption>Fig. 13. Computational time interval referred to the fractional derivatives calculus.</figcaption></figure>
allows to define the fractional derivative value at the first point of the hysteresis loop $D^\alpha \epsilon_{12}^e(t_M)$. In other words, to reproduce the hysteresis loop on the “unloading-reloading” time interval $[t_M, t_N]$, the fractional derivatives are computed on the “loading-unloading-reloading” interval $[t_L, t_N]$ to avoid the mathematical problems.

The fractional operators are approximated by using the M1-method expressed by the formulas (8), (10) which can be rewritten as following:

$$ (I_0^{1-\alpha} \epsilon_{12}^e(t))_{M1} = \frac{t^{1-\alpha}}{\Gamma(2-\alpha)} \int_0^1 \epsilon_{12}^e(t(1-\nu^{1-\alpha})) dv \quad (32) $$

$$ (D_0^\alpha \epsilon_{12}^e(t))_{M1} = \frac{(I_0^{1-\alpha} \epsilon_{12}^e(t + \Delta t))_{M1} - (I_0^{1-\alpha} \epsilon_{12}^e(t - \Delta t))_{M1}}{2\Delta t} \quad (33) $$

To use the formulas (32) and (33), the analytical expression of the $\epsilon_{12}^e(t)$ function has to be reconstructed by using the piecewise function.

# 4. Parameters identification

In the following sections the identification procedure for the model parameters is described using the experimental data for the carbon/epoxy woven composite (Figs. 5 and 6). The in-ply damage, the isotropic strain hardening and the hysteresis loops mod-elling are taken into consideration.

## 4.1. Materials characterization in fibres directions

The $[0^\circ/90^\circ]_8$ tensile test (Fig. 5) is used to identify the longitudinal (or transverse) material behaviour. As the woven fabrics have commonly a brittle linear response the in fibre directions, the elastic parameters such as elastic modulus $E_{11(22)}^0$, Poisson's ratio $v_{12(21)}^0$ and rupture stress $\sigma_{11(22)}^R$ are easily identified and presented in the Table 2.

### 4.1.1. Damage identification in fibre directions

Due to the brittle elastic behaviour, the material failure appears instantly. The Heaviside step function (22) is used to describe the damage propagation. The initial damage threshold coincides with the failure threshold: $\bar{Y}_{11(22)}^0 = \bar{Y}_{11(22)}^R$. Damage parameters are provided in the Table 3. The damage propagation is illustrated in the Fig. 14.

## 4.2. Materials characterization in shear

The material parameters are identified from the shear experimental curve for thermoset composite (Fig. 6). The shear undam-

<table><caption>**Table 2**
Longitudinal (or transverse) identification of the carbon/epoxy composite.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>E<sup>0</sup><sub>11</sub> = E<sup>0</sup><sub>22</sub></td><td>57.3 GPa</td></tr><tr><td>v<sup>0</sup><sub>12</sub> = v<sup>0</sup><sub>21</sub></td><td>0.07</td></tr><tr><td>σ<sup>R</sup><sub>11</sub> = σ<sup>R</sup><sub>22</sub></td><td>687 MPa</td></tr></tbody></table>

<table><caption>**Table 3**
Longitudinal (or transverse) damage parameters for the carbon/epoxy composite.</caption><thead><tr><td>Parameter</td><td>Carbon/epoxy</td></tr></thead><tbody><tr><td>√Ȳ<sup>0</sup><sub>11</sub> = √Ȳ<sup>R</sup><sub>11</sub></td><td>2 MPa<sup>1/2</sup></td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Function of the longitudinal damage evaluation for the carbon/epoxy composite.</figcaption></figure>

aged modulus $G_{12}^0$ and the yield stress $R_0$ are identified from the linear regression of the elastic part in the experimental stress-strain curve. The failure stress $\sigma_{12}^R$ corresponds to the maximum reached stress during the test. The elastic parameters in shear are summarized in the Table 4.

### 4.2.1. Damage and plastic strain evaluation in shear

Typical cyclic shear tests show a reduction in the material stiffness. To model the material degradation, the damage variable $d_{12}$ is linked with the associated thermodynamic force $\sqrt{\bar{Y}_{12}}$. The damage level is measured by the shear modulus reduction:

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

where $G_{12}^i$ is a degraded shear modulus.

Using the experimental data we state that damage grows linearly (Fig. 15) in the carbon/epoxy composite in accordance with the Eq. (21). In this case, the coefficient $\sqrt{\bar{Y}_{12}}$ is the speed of damage propagation. The damage parameters are presented in the following Table 5.

The cumulative plastic strain $p$ is obtained using the identified damage parameters. The strain hardening function $R(p)$ is defined by the Eq. (35).

$$ R(p) = |\bar{\sigma}_{12}| - R_0 \quad (35) $$

The experimental data are fitted using the power law (26). The approximation is illustrated in the Fig. 16. The identified parameters are presented in the Table 6.

## 4.3. Identification of fractional derivative model parameters

The fractional model (29) parameters A, B and $\alpha$ are identified by solving a non-linear optimization problem with constrains. The main constrain is imposed on the fractional derivative order $\alpha$: $0 < \alpha < 1$. The objective function is defined as a relative error:

<table><caption>**Table 4**
Shear parameters identification for the carbon/epoxy composite.</caption><thead><tr><td>Parameter</td><td>Carbon/epoxy</td></tr></thead><tbody><tr><td>G<sup>0</sup><sub>12</sub></td><td>3.36 GPa</td></tr><tr><td>R<sub>0</sub></td><td>12 MPa</td></tr><tr><td>σ<sup>R</sup><sub>12</sub></td><td>80 MPa</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>Fig. 15. Shear damage function for the carbon/epoxy composite.</figcaption></figure>

<table><caption>**Table 5**
Shear damage parameters for the carbon/epoxy composite.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>√Ȳ<sub>12</sub><sup>C</sup></td><td>2.25 MPa<sup>1/2</sup></td></tr><tr><td>√Ȳ<sub>12</sub><sup>O</sup></td><td>0.15 MPa<sup>1/2</sup></td></tr><tr><td>√Ȳ<sub>12</sub><sup>R</sup></td><td>0.8 MPa<sup>1/2</sup></td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Hardening function R(p) for the carbon/epoxy composite.</figcaption></figure>

<table><caption>**Table 6**
Hardening parameters for the carbon/epoxy composites.</caption><thead><tr><th>Parameter</th><th>Carbon/epoxy</th></tr></thead><tbody><tr><td>β</td><td>266.6 MPa</td></tr><tr><td>k</td><td>0.36</td></tr></tbody></table>

$$ \delta = \frac{\sqrt{\sum_{i=1}^{N} (\varepsilon_{12}^{\text{exp}}(t_i) - \varepsilon_{12}^{\text{loop}}(t_i))^2}}{N} \quad (36) $$

where $\varepsilon_{12}^{\text{exp}}$ is an experimental strain determined by the Eq. (37), $\varepsilon_{12}^{\text{loop}}$ is a shear strain calculated by the fractional model (29) and N is a number of time-points inside the considered interval.

$$ \varepsilon_{12}^{\text{exp}} = \frac{\sigma_{12}^{\text{exp}}}{G_{12}^{0}(1 - d_{12})} \quad (37) $$

Three fractional parameters A, B and $\alpha$ are defined for each hysteresis loops. Their values stay constant within the hysteresis loop but they are altered in every loop. How it has been mentioned above, the fractional parameters are the functions of damage.

For the carbon/epoxy composite, the fractional parameters are represented by linear functions of damage as following:

$$ A = m_A d_{12} + b_A \quad (38) $$

$$ B = m_B d_{12} + b_B \quad (39) $$

$$ \alpha = m_{\alpha} d_{12} + b_{\alpha} \quad (40) $$

The results are illustrated in the Figs. 17–19. The coefficients of fittings (38), (39) and (40) are collected in the Table 7.

<figure><img src="image_18.png" /><figcaption>Fig. 17. Fractional parameter A evaluation with damage for the carbon/epoxy composite.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Fractional parameter B evaluation with damage for the carbon/epoxy composite.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. Fractional parameter $\alpha$ evaluation with damage for the carbon/epoxy composite.</figcaption></figure>
<table><caption>**Table 7**
Coefficients of fractional parameters fitting for the carbon/epoxy composite.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>m<sub>A</sub></td><td>0.0189</td><td>b<sub>A</sub></td><td>-0.0013</td></tr><tr><td>m<sub>B</sub></td><td>1.795</td><td>b<sub>B</sub></td><td>0.8605</td></tr><tr><td>m<sub>z</sub></td><td>0.6944</td><td>b<sub>z</sub></td><td>-0.0035</td></tr></tbody></table>

## 5. Results

The carbon/epoxy composite behaviour is represented numerically in fibre directions by the collaborative model using the identified parameters in Fig. 20 and the results of simulation are more than satisfying since the exact replica of the elastic brittle behaviour is obtained.

The collaborative model is applied to represent the shear strain-stress curve. The results are in good agreement with the experimental data (Fig. 21). Small errors are observed at the first and the last points of the hysteresis loop. The error at the first point is related to the mathematical particularity of fractional calculus [43]. The error at the last point is related to the uncertainty of the optimal solution. Furthermore, the derivation of the piecewise function provides an inaccuracy at the inflection point (the transition from the unloading to the loading path).

The identification of parameters is also made for the thermoplastic woven composite: carbon/PA66 in shear direction. The damage propagation has a non-linear character and is expressed by the third order polynomial law (23). The power law (27) has

been chosen to fit the hardening function. The symmetrical hysteresis loops are transformed into the “bean” shaped loops due to the fibre reorientation (Fig. 22) and this increases the error in damage identification. To improve the simulation results, the reorientation of fibre angle has to be taken into account. Due to the non-linear damage law, the fractional parameters A and B also has a non-linear evolution and are approximated by the third order polynomial law. The fractional derivative order $\alpha$ has a linear dependence on damage but its values rises because of the loop size increments.

The collaborative model is applied to simulate thermoplastic composite response. The difference between the experimental and simulation results is observed towards the end of the resulting stress-strain curve (Fig. 22) when the level of damage is high. Thus the inexact damage identification provides a supplementary error in addition to the previously cited fractional calculus problems. The asymmetric hysteresis loops cannot be treated by the collaborative model. Despite the unsymmetrical loops form, their area measures correctly and the amount of dissipated energy can be determined precisely.

In conclusion we state that the collaborative model provides more than satisfactory results for the composite materials with different matrixes.

An additional simulation has been made to predict the appearance of hysteresis loops during cyclic loading. A six cycles' experimental shear test is replaced by the fictitious eight cycles loading for the thermoplastic composite (Fig. 23). The simulation is made with the previously identified parameters. The results are shown

<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Experimental and numerical behaviour comparison in longitudinal direction for the carbon/epoxy composite.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Comparison of experimental and simulated behaviour by collaborative model in shear for the carbon/epoxy composite.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 22.** Comparison of experimental and simulated behaviour by collaborative model in shear for the carbon/PA66 composite.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 23.** Fictitious cyclic loading for the carbon/PA66 composite.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. 24.** Simulation of a fictitious loading in shear for the carbon/PA66 composite.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 25.** Stress-strain dynamic response for the carbon/epoxy composite.</figcaption></figure>

in the Fig. 24. The size and the position of the additional loops are predicted exactly. It opens a lot of interesting perspectives for the future works.

The composites with polymer matrix are sensitive to the strain rate. The ability of the collaborative model to represent hysteresis loops in dynamic loading is shown with the example of thermoset composite for two different strain rates: 5 mm/min and 50 mm/min. The material response for strain rate of 50 mm/min is predicted by the collaborative model from the quasi-static (5 mm/min) test identification (Fig. 25). The collaborative model is a promising approach to predict the hysteresis loops' appearance and to characterize the damage propagation at different strain rates.

## 6. Conclusion

The collaborative model is developed to simulate a complete dynamic behaviour of woven composites. This model consists of two sub-models: the first one is the elastoplastic damage model [3] and the second is the fractional derivative model. The elastoplastic damage model is able to represent the composite material response except the hysteresis loops. To fill this gap, the fractional model is introduced in order to describe the history-depended phenomena. The fractional constitutive law is expressed by a non-standard Kelvin-Voigt model using Riemann-Liouville fractional derivatives. Just a few parameters are required to represent the hysteresis loop in this model. The fractional model is applied during the unloading path which significantly reduces the computation time. Moreover, the fractional model is an optional tool to the elastoplastic damage model and can be used only if the hysteresis modelling is required.

The collaborative model is validated for thermoset and thermoplastic carbon woven composites. Both materials have an elastic brittle behaviour in fibres directions. The simplified identification procedure is used. The material response in shear is characterized by the damage propagation, the irreversible strains appearance and the hysteresis mechanisms. The shear behaviour was simulated by the collaborative model. The results of simulation are in good agreement with the experimental data. Even if the collaborative model is not able to represent unsymmetrical hysteresis loops for thermoplastic composites, their area is computed correctly and the dissipated energy can be estimated precisely. The fibre reorientation during the test has to be taken into account to increase the model performance.

The collaborative model is also able to predict the appearance of additional hysteresis loops. It is an important result which can lead to the simplification of experimental characterization campaign. The model is able to represent the hysteresis loops at different strain rates from the quasi-static experimental data. It is a promising approach to signify damage propagation in crash analysis.

The collaborative model is a simple and efficient approach to make a link between the mechanical and the thermodynamic simulations. The amount of dissipated energy during the cyclic loading can be easily calculated. Therefore, the collaborative model can be used as a numerical support for the self-heating experimental method [48]. A method based on the self-heating tests provides fast fatigue limit identification for the composite materials compared to conventional methods (S-N curves). Using the proposed behaviour model, the thermo-mechanical finite element simulation can be developed in order to compute accurately the dissipated energy.

## References

[1] Ladevèze P. Sur la mécanique de l'endommagement des composites. In: Proceedings of the Journées Nationales sur les Composites (JNC/5), Paris, 1986.

[2] Allix O, Ladevèze P, Le Dantec E, Vittecoq E. Damage mechanics for composite material under complex loading. In: Boehler JP, editor. Yielding, damage and failure of anisotropic solids. London: EG, Mechanical Engineering Publications; 1990. p. 551-96.

[3] Ladevèze P, Le Dantec E. Damage modelling of the elementary ply for laminated composites. Compos Sci Technol 1992;43:257-67.

[4] Ladevèze P, Allix O, Gornet L, Lévêque D, Perret L. A computational damage mechanics approach for laminates: identification and comparison with experimental results. Stud Appl Mech 1998;46:481-500.

[5] Hochard C, Aubourg P, Charles J. Modelling of the mechanical behaviour of woven-fabric CRFP laminates up to failure. Compos Sci Technol 2001;61 (2):221-30.

[6] Hochard C, Payan J, Montagnier O. Design and computation of laminated composite structures. Compos Sci Technol 2005;65:467-74.

[7] Ladevèze P, Lubineau G. On a damage mesomodel for laminates: micro-meso relationships, possibilities and limits. Compos Sci Technol 2001;61 (15):2149-58.

[8] Ladevèze P, Daghia F, Abisset E, Le Mauff C. A micromechanics-based interface mesomodel for virtual testing of laminated composites. Adv Model Simul Eng Sci 2014;1(29):1-16.

[9] Rozycki P. Contribution au développement de lois de comportement pour matériaux composites soumis à l'impact [Ph.D. thesis]. University of Valenciennes; 2000.

[10] Rozycki P, Coutellier D. 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 2002;11(2-4):263-76.

[11] Coutellier D, Rozycki P. 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 2004;14 (1):31-50.

[12] Johnson AF, Pickett AK, Rozycki P. Computational methods for predicting impact damage in composite structures. Compos Sci Technol 2001;61 (15):2183-92.

[13] Marguet S, Rozyck P, Gornet L. A rate dependence constitutive model for carbon-fiber reinforced plastic woven fabrics. Mech Adv Mater Struct 2007;14 (8):619-31.

[14] Marguet S. Contribution to the modelling of the mechanical behaviour of sandwich structures subjected to dynamic loadings [Ph.D. thesis]. Ecole Centrale de Nantes; 2007.

[15] Scott Blair G, Caffyn JE. An application of the theory of quasiparticles to the treatment of anomalous stress-strain relations. Phil Mag 1949;40:80-94.
[16] Gorenflo F, Mainardi F. Fractional calculus, integral and differential equations of fractional order. In: Carpinteri A, Mainardi F, editors. Fractals and Fractional Calculus in Continuum Mechanics. Vienna: Springer Verlag; 1997. p. 223-76.

[17] Friedrich C. Relaxation and retardation functions of the Maxwell model with fractional derivatives. Rheol Acta 1991;30:151-8.

[18] Koeller R. Applications of fractional calculus to the theory of viscoelasticity. J Appl Mech 1984;51(12):299-307.

[19] Nonnenmacher T, Metzler R. On the Riemann-Liouville fractional calculus and some recent applications. Fractals 1995;3(3):557-66.

[20] Schiessel H, Blumen A. Hierarchical analogues to fractional relaxation equations. J Phys A: Math Gen 1993;26:5057-69.

[21] Alemany P, Blumen A. Dynamics in disordered systems. Prog Colloid Polym Sci 1994;96:15-21.

[22] Schiessel H, Blumen A. Mesoscopic pictures of the sol-gel transition: ladder models and fractal networks. Macromolecules 1995;28:4013-9.

[23] Heymans N. Fractal rheological models and fractional differential equations for viscoelastic behavior. Rheol Acta 1994;33(13):210-9.

[24] Caputo M, Mainardi F. A new dissipation model based on memory mechanism. Pure Appl Geograph 1971;91(1):134-47.

[25] Rabotnov YN. Ravnovesie uprugoj sredi s posledstviem. Prikladnaya Matematika i Mechanika 1948;12:53-62.

[26] Rabotnov YN. Polzuchest elementov konstrucij. Moscou: Nauka; 1966.

[27] Rabotnov YN. Elementy naledstvennoj mehaniki tverdyh tel. Moscow: Nauka; 1977.

[28] Bagley R, Torvik P. On the fractional calculus model of viscoelastic behaviour. Rheology 1986;30(11):133-55.

[29] Bagley R, Torvik P. A theoretical basis for the application of fractional calculus to viscoelasticity. Rheology 1983;27:201-10.

[30] Ferry JD, Landel RF, Williams M. Extensions of the Rouse theory of viscoelastic properties to undiluted linear polymers. Appl Phys 1955;26(4):359-62.

[31] Schiessel H, Friedrich C, Blumen A. Applications to problems in polymer physics and rheology. Appl Fract Calcul Phys 2000;331-76.

[32] Rouse RE. The theory of the linear viscoelastic properties of dilute solutions of coiling polymers. Chem Phys 1953;21:1272-80.

[33] Caputo M, Carcione J. Hysteresis cycles and fatigue criteria using anelastic models based on fractional derivatives. Rheol Acta 2011;50:107-15.

[34] Mateos MH, Gornet L, Zabala H, Aretxabaleta L, Rozycki P. Hysteretic behaviour of fiber-reinforced composites Italy. In: Proceedings of the 15th European Conference on Composite Materials (ECCM/15), Venice, 2012.

[35] Mateos MH. Hysteretic and damage behaviour modelling of composite materials by fractional operators [Ph.D. thesis]. Ecole Centrale de Nantes; 2014.

[36] Samko SG, Kilbas AA, Marichev Ol. Fractional integrals and derivatives: theory and applications. Amsterdam: Gordon and Breach Science; 1993.

[37] Kilbas AA, Srivastava HM, Trujillo JJ. Theory and applications of fractional differential equations. Amsterdam: Elsevier; 2006.

[38] Cosson P, Michon J. Identification by a non-integer order model of the mechanical behaviour of an elastomer. Chaos, Solitons Fractals 1996;7 (11):1807-24.

[39] Uchaikin VV. Metod drobnyh proizvodnyh. Ulianovsk: Artishok; 2008.

[40] Liouville J. Mémoire sur quelques questions de géometrie et de mécanique et sur Mémoire sur quelques questions de géométrie et de mécanique et sur un nouveau genre de calcul pour resoudre ces questions. Journal de l'École Royale Polytechnique 1832;13(21):1-69.

[41] Riemann B. Versuch einer allgemeinen auffassung der integration und differentiation, (14 janvier 1847). In: Teubner Verlag, editors. Bernhard Riemanns gesammelte mathematische Werke und wissenschaftlicher Nachlass. Leipzig, 1876. p. 353-66.

[42] Oldham K, Spanier J. The fractional falculus. New York and London: Academic Press; 1974.

[43] Krasnobrizha A. Modélisation des mécanismes d'hystérésis des composites tissés à l'aide d'un modèle collaboratif élasto-plastique endommageable à dérivées fractionnaires [Ph.D. thesis]. Ecole Centrale de Nantes; 2015.

[44] Ladevèze P. Inelastic strains and damage. In: Talreja R, editor. Damage mechanics of composite materials. Amstardam: Elsevier; 1994. p. 117-38.

[45] Lemaitre J, Chaboche J-L. Mecanique des matériaux molides. Paris: Dumond; 1996.

[46] Herakovich CT. Mechanics of fibrous composites. New-York: Wiley; 1998.

[47] Rozycki P. Loi de comportement pour composite et outils numériques en dynamique rapide. Habilitation thesis, Nantes, 2015.

[48] Westphal O. Analyse thermomécanique de l'endommagement en fatigue de stratifié de carbone/époxy: détermination de la limite de d'endurance à partir d'essais d'auto-échauffement [Ph.D. thesis]. Ecole Centrale de Nantes; 2014.