Aimant permanent : utilisation des intégrales elliptiques

1. Introduction

Ce document présente le calcul du champ magnétique créé par un ou plusieurs aimants permanent. La géométrie est axisymétrique et les aimants sont coaxiaux.

Chaque aimant, possédant une aimantation uniforme, est équivalent à deux faces possédant chaune une charge magnétique surfacique uniforme. Le calcul de l'excitation magnétique $\overrightarrow{H}$ créé par une face d'un aimant se fait par une intégrale double sur cette face, comme expliqué dans Aimant permanent.

Dans Aimant permanent, le calcul de cette intégrale double est fait entièrement de manière numérique, à partir d'un maillage de la face de l'aimant.

Nous présentons ici un calcul de l'intégrale reposant sur une expression de l'intégration par rapport à l'angle, qui fait intervenir des intégrales elliptiques complètes, ce qui permet de ramener le calcul de l'intégrale double à une intégrale simple. L'algorithme de Carlson [1] pour le calcul numérique des intégrales elliptiques permet réduire le temps de calcul par rapport au calcul numérique de l'intégrale double.

2. Calcul du champ magnétique

On considère un ou deux aimants en géométrie axisymétrique, chacun étant caractérisé par une aimantation uniforme. Un point P de l'espace est repéré par ses coordonnées cylindriques (r,z).

Figure pleine page

L'excitation magnétique en un point P est calculé par l'intégrale suivante pour chaque face des aimants :

$$\overrightarrow{H}\left(P\right)=\frac{{\sigma }_m}{4\pi }{\iint }_{disque}\frac{\overrightarrow{QP}}{\left|\right|\overrightarrow{QP}|{|}^3}dS\text{(1)}$$

où σ_m est la densité de charge magnétique de la face, égale à σ_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.

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

On a :

$$\overrightarrow{QP}=(r-r_{1}cos\left({\theta }_1\right)){\overrightarrow{u}}_x-r_{1}sin\left({\theta }_1\right){\overrightarrow{u}}_y+(z-z_1){\overrightarrow{u}}_z\text{(2)}$$

L'élément de surface s'écrit dS=r₁dθ₁dr₁. Les deux composantes de l'excitation magnétique s'écrivent donc :

$$\begin{array}{ccc}\hfill & H_r(r,z)=\frac{{\sigma }_m}{4\pi }{\int }_0^{a}d{r}_1{\int }_0^{2\pi }d{\theta }_1\frac{r_1(r-r_{1}cos\left({\theta }_1\right))}{(r^2+r_1^2-2r{r}_{1}cos\left({\theta }_1\right)+{(z-z_1)}^2{)}^{\frac{3}{2}}}\hfill & \hfill \text{(3)}\\ \hfill & H_z(r,z)=\frac{{\sigma }_m}{4\pi }{\int }_0^{a}d{r}_1{\int }_0^{2\pi }d{\theta }_1\frac{r_1(z-z_1)}{(r^2+r_1^2-2r{r}_{1}cos\left({\theta }_1\right)+{(z-z_1)}^2{)}^{\frac{3}{2}}}\hfill & \hfill \text{(4)}\end{array}$$

L'intégrale par rapport à l'angle θ₁ peut s'exprimer à l'aide des intégrales elliptiques de première et seconde espèces, définies par :

$$\begin{array}{cc}\hfill & K\left(x\right)={\int }_0^{\pi /2}\frac{d\theta }{\sqrt{1-x{sin}^2\theta }}\hfill \\ \hfill & E\left(x\right)={\int }_0^{\pi /2}\sqrt{1-x{sin}^2\theta }d\theta \hfill \end{array}$$

Pour exprimer ces intégrales, définissons les expressions suivantes :

$$\begin{array}{ccc}\hfill & b=z-z_1\hfill & \hfill \text{(5)}\\ \hfill & A={(r-r_1)}^2+b^2\hfill & \hfill \text{(6)}\\ \hfill & B={(r+r_1)}^2+b^2\hfill & \hfill \text{(7)}\\ \hfill & C=r^2-r_1^2-b^2\hfill & \hfill \text{(8)}\\ \hfill & x=-\frac{4r{r}_1}{A}\hfill & \hfill \text{(9)}\end{array}$$

Voici les expressions des intégrales sur la variable θ₁, obtenues grace au logiciel de calcul formel Mathematica :

$$\begin{array}{ccc}\hfill & f_r=\frac{2r_1}{r}\frac{CE\left(x\right)+BK\left(x\right))}{B\sqrt{A}}sir>0\hfill & \hfill \text{(10)}\\ \hfill & f_r=0sir=0\hfill & \hfill \text{(11)}\\ \hfill & f_z=\frac{4b{r}_{1}E\left(x\right)}{B\sqrt{A}}\hfill & \hfill \text{(12)}\end{array}$$

Les intégrales elliptiques sont calculées au moyen des fonctions ellipe et ellipk du module scipy.special. Ces deux fonctions implémentent l'algorithme itératif de Carlson, qui converge nettement plus vite qu'une méthode de calcul de type méthode des trapèzes.

import numpy as np
from scipy.special import ellipk, ellipe
from matplotlib.pyplot import *

def frfz(r,z,r1,z1):
    b = z-z1
    b2 = b*b
    A = (r-r1)**2+b2
    B = (r+r1)**2+b2
    C = r**2-r1**2-b2
    x = -4*r*r1/A
    E = ellipe(x)
    K = ellipk(x)
    D = B*np.sqrt(A)
    if r!=0:
        fr =  2*r1/r*(C*E+B*K)/D
    else:
        fr=0
    fz = 4*b*r1*E/D
    return fr,fz

frfz = np.frompyfunc(frfz,4,2)
                   

Voici un exemple :

z1 = 0
a = 1
z = 0.1
r = 0.1
N = 100
r1 = np.linspace(0,a,N)
fr,fz = frfz(r,z,r1,z1)
figure()
subplot(211)
plot(r1,fr)
ylabel('fr')
grid()
subplot(212)
plot(r1,fz)
ylabel('fz')
xlabel('r1')
grid()
                   

frfz-1.pdf

et un autre à une distance plus petite du disque :

z1 = 0
a = 1
z = 0.01
r = 0.1
N = 500
r1 = np.linspace(0,a,N)
fr,fz = frfz(r,z,r1,z1)
figure()
subplot(211)
plot(r1,fr)
ylabel('fr')
grid()
subplot(212)
plot(r1,fz)
ylabel('fz')
xlabel('r1')
grid()
                   

frfz-2.pdf

Le calcul de l'intégrale sur la variable r₁ peut se faire par la méthode des trapèzes :

dr = r1[1]-r1[0]
Iz = (np.sum(fz[1:N-1])+0.5*(fz[0]+fz[N-1]))*dr
Ir = (np.sum(fr[1:N-1])+0.5*(fr[0]+fr[N-1]))*dr
M = 1
Hz = M/(4*np.pi)*Iz
Hr = M/(4*np.pi)*Ir
                   

print((Hz,Hr))
--> (0.49496075263083483, 0.02507405131733448)

Voici une fonction qui effectue le calcul du champ en un point de coordonnées (r,z) :

def HrHz(r,z,z1,a,M,N):
    r1 = np.linspace(0,a,N)
    dr = r1[1]-r1[0]
    fr,fz = frfz(r,z,r1,z1)
    Hz = M/(4*np.pi)*(np.sum(fz[1:N-1])+0.5*(fz[0]+fz[N-1]))*dr
    Hr = M/(4*np.pi)*(np.sum(fr[1:N-1])+0.5*(fr[0]+fr[N-1]))*dr
    return Hr,Hz
                   

Voici par exemple le tracé du champ sur l'axe, qui peut être comparé à l'expression exacte :

P = 100
z = np.linspace(1e-3,2,P)
Hz = np.zeros(P)
r = 0
z1 = 0
a = 1
M = 1
N = 500
for k in range(P):
    hr,hz = HrHz(r,z[k],z1,a,M,N)
    Hz[k] = hz
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_exact = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz,'b')
plot(z,Hz_exact,'r')
grid()
xlabel('z')
ylabel('Hz')
                   

Hz-1.pdf

Le calcul numérique de l'intégrale est incorrect lorsque le point est très proche du disque, car la fonction à intégrer présente dans ce cas des variations très rapides au voisinage de r₁=r. Dans le cas ci-dessus, la valeur N=500 n'est pas suffisante pour z-z₁=10^-3. À l'inverse, pour les points plus éloignés du disque il devrait être possible de réduire la valeur de N. Il est donc judicieux de mettre en place une méthode d'ajustement de N en fonction d'une tolérance. La fonction suivante répète le calcul en augmentant N d'un facteur deux jusqu'à ce que la variation relative de la norme du champ soit inférieure à la tolérance :

def HrHz_adapt(r,z,z1,a,M,Nmin,tol):
    N = Nmin
    Hr1,Hz1 = HrHz(r,z,z1,a,M,N)
    H1 = np.sqrt(Hr1*Hr1+Hz1*Hz1)
    N *= 2
    Hr,Hz = HrHz(r,z,z1,a,M,N)
    H = np.sqrt(Hr*Hr+Hz*Hz)
    while abs((H-H1)/H)>tol:
        H1 = H
        N *= 2
        Hr,Hz = HrHz(r,z,z1,a,M,N)
        H = np.sqrt(Hr*Hr+Hz*Hz)
    return Hr,Hz,N
                    

P = 100
z = np.linspace(1e-3,1,P)
Hz = np.zeros(P)
list_N = np.zeros(P)
r = 0
z1 = 0
a = 1
M = 1
Nmin = 10
tol=1e-2
for k in range(P):
    hr,hz,N = HrHz_adapt(r,z[k],z1,a,M,Nmin,tol)
    Hz[k] = hz
    list_N[k] = N
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_exact = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz,'b')
plot(z,Hz_exact,'r')
grid()
xlabel('z')
ylabel('Hz')
                   

Hz-2.pdf

Voici les valeurs de N :

figure()
plot(z,list_N)
xlabel('z')
ylabel('N')
yscale('log')
grid()
                   

N-2.pdf

L'algorithme de calcul de l'intégrale peut être amélioré en remarquant que la somme des N-2 termes peut être réutilisée après avoir multiplié N par 2. Pour cela, il faut définir N comme étant le nombre de sous-intervalles de [0,a] (le nombre de points est N+1). On pose donc :

$$\begin{array}{ccc}\hfill & r_k^N=k\frac{a}{N}k=0,\dots N\hfill & \hfill \text{(13)}\\ \hfill & S^N=\sum _{k=1}^{N-1}f\left(r_k^N\right)\hfill & \hfill \text{(14)}\\ \hfill & I^N=\frac{a}{N}(\frac{1}{2}f\left(0\right)+\frac{1}{2}f\left(a\right)+S^N)\hfill & \hfill \text{(15)}\end{array}$$

I^N est l'approximation de l'intégralle pour N sous-intervalles. La somme pour 2N sous-intervalles s'écrit :

$$S^{2N}=S^N+\sum _{p=0}^{N-1}{r}_{2p+1}^{2N}\text{(16)}$$

Pour implémenter cet algorithme, il faut écrire une fonction qui prend en argument la somme S^N/2 et renvoie la somme S^N et l'approximation de l'intégrale I^N :

                     
def HrHz_partiel(r,z,z1,a,M,N,sum_fr,sum_fz):
    dr = a/N
    r1 = np.arange(1,N,2)*dr
    fr,fz = frfz(r,z,r1,z1)
    sum_fr += np.sum(fr)
    sum_fz += np.sum(fz)
    Hz = M/(4*np.pi)*(sum_fz*dr)
    Hr = M/(4*np.pi)*(sum_fr*dr)
    return Hr,Hz,sum_fr,sum_fz
                     

Pour démarrer l'itération à un N quelconque, il faut aussi une fonction qui fait le calcul complet de la somme :

def HrHz_complet(r,z,z1,a,M,N):
    fr0,fz0 = frfz(r,z,0,z1)
    fra,fza = frfz(r,z,a,z1)
    dr = a/N
    r1 = np.arange(1,N)*dr
    fr,fz = frfz(r,z,r1,z1)
    sum_fr = np.sum(fr)+0.5*(fr0+fra)
    sum_fz = np.sum(fz)+0.5*(fz0+fza)
    Hz = M/(4*np.pi)*(sum_fz*dr)
    Hr = M/(4*np.pi)*(sum_fr*dr)
    return Hr,Hz,sum_fr,sum_fz
                     

La fonction suivante effectue les calculs intératifs en partant de N=10.

def HrHz_iter(r,z,z1,a,M,tol):
    N = 10
    Hr1,Hz1,sum_fr,sum_fz = HrHz_complet(r,z,z1,a,M,N)
    H1 = np.sqrt(Hr1*Hr1+Hz1*Hz1)
    N *= 2
    Hr,Hz,sum_fr,sum_fz = HrHz_partiel(r,z,z1,a,M,N,sum_fr,sum_fz)
    H = np.sqrt(Hr*Hr+Hz*Hz)
    while abs((H-H1)/H)>tol:
        H1 = H
        N *= 2
        Hr,Hz,sum_fr,sum_fz = HrHz_partiel(r,z,z1,a,M,N,sum_fr,sum_fz)
        H = np.sqrt(Hr*Hr+Hz*Hz)
    return Hr,Hz,N
    
                     

Voici la reprise du calcul précédent :

P = 100
z = np.linspace(1e-3,1,P)
Hz = np.zeros(P)
list_N = np.zeros(P)
r = 0
z1 = 0
a = 1
M = 1
Nmin = 10
tol=1e-2
for k in range(P):
    hr,hz,N = HrHz_iter(r,z[k],z1,a,M,tol)
    Hz[k] = hz
    list_N[k] = N
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_exact = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz,'b')
plot(z,Hz_exact,'r')
grid()
xlabel('z')
ylabel('Hz')
                   

Hz-3.pdf

et les valeurs de N :

figure()
plot(z,list_N)
xlabel('z')
ylabel('N')
yscale('log')
grid()
                   

N-3.pdf

3. Lignes de champ d'un aimant

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

def champAimant(r,z,a,zs,zn,M,mu0=1,tol=1e-2):
    Hr1,Hz1,N1 = HrHz_iter(r,z,zs,a,-M,tol)
    Hr2,Hz2,N2 = HrHz_iter(r,z,zn,a,M,tol)
    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 μ₀ 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-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-ligneB.pdf