Équation de diffusion à une dimension avec Python

1. Fonctions de calcul

La méthode numérique est expliquée dans Équation de diffusion à une dimension. Dans l'implémentation détaillée ci-dessous, le schéma de Cranck-Nicolson est utilisé, 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 fonctions sont incluses dans une classe. Le constructeur __init__ crée les différents tableaux. Son unique argument est le nombre N de points sur le domaine [0,1]. La fonction config définit le schéma numérique avec les conditions limites. Les tableaux a,b,c,d,e,f et s sont remplis dans cette fonction. Les tableaux β et γ utilisés dans l'élimination de Gauss-Jordan sont aussi remplis.

Les arguments de la fonction config sont :

dt : pas de temps
coef : coefficients de diffusion et frontières sous la forme [[x1,D1],[x2,D2],...[1,Dk]]
type0 : type de la condition limite en x=0, "neumann" ou "dirichlet"
lim0 : valeur de la condition limite en x=0
type1 : type de la condition limite en x=1
lim1 : valeur de la condition limite en x=1
S : source, tableau numpy à N éléments

Diffusion1D.py

import numpy

class Diffusion1D:
    def __init__(self,N):
        self.N = N
        self.a = numpy.zeros(N)
        self.b = numpy.zeros(N)
        self.c = numpy.zeros(N)
        self.d = numpy.zeros(N)
        self.e = numpy.zeros(N)
        self.f = numpy.zeros(N)
        self.s = numpy.zeros(N)
        self.alpha = numpy.zeros(N)
        self.beta = numpy.zeros(N)
        self.gamma = numpy.zeros(N)
        self.x = numpy.zeros(N)
    def config(self,dt,coef,type0,lim0,type1,lim1,S):
        N = self.N
        dx = 1.0/(N-1)
        self.dt = dt
        v = dt/(dx*dx)
        i1 = 0
        for k in range(len(coef)): # coefficients de diffusion
            x = coef[k][0]
            diff = coef[k][1]
            i2 = int(x*N)
            for i in range(i1,i2):
                self.alpha[i] = diff*v
            i1 = i2
        for i in range(1,N-1): # Cranck-Nicolson
            self.a[i] = -self.alpha[i]/2
            self.b[i] = (1+self.alpha[i])
            self.c[i] = self.a[i]
            self.d[i] = self.alpha[i]/2
            self.e[i] = 1-self.alpha[i]
            self.f[i] = self.alpha[i]/2
            self.s[i] = dt*S[i]
        if len(coef)>1: # frontieres entre deux milieux differents
            for k in range(len(coef)-1):
                i = int(coef[k][0]*N)
                self.a[i] = -coef[k][1]
                self.b[i] = coef[k][1]+coef[k+1][1]
                self.c[i] = -coef[k+1][1]
                self.d[i] = 0
                self.e[i] = 0
                self.f[i] = 0
                self.s[i] = 0
        if type0=="dirichlet":
            self.b[0] = 1.0
            self.c[0] = 0.0
            self.e[0] = 0
            self.f[0] = 0
            self.s[0] = lim0
        if type1=="dirichlet":
            self.a[N-1] = 0
            self.b[N-1] = 1.0
            self.d[N-1] = 0
            self.e[N-1] = 0
            self.s[N-1] = lim1
        if type0=="neumann":
            self.b[0] = 1.0+self.alpha[0]
            self.c[0] = -self.alpha[0]
            self.e[0] = 1.0-self.alpha[0]
            self.f[0] = self.alpha[0]
            self.s[0] = -2*self.alpha[0]*lim0*dx+dt*S[0]
        if type1=="neumann":
            self.a[N-1] = -self.alpha[N-1]
            self.b[N-1] = 1.0+self.alpha[N-1]
            self.d[N-1] = self.alpha[N-1]
            self.e[N-1] = 1.0-self.alpha[N-1]
            self.s[N-1] = 2*self.alpha[N-1]*dx*lim1+dt*S[N-1]
        
        self.beta[0] = self.b[0]
        self.gamma[0] = self.c[0]/self.beta[0]
        for k in range(1,N):
            self.beta[k] = self.b[k]-self.a[k]*self.gamma[k-1]
            if self.beta[k]==0:
                raise Exception("Impossible de resoudre le systeme ligne "+str(k))
            self.gamma[k] = self.c[k]/self.beta[k]

La fonction config sera généralement appelée plusieurs fois au cours d'une simulation, lorsque le pas de temps doit être changé. Les phénomènes de diffusion sont en effet très rapides au début (lorsque les gradients sont grands) et beaucoup plus lents à l'approche du régime stationnaire. La présente implémentation n'effectue pas d'adaptation du pas de temps; c'est donc à l'utilisateur de le faire.

La fonction iterations implémente le schéma de Cranck-Nicolson. Ses arguments sont

U0 : tableau numpy à N éléments
ti : temps initial
tf : temps final

La fonction renvoit le tableau U et le dernier instant calculé.

    def iterations(self,U0,ti,tf):
        if U0.size!=self.N:
            raise Exception("Taille de U incorrecte")
        U = U0.copy()
        t = ti
        while t<tf:
            self.x[0] = (self.e[0]*U[0]+self.f[0]*U[1]+self.s[0])/self.beta[0]
            for k in range(1,self.N-1):
                self.x[k] = (self.d[k]*U[k-1]+self.e[k]*U[k]+self.f[k]*U[k+1]+self.s[k]-self.a[k]*self.x[k-1])/self.beta[k]
            k = self.N-1
            self.x[k] = (self.d[k]*U[k-1]+self.e[k]*U[k]+self.s[k]-self.a[k]*self.x[k-1])/self.beta[k]
            U[self.N-1] = self.x[self.N-1]
            for k in range(self.N-2,-1,-1):
                U[k] = self.x[k]-self.gamma[k]*U[k+1]
            t += self.dt
        return [U,t]

2. Exemple : diffusion thermique dans une plaque

2.a. Température fixée aux bords

On considère une plaque (perpendiculaire à l'axe x) de conductivité thermique uniforme, soumise en x=0 à une température constante U=1 et en x=1 à une température constante U=0. Il n'y a aucune source thermique dans la plaque. Initialement la température est nulle sur l'intervalle [0,1].

Plusieurs groupes d'itérations sont effectuées, avec à chaque fois une augmentation d'un facteur 10 du pas de temps. Cela permet de mettre en évidence à la fois les variations de la température à court terme et à long terme.

import numpy
from Diffusion1D import Diffusion1D
from matplotlib.pyplot import *
               
N=100
U=numpy.zeros(N)
S=numpy.zeros(N)
x=numpy.arange(N)*1.0/N
coef = [[1,1]]
t=0
diffusion = Diffusion1D(N)
diffusion.config(0.0001,coef,"dirichlet",1,"dirichlet",0,S)
[U1,t]=diffusion.iterations(U,t,0.001)
diffusion.config(0.001,coef,"dirichlet",1,"dirichlet",0,S)
[U2,t]=diffusion.iterations(U1,t,0.01)
diffusion.config(0.01,coef,"dirichlet",1,"dirichlet",0,S)
[U3,t]=diffusion.iterations(U2,t,0.1)
diffusion.config(0.1,coef,"dirichlet",1,"dirichlet",0,S)
[U4,t]=diffusion.iterations(U3,t,1)
figure(figsize=(8,6))
plot(x,U1,label="t=0.001")
plot(x,U2,label="t=0.01")
plot(x,U3,label="t=0.1")
plot(x,U4,label="t=1")
legend(loc="upper right")
xlabel("x")
ylabel("T")
axis([0,1,0,1])
grid()

figA.pdf

2.b. Système à trois plaques

On considère un système isolé formé de deux plaques initialement à deux températures différentes, mises en contact thermique par une troisième plaque mince de conductivité plus faible.

N=100
U=numpy.zeros(N)
S=numpy.zeros(N)
x=numpy.arange(N)*1.0/N
coef = [[0.45,1],[0.55,0.05],[1,1]];
for i in range(int(N/2)):
    U[i] = 1
t=0
diffusion = Diffusion1D(N)
diffusion.config(0.0001,coef,"neumann",0,"neumann",0,S)
[U1,t]=diffusion.iterations(U,t,0.001)
diffusion.config(0.001,coef,"neumann",0,"neumann",0,S)
[U2,t]=diffusion.iterations(U1,t,0.01)
diffusion.config(0.01,coef,"neumann",0,"neumann",0,S)
[U3,t]=diffusion.iterations(U2,t,0.1)
diffusion.config(0.1,coef,"neumann",0,"neumann",0,S)
[U4,t]=diffusion.iterations(U3,t,1)
[U5,t]=diffusion.iterations(U4,t,4)
figure(figsize=(8,6))
plot(x,U1,label="t=0.001")
plot(x,U2,label="t=0.01")
plot(x,U3,label="t=0.1")
plot(x,U4,label="t=1")
plot(x,U5,label="t=4")
legend(loc="upper right")
xlabel("x")
ylabel("T")
axis([0,1,0,1])
grid()

figB.pdf

3. Exemple : diffusion de particules

On considère la diffusion de particules dans un récipient. Initialement, les particules sont concentrées au milieu du domaine (sur un dizième de la largeur). Les conditions limites sont de type Neumann avec un flux nul.

N=100
U=numpy.zeros(N)
S=numpy.zeros(N)
x=numpy.arange(N)*1.0/N
coef = [[1,1]]
t=0
i1 = int(0.45*N)
i2 = int(0.55*N)
for i in range(i1,i2):
    U[i] = 1.0
diffusion = Diffusion1D(N)
diffusion.config(0.00001,coef,"neumann",0,"neumann",0,S)
[U1,t]=diffusion.iterations(U,t,0.0001)
diffusion.config(0.0001,coef,"neumann",0,"neumann",0,S)
[U2,t]=diffusion.iterations(U1,t,0.001)
diffusion.config(0.001,coef,"neumann",0,"neumann",0,S)
[U3,t]=diffusion.iterations(U2,t,0.01)
diffusion.config(0.01,coef,"neumann",0,"neumann",0,S)
[U4,t]=diffusion.iterations(U3,t,0.1)
figure(figsize=(8,6))
plot(x,U1,label="t=0.0001")
plot(x,U2,label="t=0.001")
plot(x,U3,label="t=0.01")
plot(x,U4,label="t=0.1")
legend(loc="upper right")
xlabel("x")
ylabel("T")
axis([0,1,0,1])
grid()

figC.pdf

Reprenons la première étape du calcul avec un pas d'espace Δx dix fois plus petit :

N=1000
U=numpy.zeros(N)
S=numpy.zeros(N)
x=numpy.arange(N)*1.0/N
coef = [[1,1]]
t=0
i1 = int(0.45*N)
i2 = int(0.55*N)
for i in range(i1,i2):
    U[i] = 1.0
diffusion = Diffusion1D(N)
diffusion.config(0.00001,coef,"neumann",0,"neumann",0,S)
[U1,t]=diffusion.iterations(U,t,0.0001)
figure(figsize=(8,6))
plot(x,U1,label="t=0.0001")
legend(loc="upper right")
xlabel("x")
ylabel("T")
axis([0.4,0.6,0,1])
grid()

figD.pdf

On voit apparaître un artefact au voisinage de la discontinuité initiale de concentration. Réduisons le pas de temps Δt d'un facteur 5 :

N=1000
U=numpy.zeros(N)
S=numpy.zeros(N)
x=numpy.arange(N)*1.0/N
coef = [[1,1]]
t=0
i1 = int(0.45*N)
i2 = int(0.55*N)
for i in range(i1,i2):
    U[i] = 1.0
diffusion = Diffusion1D(N)
diffusion.config(0.000002,coef,"neumann",0,"neumann",0,S)
[U1,t]=diffusion.iterations(U,t,0.0001)
figure(figsize=(8,6))
plot(x,U1,label="t=0.0001")
legend(loc="upper right")
xlabel("x")
ylabel("T")
axis([0.4,0.6,0,1])
grid()

figE.pdf

Le problème a disparu. Bien que le schéma de Cranck-Nicolson soit toujours stable (quelque soient Δt et Δx), il peut être nécessaire de réduire le pas de temps lorsque le pas d'espace est très petit.