Table des matières Python

Aimant permanent

1. Introduction

Ce document présente le calcul du champ magnétique créé par un ou plusieurs aimants permanent et le calcul de la force d'interaction entre deux aimants.

2. Champ magnétique créé par un aimant

2.a. Charges magnétiques

Dans cette partie, on considère le problème général du champ magnétique créé par un aimant permanent, en l'absence de tout courant électrique.

Les matériaux ferromagnétiques peuvent posséder une aimantation en l'absence de source de champ extérieur. Un aimant permanent est défini comme un volume (V) en tout point duquel la densité volumique de moment magnétique M , appelée aimantation, est supposée connue et constante.

Le champ d'aimantation est équivalent à un vecteur densité de courant électrique égal au rotationnel de l'aimantation ([1]). En régime quasi stationnaire, le champ magnétique créé par l'aimant obéit donc à l'équation de Maxwell-Ampère à laquelle ce courant équivalent est ajouté :

B=μ0(J+rotM)(1)

J est la densité de courant électrique, qui est nulle dans le problème considéré ici.

L'équation (1) est une version de l'équation de Maxwell-Ampère (en régime quasi stationnaire) valable dans la matière et faisant intervenir des champs moyens, c'est-à-dire des champ moyennés dans l'espace sur une échelle beaucoup plus grande que l'échelle atomique

La propriété de conservation du flux magnétique reste valable pour le champ moyen :

divB=0(2)

On définit le vecteur excitation magnétique par :

H=Bμ0-M(3)

L'équation (1) s'écrit alors :

rotH=J(4)

et la conservation du flux magnétique :

divH=divM(5)

Le vecteur excitation magnétique H s'exprime en A/m. Il est parfois nommé champ magnétique et dans ce cas le vecteur B est désigné comme le vecteur densité de flux magnétique, qui est une dénomination tout à fait correcte, ou bien induction magnétique, qui est plus discutable. Dans les problèmes de magnétostatique dans le vide, le vecteur H n'a pas d'intérêt et il est d'usage de désigner le vecteur B comme étant le champ magnétique. Nous préférons garder cette dénomination et désigner H par excitation magnétique.

En l'absence de courant électrique, les deux équations vérifiées par l'excitation magnétique sont :

rotH=0divH=divM

Ces deux équations sont similaires à celles vérifiées par le champ électrostatique. Cela nous amène à définir une densité de charge magnétique (par analogie avec la densité volumique de charge électrique) par :

ρm=divM(6)

Plus précisément, si le champ électrostatique est analogue à l'excitation magnétique, ρm et analogue à ρ/ε0.

Soit S la surface fermée qui délimite le volume V de l'aimant. En principe, l'aimantation, qui est un champ moyen défini à l'échelle mésoscopique, devrait varier continûment lorsqu'on passe de l'intérieure de l'aimant à l'extérieur. Cependant, la distance sur laquelle se fait la variation de l'aimantation au voisinage de la surface est négligeable à l'échelle macroscopique. Il peut donc être intéressant, pour des raisons de simplicité du modèle, de considérer que l'aimantation est non nulle dans le volume V et nulle en dehors de ce volume (dans l'air qui entoure l'aimant). En conséquence, il subit une discontinuité sur la surface S et sa divergence n'est donc pas définie sur cette surface. Ce problème est résolu en considérant l'intégrale de la densité de charge magnétique sur le volume :

Vρmdv=VdivMdv=SMndS(7)

n est le vecteur normal sortant. Une densité surfacique de charge magnétique doit donc être introduite sur la surface de l'aimant, définie par :

σm=Mn(8)

Par analogie avec le champ électrostatique, nous pouvons écrire l'excitation magnétique en un point P de l'espace sous la forme d'une intégrale sur le volume de l'aimant plus une intégrale sur sa surface :

H(P)=14πVQP||QP||3ρmdv+14πSQP||QP||3σmdS(9)

Q désigne un point quelconque du volume ou de la surface.

Pour un barreau aimanté, on considère généralement que l'aimantation est uniforme dans le barreau. On supposera donc que l'aimantation est uniforme dans l'aimant. Dans ce cas, il reste seulement l'intégrale sur la surface :

H(P)=14πSQP||QP||3Mn dS(10)

L'excitation en un point P étant connue, le champ magnétique s'en déduit par :

B=μ0(H+M)dansl'aimant(11)B=μ0Hdansl'air (12)

L'équation (10) permet de faire un calcul direct de l'excitation magnétique par un calcul d'intégrale de surface.

Une autre manière d'aborder le problème consiste à introduire un potentiel magnétique. En effet, le rotationnel de l'excitation magnétique étant nul en l'absence de courants électriques, il existe un potentiel Φm tel que :

H=-gradΦm(13)

L'opérateur divergence appliquée à cette équation conduit à l'équation :

ΔΦm=-ρmdansl'aimant, ΔΦm=0dansl'air(14)

Δ est l'opérateur laplacien. Il s'agit de l'équation de Poisson, similaire à l'équation de Poisson pour le potentiel électrostatique. Cette équation n'a de sens que si la charge magnétique est définie, autrement dit si la divergence de l'aimantation est définie en tout point. Il faudra donc adopter un modèle qui respecte cette condition, ce qui exclue l'introduction d'une charge magnétique surfacique.

Il faut aussi écrire les conditions limites à la frontière entre l'air et l'aimant. L'équation divB=0 implique la continuité de la composante normale du champ magnétique. Dans le cas où l'aimantation est discontinue à l'interface aimant-air, si n désigne le vecteur normal dirigé de l'aimant vers l'air, cette condition s'écrit :

(Haimant+M)n=Hairn(15)

La composante normale de l'excitation magnétique a donc une discontinuité égale à l'aimantation sur la surface.

Si l'aimantation varie continûment entre l'aimant et l'air, la continuité de la composante normale du champ magnétique et de l'aimantation implique la continuité de la composante normale de l'excitation magnétique.

L'équation rotH=0 implique la continuité de la composante tangentielle :

n(Haimant-Hair)=0(16)

L'équation de Poisson peut être résolue numériquement avec une condition limite de potentiel nul sur les bords du domaine, qui doit être de très grande taille par rapport à l'aimant. Une méthode de discrétisation par différence finie associée à une méthode de résolution multigrille devrait constituer un moyen très efficace de résoudre numériquement cette équation.

2.b. Barreau aimanté cylindrique

Le barreau aimanté est cylindrique, d'axe (Oz). Sa section droite peut être circulaire ou rectangulaire. On suppose que l'aimantation dans le volume du barreau est uniforme et dirigée selon l'axe (Oz) :

M=Muz(17)

On s'intéresse aux faces sud et nord du barreau, respectivement situées en z=zs et z=zn car c'est dans leur voisinage qu'une charge magnétique est présente.

Calcul direct

Si l'on souhaite effectuer un calcul direct de l'excitation magnétique par la relation (10), il est intéressant de supposer que lorsqu'on traverse la face sud ou la face nord de l'aimant, l'aimantation subit une discontinuité en passant de 0 à M ou inversement. Comme expliqué plus haut, cette discontinuité introduit une charge de surface, ce qui permet de ramener le calcul direct de l'excitation magnétique à un calcul d'intégrale sur une surface. La densité surfacique de charge magnétique sur la face nord est σo=M et sur la face sud o.

aimant-1.svgFigure pleine page

Le calcul de l'excitation magnétique se fait pour chaque face. Par exemple pour la face nord :

H(P)=14πNQP||QP||3M dS(18)

Équation de Poisson

Pour appliquer la méthode de l'équation de Poisson (14), il faut considérer que l'aimantation varie continument à la traversée de la face sud ou nord. Le modèle le plus simple consiste à supposer que la divergence, c'est-à-dire la dérivée de Mz par rapport à z est uniforme sur une épaisseur δ est nulle en dehors.

aimant-2.svgFigure pleine page

La densité volumique de charge magnétique dans la couche d'épaisseur δ de la face nord est :

ρo=Mδ(19)

et l'opposée dans la couche de la face sud. La densité volumique est nulle en dehors de ces deux couches. L'épaisseur δ doit être très petite devant la longueur du barreau.

L'aimantation varie continûment entre l'intérieur de l'aimant et l'extérieur donc l'excitation magnétique varie continûment. La résolution de l'équation de Poisson se fait donc dans un seul domaine.

L'équation de Poisson doit être discrétisée sur une grille. L'épaisseur minimale (et suffisante) de la couche chargée se réduit au nœuds de la grille qui sont sur la face, ce qui finalement revient à considérer que la charge est surfacique.

2.c. Force magnétique

L'énergie magnétique est, par définition, l'énergie qu'il faut fournir pour établir le champ magnétique à partir d'une situation où aucune champ n'est présent. Si la variation du champ magnétique en un point de l'espace est δB , la variation d'énergie magnétique s'écrit ([1]) :

δW=HδBdv(20)

où l'intégrale est étendue à tout l'espace.

Considérons deux aimants cylindriques, notés 1 et 2, et supposons qu'ils sont alignés sur le même axe (Oz). Dans ce cas, ils exercent l'un sur l'autre une force dirigée selon cet axe. Notons δz le déplacement virtuel de l'aimant 2 par rapport à l'aimant 1, ce dernier restant immobile. Au cours de ce déplacement, chaque aimant garde son aimantation. et le champ magnétique varie de δB à cause de ce déplacement, ce qui implique une variation d'énergie δW. Dans l'air qui entoure les aimants, nous avons :

δWair=airHμ0δHdv(21)

Dans un aimant, par exemple l'aimant 1, on a :

δH-δBμ0=δM1=0(22)

puisque son aimantation ne change pas. Dans les aimants, on a donc la même expression de la variation d'énergie que dans l'air :

δWaimant=aimantHμ0δHdv(23)

Finalement, on obtient par intégration l'expression de l'énergie magnétique totale :

Wm=μ02HHdv=μ02||H||2dv(24)

Cette énergie dépend de la position z de l'aimant 2 par rapport à l'aimant 1. La force agissant sur l'aimant 2 s'écrit :

Fz=-dWmdz(25)

Pour calculer cette force, il faudra calculer l'énergie magnétique pour deux positions voisines de l'aimant 2, puis appliquer une différence finie :

Fz-Wm(z+δz)-Wm(z)δz(26)

Plus généralement, il sera possible de calculer la résultante et le moment des actions mécaniques d'un aimant sur l'autre, quel que soit leurs positions et orientations relatives.

Pour effectuer les calculs numériques, il est intéressant d'utiliser des grandeurs sans dimensions. On conviendra que M désigne en fait l'aimantation divisée par une aimantation de référence M0, qui pourrait être par exemple celle d'un des aimants. Le vecteur H est proportionnel à M0. Supposons aussi que les distances soient définies par rapport à une distance de référence 0 (par exemple 1 mm). On peut aussi ne pas prendre en compte la constante μ0 dans le calcul. La valeur de la force obtenue par le calcul numérique est alors sans dimensions et devra être multipliée par :

μ0M0202(27)

Prenons l'exemple d'aimants en Alnico. Le champ magnétique rémanent pour cet alliage est B0=1,3T. Son aimantation est donc M0=B00. Pour 0=1mm on a donc :

μ0M0202=1μ0B0202=1,34N(28)

3. Calcul direct

3.a. Hypothèses

On s'intéresse au cas particulier d'un aimant ou de deux aimants cylindriques de section circulaire, alignés sur le même axe.

deuxAimants.svgFigure pleine page

Un point P est repéré par ses coordonnées cylindriques (r,θ,z). Tout plan contenant l'axe (Oz) est un plan de symétrie de l'aimantation (vecteur axial). En conséquence, c'est aussi un plan de symétrie de l'excitation magnétique (vecteur axial). Il y a par ailleurs invariance par rotation autour de l'axe. Il s'en suit que :

H=Hr(r,z)ur+Hz(r,z)uz(29)Φm(r,z) (30)

L'énergie magnétique s'écrit :

Wm=μ02-dz0dr02πdθ||H||2r=μ0π-dz0dr ||H||2r(31)

Le calcul numérique de cette intégrale doit se faire sur un domaine fini de l'espace autour des deux aimants, assez grand pour que l'énergie en dehors de ce domaine soit négligeable.

3.b. Calcul du champ magnétique

Le calcul direct de l'excitation magnétique consiste, pour chaque face des aimants, à calculer l'intégrale (18). Il faut calculer l'intégrale pour chaque face (2 faces pour un aimant, 4 faces pour deux aimants) et sommer les champs obtenus pour obtenir le champ complet. Voyons comme se fait le calcul de l'intégrale pour une face d'un aimant de rayon a, c'est-à-dire un disque de rayon a. La densité de charge magnétique sur cette face est σm=M ou σm=-M selon qu'il s'agit de la face nord ou de la face sud de l'aimant d'aimantation M.

L'intégrale s'écrit :

H(P)=σm4πdisqueQP||QP||3dS(32)

Pour expliciter cette intégrale, on peut considérer que θ=0. Notons (r11,z1) les coordonnées cylindriques du point Q situé sur une face.

On a :

QP=(r-r1cos(θ1))ux-r1sin(θ1)uy+(z-z1)uz(33)

L'élément de surface s'écrit dS=r11dr1. Les deux composantes de l'excitation magnétique s'écrivent donc :

Hr(r,z)=σm2π0adr10πdθ1 r1(r-r1cos(θ1))(r2+r12-2rr1cos(θ1)+(z-z1)2)32(34)Hz(r,z)=σm2π0adr10πdθ1 r1(z-z1)(r2+r12-2rr1cos(θ1)+(z-z1)2)32(35)

Pour calculer numériquement ces intégrales, on doit effectuer un maillage du demi-disque de rayon a. L'angle θ1 est échantillonné de la manière uniforme :

Δθ= πNθ(36)θ1,i=iΔθi=0,1,Nθ-1(37)

Pour le rayon r1, il est préférable de choisir un échantillonnage qui permette d'avoir des mailles d'aire constante entre le centre et les bords du disque. Celui-ci devrait convenir :

Δr = aNr(38)r1,j=jΔrj=0,1,Nr-1(39)

La figure suivante représente la maille (i,i+1,j,j+1) :

maille.svgFigure pleine page

Son aire est :

ΔS = Δθ2(r1,j+12-r1,j2)=Δθ(Δr)22(j+1-j)=Δθ(Δr)22(40)

Toutes les mailles ont bien la même aire.

Les approximations des intégrales sont obtenues en considérant que l'intégrale sur une maille est égale à l'aire de la maille multipliée par la valeur de la fonction intégrée au centre de la maille.

Hr(r,z)σm2πΔθ(Δr)22i=0Nθ-1j=0Nr(r-r1,j+12cos(θ1,i+12))(r2+r1,j+122-2rr1,j+12cos(θ1,j+12)+(z-z1)2)32(41) Hz(0,z)σm2πΔθ(Δr)22i=0Nθ-1j=0Nr(z-z1)(r2+r1,j+122-2rr1,j+12cos(θ1,j+12)+(z-z1)2)32(42)

Il peut être intéressant de valider ce calcul en comparant Hz à son expression sur l'axe (Oz), qui admet une forme analytique :

Hz(0,z)=2πσm4π(z-z1)0ar1(r12+(z-z1)2)3/2dr1=σm2[z-z1|z-z1|-z-z1(a2+(z-z1)2)1/2](43)

De part et d'autre de la surface, on a :

Hz(0,z1+)=σ2=M2(44)Hz(0,z1-)=-σ2=-M2(45)

ce qui confirme la discontinuité de la composante normale de l'excitation magnétique, égale à l'aimantation M.

Numpy permet de faire des calculs de somme double très rapidement. Pour cela, on utilise tout d'abord la fonction numpy.meshgrid, qui permet de crééer une grille à partir de deux tableaux. Par exemple, pour créer une grille définie par xi=i et yj=j avec 3 points sur chaque axe, on écrit :

import numpy as np
x = [0,1,2]
y = [0,1,2]
xx,yy = np.meshgrid(x,y)
                   
print(xx)
--> array([[0, 1, 2],
       [0, 1, 2],
       [0, 1, 2]])
print(yy)
--> array([[0, 0, 0],
       [1, 1, 1],
       [2, 2, 2]])

Le tableau xx contient les valeurs de x pour les points de la grille, le tableau yy contient les valeurs de y. L'évaluation d'une fonction a deux variables sur cette grille se fait alors de la manière suivante (pour une fonction universelle) :

def f(x,y):
    return x+y
z = f(xx,yy)
                   
print(z)
--> array([[0, 1, 2],
       [1, 2, 3],
       [2, 3, 4]])

La somme des valeurs sur les nœuds de la grille se calcule aisément :

print(z.sum())
--> 18

Les fonctions suivantes effectuent le calcul des sommes (41) et (42) :

def HrHz(r,z,z1,rr,tt,deltaS,M):
    b2 = (z-z1)**2
    cos = np.cos(tt);
    D = (r**2+rr**2-2*r*rr*cos+b2)**1.5
    Hr = (r-rr*cos)/D
    Hz = (z-z1)/D
    Hr = Hr.sum()*deltaS/(2*np.pi)*M
    Hz = Hz.sum()*deltaS/(2*np.pi)*M
    return (Hr,Hz)
    
def grille(Nr,Ntheta,a):
    delta_theta = np.pi/Ntheta
    delta_r = a/np.sqrt(Nr)
    theta1 = (0.5+np.arange(Ntheta))*delta_theta
    r1 = np.sqrt(0.5+np.arange(Nr))*delta_r
    deltaS = delta_theta*delta_r**2*0.5
    rr,tt = np.meshgrid(r1,theta1)
    return (rr,tt,deltaS)
               
def champH(Nr,Ntheta,r,z,z1,a,M):
    (rr,tt,deltaS) = grille(Nr,Ntheta,a)
    return HrHz(r,z,z1,rr,tt,deltaS,M)
                 

Voici une exemple :

Nr = 100
Ntheta = 50
a = 1
z = 1
z1 = 0
r = 0
M = 1
(Hr,Hz) = champH(Nr,Ntheta,r,z,z1,a,M)

                 
print((Hr,Hz))
--> (0.0, 0.14644532315842412)

Pour comparaison, la fonction suivante calcule le champ sur l'axe, défini par (43) :

def Hz_axe(z,z1,a,M):
    return M/2*((z-z1)/abs(z-z1)-(z-z1)/(a**2+(z-z1)**2)**0.5)

Hz = Hz_axe(z,z1,a,M)
                 
print(Hz)
--> 0.14644660940672627

Pour faire une comparaison plus complète, traçons Hz en fonction de z pour r=0 :

from matplotlib.pyplot import *
z = np.linspace(0.01,5,100)
Nr = 100
Ntheta = 50
a = 1
z1 = 0
r = 0
M = 1
(rr,tt,deltaS) = grille(Nr,Ntheta,a)
Hz1 = np.zeros(len(z))
for k in range(len(z)):
    Hr,Hz1[k] = HrHz(r,z[k],z1,rr,tt,deltaS,M)
Hz2 = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz1,"r-",label='approché')
plot(z,Hz2,"k",label='exact')
grid()
xlabel('z')
ylabel('Hz')
legend(loc='upper right')
ylim(0,1)
                 
Hz-axe-1faceHz-axe-1face.pdf

Pour cette taille de grille, le calcul approché de l'intégrale donne un résultat convenable sauf à petite distance du disque. À très petite distance du disque, la fonction à intégrer prend des valeurs très grandes lorsque le point Q se trouve proche de P et elle varie beaucoup lorsqu'il s'éloigne de cette position. Pour améliorer la précision du calcul de l'intégrale, il faut augmenter la finesse de la subdivision de r :


z = np.logspace(-4,-2,10)
Nr = 1000
Ntheta = 50
a = 1
z1 = 0
r = 0
M = 1
(rr,tt,deltaS) = grille(Nr,Ntheta,a)
Hz1 = np.zeros(len(z))
for k in range(len(z)):
    Hr,Hz1[k] = HrHz(r,z[k],z1,rr,tt,deltaS,M)
Hz2 = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz1,"r-",label='approché')
plot(z,Hz2,"k",label='exact')
grid()
xlabel('z')
ylabel('Hz')
legend(loc='upper right')
xscale('log')
                 
Hz-axe-1face-2Hz-axe-1face-2.pdf

Une augmentation d'un facteur 100 du nombre de points selon r permet bien d'augmenter la précision, mais à très courte distance, il y a toujours une divergence de l'intégrale.

Il faut donc adapter le nombre Nr à la distance z et ne pas faire le calcul à une distance trop petite. Il faut trouver un critère pour choisir ce nombre. On peut par exemple considérer que le rayon de la maille centrale doit être de l'ordre de z/p, où p est un paramètre à ajuster. Cette condition conduit à :

Nr=E[(pa|z|)2](46)

mais Nr ne doit pas descendre en dessous d'une certaine valeur.

z = np.linspace(0.01,5,100)
Ntheta = 50
a = 1
z1 = 0
r = 0
M = 1
Hz1 = np.zeros(len(z))
p = 3
Nr_min = 10
for k in range(len(z)):
    Nr = max(int((p*a/abs(z[k]))**2),Nr_min)
    (rr,tt,deltaS) = grille(Nr,Ntheta,a)
    Hr,Hz1[k] = HrHz(r,z[k],z1,rr,tt,deltaS,M)
Hz2 = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz1,"r-",label='approché')
plot(z,Hz2,"k",label='exact')
grid()
xlabel('z')
ylabel('Hz')
legend(loc='upper right')
ylim(0,1)
                 
Hz-axe-1face-3Hz-axe-1face-3.pdf

Il reste à savoir si la méthode fonctionne aussi lorsque le point P est proche du disque mais hors de l'axe. L'aire de la maille qui se trouve au plus proche du point P a la même aire quelle que soit la position de P. On peut donc raisonnablement s'attendre à un calcul approché correct également dans ce cas.

Traçons les composantes Hr et Hz ainsi que la norme de H, en fonction de r, pour une distance z petite :

r = np.linspace(0,2,50)
z = 1e-1
a = 1
z1 = 0
M = 1
Hz1 = np.zeros(len(r))
Hr1 = np.zeros(len(r))
Nr_min = 10
p = 2
Nr = max(int((p*a/abs(z))**2),Nr_min)
(rr,tt,deltaS) = grille(Nr,Ntheta,a)
for k in range(len(r)):    
    Hr1[k],Hz1[k] = HrHz(r[k],z,z1,rr,tt,deltaS,M)
H1 = np.sqrt(Hr1**2+Hz1**2)
Hz2 = np.zeros(len(r))
Hr2 = np.zeros(len(r))
p = 3
Nr = max(int((p*a/abs(z))**2),Nr_min)
(rr,tt,deltaS) = grille(Nr,Ntheta,a)
for k in range(len(r)):    
    Hr2[k],Hz2[k] = HrHz(r[k],z,z1,rr,tt,deltaS,M)
H2 = np.sqrt(Hr2**2+Hz2**2)
figure()
subplot(311)
plot(r,Hr1,'r-')
plot(r,Hr2,'b-')
ylabel('Hr')
grid()
subplot(312)
plot(r,Hz1,'r-')
plot(r,Hz2,'b-')
ylabel('Hz')
grid()
subplot(313)
plot(r,H1,'r')
plot(r,H2,'b-')
ylabel('H')
xlabel('r')
grid()
                  
HrHz-prox-1faceHrHz-prox-1face.pdf

Comme souvent en calcul numérique, le moyen le plus simple de vérifier que le calcul est assez précis est de faire varier le paramètre qui contrôle la précision, ici p, et de comparer les résultats. Les deux courbes précédentes sont calculées pour p=2 et p=3. La valeur p=2 semble suffisante.

Dans l'objectif de tracer des lignes de champ, la méthode de calcul directe semble intéressante, à condition de démarrer les lignes de champ à une distance pas trop petite des faces des aimants. À grande distance, le calcul précis se fait avec une faible complexité.

Le calcul du champ créé par un aimant est fait par la fonction suivante, qui renvoie les composantes de H et de B . L'aimant est caractérisé par son rayon a, les côtes de ses faces zs et zn et son aimantation M.

def champAimant(r,z,a,zs,zn,M,mu0=1):
    p = 3
    Nr_min = 10
    Ntheta = 50
    Nr = Nr_min
    if abs(r)<=a:
        Nr = max(int((p*a/abs(z-zs))**2),Nr_min)
    (rr,tt,deltaS) = grille(Nr,Ntheta,a)
    (Hr1,Hz1) = HrHz(r,z,zs,rr,tt,deltaS,-M)
    Nr = Nr_min
    if abs(r)<=a:
        Nr = max(int((p*a/abs(z-zn))**2),Nr_min)
    (rr,tt,deltaS) = grille(Nr,Ntheta,a)
    (Hr2,Hz2) = HrHz(r,z,zn,rr,tt,deltaS,M)
    Br1 = mu0*Hr1
    Br2 = mu0*Hr2
    if abs(r)<a and z>min(zs,zn) and z<max(zs,zn):
        Bz1 = mu0*(Hz1+M)
        Bz2 = mu0*(Hz2+M)
    else:
        Bz1 = mu0*Hz1
        Bz2 = mu0*Hz2
    return (Hr1+Hr2,Hz1+Hz2,Br1+Br2,Bz1+Bz2)
                    

La fonction suivante permet de dessiner une flèche au milieu d'une ligne.

def fleche(x,y,sens,long,style='k-'):
    n = len(x)//2
    xa = x[n]
    xb = x[n+sens]
    ya = y[n]
    yb = y[n+sens]
    z = (xb-xa)+1j*(yb-ya)
    phi = np.angle(z)
    a = np.pi/2+np.pi/3
    xc = xb+long*np.cos(phi-a)
    yc = yb+long*np.sin(phi-a)
    xd = xb+long*np.cos(phi+a)
    yd = yb+long*np.sin(phi+a)
    plot([xb,xc],[yb,yc],style)
    plot([xb,xd],[yb,yd],style)
                    

Le tracé d'une ligne de champ de l'excitation magnétique se fait en partant d'un point à proximité d'une face, jusqu'à un point à proximité d'une face ou bien jusqu'à sortie d'un rectangle centré sur l'aimant. Il faut remarquer que le choix de la valeur de M et de μ0 n'a aucune influence sur les lignes de champ.

def ligneH(a,zs,zn,M,ri,zi,sens,rmax,zmax,dmin):
    ds = 0.01*sens
    r = ri
    z = zi
    ligne_z = []
    ligne_r = []
    continuer = True
    while continuer:
        ligne_z.append(z)
        ligne_r.append(r)
        (Hr,Hz,Br,Bz) = champAimant(r,z,a,zs,zn,M)
        H = np.sqrt(Hr**2+Hz**2)
        r += Hr/H*ds                                    
        z += Hz/H*ds       
        if (abs(z-zs) < dmin and abs(r)<a) or (abs(z-zn)< dmin and abs(r)<a) or (abs(r)>rmax) or (abs(z)>zmax) : continuer=False
    return (np.array(ligne_r),np.array(ligne_z))
a=1
M = 1 
zs=-2
zn=2
ds = 0.05
figure(figsize=(8,8))
zmax = 10
rmax = 10
longfleche = 0.2
dmin = 0.1 # distance minimale d'approche des faces
for ri in [0,0.2,0.4,0.6,0.8,1.0]:
    sens = 1
    ligne_r,ligne_z = ligneH(a,zs,zn,M,ri,zn+dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
    sens = -1
    ligne_r,ligne_z = ligneH(a,zs,zn,M,ri,zs-dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
    sens = 1
    ligne_r,ligne_z = ligneH(a,zs,zn,M,ri,zn-dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
axis('equal')
xlabel('z')
ylabel('r')
xlim(-zmax,zmax)
ylim(-rmax,rmax)
plot([-zs,zs,zs,-zs,-zs],[a,a,-a,-a,a],'r-')
grid()
title('Lignes de H')
                     
aimant-ligneHaimant-ligneH.pdf
def ligneB(a,zs,zn,M,ri,zi,sens,rmax,zmax,dmin):
    ds = 0.01*sens
    r = ri
    z = zi
    ligne_z = []
    ligne_r = []
    continuer = True
    while continuer:
        ligne_z.append(z)
        ligne_r.append(r)
        (Hr,Hz,Br,Bz) = champAimant(r,z,a,zs,zn,M)
        B = np.sqrt(Br**2+Bz**2)
        r += Br/B*ds
        z += Bz/B*ds
        if (abs(z-zs) < dmin and abs(r)<a) or (abs(z-zn)< dmin and abs(r)<a) or (abs(r)>rmax) or (abs(z)>zmax) : continuer=False
    return (np.array(ligne_r),np.array(ligne_z))
    
figure(figsize=(8,8))
for ri in [0,0.2,0.4,0.6,0.8,1.0]: 
    sens = 1
    ligne_r,ligne_z = ligneB(a,zs,zn,M,ri,zn+dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
    sens = -1 
    ligne_r,ligne_z = ligneB(a,zs,zn,M,ri,zs-dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
    sens = -1
    ligne_r,ligne_z = ligneB(a,zs,zn,M,ri,zn-dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
axis('equal')
xlabel('z')
ylabel('r')
xlim(-zmax,zmax)
ylim(-rmax,rmax)
plot([-zs,zs,zs,-zs,-zs],[a,a,-a,-a,a],'r-')
grid()
title('Lignes de B')
                     
aimant-ligneBaimant-ligneB.pdf

3.c. Calcul de l'énergie magnétique

Le calcul de l'énergie magnétique par l'intégrale (31) nécessite de définir une grille de points sur un domaine assez large centré sur l'aimant. La difficulté de ce calcul approché d'intégrale vient du fait que la norme du champ varie beaucoup lorsqu'on s'éloigne de l'aimant. Elle est d'ailleurs maximale à proximité des faces sud et nord, là où le calcul du champ pose problème.

Pour réduire le nombre de nœuds sans compromettre la précision du résultat, nous devons définir une grille dont les mailles sont d'autant plus fines qu'on est proche de l'aimant. Soit R le rayon du domaine d'intégration. L'échantillonnage de r peut être fait de la manière suivante :

Δr=RNrn(47)rj = jnΔrj=0,1,Nr-1(48)

n est un nombre strictement supérieur à 1, par exemple 2 ou 3.

Il est beaucoup plus difficile de faire un échantillonnage irrégulier de z, car celui-ci devrait dépendre de la position de l'aimant sur l'axe et serait complexe à définir avec deux aimants. On opte donc, dans un premier temps, pour l'échantillonnage suivant, où (Zmin,Zmax) désigne l'intervalle centré sur l'aimant où on fera le calcul de l'intégrale.

Δz=Zmax-ZminNz(49)zi = Zmin+iΔzi=0,1,Nz-1(50)

L'aire de la maille comprise entre les rayons ri et ri+1, et entre zj et zj=1 est :

ΔSj=(rj+1-rj)Δz=((j+1)n-jn)ΔrΔz(51)

Cette aire est d'autant plus grande que j est grand et que l'exposant n est grand.

Le module de H doit être évalué au milieu de cette maille, c'est-à-dire en :

zm,i=zi+Δz2(52)rm,j=rj+rj+12=(jn+(j+1)n)Δr2(53)

L'approximation de l'intégrale (31) s'écrit :

Wmμ0πi=0Nzj=0Nrrm,i||H(rm,i,zm,i)||2ΔSj(54)

Nous avons vu qu'il est en pratique impossible (pour des raisons de temps de calcul) de calculer le champ à une distance des faces inférieure à δz, estimée à 0,01 pour cette implémentation en Python. Sachant que le champ est évalué au milieu des mailles, il est judicieux de définir la subdivision de l'axe z de telle sorte que les faces soient précisément sur des nœuds de la grille et de choisir Δz au moins égal à 2δz. Pour y parvenir simplement, nous supposons que zs,zn,Zmin,Zmax sont des nombres entiers, ce qui est une convention facile à respecter si les distances sont exprimées en mm (il suffira de diviser par un déplacement en mm pour calculer la force). Dans ce cas, si nous choisissons Nz multiple de Zmax-Zmin, l'abscisse (z) de chaque face coïncide bien avec une des valeurs zj. Si zs,zn multipliés par 10 sont des nombres entiers (valeurs définies aux dixième de mm), il faut que Nz soit multiple de 10(Zmax-Zmin).

La fonction suivante calcule l'énergie magnétique pour un aimant :

def energieAimant(a,zs,zn,M,Zmin,Zmax,R,Nz,Nr,n,mu0=1):
    Delta_z = (Zmax-Zmin)/Nz
    Delta_r = R/Nr**n
    W = 0
    for i in range(Nz):
        for j in range(Nr):
            rm = (j**n+(j+1)**n)*Delta_r/2
            zm = Zmin + (i+0.5)*Delta_z
            DeltaS = ((j+1)**n-j**n)*Delta_r*Delta_z
            (Hr,Hz,Br,Bz) = champAimant(rm,zm,a,zs,zn,M,mu0) 
            W += rm*(Hz**2+Hr**2)*DeltaS
    W *= mu0*np.pi
    return W
    
                          

Voici des exemples de calcul. La difficulté est de choisir Nr et Nz assez grand tout en gardant un temps de calcul raisonnable. Nous essayons plusieurs valeurs de Nz car c'est le paramètre critique : une grande valeur est nécessaire pour prendre en compte les zones proches des faces, là où le champ est le plus intense.

a=1
M = 1 
zs=-2
zn=2
Nz = 300
Nr = 100
W = energieAimant(a,zs,zn,M,-10,10,10,Nz,Nr,3)
                          
print(W)
--> 1.1331510305163626
Nz = 400
Nr = 100
W = energieAimant(a,zs,zn,M,-10,10,10,Nz,Nr,3)                          
                         
print(W)
--> 1.130434482130545
Nz = 500
Nr = 150
W = energieAimant(a,zs,zn,M,-10,10,10,Nz,Nr,3)                          
                         
print(W)
--> 1.1264871005928432
Nz = 500
Nr = 200
W = energieAimant(a,zs,zn,M,-10,10,10,Nz,Nr,3)                          
                         
print(W)
--> 1.126523385572853
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.