Hypothèses

On considère ƒ une fonction continue sur Ω de carré sommable et u la solution de l'équation aux dérivées partielles suivante sur Ω (Δ est l'opérateur laplacien) :

${\displaystyle \displaystyle -\Delta u+k^{2}u=f}$

Avec la condition au bord u = 0 sur δΩ. Ceci peut également se réécrire u ∈ V₀. Cette condition au bord s'appelle la condition de Dirichlet.

On démontre qu'il existe une solution unique à ce problème d'EDP à l'aide du théorème de Lax-Milgram.

Formulation faible

Soit v ∈ V₀ quelconque. Multiplions les deux parties de l'équation précédente par v puis sommons sur le domaine Ω, puisque v et ƒ sont tous deux de carré sommable sur ce domaine. On obtient l'équation :

${\displaystyle -\int _{\Omega }{v\Delta u\,\mathrm {d} \omega }+k^{2}\int _{\Omega }{vu\,\mathrm {d} \omega }=\int _{\Omega }{vf\,\mathrm {d} \omega }\,}$

On utilise pour le premier terme une intégration par parties^[3] :

${\displaystyle -\int _{\Omega }{v\Delta u\,\mathrm {d} \omega }=-\int _{\partial \Omega }{{\frac {\partial u}{\partial n}}v\,\mathrm {d} s}+\int _{\Omega }{(\nabla u\cdot \nabla v)\,\mathrm {d} \omega }}$

Dans cette formulation, v est nulle sur le bord (v ∈ V₀) ce qui permet d'obtenir la formulation faible du problème :

${\displaystyle \int _{\Omega }{\nabla u\cdot \nabla v\,\mathrm {d} \omega }+k^{2}\int _{\Omega }{vu\,\mathrm {d} \omega }=\int _{\Omega }{vf\,\mathrm {d} \omega }}$

Si u est deux fois différentiable, il y a équivalence entre cette formulation et celle du problème initial donné dans la section hypothèse et alors la solution de la formulation faible est la même que la solution initiale. On peut donc résoudre la formulation faible au lieu de résoudre le problème initial.

La question de savoir s'il y a équivalence entre la formulation faible et la formulation initiale donnée dans les hypothèses peut être particulièrement délicate dans les cas limites où l'ouvert Ω n'est pas suffisamment régulier (par exemple s'il y a des points singuliers) ou si ƒ n'est pas suffisamment dérivable (si l'on ne suppose pas que ƒ est au moins L²). Il faut alors souvent se ramener à une étude au cas par cas et rien ne dit que la formulation faible aura les mêmes solutions que l'équation de départ. Dans la majorité des problèmes physiques la solution est souvent ${\displaystyle {\mathcal {C}}^{\infty }}$ et l'on ne se pose pas de tels problèmes. Néanmoins, pour des domaines avec des points singuliers cette équivalence peut poser problème. Ceci peut être gênant pour l'étude de fissures en mécanique des milieux continus par exemple.

Choix d'un maillage

Un exemple de maillage triangulaire.

Raffinement de maillage par diminution de la taille d'éléments (type h, à gauche) ou par augmentation du degré des éléments (type p, à droite).

La méthode des éléments finis repose sur un découpage de l'espace selon un maillage. D'habitude l'on choisit un maillage carré ou triangulaire mais rien n'interdit de choisir des maillages plus complexes. Il n'est pas non plus nécessaire que le maillage soit régulier et l'on a tendance à resserrer le maillage près des endroits d'intérêt (par exemple aux endroits où l'on pense que la solution va beaucoup varier) ; cependant, il faut veiller à avoir des éléments faiblement distordus (se rapprocher d'un polygone régulier). Plus ce maillage est resserré, plus la solution que l'on obtient par la méthode des éléments finis sera précise et proche de la « vraie » solution de l'équation aux dérivés partielles.

On appelle traditionnellement h la plus grande dimension d'un élément (le diamètre de la sphère dans laquelle s'inscrit l'élément), et p le degré du polynôme décrivant le côté ou l'arête (p = 1 pour des côtés/arêtes droits, p = 2 pour des côtés/arête présentant une courbure).

Dans l'idéal, le maillage doit donc épouser les contours δΩ du domaine. Si δΩ est courbe, alors on peut :

• soit utiliser des éléments plus petits, on parle de raffinement de type h ;
• soit utiliser des éléments dont les côtés (en 2D) ou arêtes (en 3D) sont courbes, on parle de raffinement de type p.

Lorsque le degré des polynômes p est élevé, on parle de méthode des éléments spectraux.

Définition formelle

Le cadre mathématique des éléments finis permet de poser un système assurant la bonne définition de la solution.

Cette définition offre a priori beaucoup de liberté, mais en général, plusieurs conditions sont imposées : le domaine K est pris comme non dégénéré, l'espace des fonctions de base sera choisi de dimension finie et simples à calculer et les degrés de liberté seront pris comme vérifiant la propriété d'unisolvance :

Fonctions de base

Fonction de base en dimension 1. Les x_i sont les nœuds du réseau.

On doit après prendre une base de fonctions « adaptées » au maillage. Plusieurs choix sont alors possibles. En général, les fonctions de base utilisées pour les éléments finis sont interpolantes, c'est-à-dire que les valeurs nodales sont les valeurs des grandeurs inconnues aux nœuds.

La plus simple est l'emploi des polynômes de Lagrange. Dans cette méthode les fonctions de base valent 1 à un nœud du maillage et 0 à tous les autres. La fonction de base i est alors la fonction valant 1 au nœud i et 0 sur les autres nœuds et qui est polynomiale sur chaque élément. Un exemple de telles fonctions est représenté en dimension 1 à côté. Il y a autant de fonctions de base par élément que de nombre de nœuds.

On appelle « élément » la donnée d'une géométrie (souvent polygonale en 2D, polyédrique en 3D) et de fonctions de base associées à cette géométrie.

D'autres solutions peuvent exister pour les fonctions de base. On cite ici un seul exemple les éléments finis d'Hermite qui ont la particularité d'avoir deux fonctions de base associées à chaque nœud. Dans cette version, la valeur de la solution est ajustée avec la première fonction alors que la deuxième permet d'ajuster la valeur de la dérivée. Ce type de fonctions de base peut avoir un intérêt pour la résolution de certaines équations aux dérivées partielles (par exemple l'équation des plaques en mécanique des milieux continus), même si elle nécessite d'avoir deux fois plus de fonctions pour un maillage donné.

Quelques éléments classiques

Principaux types d'éléments utilisés en 3D.

En 2D :

• triangles :
  • triangles de degré 1 (triangles à 3 nœuds, fonctions linéaires),
  • triangles de degré 2 (triangles à 6 nœuds, polynômes de degré 2) ;
• quadrilatères :
  • quadrilatères de degré 1 (carrés à quatre nœuds, fonctions linéaires),
  • quadrilatères de degré 2 (carrés à 8 ou 9 nœuds, polynômes de degré 2).

En 3D :

• tétraèdres :
  • tétraèdres de degré 1 (quatre nœuds, fonctions linéaires)^[4],
  • tétraèdres de degré 2 (dix nœuds, polynômes de degré 2)^[5] ;
• hexaèdres^[6] :
  • hexaèdres de degré 1 (huit nœuds, fonctions linéaires),
  • hexaèdres de degré 2 (vingt nœuds, polynômes de degré 2),
  • hexaèdres triquadratiques (vingt-sept nœuds, polynômes de degré 2).

Les deux fonctions de base d'Hermite associées au nœud 0 avec les nœuds voisins en +1 et -1 en dimension 1.

Condition de Neumann

La condition qui suit est très différente de celle de Dirichlet. On pose comme condition au bord que la dérivée normale existe sur le bord ${\displaystyle {\bar {\Omega }}}$, et que la condition de Neumann ${\displaystyle {\frac {\partial u}{\partial n}}({\overline {\Omega }})=0}$ est vérifiée.

Si la fonction est supposée différentiable au bord ${\displaystyle D_{u}(X)\cdot n=0,\forall X\in {\overline {\Omega }}}$, voire si elle admet un gradient ${\displaystyle \langle \nabla u(X)|n\rangle =0,\forall X\in {\overline {\Omega }}}$.

Le résultat fonctionne de la même manière car l'élément clef de la démonstration où intervient l'hypothèse de bord est que l'intégrale ${\displaystyle \int _{\partial \Omega }{{\frac {\partial u}{\partial n}}v\,\mathrm {d} \omega }=0}$ parce que cette fois-ci ce n'est pas la fonction test mais la dérivée normale qui est nulle.

Par la suite, la différence réside surtout dans le choix des vecteurs de base pour la discrétisation : il faut conserver les fonctions tests propres aux nœuds du bord.

Définitions et notations

On cherche ici à déterminer le vecteur des déplacements ${\displaystyle \left\{U\right\}.\,}$ C'est un vecteur dont chaque composante est également appelée degré de liberté (ddl) :

• 3 ddl de translation : ${\displaystyle U_{x},U_{y},U_{z}\,}$ ;
• 3 ddl de rotation : ${\displaystyle \theta _{x},\theta _{y},\theta _{z}\,}$.

On écrit alors le tenseur des déformations ${\displaystyle \left[\epsilon \right]\,}$, qui modélise la façon dont le matériau va se déformer par rapport à sa position initiale.

Sous l'hypothèse des petites déformations, on a ${\displaystyle (A'B')^{2}=\left(dx+{\partial u \over \partial x}\ dx\right)^{2}+\left({\partial v \over \partial x}\ dx\right)^{2}\,}$

${\displaystyle \epsilon _{x}={{A'B'-AB} \over AB}\,}$

Comme ${\displaystyle AB=dx\,}$, on a ${\displaystyle A'B'=(1+\epsilon _{x})\,dx\,}$

${\displaystyle \Rightarrow 2\epsilon _{x}+\epsilon _{x}^{2}=2{\partial u \over \partial x}+\left({\partial u \over \partial x}\right)^{2}+\left({\partial v \over \partial x}\right)^{2}\,}$

Sous l'hypothèse des petites déformations, on néglige les termes d'ordre 2 :

${\displaystyle \epsilon _{x}={\partial u \over \partial x}\,}$

Remarque : ${\displaystyle \epsilon \,}$ est sans dimension.

On écrit ensuite le tenseur des contraintes) ${\displaystyle \left[\sigma \right]\,}$, qui représente les efforts internes qui s'appliquent dans la structure.

• ${\displaystyle \sigma _{11},\sigma _{22},\sigma _{33}\,}$ : contraintes normales.
• ${\displaystyle \sigma _{12},\sigma _{13},\sigma _{23}\,}$ : contraintes de cisaillement.

Une contrainte est homogène à une pression (N/m²).

Ces deux tenseurs sont liés par une loi de comportement. On considère la loi de Hooke :

${\displaystyle {\begin{Bmatrix}\sigma _{11}\\\sigma _{22}\\\sigma _{33}\\\sigma _{12}\\\sigma _{23}\\\sigma _{31}\end{Bmatrix}}={{E(1-\nu )} \over {(1+\nu )(1-2\nu )}}{\begin{bmatrix}1&{\nu \over {1-\nu }}&{\nu \over {1-\nu }}&0&0&0\\&1&{\nu \over {1-\nu }}&0&0&0\\&&1&0&0&0\\&&&{(1-2\nu ) \over {2(1-\nu )}}&0&0\\&{\mbox{(sym)}}&&&{(1-2\nu ) \over {2(1-\nu )}}&0\\&&&&&{(1-2\nu ) \over {2(1-\nu )}}\end{bmatrix}}{\begin{Bmatrix}\epsilon _{11}\\\epsilon _{22}\\\epsilon _{33}\\\gamma _{12}\\\gamma _{23}\\\gamma _{31}\end{Bmatrix}}\,}$

avec

• ${\displaystyle E\,}$: module de Young (N/m²)
• ${\displaystyle \nu \,}$ : coefficient de Poisson (sans dimension)

On utilise parfois le module de cisaillement : ${\displaystyle G={E \over {2(1+\nu )}}\,}$

Pour un matériau isotrope, il n'y a que deux paramètres indépendants. Il y en a 6 pour un matériau isotrope transverse, 9 pour un matériau orthotrope et 21 pour un matériau anisotrope.

En notation matricielle, on écrit : ${\displaystyle \left\{\sigma \right\}=\left[D\right]\left\{\epsilon \right\}\,}$

${\displaystyle \left[D\right]\,}$ est appelée matrice d'élasticité du matériau.

L'énergie de déformation W s'écrit

${\displaystyle W={1 \over 2}\int _{v}\sigma .\epsilon .\,\mathrm {d} v\,}$

Travail d'une force

Le travail d'une force est le produit de la force par le déplacement de son point d'application : ${\displaystyle \tau _{F}={1 \over 2}F.U_{F}\,}$

Moment

Le moment d'une force est la mesure d'une action entraînant une rotation, proportionnelle à la force et à la longueur du bras de levier entre le centre de rotation du solide et le point d'application de la force.

Équations fondamentales

Équations d'équilibre local

${\displaystyle {\partial \sigma _{x} \over \partial x}+{\partial \sigma _{xy} \over \partial y}+{\partial \sigma _{xz} \over \partial z}+F_{x}=0\,}$

${\displaystyle {\partial \sigma _{y} \over \partial y}+{\partial \sigma _{xy} \over \partial x}+{\partial \sigma _{yz} \over \partial z}+F_{y}=0\,}$

${\displaystyle {\partial \sigma _{z} \over \partial z}+{\partial \sigma _{zx} \over \partial x}+{\partial \sigma _{zy} \over \partial y}+F_{z}=0\,}$

Relations déformations-déplacements

${\displaystyle {\epsilon _{x}={\partial u \over \partial x}}\qquad {\epsilon _{y}={\partial v \over \partial y}}\qquad {\epsilon _{z}={\partial w \over \partial z}}\,}$

${\displaystyle \gamma _{xz}={\partial w \over \partial x}+{\partial u \over \partial z}\,}$

${\displaystyle \gamma _{xy}={\partial u \over \partial y}+{\partial v \over \partial x}\,}$

${\displaystyle \gamma _{yz}={\partial w \over \partial y}+{\partial v \over \partial z}\,}$

Symboliquement, on écrit ${\displaystyle \{\epsilon \}=[S]\{U\}\,}$

Signification du coefficient de Poisson.

Si on applique au barreau une contrainte ${\displaystyle \sigma _{x}\,}$, on observe un rétrécissement dans la direction y correspondant à une déformation ${\displaystyle -\nu \sigma _{x} \over E\,}$

Exemple de formulation : barre en traction

On suppose que le déplacement en tout point de la barre est donné par un polynôme du 1^er degré : ${\displaystyle u(x)=a_{1}+a_{2}x\,}$

On a ${\displaystyle u(0)=u_{1}\,}$ et ${\displaystyle u(L)=u_{2}\,}$

d'où ${\displaystyle u(x)=u_{1}\left(1-{x \over L}\right)+{x \over L}u_{2}\,}$

qu'on écrit symboliquement : ${\displaystyle u(x)=[N(x)]\{U\}\,}$ avec ${\displaystyle \left|{\begin{matrix}[N(x)]&=&\left[\left(1-{x \over L}\right);\left({x \over L}\right)\right]\\\{U\}&=&{\begin{Bmatrix}u_{1}\\u_{2}\end{Bmatrix}}\end{matrix}}\right.\,}$

On en déduit :

${\displaystyle \epsilon =\left[-{1 \over L};{1 \over L}\right]{\begin{Bmatrix}u_{1}\\u_{2}\end{Bmatrix}}=[B]\{U\}\,}$

${\displaystyle \sigma =E\left[-{1 \over L};{1 \over L}\right]{\begin{Bmatrix}u_{1}\\u_{2}\end{Bmatrix}}=[D]\epsilon \,}$

D'autre part, on a par définition :

${\displaystyle {\begin{Bmatrix}F_{1}\\F_{2}\end{Bmatrix}}={\begin{Bmatrix}-\sigma S\\\sigma S\end{Bmatrix}}\,}$ où S est l'aire de la section de la barre.

On pose :

${\displaystyle \{F\}={\begin{Bmatrix}F_{1}\\F_{2}\end{Bmatrix}}=[A]\sigma \qquad {\mbox{avec}}\qquad \{A\}=S{\begin{Bmatrix}-1\\1\end{Bmatrix}}\,}$

On obtient finalement :

${\displaystyle \{F\}=[A][D][B]\{U\}\,}$

Soit une relation du type :

${\displaystyle \{F\}=[K]\{U\}\qquad {\mbox{avec}}\qquad [K]=[A][D][B]\,}$

En explicitant :

${\displaystyle [K]={ES \over L}{\begin{bmatrix}1&-1\\-1&1\end{bmatrix}}\,}$

On voit que la matrice de rigidité se calcule comme le produit de trois matrices :

${\displaystyle [B]\,}$ : Transformation des déplacements aux déformations

${\displaystyle [D]\,}$ : Matrice d'élasticité du matériau

${\displaystyle [A]\,}$ : Transformation des contraintes en forces

Formulation générale (méthode directe)

La démarche est la suivante :

1. On exprime le déplacement ${\displaystyle \{u(x)\}\,}$ en tout point de l'élément en fonction des déplacements aux nœuds ${\displaystyle \{u(x)\}=[N(x)]\{U\}\,}$ ;
2. On exprime les déformations en fonction des déplacements ${\displaystyle \{\epsilon (x)\}=[S]\{u(x)\},\,}$ d'où

${\displaystyle \{\epsilon (x)\}=[S]\{N(x)\}\{U\}=[B(x)]\{U\}\,}$

On écrit la loi de comportement du matériau qui relie les contraintes aux déformations :

${\displaystyle \{\sigma (x)\}=[D]\{\epsilon (x)\}\,}$

On écrit que le travail des forces externes appliquées à la structure pour un déplacement virtuel ${\displaystyle \delta U\,}$ est égal au travail interne des contraintes pour ce même déplacement :

${\displaystyle \{\delta U\}^{T}\{F\}=\int _{V}\{\delta \epsilon \}^{T}\{\sigma \}\,\mathrm {d} v\,}$

En explicitant, on a :

${\displaystyle \{\delta U\}^{T}\{F\}=\{\delta U\}^{T}\left(\int _{V}[B]^{T}[D][B]\,\mathrm {d} v\right)\{U\}\,}$

Comme cette relation est vraie pour tout déplacement virtuel, on en déduit :

${\displaystyle \{F\}=[K]\{U\}\,}$

avec ${\displaystyle [K]\,}$ sous sa forme plus générale : ${\displaystyle [K]=\int _{V}[B]^{T}[D][B]\,\mathrm {d} v\,}$

Remarques

• La relation ci-dessus montre que ${\displaystyle [K]\,}$ est symétrique.
• Le terme courant ${\displaystyle K_{ij}\,}$ de la matrice correspond à la force qui s'exerce sur le nœud ${\displaystyle i\,}$ lorsqu'on impose un déplacement unitaire du nœud ${\displaystyle j\,}$.

La symétrie de ${\displaystyle [K]\,}$ qui s'écrit ${\displaystyle K_{ij}=K_{ji}\,}$ correspond mécaniquement au théorème de réciprocité de Maxwell-Betti.

Les types d'éléments finis sont différenciés par le choix des fonctions N(x) et la nature de l'opérateur S (reliant déformations et déplacements), qui dépend du type de théorie élastique utilisée : théorie des poutres, théorie des plaques, théorie de la contrainte ou de la déformation plane, théorie des coques, théorie des corps de révolution, théorie de l'élasticité 3D.

Remarques Le processus de formulation d'un élément fini décrit ici est celui de la méthode directe (dite aussi méthode des déplacements). Il existe d'autres approches : la méthode des résidus pondérés, l'application du Principe des travaux virtuels ou des puissances virtuelles, et la minimisation de l'énergie potentielle. Ces approches sont équivalentes et aboutissent à la construction de la même matrice de rigidité.

Interpolation de la géométrie

${\displaystyle {\begin{matrix}x&=&\ &{\tilde {N}}_{1}(x)x_{1}+{\tilde {N}}_{2}(x)x_{2}\\\ &\ &+&{\tilde {N}}_{3}(x)x_{3}+{\tilde {N}}_{4}(x)x_{4}\end{matrix}}\,}$ Idem pour les autres coordonnées.

Coordonnées locales (cas 2D)

Élément isoparamétrique

Un élément est dit isoparamétrique si on prend les mêmes fonctions d'interpolation pour le déplacement et la géométrie.

${\displaystyle {\tilde {N}}(x)=N(x)\,}$

Autres classes d'éléments

Évaluation de K

La forme générale s'écrit :

${\displaystyle [K]=\int _{x}\int _{y}\int _{z}G(x,y,z)\,\mathrm {d} x\,\mathrm {d} y\,\mathrm {d} z\,}$

On passe en variables locales ${\displaystyle \xi ,\eta ,\zeta \,}$

On a ${\displaystyle \,\mathrm {d} x\,\mathrm {d} y\,\mathrm {d} z=\det(J)\,\mathrm {d} \xi \,\mathrm {d} \eta \,\mathrm {d} \zeta \,}$, où ${\displaystyle J\,}$ est la matrice jacobienne.

On est alors amené à calculer des intégrales du type :

${\displaystyle \int _{-1}^{+1}\int _{-1}^{+1}\int _{-1}^{+1}G(\xi ,\eta ,\zeta )\det J\ \,\mathrm {d} \xi \,\mathrm {d} \eta \,\mathrm {d} \zeta \,}$

Bénéfice de l'approche

On s'est ramené à un domaine d'intégration simple et invariant pour lequel on peut appliquer les formules de quadrature de Gauss :

${\displaystyle \int _{-1}^{+1}f(x)\,\mathrm {d} x=\sum _{k=1}^{n}H_{k}f(a_{k})\,}$

les ${\displaystyle a_{k}\,}$ et ${\displaystyle H_{k}\,}$étant tabulés.

Les ${\displaystyle a_{k}\,}$ sont appelés points d'intégration de l'élément ou encore points de Gauss de l'élément.