Introducere 2 Formularea problemei de câmp electromagnetic cvasi-stationar. 4 Formularea problemei de câmp termic 13 Aspecte calitative 13… [311223]

CURPINS

Introducere 2

[anonimizat]. 4

Formularea problemei de câmp termic 13

Aspecte calitative 13

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

[anonimizat] 25

Aplicații numerice și rezultate 27

Concluzii 38

Bibliografie 42

INTRODUCERE

Încălzirea prin inducție electromagnetică se bazează pe efectul de pătrundere a [anonimizat]. Curenții turbionari determinați de tensiunile electromotoare induse conduc la încălzirea acestuia prin efectul Joule.

Procesul de încălzire prin inducție implică două fenomene cuplate: 1 – transferul de energie pe cale electromagnetică de la bobină la corpul de încălzit; 2 – [anonimizat], prin efect Joule și transferul de căldură prin conducție termică în întreaga masă a corpului. Cele două fenomene prezintă fiecare dimensiuni de timp diferite. Datorită fenomenelor complexe ce apar, o etapă importantă în realizarea echipamentelor inductive o [anonimizat], modelarea numeric a [anonimizat].

[anonimizat] [1]. [anonimizat], [anonimizat] – Mecanica reprezintă cadrul teoretic ce trebuie luat în considerare de modelele folosite de sistemele CAD în electrotehnică [2].

Ca urmare a progreselor industriale în domeniul materialelor și a schimbărilor fundamentale în structura și funcționarea sistemelor electrotehnice se impune o continuă revizuire și îmbunătățire a metodelor tradiționale de concepție și de calcul. [anonimizat], considerarea regimurilor tranzitorii au impus metode și tehnici de concepție asistată de calculator [3], care răspund eficient atât necesităților de optimizare a structurilor clasice cât și acelora de predeterminare a structurilor noi înainte de construcția unor prototipuri fiabile.

Proiectarea asistată de calculator poate avea ca efect cumulativ pe termen lung dezvoltarea creativității prin posibilitatea de a experimenta ușor prin simulare a unoi idei noi și de a [anonimizat] “laborator” de construcție de prototipuri (fără costurile și termenele unei fabricații veritabile) și un partener de concepție în luarea deciziilor.

Modelele utilizate pentru simularea funcționării dispozitivelor electrotehnice sunt descrise de ecuații cu derivate partiale în domenii adesea cu geometrie foarte complex, fac apel la tehnici de discretizare a domeniilor de calcul [5].

Ecuațiile cu derivate parțiale se transformă în modelarea numerică în sisteme de ecuații algebrice, a căror soluție furnizează o aproximație a câmpului într-un număr discret de puncte ale domeniului de calcul [6].

Metoda elementelor finite permite determinarea unei funcții globale ce reprezintã fenomenul de studiat în orice punct al domeniului de calcul. Deoarece domeniu de calcul este divizat în multiple subdomenii adiacente, funcția globalã este la rândul său un ansamblu de funcții asociate fiecãruia din aceste elemente, discretizarea domeniului de calcul fiind efectuată în general cu elemente regrupate pe familii topologice ca: segmente de dreaptã, triunghiuri, patrulatere, tetraedre, paralelipipede, prisme, etc, funcția ce modeleazã fenomenul fiind definitã printr-o interpolare polinomialã în raport cu valorile necunoscutei în fiecare din nodurile elementului respective.

În lucrarea de față s-a realizat un studiu privind analiza cu elemente finite a unui model neliniar de câmp electromagnetic exprimat în potențialul magnetic vector, cuplat cu un câmp termic în sisteme electrotemice ce utilizează inductoare cu geometrie variabilă [11]. Analiza s-a realizat prin simularea fenomenului de încălzire inductivă, aplicație caracteristică problemelor cuplate de câmp electromagnetic și termic, utilizând un program de analiză cu elemente finite.

Capitolul 1 al lucrării prezintă aspect fundamentale privind încălzirea prin inducție electromagnetică.

Capitolul 2 al lucrării prezintă formularea problemei de camp electromagnetic utilizată în soluționarea problemei de curenți turbionari.

Capitolul 3 prezintă modul de formulare și rezolvare numerică a problemei de câmp termic prin metoda elementului finit

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

Capitolul 5 prezintă modelarea numerică a fenomenului de încălzire prin inductie electromagnetică a unor semifabrica utilizând inductoare cu geometrie variabilă.

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

incĂlzirii prin inducție electromagnetică

1.1 Aspecte fundamentale

Pătrundere câmpului electromagnetic variabil în timp în materiale conductoare determină încălzirea acestora prin intermediul curenților turbionari, determinați de tensiunile electromotoare induse.

Orice corp conducător de electricitate se încalzește prin efect Joule atunci când este parcurs de un curent electric. Încălzirea se obține fie aplicând la extremitățile conductorului o diferență de potențial – încălzire clasică prin rezistență electrică – fie amplasând acest conductor într-un câmp magnetic variabil în timp – incălzirea prin inducție electromagnetică.

În cazul unei bobine (figura 2.1.a) parcursă de curent alternativ, acesta va crea în interiorul ca și în exteriorul bobinei un câmp magnetic variabil. Dacă în interiorul bobinei se introduce un corp C din material conductor, fluxul magnetic variabil în timp care parcurge materialul va induce o tensiune electromotoare ce determină apariția curenților turbionari. Dacă fluxul electromagnetic inductor 1 este alternativ, de pulsație , în piesă apar curenți induși I2, de densitate J2 (figura 2.1.b) al căror sens este astfel încât fluxul lor , se opune fluxului inductor.

Câmpul magnetic determinat de curenții de conducție conduce la reducerea câmpului magnetic rezultat în interiorul corpului. Reducerea câmpului magnetic rezultă în interiorul corpului. Reducerea câmpului magnetic rezultant este cu atât mai pronunțată cu cât frecvența este mai mare. În axul corpului, câmpul magnetic va avea valoarea cea mai mica și va creste treptat spre exterior.

Curenții electrici turbionari induși produc căldura prin efect Joule, corpul conductor în care aceștia au luat naștere se va încalzi. Bobina constituie circuitul primar sau inductor, iar materialul conductor constituie circuitul secundar sau Indus. În acest fel echipamentele încalzite prin inducție electromagnetica, sunt, ca principiu de funcționare, cu particularitațile specifice de conducție, asimilabile cu transformatoarele cu sau fara circuit magnetic, primarul fiind alimentat la frecvența industriala medie sau înaltă.

Aplicațiile încălzirii prin inducție nu se limitează la cazul cel mai frecvent al unui corp conductor plasat în interiorul unui solenoid, fiind utilizate configurații foarte variate de inductoare ( inductoare plane, inductoare liniare, inductoare de tip tunnel etc.) și poziții cariate ale corupului de încalzit în raport cu inductorul.

Avantajul solenoidului este acela că valoarea câmpului magnetic ca și inducție magnetică în interiorul acestuia pot avea valori mari. În alte configurații se face apel la circuite care intensifica câmpul magnetic.

Spre deosebire de cazul încălzirilor clasice, unde pentru încălzire se utilizează o sursa de caldură cu temperatura mult mai ridicată, în cazul încălzirii prin inducție electromagnetică căldura este produsă cu ajutorul unui inductor ce ramâne rece în raport cu corpul de încălzit, care este adus la temperaturi foarte mari datorită curenților turbionari induși în acesta.

Calculul simplificat al puterii transmise într-o piesă cilindrică

Creșterea puterii disipate în material procesat în câmpul electromagnetic produs de instalația de încălzire prin inducție poate fi realizată prin creșterea intensității câmpului magnetic și implicit a valorii curentului electric I care parcurge materialul, în sensul câmpului magnetic inductor, fiind numeric egală cu aria suprafeței cuprinsă sub curba J(x) :

I=. (1.2)

Intensitatea curentului electric care parcurge materialul pe o zonă corespunzătoare adâncimii de pătrundere este:

(1.3)

Rezultă astfel că doar 63,20% din curentul electric este concentrat în zona adâncimii de pătrundere.

Puterea dezvoltată în zona adâncimii de pătrundere , este:

, unde P- este puterea totală disipată în semispațiu.

Tipuri de inductoare utilizate în procesele de încălzire prin inducție electromagnetică

CÂMPUL ELECTROMAGNETIC CVASI-STAȚIONAR ÎN STRUCTURI CU MEDII NELINIARE

Pentru formularea corectă a problemei de câmp electromagnetic se impune considerarea unui domeniu conductor C plasat într-un mediu neconductor și nemărginit 0 (aer), în care sunt plasate surse de current (figura 2.1). Dacă presupunem C frontiera domeniului conductor și 0 ca fiind frontiera exterioară a domeniului de calcul, vor fi valabile toate afirmațiile ce vor urma atât pentru corpuri izolante cât și în mișcare. Câmpul electromagnetic în regim cvasi-staționar este descrisă de ecuațiile lui Maxwell, [5], [19], astfel că în domeniul conductor C acesta este descris astfel:

(2.1)

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

(2.3)

(2.4)

(2.5)

unde J reprezintă densitatea curreților turbionari induși.

Figura. 2.1. Configurația domeniului de calcul pentru problemele de curenți turbionari

În 0, vom avea :

(2.6)

(2.7)

unde J0 reprezintă densitatea de current impusă.

În continuare se analizează câmpul electromagnetic în domeniul C, deci se vor rezolva ecuațiile (2.1),…,(2.7). Frontiera C a domeniului C conține și suprafețele ce mărginesc corpurile conductoare. Frontiera domeniului C poate fi în mișcare deoarece și corpurile conductoare pot fi în mișcare. Mișcarea frontierei poate fi arbitrară. Se va considera că toate corpurile conductoare aflate în mișcare se deplasează în interiorul domeniului C . Se presupunem că pe frontiera C vom avem următoarele condiții de frontieră:

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

Ht=f (2.8)

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

Bn=g (2.9)

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

(2.10)

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

(') (2.11')

sau

(") (2.11")

Observații:

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

Bn= 0 (2.12)

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

c. 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

(").

Cunoașterea valorii inducției magnetice: pentru t=0 este o prblemă importantă ce nu trebuie neglijată. În cazul în care regimul este periodic, atunci la ecuațiile (2.1),…,(2.7) și condițiile de frontieră se adaugă condiția de periodicitate, astfel că, condiția inițială nu este necesară. Pentru mediilor în mișcare regimul periodic presupune dificultăți suplimentare de analiză.

Pentru a se asigura unicitatea soluției se impune ca ecuațiile (2.1), …(2.7) împreună cu condițiile de frontieră, condițiile de trecere și condițiile înițiale să definească unic problema de câmp [20]. În acest sens în procesul de modelare numerică sunt utilizate condiții de etalonare, cum ar fi de exemplu condiția de etalonare Coulomb ce impune potențialul magnetic vector A de forma: divA=0.

In caz particular se preferă utilizarea soluțiilor de potențial, în acest sens în continuare find prezentată o metodă de analiză numerică a câmpului electromagnetic cuasistaționar utilizată pentru soluționarea problemei considerate, folosind proceduara (A-V,A), procedura (A-V, A) putând fi aplicată și în cazul în care este multiplu conex. Aceasta este greu de aplicat la structuri cu corpuri în mișcare deoarece rețeaua de discretizare se modifică în timp.

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

1 (2.13)

Atunci din legea inducției rezultă:

2 (2.14)

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

Apoi, din legea conducției:

unde

În consecință avem:

3 (2.15)

0 (2..16)

Din ecuația (2.15) rezultă:

4 (2..17)

Deoarece ecuațiile (2.15) și (2.16) nu asigură unicitatea potențialelor A și V: putem adăuga la A și scădea din gradV gradientul unui potențial cu derivată după normală nulă pe S'. Din (2.13) rezultă că pe componenta tangențială a lui A diferă printr-un gradient. Putem forța conservarea componentei tangențiale a lui E impunând 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.18)

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

Apoi, din (2.15) avem:

Pentru a asigura unicitatea soluției (A,V) trebuie impusă o condiție suplimentară pentru potențialul vector A, numită condiție de etalonare. În plus, pe frontiera C (C 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.19)

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

7 (2.20)

În plus, din (2.17) 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: Potențialele A,V care verifică ecuațiile (2.15), (2.16), (2.19) și condițiile de frontieră (FRA), (2.20) sunt unic determinate în C.

Demonstrație: Pentru rezolvarea numerică a problemei se folosește tehnica Galerkin, proiectându-se ecuațiile (2.15), (2.16), (2.19) 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 va scrie:

,

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

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

8 (2.21)

Integrând prin părți primul termen rezultă:

9

(2.22)

Se proiectează ecuațiile (2.17) și (2.18) pe funcțiile Vk și se obține:

10 (2.23)

Integrând prin părți primul termen, rezultă:

11 (2.23')

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

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

12 (2.24)

Integrând prin părți rezultă:

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

Sistemul (2.21), (2.23), (2.24) 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.21), (2.23) de nE = nN+nV ecuații cu nE necunoscute.

Se poate plasa condiția de etalonare în ecuația (2.15) [9], [10], transformând-o în:

13 (2.25)

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

14 (2.26)

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

15 (2.27)

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

16

(2.28)

Integrând prin părți primii doi termeni, avem:

17

(2.29)

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.23). 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.28), (2.23) este mai agreabil decât (2.21), (2.23), (2.24). Dacă adoptăm elementele de muchie atunci condiția de etalonare este implicită și rezultă sistemul format din ecuațiile (2.21), (2.23).

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

18 (2.30)

ecuația (2.15) devine:

19 (2.31)

În aceast caz apar însă dificultăți atunci când mediul conductor nu este omogen electric.

3. SOLUȚIONAREA NUMERICĂ A CÂMPULUI TERMIC

3.1. Difuzia câmpului termic

Se consideră domeniul cu bordura . Domeniul este ocupat de dielectricul cu pierderi și este inclus în domeniul de calcul al câmpului electromagnetic . Ecuația Fourier pentru regimul staționar al câmpului termic este[3]:

-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ă. Condiția de frontieră este:

– (3.2)

unde este coeficientul de transfer termic la suprafață, iar este temperatura în exteriorul domeniului . De multe ori se poate lua = 0 și în acest caz se determină supra-temperatura față de mediul ambiant. Dacă în relația (3.2) avem = 0, obținem condiții Neumann omogene, iar dacă = 0, rezultă condiția de frontieră Dirichlet.

Teorema 3.1. Ecuația (3.1) are soluție unică în condiția de frontieră (3.2). Demonstrație. Într-adevăr dacă două câmpuri termice verifică ecuația (3.1) și condiția de frontieră (3.2), atunci câmpul diferență 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:

(3.6)

unde c este capacitatea termică volumică [12]. La ecuația (3.6) se adaugă condițiile de frontieră (3.2) și condiția inițială pentru temperatură: .

Teorema 3.2. Ecuația (3.6) are soluție unică în condiția de frontieră (3.2) și cu condiție inițială impusă.

Demonstrație. Este valabilă relația (3.3), dar în locul relației (3.4) avem:

(3.7)

dacă se ține cont că temperatura diferență verifică ecuația (3.6). Din (3.3) și (3.7) rezultă prin integrare:

de unde Td = 0.

3.2. Metoda elementelor finite

Se presupune că frontiera a domeniului este formată din 3 suprafețe , , , unde sunt îndeplinite condițiile de frontieră Dirichlet: , Neumann: și, respectiv, mixte -. Primele 2 condiții de frontieră pot fi considerate cazuri particulare ale condițiilor mixte, atunci când = 0 și respectiv = 0 (Subcap 3.1.). Din punct de vedere tehnic, condiția de frontieră Neumann apare pe suprafețe de simetrie (suprafețe de câmp), fiind întotdeauna nulă. Spre deosebire de cazul câmpului electromagnetic, în cazul câmpului termic staționar se va prezenta tehnica Ritz de soluționare numerică a ecuației (3.1).

Fie funcționala:

(3.8)

al cărei minim se caută în mulțimea funcțiilor care au condiția de frontieră Dirichlet impusă.

Teorema 3.3. Soluția ecuației (3.1), care verifică condițiile de frontieră, minimizează funcționala (3.8) și, reciproc, minimul funcționalei (3.8) este dat de soluția (slabă) a ecuației (3.1) care verifică și condițiile de frontieră Neumann și mixte.

Demonstrație. Fie minimul funcționalei (3.8) și fie o temperatură oarecare, dar care îndeplinește condițiile de frontieră Dirichlet, cu și arbitrari, având pe . Înlocuind în (3.8), se obține:

(3.9)

Se integrează prin părți ultimul termen:

(3.10)

Înlocuind în (3.16) și grupând după puterile lui , se obține:

+>0 (3.11)

Dacă minimizează funcționala , atunci inegalitatea (3.11) 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 rezultă condiția de frontieră (3.2), din al doilea, condiția Neumann iar din al treilea, ecuația (3.1). Reciproc, dacă cele două condiții de frontieră și ecuația (3.1) sunt îndeplinite, atunci termenul în este nul, inegalitatea (3.11) este îndeplinită și este minimul funcționalei .

Teorema 3.3. Funcționala (3.8) este strict convexă:

, (3.12)

Demonstrație. Folosind expresia (3.8), rezultă, prin câteva calcule simple:

Teorema 3.4 arată că funcționala (3.8) are un singur minim local. În plus, rezultă consistența metodei elementelor finite: alegem un șir de subspații finit dimensionale pentru T: … unde W este spațiul funcțiilor definite pe , cu valori reale și care îndeplinesc condiția de frontieră Dirichlet. Obținem valorile , care minimizează (3.8) și pentru care:

>>>…>

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

(3.22)

unde sunt funcții liniar independente cu condiții de frontieră Dirichlet nule, numite funcții de formă, iar o funcție arbitrară care îndeplineste condiția de frontieră Dirichlet. Se proiectează ecuația (3.1) pe funcțiile liniar independente , cu condiții de frontieră Dirichlet nule, numite funcții test:

(3.23)

Se integrează prin părți membrul stâng:

(3.24)

Dacă se ține cont de condițiile de frontieră, se obține:

(3.25)

Inlocuind (3.25) în (3.24) și apoi în (3.23), rezultă:

(3.26)

Ținând cont de (3.22), rezultă sistemul de ecuații:

+

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

Dacă se aleg funcțiile test aceleași cu funcțiile de formă, atunci ecuația (3.27) este identică cu (3.18).

Se recomandă utilizarea elementelor nodale.

Difuzia termică

Difuzia termică este descrisă de ecuația (3.6):

– div grad T + c (3.28)

Operatorul de difuzie nu se mai bucură de proprietatea de a fi pozitiv definit și simetric. Din acest motiv, rezolvarea numerică ecuației (3.28) se face prin procedura Galerkin. Temperatura se scrie sub forma:

(3.29)

Spre deosebire de expresia (3.22), în relația (3.29), coeficienții sunt funcții de timp. Ecuația (3.28) capătă forma:

++

++

+, j = 1, 2, … , n (3.30)

3.3. Discretizarea în domeniul timp

Sistemul (3.30) se mai poate scrie:

(3.31)

unde

=++

Soluționarea numerică a ecuației (3.31) se poate face prin metoda Crank-Nicholson [1]. Se integrează (3.31) pe intervalul de timp , de lungime :

(3.32)

Se aproximează integrala în domeniul timp cu:

= (3.33)

unde este un coeficient. În particular, avem:

pentru = 0, metoda Euler explicită;

pentru = 1, metoda Euler implicită;

pentru = 1/2, metoda trapezelor.

Folosind (3.33), relația (3.32) devine:

(3.34)

unde:

si

Relația (3.34) reprezintă un sistem de ecuații cu necunoscutele . Pasul de timp se corectează astfel încât norma variației de temperatură:

=

să se păstreze între niște limite impuse. Dacă la pasul de timp se obține >, atunci se înjumătățește următorul pas de timp: =, iar dacă <, atunci următorul pas de timp se dublează =.

4. CUPLAREA PROBLEI DE CÂMP ELECTROMAGNETIC CU CEA DE DIFUZIE TERMICĂ

Abordarea numerică a problemei complexe de camp electromagnetic caracterizată de ecuații cu derivate parțiale presupune parcurgea următoarelor etape succesive:

Descrierea problemei care presupune realizarea geometrie, stabilirea proprietăților fizice, realizarea rețelei de discretizare și asocierea caracteristicilor de material

Dezvoltarea metodei numerice de simulare (MDF, MEF);

Verificarea, vizualizarea și interpretarea rezultatelor simulării.

Etapele menționate mai sus corespund corespund structurii naturale a software-ului în trei procesoare:

– procesorul de introducere a datelor (preprocesorul);

– procesorul de calcul;

– procesorul de extragere a rezultatelor (postprocesorul).

Descrierea geometriei, domeniului de calcul și discretizarea cu elemente finite este asigurată de către preprocesorul geometric al pachetului software de analiză. Discretizarea domeniului de calcul reprezintă trecerea de la mediul continuu la un mediu discret, ea fiind strâns legată de metoda utilizată pentru aproximarea ecuațiilor cu derivate parțiale.

În metoda elementului finit, discretizarea domeniului divizarea domeniului de calcul într-un ansamblu de subdomenii – elemente – cu respectarea frontierelor și a suprafețelor de trecere ale domeniului initial, rezultând astfel numărul de noduri. Relativa ușurință a discretizării în elemente finite și marea generalitate a procedurilor numerice asociate fac ca această metodă să fie foarte utilizată în softurile CAD [17].

Pentru discretizarea cât mai bună a modelului considerat, trebuie să fie îndeplinite următoarele condiții:

– densitatea elementelor în zonele unde fenomenul este mai intens trebuie să fie cât mai mare;

– elementele trebuie să fie suficient de regulate (dacă avem triunghiuri, ele trebuie să fie cât mai aproape de echilateral);

Utilizatorul definește densitatea elementelor propunând discretizarea frontierelor domeniului, în timp ce sarcina divizorului automat (cu care este prevăzut programul ) este aceea de a propaga automat divizarea către interiorul domeniului, respectând densitățile alese.

Descrierea proprietăților fizice presupune două etape distincte:

– identificarea diferitelor regiuni și a diferitelor porțiuni ale frontierei domeniului de calcul;

– specificarea caracteristicilor fizice, adică descrierea efectivă a caracteristicilor de material a surselor câmpului și a tipului de condiții la limită pe fiecare frontieră. Legăturile între regiuni și caracteristicile fizice se fac prin intermediul numelor asociate regiunilor și frontierelor în faza de descriere topologică.

Rezolvarea problemei de curenți turbionari presupune luarea în calcul a parametrilor de material (rezistivitatea electrică și dependentă B-H) care sunt direct dependenți de temperatură, spre deosebire de rezolvarea problemei de câmp termic, în cazul căreia parametrii de material vor depinde de rezultatele obținute în urma rezolvării problemei de curenți turbionari, motiv pentru care, la fiecare pas de timp adoptat pentru problema termică, se impune revenirea la problema de curenți turbionari și la problema de difuzie, corectându-se parametrii de material. Când corecția nu este semnificativă, se avansează la următorul pas de timp [12]. Dacă se semnalează instabilitate în timp, atunci pasul de timp se micșorează. Tratarea neliniarității relației constitutive B-H se face iterativ, prin modelul pseudo-liniar, corectând permeabilitatea magnetică în funcție de valoarea efectivă a inducției magnetice.

MODELAREA NUMERICĂ A INCĂLZIRII PIRN INDUCȚIE A SEMIFABRICATELOR UTILIZÂND INDUCTOARE CU GEOMETRIE VARIABILĂ

Dezvoltările tehnologice din domeniu impun tot mai des utilizarea unur subansamble având o formă geometrică mai complexă, respectiv utilizarea unor materiale având caracteristici tehnice superioare. Pornind de la acest deziderat, în vederea raspunderii solicitărilor primite, se încearcă îmbunătățirea proprietăților de material prin aplicarea de tratamente termice. Datorită avantajelor ce le prezintă, incălzirea prin inducție electromagnetică, în cdrul proceselor electrotermice inductive se încearca adaptarea echipamentelor inductive funcție de tipul aplicației. Una din cele mai utilizate metode consta în adaptarea formei inductorului la forma geometrică a piesei de încălzit, astfel încât să asigure o distribuție cât mai uniformă a câmpului electromagnetic și implicit a câmpului termic în zona de interes. Din dorința de a mări randamentul de conversie electroenergetic la nivelul ansamblului inductor piesă, în cadrul modelului considerat s-a incercat adoptarea unei soluții care a presupus utilizarea unui inductor având un număr de două spire cu diametru diferit și secțiune rectangulară.

Constante de material

Se presupune cunoscută frecvența de lucru a instalației ca fiind 15 kHz, iar densitatea de curent prin inductor ca fiind 21 [A/mm2], mărimi caracteristice similare unei instalații de călire prin inducție produsă de firma LEPEL CORP USA [21].

Materialul utilizat în cazul ambelor repere este identic având următoareale proprietăți:

Materialul din care este realizat semifabricatul este oțel AISI-SAE-4135, având 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 [55]:

[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ă în Figura 5.3. Poate fi utilizată [4] relația:

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ă.

– Conductivitatea termică izotropă  [W/m/C], a carei dependență de temperatură  este dată în Figura. 5.6. Următoarele funcții [4] 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 [4]:

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.7.

Temperatura de austenitizare este θa = 770 °C.

5.1.Definirea problemei

5.1.1 Inductor cu două straturi de secțiune rectangulară utilizat pentru in încălzirea prin inducție a suprafetelor cindrice cu exterioare

Se considera un semifabricatul cilindric cu geometrie complexă având caracteristicile de mai sus, plasat într-o instalație de încălzire prin inducție electromagnetică dotat cu un inductor format din două spire având diametrul diferit pentru cele două spire, astfel încât acesta să urmărească forma geometrică a suprafeței de încălzit. Inductorul 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ă. Din motive de simetrie axială, se consideră că este suficient să se realizeze analiza numerică doar pentru un sfert din domeniul de calcul. Semifabricatul considerat, având geometrie complexă este reprezentat de un reper cu destinație specială în industria auto, reprezentată în Figura 5.1. Pentru o mai buna distribuție a câmpului termic în semifabricat se va considera inductorul cu secțiune rectangulara rectangulară cu dimensiunile laturilor de 10×10 [mm], o grosime a peretelui de 2,5 [mm] și distanță între spire de 6 [mm], alimentate de la o sursă cu frecvența tensiunii de alimentare de 15 (KHz)

Fig 5.1. Geometria modelului

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

Domeniul de calcul pentru câmpul termic este cel al piesei astfel că 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ă.

În urma simulării numerice în vederea determinării temperaturii atinse în semifabricat s-a constattat că este nevoie de 12 secunde pentru ca semifabricatul sa atinga temperatura dorita de 930 grade C.

Rezultattele obținute în urma simulării numerice scot în evidență valorile și distribuția mărimilor de câmp. Figurile 5.2-5.5 prezintă distrinuțiile fluxului magnetic pentru diferite momente de timp în ansamblul piesă inductor pentru durata tde timp necesară atingerii temperaturii dorite. Pentru o mai bună întelegere a fenomenelor electromagnetice, în figurile 5.6.-5.9 se prezintă distribuția câmpului magnetic pentru ansamblul inductor-piesă.

vor prezenta rezultattele simulărilor în vederea evidențierii distribuție mărimilor de camp electromagnetic si termic pentru o mai bună înțelegere a fenomenelor ce apar pe durata încălzirii.

Constante de material și condiții de frontieră utilizate în problema de câmp electromagnetic

Materialul din care este realizat semifabricatul este oțel AISI-SAE-4135, având 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 [55]:

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

pt. [20, 800]C,

și

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

pt. [800, 1400]C

Relatia B-H, a cărei dependență de temperatură este descrisă în Figura 5.3. Poate fi utilizată [4] relația:

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ă.

Condițiile de periodicitate și sensul curentului impun utilizarea următoarelor condiții de frontieră: Dirichlet nulă pentru potențialul vector și axa de simetrie așa cum se poate observa în Figura 5.4

Domeniul de calcul pentru câmpul termic este cel al piesei astfel că 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.8).

Constante de material și condiții de frontieră utilizate în problema termică

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

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

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

și

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

Figura 5.4. Condiție de frontieră pentru problema de câmp electromagnetic

Domeniul de calcul pentru câmpul termic este cel al axului. 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ă.

Resultatele obținute prin simulare numerică scot în evidență influența formei geometrice a inductorului asupra distribuției câmpului electromagnetic, termic și timpului necesar pentru atingerea temperaturii dorite în procesul de încălzire inductivă așa cum se poate onbserva din rezul.

Rezultatele simulării prezentate în figurile 5.9 și 5.10 scot în evidență infuența frecvenței asupra densității de current indusă și densității de putere indusă în semifabricat pe durata procesului de încălzire inductivă. Deasemenea în figurile 5.11, 5.12 și 5.13 sunt prezintă influența frecvenței asupra distribuției câmpului termic și mișcare conform pașilor menționați în tabelul 1. Se observă că, odată cu creșterea frecvenței se obține reducerea semnificativă a timpului de încălzire prin inducție. Deasemenea distribuția câmpului termic în interiorul semifabricatului este similar pentru cele trei cazuri studiate.

Figura 5.13 prezintă distribuția câmpului termic pentru fiecare pas de timp prezentat în tabelul nr 1 în cel de-al treilea caz de analiză considerat, când frecvența tensiunii de alimentare a fost de 10 KHz

CONCLUZII.

Simularea prin metoda elementului finit da posibilitatea calculului mărimilor specifice câmpului magnetic, posibilitatea integrării pe o curbă, pe o suprafață sau un volum bine definit, care este mult mai avantajoasă față de metoda clasică, de proiectare, realizare și încercare a prototipurilor, fiind și mult mai ieftină, oferind o bogăție de informații concludente pentru utilizator.

Pentru analiza geometrică ale celor trei configurații am folosit programul de analiză MAXWELL 2D

Analizad cele doua metode se observa ca cea care da mai multe informatii este metoda de calcul analitica care are urmatoarele caracteristici:

este mai complexa ,

necesita un nivel de cunostinte mai ridicat pentru a rezolva ecuatiile matematice

necesita un timp de efort mai mare pentru obtinerea rezultatului

timpul pentru obtinerea rezultatului este considerabil si creste odata cu complexitatea modelului de analiza adoptat

este o metoda care poate sa obtina rezultate fine si precise, reflectand un rezultat cat mai exact fata de conditiile reale

Metoda de calcul numerica are urmatoarele caracteristici:

este o metoda de calcul simplificata

necesita un timp mai scurt pentru obtinerea rezultatelor

resursele folosite sunt mult diminuate

datele obtinute prin aceasta metoda reflecta rezultate obtinute in conditii particulare in compartie cu alte date obtinute in conditii analitic generalizate

Metoda de calcul folosind programul informatic MAXWELL are urmatoarele caracteristici:

ofera posibilitatea de a modela modelul matematic si fizic privind analiza si calculul complex pornind de la un model simplificat la un model analitic complex

are limite privind implementarea modelului matematic, fiind un program informatic dezvoltat pe o anumita structura nu permite abstractizarea metodei de calcul, altfel spus are la baza sabloane predefinite privind realizarea calculelor

este o solutie rapida prin care se pot obtine rezultate reale, se pot face simulari in timp scurt, fapt ce poate contribui la luarea unei decizii rapide privind solutia adoptata

pune la dispozitie reprezentari grafice a modelului realizat si a marimilor calculate

necesita o dotare privind tehnica de calcul, instruirea personalului, licentierea aplicatiilor informatice, toate acestea insemnand costuri.

Concluzionand putem face o comparatie intre metodele prezentate anterior dupa cum urmeaza:

gradul de detaliere si precizie a rezultatelor

timpul de obtinere a rezultatelor

costurile de obtinere a rezultatelor

Metota analitica ne da cea mai clara imagine asupra rezultatului, dar implica un efort uman (nivel de cunostinte si timp) mai ridicat, posibil ca in anumite cazuri sa nu fie eficient din punct de vedere al raportului rezultate/timp.

Metoda numerica este cea care poate obtine un rezultat bun in cadrul proiectelor in care nu se doreste a se investi multi bani sau timp, obtinand rezultatele dorite.

Metoda de calcul folosind programe informatice este cea care in conditiile in care se doreste obtinerea in timp scurt a unor rezultate foarte bune este cea mai buna solutie privind raportul rezultate/timp, dar are dezavantajul ca metodele de calcul nu pot fi abstractizate la acelasi nivel ca si in metoda analitica.

Astfel in functie de rezultatul dorit privind rapoartele rezultate/timp/costuri se poate opta in functie de caz pentru oricare din metodele prezentate anterior.

Deci fiecare metoda folosita, are si avantajele ei economice ,in functie de timpul necesar pentru a obtine un rezultat cat mai precis, de bani investiti pentru a obtine rezultate cat mai precise si cat mai reale, de implicarea factorului uman in obtinerea acestor rezultate.

Alegerea celei mai viabile metode se face in functie de asteptarile pe care le au cei care decid acest lucru si in functie de interesele economice ale fiecarei companii, dar tinandu- se cont de :rezultate,timp si costuri.

Calculul problemei de câmp electromagnetic cuplat cu cel termic utilizând metoda aleasă permite realizarea analizei numerice pentru problema considerată fără modificarea rețelei de discretizare, datele obținute pentru fiecare etapă fiind considerate ca date inițiale pentru etapa următoare sporind astfel precizia.

În lucrare se tratează o variantă eficientă pentru aplicațiile industriale: se aleg pozițiile pe care le ocupă porțiunea activă a bobinei fictive (în lucrare sunt echidistante). La fiecare poziție se “încălzește” piesa până ajunge la temperatura dorită, apoi se trece la poziția următoare. Rezultă deci timpul petrecut de piesă în fiecare poziție și, în consecință, viteza de deplasare. Se poate determina astfel viteza de deplasare optimă a piesei în raport cu inductorul pentru a se obține parametrii termici indicați pentru procesul de tratare termică de suprafață prin inducție, parametri fundamentali pe care-i dorește și tehnologul pentru realizarea procesului tehnologic.

Importanța economică a problemei studiate rezultă din utilizarea în industrie a tehnologiilor de tratare termică de suprafață a axelor, bolțurilor și, în general, a pieselor care au cel puțin o dimensiune mare, astfel încât nu permit acoperirea întregii suprafețe ce urmează a fi tratată termic de bobina inductoare. Pentru tehnolog este deosebit de util să poată stabili evoluția în timp a vitezei de deplasare a piesei, și, eventual, a intensității curentului din bobină, astfel încât să realizeze un tratament optim pentru suprafața de interes. Cunoașterea acestei evoluții permite automatizarea eficientă a proceselor specifice tratamentelor termice.

BIBLIOGRAFIE

[1] T. Leuca, Câmpul electromagnetic și termic cuplat – Curenți turbionari, Editura Mediamira, Cluj-Napoca 1996

[2] V. Firețeanu, T. Leuca, Inducția electromagnetică și tehnologii specifice, Ed. Mediamira, Cluj-Napoca, 1997

[3] www.ansys.com

[4] Maxwell 2D User Guide

[5] F. I. Hănțilă, E. Demeter, Rezolvarea numerică a problemelor de câmp electromagnetic, Ed. Ari Press, ICPE-ME, 1995

[6] F.I. Hănțilă, t. Leuca, C. Ifrim, electrotehnică teoretică, Ed. Electra București, 2002

[7] R.Pascal, P.Conraux, J.M.Bergheau Coupling Between Finite Elements and Boundary Elements for the Numerical Simulation of Induction Heating Process Using a Harmonic Balance Method T-MAG may 03 1535-1538

[8] G. Szymanski, M. Waszak Calculation of Eddy Current Energy Losses in Thin Sheets under Saturation Proceedings of COMPUMAG’04 Saratoga Springs, USA, 2003

[9] I. Sebestyén, S. Gyimóthy, J. Pávó, O. Bíró Calculation of losses in laminated ferromagnetic materials T-MAG march 04 924-927

[10] J. Zou, Y.Q. Xie, J.S. Yuan, X.S. Ma. Analysis of the thin plate eddy current problem by Finite Volume Method T-MAG march 04 1370-1373

[11] K. Yamazaki, S. Watari, A. Egawa. Adaptive Finite Element Meshing for Eddy Current Analysis of Moving Conductor T-MAG march 04 993-9

[12] M.Arion, T. Leuca, F. I. Hantila, Numerical analysis method for solving the coupled electromagnetic and thermal field Questions for induction heating systems with moving parts, Journal of Optoelectronics and Adavanced Materials 2008, Vol. 10, No.5, ISSN 1454-4164, p 1213-1217

[13] T Leuca, M. Arion, About Numerical Analysis of Electromagnetic and Thermal Field in Induction Equipment with Moving Bodies, Revue roumaine des sciences techniques
SérieÉlectrotechnique et Énergétique, No: 3, 2009, Publisher : Romanian Academy, Publishing House of the Romanian Academy,Format print and electronic ISSN: 0035-4066, p 271-280

[14] T. Leuca, Teoria câmpului electromagnetic. Aplicații utilizând tehnici informatice. Editura Universității din Oradea, pg 267, 2001, (ISBN 973-613-071-1).

[15] T. Leuca, F.I. Hănțilă, Livia Bandici, Carmen Molnar, Bazele electrotehnicii, Editura Mediamira, Cluj-Napoca, 211 pg., 2007, (ISBN 978-973-713-189-8).

[16]. Șt. Nagy, T. Leuca, A. Nagy – Conceptul CAD/CAM. Teorie și aplicații, Editura Universității din Oradea, 164 pg, 2006, ISBN (10) 973-759-054-6, ISBN (13) 978-973-759-054-1.

[17]. Șt. Nagy, T. Leuca, T. Maghiar – Procesarea materialelor în câmp electromagnetic. Aplicații utilizând tehnici informatice, Editura Universității din Oradea, 155 pg., 2002, (ISBN 973-613-120-3).

[18] 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

[19] F. I. Hănțilă, A method for solving 3-D eddy current problems in non-linear media, Rev. Roum. Sci. Techn.-Electrotechn. et Energ., no.3, 1992, p. 267-281

[20] F.I. Hănțilă, Existences and stability theorems of the stationary field without sources and with null boundarz conditions in nonlinear media, Revue Roum. Sci. Techn. Ser. Electrotechnique et Energ., no.3, 1981