Uporaba PC preglednic s poudarkom na reševanju temperaturnih polj in polj mešanja taline Application of Computer Spreadsheets vvith Emphasis on the Solution of Temperature Fields and Fields of Melt Stirring M. Bolčina, Železarna Štore 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 nelinearnimi parcialnimi diferencialnimi enačbami. Rezultati, ki jih na ta način dobimo, bistveno bolj ustrezajo dejanskemu stanju, kot rezultati, ki smo jih dobili po klasičnih postopkih, kjer smo morali vrsto vplivnih faktorjev zanemariti, da smo določene obrazce oziroma procedure lahko uporabili. Edina ovira je, da moramo imeti na razpolago ustrezno programsko in aparaturno opremo. Application and adaptation of the so called spreadsheets is presented. They can be satisfactorially used in solving partial differential equations. Nowadays they are available practically for each home computer or PC. The method of final differences and iterations are used till in any field segment the desired accuracy is achieved. Use of this method is extremely simple for the solution of a suitable form of Laplace or Poisson differential equation. It must not be neglected that always and immediately also corresponding graphical presentation of the system solution is available. The procedure was illustrated by two examples, i.e. by the solution of Fourier heat transfer differential equation, and by the somewhat more demanding solution of Navier-Stokes differential equation vvhich was applied in estimating the stirring intensity of melt in an induction furnace. 1 Uvod Prikazali bomo uporabo in prilagoditev tabelaričnih kalkulatorjev oziroma t.i. preglednic (spreadsheet), ki jih lahko v ta namen s pridom uporabimo in so na voljo praktično za vsak hišni ali osebni računalnik. Pri tem je uporabljena metoda končnih diferenc in iterativni postopek do željenega pogreška v poljubnem segmentu polja. Uporaba za reševanje ustrezne oblike Lapacove oziroma Poissonove diferencialnih enačb za 2D je skrajno enostavna, pri novejših preglednicah pa tudi za tri in več dimenzional-ne probleme. Ne smemo zanemariti, da imamo vedno in takoj na voljo ustrezno grafično ponazoritev rešitve sistema. Postopek bomo ponazorili z dvema primeroma in sicer Fourierjeve diferencialne enačbe prenosa toplote in nekoliko zahtevnejšega postopka reševanja Navier-Stokesove diferencialne enačbe, s pomočjo katere smo ocenjevali intenzivnost mešanja taline v indukcijski peči. Tako lahko npr. s formulo SUM(B2..B133) v trenutku dobimo vsoto vseh vrednosti, ki so v stolpcu B na vrsticah od 2 do 133. Rezultat se izpiše v celici, v kateri smo napisali zgornjo formulo. Dobljeni rezultat lahko ponovno koristimo pri nadaljni obdelavi. Poleg vseh bistvenih matematičnih, logičnih, časovnih itd. funkcij in izrazov, ki so implementirane v posamezne preglednice, omenimo še možnost uporabe iteracij, ki jo uporabimo pri numeričnem reševanju parcialnih diferencialnih enačb, do željenega pogreška. Za uporabnika je zelo pomembno dejstvo, da je uporaba hitra, priročna, rezultati pa so ponazorjeni v grafični obliki. 3 Uporaba pri reševanju temperaturnega polja Kot prvi primer prikažimo način uporabe preglednic pri numeričnem reševanju Fourierjeve parcialne diferencialne enačbe, ki jo za stacionarno stanje preoblikujemo v Laplaceovo obliko, 2 Preglednice — tabelarični kalkulatorji (spreadsheet) Pri delu lahko uporabljamo preglednice, kot so npr. Microsoft Works, Excel, Lotus 123, QuatroPro, Symfony itd. Ker predpostavljamo poznavanje vsaj ene od omenjenih preglednic, podajamo samo njihov kratek splošni opis. Posamezne preglednice imajo specifične prednosti. Namenjene so predvsem hitri in enostavni obdelavi podatkov po znanem ključu. Posamezne celice preglednice imajo svoje naslove sestavljene iz črk, ki definirajo stolpec in številk, ki definirajo vrstico nahajanja celice v tabeli, npr.: Al, DF205, CY5698 itd. Vse celice so na začetku prazne in med seboj enakovredne. Vanje lahko pišemo besedila, številčne vrednosti ali formule, s katerimi povezujemo posamezna polja. div(fc grad T) + gc PT 6T Ji T 62T 6y2 8q_ SV T(x,y, z,t) —TT + -TV = <5: d) (2) kjer so parametri k, q, č lahko krajevno, časovno in temperaturno odvisni. Ker pri tem uporabljamo diferenčno metodo, podajamo kratek opis izpeljave za omenjeni primer dvodimenzionalnega polja, kot kaže slika 1: A Tn A y Tn — To sledi: A 7 V A j-A Tu> Ax A Ts A y Ar _ To - T t Tw — To h To - Ts B T ATw Ar A 7V A r Mt Aa7 A y ATn Ay ATs Kot rezultat dobimo iskano temperaturo: Tw + Te + Tn + Ts + Qh2 To = - (3) (4) za poljubno vozlišče (celico) končnih difrenc znotraj obravnavanega polja. Podoben pristop lahko uporabimo tudi za časovno odvisna in večdimenzionalna polja. Na osnovi te izpeljave lahko z. pomočjo iteracij uporabljamo to metodo v prej omenjenih preglednicah. V splošnem nimamo težav z divergiranjem h končni rešitvi (velja za eliptični tip diferencialnih enačb). Tn -« T« . i- To h > —i Te Ts 10 n 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 STACIONARNO 2D TEMPERATURNO POLJE 500.00 500.00 500.00 500.00 500.00 1100.00 jlOO.OO i 100.00 |100.00 ■ 100.00 369.15 401.37 416.34 422.46 280.22 319.98 341.53 351.04 295.01 210.91 168.39 220.84 256.79 278.76 289.12 141.81 177.98 207.57 227.62 237.56 120.87 141.69 167.89 186.59 100.00 100.00 100.00 100.00 100.00 135.72 154.90 120.09 133.59 111.05 119.18 104.94 106.70 110.69 195.95 163.68 140.21 123.36 422.46 351.04 269.12 237.56 195.95 163.68 140.21 123.36 110.69 100.00 100.00 100.00 USTREZNI ALGORITEM (50 iteracij) 500.00 500.00 500.00 500.00 500.00 jlOO.OO i 100.00 i 1100.00 i 100.00 !100.00 (C17+B18+C19+018) /4 obrazec za C20 velja za vse celice v polju 100.00 100.00 100.00 100.00 100.00 11AJ.UU 1UU.UU 1UJ.UU F17 F18 F19 F20 F21 F22 F23 F24 F 25 Slika 1. Osnovna mreža končnih diferenc. Figure 1. Basic net of final differencies. Vzeli smo primer, kjer ploščinskemu homogenemu elementu specifične oblike (npr. kovinski plošči v obliki črke L), ki jo vizuelno ponazorimo v preglednici z. podajanjem robnih in začetnih pogojev in potrebnim številom celic glede na potrebno natančnost. Naša naloga je določiti pripadajoče temperaturno polje. V našem primeru naj bodo robni pogoji na vseh stranicah 100°C, razen na zgornjem robu, kjer predpostavimo konstantno temperaturo 500°C in na desni stranici, kjer predpostavljamo idealno toplotno izolacijo (slika 2). Po preoblikovanju parcialne časovno neodvisne diferencialne enačbe v ustrezno numerično obliko (po metodi končnih diferenc s korakom h), vidimo, da je temperatura v posamezni točki, ki jo ponazarja vrednost v celici, za homogene kovinske plošče enaka povprečju ob-dajajočih temperatur in toplotnemu izvoru v tej točki. Temperature na zunanjih robovih so konstante, razen na desnem robu, kjer je temperatura enaka vrednosti v sosednji točki (celici) v notranjosti elementa. Po večih interacijah pridemo do zadovoljivega pogreška, ki pa je poleg števila iteracij odvisen tudi od kvadrata koraka delitve h, kot kaže slika 2. Tako v poljubno celico preglednice napišemo splošno enačbo 4, ki jo nato prekopiramo v druge celice glede na geometrično obliko opazovanega polja. (Pri tem preglednice Slika 2. Primer uporabe preglednic za določevanje dvodimenzionalnega temperaturnega polja. Figure 2. Example of use of spreadsheets in determining a two-dimensional temperature field. avtomatično prevzamejo relativne naslove obrazcev, torej ustrezne sosednje celice posamezne prekopirane celice). Če imamo opravka s konstantnimi zunanjimi temperaturami vnesemo kot robne parametre ustrezne številčne vrednosti, drugače pa ustrezne obrazce, ki opisujejo robne efekte toplotne izolacije, toplotne prestopnosti, sevanja itd. Tudi tu si pomagamo s kopiranjem istosmiselnih obrazcev. 4 Uporaba pri reševanju polj mešanja taline V električnih prevodnikih, ki so v izmeničnem magnetnem polju, se inducira električni tok. Komponenta magnetnega polja, ki je pravokotna na ta električni tok, povzroči t.i. Lorentzovo silo. V tekočih prevodnikih je ta sila vzrok za generiranje toka te tekočine in s tem mešanje fluida. Ta osnovni princip izkoriščamo pri elektromagnetnem mešanju taline. Pri tem uporabljamo izmenično (rotirajoče ali utripajoče) magnetno polje, ki ga ustvari ustrezen zunanji elek-tromagnet. Tema tega primera je obdelati elektromagnetno mešanje taline v indukcijski lončni peči s pomočjo omenjenih preglednic. Ker smo predpostavljali osnosimetrični problem (lončna peč), smo s tem zanemarili komponente hitrosti v smeri, ki bi lahko nastale kot posledica turbolenčnih efektov. Torej je problem najlažje rešljiv v cilindričnih koordinatah. Čeprav je vzbujevalno magnetno polje časovno odvisno s frekvenco napajalne napetosti na induktorju, je komponenta Lorenzeve sile, ki premika fluid, stalna (Mof-fatt je dokazal, da je časovno odvisna komponenta kompenzirana s časovno odvisnostjo v pritisku fluida, kar pa presega okvir tega dela). Polje sil v talini Slika 3. Prikaz sil na enoto taline. Figure 3. Presentation of forces per unit melt. Razlika sil na robovih poljubno majhnega geometrijskega telesa v talini povzroči njegovo vrtinčenje. Zalo po numerični poti določimo rot F(r, r) v cilindričnih koordinatah, ki ima smer 1 ip. RF(r,z) = rot F(r, z) /01 > l
.
(7)
(8)
Operator nabla V in Laplacov operator sta v tem primeru dvodimenzionalna operatorja in sta definirana kot:
0 0
v =
V2 = +
Or Or Oz2
V2 « + (9)
Or Oz2
Končno dobimo v primeru konstantne viskoznosti (nestisljivega fluida) s pomočjo Navier-Stokesove enačbe 6 izraz:
RF{r^ + T?V2WC = -
0V 0u!r OV du>
Oz
Or Or
■sr- 1101
6 Tehnika računanja Navier-Stokesove enačbe s pomočjo preglednic
V tem primeru moramo biti pri sami tehniki računanja po opisanih metodah posebej pazljivi. Razlike pri velikih številih, ki pri tem izračunu nastopajo v posameznih celicah preglednice, so relativno majhne in velika nevarnost je, da reševanje po iterativnem postopku ne konvergira. (Paziti torej moramo, da v lastni matriki računanega sistema prevladujejo diagonalni koeficienti.) To velja predvsem za preglednice starejše generacije, ki pri računanju upoštevajo manjše število decimalnih mest.
V nadaljevanju izhajamo iz sistema naslednjih dveh parcialnih diferencialnih enačb:
Or2 Oz2
<92uc 02uic ~ +
Or- Oz2
Or Oz
d V 0uc Oz Or
5 Dinamika fluida
Osnovna enačba za stalen tok nestisljivega fluida, kar opisuje obrazec V V = 0 je:
(H) - RF(r,z) (12)
z. dvema neznankama uic in Enačbi lahko rešujemo po metodi, ki smo jo opisali v prejšnjem poglavju. Če predpostavimo 1} — oc , potem je rešitev spodnje enačbe relativno enostavna:
(6)
UcO
+ U!ce + U cn +
(13)
S to rešitvijo pa sproti iterativno rešujemo tudi zgornjo enačbo:
pri čemer velja:
+ l-e + ^n + V, +u>c0h2
Z rešitvami polja tokovne funkcije \p0 =
zadostnem številu iteracij (Vr" 1 — ^o"
-D
(14)
po
< emax), lahko
določimo hitrostno polje taline po obrazcu:
V = V^ = lr—— + lz—. (15) Ar Ar
Seveda pa moramo upoštevati dejansko vrednost za viskoznost t?(t). ki bistveno vpliva na končno rešitev. Zato v postopek iterativnega računanja preko ustreznih celic v preglednici in pripadajočih matematičnih povezav vključimo medsebojni vpliv med obema enačbama, kakor tudi vpliv polja elektromagnetnih sil, kar definira izraz:
l fdV dioc 7] \ dr 8z
OV DiJr , ,
(16)
V tem primeru postane sistem za realne vrednosti i/ večkrat nestabilen, oziroma rešitve v splošnem ne konver-girajo h končni vrednosti. Vzrok ni samo v matematičnem, ampak tudi fizikalnem ozadju opisanega primera (efekt tur-bolenc). Rešitev problema smo našli po dveh poteh. Ker vemo, da lahko imamo laminami tok tudi preko meje Re > 10'\ če eksperiment oz. povečevanje hitrosti fluida opravljamo silno previdno (zelo počasi, brez zunanjih tresljajev itd.), smo prvi način rešitve našli po tej poti. Postopek iterativnega računanja smo namreč pričeli z zadosti veliko viskoznostjo r/. Ko smo pri tej vrednosti r/ dobili dovolj majhen pogrešek e — 0, smo ?; ponovno zmanjšali. Ta postopek smo avtomatizirali tako, da je bil novi rjn+1 pri
e(r,z) —" 0:
Vn + 1 = Vn - 0.01(»7„ - '/dejanski)-
17;
Postopek računanja je zelo dolgotajen, saj je potrebno izvesti več tisoč iteracij (v preglednici Excel ob mikroprocesorju 386 s koprocesorjem je potrebno nekaj ur računanja). Pri tem si lahko pomagamo s t.i. makro ukazi.
Te težave presežemo s takoimenovano "forvvard" oz. "backward" diferenčno metodo, namesto centralnih diferenc za prve odvode. Ta način nam zagotovi pozitivni prispevek k diagonalnim členom matrike koeficientov. Napaka pri diskretizaciji na razdalji h je tako O(h), pri uporabi centralnih diferenc pa bi bila 0(h2). Opisani postopek računanja sta prva predlagala Richards in Crane. Računanje smo izvedli v preglednici Excel po naslednjih splošnih obrazcih:
+¥em) + - iVo)
= n - «>2)^r+
+cM?] + Cs^T] + (cw - +
+(Cn-Unyc^ + h2RF/i1), (18)
Cn = 1 -p C, = 1 + p
Ce = 1 + 9 C w = 1-9
Un = Cn + (p) u, = C,+(p)
Ue = Ce + (q) uw = Cw+(q) (19)
vir _ vi/ ^n -
P 4r/ 1 4 n
Wi = 0.5 IV 2 = 0.5.
Konkretna razporeditev celic je prikazana v tabeli 1. Pri tem načinu uporabe preglednic pridemo veliko hitreje do rešitve tokovne funkcije V opisanem primeru so robni pogoji za na stenah lonca konstantni (npr. 0), saj je tudi hitrost taline na teh mestih enaka 0. Konkretni izračun za naš primer je podan spodaj. Slika 4 prikazuje tipični način podajanja rezultatov s tokovnicami in tangencialnimi vektorji hitrosti fluida.
■ V ............ I \ ............. J :
P/ ♦ T.-*?.;
.................. ................*
» * .........<
• S . - ......N • ••• • ........•• \ :
# • %
\ : :: ; /
li;:?/ ~ *S!!!!i» / * \Vin
II' SM\M{* ' lil
HU-V- * * jjM.
v.l;. •. ___// .• : / : •.: •. \ N___.• :i:>:
\j h
• • \ (i .. i, v«./ .......... x • •
• .........•• ................. :
k ••••••••••••••••••• .( |a *"•••••••••••••**** M
....................milil"---•,,|lHIIII|llllll.,illlllHlrf- —
Slika 4. Tokovnice taline pri elektromagnetnem mešanju v indukcijski peči.
Figure 4. Melt tlowlines in electromagnetic stirring in an induetion furnace.
Ker lahko zanemarimo povezavo vpliva gibanja taline na elektromagnetne razmere (enačba V x B —» 0), lahko Maxwellove enačbe (oziroma električni model indukcijske peči) rešujemo na enak način neodvisno od problematike fluida. Med posameznimi točkami, kjer nimamo definiranih razmer, lahko uporabimo postopke interpolacije.
rot F = RF vhodni podatki
F15
A17 tokovna funkcija T 'trom+,=< l-w)T0n,+ W./4 { m*l +
1 t w ■ n ,
+ Hr rn + 4, m_h2
e s cO J
w,=0.5, w.,=0.5 kinematiCna viskoznost 4 10 ''Itn^/sl
TT7~
p = ( fe - >Cw)/4t)
_L32
Ni 7
q
q = ( fn - r»)/4n
032
A33
C n= 1 - p
f47
A63
Ca- 1 + p F92
-
Uw=Cw+abs(ql _F107
rxro5-
Us=Cs+abs( p) _FI22I
a4fe
Cw= 1 - q
Ffil-
ca a
(>= 1 + q
1.1)2
"T7S-
lIn=Cn+abs< p)
_F77_
77H-
Ue=C<»+abs(q)
_L107
itm-
Uo
Uo=( Un+Uw+Ue+Us)/4 _F137
75T3S
cirkulacija uco _F154
-4o1 + 1) = (1 - + > + UnJ™+» + C>