Table des matières Python

Force entre deux aimants permanents

1. Introduction

Ce document présente le calcul de la force d'interaction entre deux aimants droits coaxiaux. La méthode de calcul de l'énergie magnétique du système de deux aimants est présentée en détail dans Aimant permanent.

L'implémentation présentée est en langage C++.

2. Méthode de calcul

Cette partie donne succinctement la méthode de calcul, détaillée dans Aimant permanent.

Le problème est axisymétrique. Un point P de l'espace est repéré par ses coordonnées cylindriques (r,z).

deuxAimants.svgFigure pleine page

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

σ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.

La méthode de calcul de cette intégrale double est détaillée dans Aimant permanent : utilisation des intégrales elliptiques. L'intégrale sur l'angle θ s'exprime en fonction d'intégrales elliptiques complètes de première et deuxième espèces, lesquelles sont calculées par l'algorithme de Carlson [1]. L'intégration sur r se fait par la méthode des trapèzes avec une augmentation itérative du nombre de points jusqu'à satisfaire une tolérance donnée.

L'énergie magnétique est calculée pour un domaine de l'espace compris entre zmin et zmax et de rayon inférieur à rmax par l'intégrale suivante :

Le maillage du plan (z,r) est défini par :

n est un nombre entier strictement supérieur à 1. Nous avons choisi n=2.

Le domaine d'espace de calcul de l'énergie doit être assez grand pour que l'énergie en dehors de ce domaine soit négligeable, ce qui est facile à vérifier en faisant varier la taille du domaine.

3. Implémentation

3.a. Code source

Code source : aimant.cpp.

Documentation du code source.

3.b. Compilation

Fichier de configuration de CMake : CMakeLists.txt.

Compilation sur Windows avec CMake et MSVC, exécuter dans le répertoire contenant aimant.cpp et CMakeLists.txt :

cmake -G "Visual Studio 14 2015" -B ./
cmake --build ./
                

Compilation sur Linux avec CMake et GCC :

cmake -B ./
cmake --build ./
                

3.c. Exécution

Les paramètres de la simulation et les directives sont définies dans un fichier de configuration config.txt. Pour effectuer la simulation :

aimant.exe config.txt

4. Résultats

4.a. Définition des aimants

Les longueurs sont définies par rapport à une longueur de référence . Les aimantations sont définies par rapport à une aimantation de référence M0. Pour définir deux aimants identiques d'aimantation M0, de longueur , de rayon et distants de , on écrit dans le fichier de configuration config.txt :

La=4
Lb=4
ra=1
rb=1
Ma=1
Mb=1
distance=1
              

4.b. Lignes de champ

Pour obtenir les lignes de champ (H et B) il faut ajouter dans le fichier de configuration :

lignes=lignes # active le tracé de lignes et définit le préfixe des noms des fichiers
lignes_rmax=10 # limite en r du domaine
lignes_zmax=10 # limite en z du domaine 
lignes_dr = 0.2 # espacement radial des lignes qui partent des faces
lignes_ds = 0.01 # pas d'abscisse curviligne entre deux points
lignes_dmin = 0.1 # distance minimale d'approche des faces
lignes_nb_aimants=2 # nombre d'aimants
            

La méthode numérique utilisée pour calculer les points d'une ligne est une méthode d'Euler contrôlée par le pas d'abscisse curviligne. Si lignes_nb_aimants=1, seules les lignes de champ de l'aimant A sont tracées.

Deux fichiers textes sont créés, un pour les lignes de H, l'autre pour les lignes de B. Dans un fichier, chaque ligne de champ débute par #Ligne suivi du numéro de la ligne de champ. Les lignes suivantes comportent les points sous la forme z r. La fonction Python suivante permet de lire un fichier contenant des lignes de champ et d'extraire une liste dont chaque élément est un tableau numpy.ndarray à deux dimensions définissant une ligne de champ.

import numpy as np
from matplotlib.pyplot import *

def lire_lignes(fichier):
    liste_lignes = []
    with open(fichier) as f:
        for line in f:
            line = line.strip()
            if line[0:7]=="#Ligne ":
                line = line.replace("#Ligne ","")
                n = int(line)
                if n>1:
                    liste_lignes.append(np.array(ligne))
                ligne = []
            else:
                line = line.split(" ")
                z = float(line[0])
                r = float(line[1])
                ligne.append([z,r])
    return liste_lignes
               

Voici le tracé des lignes de H pour les deux aimants définis ci-dessus :

liste_lignes = lire_lignes('lignes-H.txt')
figure(figsize=(8,8))
for ligne in liste_lignes:
    plot(ligne[:,0],ligne[:,1],'b-')
    plot(ligne[:,0],-ligne[:,1],'b-')
grid()
axis('equal')
xlabel('z')
ylabel('r')
xlim(-10,10)
ylim(-10,10)
               
lignesH-1lignesH-1.pdf

Voici le tracé des lignes de B:

liste_lignes = lire_lignes('lignes-B.txt')
figure(figsize=(8,8))
for ligne in liste_lignes:
    plot(ligne[:,0],ligne[:,1],'b-')
    plot(ligne[:,0],-ligne[:,1],'b-')
grid()
axis('equal')
xlabel('z')
ylabel('r')
xlim(-10,10)
ylim(-10,10) 
               
lignesB-1lignesB-1.pdf

4.c. Calcul de l'énergie magnétique et de la force

Les directives suivantes permettent d'activer le calcul de l'énergie et définissent le domaine de l'espace pris en compte :

energie=energie # activation du calcul de l'énergie et préfixe du nom de fichier
energie_Rmax = 10
energie_zmin = -20
energie_zmax = 20
            

Le domaine doit être assez large pour contenir la quasi totalité de l'énergie (au centième près).

La définition du maillage du plan (z,r) pour le calcul de l'intégrale doit se faire afin que les faces des aimants soient précisément sur des nœuds de ce maillage. Par ailleurs, le nombre de points dans la direction z doit être divisible par 4 car le calcul de l'énergie est fait par quatre tâches parallèles, qui effectuent chacune le calcul de l'énergie d'un quart du domaine. Dans le cas présent, la taille suivante convient :

energie_Nz = 400
energie_Nr = 100             
              

La boucle de calcul de l'énergie consiste à faire varier la position de l'aimant B et à calculer l'énergie pour chaque distance. Les résultats sont enregistrés dans un fichier texte. Voici par exemple comment définir 20 déplacement de 0.1 :

energie_deplac = 0.1
energie_ndeplac = 20                
                

Les résultats ci-après ont été obtenus avec une distance initiale entre les deux aimants définie par :

distance=0.1

ce qui fait que la distance varie de 0.1 jusqu'à 2 par pas de 0.1.

Il est possible d'attribuer une valeur à la perméabilité magnétique du vide mais nous préférons faire le calcul avec la valeur 1 (qui est la valeur par défaut) :

energie_mu0=1                
                

Voici l'énergie en fonction de la distance :

[z1,W1] = np.loadtxt("aimants-1-Wm.txt",unpack=True,skiprows=1)
figure()
plot(z1,W1)
grid()
xlim(0,2) 
ylim(0,3)
xlabel(r"$d/\ell_0$")
ylabel(r"$W_m/W_0$")
                
energie-1energie-1.pdf

Il s'agit en fait de l'énergie par rapport à une énergie de référence .

La force est obtenue par dérivation de l'énergie par rapport à la distance. On utilise pour cela une différence finie à droite, ce qui implique que le point correspondant à la plus petite distance est perdu (0,1 dans cet exemple).

F1 = []
dz = z1[1]-z1[0]
for i in range(1,len(W1)):
    F1.append(-(W1[i]-W1[i-1])/dz)
F1 = np.array(F1) 
figure()
plot(z1[1:len(z1)],F1,'b')
grid()
xlim(0,2)
ylim(-1.5,0)
xlabel(r"$d/\ell_0$")
ylabel(r"$F_z/F_0$")
                
force-1force-1.pdf

Pour vérifier que le maillage du plan (z,r) est assez-fin, on refait le calcul avec un maillage deux fois plus fin dans la direction z et avec un pas de déplacement de l'aimant deux fois plus petit :

energie_deplac = 0.05
energie_ndeplac = 40
energie_Nz = 800
                
[z1b,W1b] = np.loadtxt("aimants-1b-Wm.txt",unpack=True,skiprows=1)
F1b = []
dz = z1b[1]-z1b[0]
for i in range(1,len(W1b)):
    F1b.append(-(W1b[i]-W1b[i-1])/dz)
F1b = np.array(F1b)
plot(z1b[1:len(z1b)],F1b,'r')
                
force-1bforce-1b.pdf

Pour compléter la courbe aux très courtes distances, on fait une simulation avec les paramètres suivants :

distance=0.01
energie_deplac = 0.01
energie_ndeplac = 10
energie_Nz = 4000
                
[z1c,W1c] = np.loadtxt("aimants-1c-Wm.txt",unpack=True,skiprows=1)
F1c = []
dz = z1c[1]-z1c[0]
for i in range(1,len(W1c)):
    F1c.append(-(W1c[i]-W1c[i-1])/dz)
F1c = np.array(F1c)
plot(z1c[1:len(z1c)],F1c,'r') 
                
force-1cforce-1c.pdf

Les deux courbes (en rouge) ne se raccordent pas tout à fait, ce qui montre qu'une diminution de la taille des mailles est encore profitable (mais très coûteuse en temps de calcul).

Cette force est définie par rapport à une force de référence . Par exemple pour deux aimants en Alnico, matériau dont le champ magnétique rémanent est B0=1,3T, on a M0=B00. Si les deux aimants ont une longueur réelle de 40 mm, on a et :

Dans ce cas présent, la force d'attraction lorsque les deux aimants sont en contact (face nord contre face sud) peut se deviner par extrapolation de la courbe ci-dessus : elle est d'environ 1,5 par rapport à la force de référence, ce qui fait une force de 200N.

En divisant par l'aire des faces de l'aimant, on obtient une force surfacique . Il est intéressant de comparer cette force surfacique à la densité volumique d'énergie magnétique entre les deux aimants. Lorsque les deux aimants sont très proches, les deux faces en vis à vis forment l'équivalent d'un condensateur plan de densité surfacique de charge magnétique σm=M. L'excitation magnétique entre ces faces a donc une intensité H=M. Il s'en suit que B=μ0M et que la densité d'énergie magnétique s'écrit :

La valeur de cette densité d'énergie dans le cas de l'Alnico est .

La force par unité de surface est donc du même ordre de grandeur que la densité d'énergie magnétique entre les deux faces en vis à vis (peut-être sont-elles théoriquement égales ?).

Si ce résultat se généralise, la force par unité de surface ne doit pas dépendre du rayon des aimants. Il est donc intéressant de refaire la simulation avec des aimants de plus petit rayon :

ra=0.5
rb=0.5
                    
[z2,W2] = np.loadtxt("aimants-2-Wm.txt",unpack=True,skiprows=1)
F2 = []
dz = z2[1]-z2[0]
for i in range(1,len(W2)):
    F2.append(-(W2[i]-W2[i-1])/dz)
F2 = np.array(F2) 
[z2b,W2b] = np.loadtxt("aimants-2b-Wm.txt",unpack=True,skiprows=1)
F2b = []
dz = z2b[1]-z2b[0]
for i in range(1,len(W2b)):
    F2b.append(-(W2b[i]-W2b[i-1])/dz)
F2b = np.array(F2b)
figure()
plot(z2[1:len(z2)],F2,'b')
plot(z2b[1:len(z2b)],F2b,'r')
grid()
xlim(0,2)
ylim(-1.5,0)
xlabel(r"$d/\ell_0$")
ylabel(r"$F_z/F_0$")
                    
force-2force-2.pdf

L'aire des faces est 4 fois plus petite que précédemment. La force à distance nulle peut être estimée à -0,4, ce qui correspond approximativement à 1/4 de la force entre les deux aimants de rayon deux fois plus grand. Ces résultats sont donc en accord avec le fait que la force par unité de surface des faces ne dépend pas du rayon de celles-ci (lorsque les aimants sont en contact).

Traçons le rapport des forces pour les deux rayons d'aimants en fonction de la distance :

rapport = F1b/F2b 
figure()
plot(z2b[1:len(z2b)],rapport)
grid()
xlabel(r"$d/\ell_0$")
ylabel('F1/F2')
ylim(4,12)

                    
rapport-forcesrapport-forces.pdf

Le rapport des deux forces est égal à 4 au contact (rapport des aires) puis augmente fortement avec la distance.

Le tracé du logarithme de la norme de la force en fonction du logarithme de la distance D entre les centres des deux aimants permet de savoir si la force suit une loi en puissance en fonction de cette distance. Voici tout d'abord le cas des aimants de rayon ra=rb=1.

figure()
La = Lb = 4
D1b = z1b[1:len(z1b)]+La
plot(np.log(D1b),np.log(abs(F1b)))
grid()
xlabel(r"$\ln(D/\ell_0)$")
ylabel(r"$\ln(|Fz|/F_0)$")

                    
force-1-logforce-1-log.pdf

Voici la pente :

coef= np.polyfit(np.log(D1b),np.log(abs(F1b)),deg=1)
                    
print(coef)
--> array([-6.06799564,  8.66921359])

La pente est d'environ -6, ce qui signifie que |Fz| diminue comme l'inverse de la distance D à la puissance 6. C'est une décroissance plus rapide que la décroissance de la force d'interaction entre deux dipôles (puissance 4).

Voyons ce qu'il en est pour les aimants de rayon ra=rb=0.5 :

D2b = z2b[1:len(z2b)]+La
coef = np.polyfit(np.log(D2b),np.log(abs(F2b)),deg=1)
                    
print(coef)
--> array([-8.62943822, 10.6309995 ])

Plus le rayon des aimants est petit, plus la force décroît rapidement lorsqu'on les éloigne. Ce résultat n'est pas étonnant compte tenu de l'évolution du rapport F1/F2 représenté plus haut.

En conclusion, les résultats de cette simulation sont en accord avec le fait que la force entre les deux aimants en contact est proportionnelle à l'aire des faces. Par ailleurs, plus cette aire est grande, plus la force décroît lentement avec la distance entre les deux centres. La simulation permet de déterminer la loi de décroissance.

Il est intéressant d'augmenter encore la distance entre les deux aimants afin de déterminer à partir de quelle distance on obtient une interaction dipolaire. Lorsque la distance est grande, le domaine de calcul doit être élargi, y compris dans la direction radiale. Voici les paramètres de cette simulation :

La=4
Lb=4
ra=0.5
rb=0.5
Ma=1
Mb=1
distance=0.5
energie=aimants-2c
energie_zmin = -30
energie_zmax = 60
energie_deplac = 0.5
energie_ndeplac = 50
energie_Rmax = 20
energie_Nz = 1800
energie_Nr = 200
                       

La force d'interaction entre deux dipôles identiques de moment magnétique m s'écrit, en fonction de la distance D entre les dipôles :

Le moment magnétique est égal à l'aimantation multipliée par le volume de l'aimant :

 
[z2c,W2c] = np.loadtxt("aimants-2c-Wm.txt",unpack=True,skiprows=1)
F2c = []
dz = z2c[1]-z2c[0]
for i in range(1,len(W2c)):
    F2c.append(-(W2c[i]-W2c[i-1])/dz) 
F2c = np.array(F2c)
La = Lb = 4
ra = rb = 0.5
D = z2c[1:len(z2c)]+La
m = La*np.pi*(ra)**2
Fdipole = 6/(4*np.pi)*m**2/D**4
figure()
plot(D,np.abs(F2c),'b',label='aimants')
plot(D,Fdipole,'r--',label='dipoles')
grid()
xlabel(r"$D/\ell_0$")
ylabel(r"$|Fz|/F_0$")
yscale('log') 
legend(loc='upper right')
                        
force-3force-3.pdf

Comme déjà constaté plus haut, la décroissance de la force est plus rapide que l'interaction dipolaire pour les faibles distances. La force donnée par la relation semble correcte pour une distance , soit une distance supérieure à 20 fois le diamètre des aimants. En dessous de cette distance, la relation conduit à une sous-estimation de la force, qui atteint un facteur 10 à très courte distance.

Références
[1]  W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery,  Numerical recipes, the art of scientific computing,  (Cambridge University Press, 2007)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.