Table des matières Python

Chauffage par induction

1. Introduction

Ce document décrit un modèle de chauffage par induction. Le système étudié comporte une bobine inductrice d'axe (Oz), et une pièce métallique (l'induit) ayant la forme d'un cylindre de révolution de même axe.

bobine.svgFigure pleine page

Le matériau de l'induit est un métal non magnétique (cuivre ou aluminium par exemple) ou un métal ferromagnétique doux dont la relation constitutive est linéaire. On introduit donc une perméabilité magnétique relative μr pour le matériau de l'induit et la perméabilité magnétique μ=μrμ0.

L'objectif est de déterminer le courant électrique dans l'induit en régime sinusoïdal puis la puissance dissipée.

2. Mise en équation

Le système est étudié en coordonnées cylindriques et les effets de bord aux extrémités sont négligés (inducteurs et induits de longueur infinie). Le vecteur densité de courant dans l'inducteur, d'épaisseur e, est de la forme :

Ji=Ji(t)uθ(1)

La densité de courant dans l'induit (cylindre de rayon a) est de la forme :

J=Jθ(r,t)uθ(2)

Dans celui-ci, la loi d'Ohm s'applique :

J=γE(3)

Les plans perpendiculaires à l'axe étant des plans de symétrie du courant, on a pour le champ magnétique :

B=Bz(r,t)uz(4)

Par ailleurs, ce champ est nul à l'extérieur de la bobine. En conséquence, le champ magnétique dans l'espace qui sépare l'induit de l'inducteur (l'entrefer) est :

Bz(entrefer)=μ0eJi(t)(5)

Sur la surface de l'induit (en r=a), il y a continuité de la composante tangentielle de l'excitation magnétique H=Bμ . On a donc la condition limite sur la surface de l'induit :

Bz(r=a-)=μeJi(t)(6)

Les équations de Maxwell pour le champ électromagnétique dans l'induit sont :

divE=0(7) rotE=-Bt(8) divB=0(9) rotB=μJ (10)

La dernière s'écrit aussi rotB=μγE , ce qui implique que la divergence de E est nulle et montre que la densité de charge dans l'induit est nulle.

On obtient avec ces équations :

ΔE=μγEt(11) ΔB=μγBt (12)

On peut chercher à déterminer en premier indifféremment le champ magnétique ou le champ électrique. Le passage de l'un à l'autre se fait avec la relation :

E=1μγrotB=-1μγBzruθ(13)

Cette relation n'est valable que dans l'induit. Dans l'entrefer et dans l'induit, on a d'après l'équation de Maxwell-Faraday :

1r(rEθ)r=-dBzdt(14)

Pour le champ électrique nous devons cependant écrire la condition limite en r=a. Celle-ci est donnée par la relation ci-dessus appliquée en r=a :

1r(rEθ)r(r=a)=-μedJidt(15)

Il s'agit donc d'une condition de type Neumann. Lorsque le champ électrique dans l'induit est déterminé, cette même équation permettra de déterminer le champ électrique dans l'entrefer, sachant que la composante tangentielle de ce champ est continu à la surface de l'induit. Remarquons qu'il est impossible de déterminer le champ électrique dans l'entrefer dès le début car on ne sait rien a priori du flux magnétique dans l'induit.

Nous avons aussi besoin d'une condition limite en r=0. Pour le champ électrique, il s'agit de :

Eθ(0,t)=0(16)

Pour le champ magnétique, si nous considérons que la fonction Bz(r) est définie sur un intervalle [-a,a], cette fonction est paire et dérivable en r=0, ce qui conduit à la condition limite suivante :

Bzr(r=0)=0(17)

Finalement, les déterminations de E et de B sont similaires puisque l'équation aux dérivées partielles est la même et dans les deux cas il y a une condition limite de Dirichlet et une condition limite de Neumann. Nous préférons effectuer d'abord la détermination de E puisque la détermination du courant dans l'induit est l'objectif premier.

En explicitant l'opérateur laplacien, et en se plaçant en régime sinusoïdale de pulsation ω, l'équation différentielle à résoudre pour le champ complexe Eθ̲(r) est :

dd r(1r(d (rEθ̲)d r))=iωμγEθ̲(18)

étant convenu que les grandeurs complexes sont proportionnelles à eiωt .

Les conditions limites sont :

Eθ̲(0)=0(19)1r(rEθ̲)r(a)=-iωμe J0(20)

J0 est l'amplitude de la densité de courant dans l'inducteur.

Pour le champ magnétique, l'équation différentielle est :

1rdd r(rd Bz̲d r)=iωμγBz̲(21)

On s'intéressera aussi à la puissance moyenne dissipée dans l'induit (par unité de longueur) :

P=12γEθ̲Eθ̲*ds=γ20REθ̲Eθ̲*2πrdr(22)

L'intensité du courant dans l'induit (par unité de longueur) est :

I2̲=0aγEθ̲dr(23)

L'intensité du courant dans l'inducteur (par unité de longueur) est :

I1=J0e(24)

Si L désigne le coefficient d'auto-inductance de l'induit, M le coefficient d'inductance mutuelle et R la résistance de l'induit, on a la relation :

I2̲=iωMI1R+iLω(25)

Remarque : il n'est pas attendu que L et M soient indépendants de la fréquence.

3. Modèle bidimensionnel rectiligne

Lorsque la profondeur de pénétration du champ dans le métal est faible devant le rayon de l'induit, une géométrie bidimensionnel rectiligne devrait donner une bonne approximation des champs. Celle-ci est constituée de deux plaques de courant de taille infinie et d'épaisseur e, disposées parallèlement. L'induit est une plaque d'épaisseur 2a. La figure suivante définit le système :

plaques.svgFigure pleine page

L'expression du champ magnétique dans l'entrefer (5) est inchangée.

Le champ électrique est dans la direction de uy et obéit à l'équation :

d2Ey̲dx2=iωμγEy̲(26)

avec les conditions limites :

Ey̲(0)=0(27) dEy̲dx(a)=-iωμeJ0 (28)

La résolution de cette équation se fait sans difficulté et conduit à :

Ey̲(x)=-1+i2δωμeJ0sinh(1+iδx)cosh(1+iδa)(29)

avec :

δ=2ωγμ(30)

On peut aussi écrire μeJ0=μrB0 B0 est l'intensité du champ magnétique à l'intérieur de l'inducteur en l'absence d'induit (égal au champ dans l'entrefer).

Cette solution devrait s'appliquer à la géométrie cyclindrique (en remplaçant x par r) lorsque δa . Dans ce cas, on a :

Ey̲(x)-1+i2δωμeJ0 exp(-1+iδ(a-x))(31)

La puissance moyenne dissipée par unité de surface de l'induit est :

P=20a12γ|Ey̲(x)|2dx=14γδ(δωμeJ0)2=12ω12γ-12μ12e2J02(32)

La puissance dissipée est donc proportionnelle à la racine carrée de la fréquence et à la racine carrée de la perméabilité magnétique.

L'équation de Maxwell-Faraday permet de calculer le champ magnétique dans l'induit :

Bz̲=-1iωdEy̲dx=μeJ0cosh(1+iδx)cosh(1+iδa)(33)

Lorsque δa on a :

Bz̲μrB0exp(-1+iδ(a-x))(34)

La pertinence de ce modèle peut sembler à première vue compromise par l'hypothèse d'une densité de courant uniforme dans l'inducteur. En réalité, la densité de courant dans l'inducteur ne peut être uniforme et présente évidemment un effet de peau similaire à celui qu'on obtient dans l'induit. Cependant, le champ magnétique dans l'entrefer

B0=μ0eJ0(35)

ne dépend que de l'intensité du courant dans l'inducteur et pas de la répartition du courant. En conséquence les résultats obtenus dans l'induit sont toujours valable si le vecteur densité de courant n'est pas uniforme dans l'inducteur. Il faut remarquer que cette conclusion est valable dans cette géométrie particulière et aussi dans la géométrie cylindrique (cylindre infini) mais qu'elle ne peut se généraliser à tout type de géométrie : dans le cas général, le champ magnétique créé par un courant dans un conducteur dépend, au moins au voisinage de celui-ci, de la répartition du courant dans le conducteur.

4. Résolution par différences finies

Développons le laplacien dans l'équation (18) :

1rdEθdr-Eθr2+d2Eθdr2=iωμγEθ(36)

où la grandeur complexe est simplement notée Eθ.

L'introduction de la variable réduite r'=r/a conduit à :

1r'dEθdr'-Eθr'2+d2Eθdr'2=iωμγa2Eθ(37)

Les conditions limites sont :

Eθ(0)=0(38)

et, d'après (15) :

Eθ(a)a+dEθdr(a)=-iωμeJ0(39)

ou bien avec la variable réduite :

Eθ(1)+dEθdr'(1)=-iωμeaJ0(40)

J0 désigne l'amplitude de la densité de courant dans l'inducteur.

On pose :

A=ωμγa2(41)C=ωaμrμ0eJ0=ωaμrB0(42)

B0 est l'amplitude du champ magnétique dans l'entrefer, qui est aussi celle du champ dans la bobine en l'absence de l'induit.

La discrétisation de l'intervalle [0,1] est effectuée de la manière suivante :

r'n=nh,n=0,1N-1(43)h=1N-1(44)

L'écriture des dérivées sous forme de différences finies conduit à (pour n=1,N-2 ) :

1nhEn+1-En-12h-En(nh)2+En+1+En-1-2Enh2=iA En(45)

qui s'écrit aussi

(1h2-12nh2)En-1-(2h2+1(nh)2+iA)En+(1h2+12nh2)En+1=0(46)

ou encore :

anEn-1+bnEn+cnEn+1=0(47)

avec :

an=1-12n(48)bn=-(2+1n2+iAh2)(49)cn=1+1n2(50)

et les conditions limites s'écrivent :

E0 = 0(51)EN-1+1h(EN-1-EN-2)=-iC(52)

Nous obtenons N-2 équations et 2 conditions limites, soit N équations linéaires pour N inconnues. Il s'agit d'un système tridiagonal :

(10a1b1c1a2b2c2aN-2bN-2cN-2-1h+1)(E0E1E2 EN-2EN-1)=(000 0-iCh)(53)

Les trois triagonales sont :

0,0,c1,c2,cN-2(54)1,b1,b2,bN-2,h+1(55)a1,a2,aN-2,-1,0(56)

Afin que les trois diagonales aient N éléments, on a ajouté un zéro au début de la première et un zéro à la fin de la dernière.

Un système tridiagonal peut être résolu efficacement par un algorithme d'élimination spécifique en O(N) (voir Équation de diffusion à une dimension). La fonction scipy.linalg.solve_banded permet de résoudre efficacement un système avec une matrice bande, donc en particulier une matrice tridiagonale. La matrice doit être définie avec les trois diagonales comme écrites ci-dessus. On doit aussi préciser le nombre maximal de termes non nul au dessus de la diagonale principale (ici 1) et le nombre maximal en dessous (ici 1).

from scipy.linalg import solve_banded
import numpy as np
from matplotlib.pyplot import *
	
	
def champE(A,C,N,cyl=1):
	#cyl = 1 : géométrie cylindriques
	#cyl = 0 : géométrie rectiligne
	P = np.zeros(N,dtype=np.complex64)
	h = 1/(N-1)
	h2 = h**2
	P[N-1] = -1j*C*h
	D1 = np.zeros(N,dtype=np.complex64)
	D2 = np.zeros(N,dtype=np.complex64)
	D3 = np.zeros(N,dtype=np.complex64)
	for n in range(1,N-1):
		D1[n+1] = 1+cyl/(2*n)
	D2[0] = 1
	D2[N-1] = cyl*h+1
	for n in range(1,N-1):
		D2[n] = -2-1j*A*h2-cyl/(n**2)
	for n in range(1,N-1):
		D3[n-1] = 1-cyl/(2*n)
	D3[N-2] = -1
	M = np.array([D1,D2,D3])
	E = solve_banded((1,1),M,P)
	r = np.arange(N)*h
	dr = r[1]-r[0]
	p = 0.5*E*np.conjugate(E)
	P = np.sum(p*r)*dr*2*np.pi
	I2 = np.sum(E)*dr
	return r,E,np.real(p),np.real(P),I2
			

La fonction champE permet aussi de résoudre par différences finies le problème bidimensionnel rectiligne (avec cyl=0), ce qui permet de valider la méthode par comparaison avec la solution formelle.

Voici un exemple : un tube en fer de rayon 1 cm, à la fréquence 100 Hz. Les parties réelle et imaginaire du champ électrique sont tracées pour la géométrie cylindrique résolue par différence finie et les mêmes grandeurs pour le modèle rectiligne.

a = 0.01
gamma = 1e7
mu0 = 4*np.pi*1e-7
mur=1000
mu = mur*mu0
B0 = 0.01
I1 = B0/mu0
f = 100
w = 2*np.pi*f 
A = w*mu*gamma*a**2
C = w*a*B0*mur
delta = np.sqrt(2/(mu*gamma*w))
N = 1000
r,E,p,P,I2 = champE(A,C,N)
P *= gamma
p*=gamma
E_rect = -(1+1j)/2*delta*w*mur*B0*np.sinh((1+1j)/delta*r*a)/np.cosh((1+1j)/delta*a)
figure(figsize=(12,6))
plot(r*a*1000,np.real(E),'r',label='Re')
plot(r*a*1000,np.imag(E),'b-',label='Im')
plot(r*a*1000,np.real(E_rect),'r--',label='Re rect')
plot(r*a*1000,np.imag(E_rect),'b--',label='Im rect')
xlabel("r (mm)",fontsize=16)
ylabel(r"$E_{\theta} (V/m)$",fontsize=16)
grid()
legend(loc='lower left')
		
figAfigA.pdf

Les courbes obtenues sont quasi identiques à celles de la géométrie rectiligne.

La profondeur de pénétration (en mm) est :

print(delta*1e3)
--> 0.5032921210448703

Puissance dissipée par unité de longueur (en kW) :

print(P*1e-3)
--> 4003.0711510918636
print(I2/I1)
--> (-1.0189114214780407e-05-2.023209178338232e-07j)

Le rapport I2/I1 est quasi réel, ce qui montre que, à cette fréquence, RLω . La valeur négative de la partie réelle indique que le courant dans l'induit est de sens opposé au courant dans l'inducteur, ce qui est bien attendu d'après la loi de Lenz.

Voici le résultat pour une fréquence de 10 kHz :

f = 10e3
w = 2*np.pi*f 
A = w*mu*gamma*a**2
C = w*a*B0*mur
delta = np.sqrt(2/(mu*gamma*w))
N=2000
r,E,p,P,I2 = champE(A,C,N)
P *= gamma
p*=gamma
I2 *=gamma
E_rect = -(1+1j)/2*delta*w*mur*B0*np.sinh((1+1j)/delta*r*a)/np.cosh((1+1j)/delta*a)
figure(figsize=(12,6))
plot(r*a*1000,np.real(E),'r',label='Re')
plot(r*a*1000,np.imag(E),'b-',label='Im')
plot(r*a*1000,np.real(E_rect),'r--',label='Re rect')
plot(r*a*1000,np.imag(E_rect),'b--',label='Im rect')
xlim(9.0,10)
xlabel("r (mm)",fontsize=16)
ylabel(r"$E_{\theta} (V/m)$",fontsize=16)
grid()
legend(loc='lower left')
		
figBfigB.pdf
print(delta*1e3)
--> 0.05032921210448704
print(P*1e-3)
--> 47953.4184474598
print(I2/I1)
--> (-109.8601127176706-10.94645337769677j)

Le champ pénètre 10 fois moins dans le métal mais son amplitude est beaucoup plus grande. La puissance dissipée est environ 15 fois plus grande.

La résistance de l'induit a beaucoup augmenté puisque la partie imaginaire de I2 n'est plus petite devant la partie réelle alors que Lω a été multiplié par 100. L'augmentation de la résistance avec la fréquence est une conséquence de la diminution de la profondeur de pénétration.

Voici la densité volumique de puissance en fonction de r au voisinage de la surface :

figure(figsize=(12,6))
plot(r*a*1000,p)
grid()
xlabel("r (mm)",fontsize=16)
ylabel(r"$p\ (\rm W\cdot m^{-3})$",fontsize=16)
xlim((a-3*delta)*1e3,a*1e3)
		
figCfigC.pdf

L'effet thermique de cette puissance doit être évalué au moyen de l'équation de diffusion thermique :

ρcpTt=λ1rr(rTr)+p(57)

Dans les deux exemples précédents, la profondeur de pénétration est beaucoup plus petite que le rayon de l'induit, ce qui fait que la solution bidimensionnelle rectiligne est valable (puisque le champ est non nul seulement au voisinage de la surface). Pour un matériau non magnétique à basse fréquence, voyons une situation où la profondeur de pénétration est du même ordre de grandeur que a :

a = 0.01
gamma = 1e7
mu0 = 4*np.pi*1e-7
mur=1
mu = mur*mu0
B0 = 0.01
I1 = B0/mu0
f = 100
w = 2*np.pi*f 
A = w*mu*gamma*a**2
C = w*a*B0*mur
delta = np.sqrt(2/(mu*gamma*w))
N = 1000
r,E,p,P,I2 = champE(A,C,N)
P *= gamma
p*=gamma
E_rect = -(1+1j)/2*delta*w*mur*B0*np.sinh((1+1j)/delta*r*a)/np.cosh((1+1j)/delta*a)
figure(figsize=(12,6))
plot(r*a*1000,np.real(E),'r',label='Re')
plot(r*a*1000,np.imag(E),'b-',label='Im')
plot(r*a*1000,np.real(E_rect),'r--',label='Re rect')
plot(r*a*1000,np.imag(E_rect),'b--',label='Im rect')
xlabel("r (mm)",fontsize=16)
ylabel(r"$E_{\theta} (V/m)$",fontsize=16)
grid()
legend(loc='lower left')
		
figDfigD.pdf

Dans ce cas, la solution en géométrie cylindrique est très différente de la solution en géométrie rectiligne.

print(delta*1e3)
--> 15.915494309189533
print(P*1e-3)
--> 7.634946679903085

La puissance dissipée est 500 fois moins grande que celle pour le matériau ferromagnétique (de perméabilité relative 1000) avec la même conductivité et la même fréquence (voir plus haut) et sa valeur est insuffisante pour obtenir un effet thermique notable. Pour chauffer un matériau non magnétique, il faut absolument utiliser une grande fréquence :

a = 0.01
gamma = 1e7
mu0 = 4*np.pi*1e-7
mur=1
mu = mur*mu0
B0 = 0.01
I1 = B0/mu0
f = 100e3
w = 2*np.pi*f 
A = w*mu*gamma*a**2
C = w*a*B0*mur
delta = np.sqrt(2/(mu*gamma*w))
N = 1000
r,E,p,P,I2 = champE(A,C,N)
P *= gamma
p*=gamma
E_rect = -(1+1j)/2*delta*w*mur*B0*np.sinh((1+1j)/delta*r*a)/np.cosh((1+1j)/delta*a)
figure(figsize=(12,6))
plot(r*a*1000,np.real(E),'r',label='Re')
plot(r*a*1000,np.imag(E),'b-',label='Im')
plot(r*a*1000,np.real(E_rect),'r--',label='Re rect')
plot(r*a*1000,np.imag(E_rect),'b--',label='Im rect')
xlabel("r (mm)",fontsize=16)
ylabel(r"$E_{\theta} (V/m)$",fontsize=16)
grid()
legend(loc='lower left')
		
figEfigE.pdf
print(delta*1e3)
--> 0.5032921210448703
print(P*1e-3)
--> 4003.0711510918636

À cette fréquence de 100 kHz on obtient une puissance à peine plus grande que pour le matériau ferromagnétique à 100 Hz. D'après le modèle bidimensionnel rectiligne, lorsque δa , La puissance dissipée devrait être proportionnelle à la racine carrée de la perméabilité relative et à la racine carrée de la fréquence. On retrouve bien cette dépendance (à peu près) : il faut multiplier la fréquence par 1000 afin de compenser une diminutation de la permabilité magnétique d'un facteur 1000.

5. Induit creux

On s'intéresse au cas où l'induit est un tube :

bobine-tube.svgFigure pleine page

Dans l'induit (b<r<a), on a :

ΔE=μγEt(58) ΔB=μγBt (59)

Dans le creux du tube (r<b), on a :

ΔE=0(60) ΔB=0 (61)

Il est impossible d'obtenir a priori une condition limite pour le champ électrique en r=b. En effet cette condition serait de la forme (15) et nécessiterait donc la connaissance du champ magnétique en r=b; or le champ magnétique dans le creux ne peut être calculé a priori. Il faut donc résoudre l'équation aux dérivées partielles sur l'intervalle [0,a], sachant que le champ électrique est partout continu.

La résolution du problème bidimensionnel rectiligne se fait sans difficulté. Dans le creux on a :

Ey̲(x)=C1x(62)

C1 est une constante. Dans l'induit, on a :

Ey̲(x)=C2exp(1+iδx)+C3exp(-1+iδx)(63)

La fonction Ey̲ est impaire donc C2=-C3 :

Ey̲(x)=C2(exp(1+iδx)-exp(-1+iδx))(64)

En conséquence, il ne reste qu'une constante inconnue pour le champ dans le métal et celle-ci est déterminée, comme précédemment, avec la condition limite en r=a. Il s'en suit que la solution établie précédemment est encore valable lorsque l'induit est creux. On a donc pour b<r<a :

Ey̲(x)=-1+i2δωμbJ0sinh(1+iδx)cosh(1+iδa)(65)

La continuité du champ électrique permet d'obtenir l'autre constante :

C1=-1+i2δωμeJ0bsinh(1+iδb)cosh(1+iδa)(66)

Dans le creux, l'équation de Maxwell-Faraday s'écrit :

dEy̲dx=-iωBz̲(67)

Le champ magnétique dans le creux est donc uniforme :

Bz̲=iC1ω=1-i2δμeJ0bsinh(1+iδb)cosh(1+iδa)(68)

L'amplitude du champ dans le creux est :

B1=δμr B0b2e-a-bδ(69)

B0 est l'amplitude en l'absence d'induit. Ce résultat est intéressant car il peut faire l'objet d'une vérification expérimentale. Il est en effet aisé de mesurer B0 et B1. Si la profondeur de pénétration est petite devant l'épaisseur du tube, le champ magnétique dans l'induit est négligeable (phénomène d'écrantage). L'étude expérimentale de B1 en fonction de la fréquence devrait être un bon moyen de valider ce modèle.

Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.