Équation de diffusion à une dimension

1. Équation de diffusion

Soit une fonction u(x,t) représentant la température dans un problème de diffusion thermique, ou la concentration pour un problème de diffusion de particules. L'équation de diffusion est :

\frac{\partial u}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}} + s (x, t)

où D est le coefficient de diffusion et s(x,t) représente une source, par exemple une source thermique provenant d'un phénomène de dissipation.

On cherche une solution numérique de cette équation pour une fonction s(x,t) donnée, sur l'intervalle [0,1], à partir de l'instant t=0. La condition initiale est u(x,0). Sur les bords (x=0 et x=1) la condition limite est soit de type Dirichlet :

u (0, t) = U_{0}

soit de type Neumann (dérivée imposée) :

\frac{\partial u}{\partial x} (0, t) = J_{0}

2. Méthode des différences finies

2.a. Définitions

Soit N le nombre de points dans l'intervalle [0,1]. On définit le pas de x par

Δ x = \frac{1}{N - 1}

On définit aussi le pas du temps $Δ t$ . La discrétisation de u(x,t) est définie par :

U_{j}^{n} = u (j Δ x, n Δ t)

où j est un indice variant de 0 à N-1 et n un indice positif ou nul représentant le temps.

Figure pleine page

La discrétisation du terme de source est

S_{j}^{n} = u (j Δ x, n Δ t)

On pose

α = \frac{D Δ t}{{(Δ x)}^{2}}

2.b. Schéma explicite

Pour discrétiser l'équation de diffusion, on peut écrire la différence finie en utilisant les instants n et n+1 pour la dérivée temporelle, et la différence finie à l'instant n pour la dérivée spatiale :

\frac{U_{j}^{n + 1} - U_{j}^{n}}{Δ t} = D \frac{U_{j + 1}^{n} - 2 U_{j}^{n} + U_{j - 1}^{n}}{{(Δ x)}^{2}} + S_{j}^{n}

Avec ce schéma, on peut calculer les U_jⁿ⁺¹ à l'instant n+1 connaissant tous les U_jⁿ à l'instant n, de manière explicite. Ce schéma est précis au premier ordre ([1]). Comme montré plus loin, sa stabilité n'est assurée que si le critère suivant est vérifié :

\frac{D Δ t}{{(Δ x)}^{2}} \leq \frac{1}{2}

En pratique, cela peut imposer un pas de temps trop petit.

L'implémentation de cette méthode est immédiate. Voici un exemple :

import numpy
from matplotlib.pyplot import *

N=100
x=numpy.linspace(0,1,N)
dx=x[1]-x[0]
dx2=dx**2
U=numpy.zeros(N)
laplacien=numpy.zeros(N)
dt = 3e-5
U[0]=1
U[N-1]=0
D=1.0
for i in range(1000):
    for k in range(1,N-1):
        laplacien[k] = (U[k+1]-2*U[k]+U[k-1])/dx2
    for k in range(1,N-1):
        U[k] += dt*D*laplacien[k]

figure()
plot(x,U)
xlabel("x")
ylabel("U")
grid()

alpha=D*dt/dx2

print(alpha)
--> 0.29402999999999996

figA.pdf

Le nombre de points N et l'intervalle de temps sont choisis assez petits pour satisfaire la condition de stabilité. Pour ces valeurs, l'atteinte du régime stationnaire est très longue (en temps de calcul) car l'intervalle de temps Δt est trop petit. Si on augmente cet intervalle, on sort de la condition de stabilité :


N=100
x=numpy.linspace(0,1,N)
dx=x[1]-x[0]
dx2=dx**2
U=numpy.zeros(N)
laplacien=numpy.zeros(N)
dt = 6e-5
U[0]=1
U[N-1]=0
D=1.0
for i in range(1000):
    for k in range(1,N-1):
        laplacien[k] = (U[k+1]-2*U[k]+U[k-1])/dx2
    for k in range(1,N-1):
        U[k] += dt*D*laplacien[k]

figure()
plot(x,U)
xlabel("x")
ylabel("U")
grid()

alpha=D*dt/dx2

print(alpha)
--> 0.58805999999999992

figB.pdf

2.c. Schéma implicite de Crank-Nicolson

La dérivée seconde spatiale est discrétisée en écrivant la moyenne de la différence finie évaluée à l'instant n et de celle évaluée à l'instant n+1 :

\frac{U_{j}^{n + 1} - U_{j}^{n}}{Δ t} = \frac{D}{2} \frac{(U_{j + 1}^{n + 1} - 2 U_{j}^{n + 1} + U_{j - 1}^{n + 1}) + (U_{j + 1}^{n} - 2 U_{j}^{n} + U_{j - 1}^{n})}{{(Δ x)}^{2}} + S_{j}^{n}

Ce schéma est précis au second ordre. Contrairement au schéma explicite, il est stable sans condition. En revanche, les $U_{j}^{n + 1}$ à l'instant n+1 sont donnés de manière implicite. Il faut donc à chaque instant n+1 résoudre le système à N équations suivant :

- \frac{α}{2} U_{j - 1}^{n + 1} + (1 + α) U_{j}^{n + 1} - \frac{α}{2} U_{j + 1}^{n + 1} = \frac{α}{2} U_{j - 1}^{n} + (1 - α) U_{j}^{n} + \frac{α}{2} U_{j + 1}^{n} + Δ t S_{j}^{n}

Ce système est tridiagonal. On l'écrit sous la forme :

A_{j, j - 1} U_{j - 1}^{n + 1} + A_{j, j} U_{j}^{n + 1} + A_{j, j + 1} U_{j + 1}^{n + 1} = R_{j}

À chaque étape, on calcule la matrice colonne R et on résout le système $A U^{n + 1} = R$ . Pour j=0 et j=N-1, l'équation est obtenue par la condition limite.

On peut aussi écrire le membre de droite sous la forme :

R_{j} = B_{j, j - 1} U_{j - 1}^{n} + B_{j, j} U_{j, n} + B_{j, j + 1} U_{j + 1}^{n} + C_{j}

ce qui donne la forme matricielle

A U^{n + 1} = B U^{n} + C

2.d. Analyse de stabilité de von Neumann

L'analyse de stabilité de von Neumann ([2][3]) consiste à ignorer les conditions limites et le terme de source, et à rechercher une solution de la forme suivante :

U_{j}^{n} = σ^{n} exp (i β j)

Il s'agit d'une solution dont la variation spatiale est sinusoïdale, avec un nombre d'onde β. Toute solution de l'équation de diffusion sans source et sans condition limite doit tendre vers une valeur uniformément nulle au temps infini. La méthode numérique utilisée est donc stable si |σ|<1 quelque soit la valeur de β.

En reportant cette solution dans le schéma explicite, on obtient :

σ = 1 + 2 α (cos β - 1)

La valeur absolue maximale de σ est obtenue pour cos(β)=-1. On en déduit la condition de stabilité : $α < 1 / 2$ .

Pour le schéma de Crank-Nicolson, on obtient :

σ = \frac{1 - α (1 - cos β)}{1 + α (1 - cos β)}

|σ| est inférieur à 1, donc le schéma est inconditionnellement stable.

2.e. Discrétisation des conditions limites

La discrétisation de la condition de Dirichlet (en x=0) est immédiate :

U_{0}^{n} = U_{0}

On pose donc pour la première équation du système précédent :

A_{0, - 1} = 0, A_{0, 0} = 1, A_{0, 1} = 0, R_{0} = U_{0}

De même pour une condition limite de Dirichlet en x=1 on pose

A_{N - 1, N - 2} = 0, A_{N - 1, N - 1} = 1, A_{N - 1, N} = 0, R_{N - 1} = U_{1}

Une condition limite de Neumann en x=0 peut s'écrire :

U_{1}^{n + 1} - U_{0}^{n + 1} = Δ x J_{0}

ce qui donne

A_{0, - 1} = 0, A_{0, 0} = - 1, A_{0, 1} = 1, R_{0} = Δ x J_{0}

Cependant, cette discrétisation de la condition de Neumann est du premier ordre, alors que le schéma de Crank-Nicolson est du second ordre. Pour éviter une perte de précision due aux bords, il est préférable de partir d'une discrétisation du second ordre ([1]) :

U_{1}^{n + 1} - U_{- 1}^{n + 1} = 2 Δ x J_{0}

Un point fictif d'indice -1 a été introduit. Pour ne pas avoir d'inconnue en trop, on écrit le schéma de Crank-Nicolson au point d'indice 0 tout en éliminant le point fictif avec la condition ci-dessus ([1]). On obtient ainsi :

(1 + α) U_{0}^{n + 1} - α U_{1}^{n + 1} = (1 - α) U_{0}^{n} + α U_{1}^{n} - 2 α Δ x J_{0} + Δ t S_{0}^{n}

On obtient de la même manière la condition limite de Neumann en x=1 :

- α U_{N - 2}^{n + 1} + (1 + α) U_{N - 1}^{n + 1} = α U_{N - 2}^{n} + (1 - α) U_{N - 1}^{n} + 2 α Δ x J_{1} + Δ t S_{N - 1}^{n}

2.f. Milieux de coefficients de diffusion différents

On suppose que le coefficient de diffusion n'est plus uniforme mais constant par morceaux. Exemple : diffusion thermique entre deux plaques de matériaux différents. Soit une frontière entre deux parties située entre les indices j et j+1, les coefficients de diffusion de part et d'autre étant D₁ et D₂. Pour j-1 et j+1, on écrira le schéma de Crank-Nicolson ci-dessus. En revanche, sur le point à gauche de la frontière (indice j), on écrit une condition d'égalité des flux :

D_{1} \frac{\partial U}{\partial x} (x^{-}) = D_{2} \frac{\partial U}{\partial x} (x^{+})

qui se traduit par

D_{1} (U_{j}^{n + 1} - U_{j - 1}^{n + 1}) = D_{2} (U_{j + 1}^{n + 1} - U_{j}^{n + 1})

et conduit aux coefficients suivants

A_{j, j - 1} = - D_{1}, A_{j, j} = D_{1} + D_{2}, A_{j, j + 1} = - D_{1}, R_{j} = 0

2.g. Convection latérale

Un problème de transfert thermique dans une barre comporte un flux de convection latéral, qui conduit à l'équation différentielle suivante :

\frac{\partial T}{\partial t} = D \frac{\partial^{2} T}{\partial x^{2}} + s - C (T - T e)

où le coefficient C (inverse d'un temps) caractérise l'intensité de la convection et T_e est la température extérieure.

On pose β=CΔt. Le schéma de Crank-Nicolson correspondant à cette équation est :

U_{j}^{n + 1} - U_{j}^{n} = \frac{α}{2} ((U_{j + 1}^{n + 1} - 2 U_{j}^{n + 1} + U_{j - 1}^{n + 1}) + (U_{j + 1}^{n} - 2 U_{j}^{n} + U_{j - 1}^{n})) + S_{j}^{n} + \frac{β}{2} (2 U_{e} - U_{j}^{n + 1} - U_{j}^{n})

c'est-à-dire :

- \frac{α}{2} U_{j - 1}^{n + 1} + (1 + α + \frac{β}{2}) U_{j}^{n + 1} - \frac{α}{2} U_{j + 1}^{n + 1} = \frac{α}{2} U_{j - 1}^{n} + (1 - α - \frac{β}{2}) U_{j}^{n} + \frac{α}{2} U_{j + 1}^{n} + Δ t S_{j}^{n} + β U_{e}

3. Résolution du système tridiagonal

Les matrices A et B étant tridiagonales, une implémentation efficace doit stocker seulement les trois diagonales, dans trois tableaux différents. On écrit donc le schéma de Crank-Nicolson sous la forme :

a_{j} U_{j - 1}^{n + 1} + b_{j} U_{j}^{n + 1} + c_{j} U_{j + 1}^{n + 1} = d_{j} U_{j - 1}^{n} + e_{j} U_{j}^{n} + f_{j} U_{j + 1}^{n} + s_{j}^{n}

Les coefficients du schéma sont ainsi stockés dans des tableaux à N éléments a,b,c,d,e,f,s. On remarque toutefois que les éléments a₀, c_N-1, d₀ et f_N-1 ne sont pas utilisés.

Le système tridiagonal à résoudre à chaque pas de temps est :

a_{j} U_{j - 1} + b_{j} U_{j} + c_{j} U_{j + 1} = R_{j}

où l'indice du temps a été omis pour alléger la notation. Le second membre du système se calcule de la manière suivante :

\begin{matrix} R_{0} = e_{0} U_{0} + f_{0} U_{1} + s_{0} \\ R_{j} = d_{j} U_{j - 1} + e_{j} U_{j} + f_{j} U_{j + 1} + s_{j} p o u r 1 \leq j \leq N - 2 \\ R_{N - 1} = d_{N - 1} U_{N - 2} + e_{N - 1} U_{N - 1} + s_{N - 1} \end{matrix}

Le système tridiagonal s'écrit :

(\begin{matrix} b_{0} & c_{0} \\ a_{1} & b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} \\ a_{N - 1} & b_{N - 1} \end{matrix}) (\begin{matrix} U_{0} \\ U_{1} \\ U_{2} \\ U_{N - 1} \end{matrix}) = (\begin{matrix} R_{0} \\ R_{1} \\ R_{2} \\ R_{N - 1} \end{matrix}) (1)

La méthode d'élimination de Gauss-Jordan permet de résoudre ce système de la manière suivante.

Les deux premières équations sont :

\begin{matrix} b_{0} U_{0} + c_{0} U_{1} = R_{0} \\ a_{1} U_{0} + b_{1} U_{1} + c_{1} U_{2} = R_{1} \end{matrix}

b₀ est égal à 1 ou -1 suivant le type de condition limite. On divise la première équation par ce coefficient, ce qui conduit à poser :

\begin{matrix} β_{0} = b_{0} \\ γ_{0} = \frac{c_{0}}{β_{0}} \\ x_{0} = \frac{R_{0}}{β_{0}} \end{matrix}

La première élimination consiste à retrancher l'équation obtenue multipliée par $a_{1}$ à la seconde :

(b_{1} - a_{1} γ_{0}) U_{1} + c_{1} U_{2} = R_{1} - a_{1} x_{0}

On pose alors :

\begin{matrix} β_{1} = b_{1} - a_{1} γ_{0} \\ γ_{1} = \frac{c_{1}}{β_{1}} \\ x_{1} = \frac{R_{1} - a_{1} x_{0}}{β_{1}} \end{matrix}

On construit par récurrence la suite suivante :

\begin{matrix} β_{k} = b_{k} - a_{k} γ_{k - 1} & (2) \\ γ_{k} = \frac{c_{k}}{β_{k}} & (3) \\ x_{k} = \frac{R_{k} - a_{k} x_{k - 1}}{β_{k}} & (4) \end{matrix}

Considérons la kième équation réduite et la suivante :

\begin{matrix} U_{k} + γ_{k} U_{k + 1} = x_{k} \\ a_{k + 1} U_{k} + b_{k + 1} U_{k + 1} + c_{k + 1} U_{k + 2} = R_{k + 1} \end{matrix}

La réduction de cette dernière équation est :

(b_{k + 1} - a_{k + 1} γ_{k}) U_{k + 1} + c_{k + 1} U_{k + 2} = R_{k + 1} - a_{k + 1} x_{k}

ce qui justifie la relation de récurrence définie plus haut.

Pour finir, voyons les deux dernières équations :

\begin{matrix} U_{N - 2} + γ_{N - 2} U_{N - 1} = x_{N - 2} \\ a_{N - 1} U_{N - 2} + b_{N - 1} U_{N - 1} = R_{N - 1} \end{matrix}

La dernière équation réduite donne :

(b_{N - 1} - a_{N - 1} γ_{N - 2}) U_{N - 1} = R_{N - 1} - a_{N - 2} x_{N - 2}

c'est-à-dire :

U_{N - 1} = x_{N - 1}

Il reste à calculer les $U_{k}$ en partant du dernier par la relation :

\begin{matrix} U_{k} = x_{k} - γ_{k} U_{k + 1} & (5) \end{matrix}

Les coefficients des diagonales sont stockés dans trois tableaux (à N éléments) a,b et c dès que les conditions limites et les pas sont fixés. Les tableaux β et γ (relations 1 et 2) sont calculés par récurrence avant le départ de la boucle d'itération. À chaque pas de l'itération (à chaque instant), on calcule par récurrence la suite $x_{k}$ (relation 3) pour k variant de 0 à N-1, et enfin la suite $U_{k}$ (relation 4) pour k variant de N-1 à 0. En pratique, dans cette dernière boucle, on écrit directement dans le tableau utilisé pour stocker les $U_{j, n}$ .

Références

[1] J.W. Thomas, Numerical partial differential equations, (Springer-Verlag, 2010)

[2] J.H. Ferziger, M. Peric, Computational methods for fluid dynamics, (Springer, 2002)

[3] R.H. Pletcher, J.C. Tannehill, D.A. Anderson, Computational Fluid Mechanics and Heat Transfer, (CRC Press, 2013)