Table des matières Python

Équation de diffusion à deux dimensions

1. Équation de diffusion

On s'intéresse à l'équation de diffusion (ou équation de la chaleur) à deux dimensions en coordonnées cartésiennes. Soit u(x,y,t) la fonction vérifiant l'équation de diffusion suivante :

ut=D(2ux2+2uy2)+s(x,y,t)(1)

D est le coefficient de diffusion. La fonction s(x,y,t) est la source.

Les applications physiques de cette équations sont nombreuses. On peut citer :

On remarque qu'une solution stationnaire (indépendante du temps) est solution de l'équation de Poisson.

La discrétisation directe de l'équation de diffusion se fait avec la méthode des différences finies. Une autre approche consiste à mettre l'équation sous forme intégrale, c'est-à-dire sous la forme d'une loi de conservation. Considérons pour cela un volume V délimité par une surface fermée S. L'intégrale de l'équation de diffusion sur ce volume s'écrit :

Vutdv=VDdiv(gradu)dv+Vsdv(2)

Le théorème d'Ostrogradsky permet de faire apparaitre le flux sortant du gradient sur la surface fermée :

Vutdv=SD(gradu)ndS+Vsdv(3)

On se limite ici à l'équation de diffusion bidimensionnelle en coordonnées cartésiennes. Dans ce cas, le gradient est dans le plan (x,y) et l'intégrale de surface se ramène à une intégrale sur une courbe fermée C; l'intégrale de volume se ramène à une intégrale sur la surface Σ délimitée par la courbe :

ΣutdΣ=CD(gradu)ndl+ΣsdΣ(4)

Cette équation permet de discrétiser par la méthode des volumes finis.

2. Maillage et différences finies

On considère un domaine rectangulaire. Les nœuds du maillage sont définis par :

xi=iΔx(5)yj=jΔy(6)

L'indice i varie de 0 à Nx-1, l'indice j de 0 à Ny-1.

Les valeurs sont recherchées aux instants définis par :

tn=nΔt(7)

n est un entier positif ou nul. On pose :

Ui,jn=u(iΔx,jΔy,nΔt)(8)Si,jn=s(iΔx,jΔy,nΔt)(9)

Dans cette notation, l'exposant sert donc à représenter le temps.

figureA.svgFigure pleine page

On s'intéresse à un point A situé à l'intérieur du domaine (pas sur un bord). La discrétisation de la dérivée seconde par rapport à x se fait par :

2ux2(xi,yi,tn)Ui+1,jn-2Ui,jn+Ui-1,jn(Δx)2(10)

La dérivée temporelle se discrétise de la manière suivante :

ut(xi,yj,tn)Ui,jn+1-Ui,jnΔt(11)

Voici le schéma explicite déduit de ces expressions :

Ui,jn+1-Ui,jnΔt=DUi+1,jn-2Ui,jn+Ui-1,jn(Δx)2+DUi,j+1n-2Ui,jn+Ui,j-1n(Δy)2+Si,jn(12)

Ce schéma permet de calculer explicitement Ui,jn+1 en fonction des valeurs à l'instant n au point considéré et sur les 4 voisins. Il est intéressant de retrouver cette équation par la méthode des volumes finis, qui consiste à discrétiser les intégrales de l'équation (4) exprimée sur le contour (a,b,c,d) qui entoure le point A (voir figure ci-dessus) :

ΣutdΣUi,jn+1-Ui,jnΔtΔxΔy(13)abD(gradu)ndlUi,j+1n-Ui,jnΔyΔx(14)bcD(gradu)ndlUi-1,jn-Ui,jnΔxΔy(15)cdD(gradu)ndlUi,j-1n-Ui,jnΔyΔx(16)daD(gradu)ndlUi+1,jn-Ui,jnΔxΔy(17)ΣsdΣSi,jnΔxΔy(18)

En reportant ces différents termes dans la forme intégrale, on retrouve bien l'équation (12). La méthode des volumes finis est très utile pour la discrétisation des conditions limites de Neumann. Elle permet aussi de discrétiser l'équation en coordonnées non cartésiennes.

3. Schéma implicite

Comme expliqué dans Équation de diffusion à une dimension, le schéma explicite comporte une condition de stabilité très contraignante dès que le pas d'espace est très petit. On utilise donc un schéma implicite, comme le suivant :

Ui,jn+1-Ui,jnΔt=DUi+1,jn+1-2Ui,jn+1+Ui-1,jn+1(Δx)2+DUi,j+1n+1-2Ui,jn+1+Ui,j-1n+1(Δy)2+Si,jn(19)

Le schéma de Cranck-Nicolson est une variante consistant à utiliser la moyenne des dérivées secondes spatiales à l'instant n et à l'instant n+1.

Pour ce schéma implicite, chaque pas de temps nécessite la résolution d'un système linéaire à 5 inconnues par équation :

Δy2(Ui+1,jn+1+Ui-1,jn+1)+Δx2(Ui,j+1n+1+Ui,j-1n+1)-(2(Δx2+Δy2)+DΔx2Δy2DΔt)Ui,jn+1=-Δx2Δy2DΔtUi,jn-1DSi,jn(20)

Dans le cas à une dimension, le système comporte 3 inconnues par équation, ce qui permet d'appliquer un algorithme très efficace d'élimination pour une matrice tridiagonale. Dans le cas bidimensionnel, on choisira plutôt une méthode itérative, comme la méthode de Gauss-Seidel multigrille. On remarque d'ailleurs que le schéma est très proche de celui utilisé pour la résolution de l'équation de Poisson. Pour le schéma de Cranck-Nicolson, il faut modifier le terme de droite.

4. Schéma implicite à directions alternées

4.a. Définition

Le schéma implicite à directions alternées consiste à décomposer le pas de temps (n,n+1) en deux sous-pas. Pour le premier, on écrit un schéma implicite pour la dérivée par rapport à x mais explicite pour la dérivée par rapport à y :

Ui,jn+1/2-Ui,jn(Δt/2)=DUi+1,jn+1/2-2Ui,jn+1/2+Ui-1,jn+1/2(Δx)2+DUi,j+1n-2Ui,jn+Ui,j-1n(Δy)2+Si,jn(21)

Le second est implicite dans la direction y :

Ui,jn+1-Ui,jn+1/2(Δt/2)=DUi+1,jn+1/2-2Ui,jn+1/2+Ui-1,jn+1/2(Δx)2+DUi,j+1n+1-2Ui,jn+1+Ui,j-1n+1(Δy)2+Si,jn(22)

Chacune de ces deux étapes consiste à résoudre un système tridiagonal, ce qui se fait par élimination directe en utilisant l'algorithme décrit dans équation de diffusion à une dimension.

On écrira les deux systèmes à résoudre sous la forme générale suivante :

Ai,jUi-1,jn+1/2+Bi,jUi,jn+1/2+Ci,jUi+1,jn+1/2=Di,jUi,j-1n+Ei,jUi,jn+Fi,jUi,j+1n+Gi,j(23)A'i,jUi,j-1n+1+B'i,jUi,jn+1+C'i,jUi,j+1n+1=D'i,jUi-1,jn+1/2+E'i,jUi,jn+1/2+F'i,jUi+1,jn+1/2+G'i,j(24)

Le premier système tridiagonal est à résoudre pour chaque ligne. Le second est à résoudre pour chaque colonne.

Pour le point A à l'intérieur du domaine, on a donc pour le premier système :

Ai,j=DΔx2(25)Bi,j=-2DΔx2-2Δt(26)Ci,j=DΔx2(27)Di,j=-DΔy2(28)Ei,j=2DΔy2-2Δt(29)Fi,j=-DΔy2(30)Gi,j=-Si,jn(31)

et pour le second :

A'i,j=DΔy2(32)B'i,j=-2DΔy2-2Δt(33)C'i,j=DΔy2(34)D'i,j=-DΔx2(35)E'i,j=2DΔx2-2Δt(36)F'i,j=-DΔx2(37)G'i,j=-Si,jn+1/2(38)

En pratique, on multipliera tous ces coefficients par Δt et on simplifiera en posant :

αx=DΔtΔx2(39)αy=DΔtΔy2(40)

Dans la plupart des cas, on prendra Δx=Δy.

4.b. Conditions limites

Une condition limite de Dirichlet sur le point (i,j) s'écrit :

Ai,j=0(41)Bi,j=1(42)Ci,j=0(43)Di,j=0(44)Ei,j=0(45)Fi,j=0(46)Gi,j=U0(47)

Considérons une condition limite de Neumann en un point du bord gauche (point D). Cette condition consiste à fixer la dérivée suivante :

ux=wx(48)

Pour discrétiser cette condition limite, on suit le même principe que pour l'équation de diffusion à une dimension. On commence par écrire une discrétisation d'ordre 2 en faisant intervenir un point fictif situé à gauche de la frontière :

U1,jn-U-1,jn2Δx=wx(49)

La discrétisation de la dérivée seconde par rapport à x s'écrit :

2ux2U1,jn-2U0,jn+U-1,jn(Δx)2(50)

Le point fictif i=-1 est éliminé avec la relation (49) :

2ux22U1,jn-2U0,jn-2Δxwx(Δx)2(51)

Pour le schéma implicite à directions alternées, le premier demi-pas, correspondant à la direction x, est défini par :

Ai,j=0(52)Bi,j=-2DΔx2-2Δt(53)Ci,j=2DΔx2(54)Di,j=-DΔy2(55)Ei,j=2DΔy2-2Δt(56)Fi,j=-DΔy2(57)Gi,j=2DwxΔx-Si,jn(58)

Le seconde demi-pas, correspondant à la direction y, est défini par :

A'i,j=DΔy2(59)B'i,j=-2DΔy2-2Δt(60)C'i,j=DΔy2(61)D'i,j=0(62)E'i,j=2DΔx2-2Δt(63)F'i,j=-2DΔx2(64)G'i,j=2DwxΔx-Si,jn+1/2(65)

Voyons comment retrouver ce résultat avec la méthode des volumes finis, qui sera plus simple pour traiter les cas plus complexes. Le contour (o,p,q,r) est défini sur la figure. Voici comment sont discrétisées les différentes intégrales :

ΣutdΣUi,jn+1-Ui,jnΔt12ΔxΔy(66)opD(gradu)ndlDUi,j+1n-Ui,jnΔy12Δx(67)pqD(gradu)ndl-DwxΔy(68)qrD(gradu)ndlDUi,j-1n-Ui,jnΔy12Δx(69)roD(gradu)ndlDUi+1,jn-Ui,jnΔxΔy(70)ΣsdΣSi,jn12ΔxΔy(71)

En écrivant le schéma implicite à directions alternées à partir de ces équations, on retouve bien les coefficients précédents (attention à la division de Δt par deux dans cette opération).

Pour les bords du domaine rectangulaire, il y a 8 types de conditions limites de Neumann (4 côtés et 4 coins). Pour un objet situé à l'intérieur du domaine (contour GHIJKL), il a 4 autres types de coins. Il y a donc en tout 12 types de condition limite de Neumann.

Pour une condition limite de Neumann, il faut une convention pour fixer l'intérieur du domaine de calcul par rapport au sens de parcours d'un contour. On adopte la convention suivante : l'intérieur du domaine est situé à gauche si l'on s'oriente dans le sens du parcours. Par exemple, le contour rectangulaire délimitant le domaine doit être parcouru dans le sens trigonométrique. Le contour de l'objet situé à l'intérieur du domaine devra être parcouru dans le sens horaire. Compte tenu de cette convention, on peut donner une nomenclature simple pour désigner les différentes arêtes et coins. Une arête ou un coin sont définis par deux segments, chacun d'eux défini par deux incréments di et dj égaux à -1, 0 ou +1. Par exemple, le coin inférieur gauche du domaine est défini par (0,-1)-(1,0). Une arête situé sur le bord inférieur est définie par (1,0)-(1,0). Le coin G de l'onjet situé dans le domaine est défini par (-1,0)-(0,1).

À titre d'exemple, voici le résultat pour le coin inférieur gauche :

Ai,j=0(72)Bi,j=-4DΔx2-2Δt(73)Ci,j=4DΔx2(74)Di,j=0(75)Ei,j=2DΔy2-2Δt(76)Fi,j=-2DΔy2(77)Gi,j=-Si,jn+2DwxΔx+2DwyΔy(78)A'i,j=0(79)B'i,j=-2DΔy2-2Δt(80)C'i,j=2DΔy2(81)D'i,j=0(82)E'i,j=2DΔx2-2Δt(83)F'i,j=-4DΔx2(84)G'i,j=-Si,jn+1/2+2DwxΔx+2DwyΔy(85)

4.c. Algorithme

La maillage est constitué de matrices à Nx colonnes et Ny lignes. La définition des pas de temps et d'espace, des conditions limites et des sources permet de construire les matrices A,B,C,D,E,F,G,A',B',C',D',E',F',G'.

Il faut aussi construire les matrices β,γ pour la résolution des systèmes sur les lignes, et les matrices β',γ' pour la résolution des systèmes sur les colonnes. Pour chaque ligne j, la récurrence à appliquer pour i variant de 1 à Nx :

β1,j=B1,j(86)γ1,j=C1,jβ1,j(87)βi,j=Bi,j-Ai,jγi-1,j(88)γi,j=Ci,jβi,j(89)

Pour chaque colonne i, la récurrence à appliquer pour j variant de 1 à Ny :

β'i,1=B'i,1(90)γ'i,1=C'i,1β'i,1(91)β'i,j=B'i,j-A'i,jγ'i,j-1(92)γ'i,j=C'i,jβ'i,j(93)

Pour chaque intervalle de temps Δt, on doit effectuer une résolution tridiagonale pour chaque ligne, puis pour chaque colonne.

Pour chaque ligne j on commence par calculer la suite suivante pour i variant de 1 à Nx:

x1=1β1,j(D1,jU1,j-1n+E1,jU1,jn+F1,jU1,j+1n+G1,j)(94)xi=1βi,j(Di,jUi,j-1n+Ei,jUi,jn+Fi,jUi,j+1n+Gi,j-Ai,jxi-1)(95)

Pour la même ligne j, on calcule les nouvelles valeurs à l'instant n+1/2 par récurrence, en faisant varier de i de Nx à 1 :

UNx,jn+12=xNx(96)Ui,jn+12=xi-γi,jUi+1,jn+12(97)

Les valeurs à l'instant n+1/2 sont stockées dans la matrice U. On procède alors au parcours des colonnes. Pour chaque colonne i, on calcule la suite suivante pour j variant de 1 à Ny :

x1=1β'i,1(D'i,1Ui-1,1n+1/2+E'i,1Ui,1n+1/2+F'i,1Ui+1,1n+1/2+G'i,1)(98)xi=1β'i,j(D'i,jUi-1,jn+1/2+E'i,jUi,jn+1/2+F'i,jUi+1,jn+1/2+G'i,j-A'i,jxi-1)(99)

puis pour j variant de Ny à 1 :

Ui,Nyn=xNy(100)Ui,jn=xi-γ'i,jUi,j+1n(101)

Lorsqu'on traite les lignes, il faut écrire dans une seconde matrice U, pour ne pas écraser les valeurs de l'instant n. Il en est de même lorsqu'on traite les colonnes.

Le nombre de matrices utilisées est 20. Pour un maillage 1024x1024 en flottants 32 bits, la mémoire utilisée est de 80 méga-octets.

Les lignes peuvent être parcourues en parallèle (idem pour les colonnes). On peut donc avec profit utiliser un processeur graphique pour faire ces calculs (via l'API CUDA ou OpenCL). Pour le traitement des lignes, chaque unité de calcul traite une ligne. Les opérations effectuées sur les différentes lignes sont rigoureusement identiques, ce qui est le cas idéal pour une implémentation sur GPU. Bien entendu, toutes les matrices doivent être transférées sur la mémoire de la carte graphique avant le calcul.

On remarque que le traitement des bords (qui doit être identique au traitement des autres points) nécessite un accès à des éléments de U hors domaine. Pour résoudre ce problème, on utilise des matrices avec des rangées fictives sur les bords. Le nombre d'éléments dans chaque direction est une puissance de 2 (plus efficace pour l'implémentation GPU). La figure suivante montre un exemple pour Nx=23 et Ny=23 (les points sont les nœuds de la grille).

matrice.svgFigure pleine page

Sur cet exemple, les matrices sont 10x10 mais le maillage effectif est 8x8.

Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.