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.
Figure pleine pageLe 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.
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 :
La densité de courant dans l'induit (cylindre de rayon a) est de la forme :
Dans celui-ci, la loi d'Ohm s'applique :
Les plans perpendiculaires à l'axe étant des plans de symétrie du courant, on a pour le champ magnétique :
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 :
Sur la surface de l'induit (en r=a), il y a continuité de la composante tangentielle de l'excitation magnétique . On a donc la condition limite sur la surface de l'induit :
Les équations de Maxwell pour le champ électromagnétique dans l'induit sont :
La dernière s'écrit aussi , ce qui implique que la divergence de est nulle et montre que la densité de charge dans l'induit est nulle.
On obtient avec ces équations :
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 :
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 :
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 :
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 :
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 :
Finalement, les déterminations de et de 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 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 est :
étant convenu que les grandeurs complexes sont proportionnelles à .
Les conditions limites sont :
où J0 est l'amplitude de la densité de courant dans l'inducteur.
Pour le champ magnétique, l'équation différentielle est :
On s'intéressera aussi à la puissance moyenne dissipée dans l'induit (par unité de longueur) :
L'intensité du courant dans l'induit (par unité de longueur) est :
L'intensité du courant dans l'inducteur (par unité de longueur) est :
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 :
Remarque : il n'est pas attendu que L et M soient indépendants de la fréquence.
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 :
Figure pleine pageL'expression du champ magnétique dans l'entrefer est inchangée.
Le champ électrique est dans la direction de et obéit à l'équation :
avec les conditions limites :
La résolution de cette équation se fait sans difficulté et conduit à :
avec :
On peut aussi écrire où 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 . Dans ce cas, on a :
La puissance moyenne dissipée par unité de surface de l'induit est :
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 :
Lorsque on a :
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
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.
Développons le laplacien dans l'équation :
où la grandeur complexe est simplement notée Eθ.
L'introduction de la variable réduite r'=r/a conduit à :
Les conditions limites sont :
et, d'après :
ou bien avec la variable réduite :
où J0 désigne l'amplitude de la densité de courant dans l'inducteur.
On pose :
où 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 :
L'écriture des dérivées sous forme de différences finies conduit à (pour ) :
qui s'écrit aussi
ou encore :
avec :
et les conditions limites s'écrivent :
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 :
Les trois triagonales sont :
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')figA.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, . 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')figB.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 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)figC.pdf
L'effet thermique de cette puissance doit être évalué au moyen de l'équation de diffusion thermique :
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')figD.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')figE.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 , 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.
On s'intéresse au cas où l'induit est un tube :
Figure pleine pageDans l'induit (b<r<a), on a :
Dans le creux du tube (r<b), on a :
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 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 :
où C1 est une constante. Dans l'induit, on a :
La fonction est impaire donc C2=-C3 :
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 :
La continuité du champ électrique permet d'obtenir l'autre constante :
Dans le creux, l'équation de Maxwell-Faraday s'écrit :
Le champ magnétique dans le creux est donc uniforme :
L'amplitude du champ dans le creux est :
où 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.