République Algérienne Démocratique et Populaire
Ministère de l'Enseignement Supérieur et de la Recherche
Scientifique.
Université Abbes Laghroure Khenchela
Faculté de des sciences et techniques
Département de Génie Mécanique
Cours
Méthode des Éléments Finis
Présenté par
Noureddine SID
2024/2025
1
I.1 Introduction
La méthode des éléments finis (MEF) est une technique numérique largement utilisée
en ingénierie et en sciences pour résoudre des problèmes complexes de calcul impliquant des
équations différentielles et des dérivées partielles. Elle est particulièrement efficace pour les
problèmes liés à la mécanique des structures, à la mécanique des fluides, à la thermique, à
l'électromagnétisme, et bien plus encore.
I.2 Historique
Le mathématicien Courant a publié en 1943 un exposé présenté en 1941 traitant de la
solution du problème de torsion de Saint Venant à l’aide d’une méthode connue de nos jours
comme la méthode des éléments finis. Le travail de Courant a été ignoré jusqu’à ce que des
ingénieurs développent une méthode identique. A cette époque cette méthode était inutilisable
car les ordinateurs digitaux n’existaient pas. Ces derniers ont été créés vers la fin des années
50, de même que le langage de programmation FORTRAN. Après la Deuxième Guerre
mondiale l’industrie aéronautique développe des méthodes matricielles permettant de traiter
des problèmes de structures assez complexes. En 1955, Argyris présente une approche unifiée
des méthodes des déplacements et des forces puis, l’année suivante, Turner et Clough publient
une présentation systématique de la méthode des déplacements. Ces deux publications
représentent le vrai début de la méthode par éléments finis. Le terme “méthode par éléments
finis” a été défini par Clough en 1960. Dès lors, on assiste au développement de nombreux
éléments tels que les membranes, les plaques, les coques et les éléments de volume et la
méthode s’étend de la statique linéaire au transfert thermique, à la dynamique linéaire, plus
récemment, à l’analyse non linéaire. Les grands programmes d’analyse par les méthodes des
éléments finis ont été développés à partir des années 60. On citera pour exemple Nastran,
développé par la NASA, Ansys et Aska. Très rapidement ces programmes présentent trois
parties distinctes, le pré-processeur pour l’introduction des données, le solveur pour le calcul et
le post-processeur pour le traitement des résultats. Dès les années 80 ces logiciels sont
disponibles sur des micro-ordinateurs. Dix publications concernant les éléments finis ont été
présentées en 1961, 134 en 1966, 844 en 1971. En 1976 le nombre cumulé de publications était
de 7000, de 20’000 en 1986 et de 40’000 au milieu des années 90.
2
I.3 Généralités
I.3.1 Processus d’analyse
De façon générale, les différentes étapes d’analyse d’un problème physique s’organisent
suivant le processus schématisé par la figure I.1. Nous partons d’un problème physique. Le
cadre précis de l’étude est défini par les hypothèses simplificatrices qui permettent de
déterminer le modèle mathématique approprié. La difficulté pour l’ingénieur est de savoir
choisir parmi les lois de la physique, celles dont les équations traduiront avec la précision
voulue la réalité du problème physique.
Figure I.1 : Processus d’analyse utilisant un modèle numérique.
Un bon choix doit donner une réponse acceptable pour des efforts de mise en œuvre non
prohibitifs.
En résumé, les questions essentielles auxquelles l’ingénieur devra répondre, s’il veut
effectuer une analyse par un modèle numérique dans de bonnes conditions, sont les suivantes :
– quel modèle mathématique utiliser ?
3
– quel modèle numérique faut-il lui associer ?
– quelle est l’erreur d’approximation commise ?
– quelle est l’erreur numérique commise ?
– peut-on améliorer le modèle numérique ?
– faut-il changer le modèle mathématique ? Etc.
Qu’est-ce qu’un modèle ? La figure I.2 illustre sur un exemple mécanique simple trois
modélisations envisageables. Chacune correspond à modèle mathématique différent mais
quelle est la bonne ? Le choix du modèle mathématique est un compromis entre le problème
posé à l’ingénieur « quelles grandeurs veut-on calculer et avec quelle précision ? » et les moyens
disponibles pour y répondre. En fait, les équations du modèle retenu sont soumises à un certain
nombre d’hypothèses basées sur les sciences de l’ingénieur et il faut connaître leur domaine de
validité pour pouvoir vérifier que la solution obtenue est satisfaisante.
Figure I.2 : Choix d’un modèle mathématique : dimensionnement statique d’un support
d’étagère.
Si le modèle mathématique n’admet pas de solution analytique, il est alors nécessaire de
chercher une solution approchée de ce modèle. Dès lors, la discrétisation du problème
correspond au choix d’un modèle numérique permettant de traiter les équations mathématiques.
4
Il est important de savoir distinguer et hiérarchiser les différents niveaux d’hypothèse utilisés
pour modéliser un phénomène physique. En effet, la solution exacte d’un modèle mathématique
qui ne correspond pas à la réalité physique est inutile.
I.4 Méthodes d’approximation
Pour discrétiser les modèles complexes de phénomènes physiques, l’ingénieur dispose
de méthodes d’approximation permettant de résoudre la plupart des problèmes pour lesquels il
n’existe pas de solution formelle.
Toutes les méthodes d’approximation ont un même objectif, à savoir remplacer un
problème mathématique défini sur un milieu continu (équations différentielles ou intégrales)
par un problème mathématique discret (équations matricielles) de dimension finie que l’on sait
résoudre numériquement.
La classification que nous proposons sur la figure I.3 n’est pas unique. Elle permet
simplement de distinguer la méthode, en fonction de la démarche utilisée pour obtenir une
forme intégrale. Il est important de noter qu’un problème physique peut être formulé de façon
équivalente en un système d’équations différentielles ou sous une formulation variationnelle.
Nous montrons par la suite comment passer de l’une à l’autre.
I.4.1 Méthode des résidus pondérés (ou annulation d’erreur) : elle utilise comme point de
départ les équations locales, équations différentielles définies sur l’intérieur du domaine, et les
conditions aux limites du problème définies sur la frontière du domaine ;
I.4.2 Méthodes variationnelles : le point de départ de ces méthodes est un principe variationnel
qui est une formulation mathématique du problème basée sur des considérations énergétiques.
La formulation obtenue dépend bien entendu des hypothèses de modélisation du problème
physique.
5
Figure I.3 : Vue synthétique des méthodes d’approximation.
I.5 Domaines d’application de la méthode des éléments finis
Les domaines d’application de la méthode des éléments finis sont très larges et très
varies. En général pour les champs d’application de cette méthode sont :
1. La mécanique des milieux continus solides
2. La mécanique des milieux continus fluides
3. L’analyse thermique
4. L’électromagnétisme et l’électrostatique
Dans le cadre de la mécanique des milieux continus solides, on peut classer les
problèmes de la manière suivante (figure I.2) :
6
Figure I.4 : Classification des problèmes de la mécanique.
I.6 Principe de la méthode des éléments finis
La méthode des éléments finis consiste a diviser le domaine physique à traiter en
plusieurs sous domaines appelés éléments finis a dimensions non infinitésimales. La solution
recherchée est remplacée dans chaque élément par une approximation avec des polynômes
simples et le domaine peut ensuite être reconstitue avec l’assemblage ou sommation de tous les
éléments.
Etape 1 : Formulation des équations gouvernantes et des conditions aux limites.
La majorité des problèmes d'ingénierie sont décrits par des équations différentielles aux
dérivées partielles associées a des conditions aux limites définies sur un domaine et son contour.
L'application de la MEF exige une réécriture de ces équations sous forme intégrale. La
formulation faible est souvent utilisée pour inclure les conditions aux limites.
Etape 2 : Division du domaine en sous domaines.
Cette étape consiste à discrétiser le domaine en éléments et calculer les connectivites de chacun
ainsi que les coordonnées de ses nœuds. Elle constitue ainsi la phase de préparation des données
géométriques.
Etape 3 : Approximation sur un élément.
Dans chaque élément la variable tel que le déplacement, la pression, la température, est
approximée par une simple fonction linéaire, polynomiale ou autre. Le degré du polynôme
d'interpolation est relie au nombre de nœuds de l'élément. L'approximation nodale est
appropriée. C'est dans cette étape que se fait la construction des matrices élémentaires.
7
Etape 4 : Assemblage et application des conditions aux limites.
Toutes les propriétés de l'élément (masse, rigidité,...) doivent être assemblées afin de former le
système algébrique pour les valeurs nodales des variables physiques. C'est à ce niveau qu'on
utilise les connectivites calculées à l'étape 2 pour construire les matrices globales à partir des
matrices élémentaires.
Etape 5 : Résolution du système global
Le système global peut être linéaire ou non linéaire. Il peut définir soit un problème d'équilibre,
de valeurs critiques ou de propagation. Le problème d’équilibre concerne les cas statiques et
les cas stationnaires. Dans un problème de valeurs critiques, on s’intéresse aux fréquences et
aux modes propres de vibrations du système physique étudie. Les problèmes de propagations,
concernent les cas transitoires dans lesquels sont déterminées les variations dans le temps des
variables physiques. Les méthodes d'intégration pas à pas conviennent mieux pour ce type de
problème. Les plus utilisées sont : méthode des différences finies centrales, méthode de
Newmark, méthode de Wilson. A ces méthodes doivent être associées des techniques d'itération
pour traiter le cas non linéaire. La plus fréquente est la méthode de Newton-Raphson.
I.7 Formes classique des éléments
Il existe plusieurs formes d’éléments classiques correspondant a des domaines a une, deux
ou trois dimensions. Chaque type d’élément est identifié par un nom précisant sa forme et par
le nombre de nœuds géométrique qui le composent.
Eléments à une dimension
· Eléments à deux dimensions triangulaires
8
Eléments à deux dimensions quadrangulaires
Eléments à trois dimensions tétraédriques
Eléments à trois dimensions hexaédriques (briques)
Figure I.5 : Formes classique des éléments.
9
I.8 Approximation nodale (fonctions de forme)
Partant d’hypothèses de petits déplacements et petites déformations, la mécanique du
solide a permis d’établir deux types de conditions régissant l’équilibre d’un corps : les équations
d’équilibre des forces et la compatibilité des déplacements.
Une des méthodes les plus utilisées pour les résoudre est la méthode dite des éléments
finis qui revient à remplacer le système continu par un système discret. Le solide est alors divisé
en un certain nombre de sous-domaines appelés éléments, dont l’assemblage permet la
reconstitution de la géométrie initiale. Chacun des éléments est relié à ces voisins par des nœuds
dont les degrés de liberté (DDL) constituent les inconnues du problème.
Figure I.6 : (a) Solide (Poutre en I) ; (b) Modèle éléments finis.
Considérant un champ de déplacement cinématiquement admissible sur l’élément, la
méthode consiste le plus souvent à approximer celui-ci au moyen d’une fonction polynomiale
formée d’un nombre fini de paramètres et à l’exprimer en fonction des déplacements nodaux
(les déplacements associés aux degrés de liberté). On aura donc :
Où les :
– Xi désignent les inconnues en déplacement associées aux degrés de liberté considérés,
– Ni (x,y,z) désignent les fonctions de forme de l’élément permettant d’obtenir les déplacements
de celui-ci en n’importe quel point de sa géométrie et ce toujours à partir des déplacements
nodaux.
À titre d’exemple, prenons un élément « barre » travaillant uniquement en traction ou
en compression. Le matériau utilisé étant supposé linéaire, la déformation est de fait constante
sur la hauteur de sa section. La fonction de déplacement u(x) est alors forcément linéaire. On a
donc :
10
Par ailleurs, on sait qu’en :
D’où
Soit sous une autre forme (figure I.7) :
Figure I.7 : Fonctions de forme d’un élément barre.
En regroupant toutes les fonctions de forme de l’élément au sein de la matrice [N] et
les déplacements nodaux dans le vecteur {qe}, le champ de déplacement s’exprime alors :
Avec {u} correspondant au vecteur des fonctions de déplacement applicables sur
l’élément : u(x, y, z), v(x, y, z), w(x, y, z), 𝛽 x(x, y, z), etc. Le degré du polynôme retenu pour
l’élément est donc intimement lié aux nombres de nœuds et de degrés de liberté par nœud. Par
exemple, un élément unidimensionnel à trois nœuds aura une fonction d’approximation
parabolique.
Généralement, les fonctions d’interpolation retenues sont linéaires, quadratiques et plus
rarement cubiques ce qui ne se révèle pas toujours conforme aux théories visées telles que celles
11
des poutres, des plaques ou des coques. En effet, les lignes et surface élastique découlant de ces
théories sont plutôt du troisième voire du quatrième degré. Une solution revient à discrétiser le
modèle de telle manière à reconstituer le champ réel de déplacement à partir de fonctions de
degré inférieur.
La qualité et l’efficacité d’une telle reconstitution dépendront donc d’une part du type
d’élément choisi, en d’autres termes de ses capacités, et d’autre part de la densité du maillage.
Dans le cas des éléments courants et toujours dans une hypothèse de linéarité, ce choix aura
pour conséquence de ne pas garantir la continuité du champ de contrainte d’où la recherche
d’une discrétisation adaptée permettant sa bonne reconstitution. De ce fait, la taille des éléments
utilisés sera généralement petite.
I.9 Résolution
La méthode de résolution s’inspire généralement de celle de Rayleigh – Ritz qui consiste
via une minimisation de l’énergie potentielle, en une recherche des termes d’une fonction
approximant le champ inconnu et satisfaisant aux conditions aux limites.
On définit l’énergie potentielle E comme étant la différence entre l’énergie de
déformation W et le travail T des forces de volume et de surface :
Le problème revient alors à rechercher un champ de déplacement minimisant E de manière à
caractériser son équilibre. L’étude d’un barreau en traction permet de décrire simplement la
démarche suivie.
I.10 Énergie de déformation élémentaire
Notant que sous une forme plus synthétique {ε} = []{u}, le vecteur des déformations
pourra être exprimé en fonction des déplacements nodaux à partir de :
avec []: opérateur dérivation dont les dimensions dépendent de celles de l’élément étudié.
12
Soit pour l’élément barre :
Le champ de contrainte étant obtenu à partir de la relation contrainte-déformation
{}= [H]{ ε }, l’énergie de déformation We associée à l’élément e peut être déduite en posant :
matrice de rigidité élémentaire de l’élément e.
Ceci donne pour l’élément barre :
I.11 Travail élémentaire des forces de volume et de surface
En traduisant les actions des éléments connectés à e sous la forme d’un vecteur forces
{ pe}, Te s’écrit :
13
I.12 Assemblage des matrices élémentaires
Les caractéristiques globales sont obtenues par assemblage des différentes quantités We
et Te. L’énergie potentielle totale de la structure est alors égale à la somme des différences entre
énergies de déformation et travaux des forces de chacun des éléments :
I.13 Changement de repères
Cependant, les éléments n’ayant pas forcément le même repère local, la sommation des
énergies de déformation et des travaux des forces extérieures ne peut être directement
appliquée. Il est donc nécessaire de les exprimer dans un référentiel unique appelé repère global.
Pour ce faire, un changement de base doit être effectué du repère global vers le repère local de
chacun des éléments.
Quand on se limite à des problèmes plans, les formules de changement de base permettant
de passer du repère global au repère local s’écrivent :
14
À titre d’exemple, le changement de base d’un élément barre plan à deux nœuds:
D’une manière générale, ce changement de base permet de calculer les déplacements
des nœuds de l’élément e en repère local {qe} à partir de ceux du repère global {Qe} via une
matrice de passage nommée [Re].
L’énergie de déformation We devient alors :
D’où l’expression de la rigidité élémentaire en repère global :
15
Soit pour l’élément barre :
I.14 Organigramme général de résolution
16
I.15 Élément de poutre plane à 2 nœuds
Soit un élément de poutre plane de longueur L, de section S et d’inertie I constantes et
dont le matériau a un module d’élasticité longitudinal E.
Cet élément à deux nœuds destiné au calcul des réseaux de poutres chargés dans leur plan fait
appel à la théorie des poutres qui permet de ramener le problème tridimensionnel à un problème
unidimensionnel en condensant ses caractéristiques au niveau de sa fibre moyenne. Chacun de
𝑑𝑣𝑖
ses nœuds possède trois degrés de liberté ui ,vi et qui permettent de reconstituer les fonctions
𝑑𝑥
de déplacements axial u(x) et transversal v(x). Ses caractéristiques sont donc :
L : Longueur.
E : Module d’élasticité longitudinal.
S : Section axiale.
I = Iz : Moment d’inertie suivant l’axe z (perpendiculaire au plan xy).
Les inconnues étant les déplacements en i et en j, la première opération consistera à exprimer
𝑑𝑣𝑖 𝑑𝑣𝑗
ces champs de déplacement en fonction des valeurs nodales ui, vi, et uj, vj, .
𝑑𝑥 𝑑𝑗
17
I.16 Matrices de rigidité élémentaire en repère local
Champ de déplacement axial u(x)
En l’absence de flexion, l’élément du fait de son raccourcissement ou de son allongement,
a un comportement identique à celui d’une barre. u(x) est donc encore une fois une fonction
linéaire de la forme u(x) = a0 + a1x. La déformation associée est bien évidemment homogène
𝑑𝑢
sur la section et vaut comme en 𝜀𝑥𝑥 = = 𝑎1 .
𝑑𝑥
De plus et comme u(0) = ui en i et u(L) = uj en j, on obtient comme pour l’élément barre,
la même fonction de déplacement que celle établie en :
Champ de déplacement transversal v(x)
En flexion, chacun des nœuds de cet élément possède deux degrés de liberté, un en
translation dans la direction transversale v et un associé à la rotation β autour de l’axe z
perpendiculaire au plan xy. Quatre conditions aux limites (vi, βi ,vj ,βj ) peuvent donc être
utilisées pour définir sa fonction d’approximation ce qui explique l’expression de v(x) sous la
forme d’un polynôme de degré 3.
Cette approche revient en fait à intégrer l’équation générale des poutres avec une charge
répartie nulle. L’équation du moment fléchissant sera alors une fonction linéaire.
À partir des conditions aux limites de v(x) (i.e. v(0) = vi ,v_(0) = βi ,v(L) = v j ,v_(L)= β j ) , le
vecteur des déplacements nodaux s’écrit :
18
D’où finalement,
C’est d’ailleurs grâce à cette relation que certains logiciels arrivent à tracer la déformée
entre les nœuds. Ils calculent les valeurs de v(x) en certains points intermédiaires et
reconstituent ensuite la déformée. Ceci permet de se limiter à la géométrie de la structure et
donc de réduire le nombre de nœuds.
Matrice de rigidité élémentaire
Reprenant l’expression de l’énergie de déformation élémentaire,
et sachant que la matrice [H], les vecteurs contrainte et déformation se réduisent respectivement
aux seuls termes E, σxx et εxx, We s’écrit :
19
D’où l’expression finale de la matrice de rigidité élémentaire :
Finalement :
12 6𝑙 −12 6𝑙
𝐸𝐼 6𝑙 4𝑙 2 −6𝑙 2𝑙 2 ]
[𝐾𝑒 ] = 3 [
𝑙 −12 −6𝑙 12 −6𝑙
6𝑙 2𝑙 2 −6𝑙 4𝑙 2
20
II.1 Elément fini bi dimensionnel (2D)
On utilise avec une grande facilite la méthode des éléments finis pour résoudre des
problèmes à deux dimensions pour différents matériaux et même pour des géométries et des
contours complexes. Pour les problèmes à deux dimensions on utilise des éléments triangulaires
ou rectangulaires.
Figure II.1 : Exemples d’élément à deux dimensions
II.2 Grille à deux dimensions
L'élément triangulaire linéaire est composé de cotés droits et possède un nœud à chaque
coin. Ce type d'élément est généralement recommande pour les contours irréguliers. La fonction
d’interpolation d’un élément triangulaire est de la forme :
φ = 𝑎1 + 𝑎2 𝑥 + 𝑎3 𝑦 (II.1)
L’élément rectangulaire linéaire est composé de cotés droits et possède un nœud à
chaque coin. On utilise l’élément rectangulaire dans le cas contours réguliers. Sa fonction
d’interpolation est de la forme :
φ = 𝑐1 + 𝑐2 𝑥 + 𝑐3 𝑥𝑦 + 𝑐4 𝑦 (II.2)
II.3 Fonction d’interpolation d’un élément triangulaire à deux dimensions
La fonction d’interpolation d’un élément à deux dimensions triangle est de la forme :
φ = 𝑎1 + 𝑎2 𝑥 + 𝑎3 𝑦
𝜑𝑖 , 𝜑𝑗 𝑒𝑡 𝜑𝑘 représentent les valeurs nodales de la fonction (𝜑) aux nœuds (i), (j) et (k).
Dans le plan (x,y) les coordonnées des nœuds (i), (j) et (k) sont respectivement :
𝑖(𝑥𝑖 , 𝑦𝑖 ), 𝑗(𝑥𝑗 , 𝑦𝑗 ) et 𝑘(𝑥𝑘 , 𝑦𝑘 )
21
Figure II.2 : Fonction d’interpolation d’un élément triangulaire à deux dimensions
Pour 𝑥 = 𝑥𝑖 𝑒𝑡 𝑦 = 𝑦𝑖 𝜑 = 𝜑𝑖
Pour 𝑥 = 𝑥𝑗 𝑒𝑡 𝑦 = 𝑦𝑗 𝜑 = 𝜑𝑗
Pour 𝑥 = 𝑥𝑘 𝑒𝑡 𝑦 = 𝑦𝑘 𝜑 = 𝜑𝑘
Si on remplace dans l’équation (II.1) on obtient le système d’équation suivant :
𝜑𝑖 = 𝑎1 + 𝑎2 𝑥𝑖 + 𝑎3 𝑦𝑖
{ 𝑗 = 𝑎1 + 𝑎2 𝑥𝑗 + 𝑎3 𝑦𝑗
𝜑 (II.3)
𝜑𝑘 = 𝑎1 + 𝑎2 𝑥𝑘 + 𝑎3 𝑦𝑘
La résolution de ce système permet de déterminer les coefficients de la fonction d’interpolation :
1
𝑎1 = 2𝐴 [(𝑥𝑗 𝑦𝑘 − 𝑥𝑘 𝑦𝑗 )𝜑𝑖 + (𝑥𝑘 𝑦𝑖 − 𝑥𝑖 𝑦𝑘 )𝜑𝑗 + (𝑥𝑖 𝑦𝑗 − 𝑥𝑗 𝑦𝑖 )𝜑𝑘 ] (II.4)
1
𝑎2 = 2𝐴 [(𝑦𝑗 − 𝑦𝑘 )𝜑𝑖 + (𝑦𝑘 − 𝑦𝑖 )𝜑𝑗 + (𝑦𝑖 − 𝑦𝑗 )𝜑𝑘 ] (II.5)
1
𝑎3 = 2𝐴 [(𝑥𝑘 − 𝑥𝑗 )𝜑𝑖 + (𝑥𝑖 − 𝑥𝑘 )𝜑𝑗 + (𝑥𝑗 − 𝑥𝑖 )𝜑𝑘 ] (II.6)
La quantité (2A) représente l’aire du triangle (de l’élément).
1 𝑥𝑖 𝑦𝑖
2𝐴 = det [1 𝑥𝑗 𝑦𝑗 ] (II.7)
1 𝑥𝑘 𝑦𝑘
La substitution de (II.4), (II.5) et (II.5) dans (II.1) et (II.3) nous conduit à l’expression de la
fonction d’interpolation 𝜑 :
𝜑 = 𝑁𝑖 𝜑𝑖 + 𝑁𝑗 𝜑𝑗 + 𝑁𝑘 𝜑𝑘 (II.8)
Avec :
1
𝑁𝑖 = 2𝐴 [𝑎𝑖 + 𝑏𝑖 𝑥 + 𝑐𝑖 𝑦] (II.9)
1
𝑁𝑗 = 2𝐴 [𝑎𝑗 + 𝑏𝑗 𝑥 + 𝑐𝑗 𝑦] (II.10)
1
𝑁𝑘 = 2𝐴 [𝑎𝑘 + 𝑏𝑘 𝑥 + 𝑐𝑘 𝑦] (II.11)
22
Et :
𝑎𝑖 = 𝑥𝑗 𝑦𝑘 − 𝑥𝑘 𝑦𝑗 𝑎𝑗 = 𝑥𝑘 𝑦𝑖 − 𝑥𝑖 𝑦𝑘 𝑎𝑘 = 𝑥𝑖 𝑦𝑗 − 𝑥𝑗 𝑦𝑖
𝑏𝑖 = 𝑦𝑗 − 𝑦𝑘 𝑏𝑗 = 𝑦𝑘 − 𝑦𝑖 𝑏𝑘 = 𝑦𝑖 − 𝑦𝑗
𝑐𝑖 = 𝑥𝑘 − 𝑥𝑗 𝑐𝑗 = 𝑥𝑖 − 𝑥𝑘 𝑐𝑘 = 𝑥𝑗 − 𝑥𝑖
II.4 Construction des fonctions d’interpolation
La méthode des éléments finis est basé sur le principe de l’utilisation d’une
approximation nodale par sous domaines. Sur chaque sous domaine on approche la fonction
𝑈(𝑥) par :
𝑈(𝑥) = ∑𝑛𝑖=1 𝑁𝑖 (𝑥). 𝑢𝑖 = 𝑁(𝑥). 𝑢𝑛 (II.12)
Selon le type d’élément et de son nombre de nœuds, les fonctions d’interpolation
diffèrent. En général on choisit des fonctions d’interpolation sous forme de séries polynomiales.
Pour le cas des problèmes uni dimensionnels la fonction d’interpolation est généralement de la
forme :
𝑈(𝑥) = ∑𝑛𝑖=0 𝑎𝑖 𝑥 𝑖 (II.13)
Pour les problèmes bidimensionnels elle est généralement de la forme :
𝑈(𝑥) = ∑𝑝𝑖=1 𝑎𝑖 𝑥 𝑗 𝑦 𝑘 𝑗+𝑘 ≤𝑛 (II.14)
Pour les problèmes tri dimensionnels :
𝑈(𝑥) = ∑𝑝𝑖=1 𝑎𝑖 𝑥 𝑗 𝑦 𝑘 𝑧 𝑙 𝑗+𝑘+𝑙 ≤ 𝑛 (II.15)
Dans les relations (II.14) et (II.15), (n) représente le degré d’interpolation.
Pour les cas simples, la détermination des fonctions d’interpolation peut être faite
directement. En général on fait recours aux éléments de références.
II.5 Détermination directe
L’objectif est de déterminer un polynôme qui prend des valeurs connues en un certain
nombre de points (nœuds). Le nombre de termes qui composent ce polynôme doit égal au
nombre de points.
On utilise par exemple une interpolation de Lagrange de type (C0) pour un élément barre
a deux nœuds.
23
Figure II.3 : Elément barre à deux nœuds.
Au nœud (1) nous avons :
𝑢(𝑥1 ) = 𝑢1 (II.16)
Et au nœud (2) :
𝑢(𝑥2 ) = 𝑢2 (II.17)
Sachant que :
𝑈(𝑥) = 𝑁1 (𝑥)𝑢1 + 𝑁2 (𝑥)𝑢2
Avec :
𝑥2 − 𝑥
𝑁1 (𝑥) =
𝑥2 − 𝑥1
Et :
𝑥 − 𝑥1
𝑁2 (𝑥) =
𝑥2 − 𝑥1
Dans le cas uni dimensionnel, les polynômes utilises dans l’écriture des fonctions
d’interpolation sont généralement définis par :
𝑥−𝑥
𝑁𝑖 (𝑥) = ∏𝑛𝑗=1 (𝑥 −𝑥𝑗 ) (II.18)
𝑖 𝑗
𝑖≠𝑗
Cette relation vérifie bien l’équation (33.2) :
1 𝑖=𝑗
𝑁𝑖 (𝑥𝑗 ) = 𝛿𝑖𝑗 = {
0 𝑖≠𝑗
On peut procéder de la même façon pour dans le cas d’un élément bidimensionnel comme par
exemple le cas d’un élément triangle a trois nœuds :
𝑈(𝑥) = 𝑁1 . 𝑢1 + 𝑁2 . 𝑢2 + 𝑁3 . 𝑢3 (II.19)
La présence de trois nœuds implique trois coefficients dans l’écriture de la fonction de Forme
𝑁𝑖 (𝑥) :
𝑁𝑖 (𝑥𝑗 ) = 𝑎𝑖 + 𝑏𝑖 𝑥 + 𝑐𝑖 𝑦 (II.20)
Avec :
1
𝑁1 (𝑥) = 2𝐴 [(𝑦3 − 𝑦2 )(𝑥2 − 𝑥) − (𝑥3 − 𝑥2 )(𝑦2 − 𝑦)] (II.21)
1
𝑁2 (𝑥) = 2𝐴 [(𝑦1 − 𝑦3 )(𝑥3 − 𝑥) − (𝑥1 − 𝑥3 )(𝑦3 − 𝑦)] (II.22)
1
𝑁3 (𝑥) = 2𝐴 [(𝑦2 − 𝑦1 )(𝑥1 − 𝑥) − (𝑥2 − 𝑥1 )(𝑦1 − 𝑦)] (II.23)
2𝐴 = (𝑥3 − 𝑥2 )(𝑦1 − 𝑦2 ) − (𝑥1 − 𝑥2 )(𝑦3 − 𝑦2 ) (II.24)
24
II.6 Elément de référence
L’élément de référence est un élément de forme simple qui peut être transforme en
chaque élément réel par une transformation géométrique () comme montre sur la figure II.4.
Figure II.4 : Elément de référence à deux dimensions
La transformation (𝜏) définit les coordonnées (x) d’un point de l’élément réel en fonction des
coordonnées (r) de l’élément de référence :
𝜏 ∶ 𝑟 → 𝑥 = 𝑥(𝑟) (II.25)
II.6.1 Intérêt des éléments de référence
L’élément de référence permet de remplacer la détermination des fonctions
d’interpolation sur l’élément réel par celles sur l’élément de référence à cause de sa forme très
simple.
II.6.2 Eléments de référence classiques
En fonction de la dimension de l’élément, il existe plusieurs types d’élément de référence
comme par exemple :
Elément de référence à une dimension
25
Elément de référence à deux dimensions triangulaires
Eléments de référence à deux dimensions carrées
Eléments de référence à trois dimensions tétraédriques
Eléments de référence à trois dimensions prismatiques
26
Eléments de référence à trois dimensions cubiques
II.6.3 Construction des fonctions d’interpolation sur l’élément de référence
a/ Définitions
1. Un élément est iso-paramétrique si les nœuds géométriques sont confondus avec les nœuds
d’interpolation.
2. Un élément est sub-paramétrique si le nombre de nœuds géométriques est inférieur au
nombre de nœuds d’interpolation.
3. Un élément est super-paramétrique si le nombre de nœuds géométriques est supérieur au
nombre de nœuds d’interpolation.
b/ Méthode générale de construction des fonctions d’interpolation
La fonction 𝑢(𝑟) est définie sur l’élément de référence par :
𝑢(𝑟) = ∑𝑛𝑖=1 𝑁𝑖 𝑢𝑖 = 𝑁(𝑟)𝑢𝑛 (II.26)
Les fonctions 𝑁𝑖 (𝑟) sont composées de sommes de termes polynomiaux.
𝑢(𝑟) = 𝜑(𝑟)𝑎 = ∑𝑝𝑗=1 𝜑(𝑟)𝑎𝑗 (II.27)
Dans la relation (II.27) φ(r) est la base polynomiale de l’interpolation.
Pour un problème unidimensionnel (1 D), la fonction 𝑢(𝑟) est choisie comme suit :
𝑢(𝑟) = 𝑎0 + 𝑎1 𝑟 + 𝑎2 𝑟 2 + ⋯ (II.28)
Pour un problème bidimensionnel (2 D), la fonction 𝑢(𝑟) est choisie comme suit :
𝑢(𝑟) = 𝑎0 + 𝑎1 𝑟 + 𝑎2 𝑠 + 𝑎3 𝑟𝑠 + 𝑎4 𝑟 2 … (II.29)
27
Le choix du nombre de terme (p) du polynôme d’interpolation dépend de plusieurs paramètres
tels que le type de l’élément et son nombre de nœuds. Pour un triangle à 6 nœuds, le
déplacement u dépend de 𝑢1 , 𝑢2 , 𝑢3 , 𝑢4 , 𝑢5 et 𝑢6 , soit 6 termes.
Pour un élément de référence triangle a 6 nœuds, on choisi comme polynôme
d’interpolation la fonction :
𝑢(𝑟, 𝑠) = 𝑎0 + 𝑎1 𝑟 + 𝑎2 𝑠 + 𝑎3 𝑟𝑠 + 𝑎4 𝑟 2 + 𝑎5 𝑠 2
Figure II.5 : Eléments de référence triangulaire quadratique à 6 nœuds
Pour un élément de référence quadrangle à 8 nœuds, on choisi comme polynôme
d’interpolation la fonction :
𝑢(𝑟, 𝑠) = 𝑎0 + 𝑎1 𝑟 + 𝑎2 𝑠 + 𝑎3 𝑟𝑠 + 𝑎4 𝑟 2 + 𝑎5 𝑠 2 + 𝑎6 𝑠 2 𝑟 + 𝑎7 𝑟 2 𝑠
Construction de N :
Sachant que :
𝑝
𝑢(𝑟) = 𝜑(𝑟)𝑎 = ∑ 𝜑𝑗 (𝑟)𝑎𝑗
𝑗=1
On peut écrire que :
𝑢(𝑟) = 𝜑1 (𝑟)𝑎1 + 𝜑2 (𝑟)𝑎2 + ⋯ + 𝜑𝑝 (𝑟)𝑎𝑝 (II.30)
P : nombre de terme de la base polynomiale.
En écriture matricielle (II.30) devient :
28
𝑎1
𝑎2
𝑢(𝑟) = [𝜑1 (𝑟) 𝜑2 (𝑟) … 𝜑𝑝 (𝑟)] { ⋮ } = 𝜑(𝑟)𝑎 (II.31)
𝑎𝑝
En chaque nœud la condition suivante est vérifiée :
𝑢(𝑟𝑖 ) = 𝑢𝑖 (II.32)
Alors :
𝑎1
𝑎2
𝑢𝑖 = [𝜑1 (𝑟𝑖 ) 𝜑2 (𝑟𝑖 ) … 𝜑𝑝 (𝑟𝑖 )] { ⋮ } (II.33)
𝑎𝑝
Soit sous forme développée :
𝑢1 𝜑1 (𝑟1 ) 𝜑2 (𝑟1 ) … 𝜑𝑝 (𝑟1 ) 𝑎1
𝑢2 ⋮ ⋮ 𝑎2
{ }=[ ] { } (II.34)
⋮ ⋮
𝑢𝑛 𝜑1 (𝑟𝑛 ) 𝜑2 (𝑟𝑛 ) … 𝜑𝑝 (𝑟𝑛 ) 𝑎𝑝
On note :
𝜑1 (𝑟1 ) 𝜑2 (𝑟1 ) … 𝜑𝑝 (𝑟1 )
𝜑𝑛 = [ ⋮ ⋮ ]
⋮ ⋮
𝜑1 (𝑟𝑛 ) 𝜑2 (𝑟𝑛 ) … 𝜑𝑝 (𝑟𝑛 )
Alors :
𝑢𝑛 = 𝜑𝑛 𝑎 (II.35)
𝜑 représente la matrice de valeurs de la base polynomiale aux nœuds.
L’inversion du système (II.34) donne :
𝑎 = 𝜑𝑛−1 𝑢𝑛 (II.36)
Si on remplace (II.36) dans (II.27) on obtient :
𝑢(𝑟) = 𝜑(𝑟)𝑎 = 𝜑(𝑟)𝜑𝑛−1 𝑢𝑛 (II.37)
On remplace (II.37) dans (II.26) pour obtenir :
𝑢(𝑟) = 𝜑(𝑟)𝜑𝑛−1 𝑢𝑛 = 𝑁(𝑟)𝑢𝑛 (II.38)
Finalement, on trouve l’expression de la fonction d’interpolation :
𝑁(𝑟) = 𝜑(𝑟)𝜑𝑛−1 (II.39)
c/ Etapes de la construction de N.
1. Choix du polynôme (le nombre de terme dépend selon le cas du type de l’élément et de son
nombre de nœuds).
2. Détermination de la matrice 𝜑𝑛 des valeurs de la base polynomiale aux nœuds.
29
3. Inversion de 𝜑𝑛 pour l’obtention de 𝜑𝑛−1
4. Calcul de 𝑁(𝑟) = 𝜑(𝑟)𝜑𝑛−1
d/ Application
On se propose de déterminer les fonctions d’interpolation 𝑁(𝑟) pour le cas d’un élément
quadrilatère iso-paramétrique à 4 nœuds.
1. Choix de la base polynomiale : l’élément comporte 4 nœuds, on choisit donc une base
de la forme :
𝜑 = [1 𝑟 𝑠 𝑟𝑠]
2. Détermination de la matrice ϕ4 des valeurs de la base polynomiale aux nœuds :
𝜑1 (𝑟1 ) 𝜑2 (𝑟1 ) … 𝜑𝑝 (𝑟1 )
3. 𝜑𝑛 = [ ⋮ ⋮ ]
⋮ ⋮
𝜑1 (𝑟𝑛 ) 𝜑2 (𝑟𝑛 ) … 𝜑𝑝 (𝑟𝑛 )
Au nœud 1 : (r,s) = (−1,−1)
𝜑 = [1 − 1 − 1 1]
Au nœud 2 : (r,s) = (1,−1)
𝜑 = [1 1 − 1 − 1]
Au nœud 3 : (r,s) = (1,1)
𝜑 = [1 1 1 1]
Au nœud 4 : (r,s) = (−1,1)
𝜑 = [1 − 1 1 − 1]
Alors :
1 −1 −1 1
𝜑𝑛 = [1 1 −1 −1]
1 1 1 1
1 −1 1 −1
4. Inversion de 𝜑𝑛 :
30
1 1 1 1
1 −1 1 1 −1]
𝜑𝑛−1 = [
4 −1 −1 1 1
1 −1 1 −1
4. Calcul de 𝑁(𝑟) = 𝜑(𝑟)𝜑𝑛−1
1 1 1 1
1
𝑁(𝑟) = 𝜑(𝑟)𝜑𝑛−1 = [1 𝑟 𝑠 𝑟𝑠] [ −1 1 1 −1]
4 −1 −1 1 1
1 −1 1 −1
Soit :
1
𝑁1 = (1 − 𝑟)(1 − 𝑠)
4
1
𝑁2 = (1 + 𝑟)(1 − 𝑠)
4
1
𝑁3 = (1 + 𝑟)(1 + 𝑠)
4
1
𝑁4 = (1 − 𝑟)(1 + 𝑠)
4
31