Reševanje stacionarnega in nestacionarnega temperaturnega polja po metodi končnih elementov na PC računalnikih
Marjan Bolčina*1	UDK: 536.241:669.046:621.78.012.5:681.3.06
ASM/SLA: J2g, D6, U4g, U4k
S hitrim razvojem računalniške tehnologije se v inženirski praksi vse pogosteje uporabljajo numerične metode reševanja konkretnih problemov reševanja polj, opisanih z diferencialnimi enačbami. Danes je takšen način projektiranja ekonomsko in kakovostno upravičen in postaja nujen.
Naša naloga je reševanje nestacionarnih temperaturnih polj poljubnih objektov in izdelava ustreznega računalniškega programa. Pri tem imamo lahko definirane kot robne pogoje vsiljene temperature ali toplotne tokove ali notranje izvore toplotnega toka, ki jih povzroča električni tok. Tako lahko projektiramo (glede na temperaturno polje) objekte, kot na primer: hladilna telesa, segrevanje kovinskih in izolacijskih delov, temperaturno obremenljivost el. zbiralk, indukcijsko kaljenje in indukcijsko taljenje kovin.
UVOD
V nadaljevanju torej iščemo rešitev dvodimenzionalnega ali aksisimetričnega trodimenzionalnega, časovno odvisnega temperaturnega polja, ki ga na objektu poljubne geometrije poleg materiala in okolice povzročajo vsiljene temperature, toplotni tokovi in toplotni izvori (električni tokovi).
Ker ne obstaja splošna analitična metoda za reševanje problemov takšne vrste, moramo iskati rešitve z nu-meričnimi metodami ali s simulacijo konkretnega primera na analognem računalniku. Ker so ustrezni analogni računalniki zelo dragi, praktično vedno uporabljamo numerične metode, prirejene za digitalne računalnike. Numerične metode v splošnem ločimo na metodo končnih diferenc in metodo končnih elementov. Izbira metode se vse močneje nagiba v korist metode končnih elementov.
Pri uporabi metode končnih elementov objekt poljubne oblike razdelimo na poljubne končne elemente znane geometrije. Na splošno velja, da večje število končnih elementov (manjši, kot so), poveča natančnost iskane rešitve sistema. Natančnost rešitve pa je odvisna tudi od predvidene funkcije polja znotraj elementa. Izbrali smo trikotno obliko končnega elementa, ker z njo najlažje opišemo poljubno obliko obravnavanega objekta. Zaradi znanih prednosti naravnih (area) koordinat pred Katezijevimi (lažje integriranje, singularnost matrik) smo izvedli ustrezne transformacije.
STACIONARNO STANJE
S stališča stacionarnega prenosa toplote nas je zanimala temperaturna porazdelitev po posameznih elemen-
" Marjan Bolčina, dipl. ing. — Slovenske železarne — Železarna Štore
Rokopis prejet november 1990 Originalno publicirano: žzb 25 (1991) 1
tih indukcijske lončne peči, še posebej pa temperaturni gradient v ovoju induktorja. Induktor se namreč greje zaradi prevajanja izmeničnega električnega toka, ki zaradi skin efekta ni enakomerno porazdeljen po bakrenem profilu, pa tudi zaradi prevajanja toplote taline skozi ognjeodporno obzidavo.
NESTACIONARNO STANJE
Število sarž, ki jih lahko stopimo z enkrat pripravljeno ognjeodporno obzidavo, je odvisno tudi od temperature, ki jo doseže talina na meji z obzidavo. To področje prenosa toplote je v časovnem prostoru relativno neraziskano, zato smo si zadali nalogo, da to problematiko natančneje obdelamo. Ker se indukcijska lončna peč v prvi aproksimaciji obnaša kot integrator vhodne moči, moramo omenjeno problematiko reševati v nestacionarnem časovnem področju.
Kot bomo videli kasneje, se največ električne energije v vložku pretvori v toploto prav ob obzidavi. Ta je direktno odvisna od moči peči in frekvence napajalne električne napetosti. Zato lahko prihaja na teh mestih do neželjenih nadtemperatur, ki presegajo dovoljeno temperaturo neke ognjeodporne obzidave, medtem ko vložek na sredini še ni dosegel niti temperature taljenja. (Dobljene rezultate pa bi lahko uporabili tudi na področju kaljenja ali toplotnega čiščenja).
Pri reševanju električnega dela indukcijske lončne peči smo pristopili k izvornemu načinu. Rezultat analitične analize teoretičnega izhodišča, podanega s Fredhol-movo integralsko enačbo II. vrste
Es (Je, Jm) - 1/Je = E0 (J0)
Hs (Je, JJ-1/Jm = H0 (J0)
je specialni programski paket IMF-CAD, ki iz geometrijskih in snovnih parametrov konstruirane peči izračuna električne veličine (el. impedanco, el. napetosti, el. tokove itd) v posameznih točkah indukcijske lončne peči (7). Pri tem je za numerični izračun uporabljena metoda končnih elementov. Rezultati so podani numerično (tabelarično) in grafično. Izhodni podatki iz tega programa za dano peč (el. tokovi) so tudi vhodni podatki za analizo prenosa toplote, kar je cilj naše naloge. Podrobnejši opis tega programa presega namen tega dela.
REŠEVANJE DIFERENCIALNE ENAČBE PRENOSA
TOPLOTE PO METODI KONČNIH ELEMENTOV
Računalniški program je napisan na osnovi naslednjih matematičnih izhodišč: Fourierjeva diferencialna enačba širjenja toplote (1), vsiljene z induciranim električnim tokom (npr.: v indukcijski lončni peči) je:
div (k gradT) + pc 5T/St = aJ2 = 8q/8V T = T (x, y, z, t)
Za prakso so pomembni tudi naslednji robni pogoji, ki jih opišemo z obrazcem:
k 8T/8n = h (T-TJ
—	predpisana temperatura ali toplotni fluks na delu površine
—	predpisana konvekcija (radiacija)
GALERKINOV POSTOPEK
Predvidimo, da so dane diferencialne enačbe L (t) =f in R (t) = g, pri čemer sta L (t) in R (t) diferencialna operatorja, tako da je enačba polja temperature zadovoljena po površini in na robovih objekta.
L (t) = 8/8x| (k 8T/8X|) + Q = pc (8T/8t) = f R (T)--k (8T/5n) = h (T-T0) = g
Integracijska enačba Galerkinovega postopka (2, 3) se tako glasi:
J (L (t) — f}■ u■ dA + j(R (t) — g}-u-dS = 0
a	s
Funkcijo T izrazimo s polinomom, ki v polju zadošča diferencialni enačbi, na robovih pa zavzame predpisane robne vrednosti, in izenačimo njeno variacijo 81 s funkcijo u:
T=blT1 + b2T2+b3T3 + . . , + bnTn u = 3T
{{L(r)-f)-3Tk.dA+/{R (t) - g) • 8Tk ■ dS = 0 A	s	(k= 1, 2, 3 ... n)
Zgornjo enačbo lahko direktno uporabimo v metodi končnih elementov tako, da jo formiramo za vsak končni element. Polinom zapišemo zato v matrični obliki T = f b] {T}e. Pri tem je [bj interpolacijska matrika, |T)e pa vektor vozliščnih temperatur končnega elementa.
Da si olajšamo diferenciranje v zgornji enačbi in da se izognemo polinomom višjega reda, operator L (t) transformiramo:
L (t) = [8/8x| mi (t)] + L' (t) oziroma L (T) 3T = 8/8X, [M, (t) 8t\--m, (t) [8/8X; (3t)] + L'(t)dT
Po urejevanju in združevanju dobimo naslednjo va-riacijsko integralno enačbo:
f {(Q - pc 8T/8t) 3T + k 8T/8xj [8/8x; (ST)]}
A dA+fh (T-To)3TdS = 0
s
Ko vpeljemo diskretizacijo za končni element temperature in njene variacije v obliki:
T = b. T, oziroma 8J = bjdT,;
lahko dobi prejšnja enačba obliko:
Hi( T + Lij (8T/8t) + F, = 0
H;, = H, + H2 = Jk (Sbi/SK^ (5bj/5xi) dA + jh b^ dS =
a	s
Fj = F, + F2= JQ bj dA + Jh b,Tc dS
a	s
LN = Jpc b,bj dA
_	a
Da bi lahko rešili enačbo končnega elementa in pripadajoče integrale, moramo določiti matriko [b], kot interpolacijsko funkcijo polja v elementu. Nato lahko zapišemo:
Te= [b] (T)
Te = T (x, y, z)
(t} =
T, T2
t3
(vektor vozličnih temperatur)
Kot smo že omenili, smo vzeli trikotno obliko končnega elementa, z vozličnimi točkami 1, 2, 3 (i, j, k). Prav tako smo predvideli, da se temperatura znotraj končnega elementa širi linearno, torej:
Te = c, + c2x + c3y = [a] (c)
[a]	= [1, x, yj
{c) = [c1, c2, c3]T
Konstante v (c) določimo iz robnih pogojev, saj poznamo vrednosti temperatur jT) = [T,, T2, T3j v vozličnih točkah xn, yn
{T} = [akj {c}
[b]	= [a] [ak]~1 = [a] [B]
če [Bj izrazimo s katezijevimi koordinatami:
[B] = 1/2A
x2y3-x3y2, x3y, — x1y3, x,y2-x2y1
y2-y3, y3-yi, yi—y2
X2 — X2, Xt — x3, x2 — xf
= 1/2A[[B1], [B2], [B3]],
pri čemer je A površina trikotnika končnega elementa. Nato računamo integrale matrik [Hj, [L], [F]. Kot smo že omenili, bomo postopek izvedli v brezdimenzijskih (area) koordinatah (I,, l2, l3). Velja zakonitost I,+ l2 + |3= 1. Kate-zijeve in area koordinate so v linearni odvisnosti. Pri računanju integralov si pomagamo z obrazcem:
[F]:
f I?112 l3 dA = 2A (m! n! p!) (m + n + p + 2)!"1
a
Sledi konkretno določanje integralov matrik [Hj, [L], [H,] = k|j[[B]T (8 [a]V8x) (8 [a]/8x) [B]] +
+ [[B]T([5a]V8y)([8a]/8y)[B])) dA Po urejanju sledi končni izraz za H,:
[HJ = k/4 A
y23 +x
yi; + x?3
y23y3i + x32xi3, y23yi2 + x32x21 y3iyi2 + xi3x21
SIMETRIČNO
[H2] = h j'[b]T [b] dS = h f
s	s
In
Yl2 + X2
[l„
dS =
Ker gre za konvekcijo, ki je le na robovih, je l3 = 0. Po integraciji po dolžini dobimo:
[H2] = hls/6
2, 1, 0 1, 2, 0 0, 0, 0
Pri tem sta vozlišči 1 in 2 na začetku in koncu ls (na robu objekta), vozlišče 3 pa je v notranjosti objekta.
[L] = pc J[b]T[b] dA = pcA/12
a
2, 1, 1 1, 2, 1 1, 1, 2
[F,] = Q [ [b]T dA = QA/3
A
[FJ = h f [b]T dS = hls/2
s
dS
Po sumiranju enačbe H^-l-L,, (5T/5t) + F, = 0 po vseh končnih elementih dobimo matrično enačbo celotnega sistema:
[H] (O) + [L] {5<D/8t|+[F] = 0
Tako lahko računamo vektor vozliščnih temperatur jO) za predpisane časovne intervale ali končno stacionarno stanje. V primeru časovno odvisnih problemov lahko uporabimo metodo desnih diferenc /2,4/. Če predvidimo, da se temperatura v malem časovnem intervalu spreminja linearno, lahko zapišemo:
{$}, = |0}t_, + 1/2 ((SO/St),., +{60/8t}t) At
in:
(8<D/5t)t= -{60/8t(,_1 + 2/At ({<!>},-(Oj,.-,) po ureditvi zadnjih treh enačb lahko zapišemo:
|[H]+2/At [L]){0}, = [L] (80/ 8t),_1+2/At [L] {©},_,-[F]
Postopek računanja je naslednji: najprej izračunamo vektor JO}, za prvi časovni interval. Temu sledi določevanje odvoda temperature (8®/St}t na koncu prvega časovnega intervala. Postopek ponavljamo za drugi, tretji in n-ti časovni interval. Seveda moramo podati (pri raču-
nanju prvega intervala) začetne pogoje, in sicer za (O|0 in {SO/8t}0. Začetna temperatura je običajno temperatura okolice, (80/8t)o pa za sisteme višjega reda lahko predvidimo, da je O ali pa ga določimo iz enačb pc (SO/St)e = aJ2=Q za posamezni element.
RAČUNALNIŠKI PROGRAM ZA REŠEVANJE
TEMPERATURNEGA POLJA
Računalniški program smo napisali po prej navedenem Galerkinomev postopku v matrični obliki.
__ Težili smo k temu, da bi bili vhodni podatki sistema in rešitve sistema čimbolj grafično podprti, pri tem pa naj grafika ne bi po nepotrebnem zmanjševala razpoložljivega računalniškega spomina in hitrosti izvajanja programa. Zato smo napisali ustrezen CAD program. Več kot ima računalnik linearnega RAM pomnilnika na razpolago, več elementov lahko vsebuje obravnavani objekt. Tako lahko na Mega2 ATARI ST računalniku (2Mb RAM linearno) rešujemo probleme z več kot 600 elementov v relativno izredno kratkem času (nekaj minut).
Celoten program je zbir procedur (podprogramo-mov), zato je fleksibilen, odprt, objektno orientiran in ga je relativno lahko modificirati. Z pomočjo CAD (GEM) vmesnika pa smo program podprli s prijaznim okoljem za interaktivno delo na PC računalniku, saj vsebuje avtomatski generator topologije sistema in robnih pogojev. Namenjen je predvsem za reševanje Fouriejeve enačbe, kjer je temperatura posledica izvorov toplote, povzročene z električnim tokom. Takšni primeri so el. vodniki, el. zbiralke, uporovni grelci, indukcijsko kaljenje, taljenje in kuhanje.
e<X<e.968(n), y<0.21(n) 221 clc*CTi"tou, 135 uozliSč
<ESC>neu <*>run <INS>edit <CURS0R>»oue <HELP>load,saue
<=KO*FICIEHTI SISTEM*: TX= SS M= 17.5
_____________Tir>*= am
28 &S. 84	112 14« 168 t* £24 ( »0= 1-2
1.9« 2-87 4.S3 V dt9- C*
27 5S 83 111 139 167 I »t «3	mSV= 7M*
ZC	Tre T3S	T54	c= 7*
2.1 2.74 4.47	TKttW M
25 53 8 1	137 16S 193 22 1	«• .oTKVl * J - 834
Rešitev sistema prikazana na grafičnem monitorju
Rešitve obravnavanega sistema po metodi končnih elementov dobimo v tabelarični in grafični obliki. Geometrija objekta se avtomatsko skalira na velikost zaslona, v vozlišča elementov se izpišejo iskane rešitve temperature oz barvna lestvica, v elemente pa povprečne temperature elementov oziroma vrednosti električnih tokovnih gostot in pripadajočih toplotnih tokov. Prav tako dobimo izpis koeficientov sistema, pri katerih smo dobili dano rešitev. Po potrebi lahko spreminjamo samo določene koeficiente in opazujemo tendenco novih rešitev. Na ta način smo dobili dobro orodje za inženirsko projektiranje.
LITERATURA
1.	Holman J. P.: HEAT TRANSFER, McGravv-Hill, New York, 1986
2.	Prelog E.: METODA KONČNIH ELEMENTOV, Fakulteta za arhitekturo, gradbeništvo in geodezijo, Ljubljana, 1975
3.	Zienkiewics O. C.: FINITE ELEMENT METHOD, McGravv-Hill, London, 1975.
4.	Kenneth H. H.: FINITE ELEMENT METHOD FOR ENGINE-ERS, John Wiley & Sons, New York, 1975
5.	DRUGI JUGOSLOVANSKI SIMPOZIJ O METODI KONČNIH ELEMENTOV IN RAČUNALNIŠKEM PROJEKTIRANJU, VTŠ Maribor, Maribor. 1979, 186—210
6.	Baker A. J.: FINITE ELEMENT COMPUTATIONAL FLUID MECHANICS, McGraw-Hill, Washington, 1985
7.	Željeznov M,, Popovič M., ValenčičV., Sinigoj A KorezA.: RAČUNALNIŠKI PROGRAM »IMF-CAD« ZA PROJEKTIRANJE ELEKTRIČNEGA DELA INDUKCIJSKIH LONČNIH PEČI Z RAČUNALNIŠKIM SISTEMOM »ATARI ST«, Fakulteta za elektrotehniko, Ljubljana, 1987.
8.	Okorokovv N. W.: ELEKTRISCHE SCHMELZOFEN FUR DIE EISENMETALLURGIE, VEB Verlag Technik, Berlin, 1953,
9.	Benkowsky G.: INDUKTIONSERVVARMUNG, VEB Verlag Technik, Berlin, 1965.
10. Hinton E. Owen D. R. J.: FINITE ELEMENT COMPUTATI-ONS, Pineridge Press Limited, Svvansea, 1979
ZUSAMMENFASSUNG
Mit der sehnellen Entvvicklung der Rechentechnologie, vor allem auf dem Gebiet der Personenrechner vverden in der Inge-nieurpraxis immer ofter numerisehe Methoden bei der Losung konkreter Probleme der Feldlosungen angevvendet, vvelehe durch Differentialgleihungen besehrieben sind. Heutzutage ist solehe Art der Projektierung aus okonomisehen und Oualitats-grunden gerechtfertig und vvird notig.
Unsere Aufgabe vvar Losung nichtstationarer Temperaturfel-der beliebiger Objekte und die Herstellung von entsprechen-dem Rechnerprogramm. Dabei konnen als Randbedingungen eindrangende Temperaturen oder VVarmestrome oder innere Ursprunge des VVarmestromes definiert vverden, die durch elek-trisehen Strom verursacht vverden. So konnen (hinsiechtlich auf das Temperaturfeld) Objekte wie z.B. Kuhlkorper, erwarmen von Metali und Isolationsteilen Temperaturbeanspruchung elek-
triseher Sammelleitungen Induktionshartung und Induktions-schmelzung von Metallen projektiert vverden.
Die Lossungen dieser Objekte nach der Methode der Ende-lemente vverden in tabelariseher und graphischer Form erhal-ten. Die Geometrie des Objektes vvird automatiseh auf die Grosse des Monitors skaliert, in den Knotenpunkten der Elemente vverden die gesuehten Temperaturlosungen bzw. Far-benskala und in die Elemente die durchschnittlichen Temperaturen der Elemente, bzvv. die Werte elektriseher Stromdichte und der zugehorigen VVarmestrome ausgeschrieben. Ebenso erhalten vvir einem Auszug der Koeffiziente des Sistemes bei vvelehen die gegebene Losung erhalten vvorden ist. Nach Be-darf konnen nur bestimmte Koeffiziente geandert und die Ten-denz neuer Losungen beodachtet vverden. Alles das stellt ein gutes Mittel fur das Ingenieurprojektieren dar.
SUMMARY
The fast development of computer technology, especially vvith PCs, numerical methods of solving concrete problems connected to fields deseribed by differential equations are more and more used in engineering. Today sueh a way of de-signing is economically and qualitatively justified and it is be-coming a necessity.
Our task vvas to solve unsteady temperature fields of op-tional objeets and to prepare a suitable computer program. As boundary conditions there can be defined the imposed temperatures or heat flovvs or internal sources of heat flovvs vvhich are caused by electric current. Thus the objeets like cooling bod-ies, heating of metallic and insulator parts, temperature loada-
bility of electric buses, induetion hardening and induetion meiting of metals can be designed (according to temperature field).
Solutions for these objeets by the method of definite ele-ments can be obtained in tables or graphically. The geometry of object is automatically scaled on the size of sereen, the nodes of elements are marked vvith sought solutions of temperature or by colour scale vvhile in elements the average temperatures of elements or values of electric current densities and corre-sponding heat flovvs appear. Also coefficients of the system are vvritten out for the given solution. If necessary only single coefficients can be varied, and the trand of nevv solutions can be observed. Thus a good tool is obtained for engineering designing.
3AKflblOHEHPlE
BcneacTBne čbicrporo pa3BMTMR BbmMC/iMTe/ibHOM TexH0/i0rMM, npe>Kfle Bcero b 06/iacTM nepcoHanbHbix BbmncnuTe/iax, b mh-MeHepHOii npaKTHKe Bce 6o;iee ncno/ib3yiOT HyMepnHecKne MeTOflbi peujeHMfl KOHKpeTHbix npoč/ieMM peujeHMfl no/iePi, onn-caHHbix c flmtepeHunanbHbiMM ypaBHeHHflMM. TaKOH cnocoč npoeKTnpoBaHMF b HacTonmeM BpeMeHM bb/meTca SKOHOMM-4ecrto m KanecBeHHO ue;iecoo6pa3HbiM. a TaKMe Hy>KHbiM.
HaiiiMM 3aaaHneM HBM/iocb peaiMTb HecTaunoHapHbie TeM-neparypHbie no/ifiKaKMX-/in6o o6"beKTOB m noaroTOBKa coot-BeTCTByKDU4ePi BbNMC/lHTe/lbHOii npOrpaMMbl. ripu 3T0M KaK npeae/ibHbie ycnoBna MoryT 6biTb onpeaeneHbi ycTaHOB/ieHHbie TeMnepaiypbi h;im TennoBbie tokm m/im BHyTpeHHbie mcto^hhkm TennoBoro TOKa, npMMHHOPi K0T0pbix RB/ifieTCfl sneKTpkmecKMfi TOK. 3tMM CnOCOČOM MO>KeM npoeKTMpOBaTb (CMOTpfl Ha TeMnepaTypHoe no/ie) o6"beKTbi, KaK Hanp: ox/iaflMTe/ibHbie Te-
na, HarpeBaHne MeTa/i/iMHecKnx m n3onnunoHHbix HacTeii, TeMnepaTypHyio Harpy3Ky en. HaKonMTenePi, MHflyKLiMOHHyio 3a-Ka/iKy m nH/iyKunoHHyio n/iaBKy MeTan/ioB.
PeUjeHUH 3TMX 06"beKT0B no MeTOfly KOHUeBblX 3/ieMeHTOB no/iyMaeM b cJjopMe TaČJinubi m rpat})HKa. reoMeTpnfl očteKTa aBTOMaTMMecKO npHcnocoč/ifleTcn Be/iMMMHe SKpaHa, b y3Jiax 3/ieMeHT0B noflB/iflkdtcfl Tpe6oBaHMe peaieHMH TeMnepaTypbi, T. e. UBeTHan LUKana, a b 3neMeHTax cpeflHbix TeMnepaTypw 3/ieMeHT0B, t. e. pa3Mepw nnoTHoeTM 3/ieKTpnMecKnx tokob h pa3Mepbi cooTBeTCTBytomnx tokob Ten/ioTbi. Ha 3KpaHe takwe nonB^RioTCfl K03c|)ct>MuneHTbi cMCTeMbi npM KOTopofi Mbi nony4n-nn aaHHoe peujeHMe. Ilo noTpeČHoeTM MOMeM M3MeHfiTb TozibKo onpeaeneHHbie KO3<t>0nuneHTbi n CMOTpeTb TeHfleHumo hobwx peujeHMPi. Otum cnocočoM Mbi no/iy4nnn xopownii HHCTpyMeHT A/ih MH>KeHepHoro npoeKTMpoBaHun.