Table des matières MathematicaPython

Boucles circulaires coaxiales

1. Introduction

Le calcul du champ magnétique créé par une spire circulaire de courant est effectué. Un programme python est présenté, qui effectue le calcul pour un ensemble de spires coaxiales et calcule les lignes de champ. Ce programme permet également de calculer la trajectoire d'une particule chargée en coordonnées cylindriques.

2. Champ magnétique

Afin de calculer le champ magnétique créé par des bobines de section circulaire, on s'intéresse tout d'abord au champ créé par une boucle de courant de rayon a, parcourue par un courant quasi-stationnaire I(t). Le calcul du potentiel vecteur en coordonnées sphériques est fait dans [1]. Afin de superposer facilement plusieurs boucles, le calcul doit être fait en coordonnées cartésiennes.

Soit Oxyz un repère orthogonal où O est le centre de la boucle et Oz un axe perpendiculaire au plan de la boucle.

figureA.svgFigure pleine page

Soit M(x,y,z) le point où on calcule le champ et P un point de la boucle repéré par son angle θ. La tangente unitaire en ce point est :

Le potentiel vecteur au point M est donné par l'intégrale :

Les composantes du potentiel vecteur sont donc :

Le champ magnétique est obtenu en calculant le rotationnel :

Le vecteur étant contenu dans le plan méridien (plan contenant M et l'axe Oz) et compte tenu de l'invariance par rotation, on peut se limiter à calculer le champ dans le plan y=0, ce qui donne :

Avec Mathematica, on peut évaluer les deux intégrales précédentes :

ibx=Integrate[z*Cos[theta]/(x^2+z^2+a^2-2*a*x*Cos[theta])^(3/2),theta];
bx=Evaluate[(ibx/.{theta->Pi})-(ibx/.{theta->0})]
ibz=Integrate[(-x*Cos[theta]+a)/(x^2+z^2+a^2-2*a*x*Cos[theta])^(3/2),theta];
bz=Evaluate[(ibz/.{theta->Pi})-(ibz/.{theta->0})]

Les expressions obtenues font intervenir les intégrales elliptiques complètes de première et seconde espèces. Voici la définition de ces intégrales :

Les lignes de champ magnétique sont obtenues par résolution numérique du système différentiel suivant :

La variable t est l'abscisse curviligne le long de la ligne de champ.

Le potentiel vecteur en coordonnées cylindriques (r=x) est obtenu par l'intégrale suivante :

On obtient avec Mathematica :

A=Integrate[Cos[theta]/(x^2+z^2+a^2-2*a*x*Cos[theta])^(1/2),theta];
Atheta=Evaluate[(A/.{theta->Pi})-(A/.{theta->0})]

Il peut être utile de calculer la dérivée suivante :

ADx=Integrate[(a*Cos[theta]-x)*Cos[theta]/(x^2+z^2+a^2-2*a*x*Cos[theta])^(3/2),theta];
ADx=Evaluate[(ADx/.{theta->Pi})-(ADx/.{theta->0})]

3. Programme python

La classe Spire permet de calculer le champ créé par une spire. Le constructeur a pour arguments le rayon de la spire, sa position sur l'axe Oz et l'intensité du courant.

SpiresCoaxiales.py
import math 
import numpy
from scipy.integrate import odeint
from scipy.special import ellipk, ellipe
from matplotlib.pyplot import *

            
class Spire:
    def __init__(self,a,zs,i):
        self.a = a
        self.a2 = a*a
        self.a4 = self.a2*self.a2
        self.zs = zs
        self.i = i
    def champB(self,x,z):
        if abs(x) < 1e-8:
            x=0
        if x>0:
            sx = 1
            x = -x
        else:
            sx = -1
        z = z-self.zs
        x2 = x*x
        z2 = z*z
        r2 = x2+z2
        b1 = self.a2+ r2
        b2 = 2*x*self.a
        b3 = b1+b2
        b4 = b1-b2
        b5 = -2*b2/b4
        b6 = math.sqrt(b3/b4)*self.i 
        rb3 = math.sqrt(b3)
        b7 = self.a*b3*rb3
        b8 = self.a4-self.a2*(x2-2*z2)+z2*(x2+z2)
        b9 = (self.a2+z2)*b3
        e = ellipe(b5)
        k = ellipk(b5)
        bz = b6*((self.a2-r2)*e+b3*k)/b7
        if x==0:
            bx = 0.0
            Atheta = 0.0
            Adx = bz/2
        else:
            bx = -sx*z/x*b6*(b1*e-b3*k)/b7
            Atheta = -sx*b6/x*(-b4*e+(self.a2+r2)*k)/(self.a*rb3)
            Adx = b6/x2*(b8*e-b9*k)/b7
        return [bx,bz,Atheta,Adx]
            

La classe SystemeSpires effectue le calcul du champ et des lignes de champ pour plusieurs spires. Les arguments du constructeur sont les bornes du domaine utilisé pour le calcul des lignes de champ.

class SystemeSpires:
    def __init__(self,xmin,xmax,zmin,zmax):
        self.xmin = xmin
        self.xmax = xmax
        self.zmin = zmin
        self.zmax = zmax
        self.spires = []
        self.n = 0
    def bornes(self,xmin,xmax,zmin,zmax):
        self.xmin = xmin
        self.xmax = xmax
        self.zmin = zmin
        self.zmax = zmax
    def ajouter(self,spire):
        self.spires.append(spire)
        self.n += 1
    def champB(self,x,z):
        bx = 0.0
        bz = 0.0
        A = 0.0
        Adx = 0.0
        for spire in self.spires:
            b = spire.champB(x,z)
            bx += b[0]
            bz += b[1]
            A += b[2]
            Adx += b[3]
        return [bx,bz,A,Adx]
            

La fonction suivante calcul la composante axiale pour un tableau de valeurs de z

    def Bz_z(self,x,z):
        nz = len(z)
        bz = numpy.zeros(nz)
        for k in range(nz):
            b = self.champB(x,z[k])
            bz[k] = b[1]
        return bz
            

Voici des fonctions analogues pour d'autres composantes et coordonnées :

    def Bz_x(self,x,z):
        nx = len(x)
        bz = numpy.zeros(nx)
        for k in range(nx):
            b = self.champB(x[k],z)
            bz[k] = b[1]
        return bz
    def Bx_z(self,x,z):
        nz = len(z)
        bx = numpy.zeros(nz)
        for k in range(nz):
            b = self.champB(x,z[k])
            bx[k] = b[0]
        return bx
    def Bx_x(self,x,z):
        nx = len(x)
        bx = numpy.zeros(nx)
        for k in range(nx):
            b = self.champB(x[k],z)
            bx[k] = b[0]
        return bx
    def A_x(self,x,z):
        nx = len(x)
        A = numpy.zeros(nx)
        for k in range(nx):
            b = self.champB(x[k],z)
            A[k] = b[2]
        return A
    def Adx_x(self,x,z):
        nx = len(x)
        Adx = numpy.zeros(nx)
        for k in range(nx):
            b = self.champB(x[k],z)
            Adx[k] = b[3]
        return Adx
            

Les deux fonctions suivantes renvoient les composantes du champ (avec le potentiel vecteur) sur un domaine, sous forme de deux tableaux à deux dimensions :

    def B(self,x,z):
        nx = len(x)
        nz = len(z)
        bx = numpy.zeros((nx,nz))
        bz = numpy.zeros((nx,nz))
        A = numpy.zeros((nx,nz))
        for i in range(nz):
            for j in range(nx):
                b = self.champB(x[j],z[i])
                bx[j][i] = b[0]
                bz[j][i] = b[1]
                A[j][i] = b[2]
        return [bx,bz,A]
            

La fonction suivante calcule une ligne de champ en partant d'un point :

    def ligneB(self,x0,z0,sens):
        def equation(y,t):
            b = self.champB(y[0],y[1])
            nb = math.sqrt(b[0]*b[0]+b[1]*b[1])
            return [sens*b[0]/nb,sens*b[1]/nb]
        tmax = (self.xmax-self.xmin)
        te = 1.0*tmax/100
        xarray = numpy.array([x0])
        zarray = numpy.array([z0])
        x=x0
        z=z0
        t = 0.0
        while t < tmax:
            y = odeint(equation,[x,z],[0,te],rtol=1e-5,atol=1e-5)
            x = y[1][0]
            z = y[1][1]
            if x<self.xmin or x>self.xmax or z<self.zmin or z>self.zmax:
                break
            xarray = numpy.append(xarray,x)
            zarray = numpy.append(zarray,z)
            t += te
        return [xarray,zarray]
            
            

La fonction suivante trace des lignes de champ pour une liste de points initiaux. Pour chaque point, les deux lignes symétriques par rapport à l'axe sont tracées.

    def plot_lignes(self,points,style):
        for p in points:
            [x,z] = self.ligneB(p[0],p[1],1)
            plot(z,x,style)
            [x,z] = self.ligneB(p[0],p[1],-1)
            plot(z,x,style)
            [x,z] = self.ligneB(-p[0],p[1],1)
            plot(z,x,style)
            [x,z] = self.ligneB(-p[0],p[1],-1)
            plot(z,x,style)
            for s in self.spires:
                plot([s.zs,s.zs],[-s.a,s.a],linestyle=" ",marker="o",markersize=4,color="black")
            

La fonction suivante calcule la trajectoire d'une particule chargée en coordonnées cylindriques. Les équations du mouvement sont présentées dans Particule chargée dans un champ magnétique axial. Les arguments sont le rapport q/m de la particule, la position et la vitesse initiale (radiales et axiales), le temps maximal du calcul et la période d'échantillonnage. La vitesse angulaire initiale est nulle.

    def trajectoire(self,alpha,r0,z0,dr0,dz0,tmax,te):
        b = self.champB(r0,z0)
        ptheta = r0*alpha*b[2]
        def equation(y,t):
            r = y[0]
            theta = y[1]
            z = y[2]
            b = self.champB(r,z)
            if ptheta==0:
                u = -alpha*b[2]
                dpr = u*alpha*b[3]
                dpz = u*alpha*(-b[0])
                dtheta = -alpha*(b[1]-b[3])
            else:
                r2 = r*r
                u = ptheta/r-alpha*b[2]
                dpr = u*(ptheta/r2+alpha*b[3])
                dpz = u*alpha*(-b[0])
                dtheta = ptheta/r2-alpha*(b[1]-b[3])
            return [y[3],dtheta,y[4],dpr,dpz]
        rarray = numpy.array([r0])
        zarray = numpy.array([z0])
        thetarray = numpy.array([0.0])
        r = r0
        theta = 0.0
        z = z0
        dtheta = 0.0
        dr = dr0
        dz = dz0
        t = 0.0
        while t < tmax:
            y = odeint(equation,[r,theta,z,dr,dz],[0,te],rtol=1e-5,atol=1e-5)
            r = y[1][0]
            theta = y[1][1]
            z = y[1][2]
            dr = y[1][3]
            dz = y[1][4]
            if r<self.xmin or r>self.xmax or z<self.zmin or z>self.zmax:
                break
            rarray = numpy.append(rarray,r)
            zarray = numpy.append(zarray,z)
            thetarray = numpy.append(thetarray,theta)
        return [rarray,thetarray,zarray]
            

4. Exemple

On considère un exemple avec une spire.

from SpiresCoaxiales import *
import numpy
from matplotlib.pyplot import *

systeme = SystemeSpires(-5.0,5.0,-5.0,5.0)
systeme.ajouter(Spire(1.0,0.0,1.0))
z = numpy.arange(-10.0,10.0,0.01)
bz = systeme.Bz_z(0.0,z)
figure()
plot(z,bz)
xlabel("z")
ylabel("Bz")
grid()
            
figAfigA.pdf
x = numpy.arange(-3.0,3.0,0.01)
bz = systeme.Bz_x(x,0.1)
figure()
plot(x,bz)
xlabel("x")
ylabel("Bz")
grid()
            
figBfigB.pdf
bx = systeme.Bx_x(x,1.0)
figure()
plot(x,bx)
xlabel("x")
ylabel("Bx")
grid()
            
figCfigC.pdf
bx = systeme.A_x(x,0.0)
figure()
plot(x,bx)
xlabel("x")
ylabel("A")
grid()
            
fig1fig1.pdf

Tracé de lignes de champ :

figure(figsize=(7,7))
systeme.plot_lignes([[0.0,0.0],[0.1,0.0],[0.3,0.0],[0.5,0.0],[0.7,0.0],[0.9,0.0]],'b')
axis([-5.0,5.0,-5.0,5.0])
xlabel('z')
ylabel('r')
            
figDfigD.pdf

Trajectoires d'une particule chargée dans le plan méridien tournant :

figure()
tmax=200.0
te=0.1
alpha = 1.0
r0 = 0.00
z0 = -4.0
dz0 =2.0
[r,theta,z] = systeme.trajectoire(alpha,r0,z0,0.05,dz0,tmax,te)
plot(z,r,'r')
[r,theta,z] = systeme.trajectoire(alpha,r0,z0,0.02,dz0,tmax,te)
plot(z,r,'g')
[r,theta,z] = systeme.trajectoire(alpha,r0,z0,0.1,dz0,tmax,te)
plot(z,r,'b')
axis([-5.0,5.0,-1.0,1.0])
xlabel('z')
ylabel('r')
grid()
            
figEfigE.pdf
Références
[1]  J.D. Jackson,  Classical Electrodynamics,  (John Wiley Sons, 1975)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.