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

HAL
open science

# Un modèle mixte-hybride du premier ordre sans facteurs correctifs de cisaillement transversal pour les structures composites

Achraf Tafla, Rezak Ayad, Nabil Talbi

► To cite this version:

Achraf Tafla, Rezak Ayad, Nabil Talbi. Un modèle mixte-hybride du premier ordre sans facteurs correctifs de cisaillement transversal pour les structures composites. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01493132

HAL Id: hal-01493132

https://hal.science/hal-01493132

Submitted on 21 Mar 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.

Public Domain
Un modèle mixte-hybride du premier ordre
sans facteurs correctifs de cisaillement
transversal pour les structures composites

Achraf TAFLA, Rezak AYAD et Nabil TALBI

Groupe de Mécanique des Matériaux et des Structures, ESIEC, Université de Reims
Champagne-Ardenne, BP 1029, 51686 Reims cedex 2
rezak.ayad@univ-reims.fr

RESUME. Cette étude concerne la formulation d'un modèle « élément fini » mixte-hybride naturel pour l'analyse linéaire des plaques et coques composites multicouches. Le modèle proposé, NHMiSP4/ml (Natural Hybrid-Mixed Shear Projection 4-node/multilayer) est basé sur la théorie linéaire du 1<sup>er</sup>ordre (Reissner/Mindlin). Il s'agit d'une adaptation au cas des structures stratifiées et sandwich du modèle MiSP4 isotrope [Ayad,1993]. Grâce à cette adaptation, le modèle NHMiSP4/ml permet la modélisation des plaques et coques composites mono et multicouches sans recourir aux facteurs de correction du CT. Il est de forme quadrilatérale à 4 nœuds et 6 ddl/nœud. Il est libre de tout verrouillage en CT et présente des performances de précision appréciables, comparativement aux éléments d'ordre supérieur.

ABSTRACT. This study deals with the formulation of a hybrid-mixed finite element model for linear analysis of multilayered composite plates and shells. The proposed model, NHMiSP4/ml (Natural Hybrid-Mixed Shear Projection 4-node/multilayer) is based on the Reissner-Mindlin fisrt order theory. It has been derived from the isotropic hybrid-mixed model MiSP4[Ayad,1993]. It can be used to simulate the behavior of multilayered plates and shells, without using shear correction factors. The element has been formulated with a quadrilateral shape and has 4 nodes and 6 dof per node. It's free of shear locking and gives good performances in terms of accuracy, in comparison with some high order elements.

MOTS CLES: Plaque et coque composite, Mindlin, élément fini, formulation mixte-hybride

KEYWORDS : Composite plate and shell, Mindlin, finite element, hybrid-mixd formulation

1. Modélisation par éléments finis des structures composites multicouches

Une théorie applicable à tous les cas possibles (matériau composite, anisotrope,
multicouches, sandwich etc...), qui serait simple et moins coûteuse en temps de
calcul est difficile à envisager. La théorie du 1<sup>er</sup> ordre de Mindlin (Mindlin, 1951)
tient compte des effets de CT à travers l'épaisseur. Elle conduit à des contraintes de
CT constantes dans l'épaisseur, en contradiction avec une représentation
quadratique, classiquement définie pour les poutres ou les plaques en flexion/CT.
Pour corriger cette insuffisance, des facteurs k de correction y sont introduits. Les
éléments du 1<sup>er</sup> ordre formulés en déplacement donnent généralement de bons
résultats pour les structures isotropes et orthotropes. Ils deviennent moins précis une
fois appliqués aux matériaux composites à plusieurs couches. Certes, l'introduction
des facteurs ki dans un modèle déplacement du 1<sup>er</sup> ordre a permis de résoudre des
problèmes de structures multicouches, mais leur évaluation dépend du nombre de stratifications. Pour écarter à jamais cet obstacle, des théories d'ordre supérieur ont été introduites. Une synthèse bibliographique leur a été consacrée dans la référence (Tafla, 2007). Nous proposons dans ce travail de développer, à l'aide d'une formulation mixte-hybride au sens de Hellinger-Reissner, un modèle « élément fini » de 1<sup>er</sup> ordre, sans recourir aux facteurs de correction de CT.

## 2. Présentation du modèle mixte-hybride naturel multicouches NHMiSP4/ml

Le modèle mixte-hybride naturel projeté en cisaillement NHMiSP4/ml (Mixed with Shear Projection 4-node/MultiLayer) est un élément quadrilatéral à 4 nœuds et 6 ddls/nœud (figure 1). Il est enrichi en membrane par une approche hybride (améliorée) au sens de Pian. L'approximation des contraintes de CT {τ_s} (continuité C⁻¹) est quadratique dans la direction de l'épaisseur(ζ). Elle est linéaire pour les contraintes de flexion {σ_s}. Dans le repère (ξ,η), {σ_s} est bilinéaire avec 12 paramètres {α_f}, {τ_s} est liée à {σ_s} via deux équations d'équilibre. Ces approximations sont en accord avec la condition de stabilité de Babutska. Pour éliminer le verrouillage en CT, nous définissons les déformations de CT {γ_s} indépendamment de l'approximation des variables cinématiques. Nous adoptons pour cela l'approche des déformations de substitution (Assumed Natural Strain).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Présentation du modèle mixte-hybride naturel multicouches NHMiSP4/ml</figcaption></figure>

### 2.1. Formulation variationnelle mixte-hybride naturelle projetée en cisaillement.

Le principe des travaux virtuels, sous sa forme mixte-hybride, est développé en prenant en compte le caractère multicouche d'une coque à n couches. Nous écrivons

$$ W = \sum_{\text{éléments}} (W_{\text{int}}^{e} - W_{\text{ext}}^{e}) = 0 \quad \forall \vec{u}_{q}^{*} \text{ avec } \vec{u}_{q} = \bar{\vec{u}}_{q} \text{ et } \vec{u}_{q}^{*} = \bar{0} \text{ sur } S_{u} \quad [1] $$
Avec

$$
\bar{u}_q^* = \begin{cases} U^* \\ V^* \\ W^* \end{cases} + \zeta \frac{h}{2} \begin{bmatrix} 0 & n_Z & -n_Y \\ -n_Z & 0 & n_X \\ n_Y & -n_X & 0 \end{bmatrix} \begin{bmatrix} \theta_X^* \\ \theta_Y^* \\ \theta_Z^* \end{bmatrix} \quad [2]
$$

$$
W_{\text{int}}^{e} = W_{\text{mf}}^{e} + W_{c}^{e} \quad (\text{travail virtuel interne}) \qquad [3]
$$

$$
W_{mf}^{e} = \int_{V^e} \langle \epsilon_s^* > \{\sigma_s\} + \langle \sigma_s^* > \{\epsilon_s\} - \langle \sigma_s^* > [H]^{-1} \{\sigma_s\} \rangle dV \quad (\text{travail de membrane-flexion}) \quad [4]
$$

$$
W_c^e = \int_{V^e} \langle \gamma_s^* > \{\tau_s\} + \langle \tau_s^* > \{\gamma_s\} - \langle \tau_s^* > [G]^{-1} \{\tau_s\} \rangle dV \quad (\text{travail de CT}) \quad [5]
$$

En adoptant des approximations linéaire et quadratique en $\zeta$ pour $\{\sigma_s\}$ et $\{\tau_s\}$ :

$$
\{\sigma_s\} = \{\sigma_0\} + \zeta\{\sigma_1\} \ ; \ \{\tau_s\} = (1 - \zeta^2)\{\tau_0\} \qquad [6]
$$

$$
\{\sigma_0\} = \begin{Bmatrix} \sigma_{x0} \\ \sigma_{y0} \\ \sigma_{xy0} \end{Bmatrix} \quad (\text{membrane}) ; \quad \{\sigma_1\} = \begin{Bmatrix} \sigma_{x1} \\ \sigma_{y1} \\ \sigma_{xy1} \end{Bmatrix} \quad (\text{flexion}) ; \quad \{\tau_0\} = \begin{Bmatrix} \tau_{xz0} \\ \tau_{yz0} \end{Bmatrix} \quad (\text{CT au pt p}(\zeta=0)) \quad [7]
$$

Pour les champs de déformations planes {$\epsilon_s$} et de CT {$\gamma_s$}, nous avons adopté un
choix simplifié de leur approximations en $\zeta$: {$\epsilon_s$} est linéaire et {$\gamma_s$} constant :

$$
\{\epsilon_s\} = \{\epsilon_0\} + \zeta\{\epsilon_1\} \ ; \quad \{\gamma_s\} = \{\gamma_0\} \qquad [8]
$$

Une intégration explicite suivant $\zeta$ des formes mixtes $W_{mf}^e$ et $W_c^e$ est réalisée dans un premier temps. Nous obtenons, en tenant compte des approximation [6] et [8] :

$$
W_{mf}^{e} = \int_{A^{e}} \left\{ \begin{array}{l} h \langle \sigma_{0}^{*} \rangle \{ \epsilon_{0} \} + h \{ \epsilon_{0}^{*} \} \langle \sigma_{0} \rangle + \frac{h}{3} \langle \epsilon_{1}^{*} \rangle \{ \sigma_{1} \} + \frac{h}{3} \langle \sigma_{1}^{*} \rangle \{ \epsilon_{1} \} \\ \\ - \langle \sigma_{0}^{*} \rangle \langle \bar{H}_{m} \rangle \{ \sigma_{0} \} \\ \\ - \frac{h^{4}}{36} \langle \sigma_{1}^{*} \rangle \langle \bar{H}_{f} \rangle \{ \sigma_{1} \} - \langle \sigma_{0}^{*} \rangle \langle \bar{H}_{mf} \rangle \{ \sigma_{1} \} - \langle \sigma_{1}^{*} \rangle \langle \bar{H}_{mf} \rangle^{T} \{ \sigma_{0} \} \end{array} \right\} dA \quad [9]
$$

$$
W_c^e = \int_{A^e} \left( \frac{2h}{3} \langle \gamma_0^* \rangle \{ \tau_0 \} + \frac{2h}{3} \langle \tau_0^* \rangle \{ \gamma_0 \} - \frac{h^2}{16} \langle \tau_0^* \rangle \langle \bar{H}_c \rangle \{ \tau_0 \} \right) dA \quad [10]
$$

Avec les matrices de comportement en membrane, flexion, couplage membrane-flexion et CT suivantes, adaptées pour notre modèle multicouches NHMiPS4/ml :

$$
[\bar{H}_m] = \sum_{i=1}^{nc} (z_{i+1} - z_i)[H]_i^{-1} ; \quad [\bar{H}_{mf}] = \frac{1}{h} \sum_{i=1}^{nc} (z_{i+1}^2 - z_i^2)[H]_i^{-1} ; \quad [\bar{H}_f] = \left( \sum_{i=1}^{nc} \frac{z_{i+1}^3 - z_i^3}{3} \right) [H]_i^{-1} \quad [11]
$$

$$
[\bar{H}_C] = \sum_{i=1}^{nc} [(z_{i+1} - z_i) + \frac{16}{5h^4}(z_{i+1}^5 - z_i^5) - \frac{8}{3h^2}(z_{i+1}^3 - z_i^3)][G]_i^{-1} \quad [12]
$$

La partie membrane est hybride naturelle au sens de Pian. Une approximation de
type C⁻¹ à 5 paramètres α_m des contraintes {σ₀} est utilisée : {σ₀}=[P₀] {α_m}. Les
contraintes de flexion {σ_I} sont bilinéaires en ξ, η avec 12 paramètres {α_f}. Les
contraintes de CT {$\tau_0$} sont définies en fonction des contraintes {$\sigma_l$} en utilisant les
deux premières équations de l'équilibre tridimensionnel en contraintes.

**Remarque :** Le modèle NHMiSP4/ml n'utilise pas de facteurs correctifs du CT pour les structures multicouches. L'expression de $W_c^e$[10] fait apparaître une nouvelle matrice $[\overline{H}_c]$ comme une somme de termes linéaire, cubique et d'ordre 5 en z associés aux différentes stratifications.

**2.2. Matrice de rigidité élémentaire**

Pour une stratification non-symétrique, l'expression finale de $W_{\text{int}}^e$ [3] s'écrit:

$$
W_{\text{int}}^{e} = \ll \alpha_{m}^{*} \gg \ll \alpha_{f}^{*} \gg \ll u_{n}^{*} \gg \begin{bmatrix} -[k_{mm}] & [k_{mf}] & [k_{mu}] \\ [k_{mf}]^{T} & [k_{fc}] & [k_{fu}] \\ [k_{mu}]^{T} & [k_{fu}]^{T} & [0] \end{bmatrix} \left\{ \begin{matrix} \{\alpha_{m}\} \\ \{\alpha_{f}\} \\ \{u_{n}\} \end{matrix} \right\} \quad [13]
$$

$$
\text{avec ;} \qquad [k_{mm}] = \int_A [P_0]^T [\overline{H}_m] [P_0] dA \quad ; \quad [k_{mu}] = h \int_A [P_0]^T [B_0] dA \qquad [14]
$$

$$
[k_{fc}] = -\frac{h^4}{36} \int_A [p_1]^T [\overline{H}_f] [p_1] dA - \frac{h^2}{16} \int_A [p_\tau]^T [\overline{H}_c] [p_\tau] dA \quad [15]
$$

$$
[k_{mf}] = \int_A [p_0]^T [\overline{H}_{mf}] [p_1] dA ; \quad [k_{fu}] = \int_A \left( \frac{h^2}{6} [p_\tau]^T [B_c] + \frac{h}{3} [p_1]^T [B_1] \right) dA \quad [16]
$$

La condensation statique au niveau local des paramètres {$\alpha_m$} et {$\alpha_f$} se traduit par :

$$
\{\alpha_m\} = [k_{mm}]^{-1} ([k_{mu}\{u_n\} - [k_{mf}]\{\alpha_f\}) \ , \ \{\alpha_f\} = [k_f]^{-1} [k_{ufi}]^T \{u_n\} \quad [17]
$$

La matrice de rigidité élémentaire finale s'écrit ainsi

$$
[k^e] = [A_1] + [A_2][A_3]^{-1}[A_2]^T ; \text{ avec } [A_1] = [k_{fu}]^T + [k_{mu}]^T[k_{mm}]^{-1}[k_{mu}] \quad [18]
$$

$$
[A_2] = [k_{mu}]^T [k_{mm}]^{-1} [k_{mf}] ; \quad [A_3] = [k_{mf}]^T [k_{mm}]^{-1} [k_{mf}] \quad [19]
$$

Un schéma de Gauss à 4 points permet une intégration exacte de toutes les matrices.

**3. Performance numérique du modèle multicouches NHMiSP4/ml**

**3.1. Plaque composite à 9 couches sous chargement doublement sinusoïdal**

Cet exemple est proposé par (Pagano&al,1972). Le matériau de base est un
composite unidirectionnel fortement orthotrope. Une stratification à 9 couches
(0/90/0/90/0/90/0/90/0) est étudiée (fig.2). Les épaisseurs de couches à 0° et 90° sont
égales. Les couches de même orientation ont toutes la même épaisseur. Nous avons
effectué les calculs pour deux types de maillage : 6x6 et 10x10.
<figure><img src="image_3.png" /><figcaption>Figure 2. Plaque composite à 9 couches</figcaption></figure>

**Géométrie :**

L=1000, h=250,100,20,0.1

**Matériau :** (orthotrope)

E1=25, E2=1., G12=0.5, G13=0.5

G23=0.2, v=0.25

**Conditions aux limites :**

w=βₓ=0 sur AB ;

βₓ=0 sur BC ; βᵧ=0 sur CD ;

w=βᵧ=0 sur DA

**Chargement doublement sinusoïdal :**

q=q₀sin(πx/L)sin(πy/L)

Les résultats du déplacement transversal *w* au point C (tableau 1), comparés avec ceux de (Lardeur,1990) (élément DSQ) et de (Batoz&al,1992)(élément Q4γWT, sont en bon accord avec la solution 3D (Pagano&al,1972). Nous avons constaté, de la même façon que (Pagano&al,1972), que plus le nombre de couches est grand, moins la structure est sensible à l'effet de CT. L'estimation des contraintes planes σₓₓ et σᵧᵧ montre que, pour L/h≥10, elles sont plus proches des solutions de référence., contrairement aux valeurs correspondantes obtenues pour le cas épais L/h=4. Les résultats des contraintes σₓₓ et σᵧᵧ sont globalement satisfaisants. La Figure 3 montre la distribution (pour L/h=10) de σₓₓ à travers l'épaisseur. Une très bonne corrélation avec les résultats exacts (Pagano&al,1972) est en effet observée.

<table><caption>Tableau 1. Plaque carrée composite à 9 couches</caption><thead><tr><th>L/h</th><th>Modèle</th><th>σ&#773;<sub>x</sub><br>(L/2,L/2,h/2)</th><th>σ&#773;<sub>y</sub><br>(L/2,L/2,h/4)</th><th>σ&#773;<sub>xz</sub><br>(0,L/2,0)</th><th>σ&#773;<sub>yz</sub><br>(L/2,0,0)</th><th>w&#773;<br>(L/2,L/2,0)</th></tr></thead><tbody><tr><td rowspan="5">4</td><td>DSQ 6x6</td><td>0.491</td><td>0.487</td><td>0.235</td><td>0.243</td><td>4.235</td></tr><tr><td>NHMiSP4/ml 6x6</td><td>0.455</td><td>0.536</td><td>0.22</td><td>0.265</td><td>4,133</td></tr><tr><td>NHMiSP4/ml10x10</td><td>0.453</td><td>0.534</td><td>0.221</td><td>0.265</td><td>4,121</td></tr><tr><td>Q4γWT(6x6)</td><td>0.468</td><td>0.485</td><td>0.24</td><td>0.239</td><td>4.29</td></tr><tr><td>Q4γWT(10x10)</td><td>0.472</td><td>0.506</td><td>0.241</td><td>0.240</td><td>4.281</td></tr><tr><td rowspan="5">10</td><td>Elasticité</td><td>0.684</td><td>0.628</td><td>0.223</td><td>0.223</td><td>4.079</td></tr><tr><td>DSQ 6x6</td><td>0.519</td><td>0.455</td><td>0.246</td><td>0.228</td><td>1.516</td></tr><tr><td>NHMiSP4/ml 6x6</td><td>0.505</td><td>0.476</td><td>0.242</td><td>0.237</td><td>1,512</td></tr><tr><td>NHMiSP4/ml10x10</td><td>0.505</td><td>0.472</td><td>0.243</td><td>0.238</td><td>1,510</td></tr><tr><td>Q4γWT(6x6)</td><td>0.506</td><td>0.458</td><td>0.257</td><td>0.220</td><td>1,529</td></tr><tr><td rowspan="5">50</td><td>Elasticité</td><td>0.509</td><td>0.461</td><td>0.258</td><td>0.220</td><td>1,530</td></tr><tr><td>DSQ 6x6</td><td>0.551</td><td>0.477</td><td>0.247</td><td>0.226</td><td>1.512</td></tr><tr><td>NHMiSP4/ml 6x6</td><td>0.538</td><td>0.432</td><td>0.253</td><td>0.216</td><td>1.015</td></tr><tr><td>NHMiSP4/ml10x10</td><td>0.54</td><td>0.435</td><td>0.256</td><td>0.218</td><td>1,019</td></tr><tr><td>Q4γWT(6x6)</td><td>0.536</td><td>0.435</td><td>0.258</td><td>0.220</td><td>1,018</td></tr><tr><td rowspan="5">50</td><td>Elasticité</td><td>0.531</td><td>0.428</td><td>0.269</td><td>0.206</td><td>1,016</td></tr><tr><td>DSQ 6x6</td><td>0.535</td><td>0.431</td><td>0.270</td><td>0.207</td><td>1,020</td></tr><tr><td>NHMiSP4/ml 6x6</td><td>0.539</td><td>0.433</td><td>0.258</td><td>0.219</td><td>1.021</td></tr><tr><td>NHMiSP4/ml10x10</td><td>0.536</td><td>0.435</td><td>0.258</td><td>0.220</td><td>1,018</td></tr><tr><td>Q4γWT(10x10)</td><td>0.535</td><td>0.431</td><td>0.270</td><td>0.207</td><td>1,020</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 3. Plaque carrée à 9 couches. Distribution de σ<sub>xz</sub> selon l'épaisseur (L/h=10)**</figcaption></figure>

**4. Conclusion**

Dans ce travail, nous avons présenté un nouvel élément fini mixte-hybride naturel pour l'analyse des plaques et des coques composites multicouches. Il est basé sur la théorie du premier ordre de Reissner-Mindlin. Le modèle NHMiSP4/ml (quadrilatère à 4-nœuds) n'utilise pas de fonctions bulles (coûteuses en temps de calcul) pour assurer sa stabilité. Il prend en compte les effets du CT à travers l'épaisseur et est libre de tout verrouillage. Il n'utilise pas de facteurs de correction du CT et ce quelque soit le nombre de couches. Les travaux de (Ayad,2002) ont montré que, dans un formalisme variationnel mixte-hybride (Hellinger-Reissner), une approximation adéquate des contraintes de CT conduirait à une correction naturelle du CT sans pouvoir introduire le facteur de correction k. En effet, k=5/6 (cas isotrope) apparaît naturellement dans la formulation théorique; il devenait en conséquence possible d'envisager la même démarche pour des sections composites multicouches. C'est ce que nous avions développé dans le cadre de nos travaux récents (Tafla,2007).

**5. Bibliographie**

Mindlin R. « Influence of rotatory inertia and shear on flexural motion of isotropic, elastic plates », *Jou.App.Mech.*, Vol 18, 31-38 (1951).

Tafla A., Eléments finis mixtes-hybrides naturels sans facteurs correctifs du cisaillement pour les plaques et les coques composites multicouches, Thèse de doctorat, Université de Reims Champagne-Ardenne, 2007.

Ayad R., Eléments finis de plaques et coques en formulation mixte avec projection en cisaillement, Thèse de doctorat, Université de Technologie de Compiègne, 1993.

Lardeur P., Développement et évaluation de deux nouveaux éléments finis de plaques et coques composites avec influence du cisaillement transversal, Thèse de doctorat, Université de Technologie de Compiègne, 1990.

Pagano N.J. et Hatfield S.J. « Elastic behaviour of multilayered bidirectional composites », *AIAA Journal*, Vol. 10, N°7, P.931-933, 1972

Batoz JL. Et Quesnel T., « Modélisation numérique du flambage de panneaux sandwich raidis », *La construction navale en composite*, Paris 7-9 décembre 1992.