# Une hiérarchie de modèles raffinés pour plaques composites: application aux effets de bord

Michele d'Ottavio, Olivier Polit

► To cite this version:

Michele d'Ottavio, Olivier Polit. Une hiérarchie de modèles raffinés pour plaques composites: application aux effets de bord. 9e colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01419939

HAL Id: hal-01419939

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

Submitted on 20 Dec 2016

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
# Une hiérarchie de modèles raffinés pour plaques composites: application aux effets de bord

M. D'Ottavio, O. Polit

Laboratoire Energétique, Mécanique, Electromagnétisme (LEME)
Université de Paris Ouest
50, rue de Sèvres
92410 Ville d'Avray, France
{michele.d_ottavio,olivier.polit}@u-paris10.fr

**Résumé** — La simulation numérique de plaques composites est souvent difficile à cause du contraste entre la simplicité des modèles permettant un calcul rapide, et la précision nécessaire pour une analyse des contraintes transverses, premières responsables du délaminage. Cette contribution présente des éléments finis plaque basés sur une hiérarchie de modèles dont la précision varie entre celle des modèles classiques plaque 2D à celle des modèles quasi-3D. L'application aux effets de bord libre démontre l'utilité de l'approche proposée.

**Mots clés** — matériaux composites, modèles raffinés plaque, effets de bord.

## 1 Introduction

Comme une grande partie de structures composites est réalisée sous forme de structures minces (plaques, coques), de nombreux modèles bidimensionnels (2D) ont été suggérés afin d'éviter l'utilisation de modèles numériques tridimensionnels (3D), souvent trop chers d'un point de vue computationnel [5]. Ces modèles approximés peuvent prendre en compte seulement partiellement ou presque entièrement les spécificités de la réponse mécanique des composites stratifiés. Sur cette base, on peut identifier une hiérarchie de modèles 2D, en partant des modèles les plus grossiers aptes à une représentation globale de la structure, jusqu'à ceux plus raffinés permettant de reprendre les distributions locales et les effets des interfaces entre couches de matériaux différents (voir figure 1) [1]. Une notation unifiée proposée par Carrera permet de formuler toute une série de modèles 2D d'une façon très avantageuse vis-à-vis de l'implémentation dans un code de calcul [2]. La formulation unifiée permet d'évaluer les effets des approximations 2D des modèles plus simples à partir de la solution quasi-3D des modèles plus raffinés. Cet article présente la hiérarchie des modèles plaque 2D en formulation unifiée et son application aux effets de bord libres des stratifiés, ce qui peut être considéré comme le paradigme de l'initiation du délaminage.

## 2 La formulation unifiée

Les modèles 2D sont exprimés en formulation unifiée en introduisant pour certaines variables des approximations qui définissent *a priori* le comportement dans l'épaisseur. On identifie trois niveaux d'approximation :

1. Choix des variables approximées indépendamment : dans une approche aux déplacements virtuels, les approximations sont introduites pour les déplacements $u_i$. En alternative, la
<figure><img src="image_2.png" /><figcaption>Figure 1 – Schéma de la réponse mécanique 3D d’un stratifié : distributions transverses des déplacements (à gauche), des contraintes de cisaillement transverse (au milieu), et des contraintes dans le plan (à droite).</figcaption></figure>

formulation variationnelle mixte partielle de Reissner [8] permet d'introduire des approximations indépendantes pour les déplacements et les contraintes transverses. Cela permet donc de définir un champ de contraintes transverse qui est a priori continu aux interfaces.

2. Choix de la description multicouche : on peut identifier l'approche couche équivalente (ESL) et l'approche couches discrètes (LW). Dans le premier cas, les approximations des variables sont exprimées pour tout l'empilement, et le modèle résultant a un nombre d'inconnues indépendant du nombre de couches. Dans le deuxième cas, les variables sont exprimées pour chaque couche séparément, et le nombre d'inconnues du modèle dépend du nombre de couches.

3. Choix de l'ordre du polynôme pour les variables : l'évolution des variables selon l'épaisseur est représentée par des polynômes dont l'ordre peut varier entre 1 (modèle linéaire) et N (modèles raffinés). Dans cette étude, l'ordre maximum sera fixé à N = 4.

## 2.1 Les approximations polynomiales

On considère une plaque composite à $N_l$ couches, comme celle illustrée en figure 2. Pour la suite, on utilisera la variable générique $g(x,y,z)$ pour indiquer la distribution des champs de déplacements $u_i(x,y,z)$ (avec $i \in [1,3]$) et de contraintes $\sigma_q(x,y,z)$ (avec $q \in [1,6]$ en adoptant la notation compacte de Voigt). La construction du modèle 2D se base sur la séparation du composition dans le plan $\Omega = x \times y$ de celui dans l'épaisseur $z$ (voir figure 2), et sur l'introduction de l'approximation polynômiale suivante :

$$g(x,y,z) = F_r(z) \hat{g}_r(x,y) \quad \text{avec } r = 0, 1, \dots, N. \qquad (1)$$

Différents polynômes $F_r(z)$ seront choisis en fonction de la description utilisée pour le stratifié.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Stratifié avec $N_l = 4$ couches : systèmes de coordonnées utilisés et effet de la fonction Zig-Zag superposée à une interpolation parabolique ESL.</figcaption></figure>
## 2.2 Modèles couche équivalente (description ESL)

Dans une description à couche équivalente, le comportement du stratifié est défini par le comportement de son plan de référence $\Omega$ à $z = 0$, que l'on pourra identifier au plan moyen du stratifié (fig. 2). En utilisant les polynômes de Taylor $F_r(z) = z^r$ dans l'eq. (1), où $z$ est la coordonnée globale suivant l'épaisseur de tout l'empilement, l'expression générale de l'approximation pour la variable générique $g$ s'écrit :

$$g(x, y, z) = z^r \hat{g}_r(x, y) \quad \text{avec } r = 0, 1, \dots, N. \qquad (2)$$

Pour permettre la discontinuité de la pente aux interfaces, un terme de l'expansion polynômiale peut être remplacé par la fonction Zig-Zag proposée par Murakami [4] :

$$F_{ZZ}(z) = (-1)^k \zeta_k \hat{g}_{ZZ}. \qquad (3)$$

On a utilisé la coordonnée locale adimensionnée $\zeta_k = \frac{2\zeta_k}{h_k}$, qui varie entre $-1$ (bottom, face inférieure de la couche) et $1$ (top, face supérieure de la couche). On a $z_k = z - z_{0k}$, où $z_{0k}$ correspond à la côte du plan moyen $\Omega_k$ de la couche $k$ (fig. 2). On notera que la variable $\hat{g}_{ZZ}$ est unique pour tout l'empilement, ce qui est conforme à l'approche ESL. La figure 2 montre qualitativement la discontinuité de la pente obtenue en superposant la fonction Zig-Zag à une approximation définie par des polynômes de Taylor.

## 2.3 Modèles couches discrètes (description LW)

Dans la description LW, la variable générique $g$ est approximée dans chaque couche $k$ séparément. Pour cela, on introduit l'indice $k$ et la coordonnée locale adimensionnée $\zeta_k$. Le modèle pour tout le stratifié est construit en assemblant les contributions des différentes couches. Pour faciliter l'imposition de la continuité aux interfaces de chaque couche $k$, on utilisera des polynômes avec des variables définies au top, $\hat{g}_t^k$, et au bottom, $\hat{g}_b^k$ :

$$g^k(x, y, z_k) = \frac{P_0(\zeta_k) + P_1(\zeta_k)}{2} \hat{g}_t^k(x, y) + \frac{P_0(\zeta_k) - P_1(\zeta_k)}{2} \hat{g}_b^k(x, y) + F_r(\zeta_k) \hat{g}_r^k(x, y), \quad (4)$$
avec $F_r(\zeta_k) = P_r(\zeta_k) - P_{r-2}(\zeta_k)$ et $r = 2, 3, \dots, N$.

L'interpolation est construite à partir des polynômes de Legendre $P_r(\zeta_k)$ définis par récurrence :

$$P_0(\zeta_k) = 1; \quad P_1(\zeta_k) = \zeta_k; \\ (n+1)P_{n+1}(\zeta_k) = (2n+1)\zeta_k P_n(\zeta_k) - n P_{n-1}(\zeta_k) \quad \text{pour } n = 1, 2, \dots, N-1. \qquad (5)$$

## 2.4 Formulation aux déplacements

Dans la formulation aux déplacements, on introduit l'approximation eq. (1) uniquement pour le champ de déplacements, $g = [u_i]$. Le modèle 2D est ainsi obtenu en utilisant le principe variationnel aux déplacements :

$$\int_{\Omega} \int_{h} \delta\varepsilon_{pG}^{T} \sigma_{pC} + \delta\varepsilon_{nG}^{T} \sigma_{nC} \, dz \, d\Omega = \delta W_{e}. \qquad (6)$$

On indique par $W_e$ le travail des forces extérieures, par $\delta$ l'operateur de variation virtuelle et par l'exposant $T$ l'operateur de transposition. L'intégrale sur l'épaisseur $h$ du stratifié est classiquement obtenue comme la somme des intégrales sur chaque couche $k$. Conformément à la notation du paragraphe suivant, on a séparé les champs de déformation et de contraintes dans leurs composantes dans le plan (index $p$) et transverse (index $n$) :

$$\varepsilon_p = [\varepsilon_1 \ \varepsilon_2 \ \varepsilon_6]^T; \quad \varepsilon_n = [\varepsilon_4 \ \varepsilon_5 \ \varepsilon_3]^T; \quad \sigma_p = [\sigma_1 \ \sigma_2 \ \sigma_6]^T; \quad \sigma_n = [\sigma_4 \ \sigma_5 \ \sigma_3]^T. \qquad (7)$$
L'indice G indique que les déformations sont obtenues par les relations géométriques linéaires qui relient les déplacements et les déformations. L'indice C rappelle que les contraintes sont obtenues à partir des déformations via la loi de comportement. Suivant la séparation donnée en eq. (7), on utilise la relation constitutive linéaire suivante :

$$ \begin{bmatrix} \sigma_{pC} \\ \sigma_{nC} \end{bmatrix}^k = \begin{bmatrix} \tilde{C}_{pp} & \tilde{C}_{pn} \\ \tilde{C}_{pn}^T & \tilde{C}_{nn} \end{bmatrix}^k \begin{bmatrix} \epsilon_{pG} \\ \epsilon_{nG} \end{bmatrix}^k. \quad (8) $$

Les coefficients $\tilde{C}^k$ décrivent en général la loi de Hooke 3D pour un matériau anisotrope et conti- nennent l'information de l'orientation d'éventuels axes principaux par rapport au repère utilisé pour la structure [7].

## 2.5 Formulation mixte partielle

Pour introduire des approximations indépendantes pour les champs de déplacements et de contraintes transverse (donc, $g = [u_i, \sigma_{nM}]$), on utilise le principe variationnelle proposé par Reiss- ner [8] :

$$ \int_{\Omega} \int_{h} \delta \epsilon_{pG}^{T} \sigma_{pC} + \delta \epsilon_{nG}^{T} \sigma_{nM} + \delta \sigma_{nM}^{T} (\epsilon_{nG} - \epsilon_{nC}) \, dz \, d\Omega = \delta W_e. \quad (9) $$

Le subscript M montre que les contraintes transverses sont modellisées indépendemment des dé- placements. Les déformations transverses $\epsilon_n$ peuvent ainsi être obtenues à partir des relations géometriques déplacements-déformations (index G), ou à partir de la loi constitutive (index C), qui est alors écrite sous la forme mixte suivante :

$$ \begin{bmatrix} \sigma_{pC} \\ \epsilon_{nC} \end{bmatrix}^k = \begin{bmatrix} C_{pp} & C_{pn} \\ -C_{pn}^T & C_{nn} \end{bmatrix}^k \begin{bmatrix} \epsilon_{pG} \\ \sigma_{nM} \end{bmatrix}^k. \quad (10) $$

Les relations entre les coefficients mixtes C et les coefficients $\tilde{C}$ utilisés dans la loi de Hooke généralisée peuvent être trouvés dans [3].

## 2.6 Nomenclature des modèles utilisés

Chaque modèle est obtenu par la combinaison de (i) la formulation variationnelle utilisée, (ii) la description du multicouche choisie (ESL ou LW) et (iii) l'ordre N de l'approximation polynômiale. Pour la suite, on pose $N \in [1, 4]$. Sur cette base, la nomenclature illustré dans la figure 3 est utilisé pour identifier univoquement chaque modèle.

<figure><img src="image_4.png" /><figcaption>Figure 3 – Schéma définissant la nomenclature des modèles implémentés.</figcaption></figure>

On note les points suivants :

- l'ordre N du polynôme est le même pour toutes les composants des champs approximés ;
- les déplacements peuvent être approximés par un approche couche équivalente (ESL ou ESL+ZZ), ou par un approche couches discrètes (LW);

- si la fonction Zig-Zag est utilisée, l'ordre de l'expansion polynômiale est $N-1$ : pour un modèle $E\{D,M\}ZN$, la fonction est donc d'ordre $z^{N-1}$;

- les contraintes transverse sont toujours approximées par couche discrètes (LW).

Le modèle le plus simple obtenu par la formulation unifiée est ED1, qui correspond au modèle FSDT avec prise en compte d'une déformation transverse normale $\epsilon_3$ constante. Le modèle le plus raffiné est LM4, qui a une variation des déplacements et des contraintes transverse décrite par un polynôme d'ordre 4 pour chaque couche. Le modèle EMZ2 a un champ de déplacement décrit ESL avec un terme linéaire et la fonction Zig-Zag, et un champ de contraintes transverse qui est décrit par couche (LW) par un polynôme quadratique. Ce modèle est le plus simple qui représente a priori le comportement Zig-Zag et la continuité interlaminaire des contraintes transverse.

### 3 Approximation éléments finis

Les modèles présentés ci-dessus ont été implémentés dans le code commerciale ABAQUS sous la forme d'éléments finis (EF) plaque via la programmation d'une user subroutine. L'approximation EF se base sur l'approche isoparamétrique et utilise un élément à 4 neuds avec interpolations bilinaires de type Lagrange [3]. Le tableau 1 présente un récapitulatif des éléments implémentés en spécifiant le comportement Zig-Zag (ZZ), la continuité interlaminaire (CI) des contraintes transverse et leur nombre de DDL par neud (DDL).

<table><caption>Tableau 1 – Récapitulatif des éléments finis implémentés pour une plaque composé par $N_l$ couches : réprésentation du comportement Zig-Zag (ZZ), de la continuité interlaminaire des contraintes transverses (CI) et leur nombre de dégrées de liberté par neud (DDL).</caption><thead><tr><th>Modèle</th><th>EDN</th><th>EDZN</th><th>LDN</th><th>EMN</th><th>EMZN</th><th>LMN</th></tr></thead><tbody><tr><td>ZZ</td><td>-</td><td>-</td><td>✓</td><td>-</td><td>✓</td><td>✓</td></tr><tr><td>CI</td><td>-</td><td>-</td><td>-</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>DDL</td><td>3(1 + N)</td><td>3(1 + N)</td><td>3(1 + N<sub>l</sub>N)</td><td>3(2 + N + N<sub>l</sub>N)</td><td>3(2 + N + N<sub>l</sub>N)</td><td>6(1 + N<sub>l</sub>N)</td></tr></tbody></table>

Pour faciliter l'utilisation de ces éléments finis, le pré-processeur graphique d'ABAQUS a été enrichi par un plug-in qui gère graphiquement la saisie des donnés nécessaires à la définition du problème. L'évaluation des résultats peut être effectuée par des logiciels de visualisation open-source (MEDIT, Gnuplot, ...), ou en important les résultats dans le post-processeur d'ABAQUS.

### 4 Application aux effets de bord

La série hiérarchique de modèles raffinés formulés en notation unifiée est appliquée à l'analyse EF des effets de bord dans une plaque composite [6]. Les résultats de la figure 4, obtenus par le modèle LM4, montre les importants gradients des contraintes de cisaillement transverse $\sigma_3$ et $\sigma_4$ qui caractérisent le bord libre $y = b$. On trouve un très bon accord avec les résultats de référence disponibles dans la littérature. Les résultats de la figure 5 montrent la capacité de la série hiérarchique à identifier les limitations des différentes approximations 2D.

### Références

[1] E. Carrera. Historical review of zig-zag theories for multilayered plates and shells. Appl. Mech. Rev., 56:287-308, 2003.
<figure><img src="image_5.png" /><figcaption>Figure 4 – Evolution vers le bord libre $y = b$ des contraintes interlaminaires $\sigma_3$ (à gauche) et $\sigma_4$ (à droite) à l'interface $z = h_k$ de deux stratifiés cross-ply symétriques.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5 – Comparaison de différents modèles 2D aux déplacement (à gauche) et mixtes (à droite) : contrainte de cisaillement transverse à l'interface $z = h_k$ du stratifié $[0°/90°]_s$.</figcaption></figure>

[2] E. Carrera. Theories and finite elements for multilayered plates and shells : A unified compact formulation with numerical assessment and benchmarking. *Arch. Comput. Meth. Eng.*, 10 : 215-296, 2003.

[3] E. Carrera and L. Demasi. Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 2 : Numerical implementations. *Int. J. Numer. Meth. Eng.*, 55 :253-291, 2002.

[4] H. Murakami. Laminated composite plate theory with improved in-plane response. *J. Appl. Mech.*, 53 :661-666, 1986.

[5] A. K. Noor and W. S. Burton. Assessment of shear deformation theories for multilayered composite plates. *Appl. Mech. Rev.*, 42 :1-13, 1989.

[6] R. B. Pipes and N. J. Pagano. Interlaminar stresses in composite laminates under uniform axial extension. *J. Compos. Mater.*, 4 :538-548, 1970.

[7] J. N. Reddy. *Mechanics of Laminated Composite Plates and Shells : Theory and Analysis*. CRC Press, Boca Raton, FL, 2 edition, 2004.

[8] E. Reissner. On a certain mixed variational theorem and a proposed application. *Int. J. Numer. Meth. Eng.*, 20 :1366-1368, 1984.