Discrétisation de l'équation de Poisson vectorielle

1. Introduction

L'équation de Poisson vectorielle est :

\nabla^{2} \vec{A} = \vec{s} (1)

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

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

- \vec{r o t} (\vec{r o t} \vec{A}) = \vec{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 Σ :

\oint_{C} (\vec{r o t} \vec{A}) \cdot \vec{t} d l = - \iint_{Σ} s d σ (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 :

\frac{\partial^{2} A_{z}}{\partial x^{2}} + \frac{\partial^{2} A_{z}}{\partial y^{2}} = s_{z} (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 :

\frac{\partial}{\partial r} (\frac{1}{r} \frac{\partial (r A_{θ})}{\partial r}) + \frac{\partial^{2} A_{θ}}{\partial z^{2}} = 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 :

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial u}{\partial r}) + \frac{\partial^{2} u}{\partial z^{2}} = s (6)

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

A_{θ} (r) = \frac{1}{2} B r + \frac{C}{r} (7)

Pour l'équation scalaire, elles sont :

u (r) = - \frac{B}{r^{2}} + 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 :

\oint_{C} (- \frac{\partial A_{z}}{\partial y} {\vec{e}}_{x} + \frac{\partial A_{z}}{\partial x} {\vec{e}}_{y}) \cdot \vec{t} d l = \iint_{Σ} s_{z} d σ (9)

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

Figure 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 :

\int_{(a b)} (- \frac{\partial A_{z}}{\partial y} {\vec{e}}_{x} + \frac{\partial A_{z}}{\partial x} {\vec{e}}_{y}) \cdot \vec{t} d l = \int_{(a b)} \frac{\partial A_{z}}{\partial y} d l (10)

Cette intégrale est identique au flux sortant du gradient de A_z 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 :

\int_{(c d)} (- \frac{\partial A_{z}}{\partial y} {\vec{e}}_{x} + \frac{\partial A_{z}}{\partial x} {\vec{e}}_{y}) \cdot \vec{t} d l = \int_{(c d)} - \frac{\partial A_{z}}{\partial y} d l (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 :

\begin{matrix} D_{i, j} = Δ y^{2} & (12) \\ G_{i, j} = Δ y^{2} & (13) \\ H_{i, j} = Δ x^{2} & (14) \\ B_{i, j} = Δ x^{2} & (15) \\ C_{i, j} = - 2 (Δ x^{2} + Δ y^{2}) & (16) \\ F_{i, j} = S_{i, j} Δ x^{2} Δ y^{2} & (17) \end{matrix}

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 :

\oint_{C} (- \frac{1}{r} \frac{\partial (r A_{θ})}{\partial r} {\vec{e}}_{z} + \frac{\partial A_{θ}}{\partial z} {\vec{e}}_{r}) \cdot \vec{t} d l = \iint_{Σ} s_{θ} d σ (18)

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

Figure pleine page

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

I_{a b} = \int_{(a b)} (- \frac{1}{r} \frac{\partial (r A_{θ})}{\partial r} {\vec{e}}_{z} + \frac{\partial A_{θ}}{\partial z} {\vec{e}}_{r}) \cdot \vec{t} d l = \int_{(a b)} \frac{1}{r} \frac{\partial (r A_{θ})}{\partial r} d l (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 :

I_{a b} ≃ \frac{1}{r_{j + \frac{1}{2}}} \frac{r_{j + 1} A_{i, j + 1} - r_{j} A_{i, j}}{Δ r} Δ z (20)

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

I_{c d} = \int_{(c d)} - \frac{1}{r} \frac{\partial (r A_{θ})}{\partial r} d l (21)

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

I_{c d} ≃ \frac{1}{r_{j - \frac{1}{2}}} \frac{r_{j - 1} A_{i, j - 1} - r_{j} A_{i, j}}{Δ r} Δ z (22)

Pour le côté (bc) on a

I_{b c} = \int_{(b c)} - \frac{\partial A_{θ}}{\partial z} d l ≃ \frac{A_{i - 1, j} - A_{i, j}}{Δ z} Δ r (23)

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

\iint_{Σ} s_{θ} d σ ≃ S_{i, j} Δ z Δ r (24)

On obtient finalement le schéma suivant :

\begin{matrix} D_{i, j} = Δ r^{2} & (25) \\ G_{i, j} = Δ r^{2} & (26) \\ H_{i, j} = \frac{j + 1}{j + \frac{1}{2}} Δ z^{2} & (27) \\ B_{i, j} = \frac{j - 1}{j - \frac{1}{2}} Δ z^{2} & (28) \\ C_{i, j} = - (\frac{j}{j + \frac{1}{2}} Δ z^{2} + \frac{j}{j - \frac{1}{2}} Δ z^{2} + 2 Δ r^{2}) & (29) \\ F_{i, j} = S_{i, j} Δ z^{2} Δ r^{2} & (30) \end{matrix}

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 B_i,1=0.

3.b. Condition limite de Neumann

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

Figure 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 :

w_{r} = \frac{1}{r} \frac{\partial (r A_{θ})}{\partial r} (33)

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

I_{k l} = w_{r, 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 :

I_{l m} = \int_{(l m)} - \frac{\partial A_{θ}}{\partial z} d l ≃ \frac{A_{i - 1, j} - A_{i, j}}{Δ z} \frac{Δ r}{2} (35)

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

\iint_{Σ} s_{θ} d σ ≃ S_{i, j} \frac{1}{2} Δ z Δ r (36)

Le schéma obtenu pour ce point est finalement :

\begin{matrix} D_{i, j} = G_{i, j} = \frac{1}{2} Δ r^{2} & (37) \\ H_{i, j} = 0 & (38) \\ B_{i, j} = \frac{j - 1}{j - \frac{1}{2}} Δ z^{2} & (39) \\ C_{i, j} = - (\frac{j}{j - \frac{1}{2}} Δ z^{2} + Δ r^{2}) & (40) \\ F_{i, j} = S_{i, j} \frac{1}{2} Δ z^{2} Δ r^{2} - w_{r, j} Δ r Δ z^{2} & (41) \end{matrix}

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 :

w_{z} = \frac{\partial A_{θ}}{\partial z} (42)

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

Figure pleine page

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

\begin{matrix} D_{i, j} = \frac{1}{2} Δ r^{2} & (43) \\ G_{i, j} = 0 & (44) \\ H_{i, j} = 0 & (45) \\ B_{i, j} = \frac{1}{2} \frac{j - 1}{j - \frac{1}{2}} Δ z^{2} & (46) \\ C_{i, j} = - (\frac{1}{2} \frac{j}{j - \frac{1}{2}} Δ z^{2} + Δ r^{2}) & (47) \\ F_{i, j} = S_{i, j} \frac{1}{4} Δ z^{2} Δ r^{2} - w_{r, i} \frac{1}{2} Δ r Δ z^{2} + w_{z, j} \frac{1}{2} Δ z Δ r^{2} & (48) \end{matrix}