Table des matières Python

Discrétisation de l'équation de Poisson en coordonnées polaires

1. Introduction

On considère l'équation de Poisson :

2u=s(1)

s est un champ scalaire appelé la source et u le champ scalaire à déterminer, qui peut être un potentiel électrostatique, la température pour un problème de transfert thermique, le potentiel des vitesses pour un écoulement, etc.

On se place dans le cas d'un problème bidimensionnel en coordonnées polaires. L'équation de Poisson s'écrit alors :

1rrrur+1r22uθ2=s(r,θ)(2)

La méthode de discrétisation des volumes finis ([1][2]) consiste à utiliser l'équation (1) intégrée sur un volume. Avec le théorème d'Ostrogradsky, on obtient pour un volume V délimitée par une surface S :

Sgradunds=V sdv(3)

Physiquement, il s'agit d'une équation de conservation exprimée pour le volume de contrôle V. Par exemple pour un problème de transfert thermique, cette équation exprime la conservation de l'énergie.

Il est important que l'équation de conservation soit vérifiée par le schéma numérique. La méthode des volumes finis consiste à discrétiser ces deux intégrales pour des volumes de contrôle choisis en fonction du système de coordonnées utilisé.

2. Définition du maillage

Le maillage est défini par Nr cercles espacés de Δr et Nθ droites espacés de Δθ. Les rayons et angles des nœuds du maillage sont :

θi=iΔθ(4)rj=jΔr(5)

Le premier cercle (j=0) est réduit au point O. On a par ailleurs :

NθΔθ=2π(6)

On définit:

Ui,j=u(θi,rj)(7)Si,j=s(θi,rj)(8)

La périodicité de l'angle s'écrit :

UNθ,j=U0,j(9) figureA.svgFigure pleine page

3. Discrétisation de l'équation

Soit (i,j) un point à l'intérieur du domaine. Le volume de contrôle est représenté en pointillé sur la figure sous forme d'un contour (a,b,c,d). Dans la direction perpendiculaire au plan, la longueur du volume est constante, égale à 1. Le volume est donc :

Δv=Δθ2(rj+Δr2)2-(rj-Δr2)2=jΔθΔr2(10)

La méthode des volumes finis consiste à remplacer les deux intégrales de l'équation (3) par des approximations en fonction des valeurs de U aux nœuds du maillage. Pour le terme de source, l'approximation la plus simple est :

vsdvSi,jΔv(11)

Le flux du gradient se ramène ici à une intégrale sur le contour (a,b,c,d). Si on considère par exemple l'arc (a,b), on peut utiliser l'approximation suivante :

(ab)gradundlUi,j-1-Ui,jΔrrj-12Δθ=(Ui,j-1-Ui,j)(j-12)Δθ(12)

L'approximation pour l'arc (c,d) est similaire. Pour le segment (b,c), on utilise l'approximation suivante :

(bc)gradundlUi+1,j-Ui,jrjΔθΔr=Ui+1,j-Ui,jjΔθ(13)

L'équation discrétisée se met sous la forme générale suivante :

Di,jUi+1,j+Gi,jUi-1,j+Hi,jUi,j+1+Bi,jUi,j-1+Ci,jUi,j=Fi,j(14)

On obtient finalement pour un point à l'intérieur du domaine :

Di,j=Gi,j=1(15)Hi,j=(j+12)jΔθ2(16)Bi,j=(j-12)jΔθ2(17)Ci,j=-Di,j-Gi,j-Hi,j-Bi,j(18)Fi,j=Si,jj2Δr2Δθ2(19)

Lorsque le centre O se trouve dans le domaine d'intégration, il nécessite un traitement particulier. On note U0 la valeur de u au point O. Le volume de contrôle est le cercle tracé sur la figure. La discrétisation de la loi de conservation s'obtient en écrivant le gradient sur chaque arc de cercle :

i=0Nθ-1Ui,1-U0Δr12ΔrΔθ=S0πΔr22(20)

Cette équation ne correspond pas au schéma (14) car elle implique tous les points du premier cercle.

4. Discrétisation des conditions limites

4.a. Conditions de Dirichlet

Une condition limite de Dirichlet consiste à imposer une valeur de u sur un bord du domaine ou d'un objet. Par exemple, si une condition de Dirichlet est définie sur le bord du domaine par une fonction g(θ), on écrira pour chaque point du bord :

Ui,Nr-1=g(iΔθ)(21)

Exprimé avec le schéma général, la condition s'écrit :

Ci,Nr-1=1(22)Fi,Nr-1=g(iΔθ)(23)

4.b. Condition de Neumann

La condition limite de Neumann consiste à fixer le gradient de u sur une frontière, plus précisément la composante du gradient perpendiculaire à la frontière. Physiquement, cela correspond à un flux imposé (champ électrique imposé, flux thermique imposé, etc).

On considère tout d'abord le bord du domaine, constitué du cercle j=Nr-1. Pour un point du bord, le volume de contrôle (e,f,g,h) est représenté en pointillé sur la figure. le volume est :

Δθ2rj2-(rj-Δr2)2=12(j-14)ΔθΔr2(24)

La dérivée imposée est par convention la dérivée par rapport à r (flux sortant du domaine) :

ur=wr(θ)(25)

On définit donc :

wr,i=w(iΔθ)(26)

Le flux sortant du volume de contrôle par le côté (g,h) est discrétisé de la manière suivante :

(gh)gradundlwr,ijΔrΔθ(27)

Pour le côté (e,f) on écrit :

(ef)gradundlUi,j-1-Ui,jΔrrj-12Δθ=(Ui,j-1-Ui,j)(j-12)Δθ(28)

Le côté (fg) peut être traité comme le côté (bc) :

(fg)gradundlUi+1,j-Ui,jrjΔθΔr2=Ui+1,j-Ui,j2jΔθ(29)

On obtient finalement pour la condition limite de Neumann sur le bord du domaine (j=Nr-1) :

Di,j=Gi,j=12(30)Hi,j=0(31)Bi,j=(j-12)jΔθ2(32)Ci,j=-Di,j-Gi,j-Hi,j-Bi,j(33)Fi,j=Si,j12(j-14)jΔθ2Δr2-wr,ij2ΔrΔθ2(34)

Il peut arriver que le domaine soit limité par un cercle interne. Dans ce cas, le premier indice est j=J0 et il faut écrire une condition limite sur ce cercle. On obtient alors :

Di,j=Gi,j=12(35)Hi,j=(j+12)jΔθ2(36)Bi,j=0(37)Ci,j=-Di,j-Gi,j-Hi,j-Bi,j(38)Fi,j=Si,j12(j+14)jΔθ2Δr2+wr,ij2ΔrΔθ2(39)

Voyons le cas d'un objet dans le domaine (contour PQRS) avec des conditions limites de Neumann. Pour les points des arcs PS et QR, la discrétisation est similaire à celle faite ci-dessus pour les bords du domaine. Sur le bord PQ, il faut faire intervenir la dérivée suivante :

uθ=wθ(40)

On pose donc pour ce côté :

wθ,j=wθ(jΔr)(41)

Le volume de contrôle est représenté sur la figure par le contour (lmno). Le volume est

Δv=12jΔθΔr2(42)

L'approximation du flux sortant du volume de contrôle par le côté (no) est :

(no)gradundlwθ,jΔr(43)

Le côté (lm) est identique au côté (ad) traité plus haut. Le côté (nm) se traite comme le côté (cd) avec une longueur d'arc réduite de moitié. On obtient finalement :

Di,j=0(44)Gi,j=1(45)Hi,j=(j+12)jΔΘ22(46)Bi,j=(j-12)jΔθ22(47)Ci,j=-Di,j-Gi,j-Hi,j-Bi,j(48)Fi,j=Si,j12j2Δθ2Δr2-wθ,jjΔrΔθ(49)

Les coins doivent recevoir un traitement particulier. Par exemple pour le coin R, le volume de contrôle est délimité par le contour (uvwxyz). Le côté (yz) est soumis à la condition de Neumann orthoradiale alors que le côté (xy) est soumis à une condition radiale. Cela signifie qu'il y a deux flux à spécifier pour les coins. On obtient pour ce point le volume :

Δv=34j+116ΔθΔr2(50)

et les coefficients suivants :

Di,j=1(51)Gi,j=12(52)Hi,j=(j+12)jΔθ2(53)Bi,j=(j-12)jΔθ22(54)Ci,j=-Di,j-Gi,j-Hi,j-Bi,j(55)Fi,j=Si,j(34j+116)jΔr2Δθ2+wθ,jj2ΔrΔθ+wr,jj22ΔrΔθ2(56)

le coin V du contour TUWXY est différent du précédent car le domaine de calcul se trouve à l'intérieur du coin. Le volume de contrôle est (a,b,c,d). Le volume est :

Δv=j4-116ΔθΔr2(57)

On obtient pour ce point :

Di,j=12(58)Gi,j=1(59)Hi,j=0(60)Bi,j=(j-12)jΔθ22(61)Ci,j=-Di,j-Gi,j-Hi,j-Bi,j(62)Fi,j=Si,j(14j-116)jΔr2Δθ2+wθ,jj2ΔrΔθ-wr,jj22ΔrΔθ2(63)
Références
[1]  J.W. Thomas,  Numerical partial differential equations,  (Springer-Verlag, 2010)
[2]  R.H. Pletcher, J.C. Tannehill, D.A. Anderson,  Computational Fluid Mechanics and Heat Transfer,  (CRC Press, 2013)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.