Particule chargée dans le champ magnétique d'une bobine
1. Introduction
On considère le mouvement d'une particule chargée au voisinage de l'axe d'une bobine à symétrie axiale, afin de modéliser une lentille magnétique simple ([1]). Les équations différentielles du mouvement sont présentées dans Particule chargée dans un champ magnétique axial. On utilise le potentiel vecteur et le champ magnétique calculés par résolution numérique de l'équation de Poisson de la magnétostatique (Bobine épaisse).
2. Programme de calcul
import math
import numpy
from scipy.interpolate import RectBivariateSpline
from scipy.integrate import odeint
from matplotlib.pyplot import *
Le programme est constituée d'une classe MvtChampBaxial. Le constructeur effectue la lecture des fichiers textes contenant le potentiel vecteur, ses dérivées et le champ magnétique. Ces données sont transformées en fonctions d'interpolation.
class MvtChampBaxial:
def __init__(self,nom):
self.Atheta = numpy.loadtxt(nom+"-Atheta.txt")
self.AthetaDr = numpy.loadtxt(nom+"-AthetaDr.txt")
self.AthetaDz = numpy.loadtxt(nom+"-AthetaDz.txt")
self.Br = numpy.loadtxt(nom+"-Br.txt")
self.Bz = numpy.loadtxt(nom+"-Bz.txt")
self.r = numpy.loadtxt(nom+"-r.txt")
self.z = numpy.loadtxt(nom+"-z.txt")
self.f_A=RectBivariateSpline(self.z,self.r,self.Atheta.transpose())
self.f_ADr=RectBivariateSpline(self.z,self.r,self.AthetaDr.transpose())
self.f_ADz=RectBivariateSpline(self.z,self.r,self.AthetaDz.transpose())
self.f_Bz=RectBivariateSpline(self.z,self.r,self.Bz.transpose())
self.f_Br=RectBivariateSpline(self.z,self.r,self.Br.transpose())
La fonction suivante trace les lignes de champ :
def lignesB(self,color,density):
streamplot(self.z,self.r,self.Bz,self.Br,color=color,density=density)
xlabel("z")
ylabel("r")
La fonction suivante calcule la trajectoire pour un rapport α=q/m donné, une position initiale donnée, des vitesses radiales et axiales données (vitesse angulaire initiale nulle). Le temps maximal et le temps d'échantillonnage sont aussi à fournir. Les équations différentielles sont obtenues par les équations de Hamilton, et utilisent donc le moment angulaire (constante du mouvement), ce qui permet de traiter le cas de trajectoires partant de l'axe.
def trajectoire(self,alpha,r0,z0,pr0,pz0,tmax,te):
ptheta = alpha*r0*self.f_A(z0,r0)
def equation(y,t):
r=y[0]
z=y[1]
A = self.f_A(z,r)
if ptheta==0:
u = -alpha*A
dpr = u*alpha*self.f_ADr(z,r)
dpz = u*alpha*self.f_ADz(z,r)
dtheta = -alpha*(self.f_Bz(z,r)-self.f_ADr(z,r))
else:
u = ptheta/r-alpha*A
dpr = u*(ptheta/(r*r)+alpha*self.f_ADr(z,r))
dpz = u*alpha*self.f_ADz(z,r)
dtheta = ptheta/(r*r)-alpha*(self.f_Bz(z,r)-self.f_ADr(z,r))
return [y[2],y[3],dpr,dpz,dtheta]
t=numpy.arange(0.0,tmax,te)
y = odeint(equation,[r0,z0,pr0,pz0,0.0],t,rtol=1e-6,atol=1e-6)
yt = y.transpose()
rt = yt[0]
zt = yt[1]
thetat = yt[4]
return [rt,zt,thetat]
La fonction suivante réalise la même chose, avec les équations différentielles obtenues par la loi de Newton. La trajectoire ne peut partir de r=0.
def trajectoire_newton(self,alpha,r0,z0,pr0,pz0,tmax,te):
if r0==0.0:
r0 = 1e-10
def equation(y,t):
r=y[0]
z=y[1]
theta=y[2]
Br = self.f_Br(z,r)
Bz = self.f_Bz(z,r)
dpr = r*y[5]*y[5]+alpha*r*y[5]*Bz
dptheta = 1/r*(-2*y[3]*y[5]+alpha*(y[4]*Br-y[3]*Bz))
dpz = -alpha*r*y[5]*Br
return [y[3],y[4],y[5],dpr,dpz,dptheta]
t=numpy.arange(0.0,tmax,te)
y = odeint(equation,[r0,z0,0.0,pr0,pz0,0.0],t,rtol=1e-6,atol=1e-6)
yt = y.transpose()
rt = yt[0]
zt = yt[1]
thetat = yt[2]
return [rt,zt,thetat]
3. Bobine courte
from matplotlib.pyplot import *
import math
import numpy
from mvtChampBaxial import MvtChampBaxial
mvt = MvtChampBaxial("../../elecmag/bobine/bobine-1")
Tracé de la dérivée du potentiel vecteur par rapport à r en z=0 :
figure(figsize=(8,5))
plot(mvt.r,mvt.f_ADr(0.0,mvt.r)[0])
xlabel("z")
ylabel("Atheta")
grid()
figA.pdf
Cette courbe permet de repérer le domaine de validité de l'approximation paraxiale, pour laquelle cette dérivée est constante (Particule chargée dans un champ magnétique axial).
Voici des trajectoires partant d'un point de l'axe, avec une vitesse radiale initiale variable. Les valeurs du champ étant fixées, on ajuste le rapport q/m, ce qui revient à modifier l'intensité du champ. Les trajectoires sont représentées dans le plan méridien tournant (courbes (z,r)).
figure(figsize=(7,7))
mvt.lignesB(color='grey',density=1)
alpha=900.0
pz0=1.0
z0 = -0.7
r0=0.00
tmax=3.0
te=0.01
[r1,z1,theta1] = mvt.trajectoire(alpha,r0,z0,0.01,pz0,tmax,te)
plot(z1,r1)
[r2,z2,theta2] = mvt.trajectoire(alpha,r0,z0,0.02,pz0,tmax,te)
plot(z2,r2)
[r3,z3,theta3] = mvt.trajectoire(alpha,r0,z0,0.05,pz0,tmax,te)
plot(z3,r3)
[r4,z4,theta4] = mvt.trajectoire(alpha,r0,z0,0.1,pz0,tmax,te)
plot(z4,r4)
[r5,z5,theta5] = mvt.trajectoire(alpha,r0,z0,0.15,pz0,tmax,te)
plot(z5,r5)
axis([-1,1,-0.2,0.2])
grid()
figB.pdf
Il y a stigmatisme pour les trajectoires très faiblement inclinées par rapport à l'axe, conformément à la théorie de l'approximation paraxiale. Le point image se trouve à z=0.5. En revanche, les deux rayons les plus inclinés ne convergent pas vers l'image paraxiale.
On trace aussi l'angle du plan méridien en fonction de z :
figure(figsize=(8,5))
plot(z1,theta1)
plot(z2,theta2)
plot(z3,theta3)
plot(z4,theta4)
xlabel("z")
ylabel("theta")
grid()
figC.pdf
4. Bobine longue
mvt = MvtChampBaxial("../../elecmag/bobine/bobine-2")
Tracé de la dérivée du potentiel vecteur par rapport à r en z=0 :
figure(figsize=(8,5))
plot(mvt.r,mvt.f_ADr(0.0,mvt.r)[0])
xlabel("z")
ylabel("Atheta")
grid()
figD.pdf
Cette courbe permet de repérer le domaine de validité de l'approximation paraxiale, pour laquelle cette dérivée est constante (Particule chargée dans un champ magnétique axial). La zone où la dérivée est constante (champ magnétique uniforme) est bien plus large que pour la bobine courte.
Voici des trajectoires partant d'un point de l'axe, avec une vitesse radiale initiale variable. Les valeurs du champ étant fixées, on ajuste le rapport q/m, ce qui revient à modifier l'intensité du champ. Les trajectoires sont représentées dans le plan méridien tournant (courbes (z,r)).
figure(figsize=(7,7))
mvt.lignesB(color='grey',density=1)
alpha=120.0
pz0=1.0
z0 = -0.7
r0=0.00
tmax=3.0
te=0.01
[r1,z1,theta1] = mvt.trajectoire(alpha,r0,z0,0.01,pz0,tmax,te)
plot(z1,r1)
[r2,z2,theta2] = mvt.trajectoire(alpha,r0,z0,0.02,pz0,tmax,te)
plot(z2,r2)
[r3,z3,theta3] = mvt.trajectoire(alpha,r0,z0,0.05,pz0,tmax,te)
plot(z3,r3)
[r4,z4,theta4] = mvt.trajectoire(alpha,r0,z0,0.1,pz0,tmax,te)
plot(z4,r4)
[r5,z5,theta5] = mvt.trajectoire(alpha,r0,z0,0.15,pz0,tmax,te)
plot(z5,r5)
axis([-1,1,-0.2,0.2])
grid()
figE.pdf
Le stigmatisme est bien meilleur qu'avec la bobine courte.
On trace aussi l'angle du plan méridien en fonction de z :
figure(figsize=(8,5))
plot(z1,theta1)
plot(z2,theta2)
plot(z3,theta3)
plot(z4,theta4)
xlabel("z")
ylabel("theta")
grid()
figF.pdf
Voyons ce qui se passe si l'on réduit la vitesse axiale :
alpha=120.0
pz0=0.05
z0 = -0.7
r0=0.00
tmax=50.0
te=0.01
[r1,z1,theta1] = mvt.trajectoire(alpha,r0,z0,0.01,pz0,tmax,te)
figure(figsize=(7,7))
plot(z1,r1)
xlabel("z")
ylabel("r")
axis([-1,1,-0.2,0.2])
grid()
figG.pdf
figure(figsize=(8,5))
plot(z1,theta1)
xlabel("z")
ylabel("theta")
grid()
figH.pdf
Le point image se trouve beaucoup plus proche de l'objet, avant même l'entrée de la bobine. La particule suit un mouvement cyclotronique en suivant les lignes de champ. À l'intérieur de la bobine, le rayon de la trajectoire est pratiquement constant car le champ magnétique est uniforme.
Voyons une réduction de la vitesse axiale (on pourrait aussi augmenter le champ via le rapport q/m).
alpha=120.0
pz0=0.03
z0 = -0.7
r0=0.00
tmax=50.0
te=0.01
[r1,z1,theta1] = mvt.trajectoire(alpha,r0,z0,0.01,pz0,tmax,te)
figure(figsize=(7,7))
plot(z1,r1)
xlabel("z")
ylabel("r")
axis([-1,1,-0.2,0.2])
grid()
figI.pdf
La particule revient vers l'arrière à environ z=-0.3; on obtient ainsi un miroir magnétique.
