Université de Relizane

Faculté des Sciences et Technologie

Département de Physique

Polycopié
1èr année Master (LMD)

METHODE DES ELEMENTS FINIS

COURS ET APPLICATIONS
Préparé par
Dr. FARES Redouane (M.C.A)

Dans ce cours de la méthode des éléments finis, on présente les principes de base
de cette méthode. Ce cours est destiné aux étudiants en Master 1.

-2023-
Contenu
Chapitre 1 (Introduction) 1
1.1 Qu'est-ce que la méthode des éléments finis 1
1.2 Description de la MEF 1
1.3 Formulation des équations d’élément fini 2
1.3.1 Méthode Galerkin 2
1.3.2 Formulation variationnelle 5
Chapitre 2 (Équations Élément Fini pour Transfert de Chaleur) 8
2.1. Déclaration du problème 9
2.2. Discrétisation des équations de transfert de Chaleur par élément Fini 10
2.3. Différents types de problèmes 11
Chapitre 3 (MEF pour les Problèmes de Mécanique du Solide ) 12
3.1. Déclaration du problème 12
3.2. Les équations d’élément Fini 13
3.3. Assemblage de l’équation du système global 15
Chapitre 4 (Fini Éléments) 18
4.1. Elément triangulaire Bidimensionnel 18
4.2. Eléments Bidimensionnel iso-paramétrique 19
4.2.1. Fonctions de forme 19
4.2.2. Matrice déformation-déplacement 20
4.2.3. Propriétés des éléments 22
4.2.4. Intégration dans des éléments quadrilatéraux 22
4.2.5. Calcul des déformations et des contraintes 24
4.3. Éléments isoparamétriques tridimensionnels 26
4.3.1 Fonctions de formes 26
4.3.2 Matrice déformation-déplacement 27
4.3.3 Propriétés des éléments 28
4.3.4 Calcul efficace de la matrice de rigidité 29
4.3.5 Intégration de la matrice de rigidité 30
4.3.6 Calcul des déformations et des contraintes 30
4.3.7 Extrapolation des déformations et contraintes 30
Chapitre 5 (Discrétisation) 32
5.1 Modèle discret du problème 32
5.2 Génération de maillage 33
5.2.1 Générateurs de maillage 33
5.2.2 Technique de cartographie (Mapping technique) 33
5.2.3 Triangulation de Delaunay 33
Chapitre 6 (Assemblage) 35
6.1 Assemblage et solution 36
6.2 Algorithme de démontage 36
6.3 Assemblage 36
6.3.1 Algorithme d'assemblage pour les vecteurs 37
6.3.2 Algorithme d'assemblage pour les matrices 37
6.4 Conditions aux limites de déplacement 37
6.4.1 Spécification explicite du déplacement BC 38
6.4.2 Méthode des grands nombres 38
6.5 Solution des équations aux éléments finis 39
6.5.1 Méthodes de résolution 40
6.5.2 Méthode LDU directe avec matrice de profil 40
Chapitre 1

Introduction
1.1 Qu'est-ce que la méthode des éléments finis
La méthode des éléments finis (FEM) est une technique numérique pour résoudre des
problèmes qui sont décrits par des équations aux dérivées partielles ou qui peuvent être
formulés comme une minimisation fonctionnelle. Un domaine d'intérêt est représenté comme
un assemblage d'éléments finis. Les fonctions d'approximation en éléments finis sont
déterminées en termes de valeurs nodales d'un champ physique recherché. Un problème
physique continu est transformé en un problème d'éléments finis discrétisé avec des valeurs
nodales inconnues. Pour un problème linéaire, un système d'équations algébriques linéaires doit
être résolu. Les valeurs à l'intérieur des éléments finis peuvent être récupérées à l'aide de
valeurs nodales.
Deux caractéristiques du FEM méritent d'être mentionnées :
1) L'approximation par morceaux des champs physiques sur les éléments finis fournit une
bonne précision même avec des fonctions d'approximation simples (en augmentant le nombre
d'éléments, nous pouvons atteindre n'importe quelle précision).
2) La localité d'approximation conduit à des systèmes d'équations creux pour un problème
discrétisé. Cela aide à résoudre des problèmes avec un très grand nombre d'inconnues nodales.

1.2 Description de la MEF

Pour résumer le fonctionnement de la méthode des éléments finis, nous énumérons les
principales étapes de la méthode des éléments finis.
1. Discrétiser le continuum. La première étape consiste à diviser une région de solution en
éléments finis. L’élément fini engrener est typiquement généré par un programme
préprocesseur. La description consiste de plusieurs tableaux principale (les coordonnées
nodales et les connectivités des éléments).
2. Sélection des fonctions d'interpolation. Les fonctions d'interpolation sont utilisées pour
interpoler les variables sur l'élément. Souvent, les polynômes sont sélectionnés comme
fonctions d'interpolation. Le degré du polynôme dépend aux nombres de nœuds attribué à
l’élément.
3. Recherche des propriétés de l'élément. L'équation matricielle pour l'élément fini doit être
établie qui relie les valeurs nodales de la fonction inconnue à d'autres paramètres. Pour cette
tâche différente approches peut être utilisé; le plus pratique sommes: l’approche variationnel et
la méthode de Galerkin.
4. Assemblage des équations des éléments. Pour trouver le système d'équation global pour
toute la solution, nous devons assembler toutes les équations élémentaires. En d'autres termes,
nous devons combiner les équations locales d'élément pour tous les éléments utilisés pour la
discrétisation. Les connectivités des éléments sont utilisées pour l'assemblage traité. Avant
l’exécution des calculs, les conditions aux limites (qui ne sont pas prises en compte dans les
équations des éléments) doivent être imposées.

1
 Fig 1.1 : Deux éléments linéaires unidimensionnels et une fonction
d’interpolation à l'intérieur de l’élément.

5. Résoudre le l’équation du système global. L’équation global du système est typiquement

symétrique et définie positive. Des méthodes directes et itératives peuvent être utilisées pour
la solution.
6. Calcule des résultats supplémentaires. Dans de nombreux cas, nous pouvons calculer des
paramètres supplémentaires.

1.3 Formulation des équations d’élément fini

Plusieurs approches peuvent être utilisées pour transformer la formulation physique du
problème en son analogue discret par élément fini. Si la formulation physique du problème est
connue sous le nom d'équation différentielle, alors la méthode la plus populaire de sa
formulation par éléments finis est la méthode Galerkin. Si le problème physique peut être
formulé comme la minimisation d'une formulation fonctionnelle puis variationnelle des
équations aux éléments finis est généralement utilisée.

Méthode Galerkin

Nous allons utiliser l’exemple unidimensionnel pour l’explication de la formulation élément

fini en utilisant la méthode Galerkin. Supposant que nous avons besoin de résoudre
numériquement l’équation différentiel suivante:

Avec les conditions aux limites suivantes

Où u est une solution inconnue. Nous allons résoudre le problème en utilisant deux éléments
fini unidimensionnels, comme indiquée dans Figure. 1.1.
D'abord, considérons un élément fini présenté à droite de la figure. L'élément à deux nœuds et
l'approximation de la fonction u(x) peut se faire comme suit :

2
Où Ni sont appelées fonctions de formes

Lequel sont utilisées pour l’interpolation de u (x) en utilisant ses valeurs nodal. Le s valeurs
Nodales u1 et u2 sont inconnues, ou ils devraient être déterminés par la discrétisation global de
l’équation du système.
Après avoir remplacé u exprimé par ses valeurs nodales et les fonctions de formes, dans
l’équation différentielle, nous aurons la formulation approximative suivante:

Où ψ est un résidu non nul en raison de la représentation approximative d'une fonction à

l'intérieur d'un élément fini. La méthode Galerkin fournit un résiduel minimiser en multipliant
les termes de l’équation (1.5) par les fonctions de forme, en intégrant dans l’élément:

Utilisation de l'intégration par partie permet d’avoir une formulation discrète de l’équation
différentielle pour l’élément fini:

D'habitude tel relation pour un élément fini est présentée comme suite:

Dans la mécanique des solide [k] est appelé matrice de raideur et F est appelé vecteur de
charge. Dans le Cas de deux éléments fini de longueur L, les matrices de raideur et de vecteurs
de charges peut être calculé comme suite:

3
 Les relations ci-dessus fournissent des équations d'éléments finis pour les deux éléments
finis séparés. L’équation globale du système pour le domaine avec 2 éléments et 3 nœuds peut
être obtenue par un Assemblée des équations des éléments. Dans notre cas simple, il est clair
que les éléments interagissent les uns avec les autres au nœud avec le nœud numéro 2.
L’assemblage global des équations du système est:

Fig 1.2 : Comparaison de la solution élément fini et la solution exact.

Après application de la condition initial u( x=0)=0 l’apparence final de l’équation global du

système est:

Les valeurs Nodal ui s o n t obtenues comme sommes des résultats de la solution de système
d’équation algébrique linéaire. L’évaluation de u en tout point à l'intérieur d'un élément fini
peut être calculée à l'aide des fonctions de forme. La solution élément fini de l’équation
différentiel est montrée dans Figure. 1.2 pour a = 1,b = 1 ,L = 1 et R = 1 .

La solution exacte est une fonction quadratique. La solution par éléments finis avec
l'utilisation de l'élément le plus simple est par morceaux linéaire. Pour une solution élément
fini précis peut être obtenu en augmentant le nombre des éléments ou alors avec l’utilisation
d’éléments avec suite compliqué des fonctions de formes. Notant que la méthode des
éléments finis fournit des valeurs exactes de u (juste pour ce problème particulier). La MEF
avec des fonctions de forme linéaires produisent des valeurs nodales exactes si la solution
recherchée est quadratique.

4
Formulation variationnelle

L’équation différentielle

avec a = EA a la signification physique suivante en mécanique des solides. Il décrit la

tension d’un bar unidimensionnel avec une surface A transversale de Matériel ou le module
d’élasticité E et soumis à une charger b distribué et une charge concentré R à son extrémité
droite comme montré dans Figure 1.3.
Le problème peut être formulé en termes suivantes pour minimiser le potentiel d’énergie
fonctionnelle Π:

Fig 1.3 : Tension d’une bar unidimensionnel soumis à une

charge distribué et une charge concentré.

En utilisant représentation de { u} avec les fonctions de formes (1.3)-(1.4) nous pouvons

évaluer l’énergiepotentielle pour le deuxième élément fini comme:

La condition pour minimiser Π est:

5
Ainsi que (1.15) est équivalent à

Il est facile à vérifier que la différenciation de Π avec ui respecté donne l’élément fini et
coïncide avec l’équation obtenu par la méthode de Galerkin:

Exemple. Donner les fonctions de formes pour l’élément quadratique unidimensionnel

avec Trois nœuds. Utiliser le système de coordonnées local −1 ≤ ξ ≤ 1 .

La solution. Avec les fonctions de formes, tout domaine à l'intérieur de élément est
présenté comme suite:

Aux nœuds l’approximation d’une fonction devraient être égal à sa valeur nodal:

Depuis l’élément a Trois nœuds les fonctions de formes peut être un polynômes quadratique
(avec Trois coefficient). La fonction de forme N 1 peut être écrit comme:

N 1 = α 1 + α 2ξ + α 3ξ 2

Les coefficients inconnus α i sont définis à partir du système d'équations suivant :

6
Les solutions sont: α1=0, α 2 = −1 / 2 , α 3 = 1 / 2 . Ainsi le façonner une fonction N 1 est
égal à:

De la même manière il est possible d’obtenir les fonctions de formes N2 et N3 :

7
Chapitre 2
Équations Élément Fini pour Transfert de
Chaleur
2.1. Déclaration du problème
Considérant un corps isotrope avec une température dépendante. L’équation de base de
transfert de chaleur est donnée comme suite:

Ou qx , qy et qz sont des composants du flux de chaleur à travers la surface unitaire ;

Q=Q(x, y, z, t ) le taux de chaleur générer par unité de volume ; ρ est la densité du matériau ; c
est la capacité calorifique ; T est la température et t est le temps. Selon la loi de Fourier les
Composants de flux de Chaleurs peut être exprimé comme suit :

Où k est le coefficient conductivité thermique. Substitution des relations de Fourier donnent

les équations de transfert de chaleur:

Les conditions de frontière peut être assumé comme suite:

1. Température Spécifié

2. Flux de chaleur Spécifié

8
3. les conditions aux limites de la Convection

Où h est le coefficient de convection, Ts est une température de surface inconnue et Te est

Température environnementale inconnue.

4. Radiation

où σ est la constante de Stefan-Boltzmann ; ε est le coefficient d'émission surfacique ; α est la

coefficient d’absorption surfacique absorbée et qr est le flux de chaleur entrant par unité
surface.
Pour les problèmes transitoires, il est nécessaire à spécifier une Température au temps t=0 :

2.2. Discrétisation des équations de transfert de Chaleur par élément Fini

Un domaine V est divisé en éléments finis connectés aux nœuds. Nous allons écrire toutes les
relations pour un élément fini. Les équations globales pour le domaine peuvent être assemblées
à partir d'équations d'éléments finis en utilisant la connectivité.
Les fonctions de formes Ni sont utilisés pour l'interpolation de température à l'intérieur un fini
élément:

La différentiel des équations d’interpolation de la Température donnent les relations

l’interpolation Suivante pour la température dégradés :

Ici T est une vecteur de températures aux nœuds ; [ N ] est une matrice des fonctions de formes
et [ B ] est une matrice pour les gradient d’interpolation de température.
En utilisant la méthode de Galerkin, nous pouvons récrire l’équation basique de transfert de
chaleur sous la forme suivante:

Appliquant le théorème de divergence au trois premier termes, nous arrivons aux relations
suivantes:

9
où { n} est la normal sur la surface du corps. Après insertion des conditions initiales et aux
limites dans l’équation (2.8), les équations discrétisés sont comme suites :

Nous pouvons écrire comme ça

Les équations des éléments finis discrétisées pour les problèmes de transfert de chaleur ont les
formulations finales suivantes:

10
2.3. Différents types de problèmes

Des équations pour différents types de problèmes peuvent être déduites de l'équation générale
ci-dessus :
Problème linéaire stationnaire

Problème nonlinéaire stationnaire

Problème linéaire transitoire

Problème non linéaire transitoire

11
Chapitre 3

MEF pour les Problèmes de

Mécanique du Solide
3.1. Déclaration du problème
Considérons un corps élastique tridimensionnel soumis à des forces de surface et de corps et à
la température domaine. De plus, les déplacements sont spécifiés sur certaines surfaces. Pour
une géométrie donnée du corps, charges appliquées, conditions aux limites de déplacement,
champ de température et loi contrainte-déformation du matériau, il est nécessaire pour
déterminer le champ de déplacement du corps. Les déformations et contraintes
correspondantes sont aussi intéressées.
Le déplacements sur coordonner axes x , y et z sommes défini par le vecteur déplacement {u} :

{u} = {u v w} (3.1)

Six composants de déformation différents peuvent être placés dans le vecteur de déformation
{ε} :

{ ε } = {εx εy εz γxy γyz γzx} (3.2)

Pour les petites déformations, la relation entre les déformations et les déplacements est :

{ ε } = [ D ]{ u } (3.3)

où [ D ] est l’opérateur de la matrice différentielle:

Six composants de stress différents forment le vecteur de stress {σ} :

Qui sont liées aux déformations des corps élastiques par la loi de Hook :

12
Ici {εe} est la partie élastique des déformations ; {εt} est la partie thermique des déformations ;
α est le coefficient d’expansion thermique ; T est la Température. La matrice d’ élasticité [ E ] a
l’apparence Suivante:

où λ et µ sont des constantes élastiques de Lame qui peuvent être exprimées par le module
d'élasticité E et Coefficient de Poisson ν :

Le but de la solution par éléments finis du problème élastique est de trouver un tel champ de
déplacement qui fournit un minimum à la fonctionnelle de l'énergie potentielle totale Π :

Ici { } { } est le vecteur de la force corporelle et{ } { } le vecteur de

force superficielle. Les déplacements prescrits sont spécifiés sur la partie de la surface du corps
où les forces de surface sont absentes.
Les conditions aux limites de déplacement ne sont pas présentes dans la fonctionnelle de Π.
Pour cette raison, les conditions aux limites de déplacement doivent être mises en œuvre après
l'assemblage des équations aux éléments finis.
3.2. Les équations d’élément Fini
Considérons un élément fini tridimensionnel abstrait ayant le vecteur des déplacements nodaux {q}:

Les déplacements en un point à l'intérieur d'un élément fini {u} peuvent être déterminés à l'aide
des déplacements nodaux {q} et des fonctions de forme Ni :

13
Ces relations peuvent être réécrites sous forme matricielle comme suit :

Les déformations peuvent également être déterminées par des déplacements aux points
nodaux :

La matrice [B] est appelée matrice des déplacements différentielles. Il peut être obtenu par les
déplacements différentiels exprimés par des fonctions de forme et des déplacements nodaux :

Maintenant, en utilisant les relations pour les contraintes et les déformations, nous sommes en
mesure d'exprimer l'énergie potentielle totale à travers les déplacements nodaux :

Les déplacements nodaux {q} qui correspondent au minimum de la fonctionnelle Π sont

déterminés par les conditions:

La différenciation de Π par rapport aux déplacements nodaux {q} produit les équations
d'équilibre suivantes pour un élément fini :

14
qui se présente généralement sous la forme suivante :

Ici [k] est la matrice de rigidité de l'élément ; {f} est le vecteur de charge ; {p} est le vecteur
des forces réelles et {h} est le vecteur thermique qui représente les forces fictives pour
modéliser la dilatation thermique.

3.3. Assemblage de l’équation du système global

Le but de l'assemblage est de former le système d'équation global

Utilisant les équations d'éléments

Ici [ki], [qi] et [fi] sont la matrice de rigidité, le vecteur de déplacement et le vecteur de charge du
ième élément fini; [K], {Q} et {F} sont la matrice de rigidité globale, le vecteur de déplacement
et le vecteur de charge.
Afin de dériver un algorithme d'assemblage, présentons l'énergie potentielle totale pour le corps
comme une somme des énergies potentielles des éléments πi :

où E 0 esti le fraction de énergie potentielle en relation à libre dilatation thermique:

15
Introduisons les vecteurs suivants et une matrice où les vecteurs d'éléments et les matrices
sont simplement placer:

Il est évident qu'il est facile de trouver la matrice [A] telle que

L'énergie potentielle totale du corps peut être réécrite sous la forme suivante :

En utilisant la condition du minimum de l'énergie potentielle totale

on arrive au système d'équation globale suivant :

La dernière équation montre que les algorithmes d'assemblage de la matrice de rigidité globale
et le vecteur de charge globale:

Ici [A] est la matrice assurant la transformation de l'énumération globale à l'énumération locale.
Fraction non nulle (unité) entrée dans la matrice [A] est très faible. De ce fait, la matrice [A]
n'est jamais utilisée explicitement dans codes informatiques réels..

Exemple. Écrivez une matrice [A], qui relie les numéros de nœuds locaux (élément) et
globaux (domaine) pour le maillage éléments finis suivant :

Ordre des nœuds

pour un élément

16
La solution.

Pour rendre la représentation matricielle compacte, supposons que chaque nœud a un degré de
liberté (notez que dans un problème de mécanique des solides en trois dimensions, il y a trois
degrés de liberté à chaque nœud). La matrice [A] relie les valeurs élémentaires et nodales
globales de la manière suivante :

où {Q} est un vecteur global de valeurs nodales et {Qd} est un vecteur contenant des vecteurs
d'éléments. Les la réécriture explicite de la relation ci-dessus se présente comme suit :

17
Chapitre 4
Fini Éléments
4.1. Élément triangulaire Bidimensionnel

L'élément fini triangulaire a été le premier élément fini proposé pour les problèmes continus. À
cause de simplicité, il peut être utilisé comme introduction à d'autres éléments. Un élément fini
triangulaire dans le système de coordonnées xy est représenté sur la figure 4.1. Puisque
l'élément a trois nœuds, l'approximation linéaire de déplacements u et v est choisi :

Les fonctions de forme Ni(x, y) peuvent être déterminées à partir du système d'équations
suivant :

Les fonctions de forme pour l'élément triangulaire peuvent être présentées comme suit :

Figure 4.1 : L'élément fini triangulaire est l'élément bidimensionnel le

plus simple.

où ∆ est la surface de l'élément. La matrice [B] d'interpolation des déformations à l'aide des
déplacements nodaux est égale a :

18
La matrice d'élasticité [E] a l'allure suivante pour les problèmes plans :

où λ et µ sont des constantes de Lame :

Ici E est le module d'élasticité et ν est le coefficient de Poisson.

La matrice de rigidité pour l'élément triangulaire à trois nœuds peut être calculée comme suit :

Ici, il a été pris en compte que les deux matrices [B] et [E] ne dépendent pas des coordonnées.
C'était supposé que l'élément a une épaisseur unitaire. Comme la matrice [B] est constante à
l'intérieur de l'élément, les déformations et les contraintes sont également constantes à
l'intérieur de l'élément triangulaire.
4.2. Eléments Bidimensionnel iso-paramétrique
Les éléments finis iso-paramétriques sont basés sur la définition paramétrique des coordonnées et
du déplacement les fonctions. Les mêmes fonctions de forme sont utilisées pour la spécification de
la forme de l'élément et pour l'interpolation de le déplacement domaine.

4.2.1. Fonctions de forme

Les éléments finis iso-paramétriques bidimensionnels linéaires et quadratiques sont présentés à

la Figure 4.2. Les fonctions de forme Ni sont définies en coordonnées locales ξ, η (−1 ≤ ξ, η ≤
1). Les mêmes fonctions de forme sont utilisées pour les interpolations de déplacements et de
coordonnées :

où u, v sont des composantes de déplacement au point de coordonnées locales (ξ, η) ; ui, vi sont
les valeurs de déplacement aux nœuds de l'élément fini ; x, y sont les coordonnées des points et
xi, yi sont les coordonnées des nœuds des éléments. La forme matricielle des relations ci-dessus
est la suivante :

19
 Linéaire élément Quadratique élément

Figure 4.2 : Eléments finis linéaires et quadratiques et leur

représentation dans le repère local.

où la matrice d'interpolation pour les valeurs nodales est :

Les fonctions de forme pour les éléments isoparamétriques bidimensionnels linéaires et

quadratiques sont données par :
Elément linéaire (4 nœuds) :

Élément quadratique (8 nœuds):

Dans les équations ci-dessus, la notation suivante est utilisée : ξ0 = ξξi , η0 = ηηi où ξi, ηi
sont les valeurs de coordonnées locales ξ, η aux nœuds.

20
4.2.2. Matrice déformation-déplacement

Pour le problème plan, le vecteur de déformation contient trois composantes :

La matrice déformation-déplacement utilisée pour calculer les déformations en tout point à

l'intérieur de l'élément à l'aide des déplacements nodaux est :

Alors que les fonctions de forme sont exprimées par les coordonnées locales ξ, η, la matrice
déformation-déplacement contient des dérivées par rapport aux coordonnées globales x, y. Les
dérivées peuvent être facilement converties d'un système de coordonnées à l'autre au moyen
de la règle en chaîne de différenciation partielle :

où [J] est la matrice jacobienne. Les dérivées par rapport aux coordonnées globales sont
calculées à l'aide de l'inverse de la matrice jacobienne:

21
Les composantes de la matrice jacobienne sont calculées à l'aide des dérivées des fonctions de
forme Ni par rapport aux coordonnées locales ξ, η et aux coordonnées globales des nœuds
d'élément xi, yi :

Le déterminant de la matrice jacobienne |J| est utilisé pour la transformation des intégrales du
système de coordonnées global vers le système de coordonnées local :

4.2.3. Propriétés des éléments

Les matrices et vecteurs d'éléments sont calculés comme suit :
Matrice de rigidité

Vecteur de force (charges volumiques et surfaciques)

Vecteur thermique (forces fictives pour simuler la dilatation thermique)

La matrice d'élasticité [E] est donnée par la relation (4.5).

4.2.4. Intégration dans des éléments quadrilatéraux

L'intégration des expressions pour les matrices de rigidité et les vecteurs de charge ne
peut pas être effectuée analytiquement pour le cas général des éléments isoparamétriques. Au
lieu de cela, les matrices de rigidité et les vecteurs de charge sont généralement évalués
numériquement à l'aide de la quadrature de Gauss sur les régions quadrilatérales. La formule de
quadrature de Gauss pour l'intégrale de volume dans le cas bidimensionnel est de la forme :

où ξi, ηj sont des abscisses et wi sont des coefficients de pondération de la règle

d'intégration de Gauss. Les abscisses et les poids de la quadrature de Gauss pour n = 1,2,3 sont
donnés dans le tableau suivant:
Abscisses et poids de Gauss quadrature

22
Pour calculer l'équivalent nodal de la charge de surface, l'intégrale de surface est remplacée par
l'intégration linéaire le long du côté de l'élément. La fraction de la charge surfacique est évaluée
comme suit :

Ici s est une coordonnée globale le long du côté élément et ξ est une coordonnée locale le long
du côté élément.
Si la charge de charge répartie est appliquée le long de la normale au côté de l'élément, comme
illustré à la Fig. 4.3, l'équivalent nodal d'une telle charge est :

Exemple. Calculez les équivalents nodaux d'une charge répartie d'intensité constante
appliquée sur le côté d'un élément quadratique bidimensionnel :

Figure 4.3 : Charge normale répartie sur un côté d'un élément quadratique.

La solution. L'équivalent nodal de la charge répartie est calculé comme suit :

23
ou

Les fonctions de forme pour l'élément quadratique unidimensionnel sont :

Les valeurs des forces nodales au nœud 1, 2 et 3 sont définies par intégration :

L'exemple montre que l'approche physique ne peut pas être utilisée pour l'estimation des
équivalents nodaux d'une charge répartie. Cela fonctionne pour les éléments linéaires;
cependant, cela ne fonctionne pas pour les éléments plus compliqués.

4.2.5. Calcul des déformations et des contraintes

Les déformations en tout point d'un élément sont déterminées à l'aide des relations de
Cauchy (4.14) avec l'utilisation de la matrice de différenciation des déplacements (4.15) :

Figure 4.4 : Numérotation des points d'intégration et des sommets pour l'élément
isoparamétrique à 8 nœuds.

Les contraintes sont calculées avec la loi de Hook :

24
où { ε t } est le vecteur de dilatation thermique libre :

La précision des déformations et des contraintes dépend de manière significative de

l'emplacement du point où elles sont calculées. La précision la plus élevée pour les gradients
de déplacement se situe au centre géométrique pour l'élément linéaire et aux points
d'intégration réduits 2 × 2 pour l'élément quadrilatère quadratique.
Pour les éléments quadratiques à 8 nœuds, les déformations et les contraintes ont la
meilleure précision à 2 × 2 points d'intégration avec les coordonnées locales ξ, η = ±1/√3. Une
manière possible de créer un champ de contraintes continu avec une précision raisonnable
consiste à :
1) extrapoler les contraintes des points d'intégration réduits aux nœuds ;
2) faire la moyenne des contributions des éléments finis à tous les nœuds du modèle
d'éléments finis. Stress ultérieurs peut être interpolés à partir de nœuds à l'aide de
fonctions de forme quadratiques.

Considérons l'élément quadratique dans le système de coordonnées local ξ, η comme le

montre la Fig. 4.4 où les points d'intégration sont numérotés comme (1)...(4); les nœuds
d'angle ont les numéros 1...4. Pour l'extrapolation (à l'intérieur d’élément quadratique), nous
allons utiliser des fonctions de forme linéaires. Sous forme matricielle, la relation
d'extrapolation peut être présentée comme suit :

où f(m) sont des valeurs de fonction connues aux points d'intégration réduits ; fi sont les
valeurs des fonctions aux nœuds des sommets et Li(m) est la matrice d'extrapolation
symétrique

25
Élément Linéaire Elément Quadratique

Figure 4.5 : Eléments finis tridimensionnels linéaires et quadratiques

et leur représentation dans le repère local.

4.3 Éléments isoparamétriques tridimensionnels

4.3.1 Fonctions de formes

Les éléments isoparamétriques tridimensionnels linéaires hexaédriques (ou de type brique) à 8

nœuds et quadratiques à 20 nœuds sont représentés à la Fig. 4.5. Le terme "isoparamétrique"
signifie que la géométrie et le champ de déplacement sont spécifiés sous forme paramétrique et
sont interpolés avec les mêmes fonctions. Les fonctions de forme utilisées pour l'interpolation
sont des polynômes des coordonnées locales ξ, η et ζ (−1 ≤ ξ, η, ζ ≤ 1). Les coordonnées et les
déplacements sont interpolés avec les mêmes fonctions de forme :

Ici u, v, w sont les déplacements au point de coordonnées locales ξ, η, ζ ; ui, vi, wi sont des
valeurs de déplacement
aux nœuds ; x, y, z sont les coordonnées des points et xi, yi, zi sont les coordonnées des nœuds.
La matrice de forme
fonctions est :

Les fonctions de forme de l'élément linéaire sont égales à :

26
Pour l'élément quadratique à 20 nœuds les fonctions de forme peuvent s'écrire sous la forme
suivante :

Dans les relations ci-dessus ξi, ηi, ζi sont les valeurs des coordonnées locales ξ, η, ζ aux nœuds.

4.3.2 Matrice déformation-déplacement

Le vecteur de déformation {ε} contient six composantes différentes du tenseur de

déformation :

La matrice déformation-déplacement pour les éléments tridimensionnels a l'aspect suivant :

Les dérivées des fonctions de forme par rapport aux coordonnées globales sont obtenues
comme suit :

où la matrice jacobienne a l'apparence :

27
Les dérivées partielles de x, y, z par rapport à ξ, η, ζ sont trouvées par différenciation des
déplacements exprimés par des fonctions de forme et des valeurs de déplacement nodal :

La transformation des intégrales du système de coordonnées global au système de coordonnées

local est effectuée à l'aide du déterminant de la matrice jacobienne :

4.3.3 Propriétés des éléments

L'équation d'équilibre de l'élément a la forme suivante :

Les matrices et vecteurs d'éléments sont :

Matrice de rigidité

Vecteur de force (charges volumiques et surfaciques)

Vecteur thermique (forces fictives pour simuler la dilatation thermique)

28
La matrice d'élasticité [E] est :

où λ et µ sont des constantes élastiques de Lame qui peuvent être exprimées par le module
d'élasticité E et le coefficient de Poisson ν :

4.3.4 Calcul efficace de la matrice de rigidité

Le calcul de la matrice de rigidité des éléments par multiplication de trois matrices implique de
nombreuses opérations arithmétiques avec des zéros. Après avoir effectué des multiplications sous
forme fermée, les coefficients de la matrice de rigidité des éléments [k] peuvent être exprimés
comme suit :

Ici m, n sont des numéros de nœuds locaux ; i, j sont des indices liés aux axes de
coordonnées (x1, x2, x3). La règle cyclique est utilisée dans l'équation ci-dessus si les
indices de coordonnées deviennent supérieurs à 3
4.3.5 Intégration de la matrice de rigidité
L'intégration de la matrice de rigidité pour les éléments isoparamétriques tridimensionnels
s'effectue dans le repère local ξ, η, ζ :

29
L'application à trois temps de la règle de quadrature de Gauss unidimensionnelle conduit à la
procédure d'intégration numérique suivante :

Habituellement, l'intégration 2 × 2 × 2 est utilisée pour les éléments linéaires et l'intégration 3 ×

3 × 3 est appliquée à l'évaluation de la matrice de rigidité pour les éléments quadratiques. Pour
une intégration plus efficace, il existe une règle spéciale de type Gauss à 14 points, qui fournit
une précision d'intégration suffisante pour les éléments quadratiques tridimensionnels.

4.3.6 Calcul des déformations et des contraintes

Après avoir calculé les matrices d'éléments et les vecteurs, le processus d'assemblage est
utilisé pour composer le système d'équation global. La solution du système d'équations globales
fournit des déplacements aux nœuds du modèle d'éléments finis. En utilisant le désassemblage,
le déplacement nodal pour chaque élément peut être obtenu.
Les déformations à l'intérieur d'un élément sont déterminées à l'aide de la matrice de
différenciation des déplacements :

Les contraintes sont calculées avec la loi de Hook :

où {εt} est le vecteur de dilatation thermique libre :

Il convient de noter que les gradients de déplacement (et donc les déformations et les
contraintes) ont une précision assez différente à différents points à l'intérieur des éléments finis.
La précision la plus élevée pour les gradients de déplacement se situe au centre géométrique
pour l'élément linéaire et aux points d'intégration réduits 2 × 2 × 2 pour l'élément hexagonal
quadratique.

4.3.7 Extrapolation des déformations et contraintes

Pour les éléments quadratiques, les dérivées de déplacement ont la meilleure précision à
2×2×2 points d'intégration avec les coordonnées locales ξ, η, ζ = ±1/√3. Afin de construire un
champ continu de déformations ou de contraintes, il est nécessaire pour extrapoler les valeurs
de résultat de 2 × 2 × 2 points d'intégration aux sommets de l'élément à 20 nœuds (la
numérotation des points d'intégration et des sommets est illustrée à la Fig. 4.6).

30
 Figure 4.6 : Numérotation des points d'intégration et des sommets
pour l'élément isoparamétrique à 20 nœuds.

Les résultats sont calculés à 8 points d'intégration, et une extrapolation trilinéaire dans la
coordonnée locale le système ξ, η, ζ est utilisé :

où f(m) sont des valeurs de fonction connues aux points d'intégration réduits ; fi sont les
valeurs des fonctions aux nœuds des sommets et Li(m) est la matrice d'extrapolation
symétrique :

Les contraintes sont extrapolées des points d'intégration à tous les nœuds des éléments. Les
valeurs des nœuds médians peuvent être calculées comme une moyenne entre les valeurs de
deux valeurs nodales de sommet. Ensuite, la moyenne des contributions des éléments finis
voisins est effectuée pour tous les nœuds du modèle d'éléments finis. Faire la moyenne produit
un champ continu de résultats secondaires spécifiés aux nœuds du modèle avec une variation
quadratique à l'intérieur des éléments finis. Plus tard, les résultats peuvent être interpolés à
n'importe quel point à l'intérieur de l'élément ou sur sa surface en utilisant des fonctions de
forme quadratiques.

31
Chapitre 5

Discrétisation
5.1 Modèle discret du problème

Afin d'appliquer les procédures d'éléments finis, un modèle discret du problème doit être
présenté dans la forme numérique. Une description typique du problème peut contenir :
 Paramètres scalaires (nombre de nœuds, nombre d'éléments, etc.) ;
 Propriétés matérielles;
 Coordonnées des points nodaux ;
 Tableau de connectivité pour les éléments finis ;
 Tableaux de types d'éléments et de matériaux d'éléments ;
 Tableaux pour la description des conditions aux limites de déplacement ;
 Tableaux pour la description des charges de surface et concentrées ;
Champ de température.
Écrivons des informations numériques pour un problème simple décrit à la Fig. 5.1.

Figure 5.1 : Modèle discret composé de deux éléments finis.

Le modèle par éléments finis peut être décrit comme suit :

1. Paramètres scalaires
Nombre de nœuds = 6
Nombre d'éléments = 2
Nombre de contraintes = 5
Nombre de charges = 2
2. Propriétés des matériaux
Module d'élasticité = 2.0e+8 MPa
Coefficient de Poisson = 0,3
3. Coordonnées du nœud (x1, y1, x2, y2 etc.)
1) 0 0 2) 0 1 3) 1 0 4) 1 1 5) 2 0 6) 2 1
4. Tableau de connectivité des éléments (sens antihoraire)
1) 1 3 4 2 2) 3 5 6 4
5. Contraintes (nœud, direction : x = 1 ; y = 2)
11 21 12 32 52
6. Forces nodales (nœud, direction, valeur)
5 1 0,5 6 1 0,5

32
Alors que pour un exemple simple comme celui démontré ci-dessus, le modèle d'éléments finis
peut être codé à la main, ce n'est pas pratique pour les modèles réels. Divers générateurs de
maillage automatiques sont utilisés pour créer des modèles d'éléments finis pour les formes
complexes.

5.2 Génération de maillage

5.2.1 Générateurs de maillage

Les modèles d'éléments finis pour l'analyse pratique peuvent contenir des dizaines de
milliers voire des centaines de milliers de degrés de liberté. Il n'est pas possible de créer de tels
maillages manuellement. Le générateur de maillage est un outil logiciel, qui divise le domaine
de la solution en plusieurs sous-domaines – éléments finis. Les générateurs de maillage peuvent
être de différents types.
Pour les problèmes bidimensionnels, nous voulons mentionner deux types : les générateurs de
maillage par blocs et les triangulateurs..
Les générateurs de maillage de blocs nécessitent une forme initiale de partitionnement brut. Le
domaine de la solution est partitionné en un nombre relativement petit de blocs. Chaque bloc
doit avoir une forme standard.
Le maillage à l'intérieur du bloc est généralement généré par la technique de cartographie.
Les triangulateurs génèrent généralement un maillage irrégulier dans des domaines arbitraires.
Les polygones de Voronoi et la triangulation de Delaunay sont largement utilisés pour générer
un maillage. Plus tard, le maillage triangulaire peut être transformé
au maillage constitué d'éléments quadrilatéraux. La triangulation de Delaunay peut être
généralisée pour les domaines tridimensionnels.

5.2.2 Technique de cartographie (Mapping technique)

Supposons que nous voulions générer un maillage quadrilatère à l'intérieur d'un domaine qui a
la forme d'un quadrilatère curviligne. La technique de cartographie illustrée à la Fig. 5.2 peut
être utilisée à cette fin.
Si chaque côté du domaine quadrilatère curviligne peut être approximé par une parabole, le
domaine ressemble alors à un élément isoparamétrique à 8 nœuds. Le domaine est mappé sur
un carré dans la coordonnée locale

Figure 5.2 : Génération de maillage avec la technique de cartographie.

système ξ, η. Le carré de coordonnées ξ, η est divisé en éléments rectangulaires puis en

coordonnées nodales sont reconvertis dans le système de coordonnées global x, y.

33
Algorithme de calcul des coordonnées pour le nœud i
nξ = nombre d'éléments dans la direction ξ
nη = nombre d'éléments dans la direction η
Ligne : R = (i − 1)/(nξ + 1) + 1
Colonne : C = mod((i − 1), (nξ + 1)) + 1

Connectivités pour élément e

Élément ligne: R = ( e − 1) /n ξ + 1
Élément colonne: C = mod(( e − 1) , n ξ ) + 1
Connectivités (global nœud Nombres):

Le déplacement des nœuds médians plus près d'un coin du domaine aide à affiner (créer
des éléments plus petits) maille près de ce coin. Si le raffinement est fait du côté de l'élément
qui est parallèle à l'axe local ξ et la taille du plus petit élément près du nœud d'angle est ∆l alors
le nœud médian doit être déplacé vers la position :

Figure 5.3 : Polygones de Voronoï et triangulation de Delaunay.

Ici nξ est le nombre d'éléments dans la direction ξ ; lm est une distance entre le nœud d'angle et
le côté médian node et l est la longueur du côté de l'élément.

34
Chapitre 6

Assemblage et solution
6.1 Assemblage et solution

Le désassemblage est une création de vecteurs d'éléments à partir d'un vecteur global donné.
L'opération de désassemblage est donnée par la relation :

Ici [A] est la matrice assurant la correspondance entre les nombres globaux et locaux de nœuds (ou
degrés de liberté), {Q} est le vecteur global et {Qd} est le vecteur composé des vecteurs d'éléments

Pour les matrices, l'opération de désassemblage n'est pas nécessaire pour la mise en œuvre de la
procédure d'éléments finis.
L'assemblage est l'opération consistant à joindre des vecteurs d'éléments (matrices) dans un
vecteur global (matrice). Pour les vecteurs l'opération d'assemblage est donnée par la relation :

et pour les matrices l'opération d'assemblage est donnée par la relation :

Ici [K] est la matrice globale et [Kd] est la matrice suivante composée des matrices d'éléments :

La fraction d'entrées non nulles (unités) dans la matrice [A] est très petite. Pour cette raison, la
matrice [A] n'est jamais utilisée explicitement dans les codes informatiques réels.
6.2 Algorithme de démontage
L'opération de désassemblage comprend la multiplication matricielle avec l'utilisation
d'une grande matrice [A], qui donne la correspondance entre les énumérations locales et
globales. La matrice [A] est presque entièrement composée de
des zéros. Il n'a qu'une seule entrée non nulle (= 1) dans chaque ligne. Les entrées non nulles
dans [A] fournissent des informations sur les adresses globales où les entrées locales doivent
être prises. Au lieu de la multiplication matricielle, il est possible de simplement prendre des
éléments vectoriels à partir de leurs positions globales et de les mettre aux positions
correspondantes dans le vecteur d'éléments.

35
Démontage pour un vecteur d'élément
n = nombre de degrés de liberté par élément
N = nombre total de degrés de liberté
E = nombre d'éléments
C[E, n] = tableau de connectivité
f[n] = vecteur de charge de l'élément
F[N] = vecteur de charge global
e = numéro de l'élément pour lequel nous avons besoin du vecteur local
faire i = 1, n
f[i] = F [C[e, i]]
fin faire
Dans l'algorithme ci-dessus, nous avons utilisé des informations de connectivité liées aux
degrés de liberté. Si le tableau de connectivité contient des numéros de nœuds globaux et que
chaque nœud a plus d'un degré de liberté, alors au lieu d'un numéro, le bloc lié au nœud doit
être sélectionné. Par exemple, dans un problème d'élasticité tridimensionnel, chaque nœud est
associé à trois déplacements (trois degrés de liberté).

6.3 Assemblage
6.3.1 Algorithme d'assemblage pour les vecteurs
Au lieu de la matrice [A], les procédures d'assemblage sont généralement basées sur une
sommation directe avec l'utilisation du tableau de connectivité des éléments.
Supposons que nous ayons besoin d'assembler le vecteur de charge global F en utilisant les
vecteurs de charge des éléments f et le tableau de connectivité C. Un pseudocode d'algorithme
d'assemblage est le suivant :
Assemblage du vecteur global
n = nombre de degrés de liberté par élément
N = nombre total de degrés de liberté
E = nombre d'éléments
C[E, n] = tableau de connectivité
f[n] = vecteur de charge de l'élément
F [N] = vecteur de charge global
faire i = 1, N
F [i] = 0
fin faire
faire e = 1, E
générer f
faire i = 1, n
F [C[e, i]] = F [C[e, i]] + f[i]
fin faire
fin faire

36
Les vecteurs de charge des éléments sont générés lorsqu'ils sont nécessaires à l'assemblage. On
peut voir que l'entrée de connectivité C[e, i] fournit simplement l'adresse dans le vecteur global
où va la ième entrée du vecteur de charge pour l'élément e. Nous supposons que le tableau de
connectivité est écrit en termes de degrés de liberté. Dans les codes réels, le tableau de
connectivité contient des numéros de nœud qui sont transformés en degrés de liberté pour
l'élément assemblé actuel.

6.3.2 Algorithme d'assemblage pour les matrices

Un algorithme d'assemblage de la matrice de rigidité globale K à partir des
contributions des matrices de rigidité élémentaires k peut être exprimé par le pseudo-code
suivant :
Assemblage de la matrice globale
n = nombre de degrés de liberté par élément
N = nombre total de degrés de liberté dans le domaine
E = nombre d'éléments
C[E, n] = tableau de connectivité
k[n, n] = matrice de rigidité des éléments
K[N, N] = matrice de rigidité globale
faire i = 1, N
faire j = 1, N
K[je, j] = 0
fin faire
fin faire
faire e = 1, E
générer k
faire i = 1, n
faire j = 1, n
K[C[e, je], C[e, j]] = K[C[e, je], C[e, j]] + k[i, j]
fin faire
fin faire
fin faire
Ici, pour des raisons de simplicité, les matrices d'éléments sont entièrement assemblées dans la
matrice globale carrée complète. Comme la matrice de rigidité globale est symétrique et creuse,
ces faits sont utilisés pour économiser de l'espace et du temps dans codes d'éléments finis réels.
6.4 Conditions aux limites de déplacement
Les conditions aux limites de déplacement n'ont pas été prises en compte dans la
fonctionnelle de l'énergie potentielle totale.
Ils peuvent être appliqués au système d'équation global après son assemblage.
Considérons l'application de la condition aux limites de déplacement au système d'équations
globales. Deux méthodes peuvent être utilisées pour la spécification de la condition aux limites
de déplacement.

37
6.4.1 Spécification explicite du déplacement BC

Dans la méthode explicite, nous substituons la valeur connue du déplacement Qm = d dans la

meme colonne et déplaçons cette colonne vers la droite. Ensuite, nous mettons des zéros à la
meme colonne et à la meme ligne de la matrice sauf l'élément diagonal principal, qui est
remplacé par 1.
Méthode explicite:

6.4.2 Méthode des grands nombres

La méthode des grands nombres utilise le fait que les calculs informatiques ont une
précision limitée. Les résultats des calculs en double précision contiennent environ 15 à 16
chiffres. Ainsi, l'addition 1.0+1e−17 produit 1.0 comme résultat..
Méthode des grands nombres ( M >> K ij ):
Kmm = M
F m = Md
La méthode des grands nombres est plus simple que la méthode explicite de spécification des
conditions aux limites de déplacement. La solution du problème des éléments finis est la même
pour les deux méthodes.

6.5 Solution des équations aux éléments finis

6.5.1 Méthodes de résolution

Les applications pratiques de la méthode des éléments finis conduisent à de grands systèmes
d'équations algébriques linéaires simultanées.
Heureusement, les systèmes d'équations aux éléments finis possèdent certaines propriétés qui
permettent de réduire le temps de stockage et de calcul. Les systèmes d'équations aux éléments finis
sont : symétriques, définis positifs et creux. La symétrie permet de ne stocker que la moitié de la
matrice y compris les entrées diagonales. Les matrices définies positives sont caractérisées par de
grandes entrées positives sur la diagonale principale. La solution peut être effectuée sans pivotement.
Une matrice creuse contient plus d'entrées nulles que d'entrées non nulles. La parcimonie peut être
utilisée pour économiser le stockage et les calculs.
Les méthodes de résolution des systèmes d'équations linéaires peuvent être divisées en deux grands
groupes : les méthodes directes et les méthodes itératives. Les méthodes de résolution directe sont
généralement utilisées pour les problèmes de taille moyenne. Pour les grands problèmes, les
méthodes itératives nécessitent moins de temps de calcul et sont donc préférables.
Les formats de stockage matriciel sont étroitement liés aux méthodes de résolution. Ci-dessous, nous
considérons deux méthodes de résolution qui sont largement utilisées dans les codes de calcul par
éléments finis. La première méthode est la solution LDU directe avec matrice de rigidité globale du
profil. La seconde est la méthode du gradient conjugué préconditionné avec un format de stockage
matriciel clairsemé par rangées.

38
 Figure 6.1 : Stockage du profil symétrique de la matrice globale.

6.5.2 Méthode LDU directe avec matrice de profil

En utilisant le format de profil symétrique, la matrice de rigidité globale [A] = Aij d'ordre
N est stockée par colonnes, comme illustré à la Fig. 6.1. Chaque colonne commence par le
premier élément supérieur non nul et se termine au élément diagonal. La matrice est représentée
par deux tableaux :
a[pcol[N + 1]] tableau de doubles contenant des éléments de matrice pcol[N + 1] tableau de
pointeurs d'entiers pour les colonnes.
Le ième élément de pcol contient l'adresse du premier élément de colonne moins un. La
longueur de la ième colonne est donnée par pcol[i + 1]−pcol[i]. La longueur du tableau a est
égale à pcol[N + 1] (en supposant que les indices du tableau commencent à 1).
Par exemple, pour la matrice ci-dessus, le tableau pcol a le contenu suivant :
pcol = {0, 1, 3, 6, 8, 12, 15, 19, 22, 24}
Il convient de noter qu'un bon ordre des nœuds peut réduire considérablement le profil de la
matrice. Habituellement, les algorithmes d'ordonnancement des nœuds sont basés sur certaines
méthodes heuristiques, car le problème d'ordonnancement complet avec recherche du minimum
global prend trop de temps.
Nous allons présenter des algorithmes en notation matricielle complète aij. Ensuite, il est
nécessaire d'avoir des relations entre la notation à deux indices pour la matrice de rigidité
globale aij et le tableau a utilisé dans les codes FORTRAN ou C. L'emplacement du premier
élément non nul dans la ième colonne de la matrice a est donné par la fonction suivante :

Les relations de correspondance suivantes peuvent être facilement obtenues pour une transition
de la notation à deux indices à la notation FORTRAN/C pour un tableau unidimensionnel a :

La solution du système algébrique linéaire symétrique se compose de trois étapes :

Factorisation : [A] = [U]T [D][U]
Solution directe : {y} = [U]−T{b}
Substitution arrière : {x} = [U]−1[D]−1{y}

39
L'algorithme de factorisation à droite d'une matrice de profil symétrique est le suivant :
Factorisation LDU (à droite)
faire j = 2, N
Cdivt(j)
faire i = j, N
Cmod(j, je)
fin faire
fin faire
faire j = 2, N
Div(j)
fin faire
Cdivt(j) =
faire i = F N(j)), j − 1
ti = Aij/Aii
fin faire
Cmod(j, i) =
faire k = max(F N(i), F N(j)), i − 1
Aji = Aji − tkAki
fin faire
Cdiv(j) =
faire i = F N(j)), j − 1
Aij/ = Aii
fin faire
La solution directe avec matrice triangulaire et la rétrosubstitution sont données par le pseudo-
code :
Réduction avant et substitution arrière
faire j = 2, N
faire i = F N(j), j − 1
bj = bj − Aij ∗ bi
fin faire
fin faire
faire j = 1, N
bj = bj/Ajj
fin faire
faire i = N, 1, −1
faire i = F N(j), j − 1
bi = bi − Aij ∗ bj
fin faire
fin faire

40