Table des matières Python

Discrétisation de l'équation de Poisson en géométrie axiale

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.

On considère le cas d'un problème à symétrie axiale, pour lequel on utilise les coordoonées cylindriques (r,θ,z). Les champs s(z,r) et u(z,r) sont alors indépendant de θ et l'équation de Poisson s'écrit :

1rr(rur)+2uz2=s(z,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 :

Sgrad(u)nds=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 les nœuds suivants

zi=iΔz(4)rj=jΔr(5)

où l'indice i varie de 0 à Nz-1 et l'indice j de 0 à Nr-1. L'indice j=0 correspond à l'axe de symétrie (axe Oz).

On définit:

Ui,j=u(zi,rj)(6)Si,j=s(zi,rj)(7) figureA.svgFigure pleine page

3. Discrétisation de l'équation

Soit (i,j) un point à l'intérieur du domaine (point A), non situé sur l'axe (j>0). 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 du maillage, le volume est délimité par des arcs de cercles centrés sur l'axe Oz, d'angle Δθ. Les sommets du contour correspondent à des arcs de cercle de longueurs :

la=lb=(j+12)ΔrΔθ(8)lc=ld=(j-12)ΔrΔθ(9)

Les segments (da) et (bc) correspondent à des surfaces dont les aires sont :

Ada=Abc=jΔr2Δθ(10)

Les segments (ab) et (cd) correspondent à des surfaces cylindriques dont les aires sont :

Aab=laΔz=(j+12)ΔrΔzΔθ(11)Acd=lcΔz=(j-12)ΔrΔzΔθ(12)

Enfin le volume est :

Δv = AdaΔz=jΔr2ΔzΔθ(13)

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(14)

Le flux sortant du gradient sur la surface (da) est évalué de la manière suivante :

(da)gradunds=Ui+1,j-Ui,jΔzAda=Ui+1,j-Ui,jΔzjΔr2Δθ(15)

Le flux sur la surface (bc) s'écrit de manière analogue. Pour la surface (ab) :

(ab)gradunds=Ui,j+1-Ui,jΔrAab=(Ui,j+1-Ui,j)(j+12)ΔzΔθ(16)

En raison de la symétrie axiale, les flux sur les deux surfaces θ constante sont nuls.

L'équation (3) 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(17)

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

Di,j=Gi,j=jΔr2(18)Hi,j=(j+12)Δz2(19)Bi,j=(j-12)Δz2(20)Ci,j=-2j(Δr2+Δz2)(21)Fi,j=Si,jjΔr2Δz2(22)

On remarque que ces coefficients sont proportionnels à l'indice j (ou j plus ou moins 1/2). Comme l'indice d'un point dépend de la taille du maillage, on divisera ce facteur par Nr, pour qu'il soit indépendant de la taille du maillage. Cette correction est indispensable pour l'utilisation des méthodes multigrilles.

On considère à présent le point B situé sur l'axe Oz d'indices (i,0), avec 0<i<Nz-1. Le volume de contrôle est représenté sur la figure par le contour (e,f,h,h). La surface (gh) est réduite à un segment, donc le flux est nul. Les aires des trois autres surfaces utiles sont :

Ahe=Afg=18Δr2Δθ(23)Aef=12ΔrΔzΔθ(24)

Le volume est :

Δv=AheΔz=18Δr2ΔzΔθ(25)

On obtient pour ce point de l'axe (point B) :

Di,0=Gi,0=18Δr2(26)Hi,0=12Δz2(27)Bi,0=0(28)Ci,0=-14Δr2-12Δz2(29)Fi,j=Si,j18Δr2Δz2(30)

4. Discrétisation des conditions limites

4.a. Condition de Dirichlet

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

Ui,Nr-1=g(iΔz)(31)

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

Ci,Nr-1=1(32)Fi,Nr-1=g(iΔz)(33)

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).

Considérons tout d'abord le cas d'un point C situé sur le bord j=Nr-1 (pas sur un coin). La dérivée imposée sur ce bord est :

ur=wr(z)(34)

On pose donc :

wr,i=wr(iΔz)(35)

Le volume de contrôle pour ce point est représenté sur la figure par le contour (k,l,m,n). On écrit les aires des différentes faces et le volume :

Akn=Alm=12(j-14)Δr2Δθ(36)Akl=jΔrΔzΔθ(37)Amn=(j-12)ΔrΔzΔθ(38)Δv=12(j-14)Δr2ΔzΔθ(39)

Le flux sortant à travers la surface (kl) est discrétisé de la manière suivante :

(da)gradunds=wr,iAkl=wr,ijΔrΔzΔθ(40)

Les flux sur les surfaces (lm), (mn) et (nk) sont discrétisés comme précédemment et on obtient :

Di,j=Gi,j=12(j-14)Δr2(41)Hi,j=0(42)Bi,j=(j-12)Δz2(43)Ci,j=-(j-14)Δr2-(j-12)Δz2(44)Fi,j=Si,j12(j-14)Δr2Δz2-wr,ijΔrΔz2(45)

Pour un point D situé sur le bord i=0, on impose la dérivée suivante :

uz=wz(r)(46)

et on pose

wz,j=wz(jΔr)(47)

Le volume de contrôle représenté par le contour (o,p,q,r) conduit aux coefficients suivants :

Gi,j=0(48)Di,j=jΔr2(49)Hi,j=12(j+12)Δz2(50)Bi,j=12(j-12)Δz2(51)Ci,j=-j(Δr2+Δz2)(52)Fi,j=Si,j12jΔr2Δz2+wz,jjΔr2Δz(53)

Pour le point E situé sur un coin du domaine (i=0,j=Nr-1), le volume de contrôle est (s,t,u,v). Le flux sur (ts) fait intervenir la dérivée par rapport à r alors que le flux sur (tu) fait intervenir la dérivée par rapport à z. On obtient pour ce point :

Gi,j=0(54)Di,j=12(j-14)Δr2(55)Hi,j=0(56)Bi,j=12(j-12)Δz2(57)Ci,j=-12(j-14)Δr2-12(j-12)Δz2(58)Fi,j=Si,j14(j-14)Δr2Δz2-wr,i12jΔrΔz2+wz,j12(j-14)Δr2Δz(59)

Le coin situé sur l'axe (i=0,j=0) nécessite un traitement spécial en raison de l'absence de flux sur l'arête (yz). On obtient pour ce point :

Gi,j=0(60)Di,j=18Δr2(61)Hi,j=14Δz2(62)Bi,j=0(63)Ci,j=-18Δr2-14Δz2(64)Fi,j=Si,j116Δr2Δz2+wz,j18Δr2Δz(65)

Pour un objet cylindrique défini par le contour ABCDEF, la condition de Neumann s'écrit de manière différente suivant le côté. Les coins doivent recevoir un traitement particulier. Sur chaque coin, il y a deux flux imposés perpendiculaires.

Le traitement des côtés est similaire à celui des bords du domaine vu ci-dessus. Par exemple, le bord GH conduit à

Gi,j=Di,j=12(j+14)Δr2(66)Hi,j=(j+12)Δz2(67)Bi,j=0(68)Ci,j=-(j+12)Δz2-(j+14)Δr2(69)Fi,j=Si,j12(j+14)Δr2Δz2+wr,ijΔrΔz2(70)

Le coin I est similaire au coin E traité ci-dessus car le domaine de calcule se trouve à l'intérieur du coin. Pour les sommets G,H,J,K,L, le domaine se trouve à l'extérieur du coin. Par exemple, pour le sommet G, dont le volume de contrôle est (a,b,c,d,e,f), on obtient :

Gi,j=jΔr2(71)Di,j=12(j-14)Δr2(72)Hi,j=12(j+12)Δz2(73)Bi,j=(j-12)Δz2(74)Ci,j=-12(3j-12)Δz2-12(3j-14)Δr2(75)Fi,j=Si,j(12(j-14)+14(j+14))Δr2Δz2-wr,i12jΔrΔz2-wz,j12(j+14)ΔzΔr2(76)
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.