# Élément fini triangulaire de membrane avec degrés de liberté rotationnels

Serge Dubigeon, Bernard Peseux

► To cite this version:

Serge Dubigeon, Bernard Peseux. Élément fini triangulaire de membrane avec degrés de liberté rotationnels. Revue Européenne des Éléments Finis, 2003, 12 (1), pp.59-97. 10.3166/reef.12.59-97. hal-01380382

HAL Id: hal-01380382

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

Submitted on 13 Oct 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
# Elément fini triangulaire de membrane avec degrés de liberté rotationnels-élément fini triangulaire de coque

Serge Dubigeon — Bernard Peseux

Laboratoire de Mécanique et Matériaux
Division Mécanique des Structures
Ecole Centrale, 1, rue de la Noë
F-44321 Nantes cedex 3
bernard.peseux@ec-nantes.fr

RÉSUMÉ. L'assemblage de plaques non coplanaires dans l'espace pose depuis toujours des difficultés car les rotations normales au plan des plaques n'interviennent pas classiquement comme inconnues. Après un rappel de l'état de l'art, on propose un élément fini membranaire convergeant, faisant intervenir explicitement les degrés de liberté rotationnels, dits de « drilling ». Cet élément associé à un élément de plaque en flexion de Zienkiewicz et Specht possède trois déplacements et trois rotations par nœud. L'élément est validé en reprenant des tests standard de convergence de la littérature. Des résultats complémentaires de tests en statique et en dynamique présentés dans l'article montrent le bon comportement de cet élément.

**ABSTRACT.** For a long time, the coupling of non coplanar plates in space makes lot of difficulties because the normal rotations to the plane of the plate don't occur as unknown. After recalling the main approaches which are actually used to solve these difficulties, we propose a finite element which explicitely takes into account the « drilling » parameters and passes the « hard » patch test. This one, associated to a bending finite element from Zienkiewicz and Specht is caracterized by three displacements and three rotations per node. The validation of this element is shown with the results of standard convergence tests. Some complementar results in static and dynamic which indicate the good accuracy of this element are presented in this paper.

MOTS-CLÉS : élément fini de coque, assemblage de plaques, rotation normale.

KEYWORDS: Shell finite element, assembly of plates, normal rotation.
# 1. Introduction

L'étude d'assemblages de plaques non coplanaires ou de coques par la méthode des éléments finis peut se faire en utilisant des éléments tridimensionnels ou des éléments de coques à surface moyenne courbe. Dans ce cas, les analyses peuvent nécessiter des efforts de mise en œuvre importants ou conduire à des systèmes matriciels à résoudre de grande taille.

Si on veut utiliser des éléments plaques ordinaires, basés sur les modèles de Kirchhoff ou Hencky pour étudier ces assemblages, on se heurte à certaines difficultés. En effet, les éléments finis plaques ordinaires sont en général construits en considérant comme degrés de liberté nodaux les trois translations et deux rotations, la rotation normale au plan de la plaque n'intervient pas dans le calcul de l'énergie de déformation, donc n'a pas de contribution à la matrice de rigidité. Pour tourner cette difficulté on peut user d'expédients numériques, l'un des plus courants consiste à alimenter le degré de liberté rotation normale par une raideur numérique fictive. Cette solution, même si elle peut s'avérer efficace, n'est pas complètement satisfaisante. C'est pourquoi les travaux relatifs à ce problème sont nombreux et les différentes méthodes proposées pour introduire dans les éléments finis membranaires des degrés de liberté rotationnels sont basées sur l'utilisation de degrés de liberté :

- ou bien de « drilling » comme le font Bergan et Felippa [BER 85], Hughes et Brezzi [HUG 89], Ibrahimbegovic et al. [IBR 90, IBR 91, IBR 93], Liu et To [LIU 95], Hughes et al. [HUG 95],

- ou bien de « vertex » introduits par Allman [ALL 84, ALL 88a, ALL 88b], Yunnus et al. [YUN 89], Cook [COO 87, COO 91], Aminpour [AMI 92],

- ou encore de « fiber rotations » étudiés par Ayad et al. [AYA 95a, AYA 95b].

Notons enfin qu'Irons et Ahmad [IRO 80] faisaient déjà allusion à ce problème de rotations normales pour les éléments membranaires.

Dans cet article, nous proposons une méthode rigoureuse pour prendre en compte la rotation normale au plan de la plaque et construire la matrice de rigidité élémentaire. La démarche présentée est basée sur l'utilisation d'une part, de degrés de liberté rotationnels dits de « drilling », et, d'autre part, en construisant une interpolation particulière (« linked » interpolation) pour exprimer le champ des déplacements. Cette interpolation est de même nature que celle proposée par Zienckiewicz et al. [ZIE 93], Taylor et Aurichkio [TAY 93], Xu et al. [XU 94].

La rotation normale à la plaque est induite par les déformations membranaires et une première solution proposée par Allman, Felippa, Yunus, Ibrahimbegovic, Cook, Aminpour, Hugues, Liu, consiste en l'introduction de multiplicateurs de Lagrange lors des assemblages.

Nous proposons ici une autre solution aboutissant à la construction d'un élément triangulaire non conforme mais passant le patch test, et introduisant explicitement
les rotations normales. Cet élément membranaire peut être associé aux meilleurs éléments de plaque en flexion pour aboutir à un élément de coque performant

## 2. Un nouvel élément membranaire - Triangle à 3 nœuds avec degrés de liberté rotationnels (T3D)

Pour l'étude du champ des déplacements dans le plan de l'élément et la construction des fonctions d'interpolation spécifiques associées aux rotations normales, nous nous inspirons des travaux de Zienckiewicz [ZIE 91] et de Specht [SPE 88] utilisant les coordonnées barycentriques du triangle et la fonction produit $L_1L_2L_3$ de ces coordonnées, pour construire un élément fini triangulaire de flexion. Cet élément s'est avéré très performant pour l'étude des plaques basées sur le modèle de Kirchhoff.

### 2.1. Construction d'un champ de rotations normales pour le triangle à 3 nœuds

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Elément triangulaire avec degrés de liberté rotationnels</figcaption></figure>

On veut construire un champ de déplacements membranaire $\vec{u}(x, y, t) = u\vec{e}_1 + v\vec{e}_2$, associé à des rotations normales $\theta_{1z}, \theta_{2z}, \theta_{3z}$ aux nœuds 1, 2 et 3 respectivement.

Ce champ est cherché sous la forme générale :

$$ \vec{u}(x, y, t) = \langle \Phi \rangle \{A\} \vec{e}_1 + \langle \Phi \rangle \{B\} \vec{e}_2 \quad [1] $$
où <Φ> est le polynôme caractéristique de l'interpolation généralisée. Il est choisi
cubique en fonction des coordonnées barycentriques. Il est donc donné par :

$$
\langle \Phi \rangle = \langle L_1^3, L_2^3, L_3^3, L_1^2 L_2, L_1^2 L_3, L_2^2 L_1, L_2^2 L_3, L_3^2 L_1, L_3^2 L_2, L_1 L_2 L_3 \rangle \quad [2]
$$

`{A}, et {B}, sont des vecteurs de coordonnées généralisées, de composantes`
`A₁, A₂, A₃ pour {A} et B₁, B₂, B₃ pour {B}. Les coordonnées barycentriques sont :`

$$
\begin{align*}
L_1 &= \frac{1}{2A} (a_1 + b_1 x + c_1 y) \\
L_2 &= \frac{1}{2A} (a_2 + b_2 x + c_2 y) \\
L_3 &= \frac{1}{2A} (a_3 + b_3 x + c_3 y)
\end{align*}
$$

avec $a_1 = x_2 y_3 - x_3 y_2$, $b_1 = y_2 - y_3$, $c_1 = x_3 - x_2$. Les autres termes sont obtenus par
permutation circulaire.

L'aire du triangle est donnée par :

$$
2A = (x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1) = c_3 b_2 - c_2 b_3 \quad [4]
$$

On considère tout d'abord le champ induit par une rotation normale unitaire θ$_{z1}$
au nœud 1 (figure 1).

Le champ des déplacements doit satisfaire les conditions suivantes :

$$
\begin{align*}
\vec{u}(x, y, t) &= \vec{0} \quad \text{en } L_1 = 0 \text{ et } L_1 = 1 \quad (L_2 = L_3 = 0) \\
\theta_z &= \frac{1}{2} (\overrightarrow{rot} \vec{u}) \cdot \vec{z} = \frac{1}{2} (v_x - u_y) \tag{5} \\
\text{et } \theta_z &= 0 \quad \text{en } L_1 = 0, \quad \theta_z = 1 \quad \text{en } L_1 = 1
\end{align*}
$$

Les conditions $\vec{u} = \vec{0}$ en $L_1 = 0$ et $\vec{u} = \vec{0}$ en $L_1 = 1, L_2 = L_3 = 0$, limitent déjà le
choix des termes de $\langle \Phi \rangle$ à :

$$
\langle \Phi \rangle = \langle L_1^2 L_2, L_1^2 L_3, L_2^2 L_1, L_3^2 L_1, L_1 L_2 L_3 \rangle \quad [6]
$$

D'autre part en exprimant θz (équation [5]), il vient :

$$
\theta_z = \frac{1}{4A} \left( \sum_{i=1}^{3} b_i \frac{\partial v(L_1, L_2, L_3)}{\partial L_i} - c_i \frac{\partial u(L_1, L_2, L_3)}{\partial L_i} \right) \quad [7]
$$
La condition $\theta_z(0, L_2, L_3) = 0$ limite encore le choix de $\langle\Phi\rangle$ à $\langle\Phi\rangle = \langle L_1^2 L_2, L_1^2 L_3 \rangle$. Les vecteurs des coordonnées généralisées deviennent :

$$ \begin{cases} \{A\}^t = \langle A \rangle = \langle A_2, A_3 \rangle \\ \{B\}^t = \langle B \rangle = \langle B_2, B_3 \rangle \end{cases} \qquad [8] $$

et la condition $\theta_z(1,0,0)=1$ impose :

$$ b_2 B_2 + c_2 A_2 + b_3 B_3 + c_3 A_3 = -4A \qquad [9] $$

Enfin la rotation au nœud 1 est imposée, il s'ensuit que le glissement au nœud 1 est nul, soit :

$$ \gamma_{xy}(1,0,0) = u, v, w = 0 \qquad [10] $$

et nous obtenons :

$$ b_2 B_2 + c_2 A_2 + b_3 B_3 + c_3 A_3 = 0 \qquad [11] $$

Le système constitué des équations [9] et [11] a une infinité de solutions et est équivalent à :

$$ \begin{cases} b_2 B_2 + b_3 B_3 = 2A \\ c_2 A_2 + c_3 A_3 = -2A \end{cases} \qquad [12] $$

Alors, compte tenu de l'expression de $2A$ donnée par l'équation [4], nous choisissons :

$$ B_2 = c_3 \qquad B_3 = -c_2 \qquad A_2 = b_3 \qquad A_3 = -b_2 \qquad [13] $$

Finalement :

$$ u = \langle L_1^2 L_2, L_1^2 L_3 \rangle + \begin{Bmatrix} b_3 \\ -b_2 \end{Bmatrix} \quad \text{et} \quad v = \langle L_1^2 L_2, L_1^2 L_3 \rangle + \begin{Bmatrix} c_3 \\ -c_2 \end{Bmatrix} \qquad [14] $$

Ce champ fournit une rotation $\theta_{I_z}=1$ au nœud 1 et par permutation on obtient les champs $u$ et $v$ associés à des rotations unitaires aux nœuds 2 et 3. Donc les trois rotations nodales induisent un champ de déplacements donné par :
$$
\begin{align*}
u &= < L_1^2(L_2 b_3 - L_3 b_2) , L_2^2(L_3 b_1 - L_1 b_3) , L_3^2(L_4 b_2 - L_2 b_1) >
  \begin{Bmatrix} \theta_{1z} \\ \theta_{2z} \\ \theta_{3z} \end{Bmatrix} \\
v &= < L_1^2(L_2 c_3 - L_3 c_2) , L_2^2(L_3 c_1 - L_1 c_3) , L_3^2(L_4 c_2 - L_2 c_1) >
  \begin{Bmatrix} \theta_{1z} \\ \theta_{2z} \\ \theta_{3z} \end{Bmatrix}
\end{align*}
\tag{15}
$$

ou encore :

$$
\left\{ \begin{matrix} u \\ v \end{matrix} \right\} = \left[ \begin{matrix} N^{(1)u\theta} & & N^{(2)u\theta} & & N^{(3)u\theta} \\ & & & & \\ N^{(1)v\theta} & & N^{(2)v\theta} & & N^{(3)v\theta} \\ & & & & \\ u_1 & & u_2 & & u_3 \\ & & & & \\ v_1 & & v_2 & & v_3 \end{matrix} \right] \left\{ \begin{matrix} \theta_{1z} \\ \theta_{2z} \\ \theta_{3z} \end{matrix} \right\} = \left[ \begin{matrix} < N^{u\theta} > & & < \theta_{1z} \\ < N^{v\theta} > & & \theta_{2z} \\ & & < \theta_{3z} > \end{matrix} \right] \left\{ \begin{matrix} \theta_{1z} \\ \theta_{2z} \\ \theta_{3z} \end{matrix} \right\} \quad [16]
$$

2.2. Modification du triangle à 3 nœuds, modèle linéaire CST

Le champ des déplacements du triangle à 3 nœuds, (modèle linéaire) est défini
sous sa forme habituelle, soit en fonction des coordonnées barycentriques :

$$
u = < L_1, L_2, L_3 > \begin{Bmatrix} u_1 \\ u_2 \\ u_3 \end{Bmatrix} \quad \text{et} \quad v = < L_4, L_5, L_6 > \begin{Bmatrix} v_1 \\ v_2 \\ v_3 \end{Bmatrix} \tag{17}
$$

Ce champ de déplacements induit une rotation normale θ' constante sur
l'élément. Compte tenu de l'équation [7] elle est telle que :

$$
4A\theta' = < -c_1, -c_2, -c_3, b_1, b_2, b_3 > \begin{Bmatrix} u_1 \\ u_2 \\ u_3 \\ v_1 \\ v_2 \\ v_3 \end{Bmatrix} \quad [18]
$$

Finalement, en superposant le champ des déplacements linéaires de l’élément membranaire CST et le champ des déplacements associé à un champ des rotations normales θᵢz - θ' à chaque nœud i, nous obtenons :
$$
\begin{equation}
\begin{aligned}
\{u\} &= \langle L_1, L_2, L_3 \rangle \begin{Bmatrix} u_1 \\ u_2 \\ u_3 \end{Bmatrix} + \langle N^{u0} \rangle \begin{Bmatrix} \begin{Bmatrix} \theta_{1z} \\ \theta_{2z} \\ \theta_{3z} \end{Bmatrix} - \begin{Bmatrix} \theta' \\ \theta' \\ \theta' \end{Bmatrix} \end{aligned} \tag{19} \\
\{v\} &= \underbrace{\langle L_1, L_2, L_3 \rangle \begin{Bmatrix} v_1 \\ v_2 \\ v_3 \end{Bmatrix}}_{\text{Partie conforme}} + \langle N^{v0} \rangle \begin{Bmatrix} \begin{Bmatrix} \theta_{1z} \\ \theta_{2z} \\ \theta_{3z} \end{Bmatrix} - \begin{Bmatrix} \theta' \\ \theta' \\ \theta' \end{Bmatrix} \}
\end{equation}
$$

ou encore :

$$
\begin{equation}
\left\{ \begin{matrix} u \\ v \end{matrix} \right\} = \left[ \begin{matrix} L_1 + \frac{c_1}{4A} \sum_i N^{(i)u0} & L_2 + \frac{c_2}{4A} \sum_i N^{(i)v0} & L_3 + \frac{c_3}{4A} \sum_i N^{(i)v0} \\ \frac{c_1}{4A} \sum_i N^{(i)v0} & \frac{c_2}{4A} \sum_i N^{(i)u0} & \frac{c_3}{4A} \sum_i N^{(i)u0} \end{matrix} \right] \begin{Bmatrix} u_1 \\ u_2 \\ u_3 \end{Bmatrix} \\
+ \left[ \begin{matrix} -\frac{b_1}{4A} \sum_i N^{(i)u0} & -\frac{b_2}{4A} \sum_i N^{(i)v0} & -\frac{b_3}{4A} \sum_i N^{(i)v0} \\ L_1 - \frac{b_1}{4A} \sum_i N^{(i)v0} & L_2 - \frac{b_2}{4A} \sum_i N^{(i)u0} & L_3 - \frac{b_3}{4A} \sum_i N^{(i)u0} \end{matrix} \right] \begin{Bmatrix} v_1 \\ v_2 \\ v_3 \end{Bmatrix} \tag{20} \\
+ \begin{Bmatrix} \langle N^{u0} \rangle \\ \langle N^{v0} \rangle \end{Bmatrix} \begin{Bmatrix} \theta_{1z} \\ \theta_{2z} \\ \theta_{3z} \end{Bmatrix}
\end{equation}
$$

A la rotation $\theta_{iz}-\theta'$ au nœud $i$, due aux termes non classiques, s'ajoute la rotation $\theta'$ due au champ linéaire, et finalement la rotation au nœud $i$ est $\theta_{iz}$.

En regroupant les variables nodales par nœud :

$$
\{u_{nodal}\} = \{u_1, v_1, \theta_{1z}, u_2, v_2, \theta_{2z}, u_3, v_3, \theta_{3z}\} \quad [21]
$$

le champ des déplacements peut être mis sous la forme :

$$
\left\{ \begin{matrix} u \\ v \end{matrix} \right\} = \left[ \begin{matrix} \langle \Psi \rangle & 0 \\ 0 & \langle \Psi \rangle \end{matrix} \right] \left[ \begin{matrix} C_u \\ C_v \end{matrix} \right] \{ u_{nodal} \} \quad [22]
$$

avec :

$$
\langle \Psi \rangle = \langle L_1, L_1^2 L_2, L_2^2 L_3, L_3^2 L_4, L_4^2 L_5, L_5^2 L_6 \rangle \quad [23]
$$

et :
$$ [C_u] = \begin{bmatrix} C_u^{11} & C_u^{12} & C_u^{13} \\ C_u^{21} & C_u^{22} & C_u^{23} \\ C_u^{31} & C_u^{32} & C_u^{33} \end{bmatrix} \qquad [C_v] = \begin{bmatrix} C_v^{11} & C_v^{12} & C_v^{13} \\ C_v^{21} & C_v^{22} & C_v^{23} \\ C_v^{31} & C_v^{32} & C_v^{33} \end{bmatrix} \qquad [24] $$

où :

$$ [C_u^{11}] = \begin{bmatrix} 1 & 0 & 0 \\ \frac{c_1 b_3}{4A} & \frac{-b_1 b_3}{4A} & b_3 \\ -\frac{c_1 b_2}{4A} & \frac{b_1 b_2}{4A} & -b_2 \end{bmatrix}, \quad [C_v^{11}] = \begin{bmatrix} 0 & 1 & 0 \\ \frac{c_1 c_3}{4A} & \frac{-b_1 c_3}{4A} & c_3 \\ -\frac{c_1 c_2}{4A} & \frac{b_1 c_2}{4A} & -c_2 \end{bmatrix} $$

$$ [C_u^{12}] = \begin{bmatrix} 0 & 0 & 0 \\ \frac{c_2 b_3}{4A} & \frac{-b_2 b_3}{4A} & 0 \\ -\frac{c_2 b_2}{4A} & \frac{b_2 b_2}{4A} & 0 \end{bmatrix}, \quad [C_v^{12}] = \begin{bmatrix} 0 & 0 & 0 \\ \frac{c_2 c_3}{4A} & \frac{-b_2 c_3}{4A} & 0 \\ -\frac{c_2 c_2}{4A} & \frac{b_2 c_2}{4A} & 0 \end{bmatrix} $$

$$ [C_u^{13}] = \begin{bmatrix} 0 & 0 & 0 \\ \frac{c_3 b_3}{4A} & \frac{-b_3 b_3}{4A} & 0 \\ -\frac{c_3 b_2}{4A} & \frac{b_3 b_2}{4A} & 0 \end{bmatrix}, \quad [C_v^{13}] = \begin{bmatrix} 0 & 0 & 0 \\ \frac{c_3 c_3}{4A} & \frac{-b_3 c_3}{4A} & 0 \\ -\frac{c_3 c_2}{4A} & \frac{b_3 c_2}{4A} & 0 \end{bmatrix} \qquad [25] $$

Les autres matrices $[C_u^{ij}]$ et $[C_v^{ij}]$ sont obtenues par permutation circulaire sur les indices $i$, avec $j=1, 2, 3$.

## 2.3. Matrice de rigidité

Les opérateurs de dérivation sont donnés par :

$$ \frac{\partial}{\partial x} = \frac{1}{2A} \sum_{i=1}^{3} b_i \frac{\partial}{\partial L_i} \quad \text{et} \quad \frac{\partial}{\partial y} = \frac{1}{2A} \sum_{i=1}^{3} c_i \frac{\partial}{\partial L_i} \qquad [26] $$
En notant $<X> = <L_1 L_2, L_2 L_3, L_3 L_4>$ et $\{u_{nodal}^i\}$ le vecteur des variables nodales du nœud $i$ ($<u_{nodal}^i> = <u_i, v_i, \theta_{iz}>$) les dérivées des composantes du champ des déplacements sont données par :

$$
\begin{align*}
u_x &= \frac{1}{2A} \sum_{i=1}^{3} <X > b_i [D^i] [C_u^i] \{u_{nodal}^i\} \\
v_x &= \frac{1}{2A} \sum_{i=1}^{3} <X > b_i [D^i] [C_v^i] \{u_{nodal}^i\} \\
u_y &= \frac{1}{2A} \sum_{i=1}^{3} <X > c_i [D^i] [C_u^i] \{u_{nodal}^i\} \\
v_y &= \frac{1}{2A} \sum_{i=1}^{3} <X > c_i [D^i] [C_v^i] \{u_{nodal}^i\}
\end{align*}
\tag{27}
$$

où les matrices $[D^i]$ sont égales à :

$$
[D^1] = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 2 \end{bmatrix} \quad [D^2] = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 2 \\ 0 & 2 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad [D^3] = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 2 \\ 0 & 2 & 0 \end{bmatrix} \qquad [28]
$$

Dans ces expressions, compte tenu des matrices $[C_u^{ij}]$ et $[C_v^{ij}]$, les termes en $L_i^2$ dus à la dérivation disparaissent.

On obtient alors :

$$
\{\epsilon\} = \left\{ \begin{matrix} u,_{x} \\ v,_{y} \\ u,_{y} + v,_{x} \end{matrix} \right\} = \frac{1}{2A} \langle X \rangle \begin{bmatrix} [\alpha] & [0] \\ [0] & [\beta] \\ [\beta] & [\alpha] \end{bmatrix} [C_u] \{u_{nodal}\} \quad [29]
$$

avec :

$$
[\alpha] = [b_1[D^1], b_2[D^2], b_3[D^3]] \qquad [30]
$$

$$
[\beta] = [c_1[D^1], c_2[D^2], c_3[D^3]] \qquad [31]
$$
On considère un milieu homogène isotrope en état de contraintes planes. Il est caractérisé par une matrice d'élasticité. $[D] = \frac{E}{1-v^2} \begin{bmatrix} 1 & v & 0 \\ v & 1 & 0 \\ 0 & 0 & \frac{1-v}{2} \end{bmatrix}$

Le calcul de l'énergie de déformation $2E_d = h \int_S \{\epsilon\}'[D]\epsilon\}dS$ conduit à la matrice de rigidité $[K]$, telle que : $2E_d = \{u_{nodal}\}'[K]\{u_{nodal}\}$.

On obtient alors :

$$[K] = \frac{Eh}{2A(1-v^2)} [C]' \begin{bmatrix} [\alpha]' & [0] & [\beta]' \\ [0] & [\beta]' & [\alpha]' \end{bmatrix} \times \int_S \begin{bmatrix} [X]' \\ [X]' \\ [X]' \end{bmatrix} [D] \begin{bmatrix} [X] \\ [X] \\ [X] \end{bmatrix} \begin{bmatrix} [\alpha] & [0] \\ [0] & [\beta] \\ [\beta] & [\alpha] \end{bmatrix} dS[C] \quad [32]$$

soit :

$$[K] = \frac{Eh}{2A(1-v^2)} [C]'[E][C] \quad [33]$$

où la matrice $[E]$ est donnée par :

$$[E] = \begin{bmatrix} [\alpha]'[R][\alpha] + [\beta]'[R][\beta] \left(\frac{1-v}{2}\right) & v[\alpha]'[R][\beta] + [\beta]'[R][\alpha] \left(\frac{1-v}{2}\right) \\ \text{Sym.} & [\beta]'[R][\beta] + [\alpha]'[R][\alpha] \left(\frac{1-v}{2}\right) \end{bmatrix} \quad [34]$$

et la matrice $[R]$ par :

$$2A[R] = \int_S \{X\}' \langle X \rangle dS = 2A \begin{bmatrix} \frac{1}{2} & \frac{1}{4!} & \frac{1}{4!} & \frac{1}{4!} \\ & \frac{4}{6!} & \frac{2!}{6!} & \frac{2!}{6!} \\ & & \frac{4}{6!} & \frac{2!}{6!} \\ \text{Sym.} & & & \frac{4}{6!} \end{bmatrix} \quad [35]$$
## 2.4. Matrice masse

Le calcul de l'énergie cinétique $2E_c = \{\dot{u}_{nodal}\}^T[M]\{\dot{u}_{nodal}\}$, conduit à la matrice masse élémentaire :

$$[M] = \left[ [C_u]'_S \{\Psi\} \langle \Psi \rangle dS[C_u] \right] + \left[ [C_v]'_S \{\Psi\} \langle \Psi \rangle dS[C_v] \right] \quad [36]$$

$$\text{où :} \quad \int_S \{\Psi\} \langle \Psi \rangle dS = \begin{bmatrix} [M_1] & [M_2] & [M_2]' \\ [M_2]' & [M_1] & [M_2] \\ [M_2] & [M_2]' & [M_1] \end{bmatrix} \quad [37]$$

avec :

$$[M_1] = 2A \begin{bmatrix} \frac{2!}{4!} & \frac{3!}{6!} & \frac{3!}{6!} \\ \frac{2!4!}{8!} & \frac{4!}{8!} & \frac{2!4!}{8!} \\ \text{Sym.} & & \end{bmatrix} \quad \text{et} \quad [M_2] = 2A \begin{bmatrix} \frac{1}{4!} & \frac{2}{6!} & \frac{2!2!}{6!} \\ \frac{2!2!}{6!} & \frac{2!3!}{8!} & \frac{2!3!}{8!} \\ \frac{2}{6!} & \frac{2!2!2!}{8!} & \frac{2!3!}{8!} \end{bmatrix} \quad [38]$$

## 2.5. Vecteur des forces généralisées associées à une répartition de tension linéique uniforme sur un côté

On considère le triangle 1-2-3 et on applique une tension uniforme d'intensité $\vec{p}_0$ sur le coté (1-2) (figure 2a).

Le travail virtuel de cette répartition linéique est égal à :

$$\delta\tau = \int_{(1-2)} \delta\vec{u} \cdot \vec{p} dl = \int_{(1-2)} \langle \delta u \rangle|_{(1-2)} \{p\} dl \quad [39]$$

où :

$$\{u\}_{(1-2)} = \begin{bmatrix} \langle \Psi \rangle|_{(1-2)} & 0 \\ 0 & \langle \Psi \rangle|_{(1-2)} \end{bmatrix} \begin{bmatrix} [C_u] \\ [C_v] \end{bmatrix} \{u_{nodal}\} \quad [40]$$

Donc il vient :
$$ \delta\tau = \langle \delta u_{nodal} \rangle = \int_{(1-2)} \left[ [C_u]' \cdot [C_v]' \right] \left[ \begin{matrix} \Psi \end{matrix} \right]_{(1-2)} \begin{bmatrix} 0 \\ \Psi \end{bmatrix}_{(1-2)} \{p\} dl \quad [41] $$

<figure><img src="image_3.png" /><figcaption>Figure 2a. Répartition de tension linéique suivant le coté (1-2)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2b. Forces généralisées correspondant à la répartition linéique suivant le coté (1-2)</figcaption></figure>

Si on considère le cas particulier d'une répartition $\{p\} = \begin{Bmatrix} p_0 \\ 0 \end{Bmatrix}$, le vecteur des forces généralisées correspondant est égal à :

$$ \{\varphi\} = p_0 [C_v]' \int_{(1-2)} \{\Psi\} \}_{(1-2)} dl \quad [42] $$
Le coté (1-2) correspondant à $L_3 = 0$ permet d'écrire :

$$ \int_{(1-2)} <\Psi|_{(1-2)} dl = \int_{(1-2)} <L_1, L_1^2, L_2, 0, L_2, 0, L_1, L_2^2, 0, 0, 0> dl \quad [43] $$

Le vecteur des forces généralisées est (figure 2b) :

$$ \{\phi\} = <0, p_0 \frac{l_{12}}{2}, p_0 \frac{l_{12}^2}{12}, 0, p_0 \frac{l_{12}}{2}, -p_0 \frac{l_{12}^2}{12}, 0, 0, 0> \quad [44] $$

## 2.6. Validation

Pour valider l'élément membranaireT3D avec degrés de liberté rotationnels, nous montrons tout d'abord qu'il passe le patch test d'ordre zéro, donc qu'il converge. Puis nous réalisons plusieurs calculs sur des tests classiques de la littérature qui concernent des patchs tests d'ordre supérieur et l'étude du panneau fuselé de Cook [COO 74], ainsi qu'une étude de convergence sur une plaque triangulaire.

### 2.6.1. Patch test « dur » d'ordre zéro

<figure><img src="image_5.png" /><figcaption>**Figure 3. Notations**</figcaption></figure>

Pour que l'élément proposé passe le patch test « dur » d'ordre zéro [ZIEN 91], l'intégrale de contour $\int_{1-2-3} \vec{T} \cdot \delta\vec{u} dl$ doit être nulle pour un champ de tension $\vec{T}$ uniforme sur le triangle. Les composantes $u$ et $v$ jouant le même rôle, il suffit de calculer $\int_{1-2-3} \delta u dl$. Avec les notations $dl$ et $\alpha$ qui sont définies sur la figure 3, sur chaque coté on a $dl \cos \alpha = dx$.

Sur le côté (2-3), on a d'après l'équation [20] :
$L_1 = 0$

$$u = L_2 u_2 + L_3 u_3$$

$$+ L_2^2 L_3 b_1 \left( \theta_{2z} + \frac{1}{4A} (c_1 u_1 + c_2 u_2 + c_3 u_3 - b_1 v_1 - b_2 v_2 - b_3 v_3) \right) \quad [45] \\ - L_3^2 L_2 b_1 \left( \theta_{3z} + \frac{1}{4A} (c_1 u_1 + c_2 u_2 + c_3 u_3 - b_1 v_1 - b_2 v_2 - b_3 v_3) \right)$$

La partie non conforme pour la contribution de ce côté provient des termes :

$$(L_2^2 L_3 b_1 - L_3^2 L_2 b_1)(c_1 u_1 - b_1 v_1) \quad [46]$$

Ainsi :

$$\begin{aligned} \delta\tau = \frac{1}{4A} & \left( \int_{(2-3)}^{1} b_1 L_2 L_3 (L_2 - L_3) \delta(c_1 u_1 - b_1 v_1) dl \right. \\ & \left. + \int_{(3-1)}^{1} b_2 L_3 L_4 (L_3 - L_4) \delta(c_2 u_2 - b_2 v_2) dl \right. \\ & \left. + \int_{(1-2)}^{1} b_3 L_4 L_2 (L_4 - L_2) \delta(c_3 u_3 - b_3 v_3) dl \right) \end{aligned} \quad [47]$$

En fait, il suffit d'évaluer : $\int_0^1 L_2 L_3 (L_2 - L_3) dx$ à $L_1 = 0$ avec $x = (x_2 - x_3)L_2 + x_3$.

Soit $\int_0^1 L_2(1-L_2)(2L_2-1)(x_2-x_3)dL_2 = \int_0^1 (3L_2^2 - 2L_2^3 - L_2)(x_2 - x_3)dL_2 = 0$

Par permutation : $\int_{(3-1)} L_3 L_1 (L_3 - L_1) dx = 0$ et $\int_{(1-2)} L_4 L_2 (L_4 - L_2) dx = 0$.

Ainsi $\delta\tau = 0$.

L'élément passe le patch test « dur » d'ordre zéro, donc il converge.

### 2.6.1. Patch test d'ordre supérieur

#### 2.6.1.1. Poutre cantilever soumise à un effort de cisaillement

Pour compléter l'étude de convergence, nous reprenons le test proposé initialement par MacNeal et Harder [MAC 84, MAC 88]. Ce test concerne l'étude d'une poutre cantilever, de longueur $l=48$ mm, de hauteur $h=12$ mm, d'épaisseur unité. Elle est constituée d'un matériau de caractéristiques $E=30\,000$ N/mm² et $\nu=0.25$. Elle est soumise à un effort de cisaillement $P=40$ N appliqué à l'extrémité libre $x=l$ dans la direction $\vec{e}_2$ (figure 4).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Poutre cantilever-Géométrie et chargement</figcaption></figure>

On peut trouver les résultats de ce test dans les articles de nombreux auteurs et en particulier dans ceux de Allman [ALL 84], Bergan et Felippa [BER 85], Hughes et Brezzi [HUG 89], Yunus et al. [YUN 89], Ibrahimbegovic et al. [IBR 90] par exemple. Pour ce problème, à l'extrémité libre, le déplacement vertical théorique du point A situé sur l'axe longitudinal de la poutre est égal à :

$$v_A = \frac{Pl^3}{3EI} + \frac{(4 + 5v)Pl}{2Eh} \quad [48]$$

La solution analytique du problème en élasticité est donnée, par exemple, dans Timoshenko et Goodier [TIM 51].

Dans le tableau 1 suivant, nous présentons, pour différentes finesses de maillages réguliers (figure 5), les résultats donnés par ces différents auteurs et ceux obtenus avec l'élément T3D proposé.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Poutre cantilever. Maillages réguliers</figcaption></figure>
<table><caption>**Tableau 1. Déplacement vertical $v_A$ (en mm)-Maillages réguliers**</caption><thead><tr><th rowspan="2">Type d'élément</th><th colspan="3">Maillages</th></tr><tr><th>4x1</th><th>8x2</th><th>16x4</th></tr></thead><tbody><tr><td>CST</td><td>0.0903</td><td>0.198</td><td>0.295</td></tr><tr><td>T3D (présente théorie)</td><td>0.235</td><td>0.315</td><td>0.344</td></tr><tr><td>Allman [ALL 84]</td><td>0.2696</td><td>0.3261</td><td>0.3471</td></tr><tr><td>Olsen et Bearden [OLS 79]</td><td>0.2302</td><td>0.2983</td><td>0.3291</td></tr><tr><td>Ibrahimbegovic et al [IBR 90]</td><td>0.3445</td><td>0.3504</td><td>0.3543</td></tr><tr><td>Yunus et al (élém. triangul.) [YUN 89]</td><td>0.3609</td><td></td><td></td></tr><tr><td>Yunus et al (élém. Quadr.) [YUN 89]</td><td>0.3312</td><td></td><td></td></tr><tr><td>Cook (modèle E5) [COO 91]</td><td>0.349</td><td></td><td></td></tr><tr><td>Théorie [TIM 51]</td><td></td><td>0.35583</td><td></td></tr></tbody></table>

Dans le tableau 2, nous présentons les résultats obtenus pour différentes finesses de maillages irréguliers (figure 6).

<figure><img src="image_8.png" /><figcaption>Maillage 4x1</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Maillage 8x2</figcaption></figure>

<caption>**Figure 6. Poutre cantilever. Maillages irréguliers**</caption>
<table><caption>**Tableau 2. Déplacement vertical $v_A$ (en mm)-Maillages irréguliers**</caption><thead><tr><th rowspan="2">Type d'élément</th><th colspan="3">Maillages</th></tr><tr><th>4x1</th><th>8x2</th><th>16x4</th></tr></thead><tbody><tr><td>CST</td><td>0.088</td><td>0.186</td><td>0.284</td></tr><tr><td>T3D (présente théorie)</td><td>0.201</td><td>0.299</td><td>0.338</td></tr><tr><td>Ibrahimbegovic et al [IBR 90]</td><td>0.3065</td><td></td><td></td></tr><tr><td>Yunus et al (élém. triangul.) [YUN 89]</td><td>0.3509</td><td></td><td></td></tr><tr><td>Yunus et al (élém. Quadr.) [YUN 89]</td><td>0.3442</td><td></td><td></td></tr><tr><td>Allman (d'après [COO 91])</td><td>0.246</td><td></td><td></td></tr><tr><td>Cook (modèle E5) [COO 91]</td><td>0.313</td><td></td><td></td></tr></tbody></table>

### 2.6.1.2. Poutre cantilever soumise à un couple

Pour compléter cette étude de patch test, la poutre cantilever est maintenant soumise à un couple à l'extrémité libre. Les résultats donnés dans les tableaux 3 et 4, sont obtenus respectivement pour les maillages réguliers (figure 5) et irréguliers (figure 6) lorsque le couple est appliqué par l'intermédiaire de deux forces (figure 7-a). Les déplacements verticaux du point A, donnés par l'élément T3D sont comparés à ceux obtenus par Allman [ALL 84] et à ceux de Cook [COO 91].

<figure><img src="image_10.png" /><figcaption>**Figure 7. Poutre cantilever soumise à un couple**</figcaption></figure>

La valeur du couple appliqué est de 6000 Nmm.

Le tableau 5 donne les résultats obtenus avec l'élément T3D, lorsque le couple appliqué en A est un moment concentré (figure 7b).
<table><caption>**Tableau 3.** Déplacement vertical $v_A$ (en mm)-Maillages réguliers-Couple à l'extrémité</caption><thead><tr><th rowspan="2"></th><th colspan="3">Maillages réguliers</th></tr><tr><th>4x1</th><th>8x2</th><th>16x4</th></tr></thead><tbody><tr><th>CST</th><td>0.37</td><td>0.86</td><td>1.32</td></tr><tr><th>T3D (présente théorie)</th><td>1.02</td><td>1.44</td><td>1.58</td></tr><tr><th>Allman [ALL 84]]</th><td>1.23</td><td></td><td></td></tr><tr><th>Cook (modèle E5) [COO 91]</th><td>1.6</td><td></td><td></td></tr><tr><th>Théorie des poutres</th><td colspan="3">v<sub>A</sub>=1.6</td></tr></tbody></table>

<table><caption>**Tableau 4.** Déplacement vertical $v_A$ (en mm)-Maillages irréguliers-Couple à l'extrémité</caption><thead><tr><th rowspan="2"></th><th colspan="3">Maillages irréguliers</th></tr><tr><th>4x1</th><th>8x2</th><th>16x4</th></tr></thead><tbody><tr><th>CST</th><td>0.38</td><td>0.79</td><td>1.32</td></tr><tr><th>T3D (présente théorie)</th><td>0.91</td><td>1.31</td><td>1.54</td></tr><tr><th>Allman [ALL 84]</th><td>1.23</td><td></td><td></td></tr><tr><th>Cook (modèle E5) [COO 91]</th><td>1.6</td><td></td><td></td></tr><tr><th>Théorie des poutres</th><td colspan="3">v<sub>A</sub>=1.6</td></tr></tbody></table>

<table><caption>**Tableau 5.** Déplacement vertical $v_A$ (en mm) -Moment à l'extrémité</caption><thead><tr><th rowspan="2"></th><th colspan="3">Maillages</th></tr><tr><th>4x1</th><th>8x2</th><th>16x4</th></tr></thead><tbody><tr><th>T3D (Maillages réguliers)</th><td>1.08</td><td>1.41</td><td>1.56</td></tr><tr><th>T3D (Maillages irréguliers)</th><td>1.08</td><td>1.41</td><td>1.57</td></tr></tbody></table>

### 2.6.2. Plaque triangulaire plane en statique-Etude de convergence

Dans cette section, pour poursuivre la caractérisation de l'élément triangulaire T3D proposé, on veut par comparaison à des résultats obtenus avec l'élément ordinaire CST, évaluer la vitesse de convergence et montrer l'influence des
conditions aux limites sur les résultats d'un calcul. Pour ce faire on réalise l'étude élémentaire d'une plaque triangulaire. La plaque triangulaire rectangle de dimensions L=2 m et l=1 m, est discrétisée en triangles homothétiques à la plaque initiale (figure 8). Elle est encastrée suivant le côté BC de longueur L=2 m et soumise à une tension uniforme p = 4 N/m suivant le côté AB de longueur l=1 m.

L'épaisseur de la plaque est h=1 cm et les caractéristiques du matériau sont E=200 000 MPa et ν=0.3.

On étudie la convergence des résultats en fonction du nombre total N de triangles construits. Si n est le nombre de subdivisions des côtés de la plaque, le nombre d'éléments finis triangulaires N construits est $N = n^2$.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Géométrie et discrétisation de la plaque triangulaire</figcaption></figure>

Le tableau 6 donne les déplacements (en 10⁻⁸m) horizontal et vertical du point A, obtenus pour une modélisation ordinaire en éléments CST et ceux obtenus avec une modélisation en éléments T3D.

Deux types de conditions aux limites sont prises en compte :

- rotation bloquée le long de l'appui (résultats notés *Drill*\*). Ce type de conditions aux limites, correspond le mieux aux calculs faits avec le CST,

- rotation libre le long de l'appui (résultats notés *Drill*\*\*).

On vérifie la bonne convergence de l'élément T3D avec les degrés de liberté rotationnels. Cet élément est plus souple que l'élément CST. De plus la rotation normale comme degré de liberté permet de prendre en compte des conditions aux limites plus réalistes au niveau des appuis.

Les courbes des figures 9, 10 et 11 montrent les résultats de la convergence et la comparaison entre les deux types de conditions aux limites.
<table><thead><tr><th>N</th><th></th><th>CST</th><th>Drill*</th><th>Drill**</th></tr></thead><tbody><tr><td rowspan="2">4</td><td>U<sub>A</sub></td><td>0.054</td><td>0.239</td><td>0.275</td></tr><tr><td>V<sub>A</sub></td><td>0.371</td><td>0.653</td><td>0.696</td></tr><tr><td rowspan="2">9</td><td>U<sub>A</sub></td><td>0.108</td><td>0.262</td><td>0.314</td></tr><tr><td>V<sub>A</sub></td><td>0.452</td><td>0.681</td><td>0.744</td></tr><tr><td rowspan="2">16</td><td>U<sub>A</sub></td><td>0.152</td><td>0.287</td><td>0.335</td></tr><tr><td>V<sub>A</sub></td><td>0.518</td><td>0.715</td><td>0.777</td></tr><tr><td rowspan="2">25</td><td>U<sub>A</sub></td><td>0.188</td><td>0.305</td><td>0.344</td></tr><tr><td>V<sub>A</sub></td><td>0.571</td><td>0.742</td><td>0.796</td></tr><tr><td rowspan="2">81</td><td>U<sub>A</sub></td><td>0.264</td><td>0.336</td><td>0.353</td></tr><tr><td>V<sub>A</sub></td><td>0.69</td><td>0.793</td><td>0.820</td></tr><tr><td rowspan="2">196</td><td>U<sub>A</sub></td><td>0.302</td><td>0.349</td><td>0.357</td></tr><tr><td>V<sub>A</sub></td><td>0.748</td><td>0.814</td><td>0.829</td></tr><tr><td rowspan="2">361</td><td>U<sub>A</sub></td><td>0.319</td><td>0.359</td><td>0.354</td></tr><tr><td>V<sub>A</sub></td><td>0.775</td><td>0.883</td><td>0.833</td></tr></tbody></table>

*Drill*: la rotation $\theta_z$ est bloquée le long de l'appui.

**Drill**: la rotation $\theta_z$ est libre le long de l'appui.

<figure><img src="image_12.png" /><figcaption>**Figure 9. Déplacement horizontal du point A**</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Déplacement vertical du point A</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Influence des conditions limites sur le déplacement vertical du point A</figcaption></figure>

### 2.6.3. Etude du panneau fuselé de Cook

La validation de l'élément T3D est complétée par l'étude du panneau fuselé de Cook. Ce test a été proposé initialement par Cook [COO 74], puis ensuite par Bergan et Felippa [BER 85] (figure 12). Ce problème est un cas test standard qui a été étudié par de nombreux auteurs et dans le tableau 7 suivant nous donnons les résultats obtenus avec l'élément T3D ainsi que quelques résultats de la littérature.

Le panneau est encastré dans la section $x=0$, il est soumis à un effort vertical $P=1$ N dans la section $x=48$. Les caractéristiques du matériau sont $E=10\ 000$ MPa, $v=0.25$ et l'épaisseur du panneau est $e=1$ mm.

La solution de référence est celle donnée par Bergan et Felippa [BER 85], elle correspond aux résultats obtenus avec un maillage 100 x 100.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Panneau fuselé de Cook-géométrie (dimensions en mm) -maillage</figcaption></figure>

<table><caption>**Tableau 7.** Panneau fuselé de Cook-déplacement vertical du point A</caption><thead><tr><th>Maillage</th><th>1x1</th><th>2x2</th><th>4x4</th><th>5x5</th><th>8x8</th><th>16x16</th><th>32x32</th></tr></thead><tbody><tr><th>CST</th><td>5.97</td><td>12.09</td><td>18.35</td><td>19.92</td><td>22.34</td><td>23.41</td><td>23.44</td></tr><tr><th>T3D</th><td>12.36</td><td>17.49</td><td>21.53</td><td>22.10</td><td>23.098</td><td>23.70</td><td>23.88</td></tr><tr><th>Allman</th><td></td><td>21.52</td><td>23.17</td><td></td><td></td><td></td><td></td></tr><tr><th>Bergan & Filippa</th><td></td><td>20.36</td><td>22.42</td><td></td><td></td><td></td><td></td></tr><tr><th>Ibrahimbegovic</th><td>14.066</td><td>20.68</td><td>22.98</td><td>23.63</td><td></td><td></td><td></td></tr><tr><th>Yunus et al (AT)</th><td></td><td>20.93</td><td>22.94</td><td></td><td></td><td></td><td></td></tr><tr><th>Yunus et al (AQ)</th><td></td><td>21.85</td><td>23.26</td><td></td><td></td><td></td><td></td></tr><tr><th>Ayad [AYA 95a]</th><td></td><td>21.10</td><td>22.88</td><td></td><td></td><td></td><td></td></tr><tr><th>Cook (E1)</th><td></td><td>18.89</td><td>22.47</td><td></td><td></td><td></td><td></td></tr><tr><th>Cook (E5)</th><td></td><td>21.76</td><td>23.19</td><td></td><td></td><td></td><td></td></tr><tr><th>Solution de référ.</th><td></td><td></td><td></td><td>23.91</td><td></td><td></td><td></td></tr></tbody></table>

### 2.6.4. Conclusion

Pour les cas tests étudiés, les différents résultats obtenus, montrent que l'élément T3D proposé, triangulaire membranaire avec degrés de liberté rotationnels, a un comportement très satisfaisant. Même si le taux de convergence est un peu moins
bon que celui de certains éléments proposés dans la littérature, il reste cependant intéressant puisqu'il assure la convergence vers la solution exacte (il passe le patch test d'ordre zéro), l'intégration analytique la matrice rigidité est exacte et l'élément ne présente pas de mode de déformation parasite. Des tests non présentés ici, confirment ce résultat. Enfin la formulation en déplacement et l'interpolation rigoureuse proposée ne nécessite pas la mise en œuvre d'une formulation mixte ou l'introduction des coefficients de forme. On dispose donc d'un élément robuste et simple d'utilisation.

### 3. Elément triangulaire pour la flexion

Pour utiliser l'élément T3D avec degrés de liberté rotationnels pour la modélisation de coques, on lui associe un élément de flexion. Les éléments finis triangulaires en flexion proposés dans la littérature sont nombreux. Pour une revue détaillée des éléments disponibles et de l'étude de leurs performances, on pourra par exemple se référer aux travaux de Batoz et Dhatt [BAT 88, BAT 90 et BAT 92].

Pour la validation envisagée, notre choix s'est porté sur l'élément triangulaire de Specht [SPE 88], repris par Zienckiewicz [ZIE 91]. Nous rappelons ci-après les caractéristiques essentielles de cet élément.

On utilise les mêmes notations que précédemment, et on introduit les nouvelles quantités :

$$ \mu_1 = \frac{l_3^2 - l_2^2}{l_1^2}, \mu_2 = \frac{l_1^2 - l_3^2}{l_2^2}, \mu_3 = \frac{l_2^2 - l_1^2}{l_3^2} \quad [49] $$

où $l_1, l_2, l_3$ sont respectivement les longueurs des cotés (2-3), (3-1), (1-2).

#### 3.1. Interpolation

Les degrés de liberté de l'élément sont $w_i, \theta_{xi}, \theta_{yi}$ à chaque nœud $i=1, 2, 3$ du triangle. Le déplacement transverse est donné par l'approximation polynomiale :

$$ w = \sum_{i=1}^{9} P_i \alpha_i = \langle P \rangle \{\alpha\} \quad [50] $$

avec $\langle P \rangle$ donné en fonction des coordonnées barycentriques par :
$$<P> = <L_1, L_2, L_3, L_1 L_2, L_2 L_3, L_3 L_1,$$

$$L_1^2 L_2 + \frac{1}{2} L_1 L_2 L_3 \{3(1 - \mu_3)L_1 - (1 + 3\mu_3)L_2 + (1 + 3\mu_3)L_3\}, \quad [51]$$

$$L_2^2 L_3 + \frac{1}{2} L_1 L_2 L_3 \{3(1 - \mu_1)L_2 - (1 + 3\mu_1)L_3 + (1 + 3\mu_1)L_4\},$$

$$L_3^2 L_4 + \frac{1}{2} L_1 L_2 L_3 \{3(1 - \mu_2)L_3 - (1 + 3\mu_2)L_4 + (1 + 3\mu_2)L_5\} >$$

## Matrice de connectivité

Les paramètres généralisés {$\alpha$} sont exprimés en fonction des degrés de liberté via la matrice de connectivité. On a d'une part :

$$\alpha_i = w_i \quad i = 1, 2, 3$$

D'autre part compte tenu de l'hypothèse de Bernoulli, on a :

$$\theta_x = w_y \quad \text{et} \quad \theta_y = -w_x \quad [52]$$

Et on déduit :

$$\theta_{x1} = \frac{1}{2A} [c_1, c_2, c_3, c_2, 0, c_3, c_2, 0, 0] \{\alpha\}$$

$$\theta_{y1} = \frac{1}{2A} [-b_1, -b_2, -b_3, -b_2, 0, -b_3, -b_2, 0, 0] \{\alpha\}$$

$$\theta_{x2} = \frac{1}{2A} [c_1, c_2, c_3, c_1, c_3, 0, 0, c_3, 0] \{\alpha\}$$

[53]

$$\theta_{y2} = \frac{1}{2A} [-b_1, -b_2, -b_3, -b_1, -b_3, 0, 0, -b_3, 0] \{\alpha\}$$

$$\theta_{x3} = \frac{1}{2A} [c_1, c_2, c_3, c_2, c_1, 0, 0, c_1] \{\alpha\}$$

$$\theta_{y3} = \frac{1}{2A} [-b_1, -b_2, -b_3, 0, -b_2, -b_1, 0, 0, -b_1] \{\alpha\}$$

Finalement la matrice de connectivité $[T]$ telle que :

$$\{\alpha\} = [T]\{u_{nodal}\} \quad [54]$$

est donnée par :
$$
\left\{
\begin{array}{c}
\alpha_1 \\
\alpha_2 \\
\alpha_3 \\
\alpha_4 \\
\alpha_5 \\
\alpha_6 \\
\alpha_7 \\
\alpha_8 \\
\alpha_9
\end{array}
\right\}
=
\begin{bmatrix}
1 & & & & & & & \\
& & & & & & & \\
& & 1 & & & & & \\
& & & & & 1 & & \\
-1 & & & & & & & \\
& & & -1 & & & & \\
1 & b_2 & c_2 & & & & & \\
2 & -b_3 & -c_3 & -2 & -b_3 & -c_3 & & \\
& & & 2 & -b_1 & -c_1 & -2 & -b_1 - c_1 \\
-2 & -b_2 & -c_2 & & & & 2 & -b_2 - c_2
\end{bmatrix}
\begin{bmatrix}
w_1 \\
\theta_{x1} \\
\theta_{y1} \\
w_2 \\
\theta_{x2} \\
\theta_{y2} \\
w_3 \\
\theta_{x3} \\
\theta_{y3}
\end{bmatrix}
\quad [55]
$$

et l'interpolation du déplacement transverse est :

$$
w = \langle P \rangle [T] \{u_{nodal}\} = \langle N \rangle \{u_{nodal}\} \quad [56]
$$

ce qui conduit à la matrice d'interpolation $\langle N \rangle = \langle\langle N_1 \rangle\rangle \langle N_2 \rangle \langle N_3 \rangle \rangle$.

On obtient pour $\langle N_1 \rangle$ :

$$
\begin{align*}
\langle N_1 \rangle = & \langle P_1 - P_4 + P_6 + 2(P_7 - P_9), \\
& b_2(P_6 - P_9) - b_3 P_7, \\
& a_2(P_6 - P_9) - a_3 P \rangle \tag{57}
\end{align*}
$$

D'une manière plus générale, on écrit pour {i, j, k} les permutations circulaires de {1, 2, 3}, et on obtient pour i=1, 2, 3:

$$
\langle N_i \rangle = \langle P_i - P_{i+3} + P_{k+3} + 2(P_{i+6} - P_{k+6}), \\
\qquad -b_j(P_{k+6} - P_{k+3}) - b_k P_{i+6}, \\
\qquad -a_j(P_{k+6} - P_{k+3}) - a_k P_{i+6} \rangle \tag{58}
$$

3.2. Matrices élémentaires

3.2.1. Matrice de rigidité

De façon ordinaire, le vecteur des déformations généralisées est :

$$
\{\epsilon\} = \begin{bmatrix} -w_{,xx} \\ -w_{,yy} \\ -2w_{,xy} \end{bmatrix} = [B_{\alpha}]\{\alpha\} \quad [59]
$$
Et la matrice de rigidité de l'élément est :

$$[k] = [T]' \int_S [B_a]' [D] [B_a] dS [T] \quad [60]$$

La matrice $[B_a]$ a l'expression suivante :

$$[B_a] = \begin{bmatrix} 0 & 0 & 0 & P_{4,xx} & P_{5,xx} & P_{6,xx} & P_{7,xx} & P_{8,xx} & P_{9,xx} \\ 0 & 0 & 0 & P_{4,yy} & P_{5,yy} & P_{6,yy} & P_{7,yy} & P_{8,yy} & P_{9,yy} \\ 0 & 0 & 0 & 2P_{4,xy} & 2P_{5,xy} & 2P_{6,xy} & 2P_{7,xy} & 2P_{8,xy} & 2P_{9,xy} \end{bmatrix} \quad [61]$$

avec :

$$P_{4,xx} = \frac{1}{2A^2} b_1 b_2 \qquad P_{4,yy} = \frac{1}{2A^2} c_1 c_2 \qquad P_{4,xy} = \frac{1}{4A^2} (c_1 b_2 + b_1 c_2) \quad [62]$$

On obtient les dérivées secondes des termes $P_5 = L_2 L_3$ et $P_6 = L_3 L_4$ par permutations circulaires des indices.

Les dérivées secondes du terme :

$$P_7 = L_1^2 L_2 + \frac{3}{2}(1 - \mu_3)L_1^2 L_2 L_3 + \left(\frac{1 + 3\mu_3}{2}\right)(L_1 L_2 L_3^2 - L_1 L_2^2 L_3) \quad [63]$$

sont :

$$P_{7,xx} = \frac{1}{4A^2} \begin{bmatrix} b_1^2(2L_2 + 3(1-\mu_3)L_2 L_3) + (b_3^2 L_2 - b_2^2 L_3)L_1(1+3\mu_3) \\ + b_1 b_2(4L_1 + 6(1-\mu_3)L_1 L_3 + (1+3\mu_3)(L_3^2 - 2L_2 L_3)) \\ + b_2 b_3(3(1-\mu_3)L_1^2 + 2(1+3\mu_3)(L_1 L_3 - L_1 L_2)) \\ + b_3 b_1(6(1-\mu_3)L_1 L_2 + (1+3\mu_3)(2L_2 L_3 - L_2^2)) \end{bmatrix} \quad [64]$$

$$P_{7,yy} = \frac{1}{4A^2} \begin{bmatrix} c_1^2(2L_2 + 3(1-\mu_3)L_2 L_3) + (c_3^2 L_2 - c_2^2 L_3)L_1(1+3\mu_3) \\ + c_1 c_2(4L_1 + 6(1-\mu_3)L_1 L_3 + (1+3\mu_3)(L_3^2 - 2L_2 L_3)) \\ + c_2 c_3(3(1-\mu_3)L_1^2 + 2(1+3\mu_3)(L_1 L_3 - L_1 L_2)) \\ + c_3 c_1(6(1-\mu_3)L_1 L_2 + (1+3\mu_3)(2L_2 L_3 - L_2^2)) \end{bmatrix} \quad [65]$$
$$ P_{7,xy} = \frac{1}{4A^2} \left[ \begin{array}{l} c_1 b_1 (2L_2 + 3(1-\mu_3)L_2 L_3) + (c_3 b_3 L_2 - c_2 b_2 L_3)L_4 (1+3\mu_3) \\ + (c_1 b_2 + b_1 c_2) \left( 2L_4 + 3(1-\mu_3)L_4 L_3 + \left( \frac{1+3\mu_3}{2} \right) (L_3^2 - 2L_2 L_3) \right) \\ + (c_2 b_3 + b_2 c_3) \left( \frac{3}{2}(1-\mu_3)L_4^2 + (1+3\mu_3)(L_4 L_3 - L_3 L_2) \right) \\ + (c_1 b_3 - b_1 c_3) \left( 3(1-\mu_3)L_4 L_2 + \left( \frac{1+3\mu_3}{2} \right) (2L_2 L_3 - L_2^2) \right) \end{array} \right] \quad [66] $$

et là encore, les dérivées secondes de $P_8$ et $P_9$ s'obtiennent par permutations circulaires des indices.

A partir de l'expression de $[B_α]$, la matrice de rigidité élémentaire est obtenue par intégration numérique de Gauss-Radau sur le triangle.

### 3.2.2. Matrice masse

La matrice masse est donnée par :

$$ [M] = [T]' [M_α] [T] \quad [67] $$

avec :

$$ [M_α] = \int_S \{P\} \langle P \rangle h \rho \, dS \quad [68] $$

## 4. Élément de coque triangulaire

L'élément triangulaire de coque, permettant l'étude par éléments finis d'assemblages de plaques non coplanaires ou de coques discrétisées en facettes triangulaires planes est obtenu par superposition de l'élément triangulaire T3D avec degrés de liberté rotationnels pour les effets de membrane et de l'élément triangulaire de Specht pour les effets de flexion. L'élément fini obtenu a 3 nœuds et 6 degrés de liberté par nœud.

### 4.1. Validation

Pour valider l'élément de coque triangulaire obtenu, nous réalisons différents calculs sur des tests standard proposés dans la littérature et présentés en détail dans l'ouvrage de Batoz et Dhatt [BAT 92].
### 4.1.1. Panneau cylindrique soumis à son propre poids

Les bords courbes du panneau reposent sur des diaphragmes rigides dans leur plan et les bords droits sont libres (figure 13).

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Panneau cylindrique soumis à son propre poids. Géométrie-Discrétisation</figcaption></figure>

Les caractéristiques du problème sont :

$$ 
\begin{aligned}
L &= 6\text{m} ; R = 3\text{m} ; h = 0.03\text{m} ; \varphi = 40^{\circ} \\
E &= 3 \times 10^{10} \text{Pa} ; v = 0. ; f_z = -0.625 \times 10^4 \text{Pa}
\end{aligned}
 $$

Compte tenu des symétries, un quart du cylindre est discrétisé et les conditions aux limites appliquées sont :

conditions d'appuis sur le diaphragme :

$$ u = w = \theta_y = 0 \text{ sur AD} $$

conditions de symétrie :

$$ u = \theta_y = \theta_z = 0 \text{ sur CD} $$

$$ v = \theta_x = \theta_z = 0 \text{ sur CB} $$

Le tableau 8 donne (en centimètres), pour différentes finesses de maillage, les valeurs du déplacement vertical du point B obtenues avec l'élément coque proposé (Coque T3D), les résultats obtenus avec d'autres éléments [BAT 92] et les résultats analytiques donnés par la théorie des coques « profondes » et des coques surbaissées [BAT 92].

Le tableau 9 donne les résultats pour le déplacement vertical du point C.
<table><caption>**Tableau 8. Panneau cylindrique-Déplacement vertical du point B (w<sub>B</sub>)**</caption><thead><tr><th>Maillage</th><th>N=2</th><th>N=4</th><th>N=8</th><th>N=10</th><th>N=12</th></tr></thead><tbody><tr><td>Coque T3D</td><td>-4.032</td><td>-3.389</td><td>-3.517</td><td>-3.549</td><td>-3.564</td></tr><tr><td>Q4γ24</td><td>-4.285</td><td>-3.418</td><td>-3.541</td><td></td><td>-3.602</td></tr><tr><td>DKQ24</td><td>-4.392</td><td>-3.424</td><td>-3.528</td><td></td><td>-3.568</td></tr><tr><td>DKT12</td><td>-3.773</td><td>-2.499</td><td>-3.129</td><td>-3.275</td><td>-3.365</td></tr><tr><td>Théorie des coques</td><td colspan="5">-3.61 (profondes), -3.703 (surbaissées)</td></tr></tbody></table>

<table><caption>**Tableau 9. Panneau cylindrique-Déplacement vertical du point C (w<sub>C</sub>)**</caption><thead><tr><th>Maillage</th><th>N=2</th><th>N=4</th><th>N=8</th><th>N=10</th><th>N=12</th></tr></thead><tbody><tr><td>Coque T3D</td><td>0.529</td><td>0.494</td><td>0.524</td><td>0.530</td><td>0.533</td></tr><tr><td>Q4γ24</td><td>0.629</td><td>0.508</td><td>0.531</td><td></td><td>0.542</td></tr><tr><td>DKQ24</td><td>0.593</td><td>0.513</td><td>0.529</td><td></td><td>0.535</td></tr><tr><td>DKT12</td><td>0.175</td><td>0.330</td><td>0.465</td><td>0.489</td><td>0.504</td></tr><tr><td>Théorie des coques</td><td colspan="5">0.541 (profondes), 0.525 (surbaissées)</td></tr></tbody></table>

### 4.1.2. Cylindre pincé avec diaphragme

On considère le cas d'un cylindre soumis à deux charges concentrées diamétralement opposées qui repose à ses extrémités sur deux diaphragmes rigides dans leur plan (figure 14).

Les caractéristiques du problème sont :

$$L=6m ; R=3m ; h=0.03m$$

$$E=3 \times 10^{10} \text{Pa} ; v=0.3 ; P=-1.\text{Pa}$$

Un huitième du cylindre est discrétisé et les conditions aux limites appliquées sont :

conditions d'appuis sur le diaphragme :

$$u=w=\theta_y=0 \text{ sur AD}$$

conditions de symétrie :

$$w=\theta_y=\theta_x=0 \text{ sur AB}$$

$$u=\theta_y=\theta_z=0 \text{ sur CD}$$

$$v=\theta_x=\theta_z=0 \text{ sur BC}$$
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Cylindre pincé. Géométrie-discrétisation</figcaption></figure>

Le tableau 10 donne, pour différentes finesses de maillage, les valeurs du déplacement vertical du point C ($\overline{w}_c = -Eh/P * w_c$) obtenues avec l'élément coque proposé (Coque T3D), les résultats obtenus avec d'autres éléments et les résultats analytiques donnés par la théorie des coques [BAT 92]. Le tableau 11 donne les résultats pour le déplacement longitudinal du point D ($\overline{v}_D = Eh/P * v_D$).

<table><caption>**Tableau 10.** Cylindre pincé. Déplacement vertical du point C ($\overline{w}_c$)</caption><thead><tr><th>Maillage</th><th>N=2</th><th>N=4</th><th>N=8</th><th>N=10</th><th>N=12</th><th>N=16</th></tr></thead><tbody><tr><td>Coque T3D</td><td>11.13</td><td>69.33</td><td>131.06</td><td>140.51</td><td>143.74</td><td>152.33</td></tr><tr><td>Q4γ24</td><td>6.10</td><td>60.93</td><td>121.95</td><td>135.18</td><td>143.37</td><td></td></tr><tr><td>DKQ24</td><td>14.07</td><td>101.24</td><td>154.84</td><td></td><td>164.95</td><td></td></tr><tr><td>DKT12</td><td>8.91</td><td>103.03</td><td>205.53</td><td>214.31</td><td>216.88</td><td></td></tr><tr><td>Théorie des coques</td><td colspan="6">164.24</td></tr></tbody></table>

<table><caption>**Tableau 11.** Cylindre pincé. Déplacement longitudinal du point D ($\overline{v}_D$)</caption><thead><tr><th>Maillage</th><th>N=2</th><th>N=4</th><th>N=8</th><th>N=10</th><th>N=12</th><th>N=16</th></tr></thead><tbody><tr><td>Coque T3D</td><td>1.3</td><td>3.78</td><td>3.96</td><td>3.94</td><td>3.78</td><td>3.94</td></tr><tr><td>Q4γ24</td><td>0.67</td><td>3.92</td><td>4.19</td><td>4.17</td><td>4.16</td><td></td></tr><tr><td>DKQ24</td><td>1.95</td><td>5.96</td><td>4.44</td><td></td><td>4.26</td><td></td></tr><tr><td>DKT12</td><td>0.8</td><td>4.12</td><td>5.00</td><td>4.59</td><td>4.42</td><td></td></tr><tr><td>Théorie des coques</td><td colspan="6">4.114</td></tr></tbody></table>
#### 4.1.3. Conclusions

L'élément coque construit en superposant l'élément de membrane T3D avec degrés de liberté rotationnels et l'élément de flexion de Specht [SPE 88] et Zienckiewicz [ZIE 91], pour les problèmes tests présentés, donnent des résultats satisfaisants.

## 4.2. Applications numériques

Pour illustrer les performances de l'élément développé, nous proposons différents calculs sur des structures simples, en statique et dynamique avec la recherche de la base modale. Les résultats obtenus sont comparés à ceux donnés par un code classique, le code de calcul Samcef [SAM 00].

### 4.2.1. Etude d'une cornière en statique

On considère un assemblage de deux plaques en acier ($E=2.10^{11}$Pa, $\rho=7\,800\text{kg/m}^3$) identiques, de longueur $l=3$ m, de largeur $b=1$ m et d'épaisseur $e_1 = e_2 = e = 0.1$ m (figure 15).

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Géométrie de la cornière</figcaption></figure>

La cornière formée de ces plaques est encastrée dans le plan $y=0$, elle est soumise à une force ponctuelle $\vec{P}=100\vec{x}+100\vec{z}$ (en Newton) appliquée à l'autre extrémité, au point A de l'arête commune.

La cornière est discrétisée en éléments finis coque (éléments plaque tridimensionnels) montrés figure 16. Deux discrétisations sont considérées, l'une avec 24 éléments l'autre avec 96 éléments.

On compare les déplacements du point A, point d'application de la force et les déplacements des points extrêmes dans la section d'application de la charge, respectivement les points B et C (figure 15), obtenus :
– d'une part par Samcef. Les éléments coque utilisés sont obtenus par superposition d'un élément de membrane de Marguerre et d'un élément de flexion DSQ de Batoz [BAT 90] ;

– d'autre part avec l'élément proposé (Coque T3D).

<figure><img src="image_19.png" /><figcaption>Figure 16. Discrétisation de la cornière</figcaption></figure>

<table><caption>Tableau 12. Résultats pour une discrétisation de la cornière en 24 éléments</caption><thead><tr><th>Point</th><th></th><th>Samcef</th><th>Coque T3D</th><th>ε% = <sup>|Sam - 1 SD|</sup>&frasl;<sub>SD</sub></th></tr></thead><tbody><tr><td rowspan="3">B</td><td>U<sub>B</sub></td><td>5.16</td><td>3.83</td><td>25.8</td></tr><tr><td>V<sub>B</sub></td><td>-0.962</td><td>-0.892</td><td>7.3</td></tr><tr><td>W<sub>B</sub></td><td>4.13</td><td>3.87</td><td>6.3</td></tr><tr><td rowspan="3">A</td><td>U<sub>A</sub></td><td>4.22</td><td>3.98</td><td>5.7</td></tr><tr><td>V<sub>A</sub></td><td>0.99</td><td>0.94</td><td>5.05</td></tr><tr><td>W<sub>A</sub></td><td>4.22</td><td>3.98</td><td>5.7</td></tr><tr><td rowspan="3">C</td><td>U<sub>C</sub></td><td>4.13</td><td>3.87</td><td>6.3</td></tr><tr><td>V<sub>C</sub></td><td>-0.962</td><td>-0.892</td><td>7.3</td></tr><tr><td>W<sub>C</sub></td><td>5.16</td><td>3.83</td><td>25.8</td></tr></tbody></table>
<table><caption>**Tableau 13.** Résultats pour une discrétisation de la cornière en 96 éléments</caption><thead><tr><th>Point</th><th></th><th>Samcef</th><th>Coque T3D</th><th>ε% = <sup>Sam - T3D</sup>&frasl;<sub>Sam</sub></th></tr></thead><tbody><tr><td rowspan="3">B</td><td>U<sub>B</sub></td><td>4.89</td><td>5.03</td><td>2.9</td></tr><tr><td>V<sub>B</sub></td><td>-1.08</td><td>-1.18</td><td>9.3</td></tr><tr><td>W<sub>B</sub></td><td>4.72</td><td>5.19</td><td>9.9</td></tr><tr><td rowspan="3">A</td><td>U<sub>A</sub></td><td>4.85</td><td>5.26</td><td>8.45</td></tr><tr><td>V<sub>A</sub></td><td>1.13</td><td>1.26</td><td>11.5</td></tr><tr><td>W<sub>A</sub></td><td>4.85</td><td>5.26</td><td>8.45</td></tr><tr><td rowspan="3">C</td><td>U<sub>C</sub></td><td>4.72</td><td>5.19</td><td>9.9</td></tr><tr><td>V<sub>C</sub></td><td>-1.08</td><td>-1.18</td><td>9.3</td></tr><tr><td>W<sub>C</sub></td><td>4.89</td><td>5.03</td><td>2.9</td></tr></tbody></table>

### 4.2.2. Modes propres et fréquences propres d'une cornière

Dans cette section, la cornière précédente est étudiée en dynamique. Les épaisseurs des plaques sont respectivement $e_1 = 0.01 \text{ m}$ et $e_2 = 0.03 \text{ m}$, la largeur des plaques est $b = 0.25 \text{ m}$; la cornière est donc dissymétrique. L'étude porte sur le calcul des modes propres et fréquences propres de la cornière encastrée dans la section $y=0$, libre à l'autre extrémité.

Le tableau 14 donne les 9 premières fréquences propres de la cornière. Les résultats présentés, sont obtenus :

- par un calcul Samcef,

- par un calcul avec les éléments proposés (Coque T3D),

- par un calcul avec assemblage d'éléments coque construits en superposant l'élément de membrane CST et l'élément de flexion de Specht [SPE 88] et Zienckiewicz [ZIE 91], nommé Coque CST. Dans ce cas, les rotations normales au plan des deux plaques constituant la cornière sont éliminées et seuls les nœuds de l'arête possèdent 6 degrés de liberté.

Les figures 18 à 22 montrent les 5 premières formes modales.

On vérifie par ces calculs que les résultats issus des différents types de modélisation sont proches.
<table><caption>**Tableau 14.** Fréquences propres en Hertz de la cornière discrétisée en 160 éléments</caption><thead><tr><th>Fréquences en Hertz</th><th>Samcef</th><th>Coque T3D</th><th>Coque CST</th><th>ε₁% = |Sam - T3D| / Sam</th><th>ε₁% = |Sam - CST| / Sam</th></tr></thead><tbody><tr><td>f₁</td><td>19.89</td><td>20.48</td><td>20.89</td><td>2.9</td><td>5.</td></tr><tr><td>f₂</td><td>29.8</td><td>31.5</td><td>32.32</td><td>5.7</td><td>8.5</td></tr><tr><td>f₃</td><td>49.7</td><td>52.5</td><td>52.87</td><td>5.6</td><td>6.4</td></tr><tr><td>f₄</td><td>81.5</td><td>86.6</td><td>85.7</td><td>6.3</td><td>5.2</td></tr><tr><td>f₅</td><td>117.8</td><td>130.4</td><td>130.3</td><td>10.7</td><td>10.6</td></tr><tr><td>f₆</td><td>123.2</td><td>138</td><td>134.6</td><td>12.</td><td>9.3</td></tr><tr><td>f₇</td><td>149.9</td><td>177.2</td><td>168.8</td><td>18.</td><td>12.6</td></tr><tr><td>f₈</td><td>172.8</td><td>190.9</td><td>188.8</td><td>10.5</td><td>9.3</td></tr><tr><td>f₉</td><td>178.9</td><td>201.7</td><td>196.9</td><td>12.7</td><td>10.1</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Cornière dissymétrique - Mode 1</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 18. Cornière dissymétrique - Mode 2</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 19. Cornière dissymétrique - Mode 3</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 20. Cornière dissymétrique - Mode 4</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 21. Cornière dissymétrique - Mode 5</figcaption></figure>
## 5. Conclusion

On a présenté dans cet article un nouvel élément fini triangulaire de membrane avec degrés de liberté rotationnels. L'interpolation rigoureuse qui permet d'exprimer le champ des déplacements membranaires en fonction des degrés de liberté rotationnels conduit à un élément qui passe le patch test d'ordre zéro. Les résultats obtenus pour différents tests de la littérature montrent qu'il a un comportement satisfaisant. Même si le taux de convergence est un peu moins bon que celui de certains autres éléments, il est cependant intéressant puisqu'il assure la convergence vers la solution exacte. L'intégration de la matrice rigidité est exacte. La formulation en déplacement et l'interpolation rigoureuse proposée ne nécessite pas de mettre en œuvre une formulation mixte et d'introduire des coefficients de forme.

En l'associant à un élément triangulaire de plaque pour prendre en compte les effets de flexion, on dispose d'un élément de coque à six degrés de liberté. Les degrés de liberté rotationnels, rendent son utilisation très simple pour l'étude de structures constituées par l'assemblage de plaques et pour l'étude des coques. Pour les applications présentées, on le superpose à l'élément triangulaire de flexion de Specht [SPE 88] et Zienckiewicz [ZIE 91]. Cet élément passe également le patch test, on dispose donc d'un élément de coque triangulaire plan performant.

## 6. Bibliographie

[ALL 84] ALLMAN D.J., « A compatible triangular element including vertex rotations for plane elasticity analysis », *Computer & Structures*, vol. 19, n°1-2, 1984, p.1-8.

[ALL 88a] ALLMAN D.J., « The constant strain triangle with drilling rotations : a simple prospect for shell analysis », *Mathematics of finite elements and applications VI*, LR Whiteman (ed) Academic Press N.Y, 1988, p 233-240.

[ALL 88b] ALLMAN D.J., « Evaluation of the constant strain triangle with drilling rotations », *International Journal for Numerical Methods in Engineering*, vol. 26, 1988, p. 2645-2655.

[AMI 92] AMINPOUR M.A., « An assumed-stress hybrid 4-node element with drilling degrees of freedom », *International Journal for Numerical Methods in Engineering*, vol. 33, 1992, p. 19-38.

[AYA 95a] AYAD R., « Une nouvelle approche variationnelle mixte-hybride pour les coques avec effet de cisaillement transversal et un nouveau concept de rotation autour de la normale pour améliorer la membrane », 2ᵉ *Colloque National en Calcul des Structures 95*, Editions Hermès, 1995, p. 257-262.

[AYA 95b] AYAD R., BATOZ J.L., DHATT G., « Un élément quadrilatéral de plaque basé sur une formulation mixte-hybride avec projection en cisaillement », *Revue Européenne des Eléments Finis*, vol. 4, n°4, 1995, p. 415-440.

[BAT 88] BATOZ J.L., DHATT G., « Revue et bilan des éléments de plaque de type Kirchhoff discret », *Calcul des Structures et I.A.*, vol. 2, Fouet et al. Eds Pluralis, 1988, p. 137-160.
[BAT 90] BATOZ J.L., DHATT G., *Modélisation des Structures par Eléments Finis*, vol. 1 : *Solides élastiques*, vol. 2 : *Poutres et Plaques*, Editions Hermès, 1990.

[BAT 92] BATOZ J.L., DHATT G., *Modélisation des Structures par Eléments Finis*, vol. 3 : *Coques*, Editions Hermès, 1992.

[BER 85] BERGAN P.G., FELIPPA CA., « A triangular membrane element with rotational degrees of freedom », *Computer Methods in Applied Mechanics and Engineering*, vol. 50, 1985, p.25-69.

[COO 74] COOK R.D., « Improved two-dimensional finite element », *Journal Structural Division ASCE*, vol. 100, 1974, p.1851-1865.

[COO 87] COOK R.D., « A plane hybrid element with rotational d.o.f. and adjustable stiffness », *International Journal for Numerical Methods in Engineering*, vol. 24, 1987, p.1499-1508.

[COO 91] COOK R.D., « Modified formulations for nine D.O.F plane triangles that include vertex rotations », *International Journal for Numerical Methods in Engineering*, vol. 31, 1991, p. 825-835.

[HUG 89] HUGHES T.J.R, BREZZI F., « On drilling degrees of freedom », *Computer Methods in Applied Mechanics and Engineering*, vol. 72, 1989, p. 105-121.

[HUG 95] HUGHES T.J.R, MASUD A., HARAZI I., « Dynamic analysis and drilling degrees of freedom », *International Journal for Numerical Methods in Engineering*, vol. 38, 1995, p. 3193-3210.

[IBR 90] IBRAHIMBEGOVIC A., TAYLOR R., WILSON E., « A robust quadrilateral membrane finite element with drilling degrees of freedom », *International Journal for Numerical Methods in Engineering*, vol. 31, 1990, p. 445-457.

[IBR 91] IBRAHIMBEGOVIC A., WILSON E., « Thick shell and solide finite elements with independent rotation fields », *International Journal for Numerical Methods in Engineering*, vol. 31, 1991, p. 1393-1414.

[IBR 93] IBRAHIMBEGOVIC A., FREY F., REBORA B., « Une approche unifiée de la modélisation des structures complexes : les éléments finis avec degré de liberté de rotation », *Revue Européenne des Eléments finis*, vol. 2, n°3, 1993, p. 257-286.

[IRO 80] IRONS B., AHMAD S., « *Techniques of finite elements* », John Wiley and Son, 1980.

[LIU 95] LIU M.L., TO S.W.S., « Vibration of structure by hybrid strain-based three nodes flat triangular shell elements », *Journal of Sound and Vibration*, vol. 18, n°5, 1995, p. 801-821.

[MAC 84] MACNEAL R.H., HARDER R.L., « A proposed standard set of problems to test finite element accuracy », *SDM Finite Element Validation Forum, 25th AIAA/ASME/ASCE/AHS Structures, Structural Dynamics and Materials Conference*, Palm Springs, CA, 1984.

[MAC 88] MACNEAL R.H., HARDER R.L., « A refined four-noded membrane element with rotational degrees of freedom », *Computers & Structures*, vol. 18, 1988, p. 75-84.

[OLS 79] OLSEN M.D., BEARDEN T.W., « A simple flat shell element revisited », *International Journal for Numerical Methods in Engineering*, vol. 14, 1979, p. 51-68.
[SAM 00] SAMCEF, *User Manual* 8.1, Samtech, Liège (Belgique), 2000.

[SPE 88] SPECHT B., « Modified shape functions for the three node plate bending element passing the patch test », *International Journal for Numerical Methods in Engineering*, vol. 26, 1988, p. 705-715.

[TAY 93] TAYLOR R.L. AURICHKIO F., « Linked interpolation for Reissner-Mindlin plate elements- Part I: a simple triangle », *International Journal for Numerical Methods in Engineering*, vol. 36, 1993, p. 3057-3066.

[TIM 51] TIMOSHENKO S., GOODIER J.N., *Theory of Elasticity*, McGraw-hill, New York, 1951.

[XU 94] XU Z., ZIENKIEWICZ O.C., ZENG L.F., « Linked interpolation for Reissner-Mindlin plate elements - Part III: an alternative quadrilateral », *International Journal for Numerical Methods in Engineering*, vol. 37, 1994, p. 1437-1443.

[YUN 89] YUNUS S.M., SAIGAL S., COOK R.D., « On improved hybrid finite elements with rotational degrees of freedom », *International Journal for Numerical Methods in Engineering*, vol. 28, n°4, p. 785-800.

[ZIE 91] ZIENKIEWICZ O.C., TAYLOR R.L., *The finite element method*, Four Edition, vol. 1 et vol. 2, Mc Graw Hill, 1991.

[ZIE 93] ZIENKIEWICZ O.C., XU Z., ZENG L.F., SAMUELSON A., WIBERG N.E., « Linked interpolation for Reissner-Mindlin plate elements - Part I: a simple quadrilateral », *International Journal for Numerical Methods in Engineering*, vol. 36, 1993, p. 3043-3056.