# Modélisation de plaques composites stratifiés et sandwichs par séparation de variables

Philippe Vidal, Laurent Gallimard, Olivier Polit

► To cite this version:

Philippe Vidal, Laurent Gallimard, Olivier Polit. Modélisation de plaques composites stratifiés et sandwichs par séparation de variables. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717031

HAL Id: hal-01717031

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

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Modélisation de plaques composites stratifiés et sandwichs par séparation de variables

Philippe VIDAL, Laurent GALLIMARD, Olivier POLIT

LEME EA 4416, Université Paris Ouest, France, {philippe.vidal,laurent.gallimard,olivier.polit}@u-paris10.fr

**Résumé** — Cet article porte sur la modélisation de plaques composites et sandwichs en exploitant une méthode par séparation de variables (Proper Generalized Decomposition). Les déplacements sont donc décomposés en une somme de produits de fonctions du plan x, y par des fonctions de la coordonnée transverse z. Ce travail se base sur la méthode proposée dans [8] et [3], et l'étend à des expansions d'ordre supérieur qui se montrent plus adaptées pour la modélisation de telles structures. On aboutit à une approche LayerWise d'ordre élevée avec une représentation par variables d'espaces séparées. Les résultats obtenus sont comparés à des solutions de référence exactes pour des sollicitations mécaniques.

**Mots clés** — composite, plaque, séparation de variables

## 1 Introduction

Les composites stratifiés et sandwichs sont de plus en plus utilisés, en particulier dans les structures aéronautique, maritime, ... grâce à leurs excellentes propriétés mécaniques. Dans cet environnement, ils sont souvent soumis à des sollicitations sévères de nature mécanique et/ou thermique. Or, les mécanismes de rupture dépendent fortement des effets locaux aux interfaces des couches ou à proximité des bords. Il s'avère donc nécessaire de les évaluer de façon précise.

Pour la modélisation des structures composites et sandwichs soumis à des chargements multi-physiques, on peut distinguer deux grandes approches :

- l'approche couche équivalente (ESL) : le nombre d'inconnues est indépendant du nombre de couches,

- l'approche couche discrète (LW) : le nombre d'inconnues est proportionnel au nombre de couches. Ces approches permettent d'éviter le recours à des modélisations 3D qui peuvent être lourdes en terme de mise en oeuvre et en temps de calculs.

Notons également que depuis une dizaine d'années, la méthode PGD (Proper Generalized Decomposition) a montré d'intéressantes propriétés dans le cadre de réduction de modèles [1]. Elle a été proposée dans le cadre de la méthode Latin par Ladevèze en introduisant une décomposition radiale. Elle a pu être avantageusement exploitée en utilisant la séparation de variables (coordonnées x, y et z) pour résoudre des EDP multidimensionnelles [2] et pour traiter des problèmes de plaques composites [8] (Solution de type Navier) et [3] (solution EF).

En nous basant sur le travail développé dans [3, 8] pour la séparation des variables et sur les travaux de Carrera [4, 5] et Vidal [9, 10] pour l'évaluation des modèles raffinés, nous en déduisons une approche de type LW avec séparation de variables en utilisant une expansion en z d'ordre 4. Le but est de construire une approche efficace et précise pour la modélisation de structures composites stratifiés et sandwichs.
## 2 Modélisation LW d'ordre 4 par variables séparées

Considérons une plaque $\mathcal{P} = \Omega \times \Omega_z$ avec $\Omega_z = [-\frac{h}{2}, \frac{h}{2}]$ dans un système de coordonnées cartésiennes $(x, y, z)$. La coordonnée $z$ est relative à la direction de l'épaisseur de la plaque $h$. $\Omega$ est une région arbitraire du plan $(x, y)$ et définit le plan moyen de la plaque.

En s'appuyant sur la représentation du déplacement par variables séparées, on écrit $(u_1(x, y, z), u_2(x, y, z), u_3(x, y, z))$ comme la somme de $N$ produits de fonctions dans le plan par des fonctions de la coordonnée transverse :

$$ u = \begin{bmatrix} u_1(x, y, z) \\ u_2(x, y, z) \\ u_3(x, y, z) \end{bmatrix} = \sum_{i=1}^{N} \begin{bmatrix} f_1^i(z) v_1^i(x, y) \\ f_2^i(z) v_2^i(x, y) \\ f_3^i(z) v_3^i(x, y) \end{bmatrix} \quad (1) $$

où $(f_1^i, f_2^i, f_3^i)$ sont définies sur $\Omega_z$ et $(v_1^i, v_2^i, v_3^i)$ sur $\Omega$.

Dans ce travail, une approximation élément fini classique à 8 noeuds est utilisée sur $\Omega$, et une description LW d'ordre 4 est adoptée dans $\Omega_z$, cette dernière s'avérant particulièrement adaptée pour la modélisation de structures composites.

## 3 Résolution du problème

Pour un déplacement admissible $\delta\vec{u} \in \delta U$, le principe variationnelle est donné par :
Chercher $\vec{u} \in U$ (espace des déplacements admissibles) tel que

$$ - \int_{\mathcal{P}} [\varepsilon(\delta\vec{u})]^T [\sigma(\vec{u})] d\mathcal{P} + \int_{\mathcal{P}} [\delta u]^T [b] d\mathcal{P} + \int_{\partial \mathcal{P}_F} [\delta u]^T [t] d\partial \mathcal{P} = 0, \quad \forall \delta\vec{u} \in \delta U \quad (2) $$

où $[b]$ et $[t]$ sont les efforts volumiques et les efforts surfaciques appliqués sur $\partial \mathcal{P}_F$.

Les couples $(f_j^i, v_j^i)$ se calculent itérativement. On suppose que les $n$ premiers couples sont connus, et on calcule le $n+1$ème couple. On peut écrire :

$$ [u^{n+1}] = [u^n] + \begin{bmatrix} f_1(z) v_1(x, y) \\ f_2(z) v_2(x, y) \\ f_3(z) v_3(x, y) \end{bmatrix} \quad (3) $$

où $[u^n]$ est connu, et $(v_1(x, y), v_2(x, y), v_3(x, y))$, $(f_1(z), f_2(z), f_3(z))$ sont les fonctions à calculer.

En injectant ces expressions dans la formulation variationnelle Eq. (2), le problème à résoudre devient non-linéaire. Une technique classique de type point fixe est adoptée. Elle conduit à la résolution successive d'un problème 2D (associé au calcul de $[v]^T = [v_1, v_2, v_3]$) et d'un problème 1D (associé au calcul de $[f]^T = [f_1, f_2, f_3]$) jusqu'à obtenir la convergence.

L'algorithme du point fixe se déroule comme suit :

- initialisation de $[f]$, $[\tilde{f}] \leftarrow [f]$

- $f$ étant connu et noté $\tilde{f}$
  résolution du système linéaire

$$ [K_z(\tilde{f})] [q^v] = [\mathcal{R}_v(\tilde{f}, u^n)] \quad (4) $$

$[q^v]$ étant les ddls associés à $[v]$.

- $v$ étant connu et noté $\tilde{v}$
  résolution du système linéaire

$$ [K_{xy}(\tilde{v})] [q^f] = [\mathcal{R}_f(\tilde{v}, u^n)] \quad (5) $$

$[q^f]$ étant les ddls associés à $[f]$.

Ce processus est arrêté dès que la distance entre deux solutions successives est inférieure à une valeur choisie par l'utilisateur. Le dernier couple $(\tilde{f}, \tilde{v})$ de ce processus est donc le $n+1$ème couple de la somme Eq. (3).
On remarque d'une part que le problème associé à f (Eq. (5)) est d'une taille très inférieure à celui associé à v (Eq. (4)). D'autre part, le coût de ce dernier problème 2D (Eq. (4)) est lui très inférieur à une approche LW classique, le nombre de ddls inconnus se limitant à ceux du maillage 2D du plan moyen.

# 4 Résultats numériques

Dans cette section, un EF quadrilatère à 8 noeuds basé sur les fonctions d'interpolation de Serendipity est utilisé pour les inconnues dépendant des coordonnées x, y. Un schéma d'intégration numérique avec 3 × 3 points permet d'évaluer les matrices élémentaires. Pour l'intégration des fonctions dépendant de la coordonnée transverse, une intégration analytique est réalisée.

Dans un premier temps, une comparaison de deux approches basées sur une expansion linéaire [3] et du quatrième ordre dans l'épaisseur est effectuée sur un exemple homogène simple. Elle vise à montrer la nécessité de choisir un ordre élevé pour la coordonnée transverse dans l'écriture des déplacements. Puis, cette approche est évaluée sur des exemples de plaques stratifiées symétriques, ainsi que sur un cas sandwich. Différents élancements sont considérés. L'influence des conditions limites est également abordée. Pour tous les exemples, la convergence du point fixe est rapide, seulement 3 itérations suffisent.

## 4.1 Influence du degré d'expansion en z

Les approches basées sur un développement de degré 1 développé dans [3] ($d° = 1$) et de degré 4 ($d° = 4$) pour le présent article sont comparées pour un test homogène. Pour cela, une plaque carrée de dimension $0.1m \times 0.1m$ avec un élancement de $S = 10$ simplement posée, et soumise à une pression bi-sinusoïdale est considérée. Elle est constituée d'un matériau isotrope tel que $E = 73GPa$ et $v = 0.34$. L'évaluation s'effectue à la fois sur des quantités locales et globales. La solution exacte 3D est donnée par [6].

Tout d'abord, une étude de convergence est effectuée. Le nombre de couches numériques pour le cas linéaire est fixé à $N_z = 4$, de sorte que le nombre de ddls de chaque approche est le même. La figure Fig. 1 montre les résultats obtenus pour l'énergie de déformation et les déplacements et contraintes. La vitesse de convergence est identique pour les 2 approches mais la précision diffère. La présente approche permet de converger vers l'énergie de déformation de référence. D'un point de vue local, l'erreur maximale constatée est de 1.2% pour un maillage $N_x = N_y = 16$. Pour l'approche linéaire, les erreurs sont plus importantes. Elles sont de 10% pour $\bar{\sigma}_{11}$, 6% pour $\bar{\sigma}_{13}$ et 500% pour $\bar{\sigma}_{33}$.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – étude de convergence en fonction de $N_x = N_y$ de l'énergie de déformation (gauche), $\bar{u}_3$, $\bar{\sigma}_{11}$, $\bar{\sigma}_{13}$ and $\bar{\sigma}_{33}$ (droite) - S=10 - plaque isotrope - $b = a$ - $N_z = 4$ avec $d°1$, $N_z = 1$ avec $d°4$</figcaption></figure>

Cette tendance est confirmée par la figure Fig. 2 qui donne la distribution de $\bar{\sigma}_{13}$ et $\bar{\sigma}_{33}$ dans l'épaisseur. L'approche adoptée ici donne la solution de référence. En revanche, l'expansion linéaire conduit à
des discontinuités de contraintes. Des oscillations apparaissent pour $\bar{\sigma}_{33}$, et les conditions sur les bords ne sont pas respectées. Le nombre de couches numériques n'est donc pas suffisant. Une étude est donc effectuée pour mesurer la vitesse de convergence par rapport à ce paramètre $N_z$. La figure Fig. 3 montre qu'un raffinement de 32 couches numériques est nécessaire pour obtenir une erreur faible pour l'énergie de déformation, $\bar{\sigma}_{11}$ et $\bar{\sigma}_{13}$. En revanche, la convergence sur $\bar{\sigma}_{33}$ est plus lente et nécessite un raffinement encore plus important.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – distribution de $\bar{\sigma}_{13}$ (gauche) and $\bar{\sigma}_{33}$ (droite) dans l'épaisseur - S=10 - plaque isotrope - $b = a$ - $N_z = 4$ avec $d^\circ 1$, $N_z = 1$ avec $d^\circ 4$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3 – étude de convergence en fonction de $N_z$ de l'énergie de déformation, $\bar{\sigma}_{11}$, $\bar{\sigma}_{13}$ et $\bar{\sigma}_{33}$ - S=10 - plaque isotrope - $b = a$ - $d^\circ 1$</figcaption></figure>

Dans la suite de l'étude, seule l'approche $d^\circ 4$ développée ici, qui s'est montrée plus efficace et plus précise, sera conservée.

## 4.2 plaque stratifiée symétrique en flexion sous pression localisée

L'exemple considéré est une plaque carrée $a = b$ stratifiée [0°/90°/0°] soumise à une pression localisée $q_0$ en son centre sur une surface $a/10 \times b/10$ (Cf. Fig. 4 droite). L'élancement $S = a/h$ varie de 4 à 100. Le matériau est tel que $E_L = 25$ GPa , $E_T = 1$ GPa , $G_{LT} = 0.2$ GPa , $G_{TT} = 0.5$ GPa, $v_{LT} = v_{TT} = 0.25$ où L désigne la direction des fibres, et T la direction transverse. Les résultats sont adimensionnés tel que : $\bar{u} = u_1(0,b/2,z) \frac{E_T}{hq_0S^3}$, $\bar{w} = u_3(a/2,b/2,z) \frac{100E_T}{S^4hq_0}$, $\bar{\sigma}_{\alpha\alpha} = \frac{\sigma_{\alpha\alpha}(a/2,b/2,z)}{q_0S^2}$, $\bar{\sigma}_{12} = \frac{\sigma_{12}(0,0,z)}{q_0S^2}$,
$\bar{\sigma}_{13} = \frac{\sigma_{13}(0, b/2, z)}{q_0 S}$, $\bar{\sigma}_{23} = \frac{\sigma_{23}(a/2, 0, z)}{q_0 S}$, $\bar{\sigma}_{33} = \frac{\sigma_{33}(a/2, b/2, z)}{q_0}$. Le maillage utilisant un space ratio est donné sur la figure Fig. 4 gauche. Il provient d'une étude de convergence non exposée ici pour des raisons de concision. Les résultats de référence sont obtenus avec 450 termes dans la série de Fourier en suivant la démarche proposée dans [7]. Cet exemple vise à montrer la capacité de la méthode à capturer des phénomènes locaux.

<figure><img src="image_5.png" /><figcaption>Fig. 4 – maillage $N_x = N_y = 22$ sr(14) (gauche) - pression localisée (droite)</figcaption></figure>

On constate que le nombre de couples varient de 5 pour $S = 100$ à 12 pour $S = 4$. Ce nombre augmente avec l'épaisseur de la plaque. Le tableau Tab. 1 montre que la précision des résultats est très bonne. L'erreur maximale est seulement de 2.8 % que ce soit pour des plaques épaisses ou minces. La répartition des déplacements et contraintes dans l'épaisseur (Cf. Fig. 5 et Fig. 6) confirme la précision des résultats obtenus par rapport à la solution de référence. La contrainte de flexion $\bar{\sigma}_{11}$ est fortement influencée par la pression localisée, et spécifiquement dans la couche supérieure. La répartition dans l'épaisseur est fortement non-linéaire et non-symétrique. L'approche permet de bien capturer cette spécificité. On note également que les conditions aux bords haut et bas sont respectées.

<table><thead><tr><th>S</th><th></th><th>u(h/2) × 10<sup>-4</sup></th><th>w(0)</th><th>σ<sub>11</sub>(h/2)</th><th>σ<sub>22</sub>(-h/6)</th><th>σ<sub>12</sub>(h/2) × 10<sup>-4</sup></th><th>σ<sub>13</sub>(0)</th><th>σ<sub>23</sub>(0)</th><th>σ<sub>33</sub>(0)</th></tr></thead><tbody><tr><td rowspan="3">4</td><td>PGD d° 4</td><td>-3.76</td><td>0.1196</td><td>0.2563</td><td>-0.0512</td><td>-15.5</td><td>0.0089</td><td>0.0029</td><td>0.2067</td></tr><tr><td>erreur</td><td>0.2 %</td><td>0.2 %</td><td>2.8 %</td><td>2.6 %</td><td>0.1 %</td><td>0.0 %</td><td>0.5 %</td><td>0.92 %</td></tr><tr><td>exact</td><td>-3.77</td><td>0.1198</td><td>0.2639</td><td>-0.0526</td><td>-15.6</td><td>0.0089</td><td>0.0029</td><td>0.2086</td></tr><tr><td rowspan="3">10</td><td>PGD d° 4</td><td>-2.97</td><td>0.0420</td><td>0.0757</td><td>-0.0427</td><td>-8.54</td><td>0.0115</td><td>-0.0016</td><td>0.4564</td></tr><tr><td>erreur</td><td>0.2 %</td><td>0.2 %</td><td>1.1 %</td><td>0.3 %</td><td>0.1 %</td><td>0.3 %</td><td>1.0 %</td><td>1.03 %</td></tr><tr><td>exact</td><td>-2.98</td><td>0.0421</td><td>0.0748</td><td>-0.0428</td><td>-8.55</td><td>0.0116</td><td>-0.0016</td><td>0.4611</td></tr><tr><td rowspan="3">40</td><td>PGD d° 4</td><td>-2.91</td><td>0.0223</td><td>0.0425</td><td>-0.0327</td><td>-5.85</td><td>0.0112</td><td>-0.0022</td><td>0.5091</td></tr><tr><td>erreur</td><td>0.1 %</td><td>0.1 %</td><td>1.0 %</td><td>0.5 %</td><td>0.0 %</td><td>1.1 %</td><td>0.9 %</td><td>1.04 %</td></tr><tr><td>exact</td><td>-2.90</td><td>0.0223</td><td>0.0429</td><td>-0.0325</td><td>-5.85</td><td>0.0111</td><td>-0.0023</td><td>0.5039</td></tr><tr><td rowspan="3">100</td><td>PGD d° 4</td><td>-2.90</td><td>0.0210</td><td>0.0406</td><td>-0.0303</td><td>-5.65</td><td>0.0110</td><td>-0.0023</td><td>0.5028</td></tr><tr><td>erreur</td><td>0.0 %</td><td>0.0 %</td><td>0.1 %</td><td>0.0 %</td><td>0.1 %</td><td>0.3 %</td><td>1.4 %</td><td>0.56 %</td></tr><tr><td>exact</td><td>-2.90</td><td>0.0210</td><td>0.0406</td><td>-0.0303</td><td>-5.65</td><td>0.0109</td><td>-0.0023</td><td>0.5000</td></tr></tbody></table>

Tableau 1 – 3 couches (0°/90°/0°) - pression localisée - $b = a - N_x = N_y = 22$ sr(14) - $N_z = NC$

## 4.3 plaque sandwich en flexion

L'approche est maintenant évaluée sur une plaque sandwich constitué de 2 faces ayant les mêmes propriétés matériau que la structure précédente et d'un coeur avec les propriétés suivantes : $E_{xx} = E_{yy} = 0.04$ GPa, $E_{zz} = 0.5$ GPa, $G_{xz} = G_{yz} = 0.06$ GPa, $G_{xy} = 0.016$ GPa, $v_{xz} = v_{yz} = 0.02$, $v_{xy} = 0.25$. L'épaisseur de chaque face vaut $h/10$. La structure est simplement appuyée et soumise à une pression bi-sinusoidale. Un seul couple suffit pour représenter la solution. Les figures Fig. 7 et Fig. 8 donnent la répartition des déplacements et contraintes dans l'épaisseur pour une plaque très épaisse. A nouveau, les résultats sont en accord avec la référence. On remarque que l'effet zigzag du déplacement dans le plan est bien représenté. L'évolution de la flèche dans l'épaisseur est non-linéaire. Les fortes variations des contraintes de
<figure><img src="image_6.png" /><figcaption>Fig. 5 – distribution de $\bar{u}_1$ (gauche), $\bar{u}_3$ (milieu) and $\bar{\sigma}_{13}$ (droite) dans l’épaisseur - S=4 - 3 couches - pression localisée - $d^\circ 4$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6 – distribution de $\bar{\sigma}_{11}$ (gauche), $\bar{\sigma}_{22}$ (milieu) et $\bar{\sigma}_{33}$ (droite) dans l’épaisseur - S=4 - 3 couches - pression localisée - $d^\circ 4$</figcaption></figure>

cisaillement transverse dans les faces sont également bien capturées. Les considérations physiques sur les contraintes normales et transverses sont respectées.

<figure><img src="image_8.png" /><figcaption>Fig. 7 – distribution de $\bar{u}_1$ (gauche), $\bar{u}_3$ (milieu) et $\bar{\sigma}_{11}$ (droite) dans l’épaisseur - S=2 - sandwich - $d^\circ 4$</figcaption></figure>

# 5 Conclusion

Dans cet article, la décomposition de la solution sous forme radiale est appliquée à la modélisation de plaques composites en utilisant une approche layerwise d'ordre 4 pour le déplacement. Elle est tout d'abord comparée à une approche linéaire sur un exemple homogène afin de montrer la pertinence du degré d'expansion en z retenu. L'approche s'avère non seulement plus précise, mais également plus efficace en terme de coût de calculs. L'approche linéaire nécessite en effet un raffinement important des couches numériques dans l'épaisseur pour atteindre un degré de précision comparable à l'approche présentée. Les résultats sur des structures plaques composite stratifiée et sandwich ont confirmé la qualité de la méthode en terme de précision sur des quantités locales pour des chargements globaux ou localisés. La capacité de l'approche à capturer des effets locaux et représenter des distributions complexes de déplacements ou contraintes a été prouvée, ces phénomènes pouvant être produits par le chargement, le type d'empilement ou le type de plaques composites.
<figure><img src="image_9.png" /><figcaption>Fig. 8 – distribution de $\bar{\sigma}_{13}$ (gauche), $\bar{\sigma}_{23}$ (milieu) et $\bar{\sigma}_{33}$ (droite) dans l’épaisseur - S=2 - sandwich - $d^\circ 4$</figcaption></figure>

Pour des empilements complexes, l’approche devient moins coûteuse qu’une méthode LW, et elle permet de raffiner le comportement suivant l’épaisseur à moindre coût.

Références

[1] O. Allix, P. Vidal. *A new multi-solution approach suitable for structural identification problems*. , Comput. Methods Appl. Mech. Eng., 191, 2727-2758, 2002.

[2] A. Ammar, B. Mokdada, F. Chinesta, R. Keunings. *A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids*, Journal of Non-Newtonian Fluid Mechanics, 139, 153-176, 2006.

[3] B. Bognet, F. Bordeu, F. Chinesta, A. Leygue, A. Poitou. *Advanced simulation of models defined in plate geometries : 3D solutions with 2D computational complexity*, Comput. Methods Appl. Mech. Eng., 201-204, 1-12, 2012.

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

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

[6] L. Demasi. $\infty^3$ plate theories for thick and thin plates : the generalized unified formulation, Compos. Struct. 84, 256-270, 2008.

[7] N. Pagano. *Exact solutions for rectangular bidirectional composites and sandwich plates*, J. Comp. Materials 4, 20-34, 1970.

[8] M. Savoia, J.N. Reddy. *A variational approach to three-dimensional elasticity solutions of laminated composite plates*, J. Applied Mech. ASME, 59, 166-175, 1992.

[9] P. Vidal, O. Polit. *A family of sinus finite elements for the analysis of rectangular laminated beams*, Compos. Struct., 84, 56-72, 2008.

[10] P. Vidal, O. Polit. *A refined sine-based finite element with transverse normal deformation for the analysis of laminated beams under thermomechanical loads*, J. Mech. Mater. Struct., 4, 1127-1155, 2009.