Méthode des éléments finis

Un article de Wikipédia, l'encyclopédie libre.

Vous avez de nouveaux messages (diff^?).

Simulation numérique d'un crash de véhicule : les cellules rectangulaires calculées par le mailleur sont visibles sur la peau géométrique du véhicule et représentent chacune un élément fini.

En analyse numérique, la méthode des éléments finis est utilisée pour résoudre numériquement des équations aux dérivées partielles représentant souvent analytiquement le comportement dynamique de certains systèmes physiques (en général mécaniques, thermodynamiques, acoustiques, ... etc).

Sommaire

1 Introduction
2 Méthode des Eléments Finis
3 Exemples issus de problèmes physiques
4 Applications et perspectives de la Méthode des Eléments Finis
- 4.1 Un outil pour l'ingénieur
- 4.2 La méthode des Eléments Finis Etendus (X-FEM)
  - 4.2.1 Origine
  - 4.2.2 Principes de base et applications
5 Voir aussi

[modifier] Introduction

La méthode des Eléments Finis fait partie des outils de mathématiques appliquées. Il s'agit de mettre en place, à l'aide des principes hérités de la formulation variationnelle ou formulation faible, un algorithme discret mathématique permettant de résoudre une équation aux dérivées partielles ou EDP sur un domaine compact avec conditions aux bornes et/ou dans l'intérieur du compact. On parle couramment de conditions de type Dirichlet (valeurs aux bornes) ou Neumann (normales aux bornes).

Il s'agit donc avant tout de la résolution d'un problème discrétisé, où, grâce à la formulation variationnelle, les solutions du problème vérifient des conditions d'existence plus faibles que celles des solutions du problème de départ. Comme de nombreuses autres méthodes numériques, outre l'algorithme de résolution en soi, se posent les questions de qualité de la discrétisation :

existence de solutions,
unicité de la solution,
stabilité,
convergence,
et bien sur mesure d'erreur entre une solution discrète et une solution unique du problème initial.

La partie 2 va présenter le cadre général de la méthode des éléments finis, ainsi que le cas pratique le plus courant considérant des équations aux dérivées partielles linéaires dont on cherche une approximation par des fonctions affines.

La présentation en partie 3 est essentiellement physique, notamment mécanique. Elle ne doit être considérée que comme une présentation des éléments constitutifs de la modélisation discrète utilisée en résistance des matériaux via la méthode des éléments finis. C'est une approche tout à fait valide, un bon exemple pédagogique. Elle apporte un biais certain quant à une approche plus générale, du fait notamment de la linéarité supposée des matériaux. Quelques éléments de compréhension mathématique additionnels sont présentés dans le second paragraphe (à rédiger).

[modifier] Méthode des Eléments Finis

[modifier] Principe général

La méthode des éléments finis permet donc de résoudre de manière discrète une EDP dont on cherche une solution approchée "suffisamment" fiable.

De manière générale, cette EDP porte sur une fonction $u$ , définie sur un domaine. Elle comporte des conditions aux bornes permettant d'assurer existence et unicité d'une solution.

La discrétisation passe par la définition d'un espace de fonctions tests approprié, sur lequel les solutions de la (re)formulation variationnelle de l'équation est exacte. Cela nécessite la définition d'un maillage du domaine en fragments triangulaires : les éléments finis.

Ce maillage permet aussi de définir une base fonctionnelle $b$ , sur laquelle on projette la fonction inconnue $u$ . On applique en outre la formulation variationnelle pour chacunes des fonctions de $b$ .

Dans ce cadre, et modulo une attention toute particulière au problème des bornes du domaine, on obtient une formulation algrébrique dite discrétisation du problème initial. La solution, si elle existe et est unique, de ce problème algébrique donne les composantes de la solution approchée dans une base $b$ .

La solution trouvée, il reste cependant à déterminer les caractéristiques de la méthode ainsi développée, notamment l'unicité de l'éventuelle solution ou encore la stabilité numérique du schéma de résolution. Il est essentiel de trouver une estimation juste de l'erreur liée à la discrétisation et montrer que la méthode ainsi écrite converge, c'est à dire que l'erreur tend vers 0 si la finesse du maillage tend elle aussi vers 0.

Dans le cas d'une EDP linéaire avec opérateur symétrique (comme l'est $Δ$ ), il s'agit au final de résoudre une équation algébrique linéaire, inversible dans le meilleur des cas.

[modifier] Dimensions

On développe ici la méthode des éléments finis en deux dimensions, c'est à dire des applications d'un support de dimension algébrique 2, à valeurs réelles. On suppose que les équations étudiées sont de degré deux de différentiation.

La méthode est généralisable à des cadres d'espaces de dimension différente ou de différentiation supérieure :

pour l'image, la dimension ne compte pas vraiment, cependant cela nécessite des écritures plus complètes. Les cas les plus couramment rencontrés dans une vie d'étudiant sont la dimension 1 (comme ici), 2 ou 3 (pour des problèmes de mécanique), 6 ou 12 (pour des problèmes d'électromagnétisme "réels ou complexes" respectivement)
pour le support, on parle de méthode des segments finis pour un support en dimension 1, et de volumes finis en dimension 3.
le degré de différentiation supérieurs sont ramenés à un degré moindre par la méthode classique de réduction de degré : on fait intervenir des variables supplémentaires, i.e. des dérivées partielles des variables de départ (ex classique : les EDP de la mécanique statique des poutres font intervenir la dérivation partielle d'ordre 4). Il est parfois possible, pour de degrés supérieurs, d'appliquer plusieurs fois les méthodes de formulation variationelles afin d'obtenir des ordres plus faibles - en tout cas lorsque le degré de dérivation est paire.

Bien que théoriquement la méthode soit transposable en dimensions supérieures du support, techniquement la complexité de création des discrétisations croît avec la dimension... et pratiquement, on résoud rarement des problèmes en dimension supérieures à 3 - y compris des problèmes de dynamique en espace à 3 dimension qui sont traités en réalité avec une méthode mixte éléments finis "en espace"/différence finies "en temps".

[modifier] Cadre algébrique, analytique et topologique

Soit le plan $\mathbb{R}^2$ . Soit l'ensemble image $\mathbb{R}$ .

Soit un domaine (ouvert borné et connexe) $Ω$ de $\mathbb{R}^2$ , de bord $δΩ$ , et d'adhérence (compacte) $\bar{\Omega}$ . Pour simplifier les représentations on suppose le bord polygonal.

Soit les fonctions de $\bar{\Omega}$ dans $\mathbb{R}$ différentiables sur $\bar{\Omega}$ (compact) et deux fois différentiables sur $Ω$ (ouvert). Notons que de telles fonctions sont continues et différentiables sur le bord du compact. Soit $V (Ω)$ l'ensemble de ces fonctions. Notons que $V$ est un espace vectoriel de dimension infinie et $V 0$ le sous espace vectioriel de fonctions de $V$ nulles sur le bord $δΩ$ .

Soit les applications continues sur $\bar{\Omega}$ et différentiables sur $Ω$ , de carré sommables sur $\bar{\Omega}$ et de gradient de carré sommable sur $\bar{\Omega}$ (ou de dérivées partielles de carré sommable, ce qui revient au même avec le support de dimension finie). Nommons cet espace $\mathcal{H}_1(\Omega)$ . On dote cet espace vectoriel d'un produit scalaire issu de celui de $L 2$ tel que si $(u, v)$ appartiennent à cet espace alors le produit scalaire de $u$ et $v$ est :

$< u | v >_{ \mathcal{H}_1(\Omega) } = \int_\Omega{ ( \nabla u . \nabla v + u v ) d\omega }$

On note $\mathcal{H}_1^0(\Omega)$ le sous espace vectoriel de $\mathcal{H}_1(\Omega)$ dont les fonctions sont nulles sur le bord $δΩ$ .

Soit un maillage $\mathcal{T}_{n,N} = \left\{X_i\right\}_,{i \in [1,...,N]}, N>n$ du compact $\bar{\Omega}$ , tel que

$\left(X_i\right)_{i \in [1,...,n]} \in \Omega$ , i.e. des points ou nœuds intérieurs
$\left(X_i\right)_{i \in [n+1,...,N]} \in \delta \Omega$ , i.e. des points ou nœuds du bord

On défini une relation $\mathcal{R}$ en trois infice (i,j,k) qui est vraie si les nœuds $X i$ , $X j$ et $X k$ forment un triangle $T j, k, l$ non dégénéré qui ne contient pas d'autre point $X l$ .

Soit les fonctions test continues sur $Ω$ $w i$ telles que

$w i$ est continue sur $\bar{\Omega}$ ,
$w i$ est affine sur chaque triangle $T j, k, l$ où $\mathcal{R}_{j,k,l}$ est vraie
$w i (X j) = δ i j$ , i.e. 1 si les deux indices sont égaux, 0 sinon

La famille $\left\{w_i\right\}_{i \in [1...N]}$ est libre. Notons $V N$ l'espace engendré par cette famille, qu'on note désormais $b$ .

[modifier] Cas organique

[modifier] Hypothèses

Soit l'équation aux dérivées partielles suivante : Soit $u$ telle que

u est continûment différentiable sur $\bar{\Omega}$
u est deux fois différentiable sur $Ω$
f est une fonction continue sur $\bar{\Omega}$ (de carré sommable suffit)
$\Delta u(X) - k^2 u(X) = \frac{\partial^2 u}{\partial x^2}(X) + \frac{\partial^2 u}{\partial y^2}(X) - k^2 u(X)= -f(X), \forall X \in \Omega$ (1)
$u(\bar{\Omega}) = 0$ (condition au bord de Dirichlet)

On démontre qu'il existe une solution unique à ce problème d'EDP.

[modifier] Formulation faible

Soit l'ensemble $V 0$ des fonctions de $V$ nulle sur le bord de $Ω$ . Soit $v \in V_0$ , alors on écrit :

$v \Delta u(X) - k^2 v u(X) = -v f(X), \forall X \in \Omega$

On montre que tous ces termes sont sommables sur $Ω$ . On somme donc pour obtenir l'équation (2):

∫	vΔudω − k²	∫	vudω = −	∫	− vfdω
Ω		Ω		Ω

Soit une formule bien connue ( le "." de la formule est le produit scalaire) :

$\int_\Omega{ v \Delta u d\omega } = - \int_\Omega{ ( \nabla u . \nabla u ) d\omega } + \int_{\delta \Omega }{ \frac{ \partial u }{ \partial n } v ds }$

Dans cette formulation, v est nulle sur le bord donc il ne reste plus que le terme $- \int_\Omega{ ( \nabla u . \nabla u ) d\omega }$ d'où le changement dans l'équation (2)

$\int_\Omega{ \nabla u . \nabla v d\omega } + k^2 \int_\Omega{ v u d\omega } = \int_\Omega{ -v f d\omega }$

Ceci est la formulation faible. On démontre que modulo l'hypothèse que $u$ est deux fois différentiable, il y a équivalence entre cette formulation et celle du problème (1) .

[modifier] Remarque algébrique

soit $a$ l'opérateur bilinéaire symétrique (de $V 2$ dans $\mathbb{R}$ ) et $\mathcal{L}$ l'opérateur linéaire (de $V$ dans $\mathbb{R}$ ) tels que :

$a(u,v) = \int_\Omega{ \left(\nabla u . \nabla v + k^2 u v \right) d\omega }$

$\mathcal{L}(v) = \int_\Omega{ f v d\omega }$

Alors

le problème se reformule ainsi : $\forall v \in V_0, a( u, v ) = \mathcal{L}(v)$

$a$ est un produit scalaire : symétrique, défini, positif, continu selon la norme $\mathcal{H}_1^0$ (cf. espace de Sobolev).

$\mathcal{L}$ est continu selon la norme $\mathcal{H}_1^0$ c'est donc une forme linéaire

Ces deux éléments sont les clefs de la démonstration de l'équivalence des problèmes et de l'existence et de l'unicité de la solution du problème, où l'on montre aussi que $u$ est la solution unique du problème d'optimisation de la fonctionnelle suivante :

$\forall v \in V_0, \mathcal{J}(v) = \frac{1}{2} a( v, v ) - \mathcal{L}(v)$

[modifier] Discrétisation

Soit le maillage $\mathcal{T}_{n,N}$ et la base $b$ associée. Parce que la condition de Dirichlet impose des fonctions nulles aux bords, on utilise uniquement la sous-base $b'$ limitée aux points intérieurs de $Ω$ et le sous-espace vectoriel $\mathcal{T}_n$ engendré par la famille $b'$ .

On cherche la solution $\bar{u}$ du problème discrétisé ainsi :

$\exists! ? \bar{u} \in V_n^0 | \forall v \in V_n^0, a( \bar{u}, v ) = \mathcal{L}(v)$

Or dans cet espace discrétisé, dire que tout vecteur vérifie la proposition précédente est équivalent à dire que tous les vecteurs de la base vérifient la proposition. Si l'on décompose la solution $\bar{u}$ dans la base des $w i$ intérieurs, en composantes $u i$ , on obtient :

$\exists! ? { u_i } \in \mathbb{R}^n | \forall l \in [1,...,n] \sum_{i =1}^n { u_i a( w_i, w_j ) } = \mathcal{L}(w_j)$

Soit si l'on note

la matrice $A$ ayant pour composantes les $a (w i, w j)$ ,
le vecteur $U$ ayant pour composantes les $u i$
le vecteur $B$ pour composantes les $\mathcal{L}( w_j )$

alors ce problème revient à résoudre l'équation linéaire de $n$ équations à $n$ inconnues :

$A U = B$

Notons que $A$ par construction est symétrique, et que parce que $a$ est un produit scalaire et que $b' = \{ w_i \}_{ i \in [1,...,n] }$ est une base, alors $A$ est elle-même symétrique, définie, positive donc inversible, et on obtient donc l'éexistence et l'unicité de $U = A - 1 B$ !

[modifier] Astuces de calcul

Ceci fait parti de la méthode, et plus précisemment de l'algorithme d'application.

En utilisant la relation de Chasles, on scinde la somme définie plus haut sur chacun des triangles formant le domaine de calcul :

$a( u, v ) = \sum_{i,j,k | i<j<k,\mathcal{R}_{i,j,k}}{ \int_{T_{i,j,k}}{ ( \nabla u . \nabla v + k^2 u v ) d\omega} }$

$\mathcal{ L }( v ) = \sum_{i,j,k | i<j<k,\mathcal{R}_{i,j,k}}{\int_{T_{i,j,k}}{ f v d\omega }}$

Soit une fonction $w l$ , on sait que cette fonction par construction est nulle sur tout triangle $T i, j, k$ si l n'appartient pas au trio ${i, j, k}$ .(un dessin convainct !)

Ainsi, on sait écrire les éléments suivants ( $l \neq m$ ) qui semblent complexes mais simplifient les calculs des éléments de la matrice de raideur $A$ et du vecteur des conditions aux bord $B$ :

$a( w_l, w_m ) = \sum_{ i=1 | i \neq l, i \neq m, \mathcal{R}_{i,l,m}}^n { \int_{T_{i,l,m}}{ ( \nabla w_l . \nabla w_m + k^2 w_l w_m ) d\omega} }$

$a( w_l, w_l ) = \sum_{i,j | i<j, i \neq l, j \neq l, \mathcal{R}_{i,j,l}}{\int_{T_{i,j,l}}{ \left( [[\nabla w_l||^2 + k^2 w_l^2 \right) d\omega }}$

$\mathcal{ L }( w_l ) = \sum_{i,j | i<j, i \neq l, j \neq l, \mathcal{R}_{i,j,l}}{\int_{T_{i,j,l}}{ f w_l d\omega }}$

[modifier] Alternative

La condition qui suit est très différente de celle de Dirichlet et donne des résultats qui n'ont rien à voir.

On pose comme condition au bord que la dérivée normale existe sur le bord $\bar{\Omega}$ , et que la condition de Neumann $\frac{\partial u}{\partial n}(\bar{\Omega}) = 0$ est vérifiée.

Notons que si la fonction est supposée différentiable au bord $D_u(X).n = 0, \forall X \in \bar{\Omega}$ voire si elle admet un gradient $<\nabla u(X)| n> = 0, \forall X \in \bar{\Omega}$ .

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

Par la suite, la différence, comme on l'a signalée, réside surtout dans le choix des vecteurs de base our la discrétisation : il faut conserver les fonctions tests propres aux nœuds du bord.

[modifier] Exemples issus de problèmes physiques

[modifier] Historique

Analyse des structures née vers 1850
- RDM (Résistance des matériaux) ⇒ calculs « manuels »
  Maxwell, Castigliano, Mohr
- Concept d'éléments finis né vers 1940
  Newmark, Hrenikoff, Mc Henry, Courant
- Développement réel depuis 1960
  Calcul numérique sur ordinateur

[modifier] Domaines d'application

Calcul de structures, étude des contacts, électricité, électromagnétisme, hydraulique, aérodynamique, ...

[modifier] Principe

Le milieu continu est « idéalisé » par la subdivision en un nombre fini d'éléments dont le comportement est représenté par un nombre finis de paramètres.
La résolution du problème global, obtenu par assemblage des éléments, suit les règles qui régissent les structures discrètes.

[modifier] Les difficultés

D'ordre théorique : formulation des éléments
D'ordre pratique :
- Discrétisation du milieu continu (maillage)
- Qualité des résultats (convergence de la méthode)

[modifier] Exemple de problème discret : un réseau électrique

Équation locale du composant e :

$\left. \begin{matrix} I_i^e = {1 \over R^e} (U_i^e - U_j^e) \\ I_j^e = {1 \over R^e} (U_j^e - U_i^e) \end{matrix} \right\} \Longrightarrow \begin{Bmatrix} I_i^e \\ I_j^e \end{Bmatrix} = {1 \over R^e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{Bmatrix} U_i^e \\ U_j^e \end{Bmatrix} \,$

soit $\left\{ I^e \right\} = \left[ K^e \right] \left\{ U^e \right\} \,$
On écrit :
- La continuité des potentiels en chaque connexion $i \,$
- L'équilibre des courants à chaque connexion
- L'adjonction des courants externes $P^i \,$
On obtient l'équation globale du système assemblé :

$P^i = \sum_{j=1}^m \sum_{e=1}^m K_{ij}^e U_j^e \,$

soit $\left\{ P \right\} = \left[ K \right] \left\{ U \right\} \,$

[modifier] Définition d'un élément fini

* Problème modélisé :		* Problème simplifié :

En calcul de structures, un élément fini est caractérisé par deux matrices :

La matrice de raideur $\left[ K \right] \,$
La matrice de masse $\left[ M \right] \,$

[modifier] Formulation d'un élément fini

[modifier] Définitions et notations

Déplacement : $\left\{ U \right\} \,$

C'est un vecteur dont chaque composante est également appelée degré de liberté

3 ddl de translation : $U_x,U_y,U_z \,$
3 ddl de rotation : $\theta_x,\theta_y,\theta_z \,$

Déformation : (voir Tenseur des déformations) $\left[ \epsilon \right] \,$

C'est le rapport de l'allongement à la longueur initiale.

En petites déformations, on a $(A'B')^2 = \left(\ dx + {\partial u \over \partial x} \ dx\right)^2 + \left({\partial v \over \partial x} \ dx\right)^2 \,$

$\epsilon_x = {{A'B'-AB} \over AB} \,$

Comme $AB=\ dx \,$ , on a $A'B'=(1+\epsilon_x)\ dx \,$

$\Rightarrow 2\epsilon_x+\epsilon_x^2 = 2{\partial u \over \partial x}+{\partial u \over \partial x}^2+{\partial v \over \partial x}^2 \,$

On néglige les termes d'ordre 2 :

$\epsilon_x = {\partial u \over \partial x} \,$

Remarque : $\epsilon \,$ est sans dimension

Contrainte : (voir Tenseur des contraintes) $\left[ \sigma \right] \,$

Elle représente les efforts internes qui s'appliquent dans la structure.

$\sigma_{11},\sigma_{22},\sigma_{33} \,$ : contraintes normales

$\sigma_{12},\sigma_{13},\sigma_{23} \,$ : contraintes de cisaillement

Une contrainte est homogène à une pression (N/m²)
Il existe un système de coordonnées dans lequel $\left[ \sigma \right] \,$ est une matrice diagonale : $\left[ \sigma \right] = \begin{bmatrix}\sigma_X & 0 & 0 \\ 0 & \sigma_Y & 0 \\ 0 & 0 & \sigma_Z\end{bmatrix} \,$

$\sigma_X,\sigma_Y,\sigma_Z \,$ sont appelées les contraintes principales

Relations Contraintes-Déformations : (voir Loi de Hooke) ≡ lois de comportement

$\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}\,$

$E \,$ : module de Young (N/m²)

$\nu \,$ : coefficient de Poisson (sans dimension)

On utilise parfois le module de cisaillement : $G = {E \over {2(1+\nu)}} \,$
Pour un matériau isotrope, il n'y a que 2 paramètres indépendants. Il y en a 6 pour un matériau orthotrope et 21 pour un matériau anisotrope
En notation matricielle, on écrit : $\left\{ \sigma \right\} = \left[D\right]\left\{ \epsilon \right\} \,$

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

Énergie de déformation : $W \,$

$W = {1 \over 2}\int_{v}\sigma.\epsilon.\ dv \,$

Travail d'une force :

C'est le produit de la force par le déplacement de son point d'application : $\tau_F = {1 \over 2}F.U_F \,$

Moment :

C'est une force appliquée sur un ddl de type rotation

[modifier] Équations fondamentales

Équations d'équilibre local :

${\partial \sigma_x \over \partial x} + {\partial \sigma_{xy} \over \partial y} + {\partial \sigma_{xz} \over \partial z} + F_x = 0\,$

Image:Elasticite lineaire

${\partial \sigma_y \over \partial y} + {\partial \sigma_{xy} \over \partial x} + {\partial \sigma_{yz} \over \partial z} + F_y = 0\,$

${\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 :

${\epsilon_x= {\partial u \over \partial x}} \qquad {\epsilon_y= {\partial v \over \partial y}} \qquad {\epsilon_z= {\partial w \over \partial z}} \,$

$\gamma_{xz}={\partial w \over \partial x} + {\partial u \over \partial z} \,$

$\gamma_{xy}={\partial u \over \partial y} + {\partial v \over \partial x} \,$

$\gamma_{yz}={\partial w \over \partial y} + {\partial v \over \partial z} \,$

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

Signification du coefficient de Poisson :

Si on applique au barreau une contrainte $\sigma_x \,$ , on observe un rétrécissement dans la direction y correspondant à une déformation $-\nu \sigma_x \over E \,$
Quelques valeurs usuelles :

Acier : E = 2,1E11 N/m² = 210 000 MPa $\nu \,$ = 0,3
Aluminium : E = 7E10 N/m² = 70 000 MPa $\nu \,$ = 0,3

Remarques : On a toujours $\nu \le \,$ 0,5
Quand $\nu \to \,$ 0.5, Le matériau est dit incompressible.

[modifier] 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é : $u(x) = a_1 + a_2 x \,$

On a $u(0) = u_1 \,$ et $u(L) = u_2 \,$

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

qu'on écrit symboliquement : $u(x) = [N(x)]\{U\} \,$ avec $\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 :

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

$\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 :

$\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 :

$\{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 :

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

Soit une relation du type :

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

En explicitant :

$[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 3 matrices :

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

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

$[A] \,$ : Transformation des contraintes aux forces

[modifier] Formulation générale (méthode directe)

La démarche est la suivante :

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

d'où $\{\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 : $\{\sigma(x)\}=[D]\{\epsilon(x)\}\,$
On écrit que le travail des forces externes appliquées à la structure pour un déplacement virtuel $\delta U \,$ est égal au travail interne des contraintes pour ce même déplacement :

$\{\delta U\}^T\{F\}=\int_V \{\delta \epsilon\}^T\{\sigma\} dv \,$

En explicitant, on a :

$\{\delta U\}^T\{F\}=\{\delta U\}^T \left (\int_V [B]^T[D][B] dv \right ) \{U\} \,$

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

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

avec $[K] \,$ sous sa forme plus générale : $[K]=\int_V [B]^T[D][B] dv \,$

Remarques :

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

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

Qu'est-ce qui va différencier les différents types d'éléments finis ?

Le choix des fonctions N(x)
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'elasticité 3D

Remarques :
Nous avons décrit le processus de formulation d'un élément fini dans le cadre 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
La minimisation de l'énergie potentielle

Toutes ces approches sont équivalentes et aboutissent à la construction de la même matrice de rigidité.

[modifier] Éléments finis en contraintes

Au lieu de rechercher une solution approchée en déplacement, on peut aussi rechercher la solution approchée en contrainte.

Dans le cas de la mécanique, l'application du principe des puissances virtuelles donne de manière non triviale les théorèmes énergétiques. On peut aboutir au même résultat en quelques lignes en écrivant l'erreur en relation de comportement.

L'approche en contrainte consiste à rechercher dans l'espace des champs de contraintes admissibles celui qui réalise le minimum de l'énergie complémentaire.

Cette approche est plus précise que l'approche en déplacement mais elle est peu développée du fait de la difficulté que l'on a à générer des champs de contraintes de divergence donnée.

[modifier] Étude des fonctions N(x)

Dans le cas général de l'élasticité tridimensionnelle, ce sont en fait des fonctions de x, y, z.
Les fonctions les plus couramment utilisées sont des polynômes.

Polynôme de degré 1 : élément linéaire (2 nœuds par arête)
Polynôme de degré 2 : élément parabolique (3 nœuds par arête)
Polynôme de degré 3 : élément cubique (4 nœuds par arête)

etc.

Les fonctions N(x) sont appelées fonctions de forme ou fonctions d'interpolation de l'élément.

[modifier] Les éléments isoparamétriques

$[K]=\int_V [B]^T[D][B] dv \,$

[modifier] Problème

Soit on construit $[K] \,$ pour un certain nombre d'éléments de forme et de géométrie figée ⇒ nécessité, pour mailler une structure complexe, d'utiliser un grand nombre d'éléments.
Soit on utilise des éléments à géométries variables ⇒ il faut reconstruire $[K] \,$ à chaque fois.

[modifier] D'où l'idée

Pour construire la matrice de raideur d'un élément à géométrie variable, on va utiliser des fonctions d'interpolation pour décrire non seulement le champ de déplacement de l'élément mais également sa géométrie. De plus, on va travailler en coordonnées locales.

[modifier] Interpolation de la géométrie

$\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.

[modifier] Coordonnées locales (cas 2D)

[modifier] É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.
$\tilde N(x) = N(x) \,$

[modifier] Autres classes d'éléments

[modifier] Évaluation de K

La forme générale s'écrit :
$[K]=\int_x \int_y \int_z G(x,y,z) dxdydz \,$
On passe en variables locales $\xi,\eta,\zeta \,$
On a $dxdydz = \det J\ d\xi d\eta d\zeta \,$

$J \,$ s'appelle la matrice Jacobienne.

On est alors amené à calculer des intégrales du type :
$\int_{-1}^{+1} \int_{-1}^{+1} \int_{-1}^{+1} G(\xi ,\eta ,\zeta) \det J\ d\xi d\eta d\zeta \,$

[modifier] 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 :
$\int_{-1}^{+1} f(x) dx = \sum_{k=1}^n H_k f(a_k) \,$ les $a_k \,$ et $H_k \,$ étant tabulés.
Les $a_k \,$ sont appelés points d'intégration de l'élément ou encore points de gauss de l'élément.

[modifier] Cas particulier : les éléments axisymétriques

Décomposition en série de Fourier :

$u(r,\theta,z) = \sum_n u_n^s(r,z)cos (n\theta) + \sum_n u_n^s(r,z)sin (n\theta) \,$

L'axisymétrie correspond à la restriction $n=0 \,$ de cette décomposition.

$u(r,\theta,z) = u_n^s (r,z) \,$

Remarque importante : Pour utiliser ce type d'élément, le problème doit être globalement axisymétrique :

la géométrie
les conditions limites
le chargement

[modifier] Processus de calcul (cas statique)

Maillage
Construction de la matrice de raideur de chaque élément $[K^e] \,$
Assemblage de la matrice globale $[K] \,$
Construction du vecteur chargement $\{F\} \,$
Élimination de certains ddl (si besoin)
Résolution : $\{U\} = [K^{-1}]\{F\} \,$
Calcul des quantités dérivées de $\{U\} \,$

$\{\epsilon^e\} = [B^e] \{U^e\} \,$

$\{\sigma^e\} = [D^e] \{\epsilon^e\} \,$

$W^e = {1 \over 2} \{\epsilon^e\}^T \{\sigma^e\} \,$

etc.

[modifier] Références

O. C. Zienkiewicz, R. L. Taylor, J.Z. Zhu : The Finite Element Method: Its Basis and Fundamentals

 Butterworth-Heinemann; 6 edition (March 21, 2005)ISBN: 0750663200

K. J. Bathe : Numerical methods in finite element analysis

 Prentice-Hall (1976) ISBN: 0136271901

Chu-Kia Wang : Matrix methods of structural analysis

 International Textbook Co; 2d ed edition (1970) ISBN: 0700222677

R. H Gallagher : Introduction aux éléments finis

 Pluralis (1977)

N. Willems : Matrix analysis for structural engineers

 Prentice-Hall (1968) ASIN: B0006BUP26

[modifier] Approche mathématique/théorique

Durant cet exemple, si f est une fonction, la dérivée partielle de f par rapport à x sera notée f_x.

La meilleure façon d'introduire ce sujet est de donner un exemple simple (un « problème modèle »). Nous allons prendre l'équation de Laplace sur le pavé [0,1[x[0,1[. Tout d'abord, notons

$\mathrm{\Omega} := \left[ 0,1\right[ \times\left[ 0,1\right[=\{\left( x,y\right); 0 \leq x,y < 1\}.$

Soit g une fonction de $Ω$ vers F, l'ensemble des scalaires (l'ensemble des réels R ou le plan complexe C). Notre problème est de trouver une fonction u du plan R² vers F qui vérifie

u est deux fois différentiable sur R²
u_xx+u_yy=g dans $Ω$
u(x,y)=u(x+1,y)=u(x,y+1), c'est-à-dire, u est périodique en x et en y, de période 1.

Le problème peut être réécrit de la manière suivante :

u appartient à $H 2 (Ω)$ ( $H 2 (Ω)$ est un espace de Hilbert, $H^2(\Omega)=\left\{ v \in L^2(\Omega) , v,v' \in L^2(\Omega)\right\}$ )
Lu=g avec g un élément de $V=H_0^1(\Omega)=\left\{ v \in L^2(\Omega) , v \in L^2(\Omega) ,v(0)=v(1)=0 \right\}$

L étant l'opérateur différentiel défini ainsi

L f = f x x + f y y .

L est appelé opérateur de Laplace. Nous cherchons maintenant un u dans V tel que Lu=g.

[modifier] Formulation faible

Soit $φ$ une forme linéaire sur V.

$\phi\left(tu+v\right)=t\phi\left(u\right)+\phi\left(v\right)$

Notons par V^* l'ensemble de ces fonctions. Ainsi, les propriétés suivantes sont équivalents

L u = g

$\forall \phi \in \mathrm{V}^*: \phi\left(\mathrm{L}u\right)=\phi\left(g\right).$

[modifier] Fonctions tests

L'ensemble V^* est plus grand que ce dont on a besoin. Par conséquent, nous choisissons la fonction de la sorte

$\phi_v\left(f\right)=\int_0^1\int_0^1 v\left(x,y\right) f\left(x,y\right) dxdy,$

où $v \in V$

Maintenant, nous écrirons

∫

pour l'intégrale double $\int_0^1\int_0^1$ . En intégrant par parties,

$\int_\Omega v u_{xx}\;dxdy=\int_0^1\left[v u_x\right]_0^1\;dy - \int_\Omega \left(v_x u_x\right)\;dxdy$

$\int_\Omega v u_{yy}\;dxdy=\int_0^1\left[v u_y\right]_0^1\;dx - \int_\Omega \left(v_y u_y\right)\;dxdy$

et du fait des conditions aux limites, les termes $\int_0^1 \left[ v u_x \right]_0^1\;dy$ et $\int_0^1 \left[ v u_y \right]_0^1\;dx$ s'annulent.

$\int_\Omega v\left(u_{xx}+u_{yy}\right)\;dxdy = -\int_\Omega \left(v_x u_x + v_y u_y \right)\;dxdy.$

L'écriture précédente est équivalente à :

$\psi\left(u,v\right) := \phi_v\left(\mathrm{L}u\right) = -\int_\Omega \left(u_x v_x + u_y v_y\right) = \int_\Omega g\cdot v = \phi_v \left(g\right).$

La fonction ψ de u est de v est en fait bilinéaire sur $H^2(\Omega) \times V$ et c'est la forme bilinéaire associée à L. Les fonctions v sont appelées fonctions tests.

Le problème variationnel abstrait s'écrit :

Trouver $u \in V$ tel que, quelque soit $v \in V$ , on ait

$\psi \left( u,v \right) = \phi_v \left( g \right)$

[modifier] Un ensemble minimal de fonctions tests

En analyse fonctionnelle, nous savons qu'il est inutile d'essayer toutes les fonctions tests v dans V. En fait, si E={e₁,e₂,...} est un sous ensemble de V, l'ensemble des combinaisons linéaires est dense dans V. Par conséquent, on peut utiliser ces fonctions comme fonctions tests. Maintenant, nous allons rechercher une solution au problème

$\textrm{(*)}\ \ \ \forall e_j \in \mathrm{E} : \psi\left(u,e_j\right)=\int_\Omega g\cdot e_j.$

[modifier] Première discrétisation

De manière à rendre ce processus implémentable, nous devons choisir un sous ensemble fini de E. Prenons F=F_n={e₁,e₂,...,e_n} un ensemble fini, F est un sous ensemble de E, et nous voulons résoudre le problème

$\textrm{(**)}\ \ \forall e_j \in \mathrm{F}_n : \psi\left(u_n,e_j\right)=\int_\Omega g\cdot e_j$

avec l'idée que quand n tend vers l'infini (et F_n croit vers E), les solutions u_n devront converger vers la solution u de (*).

Ce problème est indéterminé, par conséquent, une autre discrétisation est nécessaire.

[modifier] Seconde discrétisation

Le problème variationnel dans F a une solution et est unique (conséquence du théorème de Lax-Milgram sur V, $ψ$ est une forme bilinéaire sur $V \times V$ , V est elliptique et $φ$ est une forme bilinéaire continue.

$u_n = \sum_{k=1}^n a_k e_k$

$g_n = \sum_{k=1}^n b_k e_k.$

Remplaçons dans (**) et en développant, on obtient :

$\textrm{(***) } \sum_{k=1}^n a_k \psi\left(e_k, e_j\right) = \sum_{k=1}^n b_k \int_\Omega e_k e_j,\; j=1,\dots,n.$

Se pose maintenant la question : comment obtenir un g_n acceptable ? Il existe plusieurs approches qui dépendent du choix de l'ensemble E. Si E est une base de Fourier, g_n est la projection de g sur l'ensemble des combinaisons linéaires de F, mais d'autres approches existent.

Une autre approche peut être de calculer directement $\int_\Omega g\cdot e_j$ à partir d'une méthode de quadrature numérique.

[modifier] Problème sous forme matricielle

Le lecteur astucieux aura noté que le problème précédent peut être formulé sous forme matricielle :

P a = Q b

où a, b sont des vecteurs colonnes

$a = \begin{bmatrix} a_1\\ a_2\\ \vdots\\ a_n \end{bmatrix}, b = \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix}.$

Les matrices P et Q sont données par (***):

$\mathrm{P}_{jk} = \psi\left(e_k, e_j\right),$

Q_jk =	∫	e_ke_j.
	Ω

P est appelé matrice de rigidité et Q, matrice de masse.

[modifier] Algorithme

La méthode des éléments finis doit être conduite ainsi

On calcule la matrice de rigidité P
On détermine le membre de droite, soit par l'intermédiaire de la matrice de masse, soit directement par une méthode de quadrature
On résout le problème Pa=Qb. Selon la base qui a été choisie et selon les données du problème, il faut choisir la méthode d'inversion la plus efficace.
On peut écrire u grâce à a.

Note : Dans notre problème,

si l'on choisit une base lagrangienne pour F, la matrice P sera tridiagonale.
étant donné que notre forme bilinéaire $ψ$ est symétrique, P est symétrique
$ψ$ étant V-elliptique, P est symétrique définie positive

on peut utiliser une méthode spécifique pour les systèmes linéaires à matrice tridiagonale ou une factorisation de Cholesky par bande.

[modifier] Applications et perspectives de la Méthode des Eléments Finis

[modifier] Un outil pour l'ingénieur

[modifier] La méthode des Eléments Finis Etendus (X-FEM)

La Méthode des Eléments Finis (MEF) est donc un outil bien maîtrisé actuellement, tant d'un point de vue recherche et développement que d'un point de vue utilisation dans l'industrie. C'est une méthode robuste qui a fait ses preuves, mais les défis d'aujourd'hui et de demain présentent de nouveaux enjeux et trouvent des limites à la MEF. L'accent est mis, dans la recherche, sur l'introduction de plus de physique dans les solutions numériques que l'on souhaite calculer. Récemment, le concept d'Eléments Finis Etendus (X-FEM) a été introduit.

[modifier] Origine

[modifier] Principes de base et applications

[modifier] Voir aussi

Récupérée de « http://fr.wikipedia.org../../../m/%C3%A9/t/M%C3%A9thode_des_%C3%A9l%C3%A9ments_finis.html »

Catégories : Mécanique des milieux continus • Équations différentielles numériques

Méthode des éléments finis

Un article de Wikipédia, l'encyclopédie libre.

Sommaire

[modifier] Introduction

[modifier] Méthode des Eléments Finis

[modifier] Principe général

[modifier] Dimensions

[modifier] Cadre algébrique, analytique et topologique

[modifier] Cas organique

[modifier] Hypothèses

[modifier] Formulation faible

[modifier] Remarque algébrique

[modifier] Discrétisation

[modifier] Astuces de calcul

[modifier] Alternative

[modifier] Exemples issus de problèmes physiques

[modifier] Historique

[modifier] Domaines d'application

[modifier] Principe

[modifier] Les difficultés

[modifier] Exemple de problème discret : un réseau électrique

[modifier] Définition d'un élément fini

[modifier] Formulation d'un élément fini

[modifier] Définitions et notations

[modifier] Équations fondamentales

[modifier] Exemple de formulation : Barre en traction

[modifier] Formulation générale (méthode directe)

[modifier] Éléments finis en contraintes

[modifier] Étude des fonctions N(x)

[modifier] Les éléments isoparamétriques

[modifier] Problème

[modifier] D'où l'idée

[modifier] Interpolation de la géométrie

[modifier] Coordonnées locales (cas 2D)

[modifier] Élément isoparamétrique

[modifier] Autres classes d'éléments

[modifier] Évaluation de K

[modifier] Bénéfice de l'approche

[modifier] Cas particulier : les éléments axisymétriques

[modifier] Processus de calcul (cas statique)

[modifier] Références

[modifier] Approche mathématique/théorique

[modifier] Formulation faible

[modifier] Fonctions tests

[modifier] Un ensemble minimal de fonctions tests

[modifier] Première discrétisation

[modifier] Seconde discrétisation

[modifier] Problème sous forme matricielle

[modifier] Algorithme

[modifier] Applications et perspectives de la Méthode des Eléments Finis

[modifier] Un outil pour l'ingénieur

[modifier] La méthode des Eléments Finis Etendus (X-FEM)

[modifier] Origine

[modifier] Principes de base et applications

[modifier] Voir aussi

Views

Navigation

Contribuer

Recherche

Autres langues