3. Formulations variationnelles#

3.1. Équilibre statique linéaire et non-linéaire#

On va exprimer les équations d’équilibre statique linéaire et non-linéaire sous forme variationnelle. On se place dans le cadre de la théorie des grandes déformations. Le tenseur des déformations de Green-Lagrange \({E}^{u}\) vaut:

(3.1)#\[{E}^{u}=\frac{1}{2}\left({F}^{T}F-g\right)\]

Avec \(F\) le tenseur gradient de déplacement classique (voir § 2.4 ) et \(g\) le tenseur métrique.

Afin de pouvoir introduire l’ensemble des méthodes de stabilisation et de correction des verrouillages, on va partir d’une formulation variationnelle mixte de Hu-Washizu [Bib2] et [Bib17], écrite sur la configuration initiale et en grandesdéformations:

(3.2)#\[{\pi}^{\mathit{HW}}(u,E,S)={\int}_{{\Omega}_{0}}W(E)d{\Omega}_{0}+{\int}_{{\Omega}_{0}}S:\left({E}^{u}-E\right)d{\Omega}_{0}-{\pi}_{\text{ext}}(u)\]

Avec \(W(E)\) l’énergie de déformation de la structure et \({\pi}_{\text{ext}}(u)\) le travail des efforts extérieurs. L’application de la méthode EAS consiste à reparamétrer le tenseur des déformations de Green-Lagrange tel que:

(3.3)#\[E={E}^{u}+{E}^{\mathit{EAS}}\]

Ce qui nous donne une nouvelle expression de () qui transforme la dépendance de \(E\) vers \({E}^{\mathit{EAS}}\) :

(3.4)#\[{\stackrel{~}{\pi}}^{\mathit{HW}}(u,{E}^{\mathit{EAS}},S)={\int}_{{\Omega}_{0}}W({E}^{u}+{E}^{\mathit{EAS}})d{\Omega}_{0}+{\int}_{{\Omega}_{0}}S:{E}^{\mathit{EAS}}d{\Omega}_{0}-{\pi}_{\text{ext}}(u)\]

Il est nécessaire de trouver des approximations correctes de ces trois champs, ce qui peut être compliqué en pratique. Cette formulation à trois champs est donc généralement transformée en formulation à deux champs en imposant la condition d’orthogonalité suivante:

(3.5)#\[{\int}_{{\Omega}_{0}}S:{E}^{\mathit{EAS}}d{\Omega}_{0}=0\]

Ce qui permet d’éliminer la contrainte. Au niveau discret, cette condition peut mener à des problèmes de convergence à cause d’un espace des contraintes \({S}^{h}\) trop restreint (voir [Bib3]). Ce problème peut être corrigé en imposant au niveau discret l’inclusion au moins des fonctions constantes par morceaux dans cet espace. Dans notre cas, nous allons procéder autrement pour permettre la prise en compte du pincement. Il y a deux situations:

  • s’il n’y a pas pincement, c’est-à-dire s’il n’y a pas de pression appliquée sur la surface inférieure ou supérieure de l’élément, on imposera la condition d’orthogonalité ();

  • s’il y a pincement et qu’on considère une pression sur la face supérieure \({P}_{u}\) et sur la surface inférieure \({P}_{l}\) , on choisit la contrainte de Piola-Kirchhoff de seconde espèce telle que

(3.6)#\[{\sigma}_{33}=-\frac{1}{2}(1-\zeta ){P}_{l}-\frac{1}{2}(1+\zeta ){P}_{u}\]

Avec \(\zeta\) la cordonnée paramétrique dans l’épaisseur.

On a donc:

(3.7)#\[\pi (u,{E}^{\mathit{EAS}})={\int}_{{\Omega}_{0}}W({E}^{u}+{E}^{\mathit{EAS}})d{\Omega}_{0}-{\pi}_{\text{ext}}(u)\]

Ou encore:

(3.8)#\[\pi (u,{E}^{\mathit{EAS}})={\pi}_{int}(u,{E}^{\mathit{EAS}})-{\pi}_{\text{ext}}(u)\]

La première variation de () permet d’obtenir la forme faible:

(3.9)#\[ \begin{align}\begin{aligned}\delta \pi (u,{E}^{\mathit{EAS}})=\delta {\pi}_{int}(u,{E}^{\mathit{EAS}})-\delta {\pi}_{\text{ext}}(u)\\\textrm{avec}\\\delta {\pi}_{int}(u,{E}^{\mathit{EAS}})={\int}_{{\Omega}_{0}}\left(\delta {E}^{u}+\delta {E}^{\mathit{EAS}}\right):Sd{\Omega}_{0}\end{aligned}\end{align} \]

Le problème est non-linéaire. On va appliquer la méthode de Newton (voir [R5.03.01]), avec les notations de ce document, \(\delta {\pi}_{int}\) sont les forces internes. Il faut donc exprimer les différentes variations en fonction des inconnues du système.

Pour la déformation \({E}^{u}\) , la variation étant petite, on ne considère que la partie linéaire. Dans le système des coordonnées paramétriques de l’élément, on a:

(3.10)#\[\delta {E}^{u}=\widehat{B}\delta U\]

Et pour la déformation EAS, sur l’élément, en retenant une approximation avec un seul paramètre \(\alpha\) :

(3.11)#\[\delta {E}^{\mathit{EAS}}=\widehat{\Gamma}\delta \alpha\]

On a donc:

(3.12)#\[\delta {\pi}_{int}={\int}_{{\Omega}_{0}}{\left(\widehat{B}\delta U\right)}^{T}Sd{\Omega}_{0}+{\int}_{{\Omega}_{0}}{\left(\widehat{\Gamma}\delta \alpha \right)}^{T}Sd{\Omega}_{0}\]

Et la force interne:

(3.13)#\[{L}_{int}=\underset{{L}_{int}^{u}}{\underset{⏟}{{\int}_{{\Omega}_{0}}{\widehat{B}}^{T}Sd{\Omega}_{0}}}+\underset{{L}_{int}^{\mathit{EAS}}}{\underset{⏟}{{\int}_{{\Omega}_{0}}{\widehat{\Gamma}}^{T}Sd{\Omega}_{0}}}\]

L’application de la méthode de Newton va faire apparaître quatre matrices tangentes:

(3.14)#\[\begin{split}\Delta {L}_{int}=\begin{array}{c}\delta {U}^{T}\frac{\partial {L}_{int}^{u}}{\partial U}\Delta U+\delta {\alpha}^{T}\frac{\partial {L}_{int}^{\mathit{EAS}}}{\partial \alpha }\Delta \alpha \\ +\delta {\alpha}^{T}\frac{\partial {L}_{int}^{\mathit{EAS}}}{\partial \alpha }\Delta U+\delta {U}^{T}\frac{\partial {L}_{int}^{u}}{\partial \alpha }\Delta \alpha \end{array}\end{split}\]

Avec les quatre matrices:

(3.15)#\[\Delta {L}_{int}=\left(\delta {U}^{T}{K}^{\mathit{uu}}+\delta {\alpha}^{T}{K}^{\alpha u}\right)\Delta U+\left(\delta {U}^{T}{K}^{u\alpha }+\delta {\alpha}^{T}{K}^{\alpha \alpha }\right)\Delta \alpha\]

Pour \({K}^{\mathit{uu}}\) est la matrice «classique» (que l’on retrouve dans les éléments isoparamétriques standards par exemple), elle comprend elle-même deux parties:

  • la matrice tangente matérielle \({K}_{m}^{\mathit{uu}}\) ;

  • la matrice tangente géométrique \({K}_{g}^{\mathit{uu}}\) dans le cas des grandes déformations.

Soit:

(3.16)#\[{K}^{\mathit{uu}}={\int}_{{\Omega}_{0}}{\widehat{B}}^{T}\frac{\partial S}{\partial E}d{\Omega}_{0}+{\int}_{{\Omega}_{0}}{\widehat{B}}^{T}\frac{\partial E}{\partial u}Sd{\Omega}_{0}={K}_{m}^{\mathit{uu}}+{K}_{g}^{\mathit{uu}}\]

3.2. Calcul de stabilité#

Ce paragraphe permet d’exposer l’expression du problème aux valeurs propres pour cvalculer la stabiltié d’une structure (flambement).La matrice géométrique intervient dans deux cas:

  • en non-linéaire quand on utilise des grandes déformations, le processus de linéarisation fait apparaître une contribution d’origine géométrique, voir ();

  • pour faire des calculs de stabilité (flambement), ce qui correspond dans code_aster à l’option RIGI_GEOM, voir [R7.05.01].

Dans les deux cas, la forme générale de cette matrice est la suivante:

(3.17)#\[{K}_{g}^{\mathit{uu}}={\int}_{{\Omega}_{0}}{\widehat{B}}^{T}\frac{\partial E}{\partial u}Sd{\Omega}_{0}\]

On va exprimer l’opérateur \(G\) tell que:

(3.18)#\[G={\int}_{{\Omega}_{0}}{\widehat{B}}^{T}\frac{\partial E}{\partial u}d{\Omega}_{0}\]

Comme l’approximation de \(B\) et \(E\) ont utilisé une décomposition en polynômes, on gardera la même logique de décomposition pour \(G\) (et le même schéma d’intégration numérique de Lobatto). De même on aura le même changement de repère entre l’espace paramétrique de l’élément et l’espace réel:

(3.19)#\[G=T\widehat{G}\]

La décomposition choisie est donc de la forme suivante:

(3.20)#\[\widehat{G}={\widehat{G}}^{0}+\zeta {\widehat{G}}^{\zeta}+{\zeta}^{2}{\widehat{G}}^{\zeta \zeta }+\xi {\widehat{G}}^{\xi}+\eta {\widehat{G}}^{\eta}+\xi \eta {\widehat{G}}^{\xi \eta }+\xi \zeta {\widehat{G}}^{\xi \zeta }+\eta \zeta {\widehat{G}}^{\eta \zeta }\]

Par ailleurs, on aura également une contribution venant de l’intégration réduite et une pour stabiliser:

(3.21)#\[ \begin{align}\begin{aligned}\widehat{G}={\widehat{G}}^{u}+{\widehat{G}}^{\mathit{STAB}}\\\textrm{avec}\\{\widehat{G}}^{u}={\widehat{G}}^{0}+\zeta {\widehat{G}}^{\zeta}+{\zeta}^{2}{\widehat{G}}^{\zeta \zeta }\\{\widehat{G}}^{\mathit{STAB}}=\xi {\widehat{G}}^{\xi}+\eta {\widehat{G}}^{\eta}+\xi \eta {\widehat{G}}^{\xi \eta }+\xi \zeta {\widehat{G}}^{\xi \zeta }+\eta \zeta {\widehat{G}}^{\eta \zeta }\end{aligned}\end{align} \]

On écrit les composantes de ce tenseur sous forme de vecteur:

(3.22)#\[\widehat{G}=\left\lbrace {\widehat{G}}_{\xi \xi },{\widehat{G}}_{\eta \eta },{\widehat{G}}_{\zeta \zeta },{\widehat{G}}_{\xi \eta },{\widehat{G}}_{\eta \zeta },{\widehat{G}}_{\xi \zeta }\right\rbrace\]

Pour chaque couple de nœud \(i\times j\) , on a les expressions explicites de ces composantes. Pour la membrane:

(3.23)#\[ \begin{align}\begin{aligned}\begin{split}{\widehat{G}}_{m,ij}^{0}=\left[\begin{array}{c}{a}_{1}^{i}{a}_{1}^{j}\\ {a}_{2}^{i}{a}_{2}^{j}\\ 0\\ {a}_{1}^{i}{a}_{2}^{j}+{a}_{1}^{j}{a}_{2}^{i}\\ 0\\ 0\end{array}\right]\end{split}\\\begin{split}{\widehat{G}}_{m,ij}^{\zeta}=\frac{1}{8}\left[\begin{array}{c}{h}_{3}^{i}{a}_{1}^{j}+{h}_{3}^{j}{a}_{1}^{i}\\ {h}_{2}^{i}{a}_{2}^{j}+{h}_{2}^{j}{a}_{2}^{i}\\ 0\\ {h}_{2}^{i}{a}_{1}^{j}+{h}_{3}^{i}{a}_{2}^{j}+{h}_{3}^{j}{a}_{2}^{i}+{h}_{2}^{j}{a}_{1}^{i}\\ 0\\ 0\end{array}\right]` :math:`{\widehat{G}}_{m,ij}^{\zeta \zeta }=\frac{1}{8}\left[\begin{array}{c}{h}_{3}^{i}{h}_{3}^{j}\\ {h}_{2}^{i}{h}_{2}^{j}\\ 0\\ {h}_{2}^{i}{h}_{3}^{j}+{h}_{3}^{i}{h}_{2}^{j}\\ 0\\ 0\end{array}\right]\end{split}\end{aligned}\end{align} \]

3.3. Calcul de dynamique#

Dans ce paragraphe, on exprime la masse du système pour calculer les problèmes de dynamique linéaire et non-linéaires transitoires ainsi que le calcul modal.