r3.07.10 Éléments finis mixtes hybrides via FENICSx – PLAQ_MITC#

Résumé :

Dans le but de compléter la bibliothèque d’éléments finis de plaque plans [R3.07.03] actuellement disponibles (DKT, DST, Q4G…), on se propose d’introduire un élément fini mixte hybride basé sur l’approche MITC (en anglais, Mixed Interpolation of Tensorial Component). Ce modèle permet d’effectuer des calculs sur des structures minces en respectant les hypothèses de Reissner-Mindlin.

Avertissement

Cette documentation est provisoire, les éléments ne sont pas suffisamment validés pour être utilisés dans les études AIP et l’élément PLAQ_MITC est hors du périmètre de qualification.

Table des matières

Équations constitutives et hypothèses#

Hypothèses du modèle de Reissner–Mindlin#

Le modèle de Reissner-Mindlin permet de représenter le comportement mécanique des plaques épaisses, en prenant explicitement en compte les effets de cisaillement transverse. Contrairement au modèle de Kirchhoff-Love, il ne suppose pas que les sections normales à la surface moyenne restent perpendiculaires après déformation.

Pour ce qui suit :

  • le déplacement transverse est noté \(w(x, y)\),

  • les rotations des fibres normales à la surface moyenne sont décrites par un champ vectoriel indépendant \(\bar{\theta}(x, y) = \left(\theta_x (x,y), \theta_y (x, y)\right)\).

L’épaisseur \(t\) de la plaque est supposée petite devant ses dimensions en plan, mais non négligeable, ce qui autorise l’introduction d’une déformation de cisaillement transverse. Celle-ci est définie par :

\[\bar{\gamma} = \nabla w - \bar{\theta},\]

où \(\nabla w = \left( \cfrac{\partial w}{\partial x}, \cfrac{\partial w}{\partial y} \right)\) représente le gradient du déplacement transverse.

Lois de comportement élastique#

Le comportement élastique du matériau est décrit à l’aide des paramètres classiques :

  • le module d’Young \(E\),

  • le coefficient de Poisson \(\nu\),

  • un facteur de correction du cisaillement \(\kappa\), généralement pris égal à \(\cfrac{5}{6}\) pour les plaques homogènes.

Tenseur de courbure#

Le tenseur de courbure \(\bar{\bar{\kappa}}\) représente les variations spatiales du champ de rotation \(\bar{\theta}\). Il est défini comme la partie symétrique du gradient de \(\bar{\theta}\) :

\[\bar{\bar{\kappa}} = \mathrm{sym}(\nabla \bar{\theta}) = \frac{1}{2} \left( \nabla \bar{\theta} + (\nabla \bar{\theta})^\top \right).\]

Énergie de flexion#

Les contraintes dues à la flexion \(\bar{\bar{M}}\) sont reliés au tenseur de courbure par une loi constitutive linéaire :

\[\bar{\bar{M}} = \mathbb{C} : \bar{\bar{\kappa}},\]

où \(\mathbb{C}\) désigne le tenseur de rigidité en flexion, exprimé par :

\[\mathbb{C} = D \left[ (1 - \nu)\, \mathbb{I}_4 + \nu\, \mathbb{I}_2 \otimes \mathbb{I}_2 \right], \qquad \text{avec} \qquad D = \frac{E t^3}{12(1 - \nu^2)}.\]

Les notations utilisées sont :

  • \(\mathbb{I}_4\) : tenseur identité d’ordre 4,

  • \(\mathbb{I}_2\) : tenseur identité d’ordre 2,

  • \(\otimes\) : produit tensoriel.

L’énergie de flexion correspondante est donnée par :

\[\psi_b = \frac{1}{2} D \left[ (1 - \nu)\, \mathrm{tr}(\bar{\bar{\kappa}}^2) + \nu\, \left( \mathrm{tr}(\bar{\bar{\kappa}}) \right)^2 \right]\]

Énergie de cisaillement transverse#

L’énergie élastique associée au cisaillement transverse s’écrit :

\[\psi_s = \cfrac{E \kappa t}{4(1 + \nu)} \left\| \bar{\gamma}_R \right\|^2,\]

où \(\bar{\gamma}_R\) désigne le champ de cisaillement « réduit », interpolé dans un espace fonctionnel adapté (par exemple de type Nédélec), afin de prévenir le phénomène de verrouillage numérique lors de la discrétisation par éléments finis.

Formulations variationnelles#

Champs inconnus et espaces d’approximation#

La modélisation adoptée repose sur la théorie de Reissner-Mindlin, enrichie par une interpolation mixte de type MITC (Mixed Interpolation of Tensorial Components), dans la variante proposée par Durán et Liberman [1].

Cette approche permet de corriger efficacement le verrouillage numérique lié au cisaillement transverse, assurant ainsi une meilleure précision dans la reproduction des déformations, même pour des plaques minces.

Les champs inconnus sont discrétisés dans les espaces d’éléments finis suivants :

\[w_h \in P_1, \quad \bar{\theta}_h \in [P_2]^2, \quad \bar{\gamma}_h, \, \bar{p}_h \in \mathrm{NED}_1,\]

où :

  • \(P_1\) désigne l’espace des éléments de Lagrange linéaires, c’est-à-dire les fonctions continues dont la restriction à chaque élément du maillage est un polynôme de degré inférieur ou égal à 1 ;

  • \([P_2]^2\) désigne un espace de champs vectoriels à deux composantes, chacune étant une fonction \(P_2\), soit un polynôme de degré inférieur ou égal à 2 par morceau ;

  • \(\mathrm{NED}_1\) correspond à l’espace de Nédélec de premier type, adapté à la discrétisation des champs vectoriels tangents.

Sur un élément quadrangulaire, cela se traduit par la répartition des degrés de liberté illustrée ci-dessous :

../../../../_images/function_space_scheme.svg

Fig. Répartition des degrés de liberté pour les champs \(w_h \in P_1\), \(\bar{\theta}_h \in [P_2]^2\) et \(\bar{\gamma}_h, \bar{p}_h \in \mathrm{NED}_1\) sur un élément quadrangulaire.

Cette discrétisation permet une représentation fidèle et stable des différents effets mécaniques mis en jeu, tout en prévenant les artefacts numériques tels que le verrouillage de cisaillement.

Énergie du système et couplage variationnel#

Soit \(\Omega \subset \mathbb{R}^2\) le domaine plan occupé par la plaque, discrétisé à l’aide d’un maillage conforme \(\mathcal{T}_h\) composé d’éléments quadrangulaires. L’ensemble des arêtes est noté \(\mathcal{E}_h\) et, pour toute arête \(E \in \mathcal{E}_h\), \(\bar{t}\) désigne le vecteur unitaire tangent.

La formulation variationnelle adoptée s’exprime par la minimisation d’une fonctionnelle intégrant l’ensemble des contributions du système :

  • l’énergie de flexion, liée aux moments fléchissants et à la courbure de la plaque ;

  • l’énergie de cisaillement transverse, calculée à partir d’un champ conforme \(\bar{\gamma}_h\) ;

  • un terme de couplage faible, imposant la cohérence entre \(\bar{\gamma}_h\) et le champ reconstruit à partir de \(w_h\) et \(\bar{\theta}_h\), via un multiplicateur de Lagrange \(\bar{p}_h\).

La fonctionnelle s’écrit pour le quadruplet discrétisé \(q_h = (w_h, \bar{\theta}_h, \bar{\gamma}_h, \bar{p}_h)\) :

\[\begin{split}\begin{aligned} \Pi(q_h) =\ & \frac{1}{2} \int_\Omega \bar{\bar{M}} : \bar{\bar{\kappa}} \, \mathrm{d}x + \cfrac{E \kappa t}{4(1+\nu)} \int_\Omega \left\| \bar{\gamma}_h \right\|^2 \, \mathrm{d}x \\ & + \sum_{E \in \mathcal{E}_h} \int_E \left( \bar{\gamma}(w_h, \bar{\theta}_h) - \bar{\gamma}_h \right)\!\cdot\!\bar{t} \, (\bar{p}_h \cdot \bar{t}) \, \mathrm{d}s - W_{\mathrm{ext}}, \end{aligned}\end{split}\]

où :

  • \(\bar{\gamma}(w_h, \bar{\theta}_h) = \nabla w_h - \bar{\theta}_h\) est le champ de cisaillement reconstruit ;

  • \(\bar{\gamma}_h\) : champ conforme utilisé dans la discrétisation ;

  • \(\bar{p}_h\) : multiplicateur de Lagrange agissant sur les arêtes ;

  • \(W_{\mathrm{ext}}\) : travail des forces extérieures.

Le terme de couplage traduit, sous forme faible, l’égalité entre les deux représentations du cisaillement transverse : le champ reconstruit \(\bar{\gamma}(w_h, \bar{\theta}_h)\) et le champ conforme \(\bar{\gamma}_h\). Cette contrainte, imposée en moyenne sur chaque arête \(E \in \mathcal{E}_h\) via \(\bar{p}_h\), s’écrit :

\[\int_E \left( \bar{\gamma}_h - \bar{\gamma}(w_h, \bar{\theta}_h) \right)\!\cdot\!\bar{t} \, (\bar{p}_h \cdot \bar{t}) \, \mathrm{d}s.\]

Elle dérive localement d’un Lagrangien :

\[\Pi_R(\gamma, \gamma_R, p) = \int_E \left( (\gamma_R - \gamma)\!\cdot\!\bar{t} \right)(p \cdot \bar{t}) \, \mathrm{d}s,\]

avec \(\gamma_R = \bar{\gamma}_h\), \(p = \bar{p}_h\) et \(\bar{t}\) tangent unitaire.

Ce mécanisme de couplage garantit la cohérence entre champs de cisaillement tout en assurant la stabilité numérique du schéma, et joue un rôle clé pour éviter le verrouillage dans le régime des plaques minces.

Motivation pour l’utilisation de MITC#

Problème de verrouillage du cisaillement#

Dans le modèle de Reissner–Mindlin, lorsque l’épaisseur \(t\) tend vers zéro, le comportement attendu est celui du modèle de Kirchhoff, pour lequel :

\[\bar{\gamma} = \nabla w - \bar{\theta} \to (0, 0).\]

Cependant, avec une discrétisation standard utilisant des éléments \(P_1\) (éléments finis continus de degré 1) pour le déplacement \(w\) et les rotations \(\bar{\theta}\), et sous conditions d’encastrement, on a :

  • \(\nabla w_h \in \mathrm{P}_0\) : champ constant par morceaux ;

  • \(\bar{\theta}_h \in [\mathrm{P}_1]^2\) : champ continu, nul sur \(\partial \Omega\).

Il devient alors impossible, sauf trivialement, d’obtenir \(\nabla w_h = \bar{\theta}_h\). Ce défaut d’approximation provoque un verrouillage numérique : la solution discrète est trop rigide et ne reproduit pas correctement les modes de flexion, en particulier pour des plaques minces.

Solution : les éléments MITC#

Les éléments MITC (Mixed Interpolation of Tensorial Components) contournent ce verrouillage en modifiant l’interpolation du cisaillement transverse.

Le principe consiste à projeter le cisaillement reconstruit \(\nabla w - \bar{\theta}\) sur un sous-espace adapté (souvent défini sur les arêtes ou à l’intérieur des éléments), assurant ainsi une compatibilité faible entre les champs.

Cette technique limite les contraintes parasites, restaure la cohérence du modèle et permet une convergence robuste et uniforme, y compris dans le régime des plaques minces, sans recours à un raffinement excessif du maillage.

Formulation variationnelle faible du problème#

On cherche les points stationnaires du fonctionnel \(\Pi\) en imposant que sa dérivée directionnelle s’annule pour toute variation admissible \(\tilde{q} = (\tilde{w}, \tilde{\bar{\theta}}, \tilde{\bar{\gamma}}, \tilde{\bar{p}})\) :

\[F(q_h; \tilde{q}) := D_{\tilde{q}}[\Pi(q_h)] = 0 \quad \forall \tilde{q}.\]

Le résidu (F) se décompose en quatre contributions bilinéaires :

\[\begin{split}\begin{aligned} \text{Flexion :} \quad a^{\mathrm{b}}(w_h, \bar{\theta}_h; \tilde{w}, \tilde{\bar{\theta}}) &= \int_\Omega \bar{M} : \bar{\kappa}(\tilde{w}, \tilde{\bar{\theta}}) \, \mathrm{d}x, \\ \text{Cisaillement :} \quad a^{\mathrm{s}}(\bar{\gamma}_h; \tilde{\bar{\gamma}}) &= \int_\Omega \mathbf{T}(\bar{\gamma}_h) \cdot \tilde{\bar{\gamma}} \, \mathrm{d}x, \\ \text{Couplage (primal) :} \quad a^{\mathrm{MITC}}_a(q_h; \tilde{\bar{p}}) &= \sum_{E \in \mathcal{E}_h} \int_E \left( \nabla w_h - \bar{\theta}_h - \bar{\gamma}_h \right)\!\cdot\!\bar{t} \, (\tilde{\bar{p}} \cdot \bar{t}) \, \mathrm{d}s, \\ \text{Couplage (dual) :} \quad a^{\mathrm{MITC}}_b(\bar{\gamma}_h; \tilde{w}, \tilde{\bar{\theta}}) &= \sum_{E \in \mathcal{E}_h} \int_E (\bar{\gamma}_h \cdot \bar{t}) \, \left( \nabla \tilde{w} - \tilde{\bar{\theta}} \right)\!\cdot\!\bar{t} \, \mathrm{d}s. \end{aligned}\end{split}\]

La formulation faible globale devient :

\[F(q_h; \tilde{q}) = a^{\mathrm{b}} + a^{\mathrm{s}} + a^{\mathrm{MITC}}_a + a^{\mathrm{MITC}}_b - {W}_{\mathrm{ext}}(\tilde{q}) = 0, \quad \forall \tilde{q}.\]

Le système linéaire issu de la linéarisation prend la forme bloc :

\[\begin{split}\begin{bmatrix} A & 0 & C \\ 0 & B & D \\ C^\top & D & 0 \end{bmatrix} \begin{bmatrix} \delta z \\ \delta \bar{\gamma} \\ \delta \bar{p} \end{bmatrix} = \begin{bmatrix} b_z \\ b_\gamma \\ b_p \end{bmatrix},\end{split}\]

avec :

  • \(\delta z = (\delta w, \delta \bar{\theta})\) : déplacement transverse et rotations ;

  • \(\delta \bar{\gamma}\) : cisaillement réduit, introduit pour améliorer la performance numérique ;

  • \(\delta \bar{p}\) : multiplicateur de Lagrange de la contrainte MITC ;

  • \(A, B, C, D\) : matrices issues de la linéarisation des contributions du résidu.

Ce système est résolu globalement pour obtenir les incréments du vecteur des inconnues. Dans le cas linéaire considéré ici, une seule itération de la méthode de Newton est suffisante.

Remarque:

Les blocks \(A\) et \(B\) sont issues de la linéarisation des formes bilinaires symmetriques \(a^{\mathrm{b}}\) et \(a^{\mathrm{s}}\), respectivement. La block \(D\) est une matrice diagonale. Le système matriciel global est donc symmetrique.

BIBLIOGRAPHIE#

[1] Hale JS, Brunetti M, Bordas SP, Maurini C. Simple and extensible plate and shell finite element models through automatic code generation tools. Computers & Structures. 2018 Oct 15;209:163-81.