Table des matières Python

Discrétisation de l'équation de Poisson vectorielle

1. Introduction

L'équation de Poisson vectorielle est :

2A=s(1)

s est un champ vectoriel donné appelé la source, et A le champ vectoriel inconnu. Cette équation est principalement utilisée pour le calcul de champs magnétiques en régime quasi stationnaire : A est le potentiel vecteur et s=-μ0J .

On se place dans le cas d'un champ A de divergence nul (comme en magnétostatique). L'équation s'écrit alors :

-rot(rotA)=s(2)

Le théorème de Stokes permet d'obtenir une forme intégrale de cette équation, sur une courbe fermée C délimitant une surface Σ :

C(rotA)tdl=-Σsdσ(3)

Il s'agit du théorème d'Ampère utilisé couramment en magnétostatique.

Pour un problème bidimensionnel en coordonnées cartésiennes, les deux vecteurs ont seulement une composante sur un axe z et l'équation devient :

2Azx2+2Azy2=sz(4)

Il s'agit de l'équation de Poisson pour un champ scalaire, dont la discrétisation est décrite dans Équation de Poisson à deux dimensions. Néanmoins, la loi de conservation ne s'exprime pas sous forme d'un flux à travers une surface fermée mais sous forme d'une circulation du rotationnel sur une courbe fermée. Nous verrons comment discrétiser directement le théorème d'Ampère.

Pour un problème à symétrie axiale en coordonnées cylindriques, les vecteurs sont orthoradiaux et l'équation s'écrit :

r(1r(rAθ)r)+2Aθz2=sθ(5)

Cette équation est tout à fait différente de l'équation de Poisson pour un champ scalaire, dont la discrétisation est expliquée dans Discrétisation de l'équation de Poisson en géométrie axiale :

1rr(rur)+2uz2=s(6)

Pour l'équation vectorielle, les solutions à symétrie axiale invariantes par translation (indépendantes de z) sont :

Aθ(r)=12Br+Cr(7)

Pour l'équation scalaire, elles sont :

u(r)=-Br2+C(8)

L'objectif de ce document est la discrétisation de l'équation de Poisson vectorielle en géométrie axiale par la méthode des volumes finis, qui consiste à discrétiser la relation (3).

2. Géométrie plane

Dans un problème bidimensionnel en coordonnées cartésiennes, la forme intégrale (3) s'écrit :

C(-Azyex+Azxey)tdl=Σszdσ(9)

La figure suivante montre un point A du maillage, non situé sur un bord :

figureA.svgFigure pleine page

Pour discrétiser la forme intégrale, on considère le contour rectangulaire (a,b,c,d). La méthode des volumes finis consiste ici à utiliser non pas un volume de contrôle mais un contour de contrôle. Sur le côté (ab) on a :

(ab)(-Azyex+Azxey)tdl=(ab)Azydl(10)

Cette intégrale est identique au flux sortant du gradient de Az sur une surface (ab) (d'extension 1 selon z), qui est utilisé pour discrétiser le laplacien dans la méthode des volumes finis (Équation de Poisson à deux dimensions).

Sur le côté opposé (cd), on a :

(cd)(-Azyex+Azxey)tdl=(cd)-Azydl(11)

qui est aussi le flux sortant du gradient sur la surface (cd) d'extension 1 selon z.

On vérifie la même analogie pour les deux autres côtés. L'intégrale de surface de la source peut aussi se voir comme l'intégrale de volume sur le volume (a,b,c,d) d'extension 1 selon z.

Finalement, la forme intégrale sur ce contour est identique à la forme intégrale sur le volume de contrôle utilisé pour discrétiser l'équation de Poisson bidimensionnelle en coordonnées cartésiennes (Équation de Poisson à deux dimensions). Le schéma de discrétisation obtenu est donc identique :

Di,j=Δy2(12)Gi,j=Δy2(13)Hi,j=Δx2(14)Bi,j=Δx2(15)Ci,j=-2(Δx2+Δy2)(16)Fi,j=Si,jΔx2Δy2(17)

On retrouve cette identité en remarquant que l'équation de Poisson (4) est identique à l'équation de Poisson pour un champ scalaire. Nous allons voir qu'il n'est plus de même en coordonnées cylindriques.

3. Géométrie axiale

3.a. Discrétisation de l'équation

Pour un problème à symétrie axiale en coordonnées cylindriques, la forme intégrale s'écrit :

C(-1r(rAθ)rez+Aθzer)tdl=Σ sθdσ(18)

Le maillage est dans le plan méridien (z,r). La figure suivante représente le maillage avec z en abscisse :

figureA.svgFigure pleine page

Le contour de contrôle est (a,b,c,d). Sur le côté (ab), on a :

Iab=(ab)(-1r(rAθ)rez+Aθzer)tdl=(ab)1r(rAθ)rdl(19)

qui n'est évidemment pas le flux du gradient (dans le plan méridien le gradient est identique au gradient en coordonnées cartésiennes).

La manière la plus simple de discrétiser cette intégrale est :

Iab1rj+12rj+1Ai,j+1-rjAi,jΔrΔz(20)

Sur le côté opposé (cd) on a :

Icd=(cd)-1r(rAθ)rdl(21)

qui se discrétise de la manière suivante :

Icd1rj-12rj-1Ai,j-1-rjAi,jΔrΔz(22)

Pour le côté (bc) on a

Ibc=(bc)-AθzdlAi-1,j-Ai,jΔzΔr(23)

Le côté (da) est similaire et pour le terme de source on peut écrire :

ΣsθdσSi,jΔzΔr(24)

On obtient finalement le schéma suivant :

Di,j=Δr2(25)Gi,j=Δr2(26)Hi,j=j+1j+12Δz2(27)Bi,j=j-1j-12Δz2(28)Ci,j=-(jj+12Δz2+jj-12Δz2+2Δr2)(29)Fi,j=Si,jΔz2Δr2(30)

Lorsque j est grand (loin de l'axe), ce schéma est très proche du schéma donné plus haut pour le problème cartésien. Au voisinage de l'axe Oz, les deux schémas sont très différents.

Un point situé sur l'axe j=0 nécessite un traitement particulier puisque l'indice j-1 n'est pas défini pour ce point. On remarque que le vecteur source présente une antisymétrie par rapport à l'axe Oz :

sθ(z,-r)=-sθ(z,r)(31)

Il en est de même de Aθ(r,z). On en déduit que le vecteur A est nul sur l'axe :

Aθ(z,0)=0(32)

Un point de l'axe obéit donc à une condition limite de Dirichlet. Cela est tout à fait différent du cas de l'équation de Poisson scalaire, pour laquelle l'axe Oz est un axe de symétrie.

On remarque néanmoins que les points de l'axe n'interviennent pas dans le calcul des autres points puisque Bi,1=0.

3.b. Condition limite de Neumann

Considérons par exemple un point C situé sur le bord j=Nr-1 correspondant au rayon maximal.

figureB.svgFigure pleine page

Le contour de contrôle est (k,l,m,n). La condition limite de Neumann consiste à imposer sur le côté (kl) la composante sur z du rotationnel, c'est-à-dire la dérivée suivante :

wr=1r(rAθ)r(33)

La discrétisation de la circulation sur le côté (kl) s'écrit donc :

Ikl=wr,iΔz(34)

Pour attribuer le bon signe, il suffit de se rappeler que la circulation est analogue à un flux sortant de gradient. Sur les côtés (lm) et (nk), la discrétisation est identique à celle vue plus haut, avec une longueur de segment réduite de moitié. Par exemple :

Ilm=(lm)-AθzdlAi-1,j-Ai,jΔzΔr2(35)

Pour le terme de source, il faut tenir compte de la réduction de l'aire de la surface :

ΣsθdσSi,j12ΔzΔr(36)

Le schéma obtenu pour ce point est finalement :

Di,j=Gi,j=12Δr2(37)Hi,j=0(38)Bi,j=j-1j-12Δz2(39)Ci,j=-(jj-12Δz2+Δr2)(40)Fi,j=Si,j12Δz2Δr2-wr,jΔrΔz2(41)

Une condition limite sur un côté où i est constant (z constant) fait intervenir la composante r du rotationnel, c'est-à-dire que l'on impose la dérivée suivante :

wz=Aθz(42)

Voyons comment traiter le point E situé sur un coin du domaine :

figureC.svgFigure pleine page

Le contour est (s,t,u,v). La dérivée wr est imposée sur le côté (st), la dérivée wz est imposée sur le côté (tu). Le schéma est :

Di,j=12Δr2(43)Gi,j=0(44)Hi,j=0(45)Bi,j=12j-1j-12Δz2(46)Ci,j=-(12jj-12Δz2+Δr2)(47)Fi,j=Si,j14Δz2Δr2-wr,i12ΔrΔz2+wz,j12ΔzΔr2(48)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.