Introducere 5 Capitolul 1. Formularea problemei de câmp electromagnetic cvasi-stationar 7 Ecuațiile câmpului electromagnetic în domeniul de calcul 7… [311667]

CURPINS

Introducere 5

Capitolul 1. [anonimizat] 7

Ecuațiile câmpului electromagnetic în domeniul de calcul 7

Condiții de trecere (CT) 8

Condiții de frontieră (FR) 8

Condițiile inițiale (CI) 9

Teoremă de unicitate 9

Capitolul 2. Soluții numerice de analiză a câmpului electromagnetic cuasistationar în probleme de curenți turbionari 14

Modelul (A-V, A) 14

Modelul (A – V,) 19

Modelul (T –

Capitolul 3. Soluționarea numerică a problemei de câmp termic 24

Aspecte calitative. Difuzia câmpului termic 24

Soluționarea problemei de câmp termic prin Metoda

elementului finit 26

Capitolul 4. Cuplarea problemelor de câmp electromagnetic și termic 32

Capitolul 5. Simularea numerică a procesului de încălzire prin inducție electromagnetică 34

Definirea problemei.

Rezultate numerice 36

Concluzii 44

Bibliografie 46

[anonimizat], cu alte cuvinte determinarea variației în timp a câmpului electromagnetic indus în materialul conductor solid. In literatura de specialitate este prezentata importanța utilizării acestui tip de incălzire și posibilitatea de integrare a sitemelor electrotermice inductive în procesele industriale [5], [6], [7], [10].

O problemă de mare interes tehnic o constitue determinarea numerică a soluției prin calculul pierderilor prin curenți turbionari [2]. De cele mai multe ori aceste pierderi nu sunt dorite și tehnicienii încearcă să le reducă cât mai mult posibil. [anonimizat]. [anonimizat] [11].

Majoritatea modelelor numerice de soluționare a acestor fenomene implică un pas de timp iterativ în care câmpul electromagnetic este determinat funcție de temperatură și de proprietățile de material [12], rezultând o distribuție a [anonimizat], [anonimizat].

Lavers [13] prezintă pentru prima dată soluționarea unei problemei de încălzire prin inducție în spațiul 2D pentru un conductor feromagnetic de secțiune paralelipipedică. [anonimizat] alții prezintă un procedeu de corectare și rafinare a pasului de timp pentru tranziția la temperatura Curie, dezvoltând ulterior un program de analiză numerică pentru soluționarea problemelor cuplate [14], incluzând și fenomenele electrotermice. Pachetul softwarea realizat de acestia au pus bazele principiilor care sunt utilizate la proiectarea și realizarea programelor soft-ware utilizate astăzi.

Literatura de specialitate prezintă numeroase lucrări de specialita având ca subiect modelarea numerică a problemelor cuplate de electrotermie pentru medii liniare și neliniare, [2], [3], [9], [15]. O inovație în domeniul modelării numerice a fost utilizarea modelului Preisach pentru aproximarea efectului de histerezis în materiale, [1]. Deasemenea publicatiile recente privind soluționarea numerică a problemelor electromagnetice și termice cuplate, prezintă lucrări orientate pe domenii de cercetare a unor aplicații. Aceste aplicații includ proiectarea sistemului de încălzire a regiunilor de capăt pentru o bară din oțel, model ce include impedanța de sarcină a circuitului în simulare, [4].

Majoritatea aplicațiilor de încălzire prin inducție electromagnetică sunt modelate cu ajutorul metodei elementului finit, (FEM – Finite Element Method). Această metodă a fost descrisă la începutul deceniului șase a secolului XX, dar nu pentru rezolvarea problemelor de câmp electromagnetic. Aplicarea în domeniul încălziri prin inducție a metodei elementului finit s-a realizat abia la începutul anilor “70 de către Silvester, Chari și alții. Această implementare a lui Silvester a condus la o creștere considerabilă a cercetării în domeniul încălziri prin inducție

Cercetările care stau la baza modelării numerice utilizate astăzi de către soft-urile comerciale au fost făcute acum mai bine de 25 de ani, unele incluzând posibilități de tratare a problemelor, complexe, de câmp electromagnetic cuplate cu proprietăți de material, a procesărilor termice și deformărilor mecanice, [15], [16].

Lucrarea de față își propune prezentarea principalelor soluții numerice utilizate în modelarea numerică a problemelor complexe de câmp electromagnetic, respectiv modelarea numerică a proceselor electromagnetice și termice ce apar în cadrul proceselor de încălzire prin inductie electromagnetică în echipamentele cu medii imobile.

Capitolul 1 al lucrării abordează formularea problemei de câmp electromagnetic find prezentate ecuațiile regimurilor, condițiile de frontieră de tip electric și magnetic și teoremele de unicitate.

Capitolul 2 al lucrării realizează o prezentare a pricipalelor tehnici de soluționare a problemelor de curenți turbionari folosind soluții de potențial

Capitolul 3 abordează rezolvarea numerică a problemei de câmp termic prin metoda elementului finit.

În Capitolului 4 se prezintă realizarea cuplajului dintre problema de câmp electromagneti și cea de câmp termic.

Capitolul 5 prezintă simularea numerică pentru o problem de ce presupune de încălzire prin curenți turbionari ai unei piese cilindrice.

În final sunt prezentate concuziile și referințele bibliografice.

1. Formularea problemei de CÂMP ELECTROMAGNETIC CUASISTAȚIONAR

1.1. Ecuațiile câmpului electromagnetic în domeniul de calcul

Formularea problemei de câmp electromagnetic asociată domeniului de calcul pentru corpuri inzolate și în mișcare are la bază ecuațiile lui Maxwell [11]. Se consideră domeniul supraconductor s, respectiv T domeniul din exteriorul lui s (fig 1.1). Domeniul T este format din medii conductoare = k , respectiv medii izolante 0: T = 0 . La rândul său, fiecare domeniu k este reprezentat de un corp conductor cu o lege de mișcare bine definită. Modelul de câmp rezentat poate fi aplicat și dacă în domeniul de calcul sunt corpuri izolante în mișcare.

Fig. 1.1. Domeniu de calcul [11]

În T:

(1.1)

H = F(B) (în particular B = (H + M) (1.2)

În k:

(1.3)

(1.4)

(1.5)

În 0:

(1.6)

(1.7)

unde, densitatea de curent J0 este impusă.

1.2. Condiții de trecere (CT)

Soluționarea problemei presupune impunerea condițiilor trecere. Pe frontiera a domeniului conductor vom avea:

(1.8)

(1.9)

(1.10)

1.3. Condiții de frontieră (FR)

Pentru analiza câmpul electromagnetic în domeniul T, se impune rezolvarea ecuațiilor (1.1),…,(1.7). Frontiera T a domeniului T conține și suprafețele ce mărginesc corpurile supraconductoare s , fără ca mediile supraconductoare să facă parte din domeniul de calcul. Frontiera domeniului T poate fi în mișcare deoarece și corpurile supraconductoare pot fi în mișcare, astfel încăt mișcarea frontierei poate fi arbitrară. Se vor considera toate corpurile conductoare aflate în mișcare că se deplasează în interiorul domeniului T . Presupunem că avem următoarele condiții de frontieră [11] pe frontiera T :

() Pe o parte S' a frontierei se dă componenta tangențială a lui H:

Ht=f (1.11)

() Pe restul frontierei S" = T \ S' se dă componenta normală a lui B:

Bn=g (1.12)

() Dacă S' este formată din n suprafețe disjuncte Si, atunci pentru n-1 din aceste suprafețe se dau fluxurile magnetice:

(1.13)

() Dacă T este multiplu conex atunci fie s holurile pentru care s S'. Fie TSk cu k=1,2,…,n tăieturile (minimale) ce transformă T în domeniu simplu conex și fie Sk, k=1,2,…,ns, curbele ce înconjoară holurile (ns este ordinul de conexitate al domeniului T , iar Sk formează "baza" de curbe ce definește acest ordin). Trebuie impusă una din următoarele condiții de frontieră:

(') (1.14')

sau

(") (1.14")

Observații:

1.1. La suprafața corpurilor supraconductoare avem condiția () omogenă:

Bn= 0 (1.15)

1.2. Spirele supraconductoare conduc la condiții de frontieră de tip (). Condiția (') impune curentul prin spira supraconductoare, în timp ce condiția (") impune fluxul magnetic al spirei supraconductoare.

1.3. Condiția () poate fi formulată și sub formă hibridă: pentru o parte din spire se dau curenții, adică condiția ('), iar pentru cealaltă parte se dau fluxurile magnetice, deci condiția (").

1.4. Condițiile inițiale (CI)

Pentru t=0 trebuie cunoscută valoarea inducției magnetice: . Evident se impune divBi = 0. Dacă regimul este periodic, atunci la ecuațiile (1.1),…,(1.7) și condițiile de frontieră (FR) se adaugă condiția de periodicitate. Condiția inițială nu este necesară. În cazul mediilor în mișcare însă, regimul periodic presupune dificultăți suplimentare de analiză. În această etapă ne mulțumim să determinăm regimul permanent (periodic) ca o soluție asimptotică de regim variabil cu condiții inițiale arbitrare.

1.5. Teoremă de unicitate

Ecuațiile (1.1),…,(1.7), împreună cu condițiile de trecere (CT), cu condițiile de frontieră (FR) și condițiile inițiale (CI) definesc unic câmpul magnetic (B,H) [12].

Demonstrație. Vom presupune că două câmpuri electromagnetice distincte verifică condițiile de mai înainte și fie (Bd, Hd, Ed, Jd) câmpul diferență. Acest câmp verifică relațiile (1.1),…,(1.7) și are condiții de frontieră și condiții inițiale nule. Datorită condițiilor (), () și (1.1) fluxul lui B pe ultima suprafață Sn e nul. Notăm:s

unde integrarea în timp se face în același punct din sistemul de coordonate fixat pe corpul în mișcare. Atunci, datorită condițiilor inițiale, legea inducției electromagnetice devine:

(1.16)

Domeniul conductor poate să facă găuri în 0 (Fig1.1). Fie nc numărul acestor găuri. Numărul de găuri în 0 (Nh) poate să fie mai mare decât în ț (ns): Nh = nc + ns . Fie Tk, k = 1,2,….., Nh numărul de tăieturi care transformă 0 într-un domeniu simplu-conex . Frontiera lui '0 este , unde T'k si T"k sunt cele doua fețe ale tăieturii Tk (fig. 1.2). Dacă e necesar, notăm prin Tsk tăietura bobinei supraconductoare.

În '0 avem:

Pe suprafețele Si, V = ct. = Vi. Pe tăieturile Tk, V are un salt idk iar V/n are un flux magnetic dk. Avem:

(1.17)

și de asemenea:

+

+

Cel de-al treilea termen e nul datorită condițiilor pe frontieră () și ().

Dată fiind conservarea componentei tangențiale a lui H la suprafețele corpurilor conductoare, potențialul magnetic scalar este continuu la trecerea prin aceste suprafețe. De asemenea componenta normală a lui B se conservă. Ca urmare prima integrală din ultimul termen poate trece pe fața dinspre conductor a suprafeței i . A doua integrală din ultimul termen este chiar fluxul magnetic de pe suprafața tăieturii Ti . Luând în considerare forma integrală a legii inducției electromagnetice și ținând cont de relația (2.16), ultimul termen devine:

unde Ci este frontiera suprafeței ițTi (fig.1.2).

Deoarece:

doar termenul se păstrează. De aceea:

Utilizând relația (1.17) , obținem:

(1.18)

De aceea Bd = 0, E = 0. QED.

Observație importantă. Din relația (1.18) rezultă că pentru 2 câmpuri electromagnetice care au aceleași condiții inițiale și aceleași condiții de frontieră este valabilă relația:

(1.19)

sau:

Fig. 1.2. Domeniul conductor

În mod explicit nu apare viteza în raționamentele de mai sus. Este necesar ca vitezele să fie impuse, pentru că mărimile de câmp sunt definite în sisteme de referință locale asociate particulelor în mișcarea lor. Pentru ca și câmpul diferență să se refere la puncte bine definite, este necesar ca cele 2 câmpuri care se scad să fie asociate la aceleași viteze de deplasare ale corpurilor.

Consecință. Presupunem date condițiile de frontieră, condițiile ințiale și densitatea de curent impusă J0 ; la o magnetizație M rezultă unic câmpul electromagnetic și deci inducția magnetică B. Este deci bine definită funcția:

MB = (M)

Problema de curenți turbionari este bine formulată dacă este asigurată și existența soluției câmpului electromagnetic în regimul cuasistaționar. Din păcate restricțiile necesare aplicării operatorului în sens clasic impune clasa C1 pentru mărimile de câmp. Este posibil ca soluția să nu existe în această clasă. În [1] se arată că este necesară considerarea claselor de funcții de tipul L2 și sunt formulate teoreme de unicitate și existență a soluției de regim cuasistaționar și pentru regimul cuasistaționar periodic. Sunt luate în considerare mediile neliniare. Sunt descrise și procedurile de construire a soluțiilor.

2. Soluții numerice de ANALIZă a CÂMPULUI ELECTROMAGNETIC CUASISTAȚIONAR IN PROBLEME DE CURENTI TURBIONARI

Literatura de specialitate a oferit numeroase rezultate privind determinarea câmpului electromagnetic în regim cuasistaționar, în medii imobile, liniare sau neliniare [2], dar în structuri 2-D. Dată fiind uriașa dezvoltare a performanțelor mijloacelor de calcul din ultima vreme, atenția cercetătorilor s-a îndreptat spre elaborarea unor algoritmi de soluționare numerică a problemei curenților turbionari în structuri 3-D. La început au fost adoptate elementele nodale vectoriale. Elementele de muchie s-au dovedit mult mai potrivite pentru problemele de câmp în care intervine operatorul rot. În prezent aproape toate colectivele de ceretare din domeniul electromagnetismului utilizează elementele de muchie. Necunoscutele utilizate putând fi potențialele magnetic scalar V și vector A sau potențialele electric scalar si vector T. Se pot utiliza deasemenea direct și mărinile de câmp: intensitatea câmpului electric E, a câmpului magnetic H sau densitatea de curent J. Pentru a asigura unicitatea soluției sunt propuse condițiile de etalonare Coulomb, Lorentz sau condiția de etalonare topologică arbore-coarbore. Dacă se folosesc potențialele atunci se poate renunța la unicitatea soluției și deci la condiția de etalonare. Sistemul de ecuații se soluționează prin tehnici de gradient conjugat și în ultimul timp sunt raportate rezultate interesante în acest sens. Dacă se folosesc ca necunoscute H sau E, atunci este necesar să se impună condiția de divergență nulă pentru inducția magnetică B = H sau densitatea curentului total J + D/t = E + ,E/t. O metodă originală în care sunt folosite elementele de muchie ale corzilor și elementelor nodale a fost dezvoltată în cadrul Facultății de Electrotehică din UPB.

În acest capitol sunt prezentate pe scurt principalele tehnici de analiză numerică a câmpului electromagnetic cuasistaționar in problemele de curenți turbionari[5], [11].

Pornim la ecuațiile câmpului electromagnetic în regimul cuasistaționar pentru domeniul prezentat în figura 1.1, avem:

În T:

, (2.1)

H = F(B). (2.2)

În :

, (2.3)

, (2.4)

. (2.5)

În 0:

, (2.6)

. (2.7)

Unde J0 este densitatea de curent impusă. Bobinele supraconductoare au o valoare impusă pentru J0 care trebuie să fie luată uniform în secțiunea transversală. Fie J densitatea curenților turbionari. Ea e localizată în toate corpurile conductoare care nu au bobine. Întregul spațiu are permeabilitatea magnetică în vid 0.

2.1. Modelul (A-V, A)

Din legea fluxului magnetic rezultă că putem utiliza potențialul vector A astfel încât

1 (2.8)

Atunci din legea inducției rezultă:

2 (2.9)

Din (2.8) și din teorema lui Ampère rezultă:

Apoi, din legea conducției:

unde În consecință avem:

3 (2.10)

0 (2.11)

Din ecuația (2.10) rezultă:

4 (2.12)

Deoarece ecuațiile (2.10) și (2.11) nu asigură unicitatea potențialelor A și V: se poate adăuga la A și scădea din gradV gradientul unui potențial cu derivată după normală nulă pe S'. [11]. Din (2.8) rezultă că pe componenta tangențială a lui A diferă printr-un gradient. Putem chiar forța conservarea componentei tangențiale. Atunci, pentru a conserva componenta tangențială a lui E se impune ca și Vt să se conserve. Adică V este continuă pe T. Dar cel mai comod este să impunem ca V să fie definită doar pe . La suprafața corpului conductor componenta normală a densității de curent este nulă, deci:

5 (2.13)

Relația (2.13) rezultă din (2.11), (2.10) și continuitatea componentei tangențiale a lui rotA. Într-adevăr, din (2.11) avem:

Apoi, din (2.10) avem:

Pentru asigurarea unicității soluției (A,V) trebuie impusă o condiție de etalonare suplimentară pentru potențialul vector A. În plus, pe frontiera T (T presupus simplu conex) se dau:

(FRA) () n rot A = g pe S'

(ß) An = 0 pe S'

() At = f pe S"

() dacă S" este formată din n suprafețe disjuncte Sk, atunci pe n-1 suprafețe se dau fluxurile nule ale lui A. O condiție de etalonare foarte des utilizată este condiția de etalonare Coulomb:

6 (2.14)

care, pentru suprafețe de discontinuitate, conservă componenta normală. Atunci A este continuu pe T. Deoarece B este unic determinat, din relațiile (2.8), (2.14) și condițiile de frontieră (FRA) rezultă că A este unic determinat. Din (2.13) rezultă că pe :

7 (2.15)

Din ecuația (2.12) rezultă valoarea lui divgradV. Putem da într-un punct din valoarea lui V și ca urmare și V este unic determinat în . Rezultă deci:

Teorema 2.1 Potențialele A,V care verifică ecuațiile (2.10), (2.11), (2.14) și condițiile de frontieră (FRA), (2.15) sunt unic determinate în T.

Demonstrație Pentru rezolvarea numerică a problemei se folosește tehnica Galerkin, proiectându-se ecuațiile (2.10), (2.11), (2.14) pe o mulțime de funcții Nk, k, Vk. Funcțiile Nk au condiția de frontieră () nulă și rotNk sunt liniar independente. Funcțiile k au valori constant arbitrare pe Si, i=1,2,…,n-1 și valori nule pe Sn. Funcțiile Nk și k sunt definite pe T, iar Vk pe . Se scriu:

,

unde N0 îndeplinește condiția de frontieră ().

Proiectând ecuațiile (2.10), (2.11), condiția de frontieră () și condiția de trecere (2.4) pe Nk, k1, avem:

8 (2.16)

Într-adevăr, integrând prin părți primul termen rezultă:

9

(2.17)

Proiectăm acum ecuațiile (2.12) și (2.13) pe funcțiile Vk și se obține:

10 (2.18)

Într-adevăr, integrând prin părți primul termen, rezultă:

11 (2.18')

Relația (2.18) se poate obține și prin proiectarea relației (2.16) pe gradVk, deoarece:

Din (2.17) rezultă conservarea componentei tangențiale a lui rotA și verificarea ecuației (2.11). Atunci, componenta normală a lui rotrotA este nulă și rezultă (2.18'). Mai rămâne să proiectăm relația (2.14) și condițiile de frontieră (ß), () pe k, de unde rezultă:

12 (2.19)

Într-adevăr, integrând prin părți rezultă:

unde ki sunt valorile arbitrare ale lui k pe Si, i=1,2,…,n-1.

Sistemul (2.16), (2.18), (2.19) de nE = nN+nV+n ecuații cu nE necunoscute este folosit pentru cazul în care adoptăm elementele nodale și dorim să fie verificată condiția de etalonare, precum și condiția de frontieră (ß), (). Funcțiile Nk obținute cu elemente nodale sau elemente de suprafață nu verifică condiția divNk=0. Elementele de muchie verifică condiția de etalonare în interiorul domeniilor, nu însă și pe interfețe, unde componentele normale ale lui Nk nu se conservă. Totuși aproximația Aa cu elemente de muchie a unei funcții A cu divA=0 are proprietatea că atunci când norma divizării domeniului tinde către zero, AaA. Este motivul pentru care de multe ori se renunță la condiția de etalonare în cazul elementelor de muchie. Verificarea (slabă) a condiției de etalonare necesită adăugarea funcțiilor grad .

În calculele numerice se poate renunța la condiția de etalonare impunând ca A să fie o combinație de funcții Nk. În acest caz, faptul că A este o combinație de funcții Nk apare ca o condiție de etalonare. Rămâne doar sistemul (2.16), (2.18) de nE = nN+nV ecuații cu nE necunoscute.

Condiția de etalonare poate fi plasată în ecuația (2.10), transformând-o în:

13 (2.20)

dar relația (2.12) nu mai rezultă și ea trebuie impusă. Relația (2.12) are semnificația fizică divJ=0. Evident din (2.12) și (2.20) rezultă:

14 (2.21)

și dacă pe avem divA = 0, rezultă că divA = 0 în tot . La fel relația (2.11) poate fi modificată în:

15 (2.22)

ținând cont că div J0 = 0. Proiectăm ecuațiile (2.20), (2.12), condiția de frontieră () și condiția de trecere (2.4) pe Nk și obținem:

16 (2.23)

Într-adevăr, integrând prin părți primii doi termeni, avem:

17

(2.24)

Al treilea termen impune continuitatea lui divA, iar al 5-lea termen impune condiția de frontieră nulă pentru divA.

Se adaugă ecuația (2.18). Obținem un sistem de nE = nN + nV ecuații cu nE necunoscute, în care este inclusă și condiția de etalonare (sub forma slabă). Evident sistemul ecuațiilor (2.23), (2.18) este mai agreabil decât (2.16), (2.18), (2.19). Dacă adoptăm elementele de muchie atunci condiția de etalonare este implicită și rezultă sistemul format din ecuațiile (2.16), (2.18).

Mai rar se folosește condiția de etalonare Lorentz:

18 (2.25)

Atunci ecuația (2.10) devine:

19 (2.26)

Apar însă dificultăți atunci când mediul conductor nu este omogen electric.

Observații:

1) Procedura (A-V, A) este aplicabilă și în cazul în care este multiplu conex (T este însă simplu conex).

2) Folosirea potențialului vector A în nu este economică în comparație cu potențialul magnetic scalar redus (§2.2).

2.2. Modelul (A – V,)

Se presupunem că 0 este simplu conex. Rezolvând problema de regim staționar putem obține în 0 câmpul magnetic (B0, H0) astfel încât rotH0 = J0, divB0=0 și B0 = H0, iar pe 0 să avem (vezi §1.3):

Câmpul din 0 este (BT, HT) = (B0, H0) + (B, H) și rotH = 0, deci H = -grad . Condițiile de frontieră de tip longitudinal dau pentru câmpul (B, H) condițiile:

Atunci avem ecuațiile:

20 (2.27)

(2.10)

cu condițiile pe :

(2.13)

(2.28)

(2.29)

Putem include condiția de etalonare în ecuația (2.10) și atunci ecuațiile (2.10) și (2.14) se înlocuiesc cu (2.20) și (2.12). În calculele numerice se folosesc funcțiile test Nn definite doar pe . Potențialul vector A se scrie:

Proiectând de exemplu ecuația (2.20) și condiția de trecere (2.29) pe Nk obținem:

21 (2.30)

Într-adevăr, integrând prin părți primii doi termeni, avem:

Ultimul termen impune condiția de frontieră nulă pentru div A pe .

La relația (2.26) adăugăm relația (2.18).

Pentru folosim funcțiile test k definite pe 0:

unde k are valori constante arbitrare pe Si, i=1,2,…,n-1 și valori nule pe Sn.

Proiectăm ecuația (2.27) și condiția de trecere (2.13) pe k, obținem:

22 (2.31)

Într-adevăr, integrând prin părți primul termen rezultă:

Ecuațiile (2.30), (2.18), (2.31) formează un sistem prin rezolvarea căruia obținem A, V, .

În domeniul 0 este convenabil să folosim elementele nodale. Suprafeței Si îi asociem o funcție test care rezultă din însumarea tuturor elementelor nodale k corespunzătoare nodurilor de pe suprafața Si, deci pentru suprafața Si avem:

unde {i} este mulțimea nodurilor de pe Si.

În se pot alege elemente nodale sau elemente de muchie. În ultimul caz condiția divA este asigurată. Atunci al doilea termen din ecuația (2.30) dispare.

Observații:

1) Procedura (A-V,) este mai economică pentru memoria calculatorului decât procedura (A-V,A), pentru că în 0 necunoscutele sunt scalare.

2) Aplicarea procedurii (A-V,) pentru domenii multiplu conexe este foarte dificilă (nu putem defini potențialul ). Este necesar să separăm în 0 câmpuri magnetice Hk care să aibă componenta normală nulă pe și circulația unitară pe holul k [12]. Câmpul magnetic pe 0 se scrie:

Construirea câmpurilor Hk este însă greoaie.

2.3. Modelul (T –

Deoarece div J = 0, rezultă că putem folosi potențialul electric vector T astfel încât:

23 (2.32)

De fapt T este un câmp magnetic care verifică doar teorema lui Ampère. Deci el diferă de câmpul magnetic real printr-un gradient. Este de fapt un raționament asemănător cu cel de la procedeul potențialului magnetic redus (§2.2). Deci:

24 (2.33)

Din legea inducției electromagnetice și din (2.32) rezultă:

25 (2.34)

În 0 T poate fi chiar H0 definit în §2.2, dar având H0t = 0 pe . Deci în 0 avem ecuația:

26 (2.27)

Relația (2.34) seamănă foarte mult cu (2.10) și raționamentul ce va urma va fi asemănător cu cel de la §2.2. Din relația (2.34) rezultă:

27 (2.35)

adică legea fluxului magnetic. Adăugăm condiția de etalonare:

28 (2.36)

În locul ecuațiilor (2.34) și (2.36) putem pune ecuația (2.35) și:

29 (2.37)

Pe frontiera unde jn= 0 condiția de frontieră impusă lui T este:

30 (2.38)

admițând că este simplu conex. Atunci continuitatea componentei tangențiale a lui H pe conduce la continuitatea potențialului scalar prin . Putem pune atunci . Continuitatea componentei normale a lui B prin este:

31 (2.39)

Unicitatea potențialului electric vector rezultă din relația (2.32) știind că J este unic determinat (Teorema de unicitate), din condiția de etalonare (2.36) și din condiția de frontieră (2.38). Ca urmare este unic determinat de relațiile (2.27), (2.33), de condițiile de trecere prin (2.39) și de condițiile de frontieră de pe T .

Pentru calculele numerice definim în funcțiile Nk care îndeplinesc condiția de frontieră (2.38) și punem:

Proiectăm (2.34) sau (2.37) pe Nk. De exemplu, pentru (2.37) avem, după ce folosim Gauss-Ostrogradski și condiția (2.38):

32 (2.40)

Al doilea termen din relația (2.40) impune condiția de frontieră nulă pentru divT.

Alegem acum funcțiile k în T cu valori constante ki pe Si, i=1,2,…,n-1 și nule pe Sn. Scriem:

Proiectând relațiile (2.35), (2.27) și (2.39) pe funcțiile k, obținem:

33,

(2.41)

Într-adevăr, integrând prin părți primul termen, avem:

34

(2.42)

Ecuațiile (2.40), (2.41) formează un sistem de ecuații prin rezolvarea căruia obținem T și .

Observații:

1) Metoda este convenabilă din punct de vedere al memoriei calculatorului.

2) Pentru domenii multiplu conexe metoda este foarte dificil de aplicat

3. SOLUȚIONAREA NUMERICĂ A PROBLEMEI DE CÂMP TERMIC

3.1. Aspecte calitattive. Difuzia câmpului termic

Se consideră domeniul cu bordura .. Ecuația Fourier pentru regimul staționar al câmpului termic este:

-div grad T = p (3.1)

unde este conductibilitatea termică iar p densitatea de volum a puterii ce se transformă din forma electromagnetică în căldură, în efectul Joule p = .

Condiția de frontieră este :

(3.2)

unde este coeficientul de convecție termică iar Te este temperatura în exteriorul domeniului .

Ecuația (3.1) are soluție unică în condiția de frontieră (3.2). Într-adevăr dacă două câmpuri termice verifică ecuația (3.1) și condiția de frontieră (3.2) , atunci câmpul diferență Td verifică forma omogenă a acestor relații. Avem:

(3.3)

În același timp,

(3.4)

Din relațiile (3.3) și (3.4) rezultă:

= 0 (3.5)

de unde: .

Difuzia câmpului termic este descrisă de ecuația:

– div grad T + c (3.6)

unde c este capacitatea termică.

Dacă se dă valoarea inițială a temperaturii, atunci ecuația (3.6) cu condiția de frontieră (3.2) are soluție unică. Într-adevăr, în locul relației (3.4) avem:

= (3.7)

Din (3.3) și (3.7) rezultă, prin integrare :

(3.8)

de unde Td = 0 .

Proprietăți ale operatorului L = – divgrad

Operatorul L = – divgrad, definit pe mulțimea H a câmpului de clasă C2() cu condiția de frontieră nulă (Te=0) este simetric și pozitiv [5]. Într-adevăr:

(3.9)

și

= (3.10)

Din (3.9) și (3.10) rezultă:

< T’ , LT” > = (3.11)

Membrul drept fiind simetric în raport cu T’ și T” , rezultă:

<T’,LT”> = < T”,LT’>

Operatorul L este pozitiv:

<T,LT> 0

egalitatea având loc dacă și numai dacă T = 0. Într-adevăr, punând în (3.11) T’=T”=T rezultă:

< T,LT > =

Considerăm că L este chiar pozitiv definit, afirmație demonstrată în multe cazuri de condiții de frontieră:

< T, LT>, cu (3.12)

Dacă L este simetric și pozitiv definit, atunci putem defini produsul scalar

Fie H închiderea lui H în norma definită de acest produs scalar.

H este subspațiu al lui L2(). Într-adevăr, din relația (3.12) rezultă că

unde cu am notat norma generată de produsul scalar. Deci șirul Cauchy din H este Cauchy și în L2() q.e.d.

definită pe H are semnificație în teoria distributiilor, adică:

(3.13)

Mai spunem că T din H este soluție slabă pentru ecuația (3.1).

Oricare ar fi pL2() există T H soluție slabă pentru ecuația (3.1).

Într-adevăr, pentru funcționala U: H R avem:

conform relației (3.13). Deci funcționala U este liniară și mărginită pe H. Atunci există T astfel încât:

Deci:

adică (3.13).

Soluționarea problemei de câmp termic prin metoda elementului finit

În practică au fost elaborate mai multe proceduri de analiză numerică, însă cele mai importante proceduri dezvoltate în Metoda elementului finit sunt Ritz și Galerkin [11].

Procedura Ritz

Ecuației:

-div grad T = LT=p (3.14)

i se asociază funcționala:

(3.15)

unde T este câmp termic cu condiții de frontieră nule (). Spunem că .

Deoarece operatorul L este pozitiv, funcționala F este convexă. Într-adevăr, să arătăm că:

, pentru (3.16)

sau:

(3.17)

Dezvoltând termenii din membrul stâng al relației (3.17), rezultă:

(3.18)

sau:

(3.19)

Minimul funcționalei (3.15) este unic (fiind convexă) și reprezintă soluția ecuației (3.14). Într-adevăr, fie minimul și fie o temperatură oarecare cu si arbitrari. Avem, ținând cont de simetria operatorului L:

=

+ (3.20)

sau:

>0 pentru orice (3.21)

Membrul drept al relației (3.21) poate fi privit ca un trinom în . Inegalitate poate fi realizată pentru orice dacă și numai dacă termenul care-l conține pe este nul:

(3.22)

Egalitatea (3.22) este valabilă pentru orice . Atunci:

(3.23)

adică este soluție a ecuației (3.14).

Reciproc, soluția ecuației (3.14) este minimul funcționalei (3.15). Într-adevăr, dacă (3.23) este satisfăcută, atunci primul termen din membrul drept al relației (3.21) este nul.

În procedura Ritz, se caută minimul funcționalei (3.15), în loc să se rezolve ecuația (3.14). Din păcate H este infinit dimensional și atunci se minimizează (3.15) pe subspații finit dimensionale ale lui H. Dacă pentru subspațiile finit dimensionale avem șirul de incluziuni:

(3.24)

atunci pentru temperaturile ce realizează minimele funcționalelor din aceste subspații este, evident, valabilă relația:

(3.25)

Relația (3.25) arată că Procedura Ritz este o metodă numerică consistentă: valorile ce minimizează în subspațiile (3.24) funcționala (3.15) se apropie de soluția ecuației (3.14). Totuși, dacă șirul subspațiilor (3.24) “tinde” către un subspațiu al lui H, diferit de H, soluțiile nu tind obligatoriu către .

Există câteva proceduri consacrate de alegere a subspațiilor . Cea mai cunoscută este cea a spațiului generat de elementele nodale. Cele mai simple elemente nodale sunt cele de ordinul 1. Se definește pe domeniul o rețea de triunghiuri (Figura 3.1). Elementul nodal de ordinul 1, asociat nodului , este funcția continuă și cu variație liniară pe , care are valoare 1 în nodul și valori nule în celelalte noduri. Graficul unei astfel de funcții arată că în Figura 3.2.

Se scrie temperatura sub forma:

(3.26)

unde sunt funcții liniar independente cu condiții de frontieră nule (eventual, elemente nodale), iar se determină astfel încât funcționala (3.15) să fie minimă. Derivând în raport cu și ținând cont de simetria operatorului L, avem:

, k=1,2,…,n (3.27)

Formularea de mai sus are două dezavantaje: funcțiile trebuie să fie de două ori derivabile și condițiile de frontieră trebuie să fie nule. Să presupunem că frontiera a domeniului este formată din 3 suprafețe , , , unde sunt îndeplinite condițiile de frontieră Dirichlet: , Neumann: și mixte – (3.15). Ultimele 2 condiții de frontieră pot fi considerate cazuri particulare ale condițiilor mixte, atunci când =0 și respectiv =0. Din punct de vedere tehnic, condiția de frontieră Neumann apare pe suprafețe de simetrie (suprafețe de câmp), fiind întotdeauna nulă. Se înlocuiește funcționala (3.15) cu:

(3.28)

iar minimul se caută în mulțimea funcțiilor care au condiția de frontieră Dirichlet impusă. Într-adevăr, fie minimul funcționalei (3.28) și fie o temperatură oarecare, dar care îndeplinește condițiile de frontieră Dirichlet, cu și arbitrari, având pe . Înlocuind în (3.28), avem:

(3.29)

Integrăm prin părți ultimul termen:

=

(3.30)

Înlocuind în (3.29) și grupând după puterile lui , avem:

+>0 (3.31)

Dacă minimizează funcționala F′, atunci inegalitatea (3.31) poate avea loc pentru orice dacă și numai dacă termenul care-l conține pe este nul. Cum T este arbitrar, din primul termen al acoladei rezulta condiția de frontieră (3.2), din al doilea, condiția Neumann, iar din al treilea, ecuația (3.14). Reciproc, dacă cele două condiții de frontieră și ecuația (3.14) sunt îndeplinite, atunci termenul în este nul, inegalitatea (3.31) este îndeplinită și este minimul funcționalei F′.

Pentru calculele numerice scriem câmpul de temperaturi T sub forma:

(3.32)

unde sunt funcții liniar independente cu condiții de frontieră Dirichlet nule iar o funcție arbitrară care îndeplinește condiția de frontieră Dirichlet. Ea poate fi pusă sub forma:

(3.33)

unde sunt nodurile de frontieră Dirichlet, – elementele nodale ale acestor noduri, iar numărul de noduri de pe frontiera Dirichlet. Înlocuind în (3.28) și derivând în raport cu necunoscutele rezultă sistemul:

,

k=1,2,…,n (3.34)

4. CUPLAREA PROBLEMELOR DE CÂMP ELECTROMAGNETIC ȘI TERMIC

Datorită neliniarității caracteristicii B-H a materialului utilizat în procesul de încălzire prin inducție electromagnetică se impune cuplararea problemei de câmp electromagnetic cu cea de difuzie termică. Soluționarea neliniarității poate fi realizată prin adoptarea unui model pseudo-liniar, pentru care relația B-H este liniară: , permeabilitatea magnetică fiind corectată funcție de valoarea efectivă a inducției magnetice. Deoarece încălzirea se face peste punctul Curie, caracteristica B-H depinde puternic de temperatură [12], aceasta modificându-se până la cea a vidului, soluția de câmp electromagetic este influențată putermic de temperatură, piederile prin curenților turbionari reprezentând sursa de caldură pentru problema de difuzie termică [5]. Problema termică este la rândul său neliniară, datorită dependenței conductibilității termice de temperatură, în timp ce pierderile specifice Joule depind de pătratul densității de curent și de rezistivitate, ce sunt dependente de tempertură.

Cuplarea problemelor de câmp electromagnetic și cea de difuziune termică se face conform schemei logice prezentate în fig. 5.1. Modelul prezentat este unul adoptat de majoritatea pachetelor de programe software dedicate calculului problemei cuplate de câmp electromagnetic și termic [16].

1) Inițializarea magnetică

Este rezolvată problema magnetică fiind calculată distribuția densității puterii inițiale, cunoscând valoarea inițiala a câmpului de temperaturi la timpul (o notăm cu ). Problema de câmp fiind pseudoliniară: se corectează iterativ valoarea permeabilității funcție de inducția magnetică. Se determină astfel densitățile pierderilor Joule necesare în problema termică.

2) Prima rezolvare termică

Se alege pasul de timp inițial (cel folosit a fost = 0,01 s.) respectiv se impune limita timpului de studiu. Cunoscând tempertura , se determină, prin procedura Crank-Nicholson câmpul de temperaturi la timpul și câmpul de temperaturi , cu care se face corecția parametrilor termici.

3) Actualizarea proprietăților magnetice

La finalizarea calculului primei hărți de temperatură pentru pasul de timp curent, se actualizează proprietățile magnetice funcție de câmpul de temperaturi și repornește calculul magnetic. În continuare, pornind de la harta densității puterii actualizată, se actualizează și proprietățile termice și se calculează noua hartă a temperaturii, obținându-se si .

Fig. 4.1. Organigrama de funcționare a programului

Caracteristica B-H se ajustează funcție de valorile inducției magnetice și de câmpul de temperaturi , după care se actualizează din nou parametrii termici și se detemină noua valoare a câmpului de temperaturi . Valoarea obținută se compară cu cea anterioară: dacă diferența relativă este prea mare, continuă procesul iterativ de la pct. 3. Când diferența este suficient de mică, se trece la pasul următor de timp, evaluând, pentru început, câmpul termic, ca la pct. 2, apoi urmând procesul iterativ de la pct. 3.

5. Simularea numerică a procesului de incălzire prin inducție electromagnetică

Instalațiile de încălzire prin inducție utilizate în vederea realizării tratamentelor termice sunt realizate într-o gamă foarte largă, funcție de parametrii tehnologici ai procesului de producție: productivitatea, proprietățile de material, forma geometrică a semifabricatului și altele. Ponderea instalațiilor de încălzire prin inducție electromagnetică în gama instalațiilor electrotermice a crescut considerabil în ultimii ani, iar la nivel mondial se evidențiază o serie de tendințe precum utilizarea într-o măsură mai mică a frecvenței industriale, cele mai multe instalații de încălzire prin inducție având frecvența de lucru în gama 8000 – 200000 Hz, datorită avantajelor pe care le are utilizarea frecvențelor medii și cererea crescândă pentru soluții moderne, eficiente din punct de vedere al procesului tehnologic. Dintre avantajele prezentatre de procesarea prin inducție electromagnetică a materialelor putem aminti: instalații relativ simple, fiind posibil lucrul în vid sau atmosferă controlată, gamă foarte largă de puteri, căldura se dezvoltă în corpul care urmează a fi încălzit, viteză mare de încălzire (timp de încălzire redus), obținerea unor temperaturi ridicate cu o inerție mică, reglarea precisă a temperaturii, posibilitate de automatizare a procesului de încălzire și prelucrare, poluare redusă a mediului ambiant, spațiul ocupat este relativ redus.

5.1.Definirea problemei.

Aplicația considerată este reprezentată de procesul tehnologic de încălzire prin inducție elctromagnetică a unui semifabricat cilindric ce urmează a fi deformat la cald. Reperul supus analizei face parte din rândul componentelor utilizate în industria auto, realizarea sa reprezentând o problem de interes deoarece acesta va trebui să indeplinească o serie de condiții ce implică siguranță sporită în funcționare. Analiza numerică va fi realizată în concordanță cu parametrii instalației de încălzire prin inducție electromagnetică, având frecvența tensiunii de alimentare de 8 KHz. Inductorul având secțiune rectangulară de 12 X 9 [mm] și un număr de trei spire, este realizat astfel încât să asigure o distribuție uniformă a curentului în piesă, respectiv o distribuție unifomă a câmpului termic obținut prin încălzirea inductivă a zonei de interes. Forma geometriă a ansamblului inductor-piesă considerat este reprezentată în figura 5.1. Rezultatele obținute vor pune în evidență influența frecvenței asupra timpului necesar pentru atingerea temperaturii dorite, respectiv influentă acesteia asupra parametrilor energetici ai sistemului. Din motive de simetrie axială, problema va fi analizată pentru un sfert din intregul domeniu de calcul.

Materialul utilizat la realizarea semifabricatuluil este oțel AISI-SAE-420, cu următoarele proprietăți fizice:

– Rezistivitatea  [m], a cărei dependentă de temperatura  [C] este data în Figura 5.2.

Dependenta () poate fi aproximată prin funcția [10]:

[m] = 8,69779674E-13 2 + 4,78235404E-10 + 1,56565488E-07,

pt. [20, 760]C,

și

[m] = – 2,01190476E-13 2 + 7,12976190E-10 + 6,54785714E-07,

pt. [760, 1400]C

Relatia B-H, a cărei dependență de temperatură este descrisă prin relația [10]:

B(H, ) = 0H + (2/)Js0arctg [(r0 – 1) 0H/(2Js0)]{1-exp[( – c ) / C]}

unde:

0 = 4E-7 H/m este permeabilitatea magnetică a vidului,

Js0 = 1,95 T este inducția magnetică de saturație la 0 C,

r0 = 100 este valoarea inițială a permeabilității relative,

c = 760 C este temperatura Curie,

C = 130 C este o constantă.

Figura 5.2. Dependența de temperatură a rezistivității oțelului

– Capacitatea calorică c[J/m 3/C], a cărei dependență de temperatură   este dată în Figura 5.3.

Figura. 5.3. Dependența de temperatură a căldurii specifice a oțelului

– Conductivitatea termică izotropă  [W/m/C], a carei dependență de temperatură  este dată în Figura. 5.4. Următoarele funcții [16] pot aproxima această dependență:

[W/m/C] = – 0,0000233133 2 – 0,0142523022 + 52,4957597802, pt. 760 C,

și

[W/m/C] = 0,00778571 + 19,59285714, pt. > 760 C

– Transferul termic la suprafața piesei este caracterizat de un coeficient de convecție cu valoarea = 20 W/m2/C și de un coeficient de radiație (caracteristic transferului termic prin radiație) = 0,75 cu ajutorul cărora dependența coeficientului de transfer termic e funcție de temperatură poate fi descrisă prin [10]:

e () = + ·Cn [( +273) + (a+273)][( + 273)2 + (a+273) 2]

unde Cn = 5,73E-08 W/m2/C4. Curba e () este reprezentată în Figura 5.5.

Figura 5.4. Dependența funcție de temperatură a conductivității termice a oțelului

Temperatura de austenitizare este θa = 760 °C.

Condițiile de periodicitate și sensul curentului impun utilizarea următoarelor condiții de frontieră pentru problema de câmp electromagnetic: Dirichlet nulă pentru potențialul vector și axa de simetrie.

Domeniul de calcul pentru câmpul termic este cel al semifabricatului. Porțiunea de frontieră corespunzătoare centrului piesei va fi condiție de frontieră Neumann nulă, în timp ce, pe restul frontierei, condiția de frontieră este mixtă.

Figura 5,5. Dependența funcție de temperatură a coeficientului de transfer termic echivalent

5.2. Rezultate numerice

Calculul numeric, utilizând soluția de potential AV-A a fost realizat pentru problema descrisă mai sus în vederea determinării evoluției în timp a temperaturii pe durata procesului de încălzire prin inducție electromagnetică. La atingerea temperaturii de 930 C, calculul numeric este opritt, astfel încât să poată fi determinaă influența frecvenței de alimentareasupra duratei de timp necesară atingerii temperaturii și distribuția mărimilor de câmp în semifabricatul supus studiului. Va fi deasemenea analizată distribuția câmpului termic în semifabricat și la suprafața acestuia.

Rezultatele experimentale reprezentate în figurile 5.6 – 5.9 pun în evidentă distribuția fluxului magnetic în ansamblul inductor piesă, fiind evidențiată distribuția și valorile sale la diferite momente de timp în domeniul de calcul considerat.

Distribuția la suprafață a vectorului Intensitate a câmpului magnetic este prezentată în figura 5.10. Pentru a pune în evidență întensitatea câmpului magnetic si valorile sale în ansamblul inductor piesă au fost reprezentate figurile 5.11 -5.13, fiind astfel pusă în evidență influența modificării proprietăților de material cu temperature asupra marimilor de câmp.

In continuare sunt prezentate rezultate ce pun in evidență distribuția câmpului termic pe durata procesului de încălzire în câmp electromagnetic și modificarea proprietăților de material odata cu temperature. Astefl în figura 5.16 se reprezintă distribuția câmpului termic la începutul procesului de încălzire prin inducție electromagnetică (t=0.1 s), respectiv, în figura 5.17 se pune în valoare conductivității termice funcție de temperatură.

Rezultatele prezentate in figurile 5.14 ;I 5.15 pun în evidență distribuția desnsității de putere indusă în piesă pentru două momente de timp, si anume la începutul procesului si la finalizarea acestuia.

Pentru a pune în evidență căldura specifică dezvoltată în semifabricat, a fost reprezentat în figurile 5.24 și 5.25 distribuția acesteia în piesă în scopul determinării zonelor ce suferă ulterior modificarea proprietăților de material.

7. CONCLUZII

Avantajul încălzirii prin inducție electromagnetică a semifabricatelor în procesele de încălzire inductive constă în faptul că sursa de încălzire este controlată, putând fi folosită și în combinații cu alte procese tehnologice, permițând atingerea unei temperaturi uniforme pe secțiunea piesei supuse încălzirii. Deasemenea procesul este intuitiv, permițând repetarea sa pentru un număr nelimitat de ori.

Modelarea numerică a procesului de încălzire a semifabricatelor cu ajutorul curenților turbionari este o problemă complexă, în care sunt rezolvate, simultan, două probleme neliniare de câmp: cea de câmp electromagnetic și cea de difuzie termică. Neliniaritatea problemei de curenți turbionari este datorată relației neliniare B-H, în timp ce neliniaritatea problemei termice provine din dependența cu temperatura a parametrilor termici (conductibilitate termică, capacitate termică, coeficient de transfer termic la suprafață).

Cuplajul celor două probleme rezultă din puternica dependență a relației B-H cu temperatura, în problema de câmp electromagnetic și din sursa de câmp termic, dată de pierderile Joule, în problema de difuzie termică. Neliniaritatea relației B-H se soluționează admițând modelul pseudolinear, în care , permeabilitatea magnetică corectându-se iterativ în functie de inducția magnetică. Marele avantaj al modelului rezultă din posibilitatea adoptării regimului sinusoidal și a imagnilor în complex, forma numerică a ecuațiilor câmpului conducând la un sistem de ecuații algebrice cu coeficienți complecși.

Prin calcul numeric se permite ca, cuplajul dintre cele două probleme de câmp să fie soluționat prin procedura de discretizare în timp, unde, la fiecare pas de timp, se face corecția iterativă a parametrilor electromagetici și termici.

Modelul prezentat realizează soluționarea problemei de câmp electromagnetic prin metoda elementului finit. Se realizează o analiză calitativă a câmpului electromagnetic cvasistaționar și cvasistaționar sinusoidal, prezentându-se ecuațiile regimurilor, condițiile de frontieră de tip electric și magnetic.

Problema de difuziune termică este discretizată în timp prin procedura Crank-Nicholson, în timp ce discretizarea în spațiu se face prin metoda elementului finit.

Simularea pe calculator a proceselor, care apar în timpul încălzirii prin inducție eelctromagnetică, reprezintă o unealtă extrem de valoroasă. Ea permite selectarea și controlul multor parametrii care influențează procesul pornind de la: valorile intensității curenților, frecvențele de lucru, timpii de acționare, astfel încât să se obțină evoluția dorită pentru câmpul de temperaturi și, deci, forma și adâncimea zonei de interes.

Analiza pe calculator face posibilă scurtarea timpului și scăderea costului care ar fi necesar pentru selectarea parametrilor în timpul experimentelor. Temperaturile și timpul necesare pentru atingerea temperaturii specifice sunt similare cu cele obținute pentru procesele fizice.

Mărirea productivității și creșterea eficienței economice a proceselor electrotermice presupun cunoașterea cu exactitate a mărimilor ce intervin. Din aceste motive proiectarea, construcția și exploatarea rațională a instalațiilor electrotermice implică cunoașterea tuturor factorilor care influențează indicatorii energetici și economici în scopul optimizării proceselor tehnologice atât din punct de vedere al calității produselor, cât și din cel al consumului de energie electrică.

Bibliografie

Adly A. A. (1996) Solution of induction heating problems involving media with hysteresis, J. Appl. Phys., 78, pp.4675 – 4677;

Albanese R., Hănțilă F. I., Rubinacci G. (1996) A Nonlinear Eddy Current Integral Formulation in Terms of a Two-Component Current Density Vector Potențial, IEEE Trans. on Magn., Vol. 32, Nr. 3, pp 784–787

4. Bay F., Labbe V., Favennec Y., Chenot J. L. (2003) A numerical model for induction heating processes coupling electromagnetism and thermomechanics, Int. J. Num. Meth. in Engng., 58, pp. 839-867;

Bensaid S., Trichet D., Fouladgar J. (2005) 3-D simulation of induction heating of anisotropic composite materials, IEEE Trans. Magnetics, 41, pp. 1568-1571;

Firețeanu V., Leuca T. (1997) Inducția electromagnetică și tehnologii specifice, Editura Mediamira, Cluj – Napoca;

Comșa D., Pantelimon L. (1979) Electrotermie, Editura Didactică și Pedagocică, București;

Comșa Dan (1986) Instalații electrotermice industriale.Volum II, Editura tehnică, București;

Firețeanu V. (1995) Procesarea electromagnetică a materialelor, Editura Politehnică, București;

34. Galunin S., Zlobina M., Blinov Yu., Nacke B., et al (2004) Numerical optimization in design of induction heating systems, Int. Symp. on Heating by EM Sources HES-04, June 23-25, Padua Italy;

Golovanov N., Șora I., și alți (1997) Electrotermie și electrotehnologii, Editura Tehnică, București;

50. Hănțilă F., Preda G., Vasiliu M., Leuca T., Della Giacomo E. (2001) Calculul numeric al curenților turbionari, Editura ICPE, 2001, ISBN 973 – 8067 – 31 – 6;

F. I. Hănțilă, Mathematical Models of the relation between B and H, Rev. Roum. Sci. Techn. – Électrotechn. et. Énerg., 19,3, pp.429-448, Bucharest, 1974, ISSN-1223-2106

Lavers J. D. (1983) Numerical solution methods for electroheat problems, IEEE Trans. Mag., 19, pp. 2566-2572;

Leuca T., Șt. Nagy, Carmen Molnar (2000) Încălzirea inductivă. Tehnici de proiectare asistată de calculator, Treira, Oradea;

Maghiar T., Hoble D., Nagy Șt. (2001) Method of increasing of the Electroenergetic Conversion Efficiency, Using Magnetodielectric Field Concentrator, HIS01;

Maxwell 2D SV: “ User Guide”