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 \(\vector{\theta}(x, y) = \left(\theta_x (x,y), \theta_y (x, y)\right)\).

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

\[\vector{\gamma} = \gradScal{w} - \vector{\theta},\]

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

Loi 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 \(\tensTwo{\kappa}\) représente les variations spatiales du champ de rotation \(\vector{\theta}\). Il est défini comme la partie symétrique du gradient de \(\vector{\theta}\) :

\[\tensTwo{\kappa} = \operGradSyme{\vector{\theta}} = \frac{1}{2} \left( \gradVector{\vector{\theta}} + \transpose{\gradVector{\vector{\theta}}} \right).\]

Énergie de flexion#

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

\[\tensTwo{M} = \tensFour{D} : \tensTwo{\kappa},\]

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

\[\tensFour{D} = D \left[ (1 - \nu)\, \tensFourUnit + \nu\, \tensTwoUnit \otimes \tensTwoUnit \right], \qquad \text{avec} \qquad D = \frac{E h^3}{12(1 - \nu^2)}.\]

Les notations utilisées sont :

  • \(\tensFourUnit\) : tenseur identité d’ordre 4,

  • \(\tensTwoUnit\) : tenseur identité d’ordre 2,

  • \(\otimes\) : produit tensoriel.

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

\[\psi_b = \frac{1}{2} \tensTwo{M} : \tensTwo{\kappa}.\]

Énergie de cisaillement transverse#

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

\[\psi_s = \frac{1}{2} S \left( \vector{\gamma}_R \cdot \vector{\gamma}_R \right)\]

où \(S = \cfrac{E \kappa h}{2(1 + \nu)}\) et \(\vector{\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 [bib1] .

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 \vector{\theta}^h \in [P_2]^2, \quad \vector{\gamma}_R^h, \, \vector{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. 160 Répartition des degrés de liberté pour les champs \(w^h \in P_1\), \(\vector{\theta}^h \in [P_2]^2\) et \(\vector{\gamma}_R^h, \vector{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 au cisaillement.

Fonctions de Nédélec#

Ce sont des fonctions vectorielles principalement associées aux arêtes des éléments. Une propriété importante est qu’elles assurent la continuité de la composante tangentielle d’un champ vectoriel d’un élément à l’autre. Dans notre cas, cela permet de garantir que la projection du champ \(\gamma_R\) est cohérente à l’interface des éléments.

../../../../_images/nedelec_functions.svg

Fig. 161 Représentations des fonctions de Nédélec#

Énergie du système et couplage variationnel#

Soit \(\Omega \subset \setR^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\), \(\vector{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 \(\vector{\gamma}_R^h\) ;

  • un terme de couplage faible, imposant la cohérence entre \(\vector{\gamma}_R^h\) et le champ reconstruit à partir de \(w^h\) et \(\vector{\theta}^h\), via un multiplicateur de Lagrange \(\vector{p}^h\).

La fonctionnelle s’écrit pour le quadruplet discrétisé \(q^h = (w^h, \vector{\theta}^h, \vector{\gamma}_R^h, \vector{p}^h)\) :

\[\begin{split}\begin{aligned} \Pi(q^h) =\ & \frac{1}{2} \int_\Omega \tensFour{D} \tensTwo{\kappa} (\vector{\theta}^h) : \tensTwo{\kappa} (\vector{\theta}^h) \measDomain + \frac{1}{2} \int_\Omega S \vector{\gamma}_R^h \cdot \vector{\gamma}_R^h \measDomain \\ & + \sum_{E \in \mathcal{E}^h} \int_E \left( \vector{\gamma}(w^h, \vector{\theta}^h) - \vector{\gamma}_R^h \right)\!\cdot\!\vector{t} \, (\vector{p}^h \cdot \vector{t}) \measBound - W_{\text{ext}} \end{aligned}\end{split}\]

où :

  • \(\tensTwo{\kappa}(\vector{\theta}^h) = \operGradSyme{\vector{\theta}^h}\) est le tenseur de courbure ;

  • \(\vector{\gamma}(w^h, \vector{\theta}^h) = \gradScal{w}^h - \vector{\theta}^h\) est le champ de cisaillement reconstruit ;

  • \(\vector{\gamma}_R^h\) : champ conforme utilisé dans la discrétisation ;

  • \(\vector{p}^h\) : multiplicateur de Lagrange agissant sur les arêtes ;

  • \(W_{\text{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 \(\vector{\gamma}(w^h, \vector{\theta}^h)\) et le champ conforme \(\vector{\gamma}_R^h\). Cette contrainte, imposée sur chaque arête \(E \in \mathcal{E}^h\) via \(\vector{p}^h\), s’écrit sous la forme variationnelle suivante :

\[\int_E \left( \vector{\gamma}(w^h, \vector{\theta}^h) - \vector{\gamma}_R^h \right) \cdot \vector{t} \, (\vector{p}^h \cdot \vector{t}) \measBound\]

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 \(h\) tend vers zéro, le comportement attendu est celui du modèle de Kirchhoff, pour lequel :

\[\vector{\gamma} = \gradScal{w} - \vector{\theta} \to \vector{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 \(\vector{\theta}\), et sous condition d’encastrement, on a :

  • \(\gradScal{w}^h \in \mathrm{P}_0\) : champ constant par morceaux ;

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

Il devient alors impossible, sauf trivialement, d’obtenir \(\gradScal{w}^h = \vector{\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 \(\gradScal{w} - \vector{\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{\vector{\theta}}, \tilde{\vector{\gamma}}_R, \tilde{\vector{p}})\) :

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

Le résidu \(F = {F}_{\text{ext}} - {F}_{\mathrm{int}}\) où \({F}_{\mathrm{int}}\) se décompose en quatre contributions bilinéaires :

\[\begin{split}\begin{aligned} \text{Flexion :} \quad & a^{\mathrm{b}}(\vector{\theta}^h; \tilde{\vector{\theta}}) = \int_\Omega \tensFour{D} \tensTwo{\kappa}(\vector{\theta}^h) : \tensTwo{\kappa}(\tilde{\vector{\theta}}) \, \mathrm{d}x, \\ \text{Cisaillement :} \quad & a^{\mathrm{s}}(\vector{\gamma}_R^h; \tilde{\vector{\gamma}}_R) = \int_\Omega S \vector{\gamma}_R^h \cdot \tilde{\vector{\gamma}}_R \, \mathrm{d}x, \\ \text{Couplage (primal) :} \quad & a^{\mathrm{MITC}}_1(q^h; \tilde{\vector{p}}) = \sum_{E \in \mathcal{E}^h} \int_E \left( \gradScal{w}^h - \vector{\theta}^h - \vector{\gamma}_R^h \right) \cdot \vector{t} \, (\tilde{\vector{p}} \cdot \vector{t}) \, \mathrm{d}s, \\ \text{Couplage (dual) :} \quad & a^{\mathrm{MITC}}_2(\vector{p}^h; \tilde{q}) = \sum_{E \in \mathcal{E}^h} \int_E (\vector{p}^h \cdot \vector{t}) \, \left( \nabla \tilde{w} - \tilde{\vector{\theta}} - \tilde{\vector{\gamma}}_R \right) \cdot \vector{t} \, \mathrm{d}s. \end{aligned}\end{split}\]

La formulation faible globale devient :

\[F(q^h; \tilde{q}) = {F}_{\text{ext}} - \left( a^{\mathrm{b}} + a^{\mathrm{s}} + a^{\mathrm{MITC}}_1 + a^{\mathrm{MITC}}_2 \right) = 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 \vector{\gamma}_R \\ \delta \vector{p} \end{bmatrix} = \begin{bmatrix} f_z \\ f_\gamma \\ f_p \end{bmatrix},\end{split}\]

avec :

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

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

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

  • \(\discMatr{A}\), \(\discMatr{B}\), \(\discMatr{C}\) et \(\discMatr{D}\) : matrices issues de la linéarisation des contributions du résidu;

  • \(\discVect{f_z}\), \(\discVect{f_\gamma}\) et \(\discVect{f_p}\) : vecteurs qui contribuent au résidu global \(F\).

Le système local sert à assembler les matrices et vecteurs globaux. Le système global sera ensuite résolue pour obtenir les incréments des inconnues.

Remarque:

Les blocs \(\discMatr{A}\) et \(\discMatr{B}\) sont issus de la linéarisation des formes bilinaires symétriques \(a^{\mathrm{b}}\) et \(a^{\mathrm{s}}\), respectivement. Le bloc \(\discMatr{D}\) est une matrice diagonale. Le système matriciel global est donc symétrique.

Condensation statique#

Au niveau élémentaire, il est possible d’éliminer les variables associées aux multiplicateurs de Lagrange. Pour ce faire, on calcule le système réduit en termes des variables primaires \(\discVect{z}\).

\[\begin{split}\begin{aligned} \discMatr{A_s} \discVect{\delta z} &= \discVect{f_s} \\ \discMatr{A_s} &= \discMatr{A} + \discMatr{C} \discMatr{D}^{-1} \discMatr{B} \discMatr{D}^{-1} \discMatr{C}^\top \\ \discVect{f_s} &= \discVect{f_z} - \discMatr{C} \discMatr{D}^{-1} \discVect{f_{\gamma}} + \discMatr{C} \discMatr{D}^{-1} \discMatr{B} \discMatr{D}^{-1} \discVect{f_p} \end{aligned}\end{split}\]

Les blocs \(\discMatr{A_s}\) et \(\discVect{f_s}\) sont ensuite assemblés au niveau global pour former le système final à résoudre. Il faut noter que le bloc \(\discMatr{D}\) dans la matrice élémentaire est diagonale par construction, ce qui facilite son inversion locale sur chaque élément.

Références#

[bib1]

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.