3. Formulations variationnelles#

3.1. 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.

3.2. É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.

3.3. Motivation pour l’utilisation de MITC#

3.3.1. 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.

3.3.2. 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.

3.4. 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.

3.5. 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.