# Modélisation, identification et simulations éléments finis des phénomènes de délaminage dans les structures composites stratifiées

Laurent Gornet, David Lévêque, Lionel Perret

► To cite this version:

Laurent Gornet, David Lévêque, Lionel Perret. Modélisation, identification et simulations éléments finis des phénomènes de délaminage dans les structures composites stratifiées. Mechanics & Industry, 2000, Numéro Spécial Matériaux Composites, 1 (3), pp.267-276. hal-01007077

HAL Id: hal-01007077

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

Submitted on 21 Jan 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 CC BY 4.0 - Attribution - International License
# Modélisation, identification et simulations éléments finis des phénomènes de délaminage dans les structures composites stratifiées

Laurent Gornet<sup>a</sup>, David Lévêque<sup>b</sup>, Lionel Perret<sup>c</sup>

<sup>a</sup> ECN, École centrale de Nantes, LMM/DMS, Laboratoire mécanique et matériaux, Division mécanique des structures, 1 rue de la Noë, BP 92101, Nantes cedex 03, France

<sup>b</sup> Onera, Office national d'études et de recherches aérospatiales, Département matériaux et systèmes composites, BP 72, 29 avenue de la division Leclerc, 92322 Châtillon cedex, France

<sup>c</sup> CNES, Centre spatial de Toulouse, 18, av. Édouard-Belin, 31401 Toulouse cedex 4, France

**Résumé** — Pour prévoir la simulation de l'initiation et de la propagation du délaminage par la méthode des éléments finis, il a été proposé dans des études antérieures de décrire un stratifié comme un assemblage de couches et d'interfaces susceptibles de se dégrader. Le présent travail concerne la modélisation, l'identification et la simulation des phénomènes de délaminage dans les structures composites stratifiées. Les simulations sont effectuées en tenant compte des non-linéarités géométriques et matérielles. Les modèles d'interface développés sont implantés dans le code éléments finis Castem 2000-CEA. Des essais classiques d'initiation et de propagation du délaminage sont analysés dans cette étude afin de montrer la pertinence de l'approche retenue. © 2000 Éditions scientifiques et médicales Elsevier SAS

matériaux composites / endommagement / délaminage / éléments finis / rupture / non linéaire

**Abstract** — Modelling, identification and finite element predictions of delamination in laminated composite structures. In order to forecast the delamination initiation and propagation in a finite element context, a previously defined damage meso-modelling of composite laminates is used. At the meso-level, the laminate is described as an assembly of damageable layers and interlaminar interfaces. The present work concerns the modelling, the identification and the finite element predictions of delamination phenomena in composite structures. Finite element predictions are conducted with non-linear geometric and material hypothesis. The interface modelling is implemented in the finite element code Castem 2000 developed by CEA. Classic edge delamination tension and propagation tests are conducted in order to improve the interface damage model approach in the finite elements context. © 2000 Éditions scientifiques et médicales Elsevier SAS

composite material / damage / delamination / finite elements / fracture / nonlinear

## 1. INTRODUCTION

Cette étude vise à maîtriser la modélisation, l'identification et la simulation éléments finis de l'initiation et de la propagation des endommagements et des délaminages dans les structures composites stratifiées. Pour étudier la propagation d'un délaminage établit dans les matériaux composites stratifiés, des essais du type mécanique de la rupture sont classiquement utilisés. L'éprouvette d'essai considérée peut alors être assimilée à une poutre préfissu-

rée. La fissure est obtenue par un film non adhérent inséré lors de la fabrication des éprouvettes. Lors des essais, l'avancée de la fissure est suivie en fonction de la charge appliquée, ce qui permet de déterminer un taux de restitution d'énergie critique correspondant à l'énergie nécessaire pour faire croître le délaminage d'une surface unitaire. Suivant la géométrie de la poutre et le type de chargement considéré, différents modes de sollicitation de l'interface prédélaminée sont obtenus. De nombreux traumas, basés sur la mécanique de l'endommagement, ont été réalisés depuis plus de quinze ans dans le domaine de la modélisation des matériaux composites [1–6]. Le travail présenté est le prolongement de l'étude sur les bords
droits [7] au cadre plus général tridimensionnel [3, 8]. Ces travaux sont à la base de la réalisation et de l'implantation dans le code de calcul éléments finis Castem 2000-CEA [8, 9] de modèles dédiés à l'analyse des structures composites stratifiées. L'utilisation ici de modèle d'évolution de l'endommagement avec effet de retard assure des résultats objectifs vis à vis du maillage dans le cas de fronts de délaminage complexes.

L'objet de cet article est de présenter la modélisation, l'identification et la simulation éléments finis d'esais d'initiation et de propagation de délaminages réalisés à l'Aérospatiale-Suresnes [10]. Les éprouvettes d'esai peuvent être considérées comme de petites structures. Elles nécessitent donc une analyse éléments finis avec prise en compte des différentes non linéarités (géométriques et matérielles). Les simulations éléments finis montrent la pertinence de l'approche retenue.

## 2. MÉSO-MODÉLISATION DES STRUCTURES STRATIFIÉES

La méso-modélisation des structures stratifiées permet de décrire à tout instant l'état de dégradation de la structure jusqu'à sa rupture ultime. La structure stratifiée est modélisée à une échelle intermédiaire dite « méso » située entre l'échelle de la structure et celle des microconstituants. À cette échelle, la structure est décrite au moyen de constituants de base appelés méso-constituants (couches et interfaces). L'interface inter-laminaire est introduite pour modéliser les phénomènes de délaminage qui apparaissent au voisinage des bords ou des macrodéfauts. Une fois déterminé le comportement mécanique de la couche et de l'interface, le comportement de n'importe quel stratifié est bâti à partir de ces constituants. Le comportement loin des bords du stratifié est obtenu à l'aide de la théorie classique des stratifiés. Au voisinage des bords, il est nécessaire de procéder à une analyse tridimensionnelle. L'intérêt de cette modélisation est de décrire de façon simple les mécanismes complexes de dégradation. Les problèmes liés à l'utilisation jusqu'à la rupture ultime de modèles d'endommagement locaux classiques sont jugulés par la méso-modélisation qui impose un état d'endommagement localement constant dans l'épaisseur de la couche [6, 20]. Des modèles d'évolution de l'endommagement avec effet de retard permettent de définir une zone de rupture dans le plan de chacune des couches et des interfaces [6, 8].

<figure><img src="image_2.png" /><figcaption>Figure 1. Interface : entité surfacique entre deux couches adjacentes.</figcaption></figure>

## 3. MODÉLISATION DE L'INTERFACE

L'interface est une entité surfacique qui assure le transfert des déplacements et des contraintes normales entre deux plis adjacents. Cette modélisation couplée à la mécanique de l'endommagement permet de prendre en compte les phénomènes de délaminage (amorçage et propagation) qui peuvent survenir lors de la mise en charge de pièces mécaniques. Des études précédentes ont permis de montrer la faisabilité de la prédiction du délaminage en utilisant cette modélisation. Les zones à forts gradients sont le lieu de naissance privilégié du délaminage. Dans ces zones, le champ des contraintes est tridimensionnel.

Le comportement mécanique de l'interface ne dépend que de l'angle relatif entre les fibres des plis adjacents. Un repère « d'orthotropie », de directions principales $N_1$, $N_2$ et $N_3$ ($N_3$ normale à la plaque), est alors défini par les bissectrices des directions de fibres (figure 1). La base $N_1$, $N_2$, $N_3$ est nécessairement orthotrope puisqu'une interface $[\theta_1/\theta_2]$ est équivalente à une interface $[\theta_2/\theta_1]$. On note le saut de déplacement dans la base d'orthotropie de l'interface :

$$[U] = U^{\text{sup}} - U^{\text{inf}} = [U_1]N_1 + [U_2]N_2 + [U_3]N_3$$

L'énergie de déformation du matériau endommagé (potentiel thermodynamique de Gibbs) s'écrit :

$$E_D = \frac{1}{2} \left[ \frac{\langle -\sigma_{33} \rangle_+^2}{k_3^0} + \frac{\langle \sigma_{33} \rangle_+^2}{k_3^0 (1 - d_3)} + \frac{\sigma_{32}^2}{k_2^0 (1 - d_2)} + \frac{\sigma_{31}^2}{k_1^0 (1 - d_1)} \right]$$

Où $\langle x \rangle_+ = (|x| + x)/2$ représente la partie positive de la quantité $x$. La détérioration de l'interface est prise en compte par trois variables internes d'endommagement. La relation de comportement entre les contraintes et les sauts de déplacement s'écrit dans la base d'orthotropie de l'interface à partir du potentiel thermodynamique :
<figure><img src="image_3.png" /><figcaption>Figure 2. Modes de propagation du délaminage.</figcaption></figure>

$$
\begin{pmatrix} \sigma_{13} \\ \sigma_{23} \\ \sigma_{33} \end{pmatrix} = \begin{pmatrix} k_1^0 (1-d_1) & 0 & 0 \\ 0 & k_2^0 (1-d_2) & 0 \\ 0 & 0 & k_3^0 (1-d_3) \end{pmatrix} \begin{pmatrix} [U_1] \\ [U_2] \\ [U_3] \end{pmatrix}
$$
avec
$$
\begin{cases}
k_1^0 = \frac{2G_{13}}{e} \\
k_2^0 = \frac{2G_{23}}{e} \\
k_3^0 = \frac{E_3}{e}
\end{cases}
$$

Où les rigidités d'interface $k_1^0$, $k_2^0$ et $k_3^0$ sont liées aux
caractéristiques élastiques du stratifié considéré et $e$ est
représentatif de l'épaisseur de l'interface. Le modèle
thermodynamique est construit en tenant compte des trois
modes de délaminage possibles. On distinguera bien dans
ce qui suit, les différentes variables d'endommagement
de l'interface inter-laminaire :

• $d_3$ est associée au mode « d'ouverture » (ou mode 1) de la liaison inter-laminaire ;

• $d_1$ et $d_2$ sont associées aux deux modes de cisaillement de l'interface (modes 2 et 3). Il faut noter qu'en compression, l'interface ne s'endommage pas.

Les forces thermodynamiques conjuguées aux va-
riables d'endommagement et associées aux trois modes
de délaminage sont :

$$
Y_{d3} = \frac{1}{2} \frac{\langle \sigma_{33} \rangle_{+}^2}{k_{3}^{0} (1 - d_{3})^{2}}, \quad Y_{d1} = \frac{1}{2} \frac{\sigma_{31}^{2}}{k_{1}^{0} (1 - d_{1})^{2}}
$$

$$
Y_{d2} = \frac{1}{2} \frac{\sigma_{32}^2}{k_2^0 (1 - d_2)^2}
$$

L'énergie dissipée dans ce modèle est :

$$
\Phi = Y_{d1} \dot{d}_1 + Y_{d2} \dot{d}_2 + Y_{d3} \dot{d}_3 \quad (\Phi \ge 0)
$$

3.1. Modèle d'évolution
de l'endommagement
sans effet de retard

Dans ce qui suit une loi d'évolution de l'endommage-
ment simple, dite « isotrope », est introduite. Il est sup-
posé en effet que les différents endommagements sont
très fortement couplés et gouvernés par un taux de res-
titution d'énergie équivalent unique qui s'écrit :

$$
\underline{Y}(t) = \mathrm{Sup}_{|\tau| \le t} \left[ (Y_{d_3})^{\alpha} + (\gamma_1 Y_{d_1})^{\alpha} + (\gamma_2 Y_{d_2})^{\alpha} \right]^{1/\alpha}
$$

où γ₁ et γ₂ sont des paramètres de couplage et α un autre
paramètre matériau permettant de décrire le lieu de rup-
ture en Mécanique de la Rupture. Par rapport au modèle
d'interface présenté dans [4], le modèle présenté ici, ap-
paremment plus complexe, a pour avantage de conduire à
une identification directe et à une interprétation simple
des différents coefficients. L'interface étant une zone
d'épaisseur très faible, il semble assez naturel de consi-
dérer que l'interface rompue dans un mode de sollicita-
tion donné l'est pour tous les autres modes. En l'absence
d'informations plus précises, il semble raisonnable d'uti-
liser le modèle d'endommagement le plus simple et com-
patible avec les informations expérimentales disponibles.
La loi d'évolution de l'endommagement est alors définie
par le choix d'une fonction matériau ω(Y) telle que :

$$
\begin{gather*}
\text{si } [(d_3 < 1) \text{ et } (\underline{Y} < Y_R)] \\
d_1 = d_2 = d_3 = \omega(\underline{Y})
\end{gather*}
$$

$$
\begin{gather*}
\text{sinon} \\
d_1 = d_2 = d_3 = 1
\end{gather*}
$$

La fonction de délaminage est choisie sous la forme
suivante :

$$
\omega(\underline{Y}) = \left[ \frac{n}{n+1} \frac{\langle \underline{Y} - Y_0 \rangle_+}{Y_c - Y_0} \right]^n
$$

où donc tous les endommagements sont supposés iden-
tiques. Certaines des hypothèses précédentes concernant
le modèle d'endommagement pourraient être abandon-
nées si ce modèle s'avère insuffisant. En particulier, on
pourrait supposer que l'évolution de d₃ soit différente de
celle de d₁ et d₂. Seules des informations expérimentales
concernant l'initiation pourront éventuellement nous ren-
seigner sur ce point. En résumé, le modèle d'interface,
hors les caractéristiques de raideurs initiales, dépend de
sept paramètres matériau :

• $Y_0$ : seuil d'endommagement,

• $Y_c$ : énergie critique d'endommagement,
• $Y_R$ : force thermodynamique à rupture associée à l'endommagement ultime $d_R$ (permet de simuler une rupture avec $d_R \le 1$)

$$Y_R = Y_0 + \frac{n+1}{n} d_c^{1/n} (Y_c - Y_0)$$

• $\gamma_1$ et $\gamma_2$ : paramètres de couplage énergies de cisaillement/énergie transverse,

• $n$ : caractérise la plus ou moins grande « fragilité » de l'interface (plus $n$ est grand, plus l'interface est fragile),

• et $\alpha$. L'interprétation de ce paramètre est moins aisée, $\alpha$ permet cependant de retrouver une forme classique du lieu de rupture en propagation.

À titre d'exemple, l'évolution de la réponse pour deux valeurs de la rigidité d'interface est présentée pour une interface [0] du matériau M55J/M18 (fibres de carbone très haut module/résine époxy).

## 3.2. Modèle d'évolution de l'endommagement avec effet de retard

Les difficultés de la caractérisation de la rupture et le risque de dépendance au maillage des résultats peuvent être résolus dans le cas du modèle d'interface en introduisant une modélisation avec effet de retard [6]. La signification physique de ce type de modèle est que l'évolution de l'endommagement n'est pas instantanée par rapport au chargement. Ce type de modèle est utilisé pour les simulations de fronts de délaminage tridimensionnels complexes [8]. Une variante est proposée après. L'évolution du délaminage est supposée gouvernée par la force $Y(t)$ :

$$Y(t) = ((Y_{d_3})^{\alpha} + (\gamma_1 Y_{d_1})^{\alpha} + (\gamma_2 Y_{d_2})^{\alpha})^{1/\alpha}$$

Les lois d'évolution de l'endommagement avec effet de retard s'écrivent :

$$\begin{gathered} \text{si } [(d_3 < 1) \text{ et } (Y < Y_R)] \\ \dot{d}_3 = k(\omega(Y) - d_3)_+^m, \quad \omega(Y) \le 1 \\ \dot{d}_1 = \dot{d}_2 = \dot{d}_3 \end{gathered}$$

sinon

$$d_3 = 1$$

La procédure utilisée pour identifier les paramètres des modèles avec et sans effet de retard est présentée dans la suite.

<figure><img src="image_4.png" /><figcaption>Figure 3. Déplacement imposé sur la face supérieur de l'interface et bloqué sur la face inférieur.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Évolution du déplacement imposé sur la face supérieur de l'interface en fonction du temps.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Évolution de la contrainte en fonction du saut de déplacement. Les valeurs des paramètres utilisés dans la simulation éléments finis sont : $Y_0 = 0 \text{ N}\cdot\text{mm}^{-1}$, $Y_c = 0,11 \text{ N}\cdot\text{mm}^{-1}$, $Y_R = Y_0 + (Y_c - Y_0)(n + 1)/n$, $\gamma_1 = 0,4$, $n = 0,5$, $\alpha = 1,6$.</figcaption></figure>
## 4. LIEN ENTRE LA MÉCANIQUE DE L'ENDOMMAGEMENT ET LA MÉCANIQUE DE LA RUPTURE

Une première idée pour déterminer les paramètres d'énergie critiques des modèles avec et sans effet de retard est d'utiliser les essais classiques de la Mécanique Élastique Linéaire de la Rupture qui permettent d'obtenir des taux de restitution d'énergie critiques en modes purs et mixtes. Une façon simple de relier la Mécanique de l'Endommagement à la Mécanique Élastique Linéaire de la Rupture est de comparer la dissipation mécanique donnée par les deux approches. Cette approche est utilisée dans [17] pour des modèles indépendants du temps. Les principaux résultats sont présentés ci-après. Pour des sollicitations en modes purs, les relations sont :

$$G_{cI}^{P} = Y_{c}, \quad G_{cII}^{P} = \frac{Y_{c}}{\gamma_{1}}, \quad G_{cIII}^{P} = \frac{Y_{c}}{\gamma_{2}}$$

Et dans un cas de chargement mixte, où $G_c^P = G_I + G_{II} + G_{III}$ représente le taux de restitution d'énergie critique global à propagation, on obtient :

$$\left(\frac{G_I}{G_{cI}^P}\right)^\alpha + \left(\frac{G_{II}}{G_{cII}^P}\right)^\alpha + \left(\frac{G_{III}}{G_{cIII}^P}\right)^\alpha = 1$$

Cette équation gouverne la forme du lieu de rupture en mode mixte déterminée par le paramètre $\alpha$. Les différentes relations montrent que les seuls paramètres influents de la loi d'évolution de l'endommagement en propagation sont $Y_c, \gamma_1, \gamma_2$ et $\alpha$.

La variante du modèle dit avec effet de retard introduit deux paramètres supplémentaires $k, m$. Il faut noter que lorsque les vitesses sont faibles, le modèle avec effet de retard coïncide avec le modèle statique pour un jeu de paramètres $k, m$ choisis. Dans la suite, nous choisissons $m=1$ afin de réduire le nombre de paramètre à identifier. Nous proposons d'identifier l'influence du temps caractéristique $k$ à partir de l'énergie de décohésion de l'interface. Il faut noter que $G_{cI}^P$ peut être interprété comme le travail fourni par l'unité de surface pour rompre l'interface en mode I. En mode I, l'énergie de décohésion est l'aire relative à la courbe présentée figure 7. Il suffit d'évaluer la différence entre les dissipations des modèles avec et sans effet de retard pour identifier l'influence du temps caractéristique $k$.

Dans la suite, on présente les méthodes numériques retenues afin d'implanter les modèles dans le cadre éléments finis.

<figure><img src="image_7.png" /><figcaption>Figure 6. Évolution de la contrainte en fonction du saut de déplacement pour différente valeur du temps caractéristique $k$. Les valeurs des paramètres utilisés dans la simulation éléments finis sont : $Y_0 = 0 \text{ N}\cdot\text{mm}^{-1}$, $Y_c = 0,11 \text{ N}\cdot\text{mm}^{-1}$, $Y_R = Y_0 + (Y_c - Y_0)(n + 1)/n$, $\gamma_1 = 0,4$, $n = 0,5$, $\alpha = 1,6$, $m = 1$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Évolution de la contrainte en fonction du saut de déplacement pour le modèle avec effet de retard. Les valeurs des paramètres utilisés dans les simulations éléments finis sont : $Y_0 = 0 \text{ N}\cdot\text{mm}^{-1}$, $Y_c = 0,11 \text{ N}\cdot\text{mm}^{-1}$, $Y_R = Y_0 + (Y_c - Y_0)(n + 1)/n$, $\gamma_1 = 0,4$, $n = 0,5$, $\alpha = 1,6$, $m = 1$.</figcaption></figure>

## 5. IMPLANTATION NUMÉRIQUE DU MODÈLE D'INTERFACE

Dans le cadre éléments finis, la description géométrique de l'interface est réalisée à l'aide d'un élément isoparamétrique tridimensionnel dégénéré sans épaisseur de type [11]. Cette modélisation permet d'associer une formulation éléments finis au modèle d'interface. L'algorithme utilisé dans le cadre éléments finis classiques pour résoudre des problèmes comportant des non linéarités géométriques et matérielles peut être scindé sous la forme suivante : la première est l'étape locale qui consiste
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Essais de propagation du délaminage utilisés classiquement en mécanique de la rupture. Ces essais sont destinés à identifier les paramètres des modèles d'interfices.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Taux de restitution d'énergie critiques moyens corrigés des phénomènes dissipés dans les couches.</figcaption></figure>

à intégrer les lois de comportement retenues pourodeliser le matériau. La seconde est l'étape globale où l'équilibre de la structure est évalué à partir des équations d'équilibre. Un résidu est construit à partir des équations d'équilibre et du champ de contrainte issu de l'étape locale. Ce résidu doit être vérifié avec une précision fixée par l'utilisateur pour arrêter le schéma itératif constitué des deux étapes de la méthode incrémentale. Un proces- sus d'accélération de convergence de type BFGS est uti- lisé pour s'affranchir du calcul de la matrice de rigidité tangente [12]. Cette méthodologie permet d'obtenir une vitesse de convergence satisfaisante et limite le coup de

calcul de la phase globale. Les algorithmes numériques retenus pour intégrer les modèles de comportement de l'interface sont maintenant présentés.

## 5.1. Algorithme d'intégration du modèle sans retard

L'équation discrétisée en temps qui représente l'évolution de la dégradation normale au plan de l'interface s'écrit sous la forme suivante : $d_{3n+1} = \omega(Y_{n+1})$. Celle-ci peut être écrite en fonction des sauts de déplacements. L'implantation numérique ne pose alors pas de difficultés particulières. L'algorithme d'intégration numérique est proposé ci-après :

1. Evaluation $\omega(Y_{n+1})$ avec :

$$ 
\begin{aligned}
\underline{Y}_{n+1}(t) &= \text{Sup}_{\tau \le t} \left[ \left( (Y_{d_3})^\alpha + (\gamma_1 Y_{d_1})^\alpha + (\gamma_2 Y_{d_2})^\alpha \right)^{1/\alpha} \right]_{n+1} \\
Y_{d_3} &= \frac{1}{2} k_3^0 [U_3]_+^2, \qquad Y_{d_2} = \frac{1}{2} k_2^0 [U_2]^2 \\
Y_{d_1} &= \frac{1}{2} k_1^0 [U_1]^2
\end{aligned}
$$

2. Calculs des endommagements :

$$ 
\begin{aligned}
d_{1n+1} &= \omega(Y_{n+1}), & d_{2n+1} &= \omega(Y_{n+1}) \\
d_{3n+1} &= \omega(Y_{n+1})
\end{aligned}
$$

3. Calculs des contraintes :

$$ 
\begin{aligned}
\sigma_{33n+1} &= k_3^0 (1 - d_{3n+1}) [U_3]_{n+1} \\
\sigma_{23n+1} &= k_2^0 (1 - d_{2n+1}) [U_2]_{n+1} \\
\sigma_{13n+1} &= k_1^0 (1 - d_{1n+1}) [U_1]_{n+1}
\end{aligned}
$$

4. Actualisations et retour à la phase d'équilibre.

## 5.2. Algorithme d'intégration du modèle avec retard

Un schéma numérique de type Euler implicite est utilisé pour intégrer le modèle de comportement élastique endommageable avec retard. L'équation discrétisée en temps qui représente l'évolution de la dégradation normale au plan de l'interface s'écrit sous la forme suivante :

$$ d_{3n+1} = d_{3n} + \Delta t k \langle \omega(Y_{n+1}) - d_{3n+1} \rangle_+^m $$

Celle-ci est non linéaire par rapport à l'endommagement, un schéma de Newton est utilisé pour la résoudre. L'algorithme d'intégration est proposé ci-après.
1. Evaluation $\omega(Y_{n+1})$.

2. Initialisation :

$$i = 0, \quad d_{3n+1}^i = d_{3n}$$

3. Evaluation du résidu :

$$R_{n+1}^i = d_{3n+1}^i - d_{3n} - \Delta t k \langle w(Y_{n+1}) - d_{3n+1}^i \rangle_+^m$$

4. Testes :

si $R_{n+1}^i < \text{Tol}$ alors $d_{3n+1} = d_{3n+1}^{i+1}$
aller en 5

sinon

$$\begin{align*}
\Delta d_{3n+1} &= -R_{n+1}^i \left( \frac{\partial R_{n+1}^i}{\partial d_{3n+1}} \right)^{-1} \\
d_{3n+1}^{i+1} &= d_{3n+1}^i + \Delta d_{3n+1} \\
i &= i + 1
\end{align*}$$

aller en 3

5. Calculs des endommagements $d_1, d_2$ :

$$d_{1n+1} = d_{1n} + (d_{3n+1} - d_{3n})$$

$$d_{2n+1} = d_{2n} + (d_{3n+1} - d_{3n})$$

6. Calculs des contraintes :

$$\sigma_{33n+1} = k_3^0 (1 - d_{3n+1}) [U_3]_{n+1}$$

$$\sigma_{23n+1} = k_2^0 (1 - d_{2n+1}) [U_2]_{n+1}$$

$$\sigma_{13n+1} = k_1^0 (1 - d_{1n+1}) [U_1]_{n+1}$$

7. Actualisations et retour à la phase d'équilibre.

## 6. SIMULATIONS ÉF DES ESSAIS DE PROPAGATION

Nous présentons ici les simulations d'essais de type mécanique de la rupture DCB (double cantilever beam), ENF (end notched flexure), MMF (mixed mode flexure), CLS (cracked lap shear) [13–16]. Ces essais sont nécessaires à l'identification et la validation des paramètres du modèle d'interface élastique endommageable. Les essais ont été réalisés à l'Aérospatiale-Suresnes [10] sur le matériau M55J/M18 pour différentes stratifications. Les essais de propagation de fissure sur composites se font sur des éprouvettes de type poutre dans lesquelles un film anti-adhésif est inséré, lors de la fabrication, à l'interface étudiée pour initier le délaminage. La géométrie

« moyenne » des éprouvettes utilisées ici est de 300 mm de long, 20 mm de large, le film anti-adhésif étant inséré à l'interface centrale sur une longueur de 40 mm environ. Le drappage des éprouvettes est de la forme [$(+\theta/-\theta)_{4s}/(-\theta/+\theta)_{4s}$], pour des angles $\theta$ de 0°, 22,5° et 45° correspondants aux trois types d'interface $\pm\theta$ envisagés. L'épaisseur moyenne d'un pli est de l'ordre de 0,1 mm. Les comparaisons proposées dans la suite entre les simulations et les résultats expérimentaux sont relatives à la campagne d'essai réalisée sur la stratification [0]32. Pour simuler le comportement de ce composite unidirectionnel, le comportement de la couche est choisi élastique orthotrope et l'interface élastique endommageable. Ce choix est guidé par l'observation d'un endommagement négligeable dans les couches pour ce type de stratification. Pour des stratifications croisées un état d'endommagement dans les couches désorientées est observé expérimentalement [22, 23], ce phénomène est alors pris en compte dans les simulations au travers d'un modèle élasto-plastique endommageable dans les couches. La géométrie est la même pour tous les essais, seul les conditions limites sont différentes. Les simulations sont effectuées à partir d'un maillage formé de deux couches constituées d'éléments tridimensionnels reliées par des éléments d'interface sans épaisseur [11]. Dans le cas tridimensionnel, on utilise la version dite avec effet de retard. Cette version du modèle permet d'assurer l'objectivité de la réponse jusqu'à la phase ultime de la rupture [3, 8].

Les caractéristiques de l'interface identifiées dans [10] sont rappelés dans le tableau I ci-dessous.

À titre d'exemple de l'influence d'un effet de retard, nous considérons l'évolution de la simulation de l'essai DCB lorsque le temps caractéristique évolue. Les coefficients optimum sont déterminés tels que les résultats soient identiques pour les modèles avec ou sans effet de retard dans le cas de la fissure « unidimensionnelle » décrit précédemment. Dans un premier temps, nous prenons $m = 1, n = 0,5$. En prenant $k = 100 \text{ min}^{-1}$, les résultats post-pic sont identiques pour les deux modèles et sont comparables à l'expérience réalisée avec une vitesse de chargement de $2 \text{ mm} \cdot \text{min}^{-1}$ (figure 11).

<table><caption>TABLEAU I</caption><thead><tr><td>Interface</td><td>Y<sub>c</sub>, N·mm<sup>-1</sup></td><td>y<sub>1</sub></td><td>α</td></tr></thead><tbody><tr><td>0°/0°</td><td>0,113 ± 0,007</td><td>0,37 ± 0,15</td><td>2,0 ± 0,3</td></tr><tr><td>±22,5°</td><td>0,167 ± 0,013</td><td>0,36 ± 0,17</td><td>1,4 ± 0,3</td></tr><tr><td>±45°</td><td>0,192 ± 0,014</td><td>0,44 ± 0,16</td><td>1,2 ± 0,1</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 10. Simulation de l'essai DCB : charge en fonction de l'ouverture en tête entre les deux bras, indication des longueur de délaminage.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Simulation de l'essai ENF : charge en fonction de l'ouverture en tête entre les deux bras.</figcaption></figure>

## 7. SIMULATION D'ESSAIS D'INITIATION

L'étude expérimentale de l'initiation d'un délaminage est souvent réalisée à l'aide d'éprouvettes lisses chargées uniaxialement de type EDT (edge delamination tension) [17]. La mécanique de la rupture n'est pas adaptée alors pour analyser un tel essai puisque le taux de restitution d'énergie devient nul pour une longueur initiale de fissure nulle. De plus, le délaminage est souvent en interaction avec les mécanismes d'endommagement des couches, particulièrement avec les fissures transverses [19] la méso-modélisation couches-interfaces endommageables présente alors tout son intérêt. Nous présentons ici la dégradation correspondante au délaminage et à l'endommagement des couches. On peut voir l'évolution du front de délaminage. Notons que la forme du front est réaliste par rapport aux essais réalisés à l'Aérospatiale [10].

<figure><img src="image_13.png" /><figcaption>Figure 12. Simulation de l'essai MMF : charge en fonction de l'ouverture en tête entre les deux bras.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Simulation de l'essai CLS : charge en fonction de l'ouverture en tête entre les deux bras.</figcaption></figure>

Nous considérons ici l'analyse d'une plaque stratifiée $[0_3/(\pm45)_2/90]_s$ constituée du matériau M55J/M18 sollicitée en traction. Le comportement des interfaces est endommageable et le comportement des plis, constitués d'éléments CUB8, est élasto-plastique endommageable.
<figure><img src="image_15.png" /></figure>

**Figure 14.** Simulation DCB, pour différentes valeurs du temps caractéristique. En prenant des valeurs de *k* supérieures 100 min⁻¹, les résultats de la simulation sont identiques à l'essai.

<figure><img src="image_16.png" /></figure>

**Figure 15.** Simulation d'un essai d'initiation du délaminage sur une éprouvette constituée de M55J/M18, la stratification est [0, (±45)₂, 90]. Contrainte dans l'axe de l'éprouvette en fonction des déformations.

L'endommagement dans les couches est principalement localisé pres des bords.

## 8. DÉLAMINAGE PAR FLAMBAGE LOCAL

Des imperfections peuvent exister dans les structures composites stratifiées aéronautiques. Elles peuvent être dus aux chutes d'outils par exemples lors de phase de maintenances. Dans cette partie, on montre les possibilités du code éléments finis développé dans cette étude sur l'exemple d'une éprouvette de compression possédant une zone interne délaminée. Sur cet exemple, les calculs sont menés sous l'hypothèse des grands déplacements. L'évolution du délaminage durant le flambage

<figure><img src="image_17.png" /></figure>

**Figure 16.** Comparaisons entre la réponse éléments finis et l'expérience. Le délaminage se produit sur l'interface centrale de l'éprouvette. La photographie rayon X confirme la forme de la propagation du délaminage. Le problème éléments finis est équivalent à l'essai que loin des fixations comme le montre la photographie.

<figure><img src="image_18.png" /></figure>

**Figure 17.** Essais de compression, il y a propagation du délaminage et flambage.

de l'éprouvette est observée. Sous l'action d'une sollicitation de compression, un flambage local de l'éprouvette peut être observé avant la phase de flambage globale.

## 9. CONCLUSIONS

L'extension aux structures tridimensionnelles des modèles d'interface déjà utilisée par [5, 7] a nécessité l'utilisation de modèles d'évolution de l'endommagement avec effet de retard [8]. Ces modèles qui ont un sens physique fort régularisent le problème initial et permettent d'obtenir jusqu'à la phase ultime de la rupture des résultats objectifs vis-à-vis de la discrétisation éléments finis. Les simulations éléments finis peuvent être conduites en tenant compte des non linéarités géométriques. La
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Maillage de la moitié de l'éprouvette. Il y a propagation du délaminage et flambage locale et globale de l'éprouvette sous l'action du chargement.</figcaption></figure>

simulation des essais de type mécanique de la rupture sur des éprouvettes stratifiées [$(+θ/-θ)_{4s}/(-θ/+θ)_{4s}$], pour des angles différents de 0° nécessite d'introduire les mécanismes d'endommagement des couches.

## RÉFÉRENCES

[1] Ladeveze P., Sur la mécanique de l'endommagement des composites, Comptes rendus des JNC, Vol. 5, Bathias C., Menkès D. (eds.), Pluralis Publication, Paris, 1986, pp. 667-683.

[2] Ladeveze P., Inelastic strains and damage, in: Damage Mechanics of Composite Materials, Talreja R. (ed.), Elsevier, 1994, pp. 117-136.

[3] Ladeveze P., A damage computational approach for composites : Basic aspects and micromechanical relations, Computational Mechanics 17 (1995) 142-150.

[4] Allix O., Ladeveze P., Interlaminar interface modelling for the prediction of laminates delamination, Composite Structure 22 (1992) 235-242.

[5] Allix O., Ladeveze P., Corigliano A., Damage analysis of interlaminar fracture specimen, J. Composite Structures 31 (1) (1995) 61-74.

[6] Ladeveze P., Le Dantec E., Damage modelling of the elementary ply for laminated composites, Comp. Sci. Tech. 43 (3) (1992) 257-268.

[7] Daudeville L., Ladeveze P., Une méthode simplifiée pour l'analyse du délaminage auprès d'un bord droit, Revue Européenne des éléments finis 3 (1) (1992) 127-150.

[8] Gornet L., Simulation des endommagements et de la rupture dans les composites stratifiés, thèse de l'université Paris 6/LMT/ENS Cachan, 1996.

[9] Gornet L., Implantation dans le code éléments finis Castem 2000 de modèles de comportement dédiés à l'analyse du délaminage de structures composites stratifiées, Rapport CEA, à paraître 1999.

[10] Leveque D., Analyse de la tenue au délaminage des composites stratifiés : identification d'un modèle d'interface interlaminaire, thèse de ENS Cachan/LMT, 1998.

[11] Beer G., An isoparametric joint-interface element for finite element analysis, Internat. J. Numer. Methods Engrg. 21 (1985) 585-600.

[12] Matthies, Strang, The solution of nonlinear equations, Internat. J. Numer. Methods Engrg. 14 (1979) 1613-1626.

[13] Russell A.J., Street K.N., Moisture and temperature effects on the mixed-mode delamination fracture of unidirectional graphite/epoxy, in: Delamination and Debonding of Materials, ASTM STP, Vol. 876, Johnson W.S. (ed.), Philadelphia, 1985, pp. 349-370.

[14] Carlsson L., Gillespie J.W., Pipes R.B., On the analysis and design of end notched flexure (ENF) specimen for mode II testing, J. Composite Materials 20 (1986) 594-604.

[15] Davies P., Measurement of GIC and GIIc in carbon/epoxy composites, Comp. Sci. Technol. 39 (1990) 193-205.

[16] Laksimi A., Benzeggagh M.L., Jing G., Hecini M., Roelandt J.M., Mode I interlaminar fracture of symmetrical cross-ply composites, Comp. Sci. Tech. 41 (1991) 147-164.

[17] Allix O., Ladeveze P., Leveque D., Perret L., Identification and validation of an interface damage model for delamination prediction, in: Proc. Complas 5, Barcelone, Espagne, 1996, pp. 1139-1147.

[18] O'Brien T.K., Characterization of delamination onset and growth in a composite laminate, in: Damage in Composite Materials, ASTM STP, Vol. 775, Reifsnider K.L. (ed.), 1982, pp. 140-167.

[19] Crossman F.W., Wang A.S.D., The dependence of transverse cracking and delamination on ply thickness in graphite/epoxy laminates, in: Damage in Composite Materials, ASTM STP, Vol. 775, Reifsnider K.L. (ed.), 1982, pp. 118-139.

[20] Ladeveze P., A damage computational method for composite structures, J. Computer Structure 44 (1/2) (1992) 79-87.

[21] Schellekens J.C., de Borst R., A non-linear finite element approach for the analysis of mode I-free edge delamination in composites, Internat. J. Solids Structures 30 (9) (1993) 1239-153.

[22] Talreja R., Transverse cracking and stiffness reduction in composite laminates, J. Composite Materials 19 (1985) 355-375.

[23] Herakovich C.T., Aboudi J., Lee S.W., Strauss E.A., Damage in composite laminates: Effects of transverse crack, Mech. Materials 7 (1987) 91-107.