Méthode des éléments finis 257600 31170911 2008-06-30T19:13:56Z Badmood 4507 url externe --> url interne {{voir homonymes|Élément}} [[Image:FEM example of 2D solution.png|thumb|right|solution bidimensionnelle d'une équation [[magnétostatique]] obtenue par éléments finis (les lignes donnent la direction du champ et la couleur son intensité)]] [[Image:Example of 2D mesh.png|thumb|right|Maillage utilisé pour l'image du haut (le maillage est plus resserré autour de la zone d'intérêt)]] [[Image:FEM car crash1.jpg|thumb|right|Simulation numérique d'un [[Essai de choc]] sur une voiture: les cellules utilisées pour le [[maillage]] sont visibles sur la surface du véhicule.]] En [[analyse numérique]], la '''méthode des éléments finis''' est utilisée pour résoudre numériquement des [[équation aux dérivées partielles|équations aux dérivées partielles]]. Celles-ci peuvent par exemple représenter [[Analyse (mathématiques)|analytiquement]] le comportement [[dynamique]] de certains systèmes physiques ([[mécanique]]s, [[thermodynamique]]s, [[acoustique]]s, etc.).En anglais :"Finite element analysis" Concrètement, cela permet par exemple de calculer ''numériquement'' le comportement d'objets même très complexes, à condition qu'ils soient continus et décrits par une équation aux dérivées partielles linéaire : mouvement d'une corde secouée par l'un de ses bouts, comportement d'un [[fluide]] arrivant à grande vitesse sur un obstacle, déformation d'une structure métallique, etc. == Introduction == La '''méthode des élé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 [[Algorithmique|algorithme]] discret mathématique permettant de rechercher une solution approchée d' une [[équation aux dérivées partielles]] (ou EDP) sur un domaine [[Compacité (mathématiques)|compact]] avec conditions aux bords et/ou dans l'intérieur du compact. On parle couramment de conditions de type [[Dirichlet]] (valeurs aux bords) ou [[Neumann]] (gradients aux bords) ou de [[Robin]] (relation gradient/valeurs sur le bord). Il s'agit donc avant tout de la résolution approchée d'un problème, 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 et où une discrétisation permet de trouver une solution approchée. 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 [[solution]]s, * [[unicité]] de la solution, * [[stabilité]], * [[convergence]], * et bien sûr: [[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. == Méthode des éléments finis == === 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 bords 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 quelconques : les éléments finis. Ces fragments peuvent par exemple être de forme quelconque mais doivent former un [[pavage]] de l'espace considéré. Usuellement les éléments finis sont de forme triangulaire ou rectangulaire. Ce maillage est associé à une base fonctionnelle ''b'', sur laquelle on projette la solution ''u''. Dans ce cadre, et moyennant une attention toute particulière au problème des bords du domaine, on obtient une formulation algébrique dite discrétisation du problème initial. La solution de ce problème algébrique, si elle existe et est unique, 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. === Dimensions === On développe ici la méthode des éléments finis en deux dimensions à valeurs réelles. On suppose que les équations étudiées sont des équations différentielles d'ordre deux. La méthode est généralisable à des cadres d'espaces de dimension différente ou pour des équations aux dérivées partielles d'ordre supérieur : * on traite ici le cas d'une solution réelle à une EDP, les cas où la dimension de la solution serait plus grande se traitent de façon similaire mais nécessitent des écritures plus complètes ; les cas les plus couramment rencontrés 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) ; * les degrés 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, c'est-à-dire des dérivées partielles des variables de départ (exemple 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 variationnelles afin d'obtenir des ordres plus faibles — en tout cas lorsque le degré de dérivation est pair. 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ésout rarement des problèmes en dimension supérieures à 3 — y compris des problèmes de dynamique en espace à 3 dimension qui pourraient être traité en quatre dimensions mais sont traités en réalité avec une méthode mixte éléments finis « en espace »/différences finies « en temps ». === Cadre algébrique, analytique et topologique === Soit un domaine (ouvert borné et connexe) Ω de <math>\mathbb{R}^2</math>, de bord δΩ, et d'adhérence (compacte) {{surligner|Ω}}. Pour simplifier les représentations, on suppose le bord polygonal. Soient les fonctions de {{surligner|Ω}} dans <math>\mathbb{R}</math> différentiables sur {{surligner|Ω}} (compact) et deux fois différentiables sur Ω (ouvert). De telles fonctions sont continues et différentiables sur le bord du compact. Soit V(Ω) l'ensemble de ces fonctions (V est un espace vectoriel de dimension infinie et V<sub>0</sub> est le sous-espace vectoriel de fonctions de V nulles sur le bord δΩ). Soient les applications continues sur {{surligner|Ω}} et différentiables sur Ω, de carré sommables sur {{surligner|Ω}} et de gradient de carré sommable sur {{surligner|Ω}} (ou de dérivées partielles de carré sommable, ce qui revient au même avec le support de dimension finie). Nommons cet espace <math>\mathcal{H}_1(\Omega)</math>. Cet espace est un [[espace de Sobolev]]. On dote cet espace vectoriel d'un produit scalaire issu de celui de L<sup>2</sup> tel que si (''u'', ''v''&thinsp;) appartiennent à cet espace alors le produit scalaire de ''u'' et ''v'' est : : <math>< u | v >_{ \mathcal{H}_1(\Omega) } = \int_\Omega{ ( \nabla u \cdot \nabla v + k^2 u v ) d\omega }</math> On note <math>\mathcal{H}_1^0(\Omega)</math> le sous espace vectoriel de <math>\mathcal{H}_1(\Omega)</math> dont les fonctions sont nulles sur le bord δΩ. L'opérateur « ⋅ » est un produit scalaire sur l'espace <math>\mathbb{R}^2</math> === Cas organique === ==== Hypothèses ==== On considère ƒ une fonction continue sur {{surligner|&Omega;}} de carré sommable et ''u'' la solution de l'équation aux dérivées partielles suivante sur Ω : : <math>- \Delta u + k^2 u = - \frac{\partial^2 u}{\partial x^2} - \frac{\partial^2 u}{\partial y^2} + k^2 u = f</math> Avec la condition au bord ''u'' = 0 sur δΩ. Ceci peut également se réécrire ''u'' ∈ V<sub>0</sub>. Cette condition au bord s'appelle la [[Conditions aux limites de Dirichlet|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<sub>0</sub> 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 : : <math>- \int_\Omega{ v \Delta u d\omega } + k^2 \int_\Omega{ v u d\omega } = \int_\Omega{v f d\omega }\,</math>. On utilise pour le premier terme une [[intégration par parties]] <ref>Pour certains cas limites lorsqu'il n'y a pas suffisamment de régularité et que l'intégration par partie n'a pas de sens on utilise aussi une formulation adaptée du [[théorème de Stokes]] mais il faut alors définir la dérivée faible puis la divergence faible ce qui entraine des complications supplémentaires.</ref> : : <math>- \int_\Omega{ v \Delta u d\omega } = - \int_{\partial \Omega }{ \frac{ \partial u }{ \partial n } v ds } + \int_\Omega{ ( \nabla u \cdot \nabla v ) d\omega }</math>. Dans cette formulation, ''v'' est nulle sur le bord (''v'' ∈ V<sub>0</sub>) ce qui permet d'obtenir la [[formulation faible]] du problème : : <math>\int_\Omega{ \nabla u \cdot \nabla v d\omega } + k^2 \int_\Omega{ v u d\omega } = \int_\Omega{ v f d\omega }</math>. 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élicat dans les cas limites où l'ouvert Ω n'est pas suffisamment régulier (par exemple s'il y a des [[point singulier|points singuliers]]) ou si ƒ n'est pas suffisamment dérivable (si l'on ne suppose pas que ƒ est au moins <math>\mathcal{L}_2</math>). Il faut alors souvent se ramener à une étude au cas par cas et rien ne dit que la formulation faible aura les mêmes solution que l'équation de départ. Dans la majorité des problèmes physiques la solution est souvent <math>\mathcal{C}_\infty</math> et l'on ne se pose pas de tels problèmes. Néanmoins il faut noter que pour des domaines avec des [[point singulier|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. === Notations et cadre général=== Pour plus de généralité et pour rendre la suite plus lisible on utilisera les notations suivantes : : <math>a(u,v) = \int_\Omega{ \left(\nabla u . \nabla v + k^2 u v \right) d\omega }</math> avec ''a'' un opérateur bilinéaire symétrique (de V<sup>2</sup> dans <math>\mathbb{R}</math>) ; : <math>\mathcal{L}(v) = \int_\Omega{ f v d\omega }</math> avec <math>\mathcal{L}</math> un opérateur linéaire (de V dans <math>\mathbb{R}</math>). On peut résoudre par la méthode des éléments finis toute équation aux dérivées partielles dont la forme faible se met sous la forme : <math>a( u, v ) = \mathcal{L}(v)</math> avec ''a'' un opérateur bilinéaire [[Coercivité|coercive]] continu selon la norme <math>\mathcal{H}_1^0</math> (cf. [[espace de Sobolev]]) et L un opérateur linéaire continue également selon la norme continu selon la norme <math>\mathcal{H}_1^0</math>. Avec ces notations le problème se reformule ainsi : : <math>\forall v \in V_0, a( u, v ) = \mathcal{L}(v)</math> On montre aussi que ''u'' est la solution unique du problème d'optimisation de la fonctionnelle suivante : : <math>\forall v \in V_0, \mathcal{J}(v) = \frac{1}{2} a( v, v ) - \mathcal{L}(v)</math>. Cette égalité peut avoir un sens physique notamment du point de vue de l'énergie pour certaines équations physiques et peut servir à montrer l'existence et l'unicité de la solution grâce aux propriétés de ''a'' et de L (linéarité, coercivité, …). === Choix d'un maillage et Discrétisation === =====Choix d'un maillage===== [[Image:Finite element triangulation.svg|right|thumb|Un exemple de maillage triangulaire]] 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êts (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. =====Fonctions de base===== [[Image:Finite element method 1D illustration2.png|right|thumb|Fonction de base en dimension 1. Les x<sub>i</sub> 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 des bases 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 [[Interpolation lagrangienne|polynômes de Lagrange]]. Dans cette méthode les fonctions de base valent 1 à un nœud du maillage et 0 a tout 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 à coté. Il y a autant de fonctions de base par éléments que de nombre 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 bases 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 de [[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 ====== En 2D * triangles de degré 1, (triangles à 3 nœuds, fonctions linéaires) * triangles de degré 2 (triangles à 6 nœuds, polynômes de degré 2) * quadrangles de degré 1 (carrés à quatre nœuds, fonctions linéaires) * quadrangles de degré 2 (carrés à 8 ou 9 nœuds, polynômes de degré 2) En 3D * tétraèdre de degré 1, (quatre nœuds, linéaires) * cube de degré 1, (huit nœuds, linéaire) [[Image:Fonctionshermite.jpg|right|thumb|Les deux fonction de base de Hermite associés au nœud 0 avec les nœuds voisins en +1 et -1 en dimension 1]] ===== Discrétisation ===== Soit le maillage <math>\mathcal{M}</math> et la base <math>b=(e_1...e_n)</math> associée. Puisque la condition de Dirichlet impose des fonctions nulles aux bords, on utilise uniquement la sous-base <math>b</math> limitée aux points intérieurs de <math>\Omega</math>. On cherche la solution <math>\bar{u}</math> du problème discrétisé ainsi : <math> \bar{u} \in V_n^0 | \forall v \in V_n^0, a( \bar{u}, v ) = \mathcal{L}(v)</math> 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 <math>\bar{u}</math> dans la base des <math>e_i</math> intérieurs, en composantes <math>u_i</math>, on obtient : <math> \forall j \in [1,...,n] \sum_{i =1}^n { u_i a( e_i, e_j ) } = \mathcal{L}(e_j)</math> L'idée est que quand le maillage se ressère et que le nombre de fonctions de base n tend vers l'infini (et que l'espace engendré par cette base <math> V_n^0</math> croit vers V<sub>0</sub>), les solutions ''u<sub>n</sub>'' devront converger vers la solution ''u'' de l'équation aux dérivées partielles de départ. ===== Éventuelle deuxième discrétisation ===== Dans certains problèmes physiques, il peut être intéressant de discrétiser une deuxième fois. Cette seconde discrétisation n'est pas nécessaire pour la méthode des éléments finis. Souvent on a comme expression de <math>\mathcal{L}</math>: <math> \mathcal{L}(e_j)=\int_\Omega f\cdot e_j </math> On projette alors f sur la base b. On obtient: :<math> f_n = \sum_{k=1}^n f_k e_k. </math> et on approche <math>\mathcal{L}(e_j)= \int_\Omega f\cdot e_j </math> par <math> \sum_{k=1}^n f_k \int_\Omega e_k e_j,\; j=1,\dots,n.</math> Le problème est d'obtenir ensuite une projection <math>f_n</math> acceptable sachant qu'il n'y a pas nécessairement de produit scalaire associé à la base qui permette de projeter de façon efficace. Dans les deux exemples de bases donnés plus haut, cette projection est aisée. Dans le cas des éléments finis de [[Lagrange]], la projection sur la fonction e<sub>i</sub> est donnée par la valeur en x<sub>i</sub>, dans le cas des éléments de Hermite, c'est la valeur de la fonction ainsi que de sa dérivée qui permettent d'obtenir la projection. Pour d'autres bases, la projection peut être plus compliquée. === Problème sous forme matricielle === Si l'on note: * la matrice <math>A</math> ayant pour composantes les <math>a( e_i, e_j )</math>, * le vecteur <math>U</math> ayant pour composantes les <math>u_i</math> qui sont les coordonnées de la solution approché sur la base b * le vecteur <math>B</math> pour composantes les <math>\mathcal{L}( e_j )</math> alors ce problème revient à résoudre l'équation linéaire de <math>n</math> équations à <math>n</math> inconnues : <math>\mathbf{A} U = B</math> La matrice A est appelée matrice de rigidité par analogie avec certain problèmes de mécanique des solides. <math>A</math> est par construction symétrique, et puisque <math>a</math> est [[Coercivité|coercitive]], alors <math>A</math> est symétrique, définie positive donc inversible. On obtient donc l'existence et l'unicité de <math>U = A^{-1} B</math>. grâce aux coordonnées de <math>\bar{u}</math> sur la base b on peut alors construire la solution approchée <math>\bar{u}</math>. Quand le maillage se resserre cette solution approchée va tendre vers la vraie solution de l'équation aux dérivées partielles de départ. Pour le cas avec une deuxième discrétisation de <math>\mathcal{L}( e_j )</math> on obtient: :<math> \mathbf{A}U=\mathrm{M}f </math> où M est appelée la matrice de masse et contient les <math>\int_\Omega e_i \cdot e_j </math>. f est un vecteur contenant les coordonnées de f dans la base. La méthode est alors la même qu'avec une seule discrétisation puisque A vérifie les mêmes propriétés. Cette méthode peut parfois être préférée quand on peut obtenir de façon simple la projection de f sur la base et la matrice M. === Algorithme === La méthode des éléments finis doit être conduite ainsi # On calcule la matrice de rigidité ''A'' # On détermine le membre de droite, en calculant les termes <math>\mathcal{L}(e_j)</math> ou alors par l'intermédiaire de la matrice de masse. # On résout le problème ''AU''=''B'' ou le problème ''AU''=''Mf'' suivant le niveau de discrétisation choisi. U est alors donné par <math>U = A^{-1} B</math>. 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 pour A. C'est l'étape la plus consommatrice en termes de puissance de calcul et l'efficacité de la méthode en termes de temps de calcul se joue principalement sur cette étape. # On peut écrire <math>\bar{u}</math> grâce au vecteur U qui contient les coordonnées de <math>\bar{u}</math> sur la base b et obtenir une solution approchée au problème. ====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 <math>\bar{\Omega}</math>, et que la [[Conditions aux limites de Neumann|condition de Neumann]] <math>\frac{\partial u}{\partial n}(\bar{\Omega}) = 0</math> est vérifiée. Si la fonction est supposée différentiable au bord <math>D_u(X).n = 0, \forall X \in \bar{\Omega}</math>, voire si elle admet un gradient <math><\nabla u(X)| n> = 0, \forall X \in \bar{\Omega}</math>. 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 <math>\int_\Omega{\frac{\partial u}{\partial n} v d\omega}=0</math> 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. == Exemples issus de problèmes physiques == ==== Historique ==== * Analyse des structures née vers [[1850]] ** [[Résistance des matériaux|RDM]] (Résistance des matériaux) ⇒ calculs « manuels »<br />&nbsp;&nbsp;&nbsp;[[James Clerk Maxwell|Maxwell]], [[Carlo Alberto Castigliano|Castigliano]], [[Christian Otto Mohr|Mohr]] ** Concept d'éléments finis né vers [[1940]]<br />&nbsp;&nbsp;&nbsp;[[Nathan Newmark|Newmark]], [[Hrenikoff]], [[Mc Henry]], [[Richard Courant|Courant]] ** Développement réel depuis 1960<br />&nbsp;&nbsp;&nbsp;Calcul numérique sur [[ordinateur]] ==== Domaines d'application ==== * [[Calcul de structures]], [[étude des contacts]], [[électricité]], [[électromagnétisme]], [[hydraulique]], [[aérodynamique]], ... ==== 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ètre]]s. # La résolution du problème global, obtenu par assemblage des éléments, suit les règles qui régissent les structures discrètes. ===== 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) === Exemple de problème discret : un réseau électrique === [[Image:Pb discret.png|center]] <ol> <li> Équation locale du composant e :<br /> :<math>\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} \,</math><br /> :::soit <math>\left\{ I^e \right\} = \left[ K^e \right] \left\{ U^e \right\} \,</math> </li> <li> On écrit :<br /> :* La continuité des potentiels en chaque connexion <math> i \,</math> :* L'équilibre des courants à chaque connexion :* L'adjonction des courants externes <math> P^i \,</math> </li> <li> On obtient l'équation globale du système assemblé :<br /> :<math>P^i = \sum_{j=1}^m \sum_{e=1}^m K_{ij}^e U_j^e \,</math> :::soit <math>\left\{ P \right\} = \left[ K \right] \left\{ U \right\} \,</math> </li> </ol> === Définition d'un élément fini === {| align=center border=0 |valign="middle" align="center" width="45%" | * <u>Problème modélisé</u> : |valign="middle" align="center" width="10%" | &nbsp; |valign="middle" align="center" width="45%" | * <u>Problème simplifié</u> : |- | [[Image:milieu continu.png]] | &nbsp; | [[Image:milieu discretise.png]] |} :En calcul de structures, un élément fini est caractérisé par deux matrices : :* La matrice de raideur <math>\left[ K \right] \,</math> :* La matrice de masse <math>\left[ M \right] \,</math> === Formulation d'un élément fini === ==== Définitions et notations ==== * <u>Déplacement</u> : <math>\left\{ U \right\} \,</math> :C'est un vecteur dont chaque composante est également appelée degré de liberté :* 3 ddl de translation : <math> U_x,U_y,U_z \,</math> :* 3 ddl de rotation : <math> \theta_x,\theta_y,\theta_z \,</math> * <u>Déformation</u> : (voir [[Tenseur des déformations]]) <math>\left[ \epsilon \right] \,</math> :C'est le rapport de l'allongement à la longueur initiale.<br /> En petites déformations, on a <math> (A'B')^2 = \left(\ dx + {\partial u \over \partial x} \ dx\right)^2 + \left({\partial v \over \partial x} \ dx\right)^2 \,</math>[[Image:deformation EF.png|right]] :<math> \epsilon_x = {{A'B'-AB} \over AB} \,</math><br /> Comme <math> AB=\ dx \,</math>, on a <math> A'B'=(1+\epsilon_x)\ dx \,</math><br /> :<math> \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 \,</math><br /> On néglige les termes d'ordre 2 :<br /> :<math> \epsilon_x = {\partial u \over \partial x} \,</math><br /> <u>Remarque</u> : <math> \epsilon \,</math> est sans dimension * <u>Contrainte</u> : (voir [[Tenseur des contraintes]]) <math>\left[ \sigma \right] \,</math> :*Elle représente les efforts internes qui s'appliquent dans la structure.<br /> ::<math>\sigma_{11},\sigma_{22},\sigma_{33} \,</math> : contraintes normales<br /> ::<math>\sigma_{12},\sigma_{13},\sigma_{23} \,</math> : contraintes de cisaillement<br /> :* Une contrainte est homogène à une pression (N/m²) :* Il existe un système de coordonnées dans lequel <math>\left[ \sigma \right] \,</math> est une matrice diagonale : <math>\left[ \sigma \right] = \begin{bmatrix}\sigma_X & 0 & 0 \\ 0 & \sigma_Y & 0 \\ 0 & 0 & \sigma_Z\end{bmatrix} \,</math> ::<math>\sigma_X,\sigma_Y,\sigma_Z \,</math>sont appelées les '''contraintes principales''' * <u>Relations Contraintes-Déformations</u> : (voir [[Loi de Hooke]]) ≡ lois de comportement :<math>\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}\,</math> ::<math>E \,</math>: module de Young (N/m²) ::<math>\nu \,</math> : coefficient de Poisson (sans dimension) :*On utilise parfois le module de cisaillement : <math>G = {E \over {2(1+\nu)}} \,</math><br /> :*Pour un matériau [[isotrope]], il n'y a que 2 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 : <math>\left\{ \sigma \right\} = \left[D\right]\left\{ \epsilon \right\} \,</math> ::<math>\left[D\right] \,</math> est appelée matrice d'élasticité du matériau. * <u>Énergie de déformation</u> : <math>W \,</math> :<math>W = {1 \over 2}\int_{v}\sigma.\epsilon.\ dv \,</math> * <u>Travail d'une force</u> : : C'est le produit de la force par le déplacement de son point d'application : <math>\tau_F = {1 \over 2}F.U_F \,</math> * <u>Moment</u> : : C'est une force appliquée sur un ddl de type rotation ==== Équations fondamentales ==== <u>Équations d'équilibre local</u> :<br /> :<math>{\partial \sigma_x \over \partial x} + {\partial \sigma_{xy} \over \partial y} + {\partial \sigma_{xz} \over \partial z} + F_x = 0\,</math>[[Image:Elasticite lineaire|right]] :<math>{\partial \sigma_y \over \partial y} + {\partial \sigma_{xy} \over \partial x} + {\partial \sigma_{yz} \over \partial z} + F_y = 0\,</math><br /> :<math>{\partial \sigma_z \over \partial z} + {\partial \sigma_{zx} \over \partial x} + {\partial \sigma_{zy} \over \partial y} + F_z = 0\,</math> <u>Relations déformations-déplacements</u> :<br /> :<math>{\epsilon_x= {\partial u \over \partial x}} \qquad {\epsilon_y= {\partial v \over \partial y}} \qquad {\epsilon_z= {\partial w \over \partial z}} \,</math><br /> :<math>\gamma_{xz}={\partial w \over \partial x} + {\partial u \over \partial z} \,</math><br /> :<math>\gamma_{xy}={\partial u \over \partial y} + {\partial v \over \partial x} \,</math><br /> :<math>\gamma_{yz}={\partial w \over \partial y} + {\partial v \over \partial z} \,</math><br /> ::Symboliquement, on écrit <math>\{\epsilon\} = [S]\{U\} \,</math> <u>Signification du coefficient de Poisson</u> :<br /> [[Image:Poisson coef.png|center]] Si on applique au barreau une contrainte <math>\sigma_x \,</math>, on observe un rétrécissement dans la direction y correspondant à une déformation <math>-\nu \sigma_x \over E \,</math><br /> <u>Quelques valeurs usuelles :</u> * Acier : E = 2,1E11 N/m² = 210 000 MPa &nbsp;&nbsp;&nbsp;&nbsp;<math>\nu \,</math> = 0,3 * Aluminium : E = 7E10 N/m² = 70 000 MPa &nbsp;&nbsp;&nbsp;&nbsp;<math>\nu \,</math> = 0,3 <u>Remarques</u> : On a toujours -1 <math>\le \nu \le \,</math>0,5 <br /> Quand <math>\nu \to \,</math>0.5, le matériau est dit incompressible. ==== Exemple de formulation : Barre en traction ==== [[Image:Barre traction.png|center]] On suppose que le déplacement en tout point de la barre est donné par un [[polynôme]] du 1<sup>er</sup> degré : <math> u(x) = a_1 + a_2 x \,</math><br /> :On a <math> u(0) = u_1 \,</math> et <math> u(L) = u_2 \,</math><br /> :d'où <math> u(x) = u_1\left(1-{x \over L} \right) + {x \over L} u_2 \,</math><br /> qu'on écrit symboliquement : <math> u(x) = [N(x)]\{U\} \,</math> avec <math>\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. \,</math><br /> On en déduit :<br /> :<math> \epsilon = \left[-{1 \over L} ; {1 \over L} \right]\begin{Bmatrix}u_1 \\ u_2\end{Bmatrix} = [B]\{U\} \,</math><br /> :<math> \sigma = E\left[-{1 \over L} ; {1 \over L} \right]\begin{Bmatrix}u_1 \\ u_2\end{Bmatrix} = [D]\epsilon \,</math><br /> D'autre part, on a par définition :<br /> :<math> \begin{Bmatrix}F_1 \\ F_2 \end{Bmatrix} = \begin{Bmatrix}-\sigma S \\ \sigma S \end{Bmatrix} \,</math> où S est l'aire de la section de la barre. On pose : :<math> \{F\} = \begin{Bmatrix}F_1 \\ F_2 \end{Bmatrix} = [A]\sigma \qquad \mbox {avec} \qquad \{A\} = S \begin{Bmatrix}-1 \\ 1 \end{Bmatrix} \,</math><br /> On obtient finalement :<br /> :<math> \{F\} = [A][D][B]\{U\} \,</math><br /> Soit une relation du type :<br /> :<math> \{F\} = [K]\{U\} \qquad \mbox {avec} \qquad [K] = [A][D][B] \,</math><br /> En explicitant : :<math> [K] = {ES \over L} \begin{bmatrix}1 & -1 \\ -1 & 1 \end{bmatrix} \,</math><br /> On voit que la matrice de rigidité se calcule comme le produit de 3 matrices :<br /> :<math> [B] \,</math> : Transformation des déplacements aux déformations<br /> :<math> [D] \,</math> : Matrice d'élasticité du matériau<br /> :<math> [A] \,</math> : Transformation des contraintes en forces ==== Formulation générale (méthode directe) ==== : La démarche est la suivante : * On exprime le déplacement <math>\{u(x)\}\,</math> en tout point de l'élément en fonction des déplacements aux nœuds <math>\{u(x)\} = [N(x)]\{U\}\,</math> * On exprime les déformations en fonction des déplacements <math>\{\epsilon(x)\}=[S]\{u(x)\}\,</math> :d'où <math>\{\epsilon(x)\}=[S]\{N(x)\}\{U\}=[B(x)]\{U\}\,</math> * On écrit la loi de comportement du matériau qui relie les contraintes aux déformations : <math>\{\sigma(x)\}=[D]\{\epsilon(x)\}\,</math> * On écrit que le travail des forces externes appliquées à la structure pour un déplacement virtuel <math>\delta U \,</math> est égal au travail interne des contraintes pour ce même déplacement : :<math>\{\delta U\}^T\{F\}=\int_V \{\delta \epsilon\}^T\{\sigma\} dv \,</math><br /> :En explicitant, on a :<br /> :<math>\{\delta U\}^T\{F\}=\{\delta U\}^T \left (\int_V [B]^T[D][B] dv \right ) \{U\} \,</math> :Comme cette relation est vraie pour tout déplacement virtuel, on en déduit :<br /> :<math>\{F\}=[K]\{U\} \,</math><br /> :avec <math>[K] \,</math> sous sa forme plus générale : <math>[K]=\int_V [B]^T[D][B] dv \,</math> <u>Remarques</u> :<br /> * La relation ci-dessus montre que <math>[K] \,</math> est symétrique. * Le terme courant <math>K_{ij} \,</math> de la matrice correspond à la force qui s'exerce sur le nœud <math>j \,</math> lorsqu'on impose un déplacement unitaire du nœud <math>i \,</math>. La symétrie de <math>[K] \,</math> qui s'écrit <math>K_{ij} = K_{ji} \,</math> correspond mécaniquement au [[Théorème de Maxwell-Betti|théorème de réciprocité de Maxwell-Betti]]. <u>Qu'est-ce qui va différencier les différents types d'éléments finis</u> ?<br /> * 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'élasticité 3D]] <u>Remarques</u> :<br /> 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 [[principe des puissances virtuelles|puissances virtuelles]] :* La minimisation de l'[[énergie potentielle]] '''Toutes ces approches sont équivalentes et aboutissent à la construction de la même matrice de rigidité.''' ===É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. === É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ôme]]s. ::*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.<br /> Les fonctions N(x) sont appelées ''fonctions de forme'' ou ''fonctions d'interpolation de l'élément''. === Les éléments isoparamétriques === ::<math>[K]=\int_V [B]^T[D][B] dv \,</math><br /> ==== <u>Problème</u> ==== * Soit on construit <math>[K] \,</math> 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 <math>[K] \,</math> à chaque fois. ==== <u>D'où l'idée</u> ==== 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 <u>mais également sa géométrie</u>. De plus, on va travailler en coordonnées locales. ==== <u>Interpolation de la géométrie</u> ==== {| border=0 cellspacing=0 cellpadding=0 | valign="middle" align="left" width="33%" | <math>\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} \,</math><br /> Idem pour les autres coordonnées. | valign="middle" align="center" width="34%" | &nbsp; | valign="middle" align="left" width="33%" | [[Image:EL isoparam champ.png]] |} ==== <u>Coordonnées locales (cas 2D) </u> ==== [[Image:EL isoparam coord.png|center]] ==== <u>Élément isoparamétrique</u> ==== 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.<br /> <math>\tilde N(x) = N(x) \,</math> ==== <u>Autres classes d'éléments</u> ==== [[Image:EL super param.png|center]] ==== <u>Évaluation de K</u> ==== La forme générale s'écrit :<br /> <math>[K]=\int_x \int_y \int_z G(x,y,z) dxdydz \,</math><br /> On passe en variables locales <math>\xi,\eta,\zeta \,</math><br /> On a <math>dxdydz = \det J\ d\xi d\eta d\zeta \,</math><br /> :<math>J \,</math> s'appelle la matrice Jacobienne.<br /> On est alors amené à calculer des intégrales du type :<br /> <math>\int_{-1}^{+1} \int_{-1}^{+1} \int_{-1}^{+1} G(\xi ,\eta ,\zeta) \det J\ d\xi d\eta d\zeta \,</math><br /> ==== <u>Bénéfice de l'approche</u> ==== On s'est ramené à un domaine d'intégration simple et invariant pour lequel on peut appliquer les formules de quadrature de gauss :<br /> <math>\int_{-1}^{+1} f(x) dx = \sum_{k=1}^n H_k f(a_k) \,</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; les <math>a_k \,</math> et <math>H_k \,</math>étant <u>tabulés</u>.<br /> Les <math>a_k \,</math> sont appelés <u>points d'intégration</u> de l'élément ou encore <u>points de [[Carl Friedrich Gauss|Gauss]]</u> de l'élément. === Cas particulier : les éléments axisymétriques === ::[[Image:EL axisymetrique.png|center]] Décomposition en série de Fourier :<br /> :<math>u(r,\theta,z) = \sum_n u_n^s(r,z)cos (n\theta) + \sum_n u_n^s(r,z)sin (n\theta) \,</math><br /> L'axisymétrie correspond à la restriction <math>n=0 \,</math> de cette décomposition.<br /> :<math>u(r,\theta,z) = u_n^s (r,z) \,</math><br /> <u>Remarque importante</u> : Pour utiliser ce type d'élément, le problème doit être <u>globalement</u> axisymétrique :<br /> * la géométrie * les conditions limites * le chargement === Processus de calcul (cas statique) === # Maillage # Construction de la matrice de raideur de chaque élément <math>[K^e] \,</math> # Assemblage de la matrice globale <math>[K] \,</math> # Construction du vecteur chargement <math>\{F\} \,</math> # Élimination de certains ddl (si besoin) # Résolution : <math>\{U\} = [K^{-1}]\{F\} \,</math> # Calcul des quantités dérivées de <math>\{U\} \,</math> ::<math>\{\epsilon^e\} = [B^e] \{U^e\} \,</math> ::<math>\{\sigma^e\} = [D^e] \{\epsilon^e\} \,</math> ::<math>W^e = {1 \over 2} \{\epsilon^e\}^T \{\sigma^e\} \,</math> ::::etc. == Références== <references/> *G. Allaire and A. Craig: Numerical Analysis and Optimization:An Introduction to Mathematical Modelling and Numerical Simulation * 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''' * J. Chaskalovic, "Méthode des éléments finis pour les sciences de l'ingénieur", (2004), Ed. Lavoisier. '''ISBN-10 : 2-7430-0708-7''' * 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''' ==Quelques exemples de logiciels d'éléments finis appliqués à la mécanique des structures== Quelques exemples de logiciels généralistes utilisant la méthode des éléments finis en mécanique des structures * [[ABAQUS]]: logiciel américain * [[CosmosWorks]] : Logiciel Franco-Américain appartenant à SolidWorks qui lui même appartient à Dassault [http://www.cosmosm.com/ CosmosWorks] * [[:en:ANSYS|ANSYS]]: logiciel américain * [[CAST3M]]: logiciel français mis à disposition gratuitement pour l'enseignement et la recherche [http://www-cast3m.cea.fr CASTEM] * [[SYSTUS]]: logiciel français, traite des calculs mécaniques, thermiques et thermo-mécaniques linéaires et non-linéaires * [[SYSWELD]]: logiciel français, basé sur SYSTUS et permettant le calcul de traitement thermique, thermo-chimique, soudage, avec couplage mécanique/thermique/métallurgique voir éléctromagnétisme (trempe par induction) * [[Code Aster]]: logiciel libre français [http://www.code-aster.org Aster] * [[Nastran]]: logiciel américain * [[PERMAS]]: logiciel allemand [http://www.intes.de PERMAS] * [[SAMCEF]]: logiciel belge [http://www.samcef.com SAMCEF] * [[Morfeo]]: logiciel belge [http://www.cenaero.be Morfeo] * [[JMAG]]: logiciel Japonais (distribué en Europe par Powersys [http://www.powersys.fr Powersys]) permet un couplage entre les analyses électromagnétiques et structurelles. Quelques exemples de logiciels de dynamique non linéaire (explicite) des structures utilisant la méthode des éléments finis en mécanique des structures * [[ABAQUS]]: logiciel américain * [[EuroPlexus]]: logiciel français [http://europlexus.jrc.it EuroPlexus] * [[:en:LS-DYNA|LS DYNA]]: logiciel américain * [[PAM]]: logiciel français [http://www.esi-group.com/SimulationSoftware/ PAM] * [[Radioss]]: logiciel français ==Voir aussi== * [[Méthode des différences finies]] * [[Méthode des volumes finis]] * [[Théorème de Lax-Milgram]] * [[Méthode des Eléments Finis Etendus]] * {{fr}} [http://norma.mas.ecp.fr/wikimas/ScientificComputingSoftware Liens vers des codes Eléments Finis téléchargeables] * {{en}} [http://www.uni-ulm.de/uni/intgruppen/fem/ La méthode des éléments finis en génie biomedicale et en biomécanique] {{Portail|mathématiques|physique|Génie mécanique}} [[Catégorie:Mécanique des milieux continus]] [[Catégorie:Équations différentielles numériques]] [[ca:Mètode d'elements finits]] [[de:Finite-Elemente-Methode]] [[en:Finite element method]] [[es:Elementos finitos]] [[it:Metodo degli elementi finiti]] [[ja:有限要素法]] [[nl:Eindige-elementenmethode]] [[pl:Metoda elemetów skończonych]] [[ru:Метод конечных элементов]] [[sv:Finita elementmetoden]] [[zh:有限元分析]]