Modélisation

Modélisation
Introduction aux équations aux dérivées
partielles et leurs résolutions numériques

Exercice
On considère l'équation de propagation d'un débit en rivière qui s'écrit :
∂Q( x , t ) ∂Q( x , t ) ∂ 2 Q( x , t )
+ c(Q( x , t )) = σ(Q( x , t )) + c(Q( x , t ))q ( x; t ) ,
∂t ∂x ∂x 2
où Q(x,t) est le débit dans la rivière
c(Q) est la célérité des ondes ([c]=LT-1),
σ(Q) est le coefficient d'atténuation ([σx]=L2T-1),
q(x,t) est la fonction d'apport intermédiaire le long du bief par mètre de longueur
([q]=L2T-1) (terme éventuellement négatif)

Cette équation est dite "onde diffusante" (car il est évident, et donc sous-entendu, qu'il y a de
la convection en rivière), mais c'est une équation de diffusion-convection.
On prendra la célérité et le coefficient d'atténuation constants, et le terme d'apport le long du
bief nul.
∂Q( x , t ) ∂Q( x , t ) ∂ 2 Q( x , t )
L'équation simplifiée devient : +c =σ
∂t ∂x ∂x 2

On applique cette équation dans un canal d'irrigation. Le débit est constant dans le canal à 10
litre/s (0.01 m3/s) et on ouvre brutalement une vanne à un débit de 100 litres/s (0.1 m3/s). On
veut donc savoir comment l'augmentation du débit dans le canal va se dérouler.

1 – Réfléchir aux conditions limites qu'on peut mettre à l'amont et à l'aval du bief pour effectuer une
propagation de débit avec cette équation. Vous pouvez au moins réfléchir au type de condition qu'on
peut mettre (Dirichlet ou Neumann) –
(Avant de faire la 2ème question, regardez la correction de cette 1ère question pour partir sur des bases
"saines")

2 – Résoudre cette équation avec la méthode des différences finies. On appliquera un schéma
explicite, implicite et pondéré
(Là aussi, regardez la correction après avoir fait vos calculs pour les vérifier avant de faire
la troisième question)

On applique cette équation dans un canal d'irrigation. Le débit est constant dans le canal à 10
litre/s (0.01 m3/s) et on ouvre une vanne à un débit de 100 litres/s (0.1 m3/s) en 30 secondes
environ.

On prend un pas d'espace de 100m sur une distance totale de 500m du canal.
On estime une célérité de 0.15m/s et un coefficient d'atténuation de 100m²/s (ces 2 dernières
valeurs se callent sur des observations ou se calcul d'après des formules théoriques).

3 – Calculer en schéma explicite (avec un pas de 100m, et sur une longueur de 500m) la
valeur du débit sur 3 pas de temps. (ce calcul peut se faire avec une simple calculette, mais
c'est plus rapide avec Excel si vous savez le manipuler).

FN CRES
Modélisation

On prendra un pas de temps de 30s.

Tracer sur un même graphique les courbes Q(x) pour chaque pas de temps

Même calcul et graphique avec un pas de temps de 90s

Correction
∂Q( x , t ) ∂Q( x , t ) ∂ 2 Q( x , t )
L'équation est : +c =σ
∂t ∂x ∂x 2

On veut simuler la propagation d'un débit dans le bief.

1 – Réfléchir aux conditions limites qu'on peut mettre à l'amont et à l'aval du bief pour effectuer une
propagation de débit avec cette équation. Vous pouvez au moins réfléchir au type de condition qu'on
peut mettre (Dirichlet ou Neumann) –

Il y a 2 conditions limites : une à l'amont, et une à l'aval.

CL amont
Cette condition est physiquement très logique : on injecte un débit Q dans le bief, donc le
débit est connu à l'entrée du bief. Il s'agit donc d'une condition de type Dirichlet. Ici, on
suppose que le débit restera constant à 0.1m3/s et ce sera notre condition limite pour t > 0.

CL aval
Cette condition est plus délicate. On ne peut pas mettre une condition de type Dirichlet
puisque cela voudrait dire que quelque soit le débit à l'amont (à l'entrée du bief), la valeur du
débit à l'aval du bief est connue.
Reste des conditions de type Neumann : comment le justifier et quelle valeur ?
(Remarque : je rappelle qu'il y a d'autres conditions qu'on peut mettre en œuvre, mais on ne
les voit pas dans le cadre de ce cours. On en restera donc à Dirichlet et Neumann, ce qui
nous suffit ici).

Un des avantages d'une condition de type Neumann est que la valeur de la dérivée sera fixée,
mais pas celle du débit. Autrement dit, l'inconvénient évoqué si on prend une condition de
Dirichlet disparaît effectivement.
La justification qu'on peut avancer consiste à dire qu'on applique une équation qui a tendance
à atténuer les débits dans l'espace : si vous propagez un débit très "pointu", cette pointe
diminuera au fur et à mesure qu'on s'éloigne de l'amont (c'est l'effet de la diffusion, qu'on
appelle ici "atténuation"). Si on se met suffisamment loin de l'amont, la pointe de débit aura
même disparue et le débit sera constant (la diffusion aura fait tout son effet).
Dans ce cas, la dérivée du débit par rapport à l'espace est nulle (c'est une droite horizontale).
On a alors ∂Q/∂x = 0. C'est la condition qu'on applique à l'aval … mais il faut être
suffisamment loin !!!! Pour s'en assurer, on simule généralement sur une distance supérieure à
la distance "réelle" de façon à ce que cette condition soit "le plus vraie possible" (mais on
s'intéresse quand même à la valeur du débit à l'aval du bief "réel", même si on mène les
calculs plus loin). En fait, l'influence de cette condition, en termes de calculs, n'est pas très
forte, c'est-à-dire que son influence ne remonte pas très loin dans le bief. En tout cas, elle est
beaucoup moins contraignante qu'une condition de type Dirichlet.

Ici, pour cet exercice, on prendra cette condition à 500m (même si on sait qu'on devrait
rallonger le bief conformément à ce qui est dit au-dessus).

FN CRES
Modélisation

2 – Résoudre cette équation avec la méthode des différences finies. On appliquera un schéma
explicite, implicite et pondéré

Discrétisation de l'espace (quelque soit le schéma)

t
Q = CLamont ∂Q/∂x = 0

t+∆t
t
x
t0
0 1 i-1 i i+1 N

Schéma explicite
Discrétisation
∂Q Q it +1 − Q it ∂Q Q it+1 − Q it−1 ∂ 2 Q Q it−1 − 2Q it + Q it+1
= ; = ; =
∂t ∆t ∂x 2∆x ∂x 2 ∆x 2

Equation discrétisée
Qit +1 − Qit Q t − Q it−1 Q t − 2Qit + Qit+1
+ c i +1 = σ i −1
∆t 2∆x ∆x 2
En regroupant les termes
1  c σ   1 2σ  t  c σ 
Q it +1 + Q it−1  − − 2 
+ Q it  − +  + Q i +1  − =0
∆t  2∆x ∆x   ∆t ∆x
2
  2∆x ∆x
2


Soit en écrivant le nœud au temps (t+1) inconnu à gauche et ceux au temps (t) connus à droite
  c σ   1 2σ   c σ 
Qit +1 = ∆t Q it−1  + 2  + Qit  − 2  + Qit+1  − + 2 
  2∆x ∆x   ∆t ∆x   2∆x ∆x 
Cette équation générique permet de calculer la valeur de débit au nœud i connaissant les 3
nœuds du pas de temps précédent.

Le premier nœud et le dernier nœud sont particuliers.

Au premier nœud (n°1)
  c σ   1 2σ   c σ 
Q1t +1 = ∆t CL amont  + 2  + Q1t  − 2  + Q 2t  − + 2 
  2∆x ∆x   ∆t ∆x   2∆x ∆x 

Au dernier nœud (n°N)

  c σ   1 2σ   c σ 
Q Nt +1 = ∆t Q Nt −1  + 2  + Q tN  − 2  + Q tN +1  − + 2 
  2∆x ∆x   ∆t ∆x   2∆x ∆x 

FN CRES
Modélisation

L'équation discrétisée introduit un nœud fictif. On discrétise la CL aval

 ∂Q  Q t − Q tN −1
  = N +1 = 0 d'où Q Nt +1 = Q tN −1
 ∂ x N 2 ∆ x
En remplaçant dans l'équation précédente
  2σ  1 2σ  
Q tN+1 = ∆t Q tN −1  2  + Q tN  − 2 
  ∆x   ∆t ∆x 

On applique donc ces différentes équations aux différents nœuds, en utilisant pour le premier
pas de temps la condition initiale dans le canal, c'est-à-dire Q init = Q i0 = 0.01m 3 / s quelque
soit l'abscisse i.

Schéma implicite
Discrétisation
∂Q Q it +1 − Q it ∂Q Q it++11 − Q it−+11 ∂ 2Q Q it−+11 − 2Q it +1 + Q it++11
= ; = ; =
∂t ∆t ∂x 2∆x ∂x 2 ∆x 2

Equation discrétisée
Q it +1 − Q it Q t +1 − Q it−+11 Q t +1 − 2Q it +1 + Q it++11
+ c i +1 = σ i −1
∆t 2∆x ∆x 2
En regroupant les termes
 c σ  1 2σ   c σ   1 
Q it−+11  − − 2 
+ Q it +1  + 2 
+ Q it++11  − 2 
+ Q it  −  = 0
 2∆x ∆x   ∆t ∆x   2∆x ∆x   ∆t 

Soit en écrivant les nœuds au temps (t+1) inconnus à gauche et celui au temps (t) connu à
droite :
 c σ  1 2σ   c σ  1 
Q it−+11  − − 2 
+ Q it +1  + 2 
+ Q it++11  − 2 
= Q it  
 2∆x ∆x   ∆t ∆x   2∆x ∆x   ∆t 
Cette équation générique écrite au nœud i relie donc 3 nœuds inconnus au temps (t+1) avec 1
nœud connu au temps (t); on écrit autant d'équations que d'inconnues et on résout le système
matriciel à chaque pas de temps.

Le premier nœud et le dernier nœud sont particuliers.

Au premier nœud (n°1)
 c σ   1 2σ   c σ   1
CL amont  − − 2  + Q1t +1  + 2  + Q 2t +1  − 2  = Q1t  
 2∆x ∆x   ∆t ∆x   2∆x ∆x   ∆t 
 1 2σ   c σ   1  c σ 
Soit Q1t +1  + 2  + Q 2t +1  − 2  = Q1t   + CL amont  + 2
 ∆t ∆x   2∆x ∆x   ∆t   2∆x ∆x 

Au dernier nœud (n°N)

 c σ   1 2σ   c σ  1 
Q tN+−11  − − 2  + Q tN+1  + 2  + Q tN++11  − 2  = Q tN  
 2∆x ∆x   ∆t ∆x   2∆x ∆x   ∆t 
L'équation discrétisée introduit un nœud fictif. On discrétise la CL aval
 ∂Q  Q Nt ++11 − Q Nt +−11
  = = 0 d'où QNt ++11 = QNt +−11
 ∂x  N 2∆x

FN CRES
Modélisation

En remplaçant dans l'équation précédente

 2σ   1 2σ   1
Q Nt +−11  − 2  + Q tN+1  + 2  = Q tN  
 ∆x   ∆t ∆x   ∆t 

On peut écrire le système matriciel obtenu à chaque pas de temps.

 c σ   1 2σ   c σ  1 
On pose : A =  − − 2 
; B= + 2 
; C= − 2 
; D= 
 2∆x ∆x   ∆t ∆x   2∆x ∆x   ∆t 

Q1t +1 Q t2+1 Q 3t +1 ... Q tN+−11 Q Nt +1

 B C   Q1t +1  DQ1t − A × CL amont 
 A   t +1   
 B C   Q2   DQ 2t 
 A B C   Q 3t +1   DQ 3t 
     
 ... ... ...   ...  =  ... 
 A B C   Q t +1   DQ i t 
   i   
 ... ... ...   ...   ... 
 A B C  Q t +1   t
DQ N −1 
   N −1   
 A+C B   Q tN+1   DQ tN 

Seule la matrice de droite change à chaque pas de temps. Ce système doit être résolu à chaque
pas de temps.

Schéma pondéré
Discrétisation

∂Q Q it +1 − Q it
= ;
∂t ∆t
∂Q Q it++11 − Q it−+11 Q it+1 − Qit−1
=θ + (1 − θ) ;
∂x 2∆x 2∆x
∂ 2Q Q it−+11 − 2Qit +1 + Q it++11 Qit−1 − 2Qit + Qit+1
= θ + (1 − θ )
∂x 2 ∆x 2 ∆x 2

Equation discrétisée
Q it +1 − Q it Q t +1 − Q it−+11 Q t − Q it−1 Q t +1 − 2Q it +1 + Q it++11 Q it−1 − 2Q it + Q it+1
+ cθ i +1 + c(1 − θ) i +1 = σθ i −1 + σ(1 − θ )
∆t 2∆x 2∆x ∆x 2 ∆x 2

En regroupant les termes

 θc θσ   1 2θσ   θc θσ 
Q it−+11  − − 2 
+ Q it +1  + 2 
+ Q it++11  − 2 
+
 2∆x ∆x   ∆t ∆x   2∆x ∆x 
 σ   1 2(1 − θ)σ   c σ 
Q it−1 (1 − θ) −  + Q i +1 (1 − θ)
c
− 2 
+ Q it  − + t
− 2 
=0
 2∆x ∆x   ∆t ∆x  2∆x ∆x 
2


FN CRES
Modélisation

Soit en écrivant les nœuds au temps (t+1) inconnus à gauche et celui au temps (t) connu à
droite :

 θc θσ   1 2θσ   θc θσ 
Q it−+11  − − 2 
+ Q it +1  + 2 
+ Q it++11  − 2 
=
 2∆x ∆x   ∆t ∆x   2∆x ∆x 
 c σ   1 2(1 − θ)σ   σ 
Q it−1 (1 − θ)  + Q i +1 (1 − θ) −
c
+ 2 
+ Q it  − t
+ 2 
 2∆x ∆x   ∆t ∆x  2∆x ∆x 
2


Cette équation générique écrite au nœud i relie donc 3 nœuds inconnus au temps (t+1) avec 3
nœuds connus au temps (t); on écrit autant d'équations que d'inconnues et on résout le système
matriciel à chaque pas de temps.
Remarque : pour θ = 0, on retrouve le schéma explicite, pour θ = 1, on retrouve le schéma
implicite.

Le premier nœud et le dernier nœud sont particuliers.

Au premier nœud (n°1)
 θc θσ   1 2θσ   θc θσ 
CL amont  − − 2 
+ Q1t +1  + 2 
+ Q 2t +1  − 2 
=
 2∆x ∆x   ∆t ∆x   2∆x ∆x 
 c σ   1 2(1 − θ)σ   σ 
CL amont (1 − θ)  + Q 2 (1 − θ) −
c
+ 2 
+ Q1t  − t
+ 2 
 2∆x ∆x   ∆t ∆x  2∆x ∆x 
2

soit

 1 2θσ   θc θσ 
Q1t +1  + 2 
+ Q 2t +1  − 2 
=
 ∆t ∆x   2∆x ∆x 
 c σ   1 2(1 − θ)σ   σ   θc θσ 
CL amont (1 − θ)  + Q 2 (1 − θ) −
c
+ 2 
+ Q1t  − t
+ 2 
+ CL amont  + 2 
 2∆x ∆x   ∆t ∆x  2∆x ∆x   2∆x ∆x 
2


Au dernier nœud (n°N)

 θc θσ   1 2θσ   θc θσ 
Q tN+−11  − − 2 
+ Q tN+1  + 2 
+ Q tN++11  − 2 
=
 2∆x ∆x   ∆t ∆x   2∆x ∆x 
 c σ   1 2(1 − θ)σ   σ 
Q tN −1 (1 − θ)  + Q N +1 (1 − θ) −
c
+ 2 
+ Q tN  − t
+ 2 
 2∆x ∆x   ∆t ∆x  2∆x ∆x 
2


L'équation discrétisée introduit 2 nœuds fictifs. On discrétise la CL aval

 ∂Q  Q t − Q tN −1
  = N +1 = 0 d'où Q Nt +1 = Q tN −1
 ∂x  N 2∆x
 ∂Q  Q tN++11 − Q tN+−11
  = = 0 d'où Q tN++11 = Q tN+−11
 ∂x  N 2∆x
En remplaçant dans l'équation précédente
 θσ   1 2θσ   σ   1 2(1 − θ)σ 
Q tN+−11  − 2 2  + Q tN+1  + 2 
= Q tN −1 (1 − θ) 2 2  + Q tN  − 
 ∆x   ∆t ∆x   ∆x   ∆t ∆x 2 

On peut écrire le système matriciel obtenu à chaque pas de temps.

FN CRES
Modélisation

 θc θσ   1 2θσ   θc θσ 
On pose A =  − − 2 
; B= + 2 
; C= − 2 
 2∆x ∆x   ∆t ∆x   2∆x ∆x 
 c σ   1 2(1 − θ)σ   σ 
D = (1 − θ)  ; F = (1 − θ) −
c
+ 2 
; E= − + 2 
 2∆x ∆x   ∆t ∆x  2∆x ∆x 
2


Q 1t +1 Q 2t +1 Q 3t +1 ... Q Nt +−11 Q tN+1

 B C   Q 1t +1 
 A   t +1 
 B C   Q2 
 A B C   Q 3t +1 
   
 ... ... ...   ...  =
 A B C   Q t +1 
   i 
 ... ... ...   ... 
 A B C  Q t +1 
   N−2 
 A+C B   Q tN+−11 
D × CL amont + EQ1t + FQ 2t − A × CL amont 
 
 DQ 1t + EQ t2 + FQ 3t 
 DQ 2t + EQ 3t + FQ t4 
 
 ... 
 DQ i −1 + EQ i + FQ i +1
t t t 
 
 ... 
 DQ N − 2 + EQ N −1 + FQ N
t t t 
 
 (D + F)Q tN −1 + EQ tN 

Seule la matrice de droite change à chaque pas de temps. Ce système doit être résolu à chaque
pas de temps.

On applique cette équation dans un canal d'irrigation. Le débit est constant dans le canal à
10 litre/s (0.01 m3/s) et on ouvre une vanne à un débit de 100 litres/s (0.1 m3/s) en 30
secondes environ.

On prend un pas d'espace de 100m sur une distance totale de 500m du canal.
On estime une célérité de 0.15m/s et un coefficient d'atténuation de 100m²/s (ces 2 dernières
valeurs se callent sur des observations ou se calcul d'après des formules théoriques).
Pour un débit de 0.1m3/s, on estime que la vitesse V de l'eau sera environ de 0.25m/s, et la
profondeur h d'environ 0.5m

4 – Calculer en schéma explicite (avec un pas de 100m, et sur une longueur de 500m) la
valeur du débit sur 3 pas de temps. (ce calcul peut se faire avec une simple calculette, mais
c'est plus rapide avec Excel si vous savez le manipuler).

FN CRES
Modélisation

On prendra un pas de temps de 30s.

Tracer sur un même graphique les courbes Q(x) pour chaque pas de temps

Même calcul et graphique avec un pas de temps de 90s

Remarque : "à la main", on ne peut appliquer raisonnablement que le schéma explicite,

puisque le schéma implicite nécessite la résolution d'un système d'équations (ce qui peut aussi
se faire à la main, mais non sans mal !)

dt = 30 s
On pense que la vanne sera relevée suffisamment rapidement pour avoir les 0,1 m3/s très
rapidement. C'est pour cela qu'on prend cette valeur à l'instant 0.

Remarque : à l’abscisse x=0, on prend en fait la condition limite quelque soit le temps,
notamment à t=0, alors qu’en toute logique, on pourrait prendre la condition initiale. C’est
un problème qui peut apparaître à "l'intersection" de 2 conditions non compatibles : il faut
choisir. Ici, à t=0, il est préférable de prendre la condition limite car on va plutôt considérer
que la vanne est ouverte un temps proche de t=0, et qu'à partir de cet instant, l'eau entre dans
le canal. Sinon, tout se passe comme si la vanne était ouverte à t = ∆t.

On voit la perturbation se propager dans le canal avec la valeur 0.1m3/s qui gagne petit à petit
tout le bief.

FN CRES
Modélisation

A titre informatif, j'ai repris ci-dessous le calcul, toujours avec un pas de temps de 30s (même
si je n'ai pas mis toutes les étapes sur le graphique) jusqu'à 1020s. Ce qui est dit dans l'alinéa
précédent apparaît clairement.

dt = 90 s

FN CRES
Modélisation

Pas besoin d'être hydraulicien pour voir qu'il y a un problème qui apparaît dès le 2ème pas de
temps! Le schéma avec ce pas de temps n'est pas convergent.

FN CRES