# Numerical Simulation of mechanical behavior of woven composite at different strain rate by a collaborative elasto-plastio-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. Numerical Simulation of mechanical behavior of woven composite at different strain rate by a collaborative elasto-plastio-damage model with fractional derivatives. ECCOMAS 2016 European Congress on Computational Methods in Applied Sciences and Engineering, Jun 2016, Crete Island, Greece. hal-01349712

HAL Id: hal-01349712

https://hal.science/hal-01349712

Submitted on 23 Feb 2018

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.
NUMERICAL SIMULATION OF MECHANICAL BEHAVIOR OF
WOVEN COMPOSITE AT DIFFERENT STRAIN RATE BY A
COLLABORATIVE ELASTO-PLASTO-DAMAGE MODEL WITH
FRACTIONAL DERIVATIVES

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

¹ 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
e-mail: alina.krasnobrizha@ec-nantes.fr, patrick.rozycki@ec-nantes.fr, laurent.gornet@ec-nantes.fr, pascal.cosson@ec-nantes.fr

**Keywords:** Composite material, Behavior law, Fractional derivative, Hysteresis, Strain-rate sensitivity.

**Abstract.** This paper deals with a collaborative model to represent hysteresis behavior at different strain rates. The model consists of two sub-models. The first one treats the behavior during loading path. The elastic and in-elastic strains are computed as well as the in-ply damages. The strain-rate sensitivity is also taken into account. The second sub-model involves a fractional derivative approach to describe viscoelastic material response during unloading path. The hysteresis loops and strain rate sensitivity are taken into account by fractional constitutive law. Fractional model involves a few parameters which are easily identified through an optimization procedure from the experimental data. The model is validated for thermoset and thermoplastic composite materials at different strain rates.
# 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. Complex models are 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, 2, 3]. These models take into account the damage propagation and the in-elastic strains appearing commonly in the shear and transverse directions. The strain rate dependent models were developed for unidirectional composites [4, 5] and for woven fabrics [6].

The previously developed models describe different physical phenomena, such as, the elastic and in-elastic strains, the damages propagation in orthotropic directions and the strain rate sensitivity very well. However, they cannot represent the hysteretic behavior of the material during its unloading path. The fractional derivative approach is a good technique to model the viscoelastic behavior for some natural structures and modern heterogenic materials such as elastomers and polymers [7]. A significant number of works deals with fractional viscoelastic constitutive equations for different materials under various types of loading. Caputo applied fractional Zener model to represent the behavior of glass and a few metals [8]. Bagley and Torvik proposed a fractional law in frequency domain to describe the response of certain polymers and elastomers [9]. The physical sense of fractional operators is given in the works [10], [11] for different types of polymers by using Rouse molecular theory [12]. Rabotanov presented a generalized rheological model to describe the behavior of hereditary medium [13], [14]. The hysteresis cycles for a few metals under fatigue loading, were represented by Caputo, in the frequency domain [15]. Mateos [16] proposed a constitutive model for composite materials under cyclic quasi-static shear loading and this was based on the fractional derivative approach. This model also includes strain-rate dependence and thus can be applied for dynamic loading [17]. The model is also able to represent hysteresis composite behavior using a few material parameters. So, the elastic and irreversible strains, damage and strain rate effects are taken into account but not material hardening. To fill this gap, a collaborative model is developed [18] which includes the elastoplastic damage behavior law [2] with strain-rate sensitivity [4] and a fractional derivative approach. In this paper the simulation is made by the collaborative model at different strain rates.

# 2 THEORETICAL MODEL FOR COMPOSITE PLY

The constitutive model is developed within the framework of thermodynamics isothermal irreversible processes for a woven composite elementary ply under a state of plane stress. 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 behaviors are considered to be equivalent. The continuum damage mechanics theory is used to describe material degradation such as the matrix micro-cracking and the fiber/matrix debonding. Also an isotropic hardening is assumed and the viscoelastic effects are expressed by the fractional derivatives.

The material parameters are obtained from the experimental campaign based on the cyclic tensile tests. The tensile test performed on the [0°/90°] composite allows us to characterize the longitudinal (transverse) material behavior. The shear response is obtained from the tensile cyclic test on [±45°] composite. Commonly, woven composite material has a brittle elastic response in the fiber direction. The elastic and the in-elastic strain, damage propagation and hysteresis loops, as well as the strain-rate sensitivity are present during shear cyclic test. The
typical material response is illustrated for thermoset woven composite (carbon/epoxy) in the fiber direction in the Figure 1 and in the shear direction, for two different strain rates, in the Figure 2.

<figure><img src="image_1.png" /><figcaption>Figure 1. Longitudinal (transverse) traction curve for the carbon/epoxy composite</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2. Shear stress/strain curve for the carbon/epoxy composite for two different strain-rates</figcaption></figure>

The developed constitutive model is composed of two sub-models. The first one deals with the elastoplastic damage behavior during material loading and the second sub-model describes the hysteresis behavior using a fractional derivative law. Both sub-models are strain-rate dependent and can be applied for the dynamic problems.

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

$$ \rho\psi = \rho\psi(\varepsilon^e, d_i, p) \qquad (1) $$

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

## 2.1 Damage model

Following the second principal of thermodynamic, constitutive equations are deduced from the elastic strain energy of the damaged material (equivalent to the Helmholtz potential $\psi$) which has a following form:

$$ W_e^d = \frac{1}{2} \{ C_{11}^0 (1 - d_{11}) (\varepsilon_{11}^e)^2 + C_{22}^0 (1 - d_{22}) (\varepsilon_{22}^e)^2 + 2\nu_{21}^0 C_{11}^0 \varepsilon_{11}^e \varepsilon_{22}^e + G_{12}^0 (1 - d_{12}) (2\varepsilon_{12}^e)^2 \} \qquad (2) $$
where $v_{12}^0$ and $v_{21}^0$ are the Poisson's ratios, $C_{11}^0$, $C_{22}^0$ and $G_{12}^0$ are components of stiffness matrix (3).

$$ \mathbf{C}^0 = \begin{pmatrix} C_{11}^0 & v_{21}^0 C_{11}^0 & 0 \\ v_{12}^0 C_{22}^0 & C_{22}^0 & 0 \\ 0 & 0 & G_{12}^0 \end{pmatrix} \quad \text{with} \quad \begin{cases} C_{11}^0 = \frac{E_{11}^0}{1 - v_{12}^0 v_{21}^0}, & 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} \quad (3) $$

The stress-strain relation is:

$$ \sigma = \frac{\partial W_e^d}{\partial \varepsilon^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} \varepsilon_{11}^e \\ \varepsilon_{22}^e \\ \sqrt{2}\varepsilon_{12}^e \end{pmatrix} \quad (4) $$

Thermodynamic forces associated with internal variable $d_{ij}$ are defined as following:

$$ Y_{ij} = -\frac{\partial W_e^d}{\partial d_{ij}} \Rightarrow \begin{gathered} Y_{11} = \frac{1}{2} C_{11}^0 (\varepsilon_{11}^e)^2 \quad ; \quad Y_{22} = \frac{1}{2} C_{22}^0 (\varepsilon_{22}^e)^2 \quad ; \quad Y_{12} = \frac{1}{2} C_{12}^0 (2\varepsilon_{12}^e)^2 \end{gathered} \quad (5) $$

These associated thermodynamic forces characterize the damage propagation. The state of damage can only grow [1, 2] and therefore, the threshold of undamaged zone is defined as a maximal thermodynamic force for all previous time ($\tau$) up to the current time ($t$) [5]:

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

The damage variables $d_{ij}$, {$i, j = 1,2$} represent a loss of material stiffness in different orthotropic directions. In shear the damage variable is defined from the shear modulus diminution during experiment:

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

where $G_{12}^i$ is the current shear modulus associated to each unloading-loading.

Damage evaluation law is chosen as the best approximation of experimental data. Different types of functions can be used such as linear, polynomial, logarithmic, Heaviside function.

## 2.2 Plasticity modelling and damage-plasticity coupling

The experimental data shows irreversible strains appearance mainly in shear [3, 4, 5, 6]. Thus, plastic flow is considered to be blocked in fiber directions:

$$ \varepsilon_{11}^p = \varepsilon_{22}^p = 0 \quad ; \quad \varepsilon_{12}^p \neq 0 \quad (8) $$

The damage and plasticity coupling is made using the effective stress notation (9).

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

The isotropic strain hardening is assumed. The elastic domain is defined by the yield function $f$:
$$f = \frac{|\sigma_{12}|}{(1 - d_{12})} - R(p) - R_0 \qquad (10)$$

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 \quad \text{with} \quad p = \int_{0}^{\varepsilon_{12}^p} (1 - d_{12}) d\varepsilon_{12}^p \qquad (11)$$

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

## 2.3 Strain-rate sensitivity

Polymer-matrix composite materials have strain rate dependence especially in the shear direction [5, 19]. As the strain rate increases it becomes important to modify the shear modulus and the yield stress values accordingly. To describe the material response on dynamic loading the model developed by [4] is used. The material parameters are modified using following relations:

$$G_{12}^{0 \text{ dyn}} = G_{12}^{0 \text{ static}} (1 + F^{G}(\dot{\varepsilon}, \dot{\varepsilon}_{0})) \qquad (12)$$

$$R_{0}^{\text{dyn}} = R_{0}^{\text{static}} (1 + F^{R}(\dot{\varepsilon}, \dot{\varepsilon}_{0})) \qquad (13)$$

where $F^G(\dot{\varepsilon}, \dot{\varepsilon}_0)$, $F^R(\dot{\varepsilon}, \dot{\varepsilon}_0)$ functions are taken into account for the modification of shear modulus and the yield stress. These functions are dependent on the current strain rate $\dot{\varepsilon}$ and on the “threshold” strain rate $\dot{\varepsilon}_0$ from which the strain rate affects the composite behavior. Functions $F^G$, $F^R$ are determined from the experimental data using traction tests at different strain rates. They can be expressed using different forms such as linear, polynomial, power, etc.

<figure><img src="image_3.png" /><figcaption>Figure 3. Numerical and experimental stress-strain curve comparison at dynamic loading</figcaption></figure>

Damage parameters should now be recalculated from a new shear modulus value. The material parameter identification procedure is described in the next sections. Thus we can obtain the stress-strain curve (Figure 3) for the woven carbon/epoxy composite for a strain rate of 50 mm/min using this method by using the experimental data for a strain rate of 5 mm/min. The first sub-model represents elastoplastic damage behavior of composite for different strain rates. The simulation results are found to be in good agreement with the experimental data. The next step is to adjoin the hysteresis loops modelling to the current model.
## 2.4 Fractional derivative model

The hysteresis loops are associated with energy dissipation of composite under cyclic loading. Hysteresis is a hereditary phenomenon, i.e. the previous loading history has to be taken into account. To describe this viscoelastic behavior, fractional derivatives are introduced in the constitutive equation. The Riemann-Liouville fractional derivative [20] is defined as following:

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

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 (15)$$

According to elastoplastic damage model during unloading the plastic strain stays constant and the elastic strain is a linear function of time. However, a non-linearity of strain is observed in the experimental curve in the unloading/reloading path (Figure 4).

<figure><img src="image_4.png" /><figcaption>Figure 4. Total strain composition</figcaption></figure>

In order to represent strain non-linearity, fractional derivatives are introduced in the governing law in the second sub-model. The total strain within the hysteresis loop $\varepsilon_{12}^{FD}$ is defined as following:

$$\varepsilon_{12}^{FD} = A + 2BD^{\alpha}\varepsilon_{12}^{e}(t) \qquad (16)$$

where $\varepsilon_{12}^e$ is the elastic strain determined by the elastoplastic damaged model, $D^\alpha$ is the Riemann-Liouville fractional derivative (14) and $A, B$ and $\alpha$ are 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}) \varepsilon_{12}^{FD}(t) \qquad (17)$$

By substitution the equation (16) in the (17), the constative law is:

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

with $G_{12}^1 = G_{12}^0(1 - d_{12})B$.
## 2.5 Collaboration of two sub-models

The collaboration between the models is performed automatically, depending on the sign of the yield function (10) and its derivative. If $f = 0$ and $\dot{f} = 0$, the elastoplastic 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 the damage and plastic strain stays constant, and thus we switch to the fractional derivative model.

# 3 PARAMETERS IDENTIFICATION

In the following section, the identification procedure of material parameters is described. The suggested methodology of experimental identification is applied to unidirectional carbon/epoxy composite under different strain rates of loading.

## 3.1 Material characterization in fiber directions

In the fiber direction composite material has a linear brittle response. Damage of material is instantaneous and is described by the Heaviside function (19) (Figure 5). The parameters are hence easily identified and their values are presented in the Table 1 and Table 2.

$$d_{ii} = \mathcal{H} \left( \sqrt{\bar{Y}_{ii}} - \sqrt{\bar{Y}_{ii}^R} \right) \quad \text{if } d_{ii} < 1 \text{ and } \bar{Y}_{ii} < \bar{Y}_{ii}^R; \text{ otherwise } d_{ii} = 1, \{i = 1,2\} \quad (19)$$

where $\bar{Y}_{ii}^R$ is a failure-damage threshold.

<figure><img src="image_5.png" /><figcaption>Figure 5. Evaluation of the longitudinal damage function for the carbon/epoxy composite</figcaption></figure>

<table><caption>Table 1. Longitudinal (transverse) identification for the carbon/epoxy composite</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>$E_{11}^0 = E_{22}^0$</td><td>57.3 GPa</td></tr><tr><td>$\nu_{12}^0 = \nu_{21}^0$</td><td>0.07</td></tr><tr><td>$\sigma_{11}^R = \sigma_{22}^R$</td><td>687 MPa</td></tr></tbody></table>

<table><caption>Table 2. Longitudinal (transverse) damage parameters for the carbon/epoxy composite</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>$\sqrt{\bar{Y}_{11}^0} = \sqrt{\bar{Y}_{11}^R}$</td><td>2 MPa<sup>1/2</sup></td></tr></tbody></table>
## 3.2 Material characterization in shear

From the shear test (Figure 2) we had observed the strain rate sensitivity even for low values: 5 mm/min and 50 mm/min. The identified elastic parameters are presented in the Table 3. Significant difference between yield stresses for two different strain-rates can be observed.

<table><caption>Table 3. Elastic parameters in shear for the carbon/epoxy composite</caption><thead><tr><td>Parameter</td><td>Value for 5 mm/min</td><td>Value for 50 mm/min</td></tr></thead><tbody><tr><td>G<sub>12</sub><sup>0</sup></td><td>3.36 GPa</td><td>3.37 GPa</td></tr><tr><td>R<sub>0</sub></td><td>12 MPa</td><td>20.2 MPa</td></tr><tr><td>σ<sub>12</sub><sup>R</sup></td><td>80 MPa</td><td>94 MPa</td></tr></tbody></table>

The damage propagation is identified, using the linear law (20) (Figure 6), for the low strain-rate of 5 mm/min. Using these damage parameters which are identified, the cumulative plastic strain *p* is obtained. The strain hardening function *R*(*p*) is fitted by the power law (11) on the cumulative plastic strain *p* (Figure 7). The material parameters are presented in the Table 4.

$$d_{12} = \frac{\sqrt{\bar{Y}_{12}} - \sqrt{\bar{Y}_{12}^{0}}}{\sqrt{\bar{Y}_{12}^{c}}} \quad \text{if } d_{12} < 1 \text{ and } \bar{Y}_{12} < \bar{Y}_{12}^{R}; \text{ otherwise } d_{12} = 1 \qquad (20)$$

where $\bar{Y}_{12}^0$ is the initial damage threshold, $\bar{Y}_{12}^R$ is the failure-damage threshold and $\bar{Y}_{12}^c$ is the speed of damage propagation.

<table><caption>Table 4. Shear damage and plasticity parameters for the carbon/epoxy composite</caption><thead><tr><td>Parameter</td><td>Value for 5 mm/min</td></tr></thead><tbody><tr><td>√<span style="text-decoration: overline;">Y</span><sub>12</sub><sup>c</sup></td><td>2.25 MPa<sup>1/2</sup></td></tr><tr><td>√<span style="text-decoration: overline;">Y</span><sub>12</sub><sup>0</sup></td><td>0.15 MPa<sup>1/2</sup></td></tr><tr><td>√<span style="text-decoration: overline;">Y</span><sub>12</sub><sup>R</sup></td><td>0.8 MPa<sup>1/2</sup></td></tr><tr><td>β</td><td>266.6 MPa</td></tr><tr><td>k</td><td>0.36</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 6. Shear damage function of the carbon/epoxy composite</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 7. Hardening function R(p) for the carbon/epoxy composite</figcaption></figure>

### 3.3 Identification of fractional derivative model parameters

In order to determine fractional derivative model parameters, an optimization problem has been resolved. Due to specification of fractional operators, the differentiable function must be vanished at the first point of calculus. Otherwise, the fractional operator will tend to infinity. To avoid this singularity, zero initial conditions are required for the elastic strain function $\epsilon_{12}^e$. That's why we include in the computational interval the “preloading” path where $\epsilon_{12}^e = 0$. Once the fractional derivatives are calculated on “preload-unload-reload” time-interval, the optimal solution is found within the “unloading-reloading” time-interval or within the hysteresis loop (Figure 8). The objective function is expressed by relative error $\delta$:

$$ \delta = \frac{\sqrt{\sum_{i=1}^{N} (\epsilon_{12}^{\text{test}}(t_i) - \epsilon_{12}^{\text{FD}}(t_i))^2}}{N} \quad (21) $$

where $\epsilon_{12}^{test}$ is the experimental strain determined by the equation (15), $\epsilon_{12}^{FD}$ is the strain calculated by the fractional model (22) and N is the number of time-points inside the considering interval.

$$ \epsilon_{12}^{test} = \frac{\sigma_{12}^{test}}{G_{12}^{0}(1 - d_{12})} \quad (22) $$

<figure><img src="image_8.png" /><figcaption>Figure 8. Elastic strain referring to a hysteresis loop</figcaption></figure>
To implement the fractional derivative in a numerical code the M1-method [18] is used. The Riemann-Liouville fractional integral of elastic deformation $\epsilon_{12}^e$ can be expressed in alternative form as:

$$ (I_0^{1-\alpha} \epsilon_{12}^e(t))_{M1} = \frac{t^{1-\alpha}}{\Gamma(\alpha)} \int_0^1 \epsilon_{12}^e \left( t \left( 1 - v \frac{1}{1-\alpha} \right) \right) dv \quad (23) $$

The fractional derivative is calculated using central difference scheme:

$$ (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 (24) $$

The M1-method provides numerical error less than 1%. It can be easily implemented in the numerical code if analytical expression of function $\epsilon_{12}^e$ is known. Thus the elastic strain $\epsilon_{12}^e$ is approximated by a piecewise function in further calculations.

The fractional model parameters A, B, and α are determined by resolving an optimization problem for each loop. Their values stay constant within one hysteresis loop but they are varied loop by loop. We consider A, B, and α as a function of damage as the damage is constant within one hysteresis loop. The parameters A, B, and α can be approximated by linear function of damage, for the carbon/epoxy composite during a shear test with a strain rate of 5 mm/min, as the following:

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

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

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

Approximations are illustrated on the Figure 9, Figure 10 and Figure 11. Coefficients of fittings (25), (26) and (27) are presented in the Table 5.

<table><caption>Table 5. 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_A$</td><td>0.0189</td><td>$b_A$</td><td>-0.0013</td></tr><tr><td>$m_B$</td><td>1.795</td><td>$b_B$</td><td>0.8605</td></tr><tr><td>$m_{\alpha}$</td><td>0.6944</td><td>$b_{\alpha}$</td><td>-0.0035</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 9. Fractional parameter A evaluation with damage for the carbon/epoxy composite</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 10. Fractional parameter *B* evaluation with damage for the carbon/epoxy composite</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 11. Fractional parameter $\alpha$ evaluation with damage for the carbon/epoxy composite</figcaption></figure>

# 4 RESULTS

Taking into account the previous assumptions and the identified material parameters, the simulation of stress-strain curves is performed. The brittle elastic behavior in fiber directions is represented exactly in the Figure 12. The shear curve at 5 mm/min is represented by collaborative model in the Figure 13. The numerical simulation is in a good agreement with the experimental data. The in-elastic strains, damage and hysteresis loops are taken into account during the numerical simulation.

<figure><img src="image_12.png" /><figcaption>Figure 12. Experimental and numerical behavior comparison in longitudinal direction for thermoset composite</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 13. Experimental and numerical behavior comparison in shear for carbon/epoxy composite for the strain-rate of 5 mm/min</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 14. Experimental and numerical behavior comparison in shear for thermoplastic composite</figcaption></figure>

The next simulation concerns the strain-rate sensitivity of carbon epoxy woven composite which is observed from the experimental data (Figure 2). To take this into account, the material parameters are identified for strain rates of 5 mm/min and 50 mm/min. The evolutions of the parameters are represented by the linear laws. The material response in shear for a strain rate of 50 mm/min is simulated by the collaborative model from the 5 mm/min shear test. The model is able to reproduce hysteresis loops from the quasi-static test at low strain rates. It is a promising approach to model the hysteresis loops at different strain rates and to signify damage propagation.

<figure><img src="image_15.png" /><figcaption>Figure 15. Stress-strain dynamic response for a woven thermoset composite</figcaption></figure>
# 5 CONCLUSIONS

The proposed collaborative model takes into account the damage, plastic strain and viscoelastic effects such as hysteresis loops and strain rate sensitivity. The hysteresis behavior is modeled by a fractional derivative approach. Few parameters are required to represent the hysteresis loops. These parameters are determined by resolving an optimization problem. The simple implementation of M1-method is proposed for fractional derivatives. The model is validated for thermoset and thermoplastic carbon fiber woven composite materials. The strain rate sensitivity is demonstrated on the example of carbon/epoxy composite material at low strain rates.

This constitutive model is able to describe the behavior of woven composite materials under cyclic quasi-static and dynamic loading. One of the main advantages of the proposed model is that the elastoplastic damage model [2] and strain rate dependent model [4] are classical and are completed by the fractional derivative model. The numerical implementation is simple and the computational cost is low. The collaborative model is a promising approach to quantify and correlate the material damage evolution in respect of the strain rate.

# REFERENCES

[1] P. Ladevèze, Sur la mécanique de l'endommagement des composites. *Journées Nationales sur les Composites (JNC5)*, Paris, France, 1986.

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

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

[4] P. Rozycki, *Contribution au développement de lois de comportement pour matériaux composites soumis à l'impact*. PhD thesis, Valenciennes, 2000.

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

[6] S. Marguet, P. Rozycki, L. Gornet, A rate dependence constitutive model for carbon-fiber reinforced plastic woven fabrics. *Mechanics of Advanced Materials and Structures*, **14** (8), 619-631, 2007.

[7] V. V. Uchaikin, *Metod Drobnyh Proizvodnyh*. Ulianovsk: Artishok, 2008.

[8] M. Caputo, F. Mainardi, A new dissipation model based on memory mechanism. *Pure and Applied Geographic*, **91** (134), 134-147, 1971.

[9] R. Bagley, P. Torvik, On the fractional calculus model of viscoelastic behavior. *Rheology*, **30** (11), 133-155, 1986.

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

[11] J. D. Ferry, R. F. Landel, M. Williams, Extensions of the Rouse theory of viscoelastic properties to undiluted linear polymers. *Applied Physics*, **26** (4), 359-362, 1955.

[12] R. E. Rouse, The theory of the linear viscoelastic properties of dilute solutions of coiling polymers. *Chemical Physics*, **21**, 1272-1280, 1953.

[13] Y. N. Rabotnov, Ravnovesie uprugoj sredi s posledstviem. *Prikladnaya Matematika i Mechanika*, **12**, 53-62, 1948.

[14] Y. N. Rabotnov, *Elementy Naledstvennoj Mehaniki Tverdyh Tel*. Moscow: Nauka, 1977.
[15] M. Caputo, J. Carcione, Hysteresis cycles and fatigue criteria using anelastic models based on fractional derivatives. *Rheologica Acta*, **50**, 107-115, 2011.

[16] M. Mateos, L. Gornet, H. Zabala, L. Aretxabaleta, P. Rozycki, Hysteretic behaviour of fibre-reinforced composites. *15th European Conference on Composite Materials (ECCM15)*, Venice, Italy, 2011.

[17] M. H. Mateos, *Hysteretic and damage behaviour modelling of composite materials by fractional operators*. PhD thesis, Nantes, 2014.

[18] A. Krasnobrizha, *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*. PhD thesis, Nantes, 2015.

[19] A. Gilat, R. K. Goldberg, G. R. Roberts, Experimental study of strain-rate-dependent behavior of carbon/epoxy composite. *Composite Science and Technology*, **62**, 1469-1476, 2002.

[20] S. G. Samko, A. A. Kilbas, O. I. Marichev, *Fractional Integrals and Derivatives: Theory and Applications*. Amsterdam: Gordon and Breach Science, 1993.