0.

INTRODUCTION

La modélisation est la représentation d'un système par un autre, plus facile à

appréhender. Il peut s'agir d'un système physique, analogique ou mathématique. Le modèle
sera alors physique, analogique ou mathématique.

Dans ce cours, nous n’étudierons que les modèles mathématiques numériques.

0.1. MODELISATION PHYSIQUE

Les modèles physiques consistent à reproduire le système étudié en modèle

réduit en laboratoire. Les phénomènes à étudier se font sur ce modèle réduit qui a les
mêmes caractéristiques physiques et géométriques que le système réel mais en petite
dimension ; d’où l’appellation, « modèle réduit ».

En hydrogéologie, par exemple, on construit un modèle réduit qui reproduit

l’aquifère à étudier en petites dimensions tout en gardant ses caractéristiques physiques
(hétérogénéité, anisotropie) et géométriques. Il suffit pour cela de remplir une boite avec
des terrains reproduisant en petit ces caractéristiques de l’aquifère, d'y installer des puits,
des "conditions aux limites" (une eau à niveau invariable, un arrosage, etc). Alors on observe
tout ce qui arrive en variant la disposition, les débits ou le nombre des puits en fonction des
conditions aux limites que l’expérimentateur peut aussi varier dans le temps et dans l’espace
à sa guise.

Ces modèles de moins en moins utilisés de nos jours pour les études de terrain
sont cependant parfois encore utilisé à des fins pédagogiques.

0.2. MODELISATION ANALOGIQUE

Analogie entre le courant électrique et l’eau. L’eau et le courant électrique ont

beaucoup de points communs. C’est pourquoi une analogie entres eux est souvent utilisée
pour l’étude du circuit électrique ou de celui de l’eau. La figure (0.2.1), ci-dessous le
démontre.

Figure 0.2.1. Similitude entre le circuit hydraulique et le circuit électrique

1
Dans le descriptif ci-dessous, les mots en rouge indiqueront les organes hydrauliques, ceux
en bleu leurs homologues électriques.

 Si nous suivons le sens de l’eau nous partons d’une pompe G (générateur, alternateur,
turbine, dynamo, pile, batterie, panneau solaire, éolienne).

 Cette pompe refoule dans une conduite, l’eau se trouve alors stoppée par une vanne
hydraulique T (transistor). Pour que l’eau poursuive son chemin il faudra ouvrir le robinet b
(base), d’un circuit annexe, pour que l’eau pousse le piston qui ouvrira la vanne et permette
l’écoulement du fluide dans le circuit principal entre e et c (émetteur et collecteur).

 L’eau traverse ensuite un clapet D (diode), qui permet son écoulement dans un unique sens
et évite ainsi son retour.

 Elle rencontre ensuite un petit réservoir C (condensateur), qui permet un stockage limité.

 Enfin elle arrive dans un grand réservoir de stockage B (batterie ou accumulateur), qui
permet une grande réserve.

 L’eau de par sa hauteur s’écoule ensuite par gravité, la hauteur entre le niveau supérieur de
l’eau et l'axe de la turbine, appelée Hauteur Elévation Totale en hydraulique est l’équivalent
en électricité à la Différence De Potentiel (d.d.p.), U (Tension).

 Le débit I (Intensité), transitant dans ce tuyau ne pourra avoir lieu que si la vanne In
(interrupteur) est ouverte. Il sera proportionnel à la consommation de la turbine Re.

 On notera avant la turbine un rétrécissement R (Résistance) qui limitera le courant.

 La turbine Re (Récepteur) actionnée par l’eau produira une puissance mécanique P

(Puissance), égale au produit de la tension U par l’intensité I.

L’hydrodynamique des aquifères invisibles a longtemps été étudiée par des

circuits électriques dans lesquels les paramètres de l’aquifère étudié sont remplacé par des
éléments correspondants du circuit de simulation.

Mais avec l’arrivée des ordinateurs les analogies électriques pour l’étude des
aquifères ont cédé la place aux modèles mathématiques numériques.

0.3. MADELISATION MATHEMATIQUE

Les modèles mathématiques peuvent être analytiques ou numériques.

0.3.1. Le modèle analytique

Les solutions analytiques sont des procédures logiques qui donnent une solution
exacte.
Les solutions analytiques ne résolvent que les problèmes qui se posent dans des
conditions idéales (milieu homogène, par exemple) et pour des formes géométriques

2
classiques (aquifères, réservoirs de pétrole, nuage etc. de forme rectangulaire,
parallélépipédique, prismatique, pyramidale, sphérique, conique, ellipsoïdal, etc.). La nature
n’offre pas des aquifères, des réservoirs de pétrole, des gisements, de nuages, des lacs etc.
présentant des conditions idéales et des formes géométriques classiques. Les objets étudiés
en science et en ingénierie sont très souvent hétérogènes et même anisotropes et leurs
formes géométriques sont quelconques. Pour trouver des solutions à ces objets, il faut
recourir aux méthodes numériques.

0.3.2. Le modèle numérique

Les solutions numériques sont des procédures d’essai et d’erreur qui sont plus
longues et donnent des solutions approximatives selon la précision voulue par l’opérateur.

Les solutions numériques résolvent par approximations, les problèmes dans les
conditions non idéales (milieu hétérogène, par exemple) et de forme géométrique
quelconque ou même variable au cours de temps. Elles ne trouvent des solutions qu’au
niveau des points choisis par l’opérateur appelés nœuds. Comme l’approche de la solution
exacte s’améliore quand le nombre de nœud augmente et quand le nombre d’essais et
erreurs s’accroît, seul le calcul automatique par ordinateur peut résoudre les problèmes qui
nécessitent des solutions numériques.

Très souvent les équations à résoudre par les méthodes numériques sont des
équations différentielles.

0.4. CONTENU DU COURS

Le cours comprend cinq chapitres : le premier chapitre porte sur l’équation

fondamentale de l’écoulement permanent, le deuxième chapitre établit les équations
fondamentales des écoulements transitoires, le troisième chapitre traite des équations
différentielles, le quatrième chapitre porte sur la modélisation par la méthode des
différences finies et le cinquième chapitre présente l’introduction à la modélisation par la
méthode des éléments finis.

3
1. EQUATION GENERALE DES ECOULEMENTS PERMANTS

L’équation générale de l’écoulement permanent est la combinaison

de l’équation de continuité et de l’équation de Darcy généralisée.

1.1. EQUATION DE CONTINUITE

Considérons l’écoulement permanent à travers un milieu poreux dans lequel

nous avons découpé le parallélépipède de la figure (1.1.1), ci-dessous.

Figure 1.1.1. Parallélépipède imaginé dans un aquifère homogène et isotrope en écoulement permanent

Comme l’eau et le terrain sont incompressibles, le débit entrant dans le

parallélépipède est égal au débit qui en sort. Les côtés dx, dy et dz du Parallélépipède sont
suffisamment petits pour que l’on puisse négliger l’erreur de troncature de la série de Taylor
après le premier terme de la dérivée de la vitesse. Ainsi, on peut écrire :

 La somme des débits entrant : V x dydz + V y dxdz + V z dxdy

 La somme des débits sortant :

(V x+
∂V x
∂x ) (
dx dydz + V y +
∂V y
∂y )
dy dxdz + V z + (
∂V z
∂z
dz dxdy )
On peut donc poser que

∂V x ∂V y ∂
V x dydz +V x dydz +V z dxdy=V x dydz+ V x dydz +V z dxdy+ dxdydz + dxdydz +
∂x ∂y

Donc

∂Vx ∂V y ∂V z
dxdydz + dxdydz+ dxdydz =0
∂x ∂y ∂z

Ou

∂ V x ∂ V y ∂V z
+ + =0 ………………………………………………………………..(1.1.1)
∂x ∂ y ∂z

C’est l’équation de continuité souvent écrite sous la forme vectorielle :

4
 ¿ V =0………………………………………………………………………………………..(1.1.2)

On arriverait à la même équation pour un liquide incompressible si on faisait le bilan sur la

masse liquide qui entre et qui sort dans le cube en multipliant le volume entrant et sortant
par la masse spécifique constante, ρ ,du liquide.

1.2. EQUATION DE DARCY GENERALISEE DE DANS UN MILIEU ISOTROPE

La généralisation de la loi de Darcy est la présentation en plus d’une

dimension de la loi de Darcy :

∂φ
V x =−K ………………………………………………………(1.2.1)
∂x
∂φ
V y =−K ………………………………………………………(1.2.2)
∂y
∂φ
V z=−K ………………………………………………………(1.2.3)
∂z

1.3. ETABLISSEMENT DE L’EQUATION GENERALE DE L’ECOULEMENT PERMANENT

1.3.1. Equation de Laplace en milieu isotrope

Remplaçons les composantes du vecteur vitesse dans l’équation de continuité

(1.1.1) par leurs expressions respectives des équations (1.2.1), (1.2.2) et (1.2.3) qui
expriment la généralisation de l’équation de Darcy en milieu isotrope. Nous aurons :

∂
∂x (
−K
∂φ
+
∂
∂x ∂ y) (
−K
∂φ
+
∂
∂ y ∂z
−K
∂φ
∂z) (
=0 )
Après multiplication par -1 et division par K, nous obtenons :
2 2 2
∂ φ ∂φ ∂ φ
2
+ 2 + 2 =∆ φ=0……………………………………………………………(1.3.1)
∂x ∂ y ∂z

Le potentiel hydraulique, φ , est donc un potentiel harmonique, c’est-à-dire, une

fonction de x, y et z satisfaisant l’équation de Laplace. Cette l’équation est très connue en
mathématique et en physique, l’on appelle, équation de Laplace. En mathématique, une
fonction obéissant à l’équation de Laplace peut être entièrement déterminée et de façon
univoque en tout point d’un domaine, un aquifère par exemple, si on connait sur les limites
du domaine les caractéristiques du potentiel hydraulique.

1.3.2. Conditions aux limites

Lorsqu’on connait les valeurs de la fonction sur la limite, on dit que les conditions
en question sont de Dirichlet.

5
 Lorsque ce sont les dérivées de la fonction qui sont connues sur la limite, on
parle des conditions de Neumann.

Les conditions aux limites qui permettent de résoudre l’équation de

Laplace sont :

a) Sur une limite imperméable

∂φ
=0 (Condition de Neumann)
∂n
n étant un direction perpendiculaire à la limité imperméable

b) Sur une surface filtrante, c’est-à-dire, surface de contact de la nappe aquifère avec
une étendue d’eau libre (un lac, un cours d’eau, l’eau dans un puits etc.)

φ=C sur une verticale (Condition de Dirichlet),

te

c) Sur une surface de suintement, c’est-à-dire, surface où la nappe aquifère est en

contact avec l’air
φ=z

d) Sur une surface piézométrique

φ=z (En effet, sur la surface piézométrique p est nulle)

On a aussi
∂φ
=0
∂n

En effet, la surface piézométrique étant une surface de courant, il est donc imperméable

L’équation générale de l’écoulement souterrain permanent, ici

l’équation de Laplace, donne le champ du potentiel hydraulique en fonction des
coordonnées des différents point d’un aquifère et des conditions aux limites et indique
qu’en milieu isotrope, ce champs n’est pas fonction de la conductivité hydraulique.

1.3.3. Equation de Laplace en milieu anisotrope

Etant donné que le terrain est anisotrope, chaque direction a sa propre

conductivité hydraulique. L’équation de Darcy généralisée s’écrit donc dans le cas d’un
aquifère anisotrope :

6
 ∂φ
V x =K x …………………………………………………………….1.3.2
∂x
∂φ
V y =K y …………………………………………………………….1.3.3
∂y
∂φ
V z=−K z …………………………………………….……………….1.3.4
∂z

L’équation de continuité (1.1.1) ne change pas :

∂ V x ∂ V y ∂V z
+ + =0…………………………………………………………….1.1.1
∂x ∂ y ∂z

En remplaçant les composantes du vecteur vitesse de l’équation de

continuité par leurs expressions de l’équation de Darcy généralisée pour un milieu
anisotrope on obtient l’équation suivante :

2 2 2
∂ φ ∂ φ ∂φ
Kx 2
+ K y 2 + K z 2 =0…………………………………………………………….1.3.5
∂x ∂y ∂z

Le potentiel hydraulique n’est plus une fonction harmonique des variables

x, y et z, car elle n’est plus, par rapport à ces variables, une équation de Laplace. Son champ
dépend, cette fois, de la conductivité hydraulique qui change selon la direction considérée.
Donc on ne peut plus déterminer ce champ à l’aide des conditions aux limites seulement, il
faut, en plus, connaitre les conductivités hydrauliques dans chaque direction.

Pour déterminer le champ du potentiel hydraulique dans un milieu

anisotrope, nous devons remplacer de façon imaginaire, l’aquifère anisotrope réel par un
aquifère isotrope fictif de conductivité hydraulique K constante dans toutes les directions.
Pour cela, il faut replacer par l’imagination l’aquifère étudié dans un système d’axes
cartésiens x’,y’ et z’ tel que :

x ' =x
√ K
Kx
…………………………………………………………….1.3.6

y'= y
√ K
Ky
…………………………………………………………….1.3.7

z ' =z
√
K
Kz
……………………………………………..……………….1.3.8

On démontre que cette transformation, appelée en mathématique,

transformation affine, permet, par la substitution des variables indépendantes x, y et z par
d’autres variables indépendantes x’, y’ et z’, de remplacer l’aquifère réel anisotrope par un
aquifère fictif isotrope de conductivité hydraulique K. Ainsi, le potentiel hydraulique devient
7
une fonction harmonique dans le système de ces nouveaux axes des coordonnées car il est
maintenant une fonction de Laplace de ces nouvelles variables et l’on écrit alors :

2 2 2
∂ φ ∂ φ ∂ φ
2
+ 2
+ 2
=Δ φ=0…………………………………………………………….1.3.9
∂ x' ∂ y ' ∂z '

Une fois le champ du potentiel hydraulique déterminé dans le milieu isotrope fictif à l’aide
de la solution mathématique de l’équation de Laplace en milieu isotrope fictif, on le replace
en milieu anisotrope réelle des coordonnées x, y et z, en faisant les transformations inverses.

Ainsi on

x=x '
√ Kx
K
…………………………………………………………….1.3.10

y= y '
√ Ky
K
…………………………………………………………….1.3.11

z=z '
√Kz
K
……………………………………………….…………….1.3.12

En ces qui concerne les caractéristiques de l’écoulement obtenues dans le milieu isotrope
fictif on les convertit en caractéristiques hydrauliques de l’aquifère réel anisotrope de la
façon suivante :

u=u'
√ Kx
K
…………………………………………………………….1.3.13

v=v '
√ Ky
K
…………………………………………………………….1.3.14

w=w '
√ Kz
K
…………………………………………………………….1.3.15

Q=Q '
√ KxK yKz
K
3
…………………………………………………………….1.3.16

2. EQUATIONS GENERALES DES ECOULEMENTS TRANSITOIRES

On rappelle qu’un écoulement transitoire (écoulement en régime de

non équilibre) est caractérisé par le fait que les caractéristiques de l’écoulement (vitesse,
potentiel hydraulique) changent à chaque instant en chaque point de la nappe aquifère.

2.1. VALIDITE DE LA LOI DE DARCY EN ECOULEMENT TRANSITOIRE

8
 De façon rigoureuse la loi de Darcy aura, en régime transitoire,
l’expression suivante :
∂V
V +α =−K gradφ..............................................(2.1.1)
∂t

Les vitesses des écoulements souterrains sont très faibles, à fortiori

∂V
leurs accélérations, . En négligeant cette accélération, on retrouve l’équation de Darcy
∂t
généralisée en régime permanent :

V =−K gradφ .......................................................(2.1.2)

La loi de Darcy sera donc considérée comme valable à chaque instant

dans un écoulement en régime de non équilibre (régime transitoire).

2.2. EQUATION GENERALE DES ECOULEMENTS A SURFACE LIBRE

2.2.1. Evolution dans le temps de la surface libre

L’écoulement transitoire à surface libre se fait par exemple dans une

digue en terre envahie par l’eau endiguée et drainée lors du vidange du réservoir ou dans
une nappe aquifère à surface libre soumis à un pompage ou à une injection.

L’équation que nous établissons est valable sous deux hypothèses :

a) La loi de Darcy est valable à chaque instant de l’écoulement transitoire

b) L’eau et le terrain qui la contient sont incompressibles.

Nous avons vu au point 1.3.1. que ces deux hypothèses conduise à

l’équation de Laplace :
2 2 2
∂ φ ∂φ ∂ φ
2
+ 2 + 2 = Δ φ=0……………………………………………………………..(1.3.3)
∂x ∂ y ∂z

La seule différence est que cette équation permet de déterminer le champ du potentiel
hydraulique à un instant donné grâce aux conditions aux limites qui existent au bord de
l’aquifère à cet instant-là. Si les conditions aux limites de l’aquifère arrivent à changer à
l’instant d’après, le champ du potentiel hydraulique est recalculé à l’aide des nouvelles
conditions aux limites. Ainsi, dans un écoulement transitoire, le champ du potentiel
hydraulique change à chaque instant du fait de changement continu des conditions aux
limites de l’aquifère dans le temps.

9
 Les conditions aux limites dont il est question sont les mêmes que
celle que nous avons vues au point (1.3.2.) Sauf qu’ici la position de la surface piézométrique
n’est connue qu’à l’instant initial. Pour connaitre la position de cette surface à l’instant
d’après, il faudrait déterminer la vitesse de monter ou de descente de cette surface en
chacun de ses points en un instant donnée.

On détermine cette vitesse par le procédé suivant :

Pour simplifier, considérons que le potentiel hydraulique ne varie que

dans une des directions horizontale, x, par exemple, et dans la direction verticale, z, mais
reste constante dans l’autre direction horizontale, y.

Considérons deux positions successives de la surface libre à l’instant

initial to et à l’instant suivant to + dt (voir figure (2.2.1), ci-dessous). Pour que la surface libre
ait pu passer de sa position au moment to à sa nouvelle position au temps to + dt, elle a dû
être traversée par un débit d’eau Q. Soit dQ le débit qui traverse une surface élémentaire dS
de la surface libre de normale n . Soient ε , V n et dn respectivement le coefficient
d’emmagasinement (porosité efficace) de l’aquifère, la composante selon n du vecteur
vitesse de l’eau qui traverse la surface libre et le segment de la normale compris entre les
deux lignes qui déterminent la position de la surface libre à ces deux instants.

La croissance dh de h pendant le temps de dt de la remontée de la

surface piézométrique est

∂h
dh= dt ....................................................................................(2.2.1)
∂t

∂h
Avec la vitesse instantanée de la remontée de la surface libre de
∂t
la nappe aquifère au point considérée de cette surface.

Figure 2.2.1. : Variation dans le temps de la surface piézométrique

10
Du fait du principe de continuité, le volume d’eau qui est entré dans la nappe par la surface
élémentaire dS pour relever la surface libre pendant le temps dt est égal au volume qui a
rempli les vides efficaces du terrain aquifère pendant ce temps :

V n dSdt=εdSdn .....................................................................................(2.2.2)

Donc

V n dt=εdn ..............................................................................(2.2.3)

Or

∂h
dn= dtcos( α )..................................................................................(2.2.4)
∂t

α étant l’angle que fait dn avec la verticale.

Par ailleurs

V n=V x sin ( α ) +V z cos(α ).....................................................................(2.2.5)

Donc

∂h
V x sin ( α ) +V z cos ( α ) =ε cos (α )............................................................(2.2.6)
∂t

Ou

∂h
V x tg ( α ) +V z=ε ....................................................................................(2.2.7)
∂t

Or

−∂ h ∂φ ∂φ
tg ( α ) = ; V x =−K ; V z =−K
∂x ∂x ∂z

L’équation en deux dimensions x et z de la vitesse instantanée de

remontée de la surface libre en chacun de ses points est finalement :

= −(
∂h K ∂φ ∂h ∂φ
∂t ε ∂x ∂x ∂z )
.................................................................................(2.2.8)

En trois dimensions x, y et z, cette équation s’écrit :

(
∂h K ∂φ ∂h ∂h ∂φ ∂φ
= + −
∂t ε ∂x ∂x ∂ y ∂ y ∂z )
.............................................................................(2.2.9)

Où

11
  K est la conductivité hydraulique de l’aquifère homogène et isotrope
 ε est le coefficient d’emmagasinement (porosité efficace) de l’aquifère
 φ est le potentiel hydraulique au point de la surface piézométrique où on voudrait
connaître la vitesse de la remontée ou de la descente de la surface piézométrique
 h e st la hauteur en ce point de la surface piézométrique.

∂φ
Etant donné qu’à la surface piézométrique φ = h, cependant qui est le gradient
∂x
∂h
hydraulique à la surface libre n’est plus égal à
qui est la pente de la surface libre. Les
∂x
deux ne sont les mêmes que si l’écoulement est horizontal, c’est-à-dire, si les surfaces
équipotentielles sont verticales.

A l’aide de cette vitesse instantanée du mouvement (remontée ou

descente) de la surface libre, on peut déterminer la position et la forme de cette surface à
une date future donnée. On doit d’abord retenir que le mouvement de la surface libre se fait
à une vitesse variable dans le temps ; donc on ne connait à chaque instant que la vitesse
instantanée qui va changer, augmenter ou diminuer au cours de temps.

2.2.2.Equation générale des écoulements dans une nappe phréatique peu épaisse

L’épaisseur de la nappe est si faible par rapport à son étendue que la

composante verticale du vecteur vitesse de l’eau est considérée comme quasi nulle (la
direction de l’écoulement est quasi horizontal et perpendiculaire aux surfaces
équipotentielles quasi verticales). Ainsi, la variation du potentiel hydraulique sur l’axe
verticale z est très faible ; les surfaces équipotentielles étant presque verticales. Alors on a

∂φ ∂h ∂φ ∂h
= et = .
∂x ∂x ∂ y ∂y

Et

Or

[ ] [ ]
h 2 h 2 2 2 2
∂φ ∂ φ ∂ φ ∂φ ∂ h ∂h
=∫ 2 dz =−∫ 2
+ 2 dz=−h 2
+ 2 ..........................................(2.2.10)
∂z 0 ∂z 0 ∂x ∂y ∂x ∂ y

L’équation de la vitesse instantanée de la surface libre s’écrit alors :

∂h K ∂
= h
∂h
+( ) ( )
∂
h
∂t ε ∂x ∂x ∂ y ∂h
∂h
......................................................................(2.2.11)

12
Comme l’épaisseur de la nappe aquifère et très faible par rapport à l’étendue de celle-ci,
nous pouvons considérer h comme l’épaisseur constante, H, de la nappe et récrire l’équation
ci-dessus :

( )
2 2
∂ h KH ∂ h ∂ h
= + ......................................................................................(2.2.12)
∂t ε ∂ x2 ∂ y 2

Ou

( )
2 2
∂h T ∂ h ∂ h
= + .........................................................................................(2.2.13)
∂ t ε ∂ x2 ∂ y2

Ou

∂h T
= ∆ h................................................................(2.2.14)
∂t ε

Avec T, la transimissivité de la nappe aquifère qui est le produit de la conductivité par

l’épaisseur de la nappe.

C’est l’équation de la chaleur très connue en physique.

2.3. EQUATION GENERALE DES ECOULEMENTS TRANSITOIRES DANS UNE NAPPE

PROFONDE SOUS PRESSION (NAPPE CAPTIVE)

Lorsque la nappe aquifère n’est plus phréatique mais profonde est

sous pression et d’épaisseur e et de porosité efficace, m, l’équation générale de
l’écoulement transitoire prend en compte la compressibilité de l’eau, β , et de la roche
magasin, α .

La variation dans le temps du potentiel hydraulique s’écrit toujours

en équation de la chaleur :

∂ φ Ke T
= ∆ φ= ∆ φ ...................................................................(2.3.1)
∂t S S

Avec

 S= ρg(α +nβ)est le coefficient d’emmagasinement de l’aquifère. Il est sans

dimension et représente la quantité d’eau qu’une colonne de section unitaire
(découpée mentalement dans l’aquifère sous pression) peut libérer lorsque le
potentiel hydraulique, φ , baisse d’une unité.

13
  T =eK est la transimissivité de l’aquifère.

Le coefficient d’emmagasinement et la transimissivité sont

déterminés par des essais de pompage.

3. EQUATIONS DIFFERENTIELLES

3.1. DEFINITION

En mathématiques on distingue, d’une part, des équations différentielles

ordinaires (parfois simplement appelée équation différentielle et abrégée en EDO) dont la ou
les fonctions inconnues ne dépendent que d'une seule variable; elle se présente sous la
forme d'une relation entre ces fonctions inconnues et leurs dérivées successives et, d’autre
part, des équations aux dérivées partielles (parfois appelées équations différentielles
partielles et abrégée en EDP) dont les solutions sont les fonctions inconnues dépendant de
plusieurs variables vérifiant certaines conditions concernant leurs dérivées partielles.

Une EDP a souvent de très nombreuses solutions, les conditions étant moins
strictes que dans le cas d'une équation différentielle ordinaire à une seule variable ; les
problèmes comportent souvent des conditions aux limites qui restreignent l'ensemble des
solutions. Alors que les ensembles de solutions d'une équation différentielle ordinaire
présentent plusieurs paramètres correspondant aux conditions supplémentaires, dans le cas
des EDP, les conditions aux limites se présentent plutôt sous la forme de fonction ;
intuitivement cela signifie que l'ensemble des solutions est beaucoup plus grand, ce qui est
vrai dans la quasi-totalité des problèmes.

Les EDP sont omniprésentes dans les sciences puisqu'elles apparaissent aussi
bien en dynamique des structures ou en mécanique des fluides que dans les théories de la
gravitation, de l'électromagnétisme (équations de Maxwell), ou des mathématiques
financières (équation de Black-Scholes). Elles sont primordiales dans des domaines tels que
la simulation aéronautique, la synthèse d'images, ou la prévision météorologique. Enfin, les
équations les plus importantes de la relativité générale et de la mécanique quantique sont
également des EDP.

3.2. CLASSIFICATION DES EQUATIONS DIFFERENTIELES

Les équations différentielles sont classées suivant plusieurs critères que nous
allons examiner.

3.2.1. Classification selon le nombre de variables indépendantes

Sous ce critère, comme nous venons de le dire, on distingue généralement deux

types d'équations différentielles :

 Les équations différentielles ordinaires (EDO) où la ou les fonctions inconnues ne

dépendent que d'une seule variable. Considérons l’écoulement vers un puits d’eau

14
 ou de pétrole et la fonction inconnue recherchée est la pression du réservoir de
pétrole ou le potentiel hydraulique d’un aquifère, si l’écoulement est permanent et
dans une roche réservoir isotrope, les deux fonctions ne dépendront que de la
distance au puit (r), si l’écoulement est transitoire, en un point, ces deux fonctions
ne dépendront que du temps.

Figure 3.2.1. Milieu isotrope, la fonction inconnue dépendant d’une seule variable indépendante, r

 Les équations différentielles partielles, plutôt appelées équations aux dérivées

partielles (EDP), où la ou les fonctions inconnues peuvent dépendre de plusieurs
variables indépendantes. Considérons l’écoulement dans un réservoir anisotrope
vers un puits de pétrole ou d’eau situé à l’origine d’un système des coordonnées
cartésiennes x et y dans un plan perpendiculaire à l’axe du puits où la fonction
inconnue recherchée est la pression de pétrole ou le potentiel hydraulique ; si
l’écoulement est permanent, la fonction recherchée dépendra de deux variables
indépendantes : la distance au puits (r) et de l’angle que fait le rayon avec un des
axes des coordonnées. Si l’écoulement est transitoire, la fonction recherchée
dépendant trois variables indépendantes : la distance au puits, l’angle que fait le
rayon avec un des axes des coordonnées et du temps de production.

Figure 3.2.2. Milieu anisotrope ; la fonction recherchée dépendant de deux variables indépendantes : la
distance, r et l’angle, a.

3.2.2. Classification selon l’ordre de dérivée

15
 L’ordre d’une équation différentielle est celui de la dérivée la plus élevée qu’elle
contient. Des exemples :

 Cette EDP est d’ordre 1 :

∂u ∂u
−b =0..................................................................(3.2.1)
∂x ∂y
 Cette EDP est d’ordre 2 :

2
∂ u ∂u
2
− =0..................................................................(3.2.2)
∂x ∂ y

 Cette EDP est d’ordre 4 :

3 4
∂ u ∂ u
3
− 4 =0 ..................................................................(3.2.3)
∂x ∂ y

Lorsque plusieurs EDPs indépendantes forment un système, l’ordre de leur

système est obtenu en combinant toutes les équations en une seule EDP. A titre d’exemple,
le système des EDPs suivantes est d’ordre 2 alors que chacune de EDPs qui le forme est
d’ordre 1 :

}
∂u ∂ v ∂u
+ =
∂x ∂ y ∂ z
2 2 2
∂w ∂ w ∂ w ∂ w
u= + = ..................................................................(3.2.4)
∂x ∂ x2 ∂ y2 ∂ x ∂ z
∂w
v=
∂y

3.2.3. Classification selon les constantes

Considérons cette EDP :

2 2 2
∂ u ∂u ∂u ∂u ∂u
A 2 +B +C 2 +D +E + Fu=G....................................................(3.2.5)
∂x ∂ x ∂ y ∂y ∂ x ∂ y

 Si les constantes A, B, C, D, E et F sont indépendantes et ne dépendent ni des

variables indépendantes x et y ni de la variable dépendante u, la EDP est dite linéaire.
 Une équation aux dérivées partielles quasi-linéaire est linéaire par rapport aux
dérivées partielles d'ordre le plus élevé de la fonction u.
 Si les constantes A, B, C, D, E et F dépendent de la variable dépendante u, l’EDP est
dite non linéaire.

16
3.2.4. Classification selon la grandeur B2−4 AC

 Si B2−4 AC <0, l’EDP est qualifiée d’elliptique A = 1 ; B = 0 et C =1

2 2
∂ u ∂ u
Exemple, l’équation de Laplace en deux dimensions : 2
+ 2 =0
∂x ∂y
A = 1 ; B = 0 et C =1
D’où : 02 −4∗1∗1< 0

 Si B2−4 AC >0, l’EDP est qualifiée d’hyperbolique

2 2
∂ u 1 ∂u
Par exemple, l’équation d’onde en une dimension : 2 − 2 2 =0
∂ x c ∂t
−1
A = 1 ; B = 0 et C= 2
c
2 4∗1∗−1 4
D’où : 0 − 2
= 2 >0
c c

 Si B2−4 AC =0 l’EDP est qualifiée de parabolique

2
∂u ∂ u
Par exemple, l’équation de la chaleur : =
∂ t ∂ x2
A = 1 ; B = 0 et C=0
D’où : 02 −4∗1∗−0=0−0=0

3.3. QUELQUES EQUATIONS DUFFERENTIELLES TRES CONNUES

3.3.1. Equation des cordes vibrantes

2 2
∂ y 2∂ y
2 ..................................................................(3.3.1)
2
=a
∂t ∂x

Elle s’applique aux petites vibrations transversales d’une corde tendue, flexible,
du type corde de la guitare ou du violon, initialement tendue le long de l’axe des x et mise en
mouvement. La fonction y(x,t) représente le déplacement d’un point x quelconque de la
τ
corde à l’instant, t. La constante, a 2, a pour valeur, , où τ est la tension, constante, de la
ρ
corde et ρ , la masse (constante, par unité de longueur de la corde. On suppose qu’il n’y a pas
de force extérieure à la corde et que celle-ci ne vibre que par élasticité.

17
 Cette équation peut être facilement généralisée à un ordre de dimensions plus
élevées comme le cas d’une membrane ou d’une peau de tambour à deux dimensions, par
exemple. A deux dimensions elle prend alors la forme suivante :

( )
2 2 2
∂ z 2 ∂ z ∂ z
2
=a 2
+ 2 ..................................................................(3.3.2)
∂t ∂x ∂ y

3.3.2. Equation de propagation de la chaleur

( )
2 2 2
∂u ∂ u ∂u ∂u
=ΚΔu=Κ 2
+ 2 + 2 ..................................................................(3.3.4)
∂t ∂x ∂ y ∂z

Ici, u est la température en un point (x,y,z) d’un solide à l’instant, t. La constante, Κ , appelée
C
constante de diffusion, et est égale à ; β et ρ sont respectivement la chaleur spécifique et
βρ
la masse volumique supposées constantes du solide.

3.3.3. Equation de Laplace

2 2 2
∂ u ∂ u ∂u
2
+ 2 + 2 =0 ..................................................................(3.3.5)
∂ x ∂ y ∂z

On retrouve cette équation dans de nombreux domaines ; par exemple, dans la

théorie de la propagation de la chaleur, u, est la température stationnaire quand, après un

long temps, elle cesse de changer avec le temps

∂u
∂t ( )
=0 ; elle se répartit alors dans le corps

uniquement suite aux conditions aux limites de celui-ci. Dans le cas des aquifères et des
réservoirs de pétrole u est respectivement le potentiel hydraulique ( φ ) ou la pression du
réservoir (p).

3.3.4. Equation de la vibration longitudinale d’une poutre

2 2
∂u 2∂ u
2
=c 2 ..................................................................(3.3.6)
∂t ∂x

Où

 u(x , t): déplacement vibratoire de la poutre le long de l’axe des x

2 E
 c=
ρ
 E : module de Young

18
  ρ ; masse volumique de la poutre

3.3.5. Equation de la vibration transversale d’une poutre

2 4
∂ y 2∂ y
4 ..................................................................(3.3.7)
2
=b
∂t ∂x

Où

 y (x , t): déplacement transversale (faible) d’un point x de la poutre par rapport à

l’axe des x
2 EI
 b=
Aρ
 E : module de Young
 I : le moment par rapport à l’axe des x d’une section droite
 A : la surface de ladite section droite
 ρ ; masse par unité de longueur de la poutre

4. PRINCIPES DE MODELISATION PAR LA METHODE DES DIFFERENCES FINIES

La méthode des différences finies consiste à déterminer les valeurs de la fonction

définie par des dérivées partielles en des points choisis par l’opérateur, appelés « nœuds ».
La détermination de ces valeurs d’un point à un autre se fait par le développement des
séries de Taylor.

En revanche, la méthode des éléments finis utilise des polynômes continus par
morceaux pour interpoler entre les points nodaux. Bien que les points ou les nœuds jouent
un rôle dans la théorie des éléments finis, l'accent est mis davantage sur les fonctions
d'interpolation. Nous introduirons les éléments de base de la méthode des éléments finis au
chapitre cinq.

En différences finies, le domaine d’étude est découpé en mailles généralement

rectangulaires dont les sommets forment les nœuds. Connaissant l’équation différentielle
permettant de calculer la valeur de la fonction dérivée en chaque point du domaine étudié,
ces dérivées partielles donnent donc la valeur de la fonction à chaque nœud. L’extrapolation
des valeurs d’un nœud à l’autre se fait par le développement de la série de Taylor.

4.1. APPROVIMATION DES EDPs PAR DES DIFFERENCES FINIES

19
4.1.1. Espace monodimensionnel

Considérons u (x), comme une fonction continue de la variable indépendante

unique x (u qui peut être le potentiel électrique, le potentiel hydraulique, la pression de
l’hydrocarbure dans un réservoir, la concentration chimique d’une substance, etc. et x peut
être le temps ou la distance). Nous discrétisons le domaine x (voir figure 4.1.1) en un
ensemble de points (ou nœuds) tels que

u ( x r ) ≡u ( rh ) ≡u r..............................................................................(4.1.1)
Avec
 r =0 ,1 , 2 , …
 h: un réel positif inférieur à 1

Figure 4.1.1. Variation de la fonction, u, dans un espace monodimensionnel

En remplaçant la localisation de xr par rh, les coordonnées nodales sont

spécifiées simplement comme le produit de l'entier r et de l'espacement de grille h (ici h est
supposé constant et positif normalisé à moins que l'unité). L'entier r désigne la position du
nœud le long du domaine x par rapport à une donnée spécifiée, généralement r = 0 lorsque
x=0. Lorsque h est une constante, u (rh) peut être représenté sans ambiguïté par u r. Voir
figure (4.1.2).

Pour résoudre numériquement les EDPs les dérivées sont remplacées par des
différences finies en tronquant la série de Taylor avec une erreur de troncature que
l’opérateur choisit.

Par exemple, on peut exprimer la dérivée première par les différences finies
comme suit à partir de séries de Taylor indiquées pour le développement de u(x) à droite
(équation (4.1.2)) et à gauche (équation (4.1.3)) :

2 2 3 3
∂u h ∂ u h ∂ u
u ( x r +h )=u ( x r ) +h + + + ¿.............................................................(4.1.2)
∂ x 2! ∂ x 2 3 ! ∂ x 3

2 2 3 3
∂u h ∂ u h ∂ u
u ( x r−h )=u ( x r )−h + − +¿............................................................(4.1.3)
∂ x 2! ∂ x 2 3 ! ∂ x 3
20
On tire la dérivée première à partir de l’équation (4.1.2) :

∂ u u ( x r +h ) −u ( x r ) h ∂ u h ∂ u
2 2 3
= − − −… ..............................................................(4.1.4)
∂x h 2! ∂ x 2 3! ∂ x3

Comme h est petit (h < 1) on peut tronquer l’équation (4.1.4) et écrire l’approximation
suivante :

∂ u u ( x r +h ) −u ( x r ) ou ∂ u ur +1−ur ....................................................(4.1.5)
= =
∂x h ∂x h

On tire la dérivée première à partir de l’équation (4.1.3)

∂ u u ( x r )−u ( x r −h ) h ∂ u h ∂ u
2 2 3
= + − −¿ ..............................................................(4.1.6)
∂x h 2! ∂ x 2 3! ∂ x3

Si h est petit on peut tronquer l’équation (4.1.6) et écrire l’approximation suivante :

∂ u u ( x r )−u ( x r −h ) ou ∂ u ur −ur−1 ............................................................(4.1.7)

= =
∂x h ∂x h

Si nous additionnons les équations (4.1.5) et (4.1.7) nous pouvons écrire :

∂ u ur +1−ur −1
= ...........................................................(4.1.8)
∂x 2h
2 2
h h
Avec le premier terme de la troncature qui est de = et nous écrivons que l’erreur de la
3! 6
2 2
troncature est d’ordre h , ou O(h ).

Si nous soustrayons l’équation (4.1.7) de l’équation (4.1.5), nous obtenons :

∂ u ur +1−2 ur +ur −1
2

2
= 2 ...........................................................(4.1.9)
∂x h
2 2
2h h
Avec le premier terme de la troncature qui est de = et nous écrivons que l’erreur de
4 ! 12
2 2
la troncature est d’ordre h , ou O(h ).

Nous pouvons continuer les opérations et générer plusieurs autres approximations mais les
opérations deviennent de plus en plus fastidieuses et longues. Nous résumons, dans le

21
tableau (4.1.1), ci-dessous, les approximations des dérivées des équations
monodimensionnelles.

∂ u ur +1−ur
= O(h) ................................................................(4.1.10)
∂x h

∂ u ur −ur−1
= O(h) .................................................................(4.1.11)
∂x h

∂ u ur +1−ur −1
= O(h2) ................................................................(4.1.12)
∂x 2h

∂ u −ur +2 +4 u r+ 1−3 ur
= O(h2) .....................................................(4.1.13)
∂x 2h

∂ u −ur +2 +8 ur +1−8 ur −1−ur−2

= O(h3) ..........................................(4.1.14)
∂x 2h
∂ u ur +1−2 ur +ur −1
2

2
= 2 O(h) ...........................................................(4.1.15)
∂x h

∂ u −ur +2 +16 ur +1−30 ur +16 ur −1−ur −2

2

2
= 2 O(h4) .............................................(4.1.16)
∂x 12 h

∂ u ur +2−2 ur +1+ 2u r−1−ur −2

3

3
= 2 O(h2) ............................................................(4.1.17)
∂x 2h

∂ u u r+2−4 ur +1 +6 ur −4 ur−1 +ur −2

4

4
= 4 O(h2) .............................................(4.1.18)
∂x h

Tableau 4.1.1. Les approximations les plus utilisées en différences finies dans un espace monodimensionnel

4.1.2. Espace bidimensionnel

La variable dépendante, u, est fonction de x et de y. Les distances entre deux

nœuds consécutives sont h selon l’axe des x et k selon l’axe des y. Les variables seront
indicées par r et par s respectivement selon l’axe des x et selon l’axe des y. Voir figure (4.1.2)

22
 Figure 4.1.2. Variation de la fonction, u, dans un espace bidimensionnel et tridimensionnel

Si nous considérons l’évolution de la fonction u selon l’axe des x seulement alors

sa dérivée selon l’axe des y est nulle et nous réécrivons l’équation (4.1.19) :

∂ u ur +1−ur
= O(h) ......................................................................(4.1.19)
∂x h

En conformité avec la figure (4.1.1 a), nous écrivons :

∂ ur , s u r+1 , s−u r , s
= O(h) ................................................................(4.1.20)
∂x h

Et l’équation (1.4.10) :

∂ ur , s u r , s−ur −1 , s
= O(h) .................................................................(4.1.21)
∂x h

L’équation bidimensionnelle équivalente de l’équation monodimensionnelle (4.1.15)

∂ u ur +1−2 ur +ur −1
2

2
= 2 O(h) ...........................................................(4.1.15)
∂x h

est
2
∂ ur, s ur +1 , s−2 ur , s +ur −1 , s
2
= 2
...........................................................(4.1.22)
∂x h

Et si u évolue seulement selon la direction y, alors r reste constant :

2
∂ ur, s ur ,s +1−2 ur , s +ur ,s−1
2
= 2
...........................................................(4.1.23)
∂y k

23
Considérons maintenant que u varie dans les deux directions x et y en même temps. Dans ce
cas, nous examinerons l’approximation en différence finie de la dérivée suivante :

[ ]
2
∂u ∂ ∂u
= .............................................................................(4.1.24)
∂ x∂ y ∂x ∂ y

Par analogie avec l’équation (4.1.12),

∂ u ur +1−ur −1
= + O(h2) ................................................................(4.1.12)
∂x 2h

l’équation (4.1.24) peut s’écrire :

[
1 ∂u r+ 1, s ∂ ur−1 , s
]
2
∂u
+O ( h ) ...........................................(4.1.25)
2
= −
∂ x ∂ y 2h ∂ y ∂y

Dans la même analogie, les dérivées de u par rapport à y peuvent s’écrire en différence finie
comme :

[ ]
3 3
1 ur +1 ,s +1−ur +1 ,s −1 k ∂ u r+1 , s ur−1 , s+1−u r−1 , s−1 k 2 ∂ ur−1 , s
2 2
∂u
+ … +O ( h )......
2
= + 3
+ …− − 3
∂ x ∂ y 2h 2k 3! ∂ y 2k 3! ∂ y
..(4.1.26)

Introduisons l’approximation de la différence centrale :

( )
4 3 3
∂ ur, s 1 ∂ ur +1 ,s ∂ ur−1 , s
+O ( h )..........................................(4.1.27)
2
= −
∂ x∂ y 3
2 h ∂ y3 ∂y 2

La substitution de l’équation (4.1.27) dans l’équation (4.1.26) donne :

[
1 ur +1 ,s +1−ur −1 , s+1−u
]
2
∂u −ur −1, s −1
+O ( h ) +O ( k ) ......................(4.1.28)
2 2
= r −1, s +1

∂ x ∂ y 2h 2k

Si h = k, alors l’équation (4.1.28) se réduit à :

2
∂u 1
=
∂ x∂ y 4h
[ ur +1 ,s +1−ur −1, s +1−u r −1, s + 1
−ur−1, s−1 ] +O ( h2 )......................(4.1.29)

Les approximations les plus fréquemment utilisées en deux dimensions sont présentées ci-
dessous dans le tableau (4.1.2) :

24
 ∂ ur , s u r+1 , s−u r , s
= O(h)......................................................................(1.4.22)
∂x h

∂ ur , s u r , s−ur −1 , s
= O(h) ......................................................................(1.4.23)
∂x h

∂ ur , s u r+1 , s−u r−1 ,s 2

= O (h ) ......................................................................(1.4.23)
∂x 2h

∂ u −ur +2 , s+ 4 ur +1 , s−3 ur ,s
= O(h2) .....................................................(1.4.24)
∂x 2h

∂ u ur +1 ,s +1−ur −1 ,s +1+u r+ 1, s−1−u r−1 ,s−1

= O(h2) ..........................................(1.4.25)
∂x 4h

∂ u ur +1 ,s −2u r , s+u r−1 , s

2

2
= 2 O(h2) ...........................................................(1.4.26)
∂x h

∂ u −ur +2 ,s +16 u r+1 , s−30 ur ,s +16 u r−1 ,s −ur−2 , s

2

2
= 2 O(h4) ..................................(1.4.27)
∂x 12 h

∂ u u r+2 , s−4 ur +1 , s+ 6u r , s−4 ur−1 , s +ur−2 , s

4

4
= 4 O(h2) .............................................(1.4.28)
∂x h

∂u
2
u r+ 1, s +1−ur +1 ,s−1 +6 ur ,s −ur−1 , s+1 +ur −1 ,s−1
= 4 O(h2) .................................(1.4.29)
∂ x∂ y h

Tableau 4.1.2. Les approximations les plus utilisées en différences finies dans un espace bidimensionnel.

4.1.3. Espace tridimensionnel

Les approximations des dérivées dans un espace en deux dimensions peuvent

facilement s’étendre à un espace tridimensionnel x, y et z. Dans ce cas, les symboles spatiaux
sont r pour la direction x, s pour la direction y et t pour la direction z. On peut faire les
mêmes approximations dans un espace à quatre dimensions dont la quatrième dimension
est le temps x, y, z et t ( t ≡ temps). Les approximations dans trois dimensions sont faites dans
la même analogie avec celles des équations (1.4.18), (1.4.20) et (1.4.22). Par exemple,
l’équation différentielle de second ordre s’écrit en trois dimensions comme suit :
2
∂ u r , s ,t ur ,s , t+ 1−2u r , s ,t +u r , s ,t −1
+O ( l ).................................................(4.1.30)
2
2
= 2
∂z l

25
4.2. CONCEPTS ADDITIONNELS

Pour un essai d’application des concepts vus dans les matières précédentes,
considérons l’équation différentielle parabolique, par exemple, l’équation de la chaleur
(3.3.4) :

( )
2 2 2
∂u ∂ u ∂u ∂u
=ΚΔu=Κ 2
+ 2 + 2 ..................................................................(3.3.4)
∂t ∂x ∂ y ∂z

En deux dimensions, nous pouvons écrire :

( )
2 2
∂u ∂ u ∂u
=ΚΔu=Κ 2
+ 2 ..................................................................(4.2.1)
∂t ∂x ∂ y

Et les équations elliptiques, par exemple, l’équation de Laplace (3.3.5) :

2 2 2
∂ u ∂ u ∂u
2
+ 2 + 2 =0 ..................................................................(3.3.5)
∂ x ∂ y ∂z

En deux dimensions, nous pouvons écrire :

2 2
∂ u ∂ u
2
+ 2 =0 ..................................................................(4.2.2)
∂x ∂y

Ou l’équation elliptique de Poisson en deux dimensions :

2 2
∂ u ∂ u
2
+ 2 =f ( x , y )..................................................................(4.2.3)
∂x ∂y

Nous ignorons les conditions aux limites et nous nous concentrons uniquement
sur les approximations en différences finies de ces équations différentielles.

4.2.1. Approximation de l’équation différentielle parabolique (équation de la chaleur)

Considérant les approximations vues dans les tableaux (4.1.1) et (4.1.2) nous
pouvons l’écrire :

[ ]
2
∂u ∂ u r , s u r+1 , s−ur ,s ur ,s +1−2 ur , s +ur , s−1
+O ( h , k ) =0 .....................................(4.2.4)
2
− 2
= − 2
∂x ∂ y h k

26
 Figure 4.2.1. Différents cas des nœuds des équations différentielles vues ci-dessus

Négligeant l’erreur de troncature, l’équation de la différence finie pour le nœud s au niveau

r+1 devient :

ur +1 ,s−ur , s u r , s+1−2 ur ,s +ur ,s−1

= 2
h k
Et finalement
h
2 ( r ,s +1
ur +1 , s= u −2 ur , s +ur ,s−1 ) +u r , s
k
Ou

ur +1 , s=
h
k
u
2 r ,s +1 ( h h
+ ur , s 1−2 2 + 2 u
k )
k r , s−1

h
Si nous posons, pour simplifier, 2
=ρ , nous pouvons écrire :
k

ur +1 , s=ρ ur , s+ 1+u r , s ( 1−2 ρ )+ ρur , s−1............................................(4.2.5)

L’erreur est de l’ordre de h et de k2 quand h et k sont faibles.

L’équation (4.2.5) exprime en différences finies l’équation de la chaleur au nœud

r+1,s de la figure (4.2.1 a).

1
Dans le cas spécial où ρ= , l’équation (4.2.5) se simplifie :
2

ur , s+1 +ur ,s−1

ur +1 , s= ..............................................................(4.2.6)
2

27
L’équation (4.2.6) montre qu’à l’erreur de troncature près, la valeur de u r+1,s est la moyenne
de ur,s+1 et ur,s-1. Dans ce cas spécial, le schéma des nœuds de la figure (4.2.1 a) se présente
comme si le nœud r,s n’existait pas, figure (4.2.1 b) :

4.2.2. Approximation de l’équation différentielle elliptique (équation de Laplace)

4.2.2.1. Equation de Laplace

L’équation de Laplace peut être présentée par approximation en différence finie

comme suit :

[ ]
2 2
∂ ur ,s ∂ ur ,s ur +1 , s−2 ur , s +ur −1 , s ur ,s +1−2 ur , s +ur ,s−1
+ O ( h + k )=0............(4.2.7)
2 2
2
+ 2
= 2
+ 2
∂x ∂y h k

Après réarrangement, l’équation (4.2.6) peut s’écrire :

1 1 2 ur , s ( h2 +k 2 )
2 ( r +1 , s
u +u r−1 , s ) + 2 ( r ,s +1
u +u r , s−1 )= 2 2
+O ( h2 +k 2 ).......................(4.2.8)
h k h k

Le diagramme des nœuds de l’équation (4.2.8) est présenté à la figure (4.2.1 c), ci-dessus.

Lorsque h = k, l’équation (4.2.8) se réduit à :

2
ur +1 , s+ ur−1 , s+ ur , s+1 +ur ,s −1=4 ur ,s +O(h ) .............................................(4.2.9)

Son erreur de troncature est de l’ordre de h2.

Considérons maintenant le cas général où les écartements des nœuds sont tous
différents comme dans la figure (4.2.1 d), ci-dessus. Nous voudrons établir une
approximation qui satisfait l’équation de Laplace. La procédure que nous adoptons est la
méthode des coefficients indéterminés dans lesquels nous choisissons le nœud pour lequel
nous voulons déterminer l’approximation par différences finies de l’équation de Laplace.
Ainsi, nous posons :

2 2
∂ u r , s ∂ ur , s i=4
2
+ 2
=∑ α i ui .......................................................(4.2.10)
∂x ∂y i=0

L’équation (4.2.10) peut être interprétée comme une combinaison linéaire pondérée de cinq
points discrets. Les pondérations α 0 , α 1 … . α 4 sont déterminées de façon explicite par le
développement des séries de Taylor. Ces développements sont :
28
 2 2 3 3
∂ u x h1 ∂ u x h1 ∂ u x
u1=u 0+ h1 + + +…........................................(4.2.11)
∂ x 2! ∂ x 2 3 ! ∂ x 3
2 2 3 3
∂ u y h2 ∂ u y h2 ∂ u y
u2=u 0+ h2 + + +…........................................(4.2.12)
∂ y 2 ! ∂ y2 3 ! ∂ y3
2 2 3 3
∂ u x h3 ∂ u x h3 ∂ u x
u3=u 0 +h3 + + +….........................................(4.2.13)
∂ x 2! ∂ x 2 3 ! ∂ x 3
2 2 3 3
∂ u y h4 ∂ u y h4 ∂ u y
u 4=u0 +h 4 + + +… .........................................(4.2.14)
∂ y 2! ∂ y 2 3! ∂ y 3

Nous avons considéré que u1, u0 et u3 ont la même coordonnée y, seules leurs coordonnées x
varient et u2, u0 et u4 possèdent la même coordonnée x, seules leurs coordonnées y varient.

La substitution de ces quatre équations dans l’équation (4.2.10) donne :

[ ∂2 u ∂ 2 u
2 ]
+ 2 =[ α 0 +α 1+ α 2 + α 3 +α 4 ] u0 + [ α 1 h1 −a3 h3 ]
∂x ∂ y o
∂u
∂x( ) 0

( ) ( )
2
∂u 1 2 ∂ u
+ [ α 2 h2−α 4 h4 ] + [ α 1 h1 +α 3 h3 ]
2
2
∂y 0 2 ∂x 0

( )
2 4
+1 2 ∂ u
2
[ α 2 h 2+ α 4 h 4 ]
2
2
∂ y 0 i=1
3
[ ]
+ ∑ O ( ai hi ) ......................................(4.2.15)

Pour que dans l’équation (4.2.15) le membre de droite soit égal au membre de gauche à
l’erreur de troncature O(h3), il faudrait que :

α 0 +α 1+ α 2 +α 3 +α 4=0 ,
α 1 h 1−a3 h3=0,
α 2 h 2−α 4 h4 =0,
2 2
α 1 h 1+ α 3 h3=2 et
2 2
α 2 h 2+ α 4 h 4=2 .

La solution de ce système d’équation est:

 a o=−2
[ 1
+
1
h1 h3 h2 h4
;
]
2
 a 1= ;
h1 ( h 1 + h3 )
2
 α 2= ;
h2 ( h2 +h4 )
2
 α 3= ;
h3 ( h1 +h3 )

29
 2
 α 4= ;
h 4 ( h2 +h 4 )

Lorsque ces valeurs de α sont substituées dans l’équation (4.2.10), on obtient les
approximations en différence finie recherchées. L’examen du dernier terme de l’équation
(4.2.15) révèle que, en général, l’erreur de troncature est O(h 3). L’approximation a une
erreur de troncature de O(h2 + k2) quand h2 = h4 = k et h1 = h3 =h. On obtient l’équation (4.2.9)
lorsque h1 = h2 = h3 = h4 = h.

Supposons maintenant que nous souhaitions utiliser les points ou nœuds illustrés
dans le diagramme discret de la figure (4.2.1). De plus, laissez tous les espacements égaux à
h. Puisque c'est une forme symétrique, nous pouvons écrire :

[ ∂2 u ∂ 2 u
+
∂ x2 ∂ y2 o ]
=α 0 u0 +α 1 [ u1 +u2 +u3 +u 4 ]..............................................(4.2.16)

Les coefficients α o et α 1 peuvent être déterminés au moyen des séries de Taylor suivantes :

[( ) ( ) ] [( ) ( ) ( ) ]
2 2 2 2
∂u ∂u h ∂ u ∂ u ∂u
ui=u o +h + + 2
+2 + +…...............(4.2.17)
∂x o ∂ y o 2! ∂x o ∂ x ∂ y o ∂ y2 o

Avec i = 1, 2, 3, 4

Sans présenter les détails des manipulations, nous affirmons que les équations suivantes
doivent être satisfaites :

α 0 +4 α 1=0 et 2 α 1 h2=1

Ainsi
−2 1
α 0= 2 et α 1= 2
h 2h

[ ]
2 2
∂u ∂ u 1
+ 2 = 2 [ u1 +u2 +u3 +u4 −4 u0 ] +O ( h )=0.................................(4.2.18)
2
2
∂ x ∂ y o 2h

Ceci est la représentation en différence finie pour l'équation de Laplace avec erreur
2
0(h ).

4.2.2.2. Equation de Poisson

Si, au lieu de l'équation de Laplace, nous avions considéré l'équation de Poisson,

l'approximation aux différences finies serait :

30
 1
u +u2 +u3 +u 4−4 u0 ]=2 f 0 +O ( h )=¿ ...............................(4.2.19)
2[ 1
2

L’ordre de l’erreur est O(h2)

L’équation (4.2.19) peut aussi s’écrire sous la forme vue ci-dessus :

1
2 [ r +1 , s+1
+ur −1 , s+1 +u r−1 ,s−1 +ur +1 , s−1−4 u r , s ]=2 f 0 +O ( h ) =¿............(4.2.19)
2
u
h

Il est évident que nous pourrions continuer cette méthode de coefficients

indéterminés pour produire de nombreuses approximations de différences finies différentes.
Bien que nous ayons utilisé la EDP elliptique comme exemple, la même approche peut
également être utilisée sur les EDPs hyperboliques et paraboliques.

5. INTRODUCTION A LA MODELISATION DES NAPPES AQUIFERES EN REGIME

TRANSITOIRE

Quel que soit le régime d’écoulement, transitoire ou permanent, uniforme ou

non uniforme, tout problème d’hydraulique générale ou d’hydraulique souterraine est
résolu si on connait, en chaque point, et le champ du potentiel hydraulique, ϕ, et le champ
de la fonction de courant, ψ .

En effet, la connaissance de φ sert par exemple :

 à déterminer le sens d’écoulement

 à déterminer la perte de charge entre deux points de l’écoulement,
 à calculer le gradient hydraulique, qui, multiplié par la conductivité hydraulique,
donne la vitesse de Darcy et le débit d’écoulement;
 La connaissant z, permet à déterminer la pression de l’eau en chaque point de
l’écoulement ;

ψ multiplié par la conductivité hydraulique donne la fonction de courant, Ψ qui s’exprime en

L2/T qui est le débit de l’aquifère sur une épaisseur unitaire. Pour trouver le débit total, on
multiplieΨ par l’épaisseur de l’aquifère.

Ce sont là les paramètres hydrauliques les plus recherchés dans la solution des problèmes
hydrauliques.

Les équations mathématiques à modéliser sont l’équation de la chaleur (2.2.14),

pour l’écoulement transitoire ou l’équation de Laplace ou l’équation de Poisson, pour
l’écoulement permanent.

31
5.1. DETERMINATION DE LA POSITION DE LA PIEZOMTRIQUE

5.1.1. Détermination de la position de la surface piézométrique d’un aquifère de grande

épaisseur à surface libre.

L’équation (2.2.9) de la surface libre que nous avons vue au sous-chapitre (2.2)
∂h
donne la vitesse instantanée, , de la montée ou de la descente de la surface
∂t
piézométrique, en chaque point de celle-ci d’un aquifère à surface libre :

(
∂h K ∂φ ∂h ∂h ∂φ ∂φ
= + −
∂t ε ∂x ∂x ∂ y ∂ y ∂z )
.............................................................................(2.2.9)

Sachant que, d’après le tableau (4.1.1) :

∂ u ur +1−ur −1
= O(h2) ................................................................(4.1.12)
∂x 2h

Nous pouvons écrire l’équation (2.2.9) comme suit :

2∆t
=
ε (
hr +1−hr −1 K φr +1−φ r−1 hr +1−hr −1 φr +1−φr−1 hr +1−h r−1 φ r+ 1−φr −1
2∆ x 2∆ x
+
2∆ y 2∆ y
−
2∆ z )
......................(4.3.1)

Comme nous l’avons dit lors de l’établissement de cette équation, cette vitesse,
comme son nom l’indique, est une vitesse instantanée. Elle change à chaque instant en
fonction de la position de la surface libre et des conditions qui règnent aux limites de
l’aquifère à cet instant-là. C’est ce qu’indique d’ailleurs l’équation (2.2.9) et son écriture en
différence finie qui permet de la calculer.

Comme on le voit, l’équation différentielle de la surface libre (équation (2.2.9)

n’est pas une équation parabolique, mais, comme sa solution numérique est similaire à celle
de l’équation de la chaleur qui est une équation parabolique et comme elle régit
l’écoulement transitoire comme l’équation de la chaleur, nous l’avons traitée dans ce sous-
chapitre relatif aux équations paraboliques.

Alors, pour connaître la position et la forme de la surface libre à une date future
donnée, on procède comme suit :

32
 a) On écrit l’équation (2.2.9) en différences finies en remplaçant ses dérivées par des
approximations des tableaux (4.1.1) et (4.1.2) selon que l’étude se fait en une ou
deux dimensions.
b) Au temps initial to, on détermine la position et la forme de la surface libre à l’aide des
conditions aux limites régnant au bord de l’aquifère à cet instant.
c) On calcule en chaque point de la surface libre la vitesse instantanée de son
mouvement à l’aide de l’équation (2.2.9) écrite en différences finies.
d) On découpe en petits temps, dt, le temps compris entre l’instant initial et la date
future à laquelle on voudrait connaître la position et la forme de la surface libre.
e) On multiplie la vitesse instantanée de chaque point par dt, ce qui permet de trouver
la position et la forme de la surface libre au temps to +dt.
f) On détermine les conditions aux limites de l’aquifère en cet instant t o +dt pour
connaitre le potentiel hydraulique, φ , en chaque point de la surface libre en ce temps
to + dt.
g) A l’aide du champ du potentiel hydraulique, on calcule la nouvelle vitesse
instantanée en chaque point de la surface libre à cet instant t o +dt à l’aide de
l’équation (2.2.9) écrite en différences finies.
h) On calcule la vitesse moyenne par la moyenne arithmétique entre celle trouvée à la
litera c et celle trouvée à la g.
i) On multiplie la vitesse moyenne trouvée à la litera h par dt pour trouver la nouvelle
position de la surface libre au temps t + dt.
j) On détermine les nouvelles conditions aux limites de l’aquifère en cet instant t o +dt
pour connaitre le potentiel hydraulique, φ , en chaque point de la surface libre en ce
temps to + dt
k) A l’aide de ce nouveau champ du potentiel hydraulique, on calcul la nouvelle vitesse
instantanée en chaque point de la surface libre à cet instant t o +dt à l’aide de
l’équation (2.2.9) écrite en différences finies.
l) On calcule la nouvelle vitesse moyenne par la moyenne arithmétique entre celle
trouvée à la litera c et celle trouvée à la k.
m) On compare cette nouvelle moyenne trouvée à la litera l avec celle trouvée à la litera
h. Si les deux moyennes ne diffèrent pas beaucoup, on retient la vitesse moyenne
calculé à la litera l comme la vitesse de la variation de h.
n) On calcule la nouvelle position de la surface piézométrique et on détermine les
conditions aux limites de l’aquifère régnant à cette position de la surface
piézométrique et on recommence le calcul comme si on était à la litera b
o) Si la comparaison faite à la litera m montre que les deux vitesses diffèrent trop, on
recommence le calcul à partir de la litera k.
p) On procède ainsi jusqu’à ce qu’on arrive à la date choisie pour déterminer la position
et la forme de la surface libre.

33
5.1.2. Détermination de la position de la surface piézométrique d’un aquifère de faible
épaisseur à surface libre

5.2. DETERMINATION DU POTENTIEL HYDRAULIQUE DANS L’AQUIFERE

5.2.1. Modélisation en différences finies de l’équation de la chaleur (équation

parabolique), ci-dessous :
2
∂u ∂ u
= ………………………………………………………………………………………………………..(4.3.2)
∂ x ∂ y2

5.2.1.1. Les Approximations

Comme nous avons différents modèles des EDP paraboliques, nous avons
également un ensemble d'approximations de ces différents modèles des EDP. Ce sera donc
ces équations aux différences finies (EDF) que nous devons résoudre. Les EDF doivent être
bien posés elles-mêmes ainsi que cohérentes, convergentes et stables (voir les sections 4.4
et 4.5).
Pour obtenir ces représentations de différences finies, nous utilisons la notation
du chapitre 1, à savoir que pour un problème bidimensionnel u (x, y), les incréments en x
seront de longueur h et ceux en y seront de longueur k. S'il y a plus d'une variable spatiale,
nous utilisons l1, l2, l3... Les indices r et s (c'est-à-dire ur,s) dénotent des valeurs de u à des
valeurs incrémentales de x (r) et y (s), avec d'autres indices utilisé pour plus de directions
spatiales. Lorsque h et k sont des valeurs constantes, nous nous référons au maillage
bidimensionnel qui en résulte (voir figure 4.3.1) comme un maillage uniforme. C'est
évidemment la procédure la plus simple que nous puissions analyser. Maintenant, nous
avons l'habitude de chercher à calculer u aux points du maillage, u r,s, où se produisent les
intersections de maillage.

34
 Figure 4.3.1. Choix du maillage et des nœuds de calcul.
Cela dépendant de conditions aux limites, cependant, il peut être avantageux
d'utiliser une grille centrée sur les blocs dans laquelle u est calculé au centre d'un bloc,
ur+1/2,s+1/2. Ce dernier cas est particulièrement utile lorsque le type d’écoulement et les
conditions aux limites sont spécifiées. Plus tard, nous étudierons également les
caractéristiques des maillages non uniformes dans lesquels h et k seront respectivement des
fonctions de r et s. Cela sera particulièrement utile lorsque les gradients sont très différents
dans diverses régions du domaine (x, y). En général, un maillage non uniforme dans la seule
direction x peut économiser du temps de calcul, surtout si la solution EDP est
fondamentalement amortie en fonction exponentielle de x. Dans ce cas, une augmentation
exponentielle de h au fur et à mesure que l'on avance dans l'espace x permet de conserver la
précision tout en minimisant le temps de calcul requis.

La forme linéaire générale d'une approximation par différence finie peut être
écrite sous forme matricielle comme suit :

[ A ] o { u }r +1=[ A ] 1 {u }r + [ A ] 2 {u }r −1+ [ A ]3 { u }r −2 … … . [ A ] m {u }r−m+ 1................(4.3.3)

Où

 [ A ] o, [ A ] 1 , [ A ]2 , … , [ A ]m sont des matrices carrés

 { u }r +1 , {u }r , {u }r−1 , {u }r −2 , … , { u }r −m+1 peuvent être des vecteurs appropriés

L'indice s a été supprimé avec les valeurs incorporées dans les vecteurs {u} r.
L'équation (4.3.2) dit que nous pouvons calculer {u} r+1 [au (r + l)ème niveau] en incorporant
les m niveaux de x en dessous. Bien sûr, [ A ] o ne peut pas être une matrice singulière. Dans la
plupart des cas, m = 1 ou m = 2, car des valeurs plus élevées de m nécessitent un temps de
calcul important. Si [ A ] o est une matrice diagonale, l’équation (4.3.3) est appelé une forme
de différence explicite. Dans tous les autres cas l’équation (4.3.3) est appelée une forme de
différence implicite. Les sections ultérieures de ce chapitre clarifient cette distinction et
montrent les différences de calcul.

Autrement, on peut adopter la convention selon laquelle la forme générale de

différence finie linéaire ou non linéaire soit donnée comme :

35
 ur +1= {u }r +h ( P1 {u }r+ 1+ P 0 { u }r ).................................................(4.3.4)

Où P1 et Po sont des opérateurs composés de combinaisons polynomiales d'opérateurs de

différences définies sur la région spatiale et agissante sur leurs arguments de manière
linéaire ou non linéaire. Cela correspond évidemment au sens large au cas m = l dans (4.3.3)
et est implicite. De nombreux exemples de (4.3.4) sont donnés plus loin.

Il est important de réaliser que lorsque notre modèle EDP est linéaire,
l'approximation par différence finie du modèle résultant est également linéaire. Cependant,
Gerber et Miranker (1973) ont proposé qu'un schéma de différence non linéaire puisse avoir
des avantages. A chaque étape x, le calcul se poursuit en déterminant le meilleur schéma de
différence possible à utiliser à cette étape sous réserve d'un critère d'optimisation. Cette
approche a une grande précision tout en étant stable et convergente; il se pourrait fort bien
que cette procédure non classique soit d'une importance capitale à mesure que de nouveaux
développements se produiront.

Nous avons souligné que notre modèle EDP doit être bien posé pour diverses
données initiales, mais il s'ensuit également que l'approximation des différences finies doit
être bien posée. En général, une approximation des différences finies bien posée doit être
continue pour une classe plus large de perturbations de données initiales que l’EDP ; cela est
directement lié à la stabilité de l'approximation des différences finies dans laquelle toutes les
perturbations dans une solution calculée doivent être limitées. Lorsque les perturbations
dans les données initiales sont réduites (dans un certain sens, comme h, k -> 0), les
perturbations résultantes dans la solution calculée doivent devenir faibles, entraînant ainsi la
convergence de l'approximation des différences finies. Nous étudierons certaines de ces
fonctionnalités plus tard.

5.2.1.2. Dérivation des approximations des différences finies

Nous avons montré à la section 2.1.2 que l'on peut utiliser des extensions de
séries de Taylor tronquées pour remplacer un EDP par une approximation aux différences
finies. Les termes tronqués conduisent à l'erreur de troncature locale, qui est utilisée, au
moins dans un sens local, pour fournir une mesure de la précision de l'approximation des
différences finies. Pour une EDP, il existe de nombreuses représentations de différences
finies possibles et nous pouvons penser que des schémas avec des erreurs de troncature
locales plus petites (c'est-à-dire un ordre supérieur en puissances de h, k, ...) doivent être
préférés.
Il arrive fréquemment, mais pas toujours, qu’un cas puisse présenter d'autres
caractéristiques, en plus de l'erreur de troncature, qui font qu'une approximation de
différence finie particulière soit un candidat préférable pour un calcul ou une simulation.

36
Dans cette section, nous prenons le modèle EDP parabolique sous la forme de l’équation
(4.3.2).
2
∂u ∂ u
= ………………………………………………………………………………………………………..(4.3.2)
∂ x ∂ y2

En supposant que u(x, y) de l’équation (4.3.2) soit continue et est suffisamment de dérivable,
il est simple de montrer par différenciation répétée (en supposant que l'ordre de
différenciation peut être interchangé), que nous pouvons écrire :

2 4 3 6 mx 2 my
∂ u ∂ u ∂ u ∂ u ∂ u ∂ u
2
= 4 , 3 = 6 , … mx = 2 my ……………………………………………………………(4.3.5)
∂x ∂ y ∂x ∂ y ∂x ∂y

Maintenant, si nous développons une série de Taylor autour de ur,s, il en résulte

2 2 3 3
∂u h ∂ u h ∂ u
ur +1 , s=ur ,s +h + + +…
∂ x 2! ∂ x 2 6 ! ∂ x 3

Si nous posons h=ρ k 2, nous pouvons écrire :

2 2
∂ u ( ρ k ) ∂2 u ( ρk ) ∂ 3 u
3
+… ………………………………………………(3.3.6)
2
ur +1 , s=ur ,s + ρ k + +
∂x 2! ∂ x 2 6 ! ∂ x 3
ou

2
ρ
ur +1 , s=ur ,s + ρB+ D+…
2!

Avec
2
2 ∂u 4 ∂ u
B=k , D=k
∂x ∂x
2

De la même manière, nous pouvons dériver tous les éléments suivants (nous supprimons les
indices r, s sur les côtés droit de l’équation) :

1 2 1 3 1 4 1 5
ur ±1 ,s =u ± ρB+ ρ D ± ρ F + ρ H ± ρ J +…………………………………..(4.3.7)
2 6 24 120
1 1 1 1
ur ,s +1+ ur , s−1=2 u+B+ D+ F+ H+ J +… …………………..(4.3.8)
12 360 20160 1814400
4 8 4 8
ur ,s +2+ ur , s−2=2 u+4 B+ D+ F+ H+ J + ……………………………..(4.3.9)
3 45 315 14175
ur ±1 ,s +1 +ur ± 1, s−1=2 u+(1± 2 ρ)B+
1
12 (± ρ+ ρ2 D+
1
) (±
360 12
1 1
ρ+ ρ2 F +
3
1
± ) (
20160 360
1 1 1 1
ρ+ ρ2 ± ρ3+ ρ4
24 6 12 )
………….(4.3.10)
37
 2
∂u ∂ u
Incidemment, si l’équation (4.3.2) avait été écrite comme =α 2 avec α >0, les mêmes
∂x ∂y
résultats se seraient produits avec le seul changement que αh=ρ k .
2

En utilisant directement les séries de Taylor ou en combinant les développements ci-dessus,

nous pouvons développer des approximations aux différences finies pour le modèle PDE
parabolique qu’est l’équation (4.3.2). En outre, l'erreur de troncature locale est
immédiatement disponible. Ensuite, nous calculons certaines des combinaisons possibles.

5.2.1.2.1. L'approximation explicite classique

Multiplions (4.3.8) par ρ et soustrayons le résultat de (4.3.7) (avec l'indice plus). Cela donne :

( )
1 1 1
ur +1 , s=( 1−2 ρ ) ur , s + ρ ( ur ,s +1+ ur , s−1 ) + ρ ρ− D+ ρ ρ2−
2 6 6
1
60 ( )
F +… …….(4.3.11)

Ainsi, si nous utilisons cette forme de différence finie, appelée forme explicite classique,

ur +1 , s=( 1−2 ρ ) ur , s + ρ ( ur ,s +1+ ur , s−1 )……………………………………………………………………(4.3.12)

l'erreur de troncature locale aura pour partie majeure ou principale

1
2 ( )
1
ρ ρ− D. Mais
6
h
puisque ρ= 2 , il s'ensuit que
k

( ) ( )
2 4
1 1 1 2 hk ∂ u
ρ ρ− D= h − ……………………………………………………………………………………(4.3.13)
2 6 2 12 ∂ y 4

C'est essentiellement le résultat obtenu en (4.2.5) et on dit que (4.3.12) est 0(h 2 + hk2) ou 0(h
+ k2) ou dans une notation plus simple 0(1,2).

ur +1 , s=ρ ur , s+ 1+u r , s ( 1−2 ρ )+ ρur , s−1............................................(4.2.5)

Le lecteur peut voir que la méthode explicite classique relie des points de maillage sur la
grille (r, s) qui a la forme :

Figure 4.3.2. Représentation du nœud de calcul de l’équation explicite classique

Cette molécule de calcul montre clairement que ur +1 , s, la sphère ombrée, peut être calculée
à partir des trois points inférieursur ,s−1 ,u r , s et ur ,s +1. Nous suivons la convention d'attacher un

38
tel diagramme symbolique à côté de chaque nouvelle approximation aux différences finies
développée.

1
Notez que si nous choisissons le cas particulier où ρ= = dans (4.3.11), alors (4.3.12) tient
6
toujours mais l'erreur de troncature locale a une partie principale qui est

( )
6 6
1 2 1 1 1 3∂u 1 6 ∂ u
ρ ρ− F= F= h − hk
6 60 3240 6 ∂ y 6 360 ∂y
6

1
Ainsi, si ρ= , l'approximation explicite classique du modèle PDE parabolique est 0(h3 + hk4)
6
2 4
ou 0(h + k ) ou 0(2,4).

Tous ces résultats auraient pu être dérivés directement des extensions de séries de Taylor
autour de ur,s. Nous illustrons la procédure pour le cas ou (4.3.12) comme

ur +1 , s−( 1−2 ρ ) ur , s− ρ ( u r , s+1 +ur ,s +1 )

[ ] [ ]
2 2 3 3 2 2 3 3
∂u h ∂ u h ∂ u ∂u k ∂ u k ∂ u
¿ u r , s+ h + 2
+ 3
+ … −( 1−2 ρ ) u r , s−ρ ur ,s +k + +
∂ x 2 ∂ x 6 ∂u ∂ y 6 ∂ y2 6 ∂ y3

[
−ρ ur ,s −k
∂u k 2 ∂2 u k 3 ∂3 u
+ −
∂ y 6 ∂ y2 6 ∂ y3
+…
]
En rassemblant les termes sur le côté droit, nous constatons que :

( )
2 2 4
∂u ∂ u 1 2 ∂ u 1 2 ∂ u
ur +1 , s−( 1−2 ρ ) ur , s− ρ ( u r , s+1 +ur ,s +1 )=h − 2 + h 2
− hk 4
+ ...
∂x ∂ y 2 ∂ x 12 ∂y
2 2 4
∂u ∂ u ∂ u ∂ u
Parce que − 2 =0, nous voyons que pour 2
= 4 , nous avons une erreur de
∂x ∂ y ∂x ∂ y
troncature locale identique à celui obtenu pour (4.3.12). Nous aurions pu aborder cela de la
manière suivante :

∂u ∂ u
2
ur +1 , s−ur ,s ur , s +1−2u r , s+ ur , s−1
− 2 =0 ≈ − 2 ………………………………………………………(4.3.14)
∂x ∂ y h k

avec la première approximation des différences finies à droite étant O(h) et la seconde 0(k 2).
Une simple manipulation du côté droit donne directement (4.3.12). Pour référence future,
nous pourrions utiliser la notation :

( δ 2 u )r , s=ur ,s +1−2 ur , s +ur ,s−1……………………………………………………………………(4.3.15 a)

( δ 2 u )r +1 ,s=ur +1, s +1−2u r+ 1, s +ur +1 , s−1…………………………………….……………….(4.3.15 b)

et écrire (4.3.14) comme

39
 2
∂u ∂ u ur +1 , s−ur ,s ( δ 2 u )r , s
− =0 ≈ −
∂ x ∂ y2 h k2

5.2.1.2.2. L'Approximation Explicite DuFort-Frankel

Si nous multiplions l’équation (4.3.6) par (1+2 ρ), l’équation (4.3.8)

par (−2 ρ) et l’équation (4.3.7) avec ses indices par −(1−2 ρ) et puis nous additionnons les
trois résultats, nous aurons :

( 1+2 ρ ) ur +1, s=2 ρ ( ur ,s +1+ ur , s−1 ) +(1−2 ρ)ur −1 , s…………………………..(4.3.16)

Figure 4.3.3. Représentation du nœud de calcul de l’équation de DuFort-Frankel

L’équation (4.3.18) est appelé l’approximation explicite de DuFort-Frankel et possède une

erreur de troncature :

( ) ( )
3
ρ ρ
3 ρ
2ρ −
D+ − F
6 3 180
En replaçant ρ , D et F comme ci-dessus, nous trouvons que l’approximation par différence

( )
3
3 2 h
finie de DuFort-Frankel possède une erreur de troncature locale de O h + h k + 2 ou
k

( ) ( )
2
2 h2 2 h
O h + k + 2 ou O 2 , 2 , . Notez que cela vaut pour le rapport et pour les habituels
k 2 k
termes indépendants h et k.

Sans le démontrer ici, nous disons que si nous avions utilisé le développement de Taylor
directement, nous aurions trouvé une forme équivalente de l’erreur de troncature.

L’équation (4.3.16) peut aussi se développer sous l’approximation suivante :

40
 2
∂u ∂ u ur +1 , s−ur −1 , s u r , s+1−( ur +1 ,s +ur −1 ,s ) +u r , s−1
− 2 =0= − 2
∂x ∂ y 2h k

( )
2
δ u
ce qui n'est pas directement évident. Ce qui a été fait, c'est de remplacer u r,s dans 2
k r,s
2
∂u
(une approximation de 2 ) par la moyenne arithmétique de (u r+1,s + ur-1,s). Alternativement,
∂y
2
∂u
on peut dériver l'approximation de DuFort-Frankel en remplaçant 2 par une
∂y
approximation précise du second ordre centrée en x et en y.
Donc

−
∂ x ∂ y2
=
2h
−
[ k2
−
h2 ( )]
∂ u ∂2 u u r+1 , s−ur −1 ,s ur ,s +1−2 ur , s +ur , s−1 u r+1 , s−2 ur ,s +ur −1 ,s h2
k2
………...……(4.3.17)

2
h
Sous cette forme, on voit plus facilement pourquoi l'erreur de troncature a un terme 2 .
k
2
Aussi, il faut noter que l'approximation DuFort-Frankel a une erreur plus petite (h ) que
1
l'approximation explicite classique avec ρ ≠ .
6

5.2.1.2.3. L'approximation explicite de Richardson

Une autre approximation explicite aux différences finies qui est plus précise que l'explicite
classique est due à Richardson et est appelée approximation saute-mouton. Il peut être écrit
comme suit :
ur +1 , s−ur −1 , s−2 ρ [ ur , s +1+u r , s−1 ] + 4 ρ ur ,s=0………………………………………………….(4.3.18)

Figure 4.3.3. Représentation du nœud de calcul de l’équation de Richardson

41
et résulte d'une substitution de la forme :

2
∂u ∂ u ur +1 , s−ur −1 , s u r , s+1−2 r , s+u r , s−1
− 2 =0= − 2 …………………………………………….(4.3.19)
∂x ∂ y 2h k

L'équation (4.3.18) est 0(h2 + k2) ou 0(2,2). Nous ne développons pas l'erreur de troncature
locale explicite car, il est inconditionnellement instable et ne doit être utilisé dans aucun
calcul.

Il est cependant important de noter que les approximations de DuFort-Frankel et de

Richardson ne diffèrent que par la manière dont u r,s est traité dans la représentation aux
2
∂u
différences finies pour 2 . Cette différence s'avère suffisante pour rendre le premier
∂y
modèle stable et le second instable ; il démontre également qu'il faut faire très attention
pour obtenir une approximation valide (calculable) des différences finies.

5.2.1.3. COHÉRENCE ET CONVERGENCE

Dans cette section, nous analysons les concepts de cohérence et de convergence des
approximations aux différences finies et dans la section suivante celle de stabilité. Nous
supposons dans ce développement que le modèle PDE parabolique est bien posé.

En termes généraux, nous pouvons dire que la stabilité (de calcul) appelle la limitation de
toutes les perturbations dans une solution calculée. En d'autres termes, lorsque les
magnitudes des perturbations dans les données initiales sont rendues arbitrairement petites
comme h , k , … → 0 , les perturbations résultantes dans la solution calculée doivent
disparaître plutôt que croître. Dans un tel cas, les solutions calculées basées sur une
approximation de différence cohérente convergent vers la solution du modèle PDE
parabolique. Autrement dit, la stabilité et la cohérence impliquent une convergence telle
que définie par Lax. En d'autres termes, le succès de l'obtention d'une solution numérique
convergente basée sur une approximation particulière des différences finies réside dans la
cohérence de la forme des différences finies avec le modèle bien posé PDE et la stabilité de
la forme des différences finies. En conséquence, ces concepts sont d'une importance
majeure dans l'évaluation de la viabilité de toute approximation aux différences finies.

La cohérence a à voir avec l'assurance que, lorsque le maillage aux différences finies est
raffiné (c'est-à-dire lorsqueh , k , … → 0), les erreurs de troncature tendent vers zéro.
Alternativement, on peut dire que le modèle aux différences finies se rapproche du PDE
souhaité et non d'un autre PDE. À titre d'exemple, considérons la forme classique des
différences finies explicites. En utilisant (4.3.11), nous pouvons développer la relation
approximative

42
5.2.1.3.1. KL

5.2.1.3.2. KL

5.2.1.4. Jk
5.2.1.5. KL

6. INTRODUCTION A LA MODELISATION DES NAPPES AQUIFERES EN REGIME

PERMANENT.

Comme nous l’avons vu au sous-chapitre 1.3, l’équation qui régit le potentiel

hydraulique en chaque point de la nappe aquifère en régime permanent est l’équation de
Laplace :
2 2 2
∂ φ ∂φ ∂ φ
2
+ 2 + 2 = Δ φ=0...............................................................(4.4.1)
∂x ∂ y ∂z

Lorsque l’écoulement est permanent avec alimentation en eau ou soustraction

d’eau en certains points de l’aquifères, alors l’équation qui régit le potentiel hydraulique en
de tels points est l’équation de Poisson :
2 2 2
∂ φ ∂φ ∂ φ
2
+ 2 + 2 =f (x , y , z )..................................................................(4.4.2)
∂x ∂ y ∂z

6.1.1. Modélisation de l’équation de Laplace

Dans cette section, nous n’allons considérer que les écoulements permanents
monodimensionnels et bidimensionnels dans une nappe aquifère homogène et isotrope.

6.1.1.1. Modélisation des écoulements monodimensionnels

Pour cet écoulement, le laplacien est

2
∂ φ
2
=Δ φ=0 .....................................................................(4.4.3)
∂x

43
 Soit, φ , le potentiel au nœud r, le potentiel, φ r+1 , au nœud r+1,
suivant sur l’axe des x et séparé du noeudr par une distance, h, est donné par plusieurs
approximations telles que celles que nous avons présentées au tableau (4.1.1). L’une de ces
approximations, la plus simple d’entre elles dont l’erreur de troncature est :

∂ u ur +1−2 ur +ur −1
2

2
= 2 O(h) ...........................................................(4.1.15)
∂x h

Nous pouvons écrire cette équation en remplaçant u par φ

∂ φ φr +1−2 φr + φr −1
2

2
= 2 ........................................................................(4.4.4)
∂x h
2
∂ φ
Comme 2
=0 (Voir équation (4.4.3)
∂x

Alors

φ r+1 −2 φr +φ r−1=0
Et

φ r+1 +φ r−1
φ r=
2

φ r est donc déterminé par la moyenne arithmétique des potentiels hydrauliques aux nœuds
voisins avec comme erreur, l’erreur de troncature qui est ici, d’ordre 1

20.1.1.2. Modélisation des écoulements bidimensionnels

44
Le potentiel hydraulique est connu à chaque nœud écrit en rouge qui sont sur la limite. On
voudrais déterminer les potentiels hydrauliques inconnus aux nœuds intérieurs (6, 7, 10 et
11). Or la série de Taylor donne pour chaque nœud, i, l’expression suivante pour son
potentiel hydraulique φ i :

2 3 4
1 ∂φ 1 ∂ φ 2 1 ∂ φ 3 1 ∂ φ 4
φ i+1=φ i+ h+ h+ h+ h +…
1! ∂x 2! ∂ x 2 3 ! ∂ x3 4 ! ∂ x4

2 3 4
1 ∂φ 1 ∂φ 2 1 ∂ φ 3 1 ∂ φ 4
φ i−1 =φi− h+ h− h+ h +…
1! ∂x 2! ∂ x 2 3 ! ∂ x3 4 ! ∂ x4

2 3 4
1 ∂φ 1 ∂ φ 2 1 ∂ φ 3 1 ∂ φ 4
φ j+ 1=φi + k+ k + k + k +…
1! ∂ y 2! ∂ y 2 3! ∂ y 3 4 ! ∂ y4

2 3 4
1 ∂φ 1 ∂ φ 2 1 ∂ φ 3 1 ∂ φ 4
φ j−1=φ i− k+ k− k + k +…
1! ∂ y 2! ∂ y 2 3 ! ∂ y3 4 ! ∂ y4

L’addition de ces quatre équations donne :

2 2 4 4
∂ φ 2 ∂ φ 2 1∂ φ 4 1∂ φ 4
φ i+1 +φi−1 +φ j+ 1+ φ j−1=4 φ i+ 2
h + 2k + h+ k +…
∂x ∂y 2 ∂ x4 2 ∂ y4

Si tous les termes dont l’ordre de la dérivée est égal ou supérieur à 4 sont négligés et si on
prend h = k, on aura

φi +1+ φi−1+ φ j+1 +φ j−1

4 ∂ φ ∂φ
2 2

2 2 i
=
φ + 2+ 2
h h ∂x ∂ y
Or le laplacien du potentiel hydraulique d’un écoulement permanent est nul en chaque point
2 2
∂ φ ∂φ
de l’écoulement ( 2 + 2 =0) ; donc
∂x ∂ y

45
 φi +1+ φi−1+ φ j+1 +φ j−1 4
2
= φ
2 i
h h
Ou
φi+1 + φi−1+ φ j +1 +φ j−1
φ i=
4
Comme pour l’écoulement monodimensionnel, le potentiel hydraulique en un nœud est la
moyenne arithmétique des potentiels aux nœuds voisins ; à l’erreur de troncature près.

20.1.1.3. Modélisation des écoulements tridimensionnels

Le même raisonnement conduit à

φi+1 + φi−1+ φ j +1 +φ j−1+ φl+ 1+ φl−1

φ i=
6

6.1.2. Méthodes de résolution

La solution consiste à trouver les potentiels hydrauliques aux nœuds intérieurs

où ces paramètres sont inconnus.

Nous allons montrer ces méthodes seulement pour l’écoulement bidimensionnel

le plus utilisé en hydrogéologie et aussi du fait que le même raisonnement permet de
résoudre les problèmes des écoulements monodimensionnels et tridimensionnels.

On distingue les méthodes directes et les méthodes indirectes. Les méthodes

directes trouvent la solution exacte tandis que les méthodes indirectes trouvent une solution
approchée de la solution exacte.

6.1.2.1. Méthodes directes

Par les méthodes directes on écrit un système d’équations linéaires et on résout

ce système par une des méthodes, notamment par la méthode d’inversion des matrices que
nous allons voir.

On a vu que les potentiels hydrauliques inconnus des nœuds intérieurs peuvent s’écrire par
les équations suivantes :

46
 φ5 φ7 φ 2 φ10
φ 6= + + +
4 4 4 4

φ6 φ8 φ 3 φ11
φ 7= + + +
4 4 4 4

φ9 φ 11 φ 6 φ14
φ 10= + + +
4 4 4 4

φ10 φ12 φ7 φ15

φ 11= + + +
4 4 4 4

Ces 4 équations peuvent s’écrire sous la forme suivante

Ce qui peut encore s’écrire en mettant les valeurs connues au second membre des
équations :

Ce système d’équation peut se simplifier en écrivant :

Les potentiels hydrauliques inconnus peuvent s’écrire sous la forme d’une matrice colonne,
x:

[]
φ6
φ
x= 7
φ10
φ11

Et les potentiels hydrauliques connus peuvent s’écrire en une matrice colonne, b :

47
 [ ]
φ2 + φ5
φ +φ
b= 3 8
φ9 +φ14
φ 12+ φ15

Les coefficients quant à eux, peuvent s’écrire sous la forme d’une matrice carrée, A

[ ]
4 −1 −1 0
−1 4 0 −1
A=
−1 0 4 −1
0 −1 −1 4

Ces matrices sont liées par la relation suivante :

[ A ] { x }={ b }
La matrice colonne,{ x } , peut être trouvée par la multiplication de la matrice inverse,
[ A ]−1 , de la matrice [ A ] par la matrice colonne, { b } :
−1
{ x }= [ A ] { b }

L’inversion des matrices n’est facile que si l’équation a un nombre d’inconnu inférieur à 100,
or il arrive que les nœuds où le potentiel hydraulique à rechercher sont en nombre
dépassant plusieurs milliers, voire, plusieurs dizaines des milliers. Alors on a recours aux
méthodes indirectes.

6.1.2.2. Méthodes indirectes.

La méthode indirecte la plus connue est la méthode itérative. Dans

notre cas, sachant que la solution de l’équation donnant la matrice colonne est unique, alors
nous pouvons imaginer une itération qui nous permet de comparer les solutions trouvées à
chaque tour et arrêter le calcul au moment où la solution ne change plus d’un tour de calcul
au tour suivant. Chaque tour de calcul se fait de la façon suivante :

a) Au premier tour, l’ordinateur calcule φ i de chaque nœud par la moyenne

arithmétique des potentiels hydrauliques des quatre nœuds voisins comme l’indique
l’équation ci-dessus. Comme aux nœuds intérieurs le potentiel hydraulique est
inconnu, l’ordinateur considère, pour ce premier tour de calcul, que leur potentiel est
nul.

48
b) Après ce premier tour l’ordinateur vérifie pour chaque nœud, i, si la valeur absolue
de la différence entre le potentiel hydraulique initial et le potentiel qu’il vient d’y
calculer est nulle.

c) Si la vérification faite à la litera b) montre que cette valeur absolue de la différence

est nulle à chaque nœud, le calcul est terminé. L’ordinateur peut, selon le désire de
l’opérateur ou la conception du logiciel, soit afficher les potentiels hydrauliques de
chaque nœud soit tracer des courbes équipotentielles de l’aquifère.

d) Si la vérification faite à la litera b) montre que la valeur absolue de la différence n’est

pas nulle à chaque nœud, l’ordinateur reprend le calcul du potentiel hydraulique de
chaque nœud par la moyenne arithmétique de ses quatre nœuds voisins.

e) Si la valeur absolue de la différence entre le potentiel du potentiel hydraulique

calculé au tour précédent et celui qu’il vient d’y calculer est nulle à chaque nœud, le
calcul est terminé et l’ordinateur peut procéder comme à la litera c) ci-dessus.

f) Si cette différence n’est pas nulle à chacun des points du maillage, le calcul de
potentiel hydraulique et la vérification des valeurs absolues des différences des
potentiels hydrauliques reprennent.

g) L’ordinateur peut faire autant de tours de calcul et de vérification jusqu’à ce que la

valeur absolue de la différence à chaque nœud entre le potentiel calculé et le
potentiel y calculé au tour précédent soit nulle.

h) En pratique, l’opérateur se donne un degré de précision en choisissant un nombre

très proche de zéro, par exemple 0,001. Il demande, à chaque tour de calcul, à
l’ordinateur de chercher sur tous les nœuds, le nœud où la valeur absolue de la
différence entre le potentiel hydraulique y calculé et celui y calculé au tour précédent
est la plus grande de tous les nœuds. Il demande ensuite, à l’ordinateur, de vérifier si
cette différence absolue la plus grande est inférieure ou égale à la précision qu’il
s’était choisi (ici la précision est 0,001). Après cette vérification, si la valeur absolue
de la plus grande différence est plus grande que la précision, l’ordinateur reprend le
calcul du potentiel hydraulique à chaque nœud. Si cette valeur maximale est égale ou
inférieure à la précision, le calcul est terminé et l’ordinateur procède comme à la
litera c) ci-dessus.

49
6.1.3. Modélisation de l’équation de Poisson

Il peut arriver qu’à un nœud il y ait alimentation. L’alimentation est positive si le

nœud reçoit de l’eau (présence au nœud d’un puits injectant de l’eau, ou infiltration d’eau
de pluie ou d’irrigation, par exemple) et est négative si le nœud perd de l’eau (présence au
nœud d’un puits pompant de l’eau, par exemple).

L’équation d’écoulement qui règne à un tel nœud n’est plus celle de Laplace mais
plutôt celle de Poisson :

∆ φ=f (x , y , z)

Prenons un nœud situé au centre d’un cube formé par le bloc de son influence ;
c’est-à-dire le cube obtenu par des plans parallèles aux axes x, y et z et passant par des
milieux de segments joignant le nœud aux nœuds voisins. Nous avons démontré, au sous-
chapitre 4.2., qu’à l’absence du puits (injectant ou pompant de l’eau), le débit qui entre dans
un tel cube est égal au début qui en sort. De cette égalité et de l’hypothèse
d’incompressibilité de l’eau, nous sommes arrivés à la conclusion que la divergence du
vecteur vitesse de l’écoulement de l’eau et nulle. En remplaçant la vitesse par ses
composantes dans la loi de Darcy généralisée, nous avons établi l’équation de Laplace.

Maintenant nous allons montrer qu’avec un puits pompant la quantité de l’eau

qui entre dans ce cube est égale à celle qui en sort par le pompage. Soit Q, le débit du
pompage. Nous allons donc écrire que la quantité d’eau qui sort est égale au débit du
pompage :

( V x+
∂V x
∂x ) (
dx dydz + V x +
∂V y
∂y )
dx dxdz+ V x + (
∂V z
∂z )
dx dxdy=Q

∂ V x ∂V y V z
V x dydz +V y dxdz +V z dxdy +( + + )dxdydz =Q
∂x ∂ y ∂ z

La divergence du vecteur vitesse étant nulle, on a :

V x dydz ++V y dxdz +V z dxdy=Q

Ce qui revient à dire que le débit qui entre dans le cube d’influence du nœud est
égale au débit qui en sort par le pompage.

Si, pour simplifier, nous posons dx = dy = dz = h, la distance séparant le nœud

d’avec ses nœuds voisins, nous pouvons écrire cette équation comme :

Q
V x +V y +V z = 2
h

50
 Donc, d’après la loi de Darcy

∂φ ∂φ ∂φ Q
K +K +K =
∂x ∂y ∂ z h2
Ou
∂φ ∂φ ∂φ Q
+ + =
∂ x ∂ y ∂ z K h2

Les dérivées premières de φ , peuvent s’écrire en différence finie, de plusieurs façons.

Prenons en deux :

∂ φ φ(x , y , z )−φ(x−1 , y, z )
a) =
∂x h
∂ φ 4 φ(x+1 , y , z )−φ(x +2 , y , z) −3 φ(x , y , z)
b) =
∂x 2h

La première approximation a pour erreur de troncature, d’ordre 1, alors que l’erreur de

troncature de la seconde est d’ordre 2. Prenons donc cette dernière.

Alors, nous écrivons :

4 φ(x+1 , y , z )−φ(x+ 2, y , z)−3 φ(x , y , z) 4 φ( x, y+1 , z) −φ(x , y +2 , z)−3 φ( x, y , z) 4 φ(x , y ,z +1)−φ(x , y , z +2)−3 φ(x , y , z) Q
+ + =
2h 2h 2h K h2

Ou

2Q
−9 φ( x, y , z) +4 φ(x+1 , y , z)−φ( x+2 , y , z )+ 4 φ( x, y+1 , z) −φ(x , y+2 , z) + 4 φ(x , y , z +1)−φ(x , y , z+ 2)=
Kh

1 2Q
φ(x , y , z)= [4 φ( x+1 , y, z )−φ ( x+2 , y , z )+ 4 φ ( x, y+1 , z )−φ( x , y+2 , z )+ 4 φ( x , y ,z +1 )−φ( x , y , z +2) + ]
9 Kh

Ou
φ(x , y , z)−
1
9 [
4 φ( x+1 , y, z )−φ ( x+2 , y , z )+ 4 φ ( x, y+1 , z )−φ (x , y+2 , z )+ 4 φ (x , y , z+1 )−φ( x , y , z +2) +
]
2Q
Kh
=0

C’est l’approximation, en écoulement tridimensionnel, du potentiel hydraulique à un nœud

où sort ou entre un débit Q. Si le débit entre, Q prend un signe positif, si le débit sort, Q est
affecté d’un signe négatif.

Pour l’écoulement bidimensionnel, l’approximation s’écrit :

1 2Q
φ(x , y)= [4 φ( x+1 , y )−φ (x +2 , y )+ 4 φ (x , y+1)−φ( x , y +2) + ]
6 Kh

Ou

51
 φ(x , y)−
1
6 [
4 φ( x+1 , y )−φ ( x+2 , y )+ 4 φ (x , y+1 )−φ( x , y +2) +
]2Q
Kh
=0

Pour l’écoulement monodimensionnel, on aura :

1 2Q
φ(x , y)= [4 φ ( x+1)−φ( x+2) + ]
3 Kh

Ou

φ(x , y)−
1
3[4 φ ( x+1)−φ( x+2 )+
2Q
Kh ]
=0

Si le calcul à chaque nœud tel qu’indiqué pour l’équation elliptique de Laplace donne zéro, le
calcul est terminé et l’ordinateur procède comme à la litera c) du sous-chapitre 16.3. En
pratique, comme au sous-chapitre 16.3. ; on se donne un degré de précision et l’on compare
les résultats à ce degré de précision.

7. INTRODUCTION A LA MODELISATION PAR LA METHODE DES ELEMENTS FINIS

7.1. INTRODUCTION

Contrairement aux schémas de différences finies dans lesquels le domaine

d'intérêt est remplacé par un ensemble de points discrets, dans la méthode des éléments
finis, le domaine est divisé en sous-domaines communément appelés « éléments finis ». La
fonction inconnue, disons-nous, est représentée dans chaque élément par un polynôme
d'interpolation qui est continu avec ses dérivés dans un ordre spécifié dans l'élément. En
général, la fonction d'interpolation est de continuité d'ordre inférieur entre les éléments
qu'à l'intérieur d'un élément. Ainsi, le bloc de construction fondamental de la méthode des
éléments finis est le sous-domaine ou l'élément.

Ainsi, alors que dans la méthode des différences finies ce sont des dérivées de la
fonction recherchées qui sont remplacées par des approximations, dans la méthode des
éléments finis, c’est la fonction inconnue recherchée qui est remplacée par une fonction
d’interpolation arbitrairement choisie par l’opérateur à condition qu’elle soit dérivable sur
plusieurs ordre et qu’elle réponde aux conditions aux limites prévues pour la fonction
recherchée.

Bien que les schémas de différences finies puissent être présentés à l'aide des
séries de Taylor d'une manière assez simple, la théorie derrière la méthode des éléments
finis est relativement abstraite. Il sera nécessaire d'introduire des concepts issus de l'analyse

52
fonctionnelle et des méthodes variationnelles afin de formuler les équations algébriques
analogues aux formules aux différences finies. De plus, il faut surmonter le problème
conceptuel de l'assemblage des informations obtenues élément par élément en une
représentation globale d'un problème.

7.2. LA METHODE DES RESIDUS PONDERES

Il existe plusieurs voies qui mèneront à la même formulation par éléments finis.
Une approche conceptuellement simple, mais mathématiquement rigoureuse, peut être
formulée en utilisant la méthode des résidus pondérés. Bien que dans notre discussion sur
les techniques des éléments finis, nous considérerons principalement deux cas particuliers
de la méthode des résidus pondérés (MWR) - les méthodes de collocation et de Galerkin,
d'autres schémas découlant de MWR se trouvent dans la littérature technique. Une
discussion sur l'application du MWR en sciences de l'ingénieur peut être trouvée dans
Finlayson (1972).

Dans la méthode des résidus pondérés, la fonction souhaitée u ( . ) est remplacée

par une approximation en série finie û ( . ):

j= N
u ( . ) ≈ û ( . )= ∑ U j ∅ j ( . ).................................................................(5.2.1)
j=1

En général, l'ensemble des fonctions ∅ j ( . ), j=1 ,2 , 3 , … , N , peuvent être définies à la fois sur
le domaine temporel et spatial et U j ( . ), j=1 ,2 , 3 , … , N sont des coefficients indéterminés.
N−1
L'équation (5.2.1) peut s'écrire û ( . )=U o ∅ o ( . ) + ∑ U j ∅ j ( . ) où ∅ j ( . ) satisfait les conditions aux
j=1

limites homogènes. Dans la méthode des éléments finis, les fonctions ∅ j ( . ) sont choisies

53
pour être des polynômes qui satisfont certaines des conditions aux limites imposées au
problème. Ces fonctions sont désignées de différentes manières comme des fonctions de
forme, des fonctions de base et des fonctions d'interpolation, selon la discipline dans
laquelle la méthode est appliquée.

Même lorsque les fonctions de base sont choisies pour satisfaire toutes les conditions aux
limites imposées à un problème, elles ne satisferont normalement pas non plus la PDE. La
substitution de û ( . ) dans la PDE, Lû−f =0, donne un résidu, R,

Lû ( . )−f =R (.) ……………………………………………………………………………………….(5.2.2)

L'objectif est de sélectionner les coefficients Ul indéterminés de telle sorte que ce résidu soit
minimisé en un sens. Un schéma simple serait de mettre l'intégrale de R(.) à zéro :

∫∫ R ( . ) dVdt =0………………………………………………………………………………………(5.2.3)
Ce schéma, cependant, ne génère qu'une seule équation pour les N coefficients inconnus U j.
Il peut être convenablement modifié en introduisant les fonctions de pondération ∅ j ( . ),
j=1 ,2 , 3 , … , N .
La mise à zéro de l'intégrale de chaque résidu pondéré donne N équations indépendantes

∫∫ R ( . ) w i ( . ) dVdt =0, i=1 , 2 ,3 , … . N …………………………………………………………………(5.2.4)

En théorie, au moins, (5.2.4) peut être résolu pour les N coefficients inconnus. L'équation
(5.2.4) est l'équation générale décrivant le MWR, et une multiplicité de schémas découlent
de cette expression unique par la définition des fonctions de pondération wi.

Parmi la famille de méthodes MWR, les schémas de Galerkin et de collocation sont les plus
couramment rencontrés dans la pratique de l'ingénierie. Il est intéressant de noter comment
le support de ces fonctions change dans chaque méthode. Comme on pouvait s'y attendre,
l'effort de calcul requis dans la formulation des matrices de coefficients est lié au support de
la fonction. Examinons maintenant la forme spécifique de l'équation générale MWR (5.2.4)
qui se pose pour chaque schéma.

7.3. EXEMPLE D’UTILISATION DE LA METHODE DES RESIDUS PONDERES

Comme nous venons de le dire au point 5.2., la méthode des résidus pondérés consiste à
approcher partiellement l’annulation du résidu d’une équation différentielle pour trouver
une solution discrète approximative. Illustrons le concept par l’exemple qui suit.

Soit l'équation différentielle et la condition initiale suivante :

54
du
=−u …………………………………………………………………………………………………….(5.3.1)
dx

avec la condition aux limites : u=1quand x=0

Supposons que l'on désire connaître la solution dans l'intervalle : 0 ≤ x ≤ 1

Posons alors une solution approchée :

2
û=c o +c 1 x +c 2 x ……………………………………………………………………………..(5.3.2)

Cette solution approchée doit satisfaire les conditions aux limites donc c 0=1 et c 1 et c 2
restent des coefficients indéterminés. Afin de les déterminer nous devons émettre un critère
qui nous permettra d'ajuster ces coefficients à la vraie solution.

Considérons le résidu :

dû
R ( x )= + û………………………………………………………………………………………..(5.3.3)
dx

Si û est la solution exacte, R(x) = 0. La solution approchée û donne

R ( x )=1+c 1 (1+ x ) +c 2 (2 x+ x)…………………………………………………………..(5.3.4)

Nous devons donc essayer de rendre ce résidu nul par un moyen quelconque. Le choix de
cette fonction approchée doit satisfaire les conditions aux limites ou initiales. Comme nous
l’avons dit au point (5.2) examinons la méthode de collocation et celle de Galerkin.

7.3.1. Méthode de collocation par points

Comme nous avons deux paramètres à rechercher, c 1 et c2, il faut alors écrire deux
équations, ce qui va nous donner un système de deux équations et deux inconnus. Pour cela,
1 2
imposons que le résidu s'annule en deux points de l'intervalle, soit à x= et à x= . Nous
3 3
écrivons alors :

R ( 13 )=1+c (1+ 31 )+c ( 23 + 19 )=1+ 43 c + 79 c =0

1 2 1 2

R ( )=1+c (1+ )+c ( + )=1+ c + c =0

2 2 4 4 5 16
1 2 1 2
3 3 3 9 3 9

La résolution de ce système donne :

c 1=−0,931
c 2=0,3103

55
7.3.2. Méthode de collocation par sous-domaines

Choisissons un autre critère : la moyenne du résidu sur deux sous intervalles (parce que nous
devons avoir un système d’équation à deux inconnus), constituant le

[ ] [ ]
1 1
domaine, ici, 0 ; et ;1 , doit s'annuler :
2 2

1
2

∫ Rdx= 12 + 58 c 1+ 24
7
c 2=0
0

∫ Rdx= 12 + 78 c 1+ 25 c =0
24 2
1
2

La résolution de ce système donne :

c 1=−0,9474
c 2=0,3158

7.3.3. Méthode de Galerkin

Le principe consiste à considérer que les moyennes pondérées sur l'ensemble du domaine
s'annulent. Les pondérations sont choisies parmi les fonctions qui ont servi à construire la
solution approchée. Ici, il y a deux fonction, à savoir, x et x 2, ce qui servira à obtenir un
système d’équation à deux inconnus :

∫ xRdx= 12 + 56 c 1+ 11 c =0
12 2
0

∫ x 2 Rdx= 13 + 12
7 9
c 1+ c2 =0
20
0
La résolution de ce système donne :

c 1=−0,9143
c 2=0,2857

7.3.4. Méthode des moindres carrés

On considère que la moyenne sur l’intervalle du carré du résidu doit être minimum.

1 1
1 ∂ ∂R 3 7 9
∫
2 ∂ c1 0
R dx=∫ R
2
R c1
dx= + c 1+ c2 =0
2 3 4
0

56
 1 1
1 ∂ ∂R 4 9 38
∫
2 ∂ c2 0
R dx=∫ R
2
R c2
dx= + c 1+ c 2=0
3 3 15
0

La résolution de ce système donne :

c 1=−0,9427
c 2=0,311

7.3.5. Solution exacte

La solution générale s’écrit par séparation de variable :

1
du=−dx
u

1
∫ u du=−dx
ln ( u ) =−x+ c

−x
u=C e

La solution particulière s’obtient par introduction des conditions aux limites :

0
1=C e
D’où
C=1

L’équation particulière pour ces conditions aux limites est :

−x
u=e

7.3.6. Représentation graphique des résultats

57
 Figure 5.3.1. Comparaison des graphiques f(x) = u

Figure 5.3.2. Comparaison des graphiques f(x) = R(u)

Figure 5.3.3. Comparaison des graphiques f(x) = erreur

7.3.7. Généralisation de la méthode

58
 1

∫ w1 Rdx=0
0
1

∫ w2 Rdx=0
0

Où w1 et w2 sont des fonctions de pondération

 Fonctions de Dirac pour la collocation par points :

w 1=δ x− ( 13 ) et w =δ ( x− 23 )
2

 Fonctions « échelon unité » pour la collocation par sous domaine :

w 1=u0 ( x )−u 1 (x ) et w 2=u 1 (x )

2 2
 Fonctions de base de l’approximation pour Galerkine :

w 1=x et w 2=x 2

 Fonctions dérivées du résidu par rapport aux coeffcicients pour les moindres carrés :

∂R ∂R
w 1= et w 2=
∂ c1 ∂ c2

7.3.8. Approximation nodale

Supposons que l'on veuille donner maintenant une signification physique aux coefficients
indéterminés.

Dans l’exemple précédent, les coefficients indéterminés c i n’ont pas de signification

physique, il serait intéressant d’exprimer ces coefficients en fonction de valeurs discrètes de
la fonction solution.

Prenons une solution approchée

2
u ( x )=c 1 +c 2 x +c 3 x

1
Nous voulons exprimer le fait qu'aux points 0 , et 1, nous aurons des valeurs discrètes de la
2
solution u , u et u .
1 2 3

u1=u ( 0 )=c 1 +c 2 0+ c 3 0
u2=u
1
2() 1
=c 1 +c 2 + c 3
2
1
4

59
 u1=u ( 1 ) =c 1+ c2 1+c 3 1
Sous forme matricielle :

Donc

Avec

D’autre part la solution approchée s’écrit :

Connaissant l’expression de { c n } , la solution approchée devient :

Enfin, en forme algébrique :

Et la dérivée première :

∂u
=(−3+ 4 x ) u1+ ( 4−8 x ) u 2+ (−1+ 4 x ) u3
∂x
Donc le résidu s’écrit :

60
 R ( x )=(−2+ x +2 x 2) u1 + ( 4−4 x−4 x 2 ) u2+ (−1+3 x+ 2 x 2 ) u3
Appliquons une méthode de collocation par points en utilisant deux points du domaine, soit
1 2
et :
3 3

En posant x1 = 1 comme condition initiale :

Les solutions approchée et exacte sont :

7.3.9. Exercice 5.3.1.

Enoncé du problème

Résoudre le problème suivant par une méthode de résidu pondéré de votre choix en
prenant une fonction approchée du second degré pour 0 < t < 5.

du
=1−u .
dt

Pour t ≥ 0 et u=0 quand t=0

Solution

Si u(t) est la solution exacte, alors

du
R ( t )= −1+u=0.
dt

Soit une solution approchée qui satisfait les conditions aux limites :

2
u ( t )=a1 t+ a2 t
Alors

61
 du
=a 1+2 a 2 t
dt

Le résidu approché s’écrit :

2
R ( t )=a1 +2 a2 t−1+a1 t+ a2 t
2
R ( t )=−1+(1+ t)a1+(2 t+ t )a 2

Avec une méthode de collocation par points, posons R(1) = 0 et R(4) = 0

Cela donne deux équations avec deux inconnus :

2 a1 +3 a2 =1
5 a1 +24 a2=1
Donc
a 1=0,6363
a 2=−0,0909

62