Table des matières Python

Bobine épaisse : différences finies

1. Introduction

Ce document présente le calcul du champ magnétique créé par une bobine à symétrie axiale, par résolution numérique de l'équation de Poisson pour le potentiel vecteur.

Le calcul du coefficient d'auto-inductance de la bobine est aussi effectué.

La résolution numérique de l'équation de Poisson par la méthode des différences finies est expliquée dans Équation de Poisson à deux dimensions. La discrétisation de l'équation est expliquée dans Discrétisation de l'équation de Poisson vectorielle.

Le module python utilisé est présenté dans Équation de Poisson : programme Python.

Un calcul similaire par la méthode des éléments finis est présenté dans Bobine épaisse : éléments finis.

2. Mise en équation

On considère une bobine à symétrie axiale, modélisée par une densité de courant volumique orthoradiale.

../../../../figures/sciphys/elecmag/bobine/bobine.svgFigure pleine page

La bobine est un tore de section rectangulaire, de longueur , d'épaisseur e et de rayon moyen a.

Ce modèle peut aussi s'appliquer à un induit de forme torique, avec une densité de courant orthoradiale.

Le potentiel vecteur vérifie l'équation de Poisson en coordonnées cylindriques :

2(Aθeθ)=-μ0Jθeθ(1) r(1r(rAθ)r)+2Aθz2=-μ0Jθ(2)

Le champ magnétique (rotationnel du potentiel vecteur) est :

Br=-Aθz(3)Bz=1r(rAθ)r(4)

La densité de courant Jθ est uniforme. L'intensité du courant à travers une section du tore est :

I=Jθe(5)

La résolution numérique sera faite avec un terme de source égal à -1 et une largeur de domaine égale à 1. Soit Ld la largeur réelle du domaine. Le champ Aθ (sans dimensions) obtenu par ce calcul devra être multiplié par :

A0 = μ0JθLd2=μ0IeLd2(6)

afin d'obtenir un potentiel vecteur en teslas mètres. Les champs (Br,Bz) devrons être multipliés par :

B0=A0Ld=μ0IeLd(7)

afin d'obtenir des valeurs en teslas.

L'énergie magnétique (énergie à fournir pour établir le champ) s'écrit :

Em=12VJAdv=12VJθAθdv(8)

V désigne le volume de la bobine (volume du tore). L'énergie s'exprime aussi à l'aide du coefficient d'auto-inductance :

Em=12LIf2(9)

If est l'intensité du courant dans le fil de la bobine. Si celle-ci est constituée de N spires, on a I=NIf.

Il s'en suit que l'auto-inductance L peut être calculée par :

L=VJθAθdv(10)

avec Jθ=Ne (densité de courant si If=1A). En toute rigueur, il s'agit d'un calcul d'auto-inductance en régime quasi stationnaire de très basse fréquence, où on suppose que la densité de courant est uniforme dans la section transversale des fils. L'auto-inductance s'écrit donc :

L=2πNeAθ(z,r)rdrdz(11)

où l'intégrale est calculée sur une section rectangulaire du tore.

La discrétisation du domaine de calcul est définie par zi=iΔz et rj=jΔr , comme expliqué dans Discrétisation de l'équation de Poisson vectorielle. Une approximation de l'intégrale peut être obtenue par la méthode des rectangles :

L2πNeμ0NeLd2i,jAi,jj(Δr)2Δz(12)

Ai,j sont les valeurs du potentiel vecteur obtenues avec un terme de source égal à -1 et une largeur de domaine égale à 1.

On obtient finalement :

L=2πμ0Ld2(e)2(Δr)2Δz N2i,jAi,jj(13)

On calculera l'auto-inductance pour N=1, valeur qui ne dépend pas du diamètre du fil pour une bobine de taille donnée, puis on multipliera par N2, qui dépend du diamètre du fil. Pour un fil de diamètre d, le nombre de spires est égal au nombre de spires par couche multiplié par le nombre de couches :

N=ded(14)

Ce calcul d'auto-inductance repose sur la modélisation du courant dans une section transversale du tore par un courant uniforme. Il ne tient pas compte du fait que la densité de courant et nulle de l'air entre les fils, qui a un effet à la fois sur la densité de courant et sur le potentiel vecteur dans les fils.

3. Bobine de faible longueur

from matplotlib.pyplot import *
import math
import numpy
import poisson.main
            

Le domaine de calcul (voir figure ci-dessus) comporte 2n+1 mailles dans la direction de l'axe de symétrie z et 2n mailles dans la direction radiale. Pour définir les objets, on utilise un maillage réduit de 2p+1 par 2p. On commence par créer le maillage, discrétiser le laplacien et fixer des conditions limites de Dirichlet sur les bords du domaine (potentiel nul sur les bords). La largeur du domaine est 1 et on considère que c'est aussi sa largeur en mètres.

n=10
p=6
levels = n-p+1
width,height = 1.0,0.5 #taille du domaine en mètres
bobine=poisson.main.PoissonAxialVect(p+1,p,levels,width,height,x0=-0.5,y0=0.0)
bobine.laplacien()
bobine.dirichlet_borders(0.0)
echelle = width/2**(p+1) # echelle de la maille réduite en mètres
            

Pour définir le courant (la source), on se place sur le maillage réduit.

ci,cj = 2**(p+1)//2,2**(p+1)//16 # centre du rectangle
w = h = 2 # demi-largeur et demi-hauteur du rectangle
bobine.source_rect(ci,cj,w,h,-1.0)
e = 2*h*echelle
l = 2*w*echelle
a = cj*echelle
            
print(e)
--> 0.03125
print(a)
--> 0.0625

Compte tenu de la taille du domaine complet (1 par 1), une unité du maillage réduit correspond à une longueur de 1/128. Si on considère que la taille du domaine est donnée en mètres, on a donc pour ces valeurs : e==1/32m et a=1/16m .

La taille du domaine doit être assez grande par rapport au rayon du tore pour que l'effet des bords du domaine soit négligeable dans la région qui nous intéresse (disons la moitié du domaine).

La discrétisation de l'équation de Poisson conduit à un système d'équations linéaires qui est résolu par la méthode itérative de Gauss-Seidel (avec sur-relaxation) ou bien par la méthode multigrilles avec itérations de Gauss-Seidel (beaucoup plus rapide).

#result=bobine.opencl_iterations_norm(100,60,omega=1.95)
result=bobine.multigrid_full_norm(2,levels,100,2,2)
            

Le tracé de la norme en fonction du nombre d'itérations permet de suivre la convergence :

figure(figsize=(8,5))
plot(result[0],result[1]) 
xlabel('niter')  
ylabel('norm')
grid() 
            
plotAplotA.pdf

On récupère le potentiel et les composantes du champ magnétique. L'argument symetry permet d'obtenir le champ sur un domaine comportant aussi les valeurs négatives de r (domaine complet). Lorsqu'il vaut 1, le champ obtenu est pair. Lorsqu'il vaut -1, le champ est impair.

Atheta=bobine.get_array(symetry=-1)
Br=-bobine.get_derivZ(symetry=-1)
Bz=bobine.get_derivRUR(symetry=1)
            

Tracé des lignes d'égale valeur de Aθ :

figure(figsize=(8,8))
extent=bobine.get_extent(symetry=1)
contour(Atheta,20,extent=extent)
imshow(~bobine.get_mask(symetry=1),extent=extent,alpha=0.5, cmap=cm.gray)
xlabel('z')
ylabel('r')
grid()
            
plotBplotB.pdf

Rermarque : ces lignes d'isovaleur de Aθ ne sont pas les lignes du champ magnétique.

Tracé du champ sur l'axe :

  
figure(figsize=(10,5))
z = bobine.get_z()
plot(z,Bz[2**n,:])
xlabel('z')
ylabel(r'$B_z/B_0$')
grid()
            
plotCplotC.pdf
print(numpy.max(Bz[2**n,:]))
--> 0.007927702

Pour donner une valeur à B0 , considérons un courant d'intensité I=1A. La largeur du domaine complet est Ld=1m. On a alors B0=μ0JθLd=1,2910-3T .

mu0 = 4*numpy.pi*1e-7
I = 1.0
J = I/(e*l)
B0 = mu0*J*width
A0 = mu0*J*width**2
			 
print(B0)
--> 0.0012867963509103793

La valeur du champ magnétique au centre de la bobine est donc (pour une intensité de 1 A) :

Bz(0)=0,00819×1,2910-3=1,0610-5T(15)

Comparons cette valeur à celle pour une boucle circulaire (fil de section nulle) :

Bz(0)=μ0I2a= 1,0110-5T(16)

La valeur du champ au centre est donc un peu plus grande pour cette bobine épaisse que pour une spire dont le rayon est égal au rayon moyen.

Voici la comparaison avec l'expression du champ sur l'axe d'une spire circulaire :

figure(figsize=(10,5))
z = bobine.get_z()
Bz_spire = mu0/2*I*a**2/(a**2+z**2)**1.5
plot(z,Bz[2**n,:])
plot(z,Bz_spire/B0,'r--')
xlabel('z')
ylabel(r'$B_z/B_0$')
grid()
			
plotC2plotC2.pdf

Tracé du champ dans le plan z=0 :

figure(figsize=(10,5))
plot(bobine.get_r(symetry=1),Bz[:,2**n])
xlabel('r')
ylabel(r'$B_z/B_0$')
grid()
            
plotDplotD.pdf

Pour le calcul de l'auto-inductance (expression (12)), on a besoin des indices de maille des sommets du rectangle qui définit la section du tore :

f = 2**(n-p)
i1 = (ci-w)*f
i2 = (ci+w)*f
j1 = (cj-h)*f
j2 = (cj+h)*f
			

La précision de ce calcul est déterminée par le nombre de nœuds dans le rectangle :

print((j2-j1)*(i2-i1))
--> 4096

Voici le calcul de l'auto-inductance :

A=bobine.get_array(symetry=0)
def autoInductance(A,i1,i2,j1,j2):
	(nj,ni) = A.shape
	L = 0
	for j in range(j1,j2):
		L += numpy.sum(A[j,i1:i2])*j
	Delta_r = height/(nj-1)
	Delta_z = width/(ni-1)
	L *= 2*numpy.pi/(e*l)**2*mu0*width**2*Delta_r**2*Delta_z
	return L

L = autoInductance(A,i1,i2,j1,j2)
			
print(L)
--> 1.294149588534613e-07

Cette valeur est l'auto-inductance par spire; il faut la multiplier par le nombre de spires au carré pour obtenir l'auto-inductance de la bobine. Si par exemple, le fil de la bobine a un diamètre d=1mm :

d = 1e-3
N = (l/d)*(e/d)
			 
print(L*N**2)
--> 0.12341972241731768

4. Bobine longue

La bobine a toujours un rayon moyen a=1/16m et une épaisseur e=1/32m . Sa longueur est 5 fois son épaisseur, soit =5/32m .

n=10
p=6
levels = n-p+1
width,height = 1.0,0.5 #taille du domaine en mètres
echelle = width/2**(p+1) # echelle de la maille réduite en mètres
bobine=poisson.main.PoissonAxialVect(p+1,p,levels,width,height,x0=-0.5,y0=0.0)
bobine.laplacien()
bobine.dirichlet_borders(0.0)
ci,cj = 2**(p+1)//2,2**(p+1)//16 # centre du rectangle
w,h = 10,2 # demi-largeur et demi-hauteur du rectangle
bobine.source_rect(ci,cj,w,h,-1.0)
e = 2*h*echelle
l = 2*w*echelle
a = cj*echelle
			
print(e)
--> 0.03125
print(l)
--> 0.15625
print(a)
--> 0.0625
#result=bobine.opencl_iterations_norm(100,60,omega=1.95)
result=bobine.multigrid_full_norm(2,levels,100,2,2)
figure(figsize=(8,5))
plot(result[0],result[1]) 
xlabel('niter')  
ylabel('norm')
grid()  
            
plotEplotE.pdf
Atheta=bobine.get_array(symetry=-1)
Br=-bobine.get_derivZ(symetry=-1)
Bz=bobine.get_derivRUR(symetry=1)
figure(figsize=(8,8))
extent=bobine.get_extent(symetry=1)
contour(Atheta,20,extent=extent)
imshow(~bobine.get_mask(symetry=1),extent=extent,alpha=0.5, cmap=cm.gray)
xlabel('z')
ylabel('r')
grid()
            
plotFplotF.pdf

Voici le champ magnétique sur l'axe :

B0 = mu0*1/(e*l)            
figure(figsize=(10,5))
z = bobine.get_z()
plot(z,Bz[2**n,:]*B0)
xlabel('z',fontsize=16)
ylabel(r'$B_z (T)$',fontsize=16)
grid()
            
plotGplotG.pdf
print(numpy.max(Bz[2**n,:]))
--> 0.02466112

On a comme précédemment :

B0=μ0JθLd(17)

Il est intéressant de comparer ce résultat au champ magnétique à l'intérieur de la bobine lorsque , qui est :

B=μ0Jθe(18)

Comparons la valeur obtenue au centre de la bobine à cette valeur :

0,025B0B=0,025Lde=0,80(19)

Le champ au centre est donc moins intense que celui créé par une bobine infinie de même rayon et de même épaisseur, ce qui est bien le résultat attendu.

On peut aussi comparer au champ sur l'axe d'une bobine constituée de spires de rayons a, dont l'expression analytique est connue :

Bz(z)=μ0 I2(z+2(a2+(a+2)2)12-z-2(a2+(a-2)2)12)(20)
mu0 = 4*numpy.pi*1e-7
I = 1.0
Bz_spires = mu0*I/(2*l)*((z+l/2)/(a**2+(z+l/2)**2)**(1/2)-(z-l/2)/(a**2+(z-l/2)**2)**(1/2))
figure(figsize=(10,5))
plot(z,Bz[2**n,:]*B0)
plot(z,Bz_spires,"r--")
xlabel('z',fontsize=16)
ylabel(r'$B_z (T)$',fontsize=16)
grid()  
            
plotG1plotG1.pdf

Les deux courbes sont très proches. Nous avons vérifié que la taille du domaine est assez grande : pour la même taille de bobine, augmenter la taille du domaine de calcul d'un facteur 2 a un effet négligeable sur le résultat. Il en est de même de la finesse du maillage (valeur de n).

figure(figsize=(10,5))
plot(bobine.get_r(symetry=1),Bz[:,2**n])
xlabel('r')
ylabel(r'$B_z/B_0$')
grid()
            
plotHplotH.pdf

Calcul de l'auto-inductance :

f = 2**(n-p)
i1 = (ci-w)*f
i2 = (ci+w)*f
j1 = (cj-h)*f
j2 = (cj+h)*f
A=bobine.get_array(symetry=0)
L=autoInductance(A,i1,i2,j1,j2)
			
print(L)
--> 5.907500933649765e-08
d = 1.0e-3
N = (l/d)*(e/d)
			 
print(N)
--> 4882.8125
print(L*N**2)
--> 1.408457978641931

Comparons cette valeur à celle qui est obtenue avec le modèle de solénoïde infini :

L1=μ0N2πa2(21)
L1 = mu0*N**2/l*numpy.pi*a**2
			
print(L1)
--> 2.3530970576022527

L'approximation du solénoïde infini donne une valeur d'auto-inductance plus grande.

Une formule empirique plus précise ([1]) est :

L2=2,210-6(2a)2N22a+2,2(22)
L2 = 2.2e-6*(2*a)**2*N**2/(2*a+2.2*l)
			
print(L2)
--> 1.7484029134114583

Cette valeur est notablement plus grande que celle que nous obtenons (23% plus grande).

bobine.close()            
            

5. Bobine longue avec noyau

Un noyau cylindrique de fer doux (ou de ferrite doux) est placé dans la bobine. Le matériau de ce noyau est supposé linéaire, de perméabilité magnétique relative μr . Dans le domaine correspondant, le potentiel vecteur vérifie l'équation de Poisson sans sources. La perméabilité magnétique apparaît dans la condition limite entre le domaine ferromagnétique doux et l'extérieur. D'après l'équation rotH=J , il y a continuité de la composante tangentielle de H , c'est-à-dire continuité de la composante tangentielle de B/μr . Pour une surface perpendiculaire à l'axe (Oz) séparant un milieu non magnétique (air ou cuivre) d'un milieu magnétique linéaire, on a :

(Aθz)air=(1μrAθz)fer(23)

Pour une surface cylindrique de rayon r, la condition limite s'écrit :

(1r(rAθ)r)air=(1μr1r(rAθ)r)fer(24)

Une condition limite à l'interface entre deux milieux différents est introduite par la fonction interface_polygon. Voici le calcul avec un noyau de perméabilité relative égal à 1000 :

n=10
p=6
levels = n-p+1
width,height = 1.0,0.5
bobine=poisson.main.PoissonAxialVect(p+1,p,levels,width,height,x0=-0.5,y0=0.0)
bobine.laplacien()
bobine.dirichlet_borders(0.0)
ci,cj = 2**(p+1)//2,2**(p+1)//16 # centre du rectangle
w,h = 10,2 # demi-largeur et demi-hauteur du rectangle
R = cj-h
bobine.interface_polygon(ci-w,0,[0,2*w,0],[R,0,-R],0,1,1e-3)
bobine.source_rect(ci,cj,w,h,-1.0)
e = 2*h*echelle
l = 2*w*echelle
a = cj*echelle
			
#result=bobine.opencl_iterations_norm(100,60,omega=1.95)
result=bobine.multigrid_full_norm(2,levels,100,2,2)
figure(figsize=(8,5))
plot(result[0],result[1]) 
xlabel('niter')  
ylabel('norm')
grid()  
            
plotIplotI.pdf
Atheta=bobine.get_array(symetry=-1)
Br=-bobine.get_derivZ(symetry=-1)
Bz=bobine.get_derivRUR(symetry=1)
figure(figsize=(8,8))
extent=bobine.get_extent(symetry=1)
contour(Atheta,20,extent=extent)
imshow(~bobine.get_mask(symetry=1),extent=extent,alpha=0.5, cmap=cm.gray)
xlabel('z')
ylabel('r')
grid()
            
plotJplotJ.pdf

Calcul de l'auto-inductance :

f = 2**(n-p)
i1 = (ci-w)*f
i2 = (ci+w)*f
j1 = (cj-h)*f
j2 = (cj+h)*f
A=bobine.get_array(symetry=0)
L=autoInductance(A,i1,i2,j1,j2)
			
print(L)
--> 2.3700145897717245e-07
d = 1.0e-3
N = (l/d)*(e/d)
			 
print(N)
--> 4882.8125
print(L*N**2)
--> 5.650555109433471

Voici le champ magnétique sur l'axe :

  
figure(figsize=(10,5))
plot(bobine.get_z(),Bz[2**n,:])
xlabel('z')
ylabel(r'$B_z/B_0$')
grid()
            
plotKplotK.pdf

et le champ magnétique en r=0 :

figure(figsize=(10,5))
plot(bobine.get_r(symetry=1),Bz[:,2**n])
xlabel('r')
ylabel(r'$B_z/B_0$')
grid()
            
plotLplotL.pdf
bobine.close()
Références
[1]  T. Wildi,  Electrotechnique,  (DeBoeck Université, 2000)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.