Bobine épaisse

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.

2. Mise en équation

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

Figure pleine page

La bobine est un tore de section rectangulaire, de longueur $\ell$, 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 :

$${\nabla }^2\left(A_{\theta }\overrightarrow{e_{\theta }}\right)=-{\mu }_0{J}_{\theta }\overrightarrow{e_{\theta }}\text{(1)}$$

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

$$\begin{array}{ccc}\hfill & B_r=-\frac{\partial A_{\theta }}{\partial z}\hfill & \hfill \text{(2)}\\ \hfill & B_z=\frac{1}{r}\frac{\partial \left(r{A}_{\theta }\right)}{\partial r}\hfill & \hfill \text{(3)}\end{array}$$

La densité de courant $J_{\theta }$ est uniforme. L'intensité du courant à travers une section du tore est :

$$I=J_{\theta }e\ell \text{(4)}$$

La résolution numérique sera faite avec un terme de source égal à -1 et une largeur de domaine égale à 1. Soit L_d la largeur rélle du domaine. Le champ $A_{\theta }$ (sans dimensions) obtenu par ce calcul devra être multiplié par :

$$A_0={\mu }_0{J}_{\theta }{L}_d^2={\mu }_0\frac{I}{e\ell }{L}_d^2\text{(5)}$$

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

$$B_0=\frac{A_0}{L_d}={\mu }_0\frac{I}{e\ell }{L}_d\text{(6)}$$

afin d'obtenir des valeurs en tesla.

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

$$E_m=\frac{1}{2}{\iiint }_V\overrightarrow{J}\cdot \overrightarrow{A}dv=\frac{1}{2}{\iiint }_V{J}_{\theta }{A}_{\theta }dv\text{(7)}$$

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

$$E_m=\frac{1}{2}L{I}_f^2\text{(8)}$$

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

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

$$L={\iiint }_V{J}_{\theta }{A}_{\theta }dv\text{(9)}$$

avec $J_{\theta }=\frac{N}{e\ell }$ (densité de courant si I_f=1 A). 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=\frac{2\pi N}{e\ell }\iint A_{\theta }(z,r)rdrdz\text{(10)}$$

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

La discrétisation du domaine de calcul est définie par $z_i=i\Delta z$ et $r_j=j\Delta 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 :

$$L\approx \frac{2\pi N}{e\ell }{\mu }_0\frac{N}{e\ell }{L}_d^2\sum _{i,j}{A}_{i,j}j{\left(\Delta r\right)}^2\Delta z\text{(11)}$$

où A_i,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\pi {\mu }_0\frac{L_d^2}{{\left(e\ell \right)}^2}{\left(\Delta r\right)}^2\Delta z{N}^2\sum _{i,j}{A}_{i,j}j\text{(12)}$$

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 N², 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=\frac{\ell }{d}\frac{e}{d}\text{(13)}$$

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 dans l'espace qui sépare les fils, ce qui doit avoir 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 2ⁿ⁺¹ mailles dans la direction de l'axe de symétrie z et 2ⁿ mailles dans la direction radiale. Pour définir les objets, on utilise un maillage réduit de 2^p+1 par 2^p. 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 la sa largeur en mètres.

n=9
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é sur la maille réduite 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=\ell =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() 
            

plotA.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()
            

plotB.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()
            

plotC.pdf

print(numpy.max(Bz[2**n,:]))
--> 0.008157453

Pour donner une valeur à $B_0$, considérons un courant d'intensité I=1 A. La largeur du domaine complet est L_d=1 m. On a alors $B_0={\mu }_0{J}_{\theta }{L}_d=1,29\cdot 10^{-3}T$.

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) :

$$B_z\left(0\right)=0,00819\times 1,29\cdot 10^{-3}=1,06\cdot 10^{-5}T\text{(14)}$$

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

$$B_z\left(0\right)=\frac{{\mu }_{0}I}{2a}=1,01\cdot 10^{-5}T\text{(15)}$$

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()
			

plotC2.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()
            

plotD.pdf

Pour le calcul de l'auto-inductance (expression $\mathrm{(11)}$), 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))
--> 1024

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.3219616387711536e-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=1 mm :

d = 1e-3
N = (l/d)*(e/d)
			 

print(L*N**2)
--> 0.12607208621703658

On sauvegarde les résultats dans des fichiers textes pour une utilisation ultérieure.

numpy.savetxt("bobine-1-Atheta.txt",bobine.get_array(symetry=-1))
numpy.savetxt("bobine-1-AthetaDr.txt",bobine.get_derivR(symetry=1))
numpy.savetxt("bobine-1-AthetaDz.txt",bobine.get_derivZ(symetry=-1))
numpy.savetxt("bobine-1-Br.txt",Br)
numpy.savetxt("bobine-1-Bz.txt",Bz)
numpy.savetxt("bobine-1-r.txt",bobine.get_r(symetry=1))
numpy.savetxt("bobine-1-z.txt",bobine.get_z())
bobine.close()
            

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 $\ell =5/64m$.

n=9
p=6
levels = n-p+1
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()  
            

plotE.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()
            

plotF.pdf

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()
            

plotG.pdf

print(numpy.max(Bz[2**n,:]))
--> 0.02499651

On a comme précédemment :

$$B_0={\mu }_0{J}_{\theta }{L}_d\text{(16)}$$

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

$$B_{\infty }={\mu }_0{J}_{\theta }e\text{(17)}$$

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

$$\frac{0,025B_0}{B_{\infty }}=0,025\frac{L_d}{e}=0,80\text{(18)}$$

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

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

plotH.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.957180540454258e-08

d = 1.0e-3
N = (l/d)*(e/d)
			 

print(N)
--> 4882.8125

print(L*N**2)
--> 1.4203025199065824

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

$$L_1={\mu }_0\frac{N^2}{\ell }\pi a^2\text{(19)}$$

L1 = mu0*N**2/l*numpy.pi*a**2
			

print(L1)
--> 2.3530970576022527

Sans surprise, l'approximation du sonénoïde infini donne une valeur d'auto-inductance très surestimée.

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

$$L_2=\frac{2,2\cdot 10^{-6}{\left(2a\right)}^2{N}^2}{2a+2,2\ell }\text{(20)}$$

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). Rappelons que le modèle de bobine consiste à considérer que le courant est parfaitement uniforme sur une section transversale du tore, alors qu'en réalité le courant est concentré dans les fils et nul dans l'espace qui sépare les fils.

numpy.savetxt("bobine-2-Atheta.txt",bobine.get_array(symetry=-1))
numpy.savetxt("bobine-2-AthetaDr.txt",bobine.get_derivR(symetry=1))
numpy.savetxt("bobine-2-AthetaDz.txt",bobine.get_derivZ(symetry=-1))
numpy.savetxt("bobine-2-Br.txt",Br)
numpy.savetxt("bobine-2-Bz.txt",Bz)
numpy.savetxt("bobine-2-r.txt",bobine.get_r(symetry=1))
numpy.savetxt("bobine-2-z.txt",bobine.get_z())
            

bobine.close()