# 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

Alina Krasnobrizha, Patrick Rozycki, Pascal Cosson, Laurent Gornet

► To cite this version:

Alina Krasnobrizha, Patrick Rozycki, Pascal Cosson, Laurent Gornet. 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. Matériaux et Techniques, 2016, Nouvelles avancées dans le domaine des matériaux composites, 104, pp.1-14. 10.1051/mattech/2016018. hal-02443161

HAL Id: hal-02443161

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

Submitted on 14 Jul 2021

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

Alina Krasnobrizha, Patrick Rozycki, Pascal Cosson et Laurent Gornetª

**Résumé** – Ces travaux portent sur la modélisation du comportement des matériaux composites tissus carbone sous chargements cycliques. On propose un modèle collaboratif élasto-plastique endommageable visqueux qui se décompose en deux sous-modèles. Le premier est consacré à caractériser le comportement pendant une phase de chargement. Les déformations élastiques et inélastiques sont calculées ainsi que les endommagements dans le pli élémentaire. La deuxième partie du modèle est formulée à partir de dérivées fractionnaires pour décrire les boucles d'hystérésis pendant une phase de décharge/charge. La sensibilité à la vitesse de déformation peut également être prise en compte dans les deux sous modèles. La théorie mathématique fractionnaire permet de remplacer les dérivées entières par des dérivées fractionnaires dans la partie visqueuse du modèle de comportement. Ces lois fractionnaires sont utilisées pour décrire l'historique du comportement des matériaux non homogènes. Le modèle fractionnaire proposé permet de décrire parfaitement les boucles d'hystérésis. L'algorithme d'implémentation numérique du modèle ainsi que son identification sont présentés. La validité de notre modèle est illustrée sur des matériaux composites tissés à matrices thermodurcissables ou thermoplastiques.

**Mots clés** : Composite tissu / loi de comportement / dérivée fractionnaire / hystérésis / caractérisation

**Abstract** – Modelling of the hysteretic behaviour of woven composites using a collaborative elastoplastic damage model with fractional derivatives. This paper deals with the behaviour modelling of carbon fabric composite materials under cyclic loading. The proposed collaborative model consists of two sub-models. The first one describes the material behaviour during a loading path. The elastic and inelastic deformations as well as in ply damages are computed. The second part deals with fractional derivative approach to describe hysteretic behaviour during unloading and reloading paths. The strain-rate sensibility is taken into account in both sub-models. Using mathematical theory of fractional calculus, the integer order derivatives are replaced by fractional ones in constitutive viscous laws. Indeed, the fractional laws are used to describe the history dependent behaviour of non-homogeneous materials. The proposed fractional model is able to perfectly describe the hysteresis loops. A numerical algorithm of the complete model, the computation methodology of fractional derivatives and the detailed parameters identification method are proposed. The developed collaborative model is tested on composite materials behaviour with thermoset and thermoplastic matrices.

**Key words:** Woven composite / behavior law / fractional derivative / hysteresis / characterisation

## 1 Introduction

De nos jours, les matériaux composites sont utilisés pour la fabrication de grandes structures de formes complexes dans les domaines navale, aéronautique ou en-core automobile. En effet, ce sont des matériaux plus

qu'intéressants lorsque les critères de légèreté et de grande résistance deviennent des facteurs déterminants et prépondérants pour le design des structures. On peut citer par exemple, les matériaux composites stratifiés tissés carboné/époxy ou verre/époxy qui possèdent de très bonnes propriétés en fatigue en regard des matériaux métalliques.

<footnote>ª Auteur de correspondance : laurent.gornet@ec-nantes.fr</footnote>

<footnote n="1">¹ Institut de recherche en Génie Civil et Mécanique (GeM) - UMR CNRS, École Centrale de Nantes, UMR 6183, 1 rue de la Noë, BP 92101, 44321 Nantes Cedex 3, France</footnote>
C'est pour cela qu'ils sont naturellement retenus pour la réalisation de pâles d'hélicoptère ou d'éolienne [1,2]. Actuellement, on tend même à utiliser des matériaux composites avec des architectures de fibres plus complexes (tridimensionnelles) que l'on retrouve dans la réalisation des aubes de la soufflante des moteurs Leap Snecma [3,4]. Pour tous ces matériaux destinés à la réalisation de structures primaires supportant de grands efforts, la phase de validation doit inclure des simulations réalistes des essais mécaniques.

Depuis de nombreuses années, plusieurs modèles ont été développés pour les matériaux composites qui permettent de prendre en compte les déformations élastiques, inélastiques, les endommagements [5-7] et la sensibilité à la vitesse de déformation [8,9]. Cependant, ces modèles ne permettent pas la représentation des boucles d'hystérésis.

La théorie des dérivées fractionnaires a été initiallement introduite au 19ᵉ siècle [10-12]. Néanmoins, ce n'est que dans la seconde partie du 20ᵉ siècle que cet outil mathématique a été employé dans les sciences de l'ingénieur [13, 14] et en particulier dans le domaine antivibratoire automobile [15-17]. Les modèles constitutifs à base de dérivées fractionnaires ont commencé à être de plus en plus utilisés pour décrire des phénomènes physiques [18, 19]. En effet, ils sont particulièrement bien adaptés à la modélisation des phénomènes héréditaires viscoélastiques comme les boucles d'hystérésis [20, 21].

Afin d'enrichir les premiers modèles de matériaux composites, nous nous sommes proposés d'intégrer les avantages des lois fractionnaires aux lois de comportement pour matériaux composites précédemment cités. Le modèle développé dit « collaboratif » se compose de deux sous-modèles : le premier est basé sur les travaux initiallement proposés par P. Ladevèze et ses collaborateurs au LMT de Cachan [5]. Ce modèle de type élasto-plastique endommageable avait déjà été modifié afin de prendre en compte la sensibilité à la vitesse de déformation [8, 9, 22]. Il agira uniquement lors de l'évolution de la surface de charge correspondant à la description des déformations irréversibles. Le second sous modèle est constitué d'une loi de comportement écrite à l'aide de dérivées fractionnaires qui permet d'intégrer les effets visqueux mis en évidence par l'existence de boucles d'hystérésis. Ce second modèle n'agira que lorsque le matériau subit une décharge et tant que l'état des déformations irréversibles précédent n'est pas atteint de nouveau.

Ce type de modèle collaboratif peut s'avérer être intéressant et avantageux notamment en ce qui concerne la détermination de la limite de fatigue des matériaux composites via des essais d'auto-échauffement. En effet, la détermination expérimentale des limites de fatigue d'un matériau composite requiert un grand nombre d'essais sur plusieurs stratifications de manière à identifier les différentes composantes du tenseur des contraintes. Cette détermination sur une stratification donnée, est faite habituellement grâce à la construction de courbes, appelées courbes de Wöhler : on applique sur une éprouvette des chargements cycliques pendant une durée assez conséquente (habituellement plusieurs semaines) afin que

les chargements se rapprochent de la limite de fatigue de l'éprouvette. Pour estimer en quelques heures seulement cette limite de fatigue, une approche alternative a été développée en se basant sur l'effet d'auto-échauffement du matériau composite [23]. Ce type d'approche a été initialement proposée pour évaluer la limite de fatigue de matériaux isotropes comme les alliages métalliques [24], les élastomères [25] et plus récemment sur des composites à fibres courtes [26]. La méthode consiste à appliquer à une éprouvette une séquence de blocs de chargement cyclique en contrainte et à mesurer la température stabilisée de l'échantillon sur chaque bloc de chargement imposé. Lorsque la valeur de la température stabilisée augmente de manière significative, on considère que la limite en fatigue est atteinte [27]. Des études tomographiques combinées à des essais traditionnels de fatigue ont permis de mettre au point la méthode. Les résultats expérimentaux de la méthode auto-échauffement sont en bon accord avec ceux issus des essais classiques de fatigue conduisant à la construction de la courbe de Wöhler [28]. Des simulations numériques ont aussi été effectuées en utilisant le modèle élasto-plastique endommageable [5] : elles ont montré que la dissipation des boucles d'hystérésis devait être prise en compte pour retrouver la température de l'éprouvette observée lors des différents blocs de chargement de la méthode « d'auto-échauffement » et c'est ce qui justifie l'introduction du modèle collaboratif.

## 2 Introduction au calcul fractionnaire

Dans le cadre de ce travail, on s'intéresse au développement d'un modèle de comportement dédié aux matériaux composites soumis à des chargements cycliques. Les modèles de comportement classiques [5, 22] ne permettent pas de représenter les boucles d'hystérésis associés aux mécanismes de viscoélasticité du matériau. Afin de décrire de tels phénomènes dissipatifs, il est important que ces modèles de comportement puissent prendre en compte l'histoire des sollicitations. Nous proposons donc d'introduire les dérivées fractionnaires dans le développement de lois afin de représenter les effets de mémoire dans le comportement mécanique des matériaux.

### 2.1 Les dérivées fractionnaires

Le calcul fractionnaire est un domaine des mathématiques qui a pour objet d'étendre les définitions des intégrales et des dérivées traditionnelles à des ordres non entiers. L'intégrale et la dérivée fractionnaires représentent respectivement la généralisation à des ordres non entiers de l'intégrale et de la dérivée, tout comme la fonction puissance à exposant réel correspond à la « prolongation » de la fonction puissance à exposant entier. Plusieurs définitions ont été proposées pour la dérivation non entière [29-33]. Il faut cependant noter que ces définitions ne mènent pas toujours à des résultats identiques mais sont globalement équivalentes pour un grand nombre de fonctions. Dans ce travail,
une définition classique de l'intégrale et de la dérivée fractionnaire proposée par Riemann-Liouville [29, 30] est utilisée.

L'intégration successive de la fonction $f(x)$ se définit par la formule de Cauchy qui s'écrit sous la forme suivante :

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

En généralisant la relation (1), l'intégrale fractionnaire d'ordre $\alpha$ au sens de Riemann-Liouville de la fonction $f(x)$ est définie par la formule (2) en utilisant la fonction Gamma (3) et sa propriété principale (4).

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

où $(I_a^\alpha f)_{RL}$ est l'intégrale fractionnaire d'ordre $\alpha$, et $\Gamma$ est la fonction Gamma :

$$ \Gamma(z) = \int_{0}^{+\infty} e^{-x} x^{z-1} dx, z \in R_{+}^{*} \quad (3) $$

$$ \Gamma(z+1) = z\Gamma(z), \quad z \in R_{+}^{*} \quad (4) $$

La dérivée fractionnaire est introduite comme une opération inverse à l'intégration fractionnaire. Soit $\alpha$ un réel strictement positif et $n$ l'entier vérifiant $n-1 < \alpha < n$. La dérivée d'ordre $\alpha$ s'obtient en dérivant $n$ fois la fonction $f$ intégrée à l'ordre $n-\alpha$. Par exemple, pour $0 < \alpha < 1$, la dérivée fractionnaire d'ordre $\alpha$ de $f$ s'obtient en dérivant 1 fois la fonction intégrée à l'ordre $1-\alpha$, ce qui peut se noter :

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

Ce qui amène à la forme alternative suivante :

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

La définition (6) correspond à la dérivée fractionnaire au sens de Riemann-Liouville.

## 2.2 Implémentation numérique des opérateurs fractionnaires

Dans cette partie, on présente les schémas numériques des intégrales et des dérivées fractionnaires au sens de Riemann-Liouville.

Soient $a = 0$ et $0 < \alpha < 1$. L'approximation des dérivées fractionnaires d'ordre $\alpha$ au sens de Riemann-Liouville (6) par l'algorithme L1 [13] s'écrit :

$$ (D_0^\alpha f(x))_{L1} = \frac{1}{\Gamma(1-\alpha)} \frac{f(0)}{x^\alpha} + \frac{1}{\Gamma(2-\alpha)} \times \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) $$

où $N$ est le nombre de points, $\Delta x = \frac{x}{N}$ et $x_j = j\Delta x$ avec $j = 0, 1, \dots, N-1$.

Soit $a = 0$. L'intégrale fractionnaire au sens de Riemann-Liouville (6) peut être réécrite sous les formes alternatives suivantes :

$$ (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^{\frac{2}{\alpha}}))dv \quad (9) $$

L'implémentation numérique des intégrales (8), (9) est simple si la forme analytique intégrable de la fonction $f$ est connue. L'approximation des intégrales est alors possible à l'aide par exemple de la méthode de quadrature de Gauss. L'intérêt de ces deux formulations consiste dans le fait que toutes les intégrales sont à prendre sur un intervalle fixe $[0, 1]$. Les abscisses des points d'intégration ainsi que leurs poids sont définis en début de calcul. Il suffit uniquement, pour chaque nouvelle valeur de $x$ où il faut calculer $I_0^\alpha f(x)$, de calculer la valeur de l'intégrant à chaque point de Gauss.

Pour obtenir les approximations des dérivées fractionnaires en accord avec (5), la dérivée d'ordre 1 doit être calculée à partir des intégrales $M1$ (8), $M2$ (9). Les approximations des dérivées fractionnaires $M1$, $M2$ sont obtenues en utilisant le schéma des différences centrées :

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

À titre d'exemple d'utilisation de l'algorithme, on considère la dérivée fractionnaire d'une fonction linéaire $f(x) = x$. Cet exemple est important car on choisit régulièrement d'interpoler les courbes expérimentales de manière linéaire par morceau. La dérivée fractionnaire d'ordre $\alpha$ de cette fonction se calcule par la formule suivante :

$$ D_0^\alpha f(x) = \frac{x^{1-\alpha}}{\Gamma(2-\alpha)} \quad (12) $$

Soit $\alpha = 0.3$ et $N = 20$ points de contrôle sur le domaine. Les approximations des dérivées fractionnaires (7), (10) et (11) sont comparées avec la solution analytique (12) sur la Figure 1.

Afin d'estimer la précision de chaque méthode les valeurs des erreurs relatives (13) sont calculées en chaque point. Les résultats sont présentés sur la Figure 2. La somme des erreurs numériques $\delta$ est calculée en suivant la formule (14). Les résultats sont présentés dans Tableau 1.

$$ \delta_i^{\text{relative}} = \frac{|D_0^\alpha f(t_i)_{\text{analyt}} - D_0^\alpha f(t_i)_{\text{approx}}|}{D_0^\alpha f(t_i)_{\text{analyt}}} \quad 100\% \quad (13) $$

$$ \delta = \frac{1}{N} \sum_{i=1}^{N} \delta_i^{\text{relative}} \quad (14) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparaison des approximations numériques des dérivées fractionnaires pour une fonction $f(x) = x$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Erreurs relatives des approximations numériques des dérivées fractionnaires pour la fonction $f(x) = x$.</figcaption></figure>

<caption>Fig. 2. Relative errors of fractional derivative numerical approximations for the function $f(x) = x$.</caption>

Dans les équations (13) $D^\alpha f(t)_{\text{analyt}}$ représente la dérivée fractionnaire analytique et $D^\alpha f(t)_{\text{approx}}$ correspond respectivement aux approximations numériques (7), (10) et (11). La méthode de Gauss a été utilisée pour calculer les intégrales dans les approximations $M1$ (10) et $M2$ (11). L'intervalle $[0, 1]$ est divisé en dix sous-intervalles et cinq points de Gauss sont utilisés dans chaque sous-intervalle. On constate que les plus grandes erreurs sont observées autour de la valeur initiale zéro. Les algorithmes $M1$ et $M2$ sont les plus précis. Ils donnent des valeurs d'erreurs relatives inférieures à 1 %. Dans la suite de ce travail on utilisera l'algorithme $M1$ (10) pour calculer la partie fractionnaire du comportement.

### 3 Modèle collaboratif élasto-plastique endommageable à dérivées fractionnaires

Dans le cadre de ce travail, nous proposons une nouvelle méthode pour représenter le comportement élasto-

<caption>Tableau 1. Comparaison des sommes des erreurs relatives pour la fonction $f(x) = x$.</caption>

<table><caption>Table 1. Comparison of relative errors sum for the linear function $f(x) = x$.</caption><thead><tr><th>δ<sub>M1</sub> (%)</th><th>δ<sub>M2</sub> (%)</th><th>δ<sub>L1</sub> (%)</th></tr></thead><tbody><tr><td>0,261</td><td>0,261</td><td>5,712</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Courbe de traction longitudinale pour un matériau composite carbone/époxy.</figcaption></figure>

<caption>Fig. 3. Longitudinal tensile curve for a carbon/epoxy composite.</caption>

plastique endommageable avec description des boucles d'hystérésis pour les matériaux composites sous un chargement cyclique. L'approche proposée se base sur la collaboration de deux modèles : le premier consiste à caractériser le comportement du composite pendant une phase de chargement. Les déformations élastiques et inélastiques sont déterminées ainsi que les endommagements du pli élémentaire. La deuxième partie du modèle utilise une approche par dérivées fractionnaires pour décrire le comportement visqueux pendant une phase de « décharge – recharge ». Notre modèle collaboratif est développé pour un matériau composite tissé à l'échelle mésoscopique d'un pli élémentaire sous l'hypothèse de contraintes planes [34, 35]. Ce modèle est écrit dans le cadre thermodynamique des processus irréversibles en utilisant la méthode d'état local [36].

Les lois de comportement établies permettent de transcrire les comportements longitudinaux, transversaux et de cisaillement. Les coefficients 1 et 2 représentent respectivement les directions chaine et trame. Le composite tissé est considéré parfaitement équilibré donc le comportement suivant la direction longitudinale est identique à celui de la direction transverse. Les résultats expérimentaux à propos d'un tissu carbone/époxy, ont montré des comportements longitudinaux et transversaux de type linéaire fragile (Fig. 3) alors que le comportement en cisaillement est non-linéaire de manière globale puisqu'on constate l'apparition d'un endommagement (chute régulière du module de cisaillement), la présence de déformations irréversibles et de boucles d'hystérésis en cas de chargement cyclique (Fig. 4). On peut remarquer que dans le cadre d'une nature de fibres différentes comme
<figure><img src="image_5.png" /><figcaption>Fig. 4. Courbe contrainte/déformation de cisaillement pour un matériau composite carbone/époxy.</figcaption></figure>

Fig. 4. Shear stress/strain curve for a carbon/epoxy composite.

le verre, les directions longitudinales et transversales ont des comportements linéaires endommageables [8, 37].

## 3.1 Aspect thermodynamique

Dans le cadre de la théorie thermodynamique irréversible, on choisit le potentiel de Helmholtz $\rho\psi$ qui dépend des variables internes [36] :

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

où $\epsilon^e$, $d_i$, $p$ sont des variables internes associées respectivement aux déformations élastiques, endommagements et à la déformation plastique cumulée.

## 3.2 Modélisation de l'endommagement

Sous une sollicitation donnée, l'endommagement peut se propager dans un pli élémentaire. Ce fait conduit à la dégradation des propriétés mécaniques puis à la rupture du matériau. Afin de décrire la dégradation du matériau, on utilise la mécanique d'endommagement des milieux continus à l'échelle mésoscopique [38]. L'état d'endommagement est supposé uniforme dans l'épaisseur du pli élémentaire. En utilisant la notation des contraintes effectives [36], on introduit la matrice de rigidité dégradée du pli élémentaire :

$$ \text{avec} \quad \left\{ \begin{array}{l} C = \begin{pmatrix} (1-d_{11})C_{11}^0 & \nu_{21}^0 C_{11}^0 & 0 \\ \nu_{12}^0 C_{22}^0 & (1-d_{22})C_{22}^0 & 0 \\ 0 & 0 & (1-d_{12})G_{12}^0 \end{pmatrix} \\ C_{11}^0 = \frac{E_{11}^0}{1-\nu_{12}^0\nu_{21}^0} \quad \text{et} \quad C_{22}^0 = \frac{E_{22}^0}{1-\nu_{12}^0\nu_{21}^0} \\ \nu_{21}^0 E_{11}^0 = \nu_{12}^0 E_{22}^0 \end{array} \right. \quad (16) $$

où $\nu_{12}^0$ et $\nu_{21}^0$ sont les coefficients de Poisson, $E_1^0$ et $E_2^0$ sont les modules d'Young longitudinal et transversal $G_{12}^0$ représente le module de cisaillement.

Les paramètres d'endommagement ($d_{11}d_{22}$, $d_{12}$) sont des variables internes qui représentent l'état dégradé du matériau dans les directions d'orthotropie. Ces paramètres sont définis entre 0 et 1 :

- Quand le matériau n'est pas dégradé on a : $d_{ij} = 0, \{i, j = 1, 2\}$

- S'il y a un endommagement « complet », on a : $d_{ij} = 1, \{i, j = 1, 2\}$.

Dans le cadre de la théorie thermodynamique, on introduit l'énergie de déformation élastique du matériau endommagé $W_e^d$ (qui est équivalente au potentiel thermodynamique d'énergie libre de Helmholtz $\psi$) :

$$ 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\nu_{21}^0 C_{11}^0 \epsilon_{11}^e \epsilon_{22}^e + G_{12}^0 (1-d_{12}) (2\epsilon_{12}^e)^2 \} \quad (17) $$

La loi de comportement mécanique d'un pli élémentaire du composite tissé est définie en utilisant les relations thermodynamiques :

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

Les forces thermodynamiques sont définies quant à elles par :

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

Les forces thermodynamiques associées caractérisent la propagation d'endommagement. L'état d'endommagement ne peut que croître [5, 34, 39], par conséquent les forces thermodynamiques sont bornées à leur valeurs maximales pour toute l'histoire du chargement [37] :

$$ \bar{Y}_{ij} = \sup_{\tau \le t} (Y_{ij}(t)) \quad (20) $$

L'évolution des variables d'endommagement est choisie à partir des observations expérimentales afin de représenter au mieux les évolutions identifiées. Les fonctions suivantes
sont souvent utilisées :

– Loi linéaire [5] :

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

– Loi polynomiale [22] :

$$d_{ij} = \sum_{n=1}^{N} a_{ij}^{n} (\sqrt{\bar{Y}_{ij}} - \sqrt{\bar{Y}_{ij}^{0}})^{n} \quad \text{si } d_{ij} < 1 \\ \text{et } \bar{Y}_{ij} < \bar{Y}_{ij}^{R}; \quad \text{sinon } d_{ij} = 1 \quad (22)$$

– Fonction de Heaviside (cette fonction est utilisée pour décrire le comportement élastique fragile) :

$$d_{ij} = H(\sqrt{\bar{Y}_{ij}} - \sqrt{\bar{Y}_{ij}^R}) \quad \text{si } d_{ij} < 1 \\ \text{et } \bar{Y}_{ij} < \bar{Y}_{ij}^R; \quad \text{sinon } d_{ij} = 1 \quad (23)$$

où $\bar{Y}_{ij}^0$ sont les taux d'énergie initiaux d'endommagement et $\bar{Y}_{ij}^R$ représentent les taux d'énergie à rupture relatifs aux directions longitudinale, transversale et de cisaillement.

## 3.3 Plasticité et couplage avec l'endommagement

Lors des observations expérimentales liées au comportement mécanique de matériaux composites tissés [9, 37, 40], on remarque l'apparition des déformations irréversibles principalement dans la direction de cisaillement. Pour satisfaire aux observations expérimentales, on considère pour les composites tissés que la rigidité des fibres ainsi que le tissage bloque l'écoulement plastique dans les directions longitudinales et transversales :

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

Remarque : dans le cas des matériaux composite carbone époxy 3D à architecture complexe cette hypothèse doit être remise en cause [4].

Le couplage en cisaillement de l'endommagement et de la plasticité est réalisé par le biais de la contrainte effective [36] :

$$\sigma_{12} = \frac{\sigma_{12}}{(1 - d_{12})} \quad (25)$$

En considérant un écrouissage isotrope, la fonction critère de plasticité $f$ est définie de la manière suivante :

$$f = \frac{|\sigma_{12}|}{(1 - d_{12})} - R(p) - R_0 \quad (26)$$

où $R_0$ est la limite élastique, $R(p)$ est la fonction d'écrouissage. Cette dernière fonction, $R(p)$, doit être choisie en fonction des observations expérimentales et est souvent prise sous une forme de loi puissance :

$$R(p) = \beta p^m \quad (27)$$

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparaison numérique/expérimental pour un essai en cisaillement sur un composite carbone/époxy.</figcaption></figure>

<caption>Fig. 5. Experimental/numerical behavior comparison in shear for a carbon/epoxy composite.</caption>

La déformation plastique cumulée s'exprimera sous la forme suivante :

$$p = \int_{0}^{\epsilon_{12}^{p}} (1 - d_{12}) d\epsilon_{12} \quad (28)$$

En utilisant les paramètres identifiés à la section suivante, la simulation de comportement du matériau composite carbone/époxy par le modèle élasto-plastique endommageable est effectuée en cisaillement (Fig. 5). Ce modèle permet de représenter correctement l'enveloppe de la charge, les déformations irréversibles et l'endommagement. Par contre, il ne peut représenter les boucles d'hystérésis consécutives à une décharge/recharge.

## 3.4 Modèle à dérivée fractionnaire

L'idée principale pour la modélisation des boucles d'hystérésis est de représenter les déformations provenant de celles-ci à l'aide de dérivées fractionnaires des déformations élastiques. Considérons une décharge-recharge qui correspond à une boucle d'hystérésis (Fig. 6). Conformément au modèle élasto-plastique endommageable [5], la déformation plastique reste constante et la déformation élastique est une fonction linéaire. Par contre, une non-linéarité de déformation est observée pendant la décharge-recharge.

De manière classique, pour prendre en compte la non-linéarité observée, une déformation viscoélastique $\epsilon^{ve}$ peut être ajoutée (Fig. 6). Dans ce cas, la déformation totale est représentée comme une somme de trois déformations : élastique, plastique et viscoélastique :

$$\epsilon_{12}^{t} = \epsilon_{12}^{e} + \epsilon_{12}^{p} + \epsilon_{12}^{ve} \quad (29)$$

Par contre, cette approche conduit à un modèle complexe associé à un grand nombre de paramètres dans le cadre standard de développement de lois de comportement. Nous proposons une approche alternative : afin de
<figure><img src="image_7.png" /><figcaption>Fig. 6. Décomposition des déformations.
Fig. 6. Strains separation.</figcaption></figure>

modéliser la non-linéarité, on exprime la déformation totale dans la boucle comme une fonction non-linéaire de déformation élastique grâce au calcul fractionnaire. La déformation dans la boucle $\epsilon_{12}^{\text{boucle}}$ est donc modélisée par une fonction à dérivée fractionnaire de la déformation élastique. Elle s'écrit :

$$ \epsilon_{12}^{\text{boucle}} = A + 2BD^{\alpha} \epsilon_{12}^{e}(t) \quad (30) $$

où $\epsilon_{12}^e$ est la déformation élastique déterminée avec le modèle élasto-plastique endommageable, $D^\alpha$ est la dérivée fractionnaire au sens de Riemann-Liouville (6) et A, B et $\alpha$ sont les paramètres du modèle fractionnaire qui seront définis dans la partie suivante.

Comme la plasticité est constante, la contrainte s'écrit sous la forme d'une loi élastique :

$$ \sigma_{12}(t) = G_{12}^{0}(1 - d_{12}) \epsilon_{12}^{\text{boucle}}(t) \quad (31) $$

En substituant l'Équation (30) dans (31) le modèle devient :

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

L'équation (32) est une forme non-standard du modèle de Kelvin-Voigt fractionnaire. On considère que A est une fonction de l'endommagement $d_{12}$ ou de la déformation élastique $\epsilon_{12}^e$ en prenant en compte que $d_{12} = d_{12}(\epsilon_{12}^e)$.

Pour essayer de comprendre ce que représentent les paramètres de la loi fractionnaire (32), on se propose d'étudier l'influence de ceux-ci sur la réponse du matériau :

- La dérivée fractionnaire de la déformation élastique $D^\alpha \epsilon_{12}^e(t)$ transforme la déformation linéaire en non-linéaire, ce qui fait apparaître la boucle d'hystérésis (Fig. 7). L'ordre de dérivée non-entier $\alpha$ représente l'ouverture de la boucle. L'augmentation de $\alpha$ détermine la croissance de l'aire de la boucle et en conséquence la dissipation de l'énergie viscoélastique.

- Le multiplicateur $G_{12}^1$ corrige la pente de la boucle en accord avec les données expérimentales (Fig. 8). Grâce

<figure><img src="image_8.png" /><figcaption>Fig. 7. Représentation graphique de la dérivée fractionnaire $\sigma_{12}(t) = 2D^\alpha \epsilon_{12}^e(t)$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Représentation graphique du terme $\sigma_{12}(t) = 2G_{12}^1 D^\alpha \epsilon_{12}^e(t)$.</figcaption></figure>

<caption>Fig. 8. Graphic representation of the term $\sigma_{12}(t) = 2G_{12}^1 D^\alpha \epsilon_{12}^e(t)$.</caption>

à ce coefficient, le niveau d'endommagement est corrigé à la valeur provenant du modèle élasto-plastique endommageable.

- Le terme $G_{12}^0(1 - d_{12})A$ correspond à une contrainte de raccordement qui lie les deux parties du modèle collaboratif (modèle fractionnaire et modèle élasto-plastique endommageable classique) (Fig. 9).

## 3.5 Couplage de deux modèles

À titre d'exemple, on représente la réponse classique d'un matériau composite carbone/époxy tissé en cisaillement plan sur la Figure 10. La collaboration entre les deux modèles est assurée via la fonction critère $f$ (26). Deux cas sont possibles :

- Si la fonction critère $f$ et sa dérivée $\dot{f}$ sont nulles ($f = \dot{f} = 0$), on considère le développement d'un
<figure><img src="image_10.png" /><figcaption>Fig. 9. Représentation graphique du modèle fractionnaire $\sigma_{12}(t) = G_{12}^0(1 - d_{12})A + 2G_{12}^1 D^\alpha \epsilon_{12}^e(t)$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Couplage entre le modèle élasto-plastique endommageable et le modèle fractionnaire.</figcaption></figure>

<caption>Fig. 10. Coupling between elastoplastic damage and fractional models.</caption>

écoulement plastique et la propagation de l'endommagement. Le modèle élasto-plastique endommageable est alors utilisé. C'est ce qui est illustré sur la Figure 10 par des triangles bleus.

Si la fonction critère $f$ est négative ou si la fonction critère $f$ est nulle et sa dérivée $\dot{f}$ est inférieure à zéro ($f < 0$ ou $f = 0$ et $\dot{f} < 0$), on considère alors une décharge viscoélastique fractionnaire. Les niveaux de plasticité et d'endommagement demeurent alors constants et égaux aux derniers calculs provenant du modèle élasto-plastique endommageable. C'est ce qui est illustré sur la Figure 10 par des cercles rouges.

## 3.6 Méthodologie de calcul numérique d'une boucle d'hystérésis

Dans cette partie, nous proposons de nous attarder sur la méthodologie de calcul d'une boucle d'hystérésis à l'aide du modèle fractionnaire. Pour clarifier la démarche,

<figure><img src="image_12.png" /><figcaption>Fig. 11. Portion de déformation pour l'identification des paramètres de dérivées fractionnaires.</figcaption></figure>

<caption>Fig. 11. Deformation part for fractional derivative parameters identification.</caption>

on considère un des intervalles « décharge –recharge » sur la courbe « contrainte-déformation » de la Figure 10 (qui est marqué par les cercles rouges). Cette zone est reproduite dans le domaine « déformation-temps » sur la Figure 11. On considère que le premier point de la décharge correspond à l'instant du temps $t_M$ (début de la boucle d'hystérésis) et le dernier point est noté par $t_N$. Sur l'intervalle temporel $[t_M, t_N]$, la fonction critère $f$ satisfait en décharge les conditions : $f < 0$ ou $f = 0$ et $\dot{f} < 0$. C'est donc le modèle fractionnaire qui prend le relais dans la modélisation du comportement.

Dans ce modèle (32), la dérivée d'ordre non-entier de la déformation élastique $D^\alpha \epsilon_{12}^e(t)$ doit être calculée sur l'intervalle temporel $[t_M, t_N]$. La particularité mathématique du calcul fractionnaire est que la fonction que l'on dérive doit être nulle à l'origine et malheureusement cette condition n'est pas vérifiée dans notre cas : $\epsilon_{12}^e(t_M) \neq 0$. Sans la précédente condition, la dérivée non-entière au sens de Riemann-Liouville tend vers l'infini au premier point de calcul : $D^\alpha \epsilon_{12}^e(t_M) \to \infty$. Afin d'éviter cette singularité mathématique, on propose d'inclure la déformation de la charge précédente ce qui revient à étendre l'intervalle temporel de calcul de la dérivée fractionnaire et de prendre comme support maintenant $[t_L, t_N]$ (Fig. 11, la déformation viscoélastique expérimentale et la déformation viscoélastique simulée sont quasiment superposées sur $[t_L, t_N]$). Dans ce nouvel intervalle, la condition $\epsilon_{12}^e(t_L) = 0$ est dorénavant satisfaite. Cette technique nous permet ainsi de définir la valeur de la dérivée fractionnaire $D^\alpha \epsilon_{12}^e(t_M)$ au premier point de la boucle. En d'autres termes, pour modéliser la boucle d'hystérésis qui correspond à l'intervalle temporel de « décharge-recharge » $[t_M, t_N]$, nous utilisons les valeurs de la dérivée fractionnaire sur l'intervalle « charge-décharge-recharge » $[t_L, t_N]$ afin d'éviter des problèmes d'ordres purement mathématiques.

Afin de calculer les dérivées fractionnaires, différentes approximations numériques peuvent être utilisées.
<caption>**Tableau 2.** Identification longitudinale (transversale) des paramètres élastiques du matériau composite carbone/époxy.</caption>

<table><caption>Table 2. Longitudinal (transverse) identification of elastic parameters of the carbon/epoxy composite.</caption><thead><tr><td>paramètre</td><td>Valeur</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>ν<sup>0</sup><sub>12</sub> = ν<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>

Comme nous l'avons décrit précédemment, l'erreur minimale du calcul numérique des dérivées fractionnaires correspond à l'approximation *M1* (8), (10) qui peut être explicitée dans notre cas de la façon suivante :

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

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

Enfin, afin d'utiliser les formules (33) et (34), il est nécessaire de reconstruire l'expression analytique de la fonction $\epsilon_{12}^e(t)$ par des fonctions décrites par morceaux.

# 4 Identification des paramètres

La procédure d'identification des paramètres est illustrée sur l'exemple du matériau composite tissé carbone/époxy à partir des données expérimentales présentées sur les Figures 3 et 4.

## 4.1 Caractérisation du matériau dans la direction des fibres

La caractérisation du matériau dans la direction longitudinale (transverse) est faite en utilisant la courbe expérimentale (Fig. 3). À partir de cet essai, les paramètres peuvent être identifiés facilement : en effet on a accès rapidement au module d'Young longitudinal $E_{11}^0$, au coefficient de Poisson $\nu_{12}^0$ ainsi qu'à la contrainte de rupture $\sigma_{11}^R$. Les paramètres identifiés de cet essai sont rappelés dans le Tableau 2.

Comme nous l'avons évoqué précédemment, le comportement du matériau composite dans le sens des fibres est de type élastique fragile. La rupture du matériau apparaît « instantanément ». Par conséquent, nous proposons d'utiliser une fonction de type Heaviside (23) afin de décrire ce comportement. Dans ce cas, le seuil initial d'endommagement coïncide avec le seuil de la rupture : $\bar{Y}_{11}^0 = \bar{Y}_{11}^R$. Les paramètres associés à l'endommagement sont présentés dans le Tableau 3. La propagation d'endommagement est illustrée sur la Figure 12.

## 4.2 Caractérisation du matériau en cisaillement

À partir des données expérimentales en cisaillement plan (Fig. 4), certains paramètres peuvent être identifiés

<caption>**Tableau 3.** Paramètres d'endommagement longitudinal du matériau composite carbone/époxy.</caption>

<table><caption>Table 3. Longitudinal damage parameters of the carbon/epoxy composite.</caption><thead><tr><td>paramètre</td><td>Valeur</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_13.png" /><figcaption>Fig. 12. Évolution de la fonction d'endommagement longitudinal pour le composite carbone/époxy.</figcaption></figure>

<caption>Fig. 12. Evolution of the longitudinal damage function for the carbon/epoxy composite.</caption>

<caption>**Tableau 4.** Identification des paramètres élastiques en cisaillement du composite carbone/époxy.</caption>

<table><caption>Table 4. Identification des paramètres élastiques en cisaillement du composite carbone/époxy.</caption><thead><tr><td>paramètre</td><td>Valeur</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>

directement. Le module de cisaillement non-endommagé $G_{12}^0$ est obtenu à partir d'une régression linéaire sur le début de la courbe. La contrainte d'écoulement initiale $R_0$ est le seuil de plasticité. La contrainte de rupture $\sigma_{12}^R$ correspond à la contrainte maximale pendant l'essai. Ces paramètres sont regroupés dans le Tableau 4.

Les courbes « contrainte-déformation » en cisaillement plan Figure 4 montrent généralement une diminution de la rigidité du matériau au cours des cycles. Ceci caractérise un accroissement de l'endommagement du module de cisaillement. Afin de décrire les dégradations supérieures par le matériau, il est nécessaire de déterminer la loi qui relie la variable d'endommagement $d_{12}$ à la force thermodynamique associée $\sqrt{\bar{Y}_{12}}$ en utilisant les équations (19) et (20). La méthode consiste à identifier pour chaque cycle, le module associé ce qui permet en connaissant le module initial de pouvoir remonter à l'endommagement. En cisaillement plan, l'endommagement est donné par :

$$ d_{12} = 1 - \frac{G_{12}^i}{G_{12}^0} \quad (35) $$
<caption>**Tableau 5.** Paramètres d'endommagement en cisaillement du composite carbone/époxy.</caption>

<table><caption>Table 5. Shear damage parameters for the carbon/composite composites.</caption><thead><tr><td>paramètre</td><td>Valeur</td></tr></thead><tbody><tr><td>√Y<sup>C</sup><sub>12</sub></td><td>2,25 MPa<sup>1/2</sup></td></tr><tr><td>√Y<sup>0</sup><sub>12</sub></td><td>0,15 MPa<sup>1/2</sup></td></tr><tr><td>√Y<sup>R</sup><sub>12</sub></td><td>0,8 MPa<sup>1/2</sup></td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Évolution de la fonction d'endommagement en cisaillement pour le composite carbone/époxy.</figcaption></figure>

<caption>Fig. 13. Damage function evolution in shear direction for the carbon/epoxy composite.</caption>

Dans le cas général, pour chacun des cycles, on peut également relever les déformations élastiques correspondantes ce qui permet d'identifier la force thermodynamique associée. Une fois cela fait, il est nécessaire d'ajuster le type de loi d'évolution (21), (22) ou (23) en fonction des observations expérimentales. L'identification de la contrainte initiale d'écoulement permet d'accéder aussi à la déformation associée. Cela permet d'identifier le seuil initial d'endommagement ($\bar{Y}_{12}^0$). Le taux maximum de restitution d'énergie ($\bar{Y}_{12}^R$) est quant à lui identifié sur le dernier cycle, proche de la rupture du matériau.

Après avoir appliqué la méthode proposée, nous obtenons l'évolution proposée par la Figure 13. Les conclusions que l'on peut obtenir sur l'aspect de cette évolution, nous pousse à proposer une évolution sous une forme linéaire (21).

Dans ce cas particulier de loi, rappelons que le coefficient $\sqrt{\bar{Y}_{12}^C}$ représente la vitesse de la propagation de l'endommagement de cisaillement. Les paramètres de cette loi sont proposés dans le Tableau 5.

Après avoir déterminé l'évolution de l'endommagement en cisaillement plan, il est possible maintenant d'évaluer la déformation plastique cumulée $p$ (28) sur la base des données identifiées. L'identification de la fonction d'écoulement plastique $R(p)$ est effectuée à l'aide de la formule suivante :

$$R(p) = \frac{|\sigma_{12}|}{(1 - d_{12})} - R_0 \quad (36)$$

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Fonction d'écrouissage $R(p)$ pour le composite carbone/époxy.</figcaption></figure>

<caption>Fig. 14. Hardening function $R(p)$ for the composite carbon/epoxy.</caption>

<caption>**Tableau 6.** paramètres de plasticité du composite carbone/époxy.</caption>

<table><caption>Table 6. Paramètres de plasticité du composite carbone/époxy.</caption><thead><tr><td>paramètre</td><td>Valeur</td></tr></thead><tbody><tr><td>β</td><td>266,6 MPa</td></tr><tr><td>m</td><td>0,36</td></tr></tbody></table>

Une fois les deux quantités calculées, il est possible de choisir la meilleure approximation pour décrire l'évolution de la fonction d'écrouissage en fonction de la déformation plastique cumulée. Au vu de notre identification (Fig. 14), la meilleure approximation correspond à une loi puisance (27). Les paramètres associés sont regroupés dans le Tableau 6.

## 4.3 Identification des paramètres du modèle fractionnaire

Les paramètres du modèle fractionnaire (30) $AB$ et $\alpha$ sont déterminés par la résolution d'un problème d'optimisation non-linéaire sous contraintes. La contrainte principale est imposée sur les valeurs de l'ordre de la dérivation fractionnaire $\alpha$ qui doit être tel que $0 < \alpha < 1$. La fonction « objectif » s'exprime sous la forme d'une fonction d'erreur relative :

$$\delta = \frac{\sqrt{\sum_{i=1}^{N} (\epsilon_{12}^{\text{exp}}(t_i) - \epsilon_{12}^{\text{modèle}}(t_i))^2}}{N} \quad (37)$$

ou $\epsilon_{12}^{\text{exp}}$ est la déformation expérimentale calculée par la formule (38), $\epsilon_{12}^{\text{modèle}}$ est la déformation calculée par le modèle fractionnaire (30) et $N$ est le nombre de points dans l'intervalle de temps considéré.

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

Les trois paramètres du modèle fractionnaire $A, B$ et $\alpha$ sont associés à chaque boucle d'hystérésis. En effet,
<caption>**Tableau 7.** Coefficients de l'évolution des paramètres du modèle fractionnaire pour le composite carbone/époxy.</caption>

<table><caption>Table 7. Evolution of the coefficients of the fractional model parameters for the carbon/epoxy composite.</caption><thead><tr><th>paramètre</th><th>Valeur</th><th>paramètre</th><th>Valeur</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>α</sub></td><td>0,6944</td><td>b<sub>α</sub></td><td>-0,0035</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Évolution du paramètre A pour le composite carbone/époxy.</figcaption></figure>

<caption>Fig. 15. Parameter A evolution for the carbon/epoxy.</caption>

l'évolution de ces paramètres au cours du chargement est exprimée en fonction de l'endommagement qui demeure constant pendant la phase de « décharge-recharge ». On identifie ainsi les meilleures approximations qui décrivent les relations entre les paramètres fractionnaires et l'endommagement.

La méthodologie proposée permet d'obtenir les évolutions de chacun des paramètres fractionnaires en fonction de l'endommagement de cisaillement obtenu à chaque cycle. Les évolutions sont représentées sur les Figures 15-17. À partir de ces données, on remarque que les approximations les plus adaptées sont de type linéaire :

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

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

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

Les coefficients correspondant aux approximations (39)-(41) sont donnés dans le Tableau 7.

# 5 Résultats

Suite à l'identification des paramètres du modèle collaboratif, nous proposons de valider les comportements et l'identification par des simulations numériques des essais expérimentaux au point d'intégration. Les résultats pour les différentes orientations sont donnés sur les Figures 18 and 19 pour le tissu composite carbone à matrice thermodurcissable (époxy).

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Évolution du paramètre B pour le composite carbone/époxy.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Évolution du paramètre α pour le composite carbone/époxy.</figcaption></figure>

<caption>Fig. 17. Parameter α evolution for the carbon/epoxy.</caption>

Les résultats de la simulation sont en bonne adéquation avec les données expérimentales pour le matériau composite carbone/époxy. De petites erreurs entre les données peuvent être observées au début et à la fin de la boucle d'hystérésis ainsi qu'au point de changement de signe de la déformation (transmission de la décharge en recharge). L'écart au début de la boucle est lié au problème mathématique du calcul des dérivées fractionnaires [41]. Le décalage entre les données et le modèle à la fin de la boucle provient du fait que l'optimisation ne peut pas être déterminée précisément sur les derniers points. La différence entre la courbe expérimentale et la simulation au bas de la boucle associée au changement du signe de la déformation provoque une erreur dans le calcul des dérivées fractionnaires ou entières pour les fonctions définies par morceaux. L'identification de l'endommagement peut être inexacte et c'est la raison supplémentaire qui provoque l'erreur sur le calcul numérique.

L'identification des paramètres a été également effectuée pour un matériau composite à matrice thermoplastique carbone/PA66. Dans ce cas la propagation de
<figure><img src="image_19.png" /><figcaption>Fig. 18. Comparaison expérimental/numérique sur le composite carbone/époxy pour une traction sens fibres.</figcaption></figure>

Fig. 18. Experimental/numerical behavior comparison in fibers direction for the carbon/epoxy composite.

<figure><img src="image_20.png" /><figcaption>Fig. 19. Comparaison expérimental/numérique sur le composite carbone/époxy en cisaillement.</figcaption></figure>

Fig. 19. Experimental/numerical behavior comparison in shear for the carbon/epoxy composite.

l'endommagement a un caractère beaucoup plus non-linéaire. Une approximation polynomiale de troisième degré a été utilisée. En effet, on remarque qu'au fur et à mesure des cycles, les boucles d'hystérésis commencent à perdre leur forme symétrique et elles prennent une forme de « haricot » (Fig. 20). Du fait de la nature de la matrice, le matériau dissipe plus d'énergie : on a effectivement une largeur de boucle beaucoup plus importante que dans le cas de la matrice thermodurcissable. Par ailleurs, cette nature de matrice conduit également à une réorientation des fibres plus forte que dans le cas du composite à matrice thermodurcissable lors de la phase de charge. En conséquence, pour ce type de matériau nous commettons une erreur lors de la mesure de l'endommagement, notamment sur les derniers cycles et cela se répercute sur la loi d'évolution obtenue. Il faudrait prendre en compte la réorientation des fibres au cours de la sollicitation afin d'obtenir une meilleure description de la propagation de

<figure><img src="image_21.png" /><figcaption>Fig. 20. Comparaison expérimental/numérique sur un composite carbone/PA66 en cisaillement.</figcaption></figure>

Fig. 20. Experimental/numerical comparison in shear for the carbon/PA66 composite.

l'endommagement. Ce point constitue une perspective importante pour les développements futurs.

En ce qui concerne l'évolution de la fonction d'écrouissage, la meilleure approximation en regard de l'exploitation des résultats correspond à une loi de type puissance tout comme pour le composite carbone/époxy.

De la même manière que précédemment, on post-traite les résultats en ce qui concerne les boucles d'hystérésis du composite à matrice thermoplastique. Tout comme pour l'endommagement, les approximations qui relient les paramètres du modèle fractionnaire à l'endommagement de cisaillement ne sont pas identiques à celles du composite à matrice thermodurcissable, de manière générale. Cela peut être imputé à la nature différente de la matrice. Les meilleures approximations des paramètres du modèle fractionnaire *A* et *B* sont donc des fonctions polynomiales de degré 3. L'ordre de dérivation fractionnaire *α* est approximé par une fonction linéaire.

Les résultats de la simulation des boucles d'hystérésis en termes de « déformation-contrainte » sont présentés sur la Figure 20 pour le composite à matrice thermoplastique. Une erreur plus importante est observée pour les résultats de la simulation numérique du modèle collaboratif pour le matériau composite T700/PA66. L'identification non-exacte des paramètres (l'endommagement et la plasticité) provoque un décalage supplémentaire entre les données en plus des erreurs associées aux problèmes numériques et mathématiques qui sont notées pour le matériau composite carbone/époxy. On peut observer néanmoins qu'avec le modèle fractionnaire, l'aire des boucles d'hystérésis semble être quant à elle identifiée correctement malgré les erreurs successives. C'est un point important notamment pour les objectifs de futurs travaux.

Afin de montrer les potentiels de la loi collaborative développée, nous proposons une dernière simulation numérique : nous allons imposer un chargement fictif qui possède 8 cycles (Fig. 21) et non 6 comme lors des essais réels expérimentaux. Cela permettra de montrer la capacité du modèle fractionnaire à prédire l'apparence
<figure><img src="image_22.png" /><figcaption>Fig. 21. Chargement cyclique fictif pour le composite carbone/PA66.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 22. Simulation avec un chargement fictif pour un essai en cisaillement pour le composite carbone/PA66.</figcaption></figure>

<caption>Fig. 22. Simulation of a fictional loading in shear for the carbon/PA66 composite.</caption>

des boucles d'hystérésis additionnelles. Le matériau que nous avons retenu est un composite à matrice thermoplastique et nous utiliserons les paramètres identifiés. Les résultats de la simulation numérique sont présentés sur la Figure 22. On peut remarquer que la position et la taille des boucles d'hystérésis supplémentaires sont tout à fait correctes en regard de la courbe expérimentale. Cela justifie les choix que nous avons faits tant d'un point de vue modèle que méthodologie d'identification expérimentale des paramètres. La modélisation du comportement élastoplastique endommageable visqueux permet d'exprimer la quantité d'énergie dissipée par les boucles d'hystérésis dans le cas d'essais de fatigue ou d'auto-échauffement.

## 6 Conclusion

Dans ce travail, nous avons proposé un modèle unique pour des matériaux composites tissés en faisant collaborer deux modèles : le premier un modèle élasto-plastique

endommageable et le second un modèle fractionnaire. La capacité du modèle élasto-plastique endommageable de reproduire correctement le comportement du matériau composite a été validée de nombreuses fois. Le modèle fractionnaire est représenté par un schéma de Kelvin-Voigt non-standard avec peu de paramètres qui utilise la dérivée fractionnaire au sens de Riemann-Liouville. Ce modèle est utilisé pour reproduire des boucles d'hystérésis et appliqué seulement pendant la « décharge-recharge », ce qui permet de diminuer le temps du calcul. En plus, le modèle fractionnaire est complémentaire du modèle élasto-plastique endommageable car il peut être utilisé de manière optionnelle.

La validation du modèle a été effectuée pour deux matériaux composites à matrice différente : thermodiscassable et thermoplastique. On constate que le modèle proposé est capable de représenter le comportement élasto-plastique endommageable aussi que les boucles d'hystérésis. Les résultats sont en bonne concordance avec les données expérimentales. Cependant, on peut remarquer que le modèle fractionnaire n'est pas capable de représenter les boucles non-symétriques qui sont observées dans la réponse du matériau composite à matrice thermoplastique. En effet, il est nécessaire de poursuivre les travaux, notamment en ce qui concerne la prise en compte de la variation des angles des fibres en cours de chargement. Un autre point qui est actuellement investigué concerne la capacité du modèle collaboratif à reproduire le comportement des composites tissés pour des vitesses de déformation allant jusqu'à une gamme moyenne. Les premiers résultats obtenus semblent forte-ment encourageants.

## Références

[1] O. Castro, M. Lennie, K. Branner, G. Pechlivanoglou, P. Brondsted, C. Nayeri, C. Paschereit, Comparing fatigue life estimations of composite wind turbine blades under various load cases and fatigue methods, 20th International Conference on Composite Materials (ICCM20), Cophenague, 2015

[2] D. Caous, B. Bois, J. Wahl, T. Palin-Luc, V.J., An approach to identify multiaxial fatigue life model dedicated to wind turbine blades, 20th International Conference on Composite Materials (ICCM20), Cophenague, 2015

[3] A. Elias, M. Kaminski, F. Laurin, L. Gornet, Study of low velocity impact defects in organic interlock woven composite, 16th European Conference on Composite Materials (ECCM16), Seville, 2014

[4] L. Gornet, J. Guldner, P. Rozycki, P. Cartraud, L. Marcin, Numerical modelling of strain rate effects on a woven composite fabric, 16th European Conference on Composite Materials (ECCM16), Seville, 2014

[5] P. Ladeveze, E. Le Dantec, *Compos. Sci. Technol.* **43** (1992) 257-267

[6] P. Ladeveze, O. Allix, J. Deu, D. Leveque, *Comput. Methods Appl. Mecha. Eng.* **183** (2000) 150-122

[7] C. Hochard, Y. Thollon, *Int. J. Fatigue* **32** (2010) 1158-165
[8] P. Rozycki, Contribution au développement de lois de comportement pour matériaux composites soumis à l'impact, These de doctorat, Valenciennes, 2000

[9] S. Marguet, P. Rozycki, L. Gornet, *Mech. Adv. Mater. Struct.* **14** (2007) 619-631

[10] G.W. Leibniz, G. L'Hopital, Leibniz ensgesammelte Werke, Lebinizens mathematische Schriften, Erste Abtheilung, Band II, Pertz G. H. et Gerhardt C.J., 1849, pp. 301-302

[11] N.H. Abel, *Ouvres complètes* **1** (1881) 11-27

[12] H. Holmgren, *Kongl. Svenska Vetenskaps Akad. Hand* **5** (1865-1866) 1-83

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

[14] R. Bagley, P. Torvik, *Rheology* **27** (1983) 201-210

[15] P. Cosson, J. Michon, *Chaos, Solitons Fractals* **7** (1996) 1807-1824

[16] M. Sjöberg, On dynamic properties of rubber isolators, Stockholm, 2002

[17] A. Lion, C. Kardelky, *Int. J. Plasticity* **20** (2004) 1313-1345

[18] Y. Rabotnov, Elementy Naledstvennoj Mehaniki Tverdyh Tel, Nauka, Moscau, 1977

[19] V.V. Uchaikin, Metod drobnih proizvodnih, Artishok Ulianovsk, 2008

[20] M.H. Mateos, Hysteretic and damage behaviour modeling of composite materials by fractional operators, These de Doctorat, École Centrale de Nantes, 2014

[21] M. Caputo, J. Carcione, *Rheologica Acta* **50** (2011) 107-115

[22] P. Rozycki, Loi de comportement pour composite et outils numériques en dynamique rapide, Habilitation à Diriger des Recherches, Nantes, 2015

[23] O. Westphal, 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, These de Doctorat, École Centrale de Nantes, 2014

[24] S. Giancane, A. Chrysochoos, V. Dattoma, B. Wattrisse, *Theoret. Appl. Fract. Mech.* **2** (2009) 117-121

[25] V. Le Saux, Y. Marcoa, S. Callocha, C. Doudarda, P. Charrierb, *Proc. Eng.* **2** (2010) 949-958

[26] L. Jegou, Y. Marco, V. Le Saux, S. Calloch, *Int. J. Fatigue* **47** (2013) 259-267

[27] L. Gornet, O. Wesphal, C. Burtin, J.-C. Bailleul, P. Rozycki, L. Stainier, *Proc. Eng.* **66** (2013) 697-704

[28] L. Gornet, O. Westphal, M. Mateos, A. Krasnobrizha, P. Rozycki, C. Peyrac, F. Lefebvre, *Rev. Compos. Mtériaux Avancés* **25** (2015) 181-200

[29] J. Liouville, *Journal de l'École Royale Polytechnique* **13** (1832) 1-69

[30] G. Riemann, Versuch Einer Allgemeinen Auffasung Der Integration und Differentiation, 1876 M. Riesz, Mathematiker in Göttingen, Leipzig, Teubner, 1876, pp. 331-344

[31] A.K. Grünwald, *Zeits. Math. Phys* **12** (1867) 441-480

[32] A. Letnikov, *Matematicheskij Sbornik* **3** (1868) 1-68

[33] M. Caputo, *Annales di Geofisica* **19** (1966) 383-393

[34] P. Ladeveze, Sur la mécanique de l'endommagement des composites, Journées Nationales sur les Composites JNC5, Paris 1986

[35] O. Allix, D. Gilletta, P. Ladeveze, Mechanical behaviour of elementary constituents of laminates, Fifth International Conference on Composite Materials, ICCM-V, San Diego, 1985

[36] J. Lemaitre, J. Chaboche, A. Benallal, R. Desmorat, Mécanique des Matériaux Solides, Dunod, 2009

[37] A. Johnson, A. Pickett, P. Rozycki, *Compos. Sci. Technol* **61** (2001) 2183-2192

[38] P. Ladeveze, Inelastic strains and damage, Vol. 9 of Damage Mechanics of Composite Materials, edited by R. Talreja, Amstardam, Elsevier, 1994, pp. 117-138

[39] C.T. Herakovich, Mechanics of Fibrous Composites, Wiley, New-York, 1998

[40] C. Hochard, J. Payan, O. Montagnier, *Compos. Sci. Technol.* **65** (2005) 467-474

[41] 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, These de doctorat, École Centrale de Nantes, 2015