Équation de Poisson : méthode des éléments finis

1. Introduction

Ce document montre comment résoudre numériquement l'équation de Poisson (à deux dimensions) par la méthode des éléments finis.

Après avoir énoncé la formulation variationnelle du problème, nous verrons comment utiliser le logiciel FreeFEM à travers quelques exemples (d'électrostatique et de magnétostatique) puis comment traiter les résultats dans un script Python.

2. Formulation variationnelle

La résolution numérique d'une équation aux dérivées partielles par la méthode des éléments finis (M.E.F.) repose sur une formulation variationnelle du problème. L'équation considérée ici est l'équation de Poisson à deux dimensions, vérifiée par une fonction inconnue u(x,y) de classe C² sur un domaine Ω (ouvert de $R^{2}$ ) :

\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} = - f (x, y) (1)

où f est une fonction donnée sur le domaine Ω.

La condition limite imposée sur le bord du domaine ( $\partial Ω$ ) est de type Dirichlet :

u (x, y) = 0 s u r \partial Ω (2)

La formulation variationnelle du problème de Poisson ([1],[2]) est obtenue en multipliant l'équation $(1)$ par une fonction v(x,y) (fonction test) puis en intégrant sur un domaine Ω :

\int_{Ω} (Δ u) v d τ = - \int_{Ω} f v d τ (3)

où Δ désigne l'opérateur laplacien et dτ l'élément de volume (ici le volume a deux dimensions)

Pour une fonction w(x,y) de classe C¹ sur Ω, nous utilisons le théorème suivant (formule de Green) :

\int_{Ω} \frac{\partial w}{\partial x} d τ = \int_{\partial Ω} w n_{x} d s (4)

où ds désigne l'élément de surface de $\partial Ω$ (ici la surface a une dimension) et n_x la composante sur l'axe x du vecteur normal unitaire sortant de Ω. La même formule est évidemment valable pour la coordonnée y.

L'application de la formule de Green à w=uv conduit à la relation suivante (formule d'intégration par partie) :

\int_{Ω} v \frac{\partial u}{\partial x} d τ = - \int_{Ω} u \frac{\partial v}{\partial x} d τ + \int_{\partial Ω} u v n_{x} d s (5)

Appliquons cette relation à $\frac{\partial u}{\partial x}$ (à la place de u) :

\int_{Ω} v \frac{\partial^{2} u}{\partial x^{2}} d τ = - \int_{Ω} \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} d τ + \int_{\partial Ω} \frac{\partial u}{\partial x} v n_{x} d s (6)

La même relation étant valable pour la coordonnée y, la somme des deux relations conduit à :

\int_{Ω} v Δ u d τ = - \int_{Ω} (\vec{g r a d} u) \cdot (\vec{g r a d} v) d τ + \int_{\partial Ω} v (\vec{g r a d} u) \cdot \vec{n} d s (7)

Reprenons la relation $(3)$ et supposons que la fonction test v s'annule sur $\partial Ω$ , comme la fonction u. L'application de la formule précédente conduit à :

\int_{Ω} (\vec{g r a d} u) \cdot (\vec{g r a d} v) d τ = \int_{Ω} f v d τ (8)

La fonction u recherchée est de classe C¹ (et non pas de classe C²). La formulation variationnelle consiste à énoncer que l'égalité $(8)$ doit être vérifiée pour toute fonction test v de classe C¹.

3. Méthode des éléments finis

L'exposé détaillé de la méthode des éléments finis ([2]) sort du cadre de ce document mais nous donnons brièvement son principe.

La méthode des éléments finis (M.E.F.) repose sur la triangulation du domaine Ω, c'est-à-dire la subdivision en triangles. Le maillage est donc constitué de triangles. Chaqu'un des trois sommets d'un triangle est aussi le sommet d'un triangle voisin, ou éventuellement appartient au bord du domaine. L'avantage d'un maillage par triangulation, par comparaison aux grilles utilisées pour la méthode des différences finies, est la facilité avec laquelle le maillage s'ajuste aux formes complexes, par exemple des arcs de cercle. Par ailleurs, un maillage par triangles peut facilement être affiné en certains endroits car il est aisé de subdiviser un triangle en plusieurs triangles.

Soit V_h l'ensemble des fonctions continues sur Ω et affines sur chaque triangle, c'est-à-dire de la forme ax+by+c sur chaque triangle. Cet ensemble admet une base de fonctions aisée à définir ([2]). Si i désigne un nœud du maillage ( $i = 1, 2, \dots N$ ), on note $ϕ_{i}$ la fonction appartenant à V_h, vallant 1 sur ce nœud et 0 sur tous les autres. L'ensemble des fonctions $ϕ_{i}$ constitue une base de V_h. En effet, pour toute fonction u_h de V_h on a :

u_{h} (x, y) = \sum_{i = 1}^{N} u_{i} ϕ_{i} (x, y) (9)

où u_i est la valeur de u_h au nœud i.

La solution approchée de l'équation de Poisson avec conditions limites est recherchée sous la forme d'une fonction u_h continue et affine par morceaux (c.a.d. affine sur chaque triangle).

La M.E.F. consiste à rechercher les composantes u_i de u_h sur la base $(ϕ_{i})$ et à énoncer que la relation $(8)$ doit s'appliquer :

\int_{Ω} (\vec{g r a d} u_{h}) \cdot (\vec{g r a d} v) d τ = \int_{Ω} f v d τ (10)

pour toute fonction v de V_h. En conséquence, cette relation doit aussi être vraie en remplaçant v par $ϕ_{i}$ quel que soit le nœud i.

On a donc :

\sum_{j = 1}^{N} u_{j} \int_{Ω} (\vec{g r a d} ϕ_{j}) \cdot (\vec{g r a d} ϕ_{i}) d τ = \int_{Ω} f ϕ_{i} d τ (11)

pour $i = 1, 2, \dots N$ . On obtient ainsi un système de N équations linéaires pour les N inconnues u_i :

\sum_{j = 1}^{N} A_{i j} u_{i} = b_{i} (12)

avec les coefficients suivants :

\begin{matrix} A_{i j} = \int_{Ω} (\vec{g r a d} ϕ_{j}) \cdot (\vec{g r a d} ϕ_{i}) d τ & (13) \\ b_{i} = \int_{Ω} f ϕ_{i} d τ & (14) \end{matrix}

Les calculs des coefficients A_ij et b_i font intervenir une intégrale qui ne s'étend que sur quelques triangles voisins car le support de $ϕ_{i}$ se limite aux triangles dont le nœud d'indice i est un sommet.

Le système est résolu numériquement soit par une méthode directe soit par une méthode itérative.

Le point important à retenir est que la fonction u_h obtenue est une fonction continue affine par morceaux (c.a.d. affine sur chaque triangle) qui constitue une approximation de la fonction u(x,y) vérifiant l'équation de Poisson avec la condition limite. Le calcul fournit les valeurs aux sommets des triangles (les nœuds du maillage), qui sont les composantes u_i de u sur la base de V_h. D'après la relation $(9)$ , la connaissance de ces valeurs suffit en effet à déterminer entièrement la fonction continue et affine sur chaque triangle.

4. Utilisation de FreeFem

Le logiciel FreeFEM permet de résoudre des problèmes par la M.E.F. Les calculs à effectuer sont définis dans un script écrit dans un langage spécifique facile à apprendre (très proche de C/C++). Les représentations graphiques des résultats sont affichées sur l'écran et peuvent être exportées au format EPS (encapsulated Postscript), qu'il est possible de convertir en format PNG et SVG avec un convertisseur en ligne ou bien un logiciel Adobe (Photoshop ou Illustrator). Cette solution d'exportation étant peu pratique, nous allons plutôt exporter des résultats dans des fichiers textes afin de les récupérer en Python.

La documentation de FreeFEM contient une description détaillée du langage, une introduction à la formulation variationnelle et de nombreux exemples.

Les auteurs de FreeFem ont écrit des fonctions pour Matlab et Octave (freefem_matlab_octave_plot) qui permettent de tracer les résultats de simulations dans ces deux logiciels. Nous utilisons le fichier de macros ffmatlib.idp qui peut être téléchargé dans freefem_matlab_octave_plot. Ce fichier contient quelques fonctions (en langage FreeFem) permettant d'exporter les données sous forme de fichier texte. Nous les reproduisons ci-dessous :

//Save finite element connectivity (j-th degree of freedom for i-th element)
macro ffSaveVh(Th, Vh, filename){
  {
  ofstream file(filename);
  file.precision(10);
  for (int i=0; i<Th.nt; i++){
    for (int j=0; j<Vh.ndofK; j++){
      file << Vh(i,j) << "\n";
    }
  }
  }
} //EOM

//Save FE Space data
macro ffSaveData(u, filename){
  {
  ofstream file(filename);
  file.precision(10);
  for (int j=0; j<u[].n; j++){
    file << u[][j] << endl;
  }
  }
} //EOM

macro ffSaveData2(u1, u2, filename){
  {
  ofstream file(filename);
  int datalen=u1[].n;
  if (u2[].n!=datalen){
    cout << "error: arguments must have same size" << endl;
    exit(1);
  }
  file.precision(10);
  for (int j=0; j<datalen; j++){
    file << u1[][j] << " " << u2[][j] << endl;
  }
  }
} //EOM

macro ffSaveData3(u1, u2, u3, filename){
  {
  ofstream file(filename);
  int datalen=u1[].n;
  if ((u2[].n!=datalen) | (u3[].n!=datalen)){
    cout << "error: arguments must have same size" << endl;
    exit(1);
  }
  file.precision(10);
  for (int j=0; j<datalen; j++){
    file << u1[][j] << " " << u2[][j] << " " << u3[][j] << endl;
  }
  }
} //EOM

macro ffSaveData4(u1, u2, u3, u4, filename){
  {
  ofstream file(filename);
  int datalen=u1[].n;
  if ((u2[].n!=datalen) | (u3[].n!=datalen) | (u4[].n!=datalen)){
    cout << "error: arguments must have same size" << endl;
    exit(1);
  }
  file.precision(10);
  for (int j=0; j<datalen; j++){
    file << u1[][j] << " " << u2[][j] << " " << u3[][j] << " " << u4[][j] << endl;
  }
  }
} //EOM

5. Exemples

5.a. Champ magnétique créé par un courant rectiligne

Considérons un exemple simple de magnétostatique : le champ magnétique créé par un courant électrique uniforme dans un fil rectiligne infini de rayon a. Le fil est supposé parallèle à l'axe (Oz).

Sachant que $d i v \vec{B} =$ il existe un potentiel vecteur $\vec{A}$ tel que $\vec{B} = \vec{r o t} \vec{A}$ . En magnétostatique, on impose la condition $d i v \vec{A} = 0$ . Le potentiel vecteur vérifie alors l'équation de Poisson vectorielle :

Δ \vec{A} = - μ_{0} \vec{J} (15)

Pour cet exemple, le potentiel vecteur s'écrit $\vec{A} = A_{z} (x, y) {\vec{u}}_{z}^{}$ et l'équation de Poisson :

Δ A_{z} = - μ_{0} J_{z} (16)

La densité de courant est J_z=J (uniforme) dans le fil, nulle en dehors du fil. On prendra comme domaine Ω un disque centré sur l'axe (Oz). Sur le bord du disque, on pose A_z=0. Ce problème modélise un fil dans l'espace infini dans la mesure où le rayon du domaine est assez grand par rapport au rayon du fil.

Le champ magnétique se calcule à partir du potentiel vecteur par :

\begin{matrix} B_{x} = \frac{\partial A_{z}}{\partial y} & (17) \\ B_{y} = - \frac{\partial A_{z}}{\partial y} & (18) \end{matrix}

Il faut remarquer que les lignes d'isovaleurs de A_z sont des lignes de champ magnétique.

En prendra le terme de source (la fonction f dans l'énoncé ci-dessus) égal à 1 dans le fil, 0 en dehors. Il faudra donc multiplier le potentiel A_z obtenu par $μ_{0} J$ .

Notons u=A_z et explicitons l'équation $(8)$ sous forme d'intégrales doubles :

\iint (\frac{\partial u}{\partial x} \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \frac{\partial v}{\partial y}) d x d y = \iint f v d x d y (19)

Voici un script FreeFem (extension EDP) définissant les calculs à réaliser :

filRectiligne.edp

include "ffmatlib.idp"

real Rdomaine = 1.0; // rayon du domaine
real rayon = 0.02; // rayon du fil

// fonction source (mu0*J=1)
func real source(real &x, real &y){
	if (sqrt(x*x+y*y)<rayon) {return 1.0;}
	else {return 0;}
}

border C(t=0,2*pi){x=Rdomaine*cos(t);y=Rdomaine*sin(t);} // bord du domaine

int N=50;
mesh Th = buildmesh(C(N)); // triangulation
fespace Vh(Th,P1);
plot(Th, wait=true);
Vh u,v;

problem Poisson(u,v,eps=1.0e-6) = int2d(Th) (dx(u)*dx(v)+dy(u)*dy(v)) -int2d(Th)(source(x,y)*v)
	+on(C,u=0);
	
Poisson;
plot(u,wait=1);
int i;
real error = 0.01;
for (i = 0; i < 4; i++) { // boucle d'adaptation du maillage
	Poisson;
	cout << u[].min << " " << u[].max << endl;
	plot(u,wait=1);
	Th = adaptmesh(Th, u, err=error);
	plot(Th, wait=true); 
	u = u;
	error = error/2;
}

Vh Bx,By;
By = -dx(u);
Bx = dy(u);

plot(Th,wait=true);
plot(u,wait=true);

savemesh(Th,"filRectiligne.msh");
ffSaveVh(Th,Vh,"filRectiligne_vh.txt");
ffSaveData3(u,Bx,By,"filRectiligne_data.txt");

Le fichier ffmatlib.idp doit être dans le même répertoire que le script. Pour exécuter FreeFEM sur ce script, ouvrir une console est exécuter :

freefem++ filRectiligne.edp

Chaque affichage graphique avec l'option wait=true attend un appuie sur Enter.

La fonction inconnue est le potentiel A_z mais elle est notée u dans le script pour se conformer aux notations générales définies plus haut. La fonction test est notée v.

Le domaine est un disque délimité par le cercle défini par la commande border. La commande buildmesh génère un maillage avec un nombre N de sommets sur le bord. La commande plot(Th, wait=true) fournit une représentation graphique de ce maillage.

La fonction fespace permet de définir un espace d'éléments finis. Dans le cas présent, on définit V_h, l'espace des fonctions affines par triangle (P1) sur le maillage défini.

La fonction problem permet de définir le problème sous sa forme variationnelle (équation $(8)$ ). Les intégrales sur le domaine Ω sont définies par l'opérateur int2d(Th). Il faut remarquer que les deux intégrales doivent être définies séparément. Dans le premier terme, dx(u) désigne la dérivée partielle de u par rapport à x. La condition limite sur le bord du domaine est définie par +on(C,u=0).

La ligne Poisson; effectue une première résolution. La commande plot(u,wait=true) effectue le tracé des lignes d'isovaleurs de u (c'est-à-dire du potentiel vecteur). Le premier maillage défini est trop grossier pour traiter correctement la zone proche du fil (dont le rayon est très faible par rapport au rayon du domaine). Il faut donc réaliser une adaptation du maillage. La fonction adaptmesh effectue le calcul d'un nouveau maillage (par subdivision du précédent) en fonction de la fonction u déjà obtenue et d'une valeur d'erreur. Plusieurs adaptations avec une erreur de plus en plus petite sont en général nécessaires (ici 4). Pour finir, les composantes du champ magnétique sont calculées, sous forme de fonctions de V_h.

Le fichier filRectiligne.msh contient la définition du maillage : nombre de sommets, nombre de triangles, coordonnées des sommets puis indices des sommets des triangles. Dans cet exemple, il y a 5253 sommets et 10402 triangles.

Le fichier filRectiligne_data.txt contient les valeurs de u,Bx,By aux nœuds du maillage, qui constituent les composantes de ces fonctions sur la base $(ϕ_{i})$ .

Afin d'importer la définition du maillage et les valeurs de u,Bx,By aux nœuds du maillage, nous utilisons la class Python ci-dessous :

FreeFemData.py

import numpy as np

class FreeFemData:
    def __init__(self,name):
        self.name = name
        with open(name+".msh") as f:
            lines = f.readlines()
        line0 = lines[0].rstrip()
        data = line0.split(" ")
        self.nb_sommets  = int(data[0])
        self.nb_triangles = int(data[1])
        self.x_array = np.zeros(self.nb_sommets)
        self.y_array = np.zeros(self.nb_sommets)
        self.z_array = np.zeros(self.nb_sommets)
        j = 1
        for i in range(0,self.nb_sommets):
            line = lines[j].rstrip()
            data = line.split(" ")
            j += 1
            self.x_array[i] = float(data[0])
            self.y_array[i]= float(data[1])
            self.z_array[i] = float(data[2])
        self.triangles = np.zeros((self.nb_triangles,3),dtype=np.int64)
        for i in range(self.nb_triangles):
            line = lines[j].rstrip()
            data = line.split(" ")
            j += 1
            self.triangles[i,0] = int(data[0])-1
            self.triangles[i,1] = int(data[1])-1
            self.triangles[i,2] = int(data[2])-1

    def readData(self,n):
        self.data = np.zeros((n,self.nb_sommets))
        with open('%s_data.txt'%self.name) as f:
            lines = f.readlines()
        for i in range(0,self.nb_sommets):
            line = lines[i].rstrip()
            data = line.split(" ")
            for k in range(n):
                self.data[k,i] = float(data[k])
        return self.data

Pour créer une instance de FreeFemData, on fournit la racine du nom des fichiers (ici filRectiligne). Les données obtenues sont :

nb_sommets : nombre de sommets.
nb_triangles : nombre de triangles.
x_array,y_array,z_array : tableaux contenant les coordonnées des sommets.
triangles : tableau contenant les triangles : sur chaque ligne un triangle défini par 3 indices.

La fonction readData permet de lire un fichier de données contenant n définitions de fonctions de V_h. Le résultat est renvoyé sous la forme d'un tableau de n lignes. La ligne i contient les valeurs de la fonction d'indice i aux nœuds du maillage.

Voyons comment traiter les données produites par le script filRectiligne.edp. Nous utilisons pour cela le module matplotlib.tri, qui comporte des fonctions permettant de gérer les maillages triangulaires.

On commence par lire le maillage :

import numpy as np
from matplotlib.pyplot import *
import matplotlib.tri as tri
from FreeFemData import FreeFemData

freeFemData = FreeFemData('filRectiligne')

La fonction matplotlib.tri.Triangulation permet d'obtenir un maillage triangulaire à partir de celui fourni par FreeFem :

triang = tri.Triangulation(freeFemData.x_array,freeFemData.y_array,freeFemData.triangles)

Voici comment obtenir une représentation graphique de ce maillage :

fig0,ax0 = subplots(figsize=(10,10))
ax0.set_aspect('equal')
ax0.triplot(triang,color='0.8')
ax0.set_xlabel('x',fontsize=16)
ax0.set_ylabel('y',fontsize=16)

plotA.pdf

L'adaptation du maillage a permis d'obtenir un maillage très fin au centre (au voisinage du fil). Voici un détail de la zone centrale, avec un tracé du cercle délimitant le fil :

rayon = 0.02
theta = np.linspace(0,2*np.pi,1000)
x = rayon*np.cos(theta)
y = rayon*np.sin(theta)
fig0,ax0 = subplots(figsize=(10,10))
ax0.set_aspect('equal')
ax0.triplot(triang,color='0.8')
ax0.plot(x,y,'r-')
ax0.set_xlim(-rayon*2,rayon*2)
ax0.set_ylim(-rayon*2,rayon*2)
ax0.set_xlabel('x',fontsize=16)
ax0.set_ylabel('y',fontsize=16)

plotB.pdf

La fonction tricontour permet de tracer des lignes d'isovaleurs de A_z :

data = freeFemData.readData(3) 
Az_array = data[0,:]
Bx_array = data[1,:] 
By_array = data[2,:]

fig1,ax1 = subplots(figsize=(10,10))
ax1.set_aspect('equal')
ct = ax1.tricontour(triang,Az_array,levels=40)
fig1.colorbar(ct)
ax1.set_xlabel('x',fontsize=16)
ax1.set_ylabel('y',fontsize=16)

plotC.pdf

Cherchons à présent à tracer le champ magnétique sur l'axe (Ox), qui en principe ne doit comporter que la composante B_y. La méthode des éléments finis fournit des valeurs sur les sommets des triangles. Il faut donc procéder à une interpolation pour obtenir les valeurs de B_z en des points répartis sur une droite. Dans le cas présent, l'espace V_h est celui des fonctions continues affines par triangle. Il faut donc procéder à une interpolation linéaire. La fonction matplotlib.tri.LinearTriInterpolator génère une fonction d'interpolation à partir d'un maillage triangulaire et d'un ensemble de valeurs sur les nœuds de ce maillage.

Bx_inter = tri.LinearTriInterpolator(triang, Bx_array)
By_inter = tri.LinearTriInterpolator(triang, By_array)
x = np.linspace(-1,1,1000)
y = np.ones(len(x))*0
By = By_inter(x,y)
figure(figsize=(8,6))
plot(x,By)
xlabel('x',fontsize=16)
ylabel('By',fontsize=16)
grid()

plotD.pdf

Comme prévu, le champ B_y présente un profil affine à l'intérieur du fil et décroît à l'extérieur proportionnellement à l'inverse de la distance à l'axe. La calcul a été fait avec $μ_{0} J_{z} = 1$ . L'intensité du champ magnétique à la surface du fil est donc théoriquement égal à a/2. Voici la courbe montrant les valeurs de A_z sur la droite (Ox).

Az_inter = tri.LinearTriInterpolator(triang, Az_array)
figure(figsize=(8,6))
plot(x,Az_inter(x,y))
xlabel('x',fontsize=16)
ylabel('Az',fontsize=16)
grid()

plotE.pdf

5.b. Champ électrique créé par un conducteur chargé

Ce deuxième exemple n'est pas un problème de Poisson à proprement parler mais un problème de Laplace.

Une plaque métallique a une épaisseur e dans la direction (Ox), une largeur L dans la direction (Oy) et une longueur infinie dans la direction (Oz). Elle est chargée et son potentiel est donc V₁ par rapport au potentiel à l'infini. Le potentiel est V(x,y). On a en tout point de l'espace en dehors du conducteur :

Δ V = 0 (20)

Le bord $\partial Ω$ est constitué de deux parties. La première est la surface du conducteur, sur laquelle on a la condition limite V=V₁. La seconde est la surface fermée qui englobe le conducteur, sur laquelle on a la condition limite V=0.

Ce problème est très voisin du problème de Poisson. La différence réside dans la condition limite V=V₁ sur la surface du conducteur.

Si u désigne le potentiel V, la forme variationnelle donnée par l'équation $(8)$ est toujours valable car la fonction test v est choisie s'annulant sur toute la frontière $\partial Ω$ , y compris sur le conducteur.

Nous choisissons un domaine circulaire (comme dans l'exemple précédent). Ce choix permet d'éviter les effets de coins que l'on rencontre avec des domaines rectangulaires. En effet, le potentiel à grande distance de la plaque doit présenter une quasi invariance par rotation autour de l'axe (Oz).

Voici un script FreeFem définissant les calculs à réaliser :

plaqueConductrice.edp

include "ffmatlib.idp"
real Rdomaine = 200;
real e = 1; //épaisseur de la plaque
real L = 50; // largeur de la plaque
real e2 = e/2;
real L2 = L/2;
int C1 = 99;

border C(t=0,2*pi){x=Rdomaine*cos(t);y=Rdomaine*sin(t);}
border C11(t=0, 1){ x=e2;  y=L2-t*L;  label=C1;}
border C12(t=0, 1){ x=e2-e*t;  y=-L2;  label=C1;}
border C13(t=0, 1){ x=-e2;  y=-L2+t*L; label=C1;}
border C14(t=0, 1){ x=-e2+e*t;  y=L2; label=C1;} 
mesh Th = buildmesh(C(100)+C11(5)+C12(20)+C13(5)+C14(20));
fespace Vh(Th,P1);
Vh u,v;
plot(Th, wait=1);

problem Laplace(u,v,eps=1.0e-6) = int2d(Th) (dx(u)*dx(v)+dy(u)*dy(v))
	+on(C,u=0)
	+on(C1,u=1);
Laplace;
int i;
real error = 0.01;
for (i = 0; i < 4; i++) {
	Laplace;
	cout << u[].min << " " << u[].max << endl;
	plot(u,wait=1);
	Th = adaptmesh(Th, u, err=error);
	plot(Th, wait=1);
	u = u;
	error = error/2;
}
Vh Ex,Ey;
Ex = -dx(u);
Ey = -dy(u);

plot(Th,wait=true); 
plot(u,wait=true);
plot(u,[Ex,Ey]);
savemesh(Th,"plaqueConductriceChargee.msh");
ffSaveVh(Th,Vh,"plaqueConductriceChargee_vh.txt");
ffSaveData3(u,Ex,Ey,"plaqueConductriceChargee_data.txt");

Les bords du conducteur sont définis segment par segment. La fonction buildmesh construit le maillage avec le cercle C et les bords du conducteur comme frontières. La formulation variationnelle donnée est la même que précédemment (sans le terme de source) mais on précise la condition limite sur le conducteur par +on(C1,u=1).

Voici une vue globale du maillage obtenu après adaptation :

freeFemData = FreeFemData('plaqueConductriceChargee')
triang = tri.Triangulation(freeFemData.x_array,freeFemData.y_array,freeFemData.triangles)
Rdomaine = 100
e = 1 # épaisseur de la plaque
L = 50 # largeur de la plaque
e2=e/2
L2=L/2					 
fig0,ax0 = subplots(figsize=(10,10))
ax0.set_aspect('equal')
ax0.triplot(triang,color='0.8')
ax0.plot([-e2,e2,e2,-e2,-e2],[-L2,-L2,L2,L2,-L2],'r-')
ax0.set_xlabel('x',fontsize=16)
ax0.set_ylabel('y',fontsize=16)

plotF.pdf

et un détail au voisinage d'un bord de la plaque :

fig0,ax0 = subplots(figsize=(10,10))
ax0.set_aspect('equal')
ax0.triplot(triang,color='0.8')
ax0.plot([-e2,e2,e2,-e2,-e2],[-L2,-L2,L2,L2,-L2],'r-')
ax0.set_xlabel('x',fontsize=16)
ax0.set_ylabel('y',fontsize=16)
ax0.set_xlim(-10,10)
ax0.set_ylim(15,35)

plotG.pdf

L'adaptation du maillage a conduit à une forte densité de triangles au voisinage du bord.

Voici un tracé de lignes équipotentielles :

data = freeFemData.readData(3)
V_array = data[0,:]
Ex_array = data[1,:]
Ey_array = data[2,:]

fig1,ax1 = subplots(figsize=(10,10))
ax1.set_aspect('equal')
ct = ax1.tricontour(triang,V_array,levels=40)
fig1.colorbar(ct)
ax1.plot([-e2,e2,e2,-e2,-e2],[-L2,-L2,L2,L2,-L2],'r-')
ax1.set_xlabel('x',fontsize=16)
ax1.set_ylabel('y',fontsize=16)

plotH.pdf

Voici les lignes équipotentielles au voisinage d'un bord :

fig1,ax1 = subplots(figsize=(10,10))
ax1.set_aspect('equal')
ct = ax1.tricontour(triang,V_array,levels=40)
fig1.colorbar(ct)
ax1.plot([-e2,e2,e2,-e2,-e2],[-L2,-L2,L2,L2,-L2],'r-')
ax1.set_xlabel('x',fontsize=16)
ax1.set_ylabel('y',fontsize=16)	
ax1.set_xlim(-10,10)
ax1.set_ylim(15,35)

plotI.pdf

Champ électrique sur l'axe (Ox) :

Ex_inter = tri.LinearTriInterpolator(triang, Ex_array)
Ey_inter = tri.LinearTriInterpolator(triang, Ey_array)
x = np.linspace(e2,Rdomaine,1000)
y = np.ones(len(x))*0
Ex = Ex_inter(x,y)
figure(figsize=(8,6))
plot(x,Ex)
xlabel('x',fontsize=16)
ylabel('Ex',fontsize=16)
grid()
ylim(0,0.1)

plotK.pdf

et sur l'axe (Oy) :

y = np.linspace(L2,Rdomaine,1000)
x = np.ones(len(y))*0
Ey = Ey_inter(x,y)
figure(figsize=(8,6))
plot(y,Ey)
xlabel('y',fontsize=16)
ylabel('Ey',fontsize=16)
grid()
ylim(0,0.1)

plotL.pdf

5.c. Conducteurs chargés en symétrie axiale

Considérons un problème d'électrostatique présentant une symétrie axiale, c'est-à-dire une invariance par rotation autour de l'axe (Oz). Les coordonnées cylindriques (r,z) permettent de définir un problème bidimensionnel puisque V(r,z).

L'équation de Poisson prend la forme suivante (le potentiel est noté u) :

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial u}{\partial r}) + \frac{\partial^{2} u}{\partial z^{2}} = - f (z, r) (21)

Les conditions limites sur $\partial Ω$ sont u=0 sur la surface englobante et u=V_c sur la surface d'un conducteur dont le potentiel est V_c (par rapport au potentiel à l'infini).

La fonction test v est supposée s'annuler sur $\partial Ω$ . La formulation variationnelle prend donc la même forme que précédemment :

\int_{Ω} (\vec{g r a d} u) \cdot (\vec{g r a d} v) d τ = \int_{Ω} f v d τ (22)

Cependant, on doit ici raisonner sur un domaine Ω de R³ présentant une symétrie axiale. L'intégrale sur ce volume s'écrit sous la forme d'intégrales doubles sur (r,z) :

\iint (\frac{\partial u}{\partial r} \frac{\partial v}{\partial r} + \frac{\partial u}{\partial z} \frac{\partial v}{\partial z}) 2 π r d r d z = \iint f v 2 π r d r d z (23)

Le domaine de calcul est un rectangle défini par $- R_{m} < r < R_{m}$ et $- Z_{m} < z < Z_{m}$ . La condition limite u=0 est imposée sur les bords de ce rectangle.

Le potentiel présente une symétrie par rapport à l'axe (Oz) : V(z,r)=V(z,-r). La résolution du même problème par la méthode des différences finies (Discrétisation de l'équation de Poisson en géométrie axiale) permet de définir un domaine rectangulaire dont un bord coïncide avec l'axe (Oz) et une condition limite spécifique est définie pour les points de cet axe. Pour la mise en œuvre de la méthode des éléments finis, il est plus simple de définir un domaine symétrique contenant aussi des valeurs négatives de r. Il faudra évidemment définir des conducteurs symétriques par rapport à l'axe (Oz). Le terme r dans l'intégrale doit être remplacé par sa valeur absolue car l'intégrale s'étend aux valeurs négatives de r.

Voici un exemple, avec deux disques formant un condensateur :

conducteursAxial.edp

include "ffmatlib.idp"
real Zm = 100;
real Rm = 100;
real e = 1; //épaisseur des disques
real R = 25; // rayon des disques
real e2 = e/2;
real z1 = -5-e2; // position du premier disque
real z2 = 5+e2; // position du second disque
int C = 10;
border Ca(t=0,1){x=-Zm+2*Zm*t;y=-Rm; label=C;}
border Cb(t=0,1){x=Zm;y=-Rm+2*Rm*t; label=C;}
border Cc(t=0,1){x=Zm-2*Zm*t;y=Rm; label=C;}
border Cd(t=0,1){x=-Zm;y=Rm-2*Rm*t; label=C;}
int C1=11;
border C1a(t=0, 1){ x=z1+e2;  y=R-t*2*R; label=C1;}
border C1b(t=0, 1){ x=z1+e2-e*t;  y=-R;  label=C1;}
border C1c(t=0, 1){ x=z1-e2;  y=-R+t*2*R; label=C1;}
border C1d(t=0, 1){ x=z1-e2+e*t;  y=R; label=C1;} 
int C2=12;
border C2a(t=0, 1){ x=z2+e2;  y=R-t*2*R; label=C2;}
border C2b(t=0, 1){ x=z2+e2-e*t;  y=-R;  label=C2;}
border C2c(t=0, 1){ x=z2-e2;  y=-R+t*2*R; label=C2;}
border C2d(t=0, 1){ x=z2-e2+e*t;  y=R; label=C2;} 
int N=20;
mesh Th = buildmesh(Ca(N)+Cb(N)+Cc(N)+Cd(N)+C1a(10)+C1b(3)+C1c(10)+C1d(3)+C2a(10)+C2b(3)+C2c(10)+C2d(3));

plot(Th, wait=1);
fespace Vh(Th,P1);
Vh u,v;

// x : variable z, y : variable r 
problem Laplace(u,v,eps=1.0e-6) = int2d(Th) ((dx(u)*dx(v)+dy(u)*dy(v))*abs(y))
	+on(C,u=0)
	+on(C1,u=0.5)
	+on(C2,u=-0.5);
Laplace;
int i;
real error = 0.01;
for (i = 0; i < 4; i++) {
	Laplace;
	cout << u[].min << " " << u[].max << endl;
	plot(u,wait=1);
	Th = adaptmesh(Th, u, err=error);
	plot(Th, wait=1);
	u = u;
	error = error/2;
}
Vh Ex,Ey;
Ex = -dx(u);
Ey = -dy(u); 

plot(Th,wait=true); 
plot(u,wait=true);
plot(u,[Ex,Ey]);
savemesh(Th,"conducteursAxial.msh");
ffSaveVh(Th,Vh,"conducteursAxial_vh.txt");
ffSaveData3(u,Ex,Ey,"conducteursAxial_data.txt");

Voici une vue globale du maillage obtenu après adaptation :

freeFemData = FreeFemData('conducteursAxial')
triang = tri.Triangulation(freeFemData.x_array,freeFemData.y_array,freeFemData.triangles)
Zm=100
Rm=100
e=1
e2=e/2
R=25
z1=-5-e2
z2=5+e2	 
fig0,ax0 = subplots(figsize=(10,10))
ax0.set_aspect('equal')
ax0.triplot(triang,color='0.8')
ax0.plot([z1-e2,z1+e2,z1+e2,z1-e2,z1-e2],[-R,-R,R,R,-R],'r-')
ax0.plot([z2-e2,z2+e2,z2+e2,z2-e2,z2-e2],[-R,-R,R,R,-R],'r-')
ax0.set_xlabel('x',fontsize=16)
ax0.set_ylabel('y',fontsize=16)

plotM.pdf

et un détail au voisinage des bords des plaques :

fig0,ax0 = subplots(figsize=(10,10))
ax0.set_aspect('equal')
ax0.triplot(triang,color='0.8')
ax0.plot([z1-e2,z1+e2,z1+e2,z1-e2,z1-e2],[-R,-R,R,R,-R],'r-')
ax0.plot([z2-e2,z2+e2,z2+e2,z2-e2,z2-e2],[-R,-R,R,R,-R],'r-')
ax0.set_xlabel('x',fontsize=16)
ax0.set_ylabel('y',fontsize=16)
ax0.set_xlim(-10,10)
ax0.set_ylim(15,35)

plotN.pdf

Voici un tracé de lignes équipotentielles :

data = freeFemData.readData(3)
V_array = data[0,:]
Ex_array = data[1,:]
Ey_array = data[2,:]

fig1,ax1 = subplots(figsize=(10,10))
ax1.set_aspect('equal')
ct = ax1.tricontour(triang,V_array,levels=20)
fig1.colorbar(ct)
ax1.plot([z1-e2,z1+e2,z1+e2,z1-e2,z1-e2],[-R,-R,R,R,-R],'r-')
ax1.plot([z2-e2,z2+e2,z2+e2,z2-e2,z2-e2],[-R,-R,R,R,-R],'r-')
ax1.set_xlabel('z',fontsize=16)
ax1.set_ylabel('r',fontsize=16)
ax1.set_xlim(-50,50)
ax1.set_ylim(-50,50)

plotO.pdf

Champ électrique sur l'axe (Ox) :

Ex_inter = tri.LinearTriInterpolator(triang, Ex_array)
Ey_inter = tri.LinearTriInterpolator(triang, Ey_array)
x = np.linspace(-Zm,Zm,1000)
y = np.ones(len(x))*0
Ex = Ex_inter(x,y)
figure(figsize=(8,6))
plot(x,Ex)
xlabel('z',fontsize=16)
ylabel('Ez',fontsize=16)
grid()

plotP.pdf

5.d. Problème de magnétostatique en symétrie axiale

L'équation de Poisson pour le potentiel vecteur s'écrit :

Δ \vec{A} = - μ_{0} \vec{J} (24)

ou bien (puisque $d i v \vec{A} = 0$ ) :

\vec{r o t} (\vec{r o t} \vec{A}) = μ_{0} \vec{J} (25)

Pour un problème à symétrie axiale d'axe (Oz), on a $\vec{A} = A_{θ} (r, z) {\vec{u}}_{θ}$ et $\vec{J} = J_{θ} (r, z) {\vec{u}}_{θ}$ . L'équation aux dérivées partielles s'écrit :

\frac{\partial}{\partial r} (\frac{1}{r} \frac{\partial (r A_{θ})}{\partial r}) + \frac{\partial^{2} A_{θ}}{\partial z^{2}} = - μ_{0} J_{θ} (26)

La discrétisation de cette équation pour la méthode de différences finies est exposée dans Discrétisation de l'équation de Poisson vectorielle.

Pour la méthode des éléments finis, nous devons établir une formulation variationnelle. Il faut faire attention au fait que l'opérateur de dérivation seconde appliqué à A_θ dans l'équation ci-desssus n'est pas l'opérateur laplacien. Cette équation est donc différente de celle considérée au paragraphe précédent.

Considérons une fonction test vectorielle $\vec{v}$ . L'application de la formule d'intégration par parties $(5)$ conduit à :

\int_{Ω} v_{x} (\frac{\partial A_{z}}{\partial y} - \frac{\partial A_{y}}{\partial z}) d τ = - \int_{Ω} (A_{z} \frac{\partial v_{x}}{\partial y} - A_{y} \frac{\partial v_{x}}{\partial z}) d τ + \int_{\partial Ω} v_{x} (A_{z} n_{y} - A_{y} n_{z}) d s (27)

La même relation est valable par permutation circulaire de (x,y,z). La somme des trois relations obtenues conduit à :

\int_{Ω} \vec{v} \cdot \vec{r o t} \vec{A} d τ = \int_{Ω} \vec{A} \cdot \vec{r o t} \vec{v} d τ + \int_{\partial Ω} \vec{v} \cdot (\vec{A} \land \vec{n}) d s (28)

Appliquons cette relation à $\vec{r o t} \vec{A}$ à la place de $\vec{A}$ :

\int_{Ω} \vec{v} \cdot (\vec{r o t} \vec{r o t} \vec{A}) d τ = \int_{Ω} \vec{r o t} \vec{A} \cdot \vec{r o t} \vec{v} d τ + \int_{\partial Ω} \vec{v} \cdot (\vec{r o t} \vec{A} \land \vec{n}) d s (29)

La fonction test $\vec{v}$ s'annule sur le bord $\partial Ω$ . On obtient donc la formulation variationnelle suivante :

\int_{Ω} \vec{r o t} \vec{A} \cdot \vec{r o t} \vec{v} d τ = μ_{0} \int_{Ω} \vec{v} \cdot \vec{J} d τ (30)

pour toute fonction $\vec{v}$ de classe C¹ sur $Ω$ qui s'annule sur $\partial Ω$ .

Ce résultat est valable pour un problème de magnétostatique en général. Pour un problème axisymétrique, le rotationnel s'écrit :

\vec{r o t} \vec{A} = - \frac{\partial A_{θ}}{\partial z} {\vec{u}}_{r} + \frac{1}{r} \frac{\partial}{\partial r} (r A_{θ}) {\vec{u}}_{z} = - \frac{\partial A_{θ}}{\partial z} {\vec{u}}_{r} + (\frac{A_{θ}}{r} + \frac{\partial A_{θ}}{\partial r}) {\vec{u}}_{z} (31)

et nous obtenons, pour $\vec{v} = v_{θ} (r, z) {\vec{u}}_{θ}$ :

\iint (\frac{\partial A_{θ}}{\partial z} \frac{\partial v_{θ}}{\partial z} + (\frac{A_{θ}}{r} + \frac{\partial A_{θ}}{\partial r}) (\frac{v}{r} + \frac{\partial v}{\partial r}) 2 π r d r d z = μ_{0} \iint v_{θ} J_{θ} 2 π r d r d z (32)

Le domaine de calcul est défini par $- R_{m} < r < R_{m}$ et $- Z_{m} < z < Z_{m}$ .

Considérons comme exemple une bobine axisymétrique d'axe (Oz). Ce problème est traité par la méthode des différences finies dans Bobine épaisse.

bobineAxial.edp

include "ffmatlib.idp"
real Zm = 100;
real Rm = 100;
real a = 8; // rayon moyen
real e = 4; // épaisseur
real L = 50; // longueur
real e2 = e/2;
real L2 = L/2;
int C = 10;
//border Ca(t=0,1){x=-Zm+2*Zm*t;y=0; label=C;}
//border Cb(t=0,1){x=Zm;y=Rm*t; label=C;}
//border Cc(t=0,1){x=Zm-2*Zm*t;y=Rm; label=C;}
//border Cd(t=0,1){x=-Zm;y=Rm*(1-t); label=C;}

border Ca(t=0,1){x=-Zm+2*Zm*t;y=-Rm; label=C;}
border Cb(t=0,1){x=Zm;y=-Rm+2*Rm*t; label=C;}
border Cc(t=0,1){x=Zm-2*Zm*t;y=Rm; label=C;}
border Cd(t=0,1){x=-Zm;y=Rm-2*Rm*t; label=C;}

int N=20;
mesh Th = buildmesh(Ca(N)+Cb(N)+Cc(N)+Cd(N));

plot(Th, wait=1);
fespace Vh(Th,P1);
Vh u,v;

func real source(real &x, real &y){
	if ((abs(x)<L2)&&(abs(y-a)<e2)) {return 1.0;} 
	else if ((abs(x)<L2)&&(abs(y+a)<e2)) {return -1.0;} 
	else {return 0;}
}
  
// x : variable z, y : variable r 
// u : A_theta
problem Poisson(u,v,eps=1.0e-6) = int2d(Th) ((dx(u)*dx(v)+(u/y+dy(u))*(v/y+dy(v)))*abs(y)) -int2d(Th) (source(x,y)*v*abs(y))
	+on(C,u=0);
	 
	
Poisson; 
int i;
real error = 0.01;
for (i = 0; i < 4; i++) {
	Poisson;
	cout << u[].min << " " << u[].max << endl;
	plot(u,wait=1);
	Th = adaptmesh(Th, u, err=error);
	plot(Th, wait=1);
	u = u;
	error = error/2;
}
Vh Br,Bz; 
Br=-dx(u);  
Bz=u/y+dy(u);

plot(Th,wait=true);  
plot(u,wait=true); 
 
savemesh(Th,"bobineAxial.msh");
ffSaveVh(Th,Vh,"bobineAxial_vh.txt");    
ffSaveData3(u,Br,Bz,"bobineAxial_data.txt");

Voici une vue globale du maillage :

freeFemData = FreeFemData('bobineAxial')
triang = tri.Triangulation(freeFemData.x_array,freeFemData.y_array,freeFemData.triangles)  
Zm=100
Rm=100 
a=8 
e=4
e2=e/2
L=50
L2=L/2
fig0,ax0 = subplots(figsize=(10,10))
ax0.set_aspect('equal')
ax0.triplot(triang,color='0.8')
ax0.plot([-L2,L2,L2,-L2,-L2],[a-e2,a-e2,a+e2,a+e2,a-e2],'r-')
ax0.plot([-L2,L2,L2,-L2,-L2],[-a-e2,-a-e2,-a+e2,-a+e2,-a-e2],'r-')
ax0.set_xlabel('z',fontsize=16)
ax0.set_ylabel('r',fontsize=16)

plotQ.pdf

et un détail du voisinage de la section de la bobine :

fig0,ax0 = subplots(figsize=(10,10))
ax0.set_aspect('equal')
ax0.triplot(triang,color='0.8')
ax0.plot([-L2,L2,L2,-L2,-L2],[a-e2,a-e2,a+e2,a+e2,a-e2],'r-')
ax0.plot([-L2,L2,L2,-L2,-L2],[-a-e2,-a-e2,-a+e2,-a+e2,-a-e2],'r-')
ax0.set_xlabel('z',fontsize=16)
ax0.set_ylabel('r',fontsize=16)
ax0.set_xlim(-L,L)  
ax0.set_ylim(-L,L)

plotR.pdf

Voici les lignes d'isovaleurs de A_θ :

data = freeFemData.readData(3)
Atheta_array = data[0,:]
Br_array = data[1,:] 
Bz_array = data[2,:] 

fig1,ax1 = subplots(figsize=(10,10))
ax1.set_aspect('equal')
ct = ax1.tricontour(triang,Atheta_array,levels=10)
fig1.colorbar(ct)
ax1.plot([-L2,L2,L2,-L2,-L2],[a-e2,a-e2,a+e2,a+e2,a-e2],'r-')
ax1.plot([-L2,L2,L2,-L2,-L2],[-a-e2,-a-e2,-a+e2,-a+e2,-a-e2],'r-')
ax1.set_xlabel('z',fontsize=16)
ax1.set_ylabel('r',fontsize=16)
ax1.set_xlim(-L,L)
ax1.set_ylim(-L,L)

plotS.pdf

et le champ magnétique au voisinage de l'axe de la bobine :

Br_inter = tri.LinearTriInterpolator(triang, Br_array)
Bz_inter = tri.LinearTriInterpolator(triang, Bz_array)
x = np.linspace(-Zm,Zm,1000)
y = np.ones(len(x))*1
Bz = Bz_inter(x,y)
figure(figsize=(8,6))
plot(x,Bz) 
xlabel('z',fontsize=16)
ylabel('Bz',fontsize=16)
grid()

plotT.pdf

et le champ en z=0 en fonction de r :

y = np.linspace(-Rm,Rm,1000)
x = np.ones(len(y))*0
Bz = Bz_inter(x,y)
figure(figsize=(8,6))
plot(y,Bz)
xlabel('r',fontsize=16)
ylabel('Bz',fontsize=16)
xlim(-40,40)
plot([a-e2,a-e2],[0,4],'r-')
plot([a+e2,a+e2],[0,4],'r-')
plot([-a-e2,-a-e2],[0,4],'r-')
plot([-a+e2,-a+e2],[0,4],'r-')
grid()

plotU.pdf

Théoriquement, pour une bobine très longue d'épaisseur e, on doit avoir au centre $B_{z} = μ_{0} J e$ soit ici $B_{z} = e$ puisque le calcul a été fait avec $μ_{0} J = 1$ .

La courbe précédente est obtenue à partir des dérivées calculées dans FreeFEM. Il est aussi possible de construire une fonction d'interpolation de A_θ et de calculer les dérivées dans le script Python :

Atheta_inter = tri.LinearTriInterpolator(triang, Atheta_array)
y = np.linspace(-Rm,Rm,500)
x = np.ones(len(y))*0
(Gz,Gr) = Atheta_inter.gradient(x,y)
Atheta = Atheta_inter(x,y)
Bz = Atheta/y + Gr
figure(figsize=(8,6))
plot(y,Bz)
xlim(-40,40)
plot([a-e2,a-e2],[0,4],'r-')
plot([a+e2,a+e2],[0,4],'r-')
plot([-a-e2,-a-e2],[0,4],'r-')
plot([-a+e2,-a+e2],[0,4],'r-')
grid()
xlabel('r',fontsize=16)
ylabel('Bz',fontsize=16)

plotV.pdf


x = np.linspace(-Zm,Zm,1000)
y = np.ones(len(x))*1
(Gz,Gr) = Atheta_inter.gradient(x,y)
Atheta = Atheta_inter(x,y)
Bz = Atheta/y + Gr
figure(figsize=(8,6))
plot(x,Bz)
xlim(-100,100)
grid()
xlabel('z',fontsize=16)
ylabel('Bz',fontsize=16)

plotW.pdf

La mauvaise qualité de ces courbes est due au fait que le calcul de la composante $B_{z} = \frac{A_{θ}}{r} + \frac{\partial A_{θ}}{\partial r}$ fait intervenir une dérivée, dont l'évaluation est soumise à une augmentation du bruit présent dans $A_{θ}$ . En comparaison, les dérivées calculées sur une grille rectangulaire dans la méthode des différences finies (Bobine épaisse) soit de bien meilleure qualité.

Voici le tracé de $A_{θ}$ en fonction de r pour z=0 :

figure(figsize=(8,6))
y = np.linspace(-Rm,Rm,500)
x = np.ones(len(y))*0
Atheta = Atheta_inter(x,y)
figure(figsize=(8,6))
plot(y,Atheta)
xlabel('r',fontsize=16)
ylabel(r"$A_{\theta}$",fontsize=16)
grid()

plotX.pdf

Calculons $\frac{\partial A_{θ}}{\partial r}$ directement à partir de cette courbe :

from scipy.signal import convolve
dy = y[1]-y[0]
dAtheta_dr = convolve(Atheta,[1,-1],mode='same') / dy
figure(figsize=(8,6))
plot(y,dAtheta_dr)
xlabel('r',fontsize=16)
ylabel(r"$\frac{\partial A_{\theta}}{\partial r}$",fontsize=16)
grid()

plotY.pdf

Bien que la courbe de $A_{θ}$ semble lisse, elle comporte un léger bruit qui est nettement amplifié dans l'opération de dérivation.

Voici la composante $B_{z}$ :

Bz = Atheta/y + dAtheta_dr
figure(figsize=(8,6))
plot(y,Bz)
xlim(-40,40)
plot([a-e2,a-e2],[0,4],'r-')
plot([a+e2,a+e2],[0,4],'r-')
plot([-a-e2,-a-e2],[0,4],'r-')
plot([-a+e2,-a+e2],[0,4],'r-')
grid()
xlabel('r',fontsize=16)
ylabel('Bz',fontsize=16)

plotZ.pdf