Méthode des éléments finis/Présentation générale

La méthode des éléments finis (MÉF) est une manière numérique de résoudre certains des problèmes de physique. C'est une méthode qui permet de déterminer une solution approchée sur un domaine spatial, c'est-à-dire qui permet de calculer un champ (de scalaires, de vecteurs, de tenseurs) qui correspond à certaines équations et à certaines conditions imposées.

La méthode consiste à découper le domaine spatial en petits éléments, également appelés mailles, et à rechercher une formulation simplifiée du problème sur chaque élément, c'est-à-dire à transformer le système d'équations quelconque en un système d'équations linéaires. Chaque système d'équations linéaires peut se représenter par une matrice. Les systèmes d'équations pour tous les éléments sont ensuite rassemblés, ce qui forme une grande matrice ; la résolution de ce système global donne la solution approchée au problème.

La quantité de calculs nécessite l'emploi d'un ordinateur.

Théorie et pratique

Les logiciels modernes utilisant la méthode des éléments finis bénéficient d'une interface graphique rendant leur utilisation relativement simple. Par ailleurs, un certain nombre de tâches sont automatisables. On peut donc quasiment lancer un calcul sur ordinateur sans connaître la méthode.

Cependant, le modèle utilisé risque d'être inadapté au problème, on aura donc un résultat très éloigné de la réalité. L'utilisateur doit avoir des connaissances suffisantes pour être en mesure de

• maîtriser le modèle, c'est-à-dire utiliser les options permettant de représenter le plus fidèlement possible la réalité ;
• contrôler la qualité du résultat, détecter les résultats manifestement erronés ;
• interpréter les résultats, en tirer des conclusions, et éventuellement les post-traiter, c'est-à-dire utiliser les résultats pour faire d'autres calculs.

D'autre part, les premières versions de ces logiciels utilisaient un langage de programmation (langage DMAP pour Nastran, MAPDL pour Ansys, Code_Aster). Si les fonctions principales — les plus couramment utilisées — ont été intégrées dans l'interface graphique, certaines opérations ne sont possible qu'en entrant des instructions. Dans le cas d'une utilisation avancée du logiciel, il faut donc s'intéresser au langage de programmation.

L'utilisation d'un logiciel de résolution par la méthode des éléments finis est donc faussement simple, ce qui n'est pas sans poser problème :

[…] l'utilisation de plus en plus fréquente de ce genre de technologie par des personnels non spécialistes ou inadéquatement formés commence à être une source d'inquiétude très sérieuse, compte tenu des enjeux de sécurité sous-jacents. De manière générale, utiliser un logiciel quelconque pour résoudre un problème d'ingéniérie sans en comprendre le fonctionnement est très dangereux […].

Jean-Christophe Cuillière, Introduction à la méthode des éléments finis, Paris, Dunod, 2011 (ISBN 978-2-10-056438-5), p. 3-4

On peut essayer de déterminer des compétences minimales requises en fonction de l'utilisation. Définissons trois niveaux d'utilisation en mécanique (le niveau 2 inclut le niveau 1, le niveau 3 inclut le niveau 2) :

1. Élémentaire :
   • calcul sur une structure de poutres (p. ex. treillis) ;
   • calcul volumique sur une pièce isolée.
2. Avancé : résolution de problèmes linéaires, surfaciques et volumiques impliquant plusieurs pièces en contact.
3. Expert : mécanique non linéaire (contact pouvant adhérer ou glisser, grandes déformations, écrouissage, viscosité, …), problèmes multiphysiques (recouvrant plusieurs domaines d'application à partir de la même géométrie), développement de méthodes de post-traitement automatisée, méthodes d'optimisation de formes, …

Les compétences requises recommandées sont :

1. Élémentaire :
   • statique : connaissance des liaisons mécaniques (degrés de liberté, efforts transmissibles), résolution de problèmes isostatiques (graphique et analytique), représentation des efforts par les torseurs ;
   • résistance des matériaux, théorie des poutres (Euler-Bernoulli) : courbe de traction uniaxiale, courbes des efforts de cohésion (effort normal, effort tranchant, moment fléchissant), contraintes et déformations en sollicitations simples (traction/compression, cisaillement, torsion, flexion simple), utilisation des formulaires de flexion des poutres, contraintes en sollicitations composées (traction/compression-flexion, flexion déviée, flexion-torsion), concentrations de contraintes, matage, flambage ;
   • mécanique des milieux continus, théorie de l'élasticité : principe de Barré de Saint Venant, contraintes équivalentes de von Mises et de Tresca ;
   • dessin assisté par ordinateur : manipulation du logiciel pour simplifier la géométrie (édition des esquisses, des fonctions d'extrusion, de révolution) ;
   • mathématiques : repère (géométrie analytique), similitudes (rotations, translations, homothéties), équations différentielles simples, intégration de fonctions polynomiales ;
   • éléments finis : choix des éléments de maillage, contrôle de la qualité du maillage, contrôle de la cohérence des résultats.
2. Avancé :
   • statique et résistance des matériaux : travaux virtuels, résolution de problèmes hyperstatiques, pression de Hertz, essais mécaniques, théorie des poutres épaisses (Timoshenko), éléments de théorie des plaques et des coques ;
   • mécanique des milieux continus : tenseurs des contraintes et des déformations, contraintes principales, coefficients élastiques, énergie de déformation, critères de ruine ;
   • mathématiques : matrices, changement de repère, intégration par partie, interpolation polynomiale, équations aux dérivées partielles ;
   • éléments finis : connaissance des différents types d'éléments et de leurs degrés de liberté, notions de matrice de rigidité et de points de Gauss.
3. Expert :
   • mathématiques : théorème de Green, méthode de Gauss-Legendre ;
   • éléments finis : formes intégrales forte et faible, connaissance des langages sous-jacents du logiciel utilisé (DMAP, MAPDL, Code-Aster).

Principe de la méthode

Équation aux dérivées partielles

Un certain nombre de problèmes physiques sont décrits par des équations aux dérivées partielles (ÉDP) sur un domaine spatial, un volume. Il s'agit d'une généralisation des équations différentielles aux fonctions de plusieurs variables. Par exemple, si l'on a une fonction de trois variables ƒ(x₁, x₂, x₃), l'équation suivante :

$a_{1}{\frac {\partial f}{\partial x_{1}}}+a_{2}{\frac {\partial f}{\partial x_{2}}}+a_{3}{\frac {\partial f}{\partial x_{3}}}+\mathrm {A} =0$

est une équation aux dérivées partielles. Cette équation est assortie de conditions aux limites : valeurs de la fonction ou de ses dérivées partielles en certains points.

Notons que

• la fonction ƒ peut être une fonction vectorielle,
• l'équation fait souvent intervenir des dérivées secondes ∂²ƒ/∂x²_i ou ∂²ƒ/∂x_ix_j (voire d'ordres plus élevés),
• et que les coefficients a_i et A ne sont pas nécessairement des constantes mais peuvent être des fonctions.

La résolution exacte, analytique, de telles équations devient vite impossible manuellement. Par contre, si l'on découpe le domaine spatial en petites cellules, appelées « éléments finis » (ÉF), on peut résoudre simplement l'ÉDP sur chaque élément.

Démarche générale

La méthode des éléments finis (MÉF) consiste donc à :

• découper le modèle spatial en éléments finis : c'est le maillage ;
• écrire une version simplifiée de l'ÉDP sur chaque élément fini ; notons que les conditions limites d'un élément ne sont pas connues (en particulier pour les éléments à l'intérieur de la pièce), on ne connaît que les conditions globales (les conditions de contact à l'extérieur de la pièce) ;
• rassembler les expressions des ÉDP locales pour appliquer les conditions aux limites du problème.

On retrouve la démarche générale analyse-synthèse (découpage-regroupement).

D'un point de vue pratique, la mise en œuvre de la méthode comporte les étapes suivantes :

1. Analyse du problème : définition de l'objectif du calcul, recherche des zones pouvant poser problème, éventuellement calcul manuel avec un modèle très simplifié pour avoir un ordre de grandeur du résultat.
2. Définition du modèle de calcul : la géométrie du système est dessinée avec un logiciel de dessin assisté par ordinateur (DAO). Si le modèle numérique (= sur ordinateur) existe déjà (par exemple pièce dessinée par un bureau d'étude), il faut simplifier la géométrie afin d'avoir un calcul :
   • plus rapide (modèle « léger ») ;
   • ne présentant pas de singularités : les singularités sont des points particuliers se comportant mal vis-à-vis de l'algorithme de calcul ; en ces point-là, les résultats sont en général éloignés de la réalité.
3. Maillage : découpage du modèle en éléments finis ; il convient de choisir des formes d'éléments adaptés, permettant d'avoir un maillage régulier, et de vérifier la qualité du maillage. Le maillage est fait par l'ordinateur suivant les paramètres définis par l'utilisateur.
4. Calcul, fait par l'ordinateur.
5. Affichage des résultats, vérification de leur cohérence et post-traitement.

Maillage

Pièce volumique maillée par des tétraèdres

Exemples d'éléments finis linéaires

Le maillage (meshing) consiste donc à découper l'espace en petits domaines appelés éléments finis. Dans le cas général, on utilise des éléments volumiques (3D), mais on peut :

• utiliser des éléments surfaciques (2D) si la pièce étudiée est une coque (shell), c'est-à-dire si l'épaisseur est petite devant les autres dimensions (typiquement un rapport 1/10) ;
• utiliser des éléments linéïques (1D) si la pièce étudiée est une poutre (beam), c'est-à-dire si ses dimensions transverses sont petites devant la longueur (typiquement un rapport 1/10).

On utilise plusieurs types d'éléments finis. Dans un premier temps, retenons qu'il y a principalement deux types d'éléments pour les coques et les volumes :

• les éléments carrés (coques) ou cubiques (volumes) : ils permettent un maillage régulier et efficace (on a une bonne précision avec peu d'éléments) ; comme il n'y a pas de raison que les formes soient régulières (en particulier que des côtés/arêtes soient parallèles), on parle plutôt de quadrilatères ou d'hexaèdres ;
• les éléments triangulaires (coques) ou tétraédriques (volumes) : ils permettent de s'adapter aux formes complexes.

Un élément fini est une maille ; il est défini par ses nœuds, c'est-à-dire les sommets pour un polygone ou un polyèdre.

La MÉF consiste à calculer les valeurs de la fonction ƒ — solution des ÉDP — aux nœuds ; on ne recherche pas la valeur en tout point de l'espace, mais uniquement en certains points.

La MÉF est donc une discrétisation du problème.

La fonction ƒ est une fonction vectorielle de dimension d. En chaque nœud, il faut donc déterminer d valeurs. On dit que le nœud a « d degrés de liberté ».

Interpolation des fonctions

Donc, on s'intéresse aux valeurs de la fonction ƒ uniquement aux nœud du maillage. La valeur aux autres points est ensuite déduite par interpolation polynomiale.

Sauf cas exceptionnel, dans une interpolation, plus on s'éloigne des points connus, plus l'écart avec la valeur théorique est important. Donc, plus le maillage est grossier, plus l'erreur commise en dehors des nœuds est grande.

Considérons un élément donné dont les nœuds sont numérotés de 1 à n ; on note ƒ_i la valeur de la fonction au nœud i. Ce sont ces valeurs ƒ_i que l'on veut déterminer par la MÉF.

Au sein de l'élément, on remplace donc la fonction par son interpolation, sous la forme

$f(x,y,z)\simeq \sum _{i=1}^{n}\mathrm {N} _{i}(x,y,z)\cdot f_{i}$

où les N_i sont des fonction qui ne dépendent que de la forme de l'élément fini.

On utilise différents types d'éléments finis ; pour chaque type, on a un élément de référence, qui est l'élément fini dont les coordonnées des nœuds valent 0, 1 ou –1. On peut dire en quelque sorte que l'élément fini réel est l'élément fini de référence déformé et placé dans l'espace. Les fonctions N_i sont tabulées pour les éléments de référence ; on applique ensuite un changement de repère.

La MÉF consiste à remplacer la fonction recherchée par son interpolation entre les nœuds. Les fonctions d'interpolation pour les éléments de référence sont connues.

On écrit souvent cette opération d'interpolation sous la forme d'un produit matriciel :

sur un élément fini donné, $f(x,y,z)\simeq {\begin{pmatrix}\mathrm {N} _{1}(x,y,z)&\mathrm {N} _{2}(x,y,z)&\ldots &\mathrm {N} _{n}(x,y,z)\end{pmatrix}}\cdot {\begin{pmatrix}f_{1}\\f_{2}\\\vdots \\f_{n}\end{pmatrix}}$.

Par la suite, nous notons simplement x le vecteur position (x, y, z). L'équation ci-dessus devient donc :

$f(x)\simeq {\begin{pmatrix}\mathrm {N} _{1}(x)&\mathrm {N} _{2}(x)&\ldots &\mathrm {N} _{n}(x)\end{pmatrix}}\cdot {\begin{pmatrix}f_{1}\\f_{2}\\\vdots \\f_{n}\end{pmatrix}}$.

Formulations intégrales

Nous avons m éléments/mailles, donc m équations. Le système d'ÉDP peut s'écrire

$\left\{{\begin{aligned}\mathrm {R} _{1}(f)=0\\\mathrm {R} _{2}(f)=0\\\vdots \\\mathrm {R} _{m}(f)=0\\\end{aligned}}\right.$

la fonction R_i étant appelée « résidu » de l'élément i. Pour une fonction ψ_i arbitraire, dite « fonction de pondération », on a ψ_i × R_i(ƒ) = 0 et donc a fortiori

$\mathrm {W} _{i}=\iiint _{\mathcal {V}}\psi _{i}\mathrm {R} _{i}(f)\mathrm {dV} =0$.

Ceci constitue la forme intégrale forte. Les fonctions sont notées W par analogie avec les travaux virtuels. On choisit les fonctions de pondération de manière judicieuse afin de faciliter la résolution ; ce ne sont pas nécessairement des fonctions connues, elles peuvent contenir la fonction ƒ.

Les ÉDP sont en général d'ordre 2. Les fonctions résidus contiennent donc des termes du type ∂²ƒ/∂x². À une dimension, on peut faire une intégration par parties :

$\int _{x_{1}}^{x_{2}}\psi _{i}{\frac {\mathrm {d} ^{2}f}{\mathrm {d} x^{2}}}\mathrm {d} x=-\int _{x_{1}}^{x_{2}}{\frac {\mathrm {d} \psi _{i}}{\mathrm {d} x}}\times {\frac {\mathrm {d} f}{\mathrm {d} x}}\mathrm {d} x+\left[\psi _{i}{\frac {\mathrm {d} f}{\mathrm {d} x}}\right]_{x_{1}}^{x_{2}}$.

Cette opération permet :

• de faire apparaître les conditions aux limites ƒ’(x₁) et ƒ’(x₂) ;
• de diminuer l'ordre de dérivation de 1.

Typiquement, si l'ÉDP est de la forme

${\frac {\mathrm {d} ^{2}f}{\mathrm {d} x^{2}}}+q(x)=0$

q étant une fonction, alors la fonction résidu est de la forme

$\mathrm {R} _{i}(f)={\frac {\mathrm {d} ^{2}f}{\mathrm {d} x^{2}}}+q_{i}(x)$.

Nous avons donc :

$\mathrm {W} _{i}=\int _{x_{1}}^{x_{2}}\psi _{i}\left({\frac {\mathrm {d} ^{2}f}{\mathrm {d} x^{2}}}+q_{i}(x)\right)\mathrm {d} x=\int _{x_{1}}^{x_{2}}\psi _{i}{\frac {\mathrm {d} ^{2}f}{\mathrm {d} x^{2}}}\mathrm {d} x+\int _{x_{1}}^{x_{2}}\psi _{i}q_{i}(x)\mathrm {d} x$

et donc

$\mathrm {W} _{i}=-\int _{x_{1}}^{x_{2}}{\frac {\mathrm {d} \psi _{i}}{\mathrm {d} x}}\times {\frac {\mathrm {d} f}{\mathrm {d} x}}\mathrm {d} x+\left[\psi _{i}{\frac {\mathrm {d} f}{\mathrm {d} x}}\right]_{x_{1}}^{x_{2}}+\int _{x_{1}}^{x_{2}}\psi _{i}q_{i}(x)\mathrm {d} x$.

Pour des intégrales de volume, l'équivalent de l'intégration par partie est le théorème de Green :

$\iiint _{\mathcal {V}}\psi _{i}{\frac {\partial ^{2}f}{\partial x^{2}}}\mathrm {dV} =-\iiint _{\mathcal {V}}{\frac {\partial \psi _{i}}{\partial x}}\times {\frac {\partial f}{\partial x}}\mathrm {dV} +\iint _{\mathcal {S}}\psi _{i}{\frac {\partial f}{\partial x}}n_{x}\mathrm {dS}$

où

• ${\mathcal {V}}$ est le domaine spatial considéré ;
• ${\mathcal {S}}$ est la frontière de ce domaine ;
• n_x est la composante selon x du vecteur ${\vec {n}}$ normal à ${\mathcal {S}}$ et pointant vers l'extérieur.

Les conditions limites apparaissent dans l'intégrale sur la surface.

Lorsque l'on applique le théorème de Green aux fonctions W_i, on obtient la « forme intégrale faible ».

La forme intégrale faible permet de diminuer d'un degré l'ordre de dérivation et fait apparaître les conditions aux limites. Le choix de fonctions de pondération ψ_i judicieuses facilite la résolution.

Matrice de rigidité

La méthode consiste donc à résoudre le système d'équations

$\left\{{\begin{aligned}\mathrm {W} _{1}=0\\\vdots \\\mathrm {W} _{m}=0\\\end{aligned}}\right.$

les fonctions W_i étant les formulations intégrales faibles des fonctions résidu pondérées (les fonctions résidu étant les fonctions originales du système d'ÉDP).

Sur un élément fini i donné, on ne travaille pas sur la fonction ƒ elle-même, mais sur son interpolation linéaire (pour le nœud k de l'élément i)

$f(x)\simeq \sum _{k}\mathrm {N} _{i,k}(x)\cdot f_{i,k}$.

La formulation intégrale faible fait de ce fait intervenir les intégrales des fonctions N_{i, k}(x) qui sont connues. Ainsi, la fonction W_i(x) est une fonction faisant intervenir (pour simplifier, on reprend la forme à une dimension) :

$\mathrm {W} _{i}=\sum _{k}f_{i,k}\left(-\int _{x_{1}}^{x_{2}}{\frac {\mathrm {d} \psi _{i}(x)}{\mathrm {d} x}}\times {\frac {\mathrm {d} \mathrm {N} _{i,k}(x)}{\mathrm {d} x}}\mathrm {d} x+\left[\psi _{i}(x){\frac {\mathrm {d} \mathrm {N} _{i,k}(x)}{\mathrm {d} x}}\right]_{x_{1}}^{x_{2}}\right)+\int _{x_{1}}^{x_{2}}\psi _{i}q_{i}(x)\mathrm {d} x$.

On obtient ainsi une équation matricielle de la forme :

F = K·U

où

• U est le vecteur contenant les valeurs de la fonction solution du système d'équations (W_i = 0) aux nœuds de l'élément (U n'est pas la fonction ƒ, en raison des fonctions de pondération choisies) ;
• F est le vecteur contenant les valeurs représentant les actions extérieures sur l'élément fini (action des éléments voisins, actions à distance, et conditions limites si l'élément est à la surface) ;
• K représente la loi de comportement de l'élément fini ; ses composantes sont connues et se calculent à partir de la forme de l'élément fini et des propriétés de l'espace (propriétés du matériau si c'est un milieu matériel).

Si l'élément fini a n degrés de liberté, alors l'équation s'écrit avec les composantes :

${\begin{pmatrix}\mathrm {F} _{1}\\\vdots \\\mathrm {F} _{n}\\\end{pmatrix}}={\begin{pmatrix}\mathrm {K} _{11}&\ldots &\mathrm {K} _{1n}\\\vdots &\ddots &\vdots \\\mathrm {K} _{n1}&\ldots &\mathrm {K} _{nn}\\\end{pmatrix}}\cdot {\begin{pmatrix}\mathrm {U} _{1}\\\ \vdots \\\mathrm {U} _{n}\end{pmatrix}}$

La méthode des éléments finis a d'abord été développée pour la résistance des matériaux. Dans ce domaine, F_i est une composante d'une force qui s'exerce sur un nœud i et U_i est une composante du déplacement d'un nœud. On a un rapport évident avec la loi des ressorts

F = k·Δl

où k est l raideur du ressort. De fait, la matrice K est appelée « matrice de rigidité. »

Comme indiqué précédemment, on ne connaît pas les conditions limites pour un élément fini quelconque. On connaît donc la matrice K, qui est caractéristique de l'ÉF, mais rien sur F ni sur U.

Lors de la synthèse, les vecteurs F locaux sont mis bout-à-bout pour former un vecteur global F_g, contenant toutes les composantes de tous les ÉF. De même, les vecteurs locaux U viennent former un grand vecteur global U_g, et les matrices de rigidité locales viennent former une matrice globale K_g. Cette opération est appelée « expansion-assemblage ».

On a ainsi un système linéaire de la forme

F_g = K_g·U_g

qui contient les conditions aux limites (valeurs imposées de F_g ou U_g).

Habituellement, dans une équation matricielle comme celle-ci, soit l'on connaît entièrement U_g et l'on cherche F_g, soit l'on connaît entièrement F_g et l'on cherche U_g. Ici, c'est différent, puisque certaines composantes de U_g sont connues et d'autres pas ; et certaines composantes de F_g sont connues et d'autres sont inconnues.

Par exemple, en résistance des matériaux, au niveau des appuis, les déplacements U sont imposés (égaux à 0), et les actions des appuis F (forces) sont inconnues. Aux points d'application des forces de contact, les actions F sont connues mais pas les déplacements U.

Mais on se retrouve globalement à résoudre un système linéaire.

La MÉF fait donc une linéarisation du problème.

Exemples d'application

Les quatre principaux domaines pour lesquels on a recours à la MÉF sont :

• la résistance des matériaux ;
• le transfert de chaleur;
• la mécanique des fluides (CFD : computational fluid dynamics) ;
• l'électromagnétisme.

Résistance des matériaux

Équilibre mécanique d'un élément de matière

Historiquement, c'est le domaine pour lequel la MÉF a été inventée.

Le système d'équations aux dérivées partielles est fourni par les équations de la statique (les efforts extérieurs s'exerçant sur un élément de matière s'annulent) :

$\left\{{\begin{aligned}{\frac {\partial \sigma _{x}}{\partial x}}+{\frac {\partial \tau _{xy}}{\partial y}}+{\frac {\partial \tau _{xz}}{\partial z}}+\mathrm {F} _{\mathrm {V} x}=0\\{\frac {\partial \tau _{xy}}{\partial x}}+{\frac {\partial \sigma _{y}}{\partial y}}+{\frac {\partial \tau _{yz}}{\partial z}}+\mathrm {F} _{\mathrm {V} y}=0\\{\frac {\partial \tau _{xz}}{\partial x}}+{\frac {\partial \tau _{yz}}{\partial y}}+{\frac {\partial \sigma _{z}}{\partial z}}+\mathrm {F} _{\mathrm {V} z}=0\\\end{aligned}}\right.$

où

• σ_i est la contrainte normale s'appliquant à la face normale à l'axe i ;
• τ_ij (i ≠ j ) est la contrainte de cisaillement s'appliquant sur la face normale à l'axe j et orientée selon l'axe i ; on a τ_ij = τ_ji ;
• F_Vi est la composante selon l'axe i de la densité de force volumique, en général le poids :
  F_Vx = F_Vy = 0, F_Vz = ρg
  ρ étant la masse volumique du matériau et g l'accélération de la gravité.

Les conditions aux limites sont de deux ordres :

• déplacements imposés : en général, on impose que le déplacement de certaines surfaces soit nul ; ce sont les endroits où la pièce étudiée est fixée à son support (appuis) ;
• forces surfaciques imposées : ce sont les efforts extérieurs.

Le déplacement d'un point sous l'effet de la charge est représenté par le vecteur

${\vec {\mathrm {U} }}{\begin{pmatrix}u_{x}\\u_{y}\\u_{z}\end{pmatrix}}$

On impose donc un champ de déplacement ${\vec {\mathrm {U} }}$ sur une surface S_u. On relie ce déplacement aux contraintes (grandeurs de l'ÉDP) par le passage aux déformations :

$\varepsilon _{ij}={\frac {1}{2}}\left({\frac {\partial u_{i}}{\partial j}}+{\frac {\partial u_{j}}{\partial i}}\right)$, par exemple $\varepsilon _{xx}={\frac {\partial u_{x}}{\partial x}}$ et $\varepsilon _{xy}={\frac {1}{2}}\left({\frac {\partial u_{x}}{\partial y}}+{\frac {\partial u_{y}}{\partial x}}\right)$

et l'utilisation de la loi de Hooke :

σ = H·ε

où H est une matrice de coefficients d'élasticité, caractéristique du matériau ; σ et ε sont des vecteurs dont les composantes sont respectivement les contraintes et les déformations.

Équilibre d'un élément de matière ayant une face libre sur laquelle s'exerce une force surfacique

On impose des efforts surfaciques ${\vec {f}}_{\mathrm {S} }{\begin{pmatrix}f_{\mathrm {S} x}\\f_{\mathrm {S} y}\\f_{\mathrm {S} z}\end{pmatrix}}$ (des pressions) sur une surface S_f. Si le vecteur normal unitaire à la surface est notée

${\vec {n}}{\begin{pmatrix}n_{x}\\n_{y}\\n_{z}\end{pmatrix}}$

alors l'équilibre de la surface impose

$\left\{{\begin{aligned}\sigma _{x}\cdot n_{x}+\tau _{xy}\cdot n_{y}+\tau _{xz}\cdot n_{z}=f_{\mathrm {S} x}\\\tau _{xy}\cdot n_{x}+\sigma _{y}\cdot n_{y}+\ +\tau _{yz}\cdot n_{z}=f_{\mathrm {S} y}\\\tau _{xz}\cdot n_{x}+\tau _{yz}\cdot n_{y}+\sigma _{z}\cdot n_{z}=f_{\mathrm {S} z}\\\end{aligned}}\right.$

Notons bien que

la solution de la méthode des éléments finis est le champ de déplacement ${\vec {\mathrm {U} }}$.

C'est ce que calcule l'ordinateur. Ce champ est transformé en champ de déformation ε (avec la définition fournie ci-dessus), puis le champ des contraintes σ est calculé avec la loi de Hooke.

• 
  Étude de la flexion d'une poutre encastrée
• 
  Champ de déplacement calculé
• 
  Champ de contraintes calculé

Transfert de chaleur

C'est le domaine pour lequel les équations sont les plus simples.

L'équation aux dérivées partielles est l'équation de la chaleur (équation de Fourier) :

${\frac {\partial \mathrm {T} }{\partial t}}=\mathrm {D} \cdot \left({\frac {\partial ^{2}\mathrm {T} }{\partial x^{2}}}+{\frac {\partial ^{2}\mathrm {T} }{\partial y^{2}}}+{\frac {\partial ^{2}\mathrm {T} }{\partial z^{2}}}\right)+{\frac {\mathrm {P} }{\rho c}}$

où

• T est la température ;
• t est le temps ;
• D est le coefficient de diffusivité thermique caractéristique du matériau ;
• P est un terme de production de chaleur volumique ;
• ρ est la masse volumique du matériau ;
• c est la chaleur massique du matériau.

En régime stationaire, on a

$\mathrm {D} \cdot \left({\frac {\partial ^{2}\mathrm {T} }{\partial x^{2}}}+{\frac {\partial ^{2}\mathrm {T} }{\partial y^{2}}}+{\frac {\partial ^{2}\mathrm {T} }{\partial z^{2}}}\right)+{\frac {\mathrm {P} }{\rho c}}=0$

Les conditions limite consistent à imposer :

• des températures sur certaines surfaces ;
• des sources chaleur ou des puits (pertes de chaleur par convection par exemple) sur d'autres surfaces.

Mécanique des fluides

Étude de la chute d'une goutte par méthode du volume de fluide (VoF).

La simulation numérique des écoulement est en générale appelée MFN (mécanique des fluides numérique) ou bien CFD (computational fluids dynamics).

On utilise parfois la méthode des éléments finis, mais la plupart du temps, c'est la méthode des volumes finis qui est utilisée : les équations différentielles ne sont pas résolues de la même manière. Dans le cas des logiciels munis d'une interface graphique, on voit peu de différence avec un logiciel d'éléments finis, mais dans tous les cas, le traitement de la géométrie du système est très différent : on ne peut pas simplifier les détails de la même manière, et la stratégie de maillage obéit à des impératifs différents.

Algorithmes de résolution

La méthode étant posée, les logiciels ont recours à des méthodes numériques — calculs approchés — pour effectuer les calculs. Nous présentons ici quelques méthodes utilisées classiquement.

Un modèle d'éléments finis peut contenir dix mille, cent mille voire parfois plus d'un million d'éléments. Le nombre de calcul est très important, la résolution est donc extrêmement gourmande en ressources et en temps : plusieurs heures voire plusieurs jours de calcul, et peut nécessiter le lancement d'un calcul parallélisé sur un superordinateur. L'utilisation d'algorithmes de calcul performants n'est donc pas inhérent à la méthode des éléments finis, mais est toutefois indispensable pour des raisons de performance.

Intégration

Exemple de points de Gauss pour un domaine d'intégration carré.

De manière intuitive, pour résoudre une équation différentielle, il faut inverser le processus de dérivation c'est-à-dire calculer une intégrale. Dans le MÉF, les intégrales sur les éléments finis sont calculées de manière approchée, numérique. On utilise en général la méthode de Gauss-Legendre :

Une fonction ƒ étant définie sur l'intervalle [-1 ; 1], on peut approcher son intégrale

$\int _{-1}^{1}f(x)\mathrm {d} x\simeq \sum _{i=1}^{n}\omega _{i}\cdot f(x_{i})$.

Les points x_i sont appelés « points de Gauss ». Les coefficients ω_i sont des constantes. Les valeurs des x_i et des ω_i dépendent du nombre de points n que l'on considère. Plus ce nombre est élevé, plus le calcul est précis, mais plus il est long.

Cette méthode se généralise aux fonctions de plusieurs variables. La figure ci-contre montre un exemple de carré avec quatre points de Gauss (il y a plusieurs possibilité) ; les coefficients de pondération ω_i sont ici tous égaux à 1. Soit ƒ(x, y ) une fonction définie sur [-1 ; 1]² ; on peut approcher son intégrale par

${\begin{aligned}\int _{-1}^{1}\int _{-1}^{1}f(x,y)\mathrm {d} x\mathrm {d} y\simeq \ &1\times f(1/{\sqrt {3}},1/{\sqrt {3}})\\&+1\times f(-1/{\sqrt {3}},1/{\sqrt {3}})\\&+1\times f(-1/{\sqrt {3}},-1/{\sqrt {3}})\\&+1\times f(1/{\sqrt {3}},-1/{\sqrt {3}})\\\end{aligned}}$

ce calcul est exact si ƒ est un polynôme de dergé inférieur ou égal à trois, c'est-à-dire comportant des termes de la forme xⁱ×y^j avec i + j ≤ 3.

Ceci se généralise pour d'autres formes (en particulier triangle) et pour des intégrations sur des volumes (cubes, tétraèdres, prismes droits, …).

Bien évidemment, il n'y a aucune raison pour qu'un élément fini ait ses frontières aux coordonnées 1 ou -1. On fait donc un changement de repère pour relier l'élément fini du maillage à un élément fini de référence. Ce changement d'espace de coordonnées dans le contexte dérivée/intégrale fait classiquement intervenir les matrices jacobiennes.

Rappelons que l'on ne travaille pas directement avec la fonction ƒ mais avec son interpolation linéaire ∑_{j = 1}^m N_j·ƒ_j (si l'élément fini a m nœuds). Le calcul de l'intégrale — à une dimension x — devient donc

$\int _{-1}^{1}f(x)\mathrm {d} x\simeq \sum _{i=1}^{n}\omega _{i}\cdot \left(\sum _{j=1}^{m}\mathrm {N} _{j}(x_{i})\cdot f_{j}\right)=\sum _{i=1}^{n}\sum _{j=1}^{m}\omega _{i}\cdot \mathrm {N} _{j}(x_{i})\cdot f_{j}$.

Rappelons que tout comme les x_i et ω_i, les N_j sont tabulées pour les éléments de référence. On peut faire ressortir les inconnues :

$\int _{-1}^{1}f(x)\mathrm {d} x\simeq \sum _{j=1}^{m}f_{j}\cdot \underbrace {\left(\sum _{i=1}^{n}\omega _{i}\cdot \mathrm {N} _{j}(x_{i})\right)} _{\alpha _{j}}$

les coefficients α_j étant connus.

Notes et références

---

Introduction < > Rappels de mécanique