5. Construction du système matriciel

Nous allons maintenant présenter les différents ingrédients menant à la construction du système matriciel qui permettra de résoudre le problème.

5.1. Nouvelle notation (notation de Voigt)

Afin de bien comprendre la construction des termes discrets dans la méthode des éléments finis, nous allons utiliser une notation plus compacte:

\(\langle V\rangle\) est un vecteur ligne

\(\left\lbrace V\right\rbrace\) est un vecteur colonne

\(\left[A\right]\) est une matrice

Ainsi l’interpolation géométrique s’écrit suivant les trois dimensions de l’espace:

(5.1)\[{x}^{e}={x}_{1}^{e}=\langle {x}_{1,i}^{e}\rangle .\left\lbrace {N}_{i}^{e}\right\rbrace =\langle {N}_{i}^{e}\rangle .\left\lbrace {x}_{1,i}^{e}\right\rbrace {y}^{e}={y}_{1}^{e}=\langle {y}_{1,i}^{e}\rangle .\left\lbrace {N}_{i}^{e}\right\rbrace =\langle {N}_{i}^{e}\rangle .\left\lbrace {y}_{1,i}^{e}\right\rbrace {z}^{e}={z}_{1}^{e}=\langle {z}_{1,i}^{e}\rangle .\left\lbrace {N}_{i}^{e}\right\rbrace =\langle {N}_{i}^{e}\rangle .\left\lbrace {z}_{1,i}^{e}\right\rbrace\]

Ou encore de manière plus compacte sous forme vectorielle:

(5.2)\[\left\lbrace {x}^{e}\right\rbrace =\left[{N}_{i}^{e}\right].\left\lbrace {x}_{i}^{e}\right\rbrace =\langle {x}_{i}^{e}\rangle .{\left[{N}_{i}^{e}\right]}^{T}\]

Avec la matrice \(N\) des fonctions de forme. En considérant un élément à deux noeuds, on obtient sous forme développée:

(5.3)\[\begin{split}\left\lbrace \begin{array}{c}{x}^{e}\\ {y}^{e}\\ {z}^{e}\end{array}\right\rbrace =\left[\begin{array}{cccccc}{N}_{1}^{e}& 0& 0& {N}_{2}^{e}& 0& 0\\ 0& {N}_{1}^{e}& 0& 0& {N}_{2}^{e}& 0\\ 0& 0& {N}_{1}^{e}& 0& 0& {N}_{2}^{e}\end{array}\right].\left\lbrace \begin{array}{c}{x}_{1}^{e}\\ {y}_{1}^{e}\\ {z}_{1}^{e}\\ {x}_{2}^{e}\\ {y}_{2}^{e}\\ {z}_{2}^{e}\end{array}\right\rbrace\end{split}\]

5.2. Système discrétisé

On se place dans le cas hyperélastique en petites déformations, le problème de mécanique à résoudre on écrit de manière plus compacte:

(5.4)\[\textrm{Trouver} u\in {E}^{h} \textrm{avec} a(u,\tilde{u})+l(\tilde{u})=0\]

Avec \(a(u,\tilde{u})\) une forme bilinéaire, symétrique qui représente l’énergie potentielle de la structure et \(l(\tilde{u})\) le potentiel [1] des efforts volumiques et surfaciques:

(5.5)\[\begin{split}\begin{array}{}a(u,\tilde{u})=\underset{{\Omega}^{h}}{\int}\varepsilon (\tilde{u}):\sigma (u).d{\Omega}^{h}\\ l(\tilde{u})=\underset{{\Omega}^{h}}{\int}f.\tilde{u}.d{\Omega}^{h}+\underset{{\Gamma}_{N}^{h}}{\int}g.\tilde{u}.d{\Gamma}^{h}\end{array}\end{split}\]

La discrétisation consiste à choisir une base de l’espace \({\Omega}^{h}\) et à calculer numériquement les termes de la matrice \(A\) et du vecteur \(L\) . Pour cela, on exprime la forme bilinéaire \(a(.,.)\) et la forme linéaire \(l(.)\) comme une somme sur des éléments, définis par découpage du domaine de base:

(5.6)\[\begin{split}\lbrace \begin{array}{cc}a({u}_{i},{\tilde{u}}_{j})& =\sum_{\text{éléments}{\Omega}^{e}}\underset{{\Omega}^{e}}{\int}{\sigma}_{kl}({u}_{i}).{\varepsilon}_{kl}({\tilde{u}}_{j}).d{\Omega}^{e}\\ l({\tilde{u}}_{i})& =\sum_{\text{éléments}{\Omega}^{e}}\underset{{\Omega}^{e}}{\int}{f}_{i}.{\tilde{u}}_{i}.d{\Omega}^{e}+\underset{{\Gamma}_{N}^{e}}{\int}{g}_{i}.{\tilde{u}}_{i}.d{\Gamma}_{N}^{e}\end{array}\end{split}\]

Les termes \({A}_{ij}\) , qui représentent l’interaction entre deux degrés de liberté \(i\) et \(j\) sont construits en assemblant ( l’opération notée \(\sum_{\text{éléments}{\Omega}^{e}}(\mathrm{...})\) ) les contributions provenant de chacun des éléments qui contiennent les nœuds correspondants ; on procède de la même façon pour construire le vecteur second membre \({L}_{i}\) . Ces contributions, appelées termes élémentaires, sont calculées lors d’une boucle sur les éléments et ne dépendent que des seules variables de l’élément \({\Omega}^{e}\) :

(5.7)\[\begin{split}\lbrace \begin{array}{cc}{a}^{e}& =\underset{{\Omega}^{e}}{\int}{\sigma}_{kl}.{\varepsilon}_{kl}.d{\Omega}^{e}\\ {l}^{e}& =\underset{{\Omega}^{e}}{\int}{f}_{i}.{w}_{i}.d{\Omega}^{e}+\underset{{\Gamma}_{N}^{e}}{\int}{g}_{i}.{w}_{i}.d{\Gamma}_{N}^{e}\end{array}\end{split}\]

La relation entre le tenseur des contraintes de Cauchy \(\sigma\) et les déplacements \(u\) est donnée par la relation de comportement, et est indépendante de l’écriture de la formulation variationnelle. Dans le cas élastique, on a:

(5.8)\[{\sigma}_{ij}({w}_{i})={\Lambda}_{ijkl}.{\varepsilon}_{\mathrm{ku}}({w}_{i})\]

\({\Lambda}_{ijkl}\) est le tenseur d’élasticité de Hooke. Cette forme tensorielle est peu pratique, on utilise préférentiellement la notation de Voigt, qui permet d’écrire:

(5.9)\[\sigma :\varepsilon =\langle \sigma \rangle .\lbrace \varepsilon \rbrace\]

En coordonnées cartésiennes, on a:

(5.10)\[\langle \sigma \rangle =\langle \begin{array}{cccccc}{\sigma}_{xx}& {\sigma}_{yy}& {\sigma}_{zz}& {\sigma}_{xy}& {\sigma}_{xz}& {\sigma}_{yz}\end{array}\rangle\]

Et une forme modifiée des composantes de déformation pour permettre d’exprimer le produit contracté, soit:

(5.11)\[\langle \varepsilon \rangle =\langle \begin{array}{cccccc}{\varepsilon}_{xx}& {\varepsilon}_{yy}& {\varepsilon}_{zz}& 2.{\varepsilon}_{xy}& 2.{\varepsilon}_{xz}& 2.{\varepsilon}_{yz}\end{array}\rangle\]

Remarque importante:

Dans l’intégration des lois de comportement, les composantes de cisaillement des contraintes et de déformations utilisés par Code_Aster sont:

\(\langle \sigma \rangle =\langle \begin{array}{cccccc}{\sigma}_{xx}& {\sigma}_{yy}& {\sigma}_{zz}& \sqrt{2}.{\sigma}_{xy}& \sqrt{2}.{\sigma}_{xz}& \sqrt{2}.{\sigma}_{yz}\end{array}\rangle\)

\(\langle \varepsilon \rangle =\langle \begin{array}{cccccc}{\varepsilon}_{xx}& {\varepsilon}_{yy}& {\varepsilon}_{zz}& \sqrt{2}.{\varepsilon}_{xy}& \sqrt{2}.{\varepsilon}_{xz}& \sqrt{2}.{\varepsilon}_{yz}\end{array}\rangle\)

Le produit de ces deux vecteurs donne bien le même résultat que le double produit contracté ( ).

Avec cette nouvelle notation, nous avons en élasticité:

(5.12)\[\lbrace \sigma \rbrace =[A].\lbrace \varepsilon \rbrace\]

Nous repartons de l’écriture EF du champ des déplacements:

(5.13)\[\left\lbrace {u}^{e}\right\rbrace =\left[{N}_{i}^{e}\right].\left\lbrace {u}_{i}^{e}\right\rbrace =\langle {u}_{i}^{e}\rangle .{\left[{N}_{i}^{e}\right]}^{T}\]

Et, de manière analogue, le champ des déplacements virtuels:

(5.14)\[\left\lbrace {\tilde{u}}^{e}\right\rbrace =\left[{N}_{i}^{e}\right].\left\lbrace {\tilde{u}}_{i}^{e}\right\rbrace =\langle {\tilde{u}}_{i}^{e}\rangle .{\left[{N}_{i}^{e}\right]}^{T}\]

Par souci de simplification des notations, on omettra la référence à l’élément. Il est nécessaire d’xprimer le tenseur des déformations (virtuelles ou réelles):

(5.15)\[\left\lbrace \varepsilon \right\rbrace =\left[B\right].\left\lbrace u\right\rbrace =\langle u\rangle .{\left[B\right]}^{T}\]

On obtient alors pour la matrice relative à la forme bilinéaire:

(5.16)\[\begin{array}{cc}[a]& =\langle \tilde{u}\rangle .\underset{{\Omega}^{e}}{\int}{\left[B\right]}^{T}.[\Lambda ].\left[B\right].d{\Omega}^{e}.\left\lbrace u\right\rbrace \end{array}\]

Les matrices \(\left[B\right]\) et \([\Lambda ]\) contiennent la non-linéarité éventuelle du comportement et dépendront des déplacements:

1. \(\left[B\right]\) est une fonction des déplacements si on se trouve dans la situation des grandes déformations ou des grandes transformations (grandes rotations et/ou grands déplacements).

2. \([\Lambda ]\) est la matrice de comportement. Elle devient dépendante des déplacements (et d’autres variables) dans le cas des comportements non-linéaires et/ou non-élastiques.

Dans ces deux cas, le processus de résolution des équations impliquera un traitement spécifique (le plus souvent, une linéarisation de type Newton-Raphson). De manière analogue, on obtiendra facilement la forme élémentaire pour le second membre.

5.3. Calcul des termes élémentaires

Les termes élémentaires à calculer sont de la forme:

(5.17)\[\underset{{\Omega}_{e}}{\int}f(u(x),\frac{\partial u(x)}{\partial x}).dx\]

Trois types d’opérations sont à effectuer:

1. la transformation des dérivées par rapport à \(x\) en dérivées par rapport à \(\xi\) ;

2. le passage d’une intégration sur l’élément réel à une intégration sur l’élément de référence,

3. la réalisation numérique de cette intégration qui est généralement faite par une formule de quadrature.

5.3.1. Transformation des dérivées

La transformation des dérivées s’effectue grâce à la matrice jacobienne \(J\) , d’après la règle de dérivation en chaîne:

(5.18)\[\frac{\partial {u}_{\alpha}}{\partial {x}_{\beta}}=\frac{\partial {\xi}_{\gamma}}{\partial {x}_{\beta}}.\frac{\partial {u}_{\alpha}}{\partial {\xi}_{\gamma}}={J}^{\text{-1}}.{\left[\frac{\partial N}{\partial \xi }\right]}^{T}.{u}_{\alpha}^{\text{nod}}\]

où \({u}_{\alpha}^{\text{nod}}\) est le vecteur des valeurs nodales de la composante \(\alpha\) du déplacement.

5.3.2. Changement de domaine d’intégration

Le passage à l’intégration sur l’élément de référence s’effectue en multipliant l’intégrande par le déterminant de la matrice jacobienne, appelé jacobien:

(5.19)\[\underset{{\Omega}_{e}}{\int}f(u(x),\frac{\partial u(x)}{\partial x}).dx=\underset{{\Omega}_{r}}{\int}f(u(\xi ),\frac{\partial u(\xi )}{\partial \xi }).\det(J(\xi )).d\xi\]

Le passage de l’élément de référence à l’élément réel implique la bijectivité de la transformation \(\tau\) . Il faut donc \(\det(J(\xi ))\ne 0\) , ce qui implique que l’élément ne doit pas se retourner ou dégénérer (par exemple il ne faut pas que le quadrangle dégénère en triangle).

5.3.3. Intégration numérique

Dans certains cas particuliers, on peut calculer analytiquement les intégrales. Par exemple, pour un triangle en deux dimensions, le Jacobien est constant sur le triangle, et les intégrandes se ramènent à des monômes que l’on sait intégrer exactement grâce à la formule d’intégration numérique dite «de Gauss» [2] :

(5.20)\[\underset{0}{\overset{1}{\int}}\underset{0}{\overset{1-\xi }{\int}}{\xi}_{1}^{\alpha}.{\xi}_{2}^{\beta}.d{\xi}_{1}.d{\xi}_{2}=\frac{\alpha !\beta !}{(\alpha +\beta +2)!}\]

Cependant, ces cas particuliers sont rares, et on préfère évaluer numériquement les intégrales en faisant appel à des formules de quadrature. Celles-ci donnent une approximation de l’intégrale sous forme d’une somme pondérée des valeurs de l’intégrande en un certain nombre de points de l’élément appelés points d’intégration :

(5.21)\[\underset{{\Omega}_{r}}{\int}g(\xi ).d\xi \approx \sum_{g=1}^{r}{\omega}_{g}.g({\xi}_{g})\]

Les scalaires \({\omega}_{g}\) sont appelés les poids d’intégration, et les coordonnées \({\xi}_{g}\) sont les coordonnées des \(r\) points d’intégration dans l’élément de référence.

Dans les méthodes d’intégration de Gauss, les points et poids d’intégration sont déterminés de manière à intégrer exactement des polynômes d’ordre donné. C’est ce type de méthode que l’on utilise dans Code_Aster , les points d’intégration s’appellent alors des points de Gauss.

Le nombre de points de Gauss choisi permet d’intégrer exactement dans l’élément de référence. En fait, à cause de la non-linéarité éventuelle de la transformation géométrique ou de la dépendance spatiale des coefficients (par exemple pour des éléments déformées ou de second ordre), l’intégration n’est pas exacte dans l’élément réel.

Pour chaque élément \({\Omega}_{e}\) , on a su calculer les termes dits élémentaires: matrice élémentaire \({\mathrm{A}}^{\mathrm{e}}\) et vecteur élémentaire \({L}^{e}\) . La matrice \(A\) et le vecteur \(L\) sont obtenus par une procédure que l’on appelle l’assemblage des termes élémentaires.

Si on reprend la forme élémentaire de rigidité:

(5.22)\[[a]=\underset{{\Omega}^{e}}{\int}\lbrace \sigma \rbrace .\langle \varepsilon \rangle .d{\Omega}^{e}\]

L’intégration numérique implique que l’on évalue les contraintes et les déformations aux points d’intégration:

(5.23)\[[a]=\underset{{\Omega}^{e}}{\int}\lbrace \sigma \rbrace .\langle \varepsilon \rangle .d{\Omega}^{e}\approx \sum_{g=1}^{r}{\omega}_{g}.\lbrace \sigma ({\xi}_{g})\rbrace .\langle \varepsilon ({\xi}_{g})\rangle\]

Ce qui signifie que les contraintes et les déformations sont les plus exactes (ou les moins fausses) aux points d’intégration (champs dits «ELGA» dans Code_Aster ). Le simple fait d’extrapoler ces valeurs aux nœuds pour l’affichage introduit une erreur. Il s’agit d’ailleurs d’une méthode d’évaluation de l’erreur, appelée indicateur d’erreur de Zhu-Zienkiewicz.

En élasticité 2D, un triangle exhibant un jacobien constant, un seul point de Gauss est suffisant pour intégrer exactement les termes de la matrice et du second membre (s’il est constant).

Le coût calcul augmente avec le nombre de points d’intégration, particulièrement pour les lois de comportement non-linéaires. Par exemple, un hexaèdre à 27 nœuds a besoin de 27 points de Gauss pour intégrer les quantités. Il arrive donc fréquemment qu’on «sous-intègre», c’est à dire qu’on utilise moins de points d’intégration que le minimum requis, commettant ainsi une erreur que l’on compensera éventuellement par un maillage plus fin. En plus de cette erreur systématique, cette sous-intégration doit être faite avec précaution car elle peut produire des défauts de rang de la matrice et rendre ainsi le système linéaire non inversible.

[1]

Le potentiel des efforts extérieurs ne dépend pas du déplacement de la structure, c’est ce qu’on appelle un chargement mort ou non-suiveur. Dans les cas des grandes déformations, les chargements de type pression ne peuvent pas respecter cette hypothèse.

[2]

Par abus de langage, on appelle fréquemment les schémas d’intégration numériques « schémas de Gauss » bien qu’il en existe de plusieurs sortes (Hammer pour le triangle, Gauss-Radau, Newton-Cotes, etc.).