SIMULACE TUHNUTÍ A CHLADNUTÍ PŘEDLITKU NA ZPO Příhoda Miroslav – Molínek Jiří – Vu Quoc Hung – Pyszko René Vysoká škola báňská – Technická univerzita Ostrava, 17. listopadu 15, 708 33 Ostrava – Poruba, ČR Abstract For solving different problems from thermodynamics there are available standard program packets, e.g. Ansys, Fluent etc. The universality of these programs brings also certain disadvantages. For example, it is hardly to describe the complexity of border conditions occurring in a continuous casting process. Otherwise, their prices are often high. Specialized software for solving thermal field of square and circular blanks based on the explicit finite difference method is described in this paper. The choice of solving method and border conditions is also discussed here. Expressions for numerical stability conditions of round profile are presented. The programs enable computation of blank temperature field from the steel level in the mould down to cutting machine on given CC and also the length of liquid substance. The created programs enable the user to provide data entry of initial conditions of solved task in communicative way and present results in graphical representation. An example of calculation for round profile of diameter 320 mm is presented. The results obtained are in good accordance with temperatures measured on real CC. 1 ÚVOD Tepelné procesy v předlitku dávají přehled o výrobnosti ZPO i o kvalitě předlitku. Znalost teplotního pole předlitku umožňuje řešit problematiku protrhnutí utuhlé kůry, vnitřní strukturu, povrchovou jakost, mechanické vlastnosti předlitku a optimalizovat jednotlivé parametry výrobní technologie. Standardní programové balíky, jakými jsou např. Ansys, Fluent aj. jsou univerzální a jejich aplikace na proces plynulého odlévání je komplikovaná. Nezanedbatelnou skutečností není ani jejich relativně vysoká cena. Z těchto důvodů byl vytvořen vlastní, specializovaný software pro řešení teplotního pole předlitků. 2 TEPLOTNÍ POLE PŘEDLITKU Kinetiku teplotního pole tuhnoucího a chladnoucího předlitku popisuje Fourierova – Kirchhoffova (F.-K.) rovnice. Naše dřívější práce (viz např. [1]) prokázaly, že konvekční proudy oceli, které se projevují do vzdálenosti 1 až 2 m od hladiny taveniny v krystalizátoru, významně neovlivňují dobu tuhnutí předlitku. Vliv konvekce lze zahrnout do hodnoty tzv. ekvivalentního součinitele tepelné vodivosti a místo rovnice F.-K. řešit klasickou Fourierovu rovnici ve tvaru ∂ (t ⋅ c p ⋅ ρ ) ∂τ kde t τ cp ρ
= div(λ ⋅ ∇t ) + qV je – – –
(W ⋅ m −3 )
teplota (°C), čas (s), měrná tepelná kapacita (J.kg-1.K-1), hustota (kg.m-3),
(1)
λ qV
– součinitel tepelné vodivosti (W.m-1 .K-1), – intenzita vnitřního tepelného zdroje (W.m-3).
2.1 Výběr metody řešení Kinetiku teplotního pole tuhnoucího a chladnoucího předlitku lze efektivně řešit pouze některou z numerických metod. S přihlédnutím k okolnosti, že řešené předlitky měly jednoduchý pravoúhlý resp. kruhový tvar, byla vybrána síťová metoda. Konkrétně to byla explicitní diferenční metoda, která je výhodná zejména v těch případech, kdy jde o rychle probíhající děje. Za těchto podmínek má přednost před metodou implicitní. Z kapitoly 2.2 vyplývá, že fyzikální a povrchové podmínky řešení úlohy jsou složité. Jejich správná volba je nesnadná a zejména na ní závisí, jak budou vypočtené výsledky korelovat s teplotním polem reálného předlitku. Není tudíž efektivní vytvářet „přesné“ modely, když vlastní simulace bude probíhat s ne zcela exaktními vstupními parametry. Konkrétně se jedná např. o rozhodnutí, zda řešit teplotní pole předlitku jako rovinnou či prostorovou úlohu. Trojrozměrné řešení je podstatně náročnější na strojový čas výpočtu, přičemž vždy nezaručuje přesnější výsledky. Naše dlouhodobé zkušenosti např. ukazují, že větší vliv na přesnost řešení má správná volba prostorového dělení sítě, než skutečnost, zda se chladnutí a tuhnutí předlitku řeší jako 2 D či 3 D problém. Z tohoto důvodu se v dalším textu uvádí řešení rovnice (1) ve dvou směrech. Při výpočtu se předpokládá symetrické ochlazování a tudíž se řeší teplotní pole v jedné polovině předlitku. 2.2 Podmínky jednoznačnosti Při výpočtu teplotního pole předlitku je důležitý správný výběr podmínek jednoznačnosti řešení rovnice (1), tedy podmínek geometrických, fyzikálních, počátečních a povrchových. Z této množiny je problematické určení podmínek fyzikálních a zejména povrchových, neboť podmínka geometrická plyne z konkrétního tvaru odlévaného předlitku a podmínka počáteční souvisí s teplotou oceli v mezipánvi. Fyzikální podmínky mj. zahrnují součinitel tepelné vodivosti, měrnou tepelnou kapacitu a hustotu oceli. K vyjádření funkční závislosti těchto vlastností na teplotě a chemickém složení ve tvaru polynomů třetího stupně se využilo metodiky, založené na experimentálních údajích BISRA. Rovnice (1) byla řešena jako kvazilineární, tj. za předpokladu, že fyzikální vlastnosti oceli jsou ve zvoleném časovém kroku řešení ∆τ konstantní. K fyzikálním podmínkám také patří vnitřní tepelný zdroj o intenzitě qV, který je při tuhnutí předlitku představován uvolňováním latentního tepla tuhnutí L. V dnes už klasickém postupu dle Mizikara se vliv latentního tepla respektuje prostřednictvím zvýšené hodnoty měrné tepelné kapacity v intervalu tuhnutí. Esser a Kruse [2] ukázali, že takový způsob řešení vede k větším délkám tekutého jádra předlitku oproti skutečnosti. Na katedře tepelné techniky VŠB–TU Ostrava se osvědčil způsob, který byl vícekrát verifikován i měřením na reálných licích strojích a spočívá v přiřazení „teplotního“ zdroje každému uzlovému bodu řešené oblasti. Vydatnost zdroje se určí z podílu měrného skupenského tepla tuhnutí a měrné tepelné kapacity. Vzhledem k charakteru úlohy byly při výpočtu používány povrchové podmínky II. resp. III. druhu. Způsob jejich určení se liší podle toho, zda je jedná o primární, sekundární či terciární oblast chlazení. Vychází se z kombinace teoretického výpočtu a experimentálního měření, přičemž v prvních dvou oblastech převažoval experiment a ve třetí oblasti se teoretické hodnoty korigovaly experimentálními výsledky. Transport tepla v krystalizátoru je teoreticky obtížně stanovitelný, takže intenzita odvodu tepla z tuhnoucí oceli, vyjádřená hustotou tepelného toku či součinitelem prostupu tepla, byla měřena na reálných ZPO. Principiálně lze aplikovat dva přístupy, z nichž každý má své
přednosti i nedostatky [3, 4]. Jeden, využívající měření parametrů chladicí vody krystalizátoru, je poměrně jednoduchý a přesný, ovšem umožní zjistit pouze střední hodnoty tepelného toku na příslušné desce krystalizátoru. Druhý, založený na měření teplotních gradientů v různých místech po výšce a obvodu pracovního povrchu krystalizátorové vložky, dovoluje vypočítat hodnoty místních tepelných toků, ale měření je náročné na přípravu i vlastní provedení experimentu. Optimální je současné použití obou postupů. V sekundární oblasti chlazení bylo pro určení hodnoty součinitele přestupu tepla ostřikem α využito zejména laboratorního výzkumu. Na katedře tepelné techniky byl navržen a postaven teplý model sekundární oblasti chlazení, pracující na principu přímožhavené sondy [5]. Pracovní část sondy, vyrobená z kanthalového pásku, se odporově vyhřívá elektrickým proudem a měří se elektrický příkon, potřebný pro udržení konstantní teploty sondy. Elektrický příkon dodávaný do měřicí plošky se rovná tepelnému toku odváděnému z této části sondy. Od měřeného příkonu sondy se odečítají ztráty, představované tepelným tokem odváděným do držáku sondy a přívodních vodičů. Ztrátový tepelný tok se určí cejchovacími měřeními, během nichž je aktivní část sondy tepelně izolovaná. Ztráty závisí na teplotě sondy a teplotě okolí a proto cejchování probíhá opakovaně před a po vlastním měření. Přímé měření elektrických veličin je relativně jednoduché a přesné. Hlavní výhodou sondy je její malá tepelná setrvačnost. Sonda se pohybuje pomocí posuvného mechanizmu ve dvou kolmých směrech, přičemž měření probíhá v n vertikálních a m horizontálních řezech. Měření se v každém řezu několikrát opakuje a naměřené hodnoty napětí a teploty jsou registrovány měřicí ústřednou. Hodnoty α (nebo analogicky hustoty tepelného toku q) jsou po zpracování tabulkovým kalkulátorem uloženy ve tvaru matice s prvky αi,j, resp. qi,j. Pro různé tlaky chladicí vody je možné získat funkční závislost lokální i globální intenzity chlazení na tlaku chladicí vody či teplotě ochlazovaného povrchu. Použitý teplý model s přímožhavenou sondou je v porovnání s ”klasickým” deskovým teplým modelem energeticky méně náročný, neboť se ohřívá jen malá měrná ploška sondy. Teplotní rozsah měření je široký a je v zásadě ohraničen pouze teplotní odolností materiálu sondy. Protože vyhodnocení se provádí na principu měření elektrických veličin, není nutné znát závislost součinitele tepelné a teplotní vodivosti materiálu sondy na teplotě. Přesnost měření je limitována přesností měřicích přístrojů pro cejchování a vlastní měření. V terciární oblasti chlazení se radiační složka součinitele přestupu tepla určuje výpočtem ze Stefanova – Boltzmannova zákona. Konvekční složka, kde dominantní roli hraje přirozené proudění, se řeší z kriteriální rovnice, udávající závislost mezi kritériem Nusseltovým a Rayleighovým. Teoretické hodnoty byly nepřímo ověřovány měřením povrchovým teplot předlitku. 2.3 Numerická stabilita O explicitní diferenční metodě je známo, že může být, na rozdíl od metody implicitní, numericky nestabilní. Taková situace nastává, není-li dodržena určitá závislost mezi časovým a prostorovým řešením V dalším textu bude ukázáno řešení teplotního pole u předlitku kruhového průřezu, tudíž se blíže zmíníme o podmínce numerické stability v cylindrické soustavě souřadnic. U kruhové oblasti rozlišujeme tři typy elementů, a sice elementy vnitřní, vnější a středové. Po obvodu se oblast dělí na výseče o středovém úhlu ∆ϕ a po poloměru na mezikruží o šířce ∆r. Výjimku tvoří vnější a středové elementy, jejichž šířka je ∆r/2. Každý element je označen indexy i ∈ <0,m> (po obvodu) a j <0,n> (ve směru poloměru). Povrchové elementy mají index j=0, u středových pak platí j=n. Pro jednoduchost se předpokládá, že pro všechny
elementy sítě platí, že během časového intervalu ∆ τ je součin ρ .cp a hodnota λ jenom funkcí souřadnic i, j. S využitím druhého zákona termodynamiky je možno odvodit vzájemnou závislost mezi časovým krokem a prostorovým dělením, zaručující numerickou stabilitu. Pro vnitřní elementy platí ∆τ ≤
ρ i , j ⋅ c i , j ⋅ ∆ϕ ⋅ R j ⋅ ∆ r λi , j ⋅ ∆r 1 1 + λi , j ⋅ ∆ϕ ⋅ + R Rj ∆ϕ ln j R j ⋅ sin ln R R j −1 2 j +1
(s )
(2)
U středových elementů vychází ∆τ str ≤
ρ i , n ⋅ c i , n ⋅ ∆ϕ ⋅ ∆ r 2 2 ∆ϕ 8 ⋅ λi ,n ⋅ + ∆ϕ ln 0,25 sin 2
(s )
(3)
U elementů povrchových lze odvodit ∆τ pov ≤
ρ i , 0 ⋅ c i , 0 ⋅ ∆ϕ ⋅ ∆ r ⋅ ( R − ∆ r / 4 ) λi , 0 ⋅ ∆r 2 + λi ,0 ⋅ ∆ϕ ⋅ + 2 ⋅ α ⋅ R ⋅ ∆ϕ R ∆ϕ ln R ⋅ sin R − ∆r 2
(s )
(4)
Časový krok stabilního numerického řešení, který odpovídá minimální hodnotě, určené z nerovností (2) až (4), je volen a kontrolován během výpočtu automaticky. Pro běžné podmínky na ZPO to bývá hodnota vypočtená z rovnice (3). Kritérium numerické stability, uvedené např. v literatuře [6], není dostačující. 3 POČÍTAČOVÉ PROGRAMY Byly vytvořeny dva programy, řešící teplotního pole jak u obdélníkových, tak u kruhových předlitků. Princip práce i obsluha obou programů je shodný, takže budou komentovány společně. 3.1 Popis programu Program vyžaduje vstupní údaje ve dvou datových souborech. Jde o hlavní datový soubor, obsahující mj. teplotu oceli lité do krystalizátoru, teplotu likvidu a solidu odlévané oceli, latentní teplo tuhnutí oceli. Následuje tloušťka předlitku a počet elementů po tloušťce, šířka předlitku a počet elementů po půlce šířky, rychlost lití, emisivita, Největší část souboru zaujímají povrchové podmínky v primární a sekundární zóně chlazení, včetně jejich délek. U terciární zóny se zadává pouze celková délka zóny. Druhý datový soubor obsahuje údaje o chemickém složení odlévané oceli. Oba datové soubory lze připravit předem resp. upravit pomocnými programy. Nejsou-li po spuštění programu nalezeny oba soubory v aktuálním
adresáři, jsou nabídnuty soubory prázdné a ty je třeba s pomocí nápovědy doplnit. Na konci zadávání vstupních údajů jsou soubory pojmenovány, uloženy a jsou k dispozici pro případné další použití. Ukázka jedné ze stránek vstupního souboru je na obr. 1.
Obr. 1 Ukázka jedné z obrazovek vstupního souboru Po zadání konkrétního složení uhlíkové nebo nízkolegované oceli program určí hodnotu tzv. ekvivalentního uhlíku (Cekv). Pokud je u uhlíkové oceli Cekv < 7 a u nízkolegované oceli Cekv < 8, vypočtou se teplotní závislosti fyzikálních vlastností oceli (součinitel tepelné vodivosti, měrná tepelná kapacita, hustota oceli) ve formě polynomů třetího stupně. Pro vyšší hodnoty Cekv a pro legované oceli neumí program teplotní funkce určit a koeficienty jednotlivých polynomů je potřeba zadat ručně. Koeficienty u lineárního až kubického členu mohou být i nulové (konstantní hodnota dané fyzikální veličiny). Na poslední obrazovce před zahájením výpočtu se znázorní síť uzlových bodů, kterou je pokryta řešená oblast. Podle nápovědy se vyberou uzlové body, v nichž se budou během výpočtu vykreslovat průběhy teplot. Je vhodné volit počet uzlů tak, aby kreslený obrázek byl přehledný. Při větším počtu vybraných uzlů mohou teplotní křivky splývat. Během výpočtu se vytváří několik výstupních souborů, které obsahují teplotní pole předlitku, měrnou entalpii předlitku, hustotu tepelného toku z jednotlivých povrchů do okolí, tloušťku utuhlé kůry na malém a velkém rádiusu. Všechny tyto údaje jsou ukládány v závislosti na vzdálenosti od hladiny oceli v krystalizátoru. Dále program vytvoří čtyři grafické soubory, což jsou obsahy čtyř posledních grafických obrazovek. Tyto soubory mají formát BMP, lze je později prohlížet standardním programem typu „picture viewer“. Jsou to obrazovky teplotního průběhu vybraných bodů, měrné entalpie předlitku, průměrné hustoty tepelného toku z povrchu předlitku do okolí a tloušťky utuhlé kůry. 3.2 Příklad výpočtu Pro demonstraci programu bylo vybráno odlévání kruhového průřezu ∅ 320 mm z oceli jakosti ČSN 95 ČSD-VK s 0,391 hmot.% C. Ze vstupních údajů lze uvést: teplota oceli v krystalizátoru 1518 °C, teplota likvidu 1493 °C, teplota solidu 1468 °C, rychlost lití 0,91 m.min-1, délka krystalizátoru 0,7 m, délka sekundární zóny 4,76 m, délka terciární zóny
29,14 m. Počítaná polovina průřezu byla rozdělena na elementy o rozměrech ∆r = 16 mm a ∆ϕ = 18°. Ukázka vypočteného teplotního pole v sedmi vybraných bodech průřezu je na obr. 2. Originální obrázek je barevný, zde uvedená černobílá kopie má sníženou vypovídací hodnotu. Přesto je zřejmý reohřev povrchové vrstvy po výstupu z krystalizátoru, kolísání povrchové teploty vlivem chladicích trysek a opětovný nárůst teploty povrchu na počátku terciární zóny. Z obrázku lze odečíst i metalurgickou délku, která je pro zadané technologické podmínky necelých 16 m.
Obr. 2 Teplotní průběh ve vybraných bodech kruhového předlitku
metalurgická délka (m)
Vypočtené povrchové teploty byly porovnávány s teplotami změřenými na provozním blokovém ZPO. Ve vzdálenosti 10,3 m od menisku byla vypočtená teplota 1015 °C, zatímco naměřená činila 1013 °C. Ve vzdálenosti 16,2 m bylo vypočteno 984 °C a naměřeno 967 °C. Je zřejmé, že zjištěný teplotní rozdíl je technicky přijatelný.
18,5 18
l = 23,851.v - 5,759 R2 = 1
17,5 17 16,5 16 15,5 0,9
0,92
0,94
0,96
0,98
1
-1
rychlost lití (m.min )
Obr. 3 Vliv licí rychlosti na metalurgickou délku
Popisovaný program je možno využít i k simulaci vlivu změn technologických parametrů na rychlost tuhnutí a chladnutí předlitku. Jako příklad je pro výše zmíněný kruhový blok uveden na obr. 3 vliv změny licí rychlosti na metalurgickou délku. Ukazuje se, že rychlost odlévání ovlivňuje metalurgickou délku lineárně a její zvýšení o 0,1 m.min-1 se projeví zvětšením délky tekutého jádra o přibližně 2,4 m.
Obdobně lze programem simulovat vliv rychlosti lití či licí teploty na tloušťku licí kůry při výstupu z krystalizátoru, vliv chemického složení oceli na metalurgickou délku apod. V součinnosti s laboratorním výzkumem program umožňuje také např. posoudit rozmístění trysek v sekundární oblasti chlazení. 4 ZÁVĚR Pro řešení teplotního pole předlitku na ZPO byla použita explicitní diferenční metoda a sestaveny programy pro řešení pravoúhlého a kruhového formátu. Výpočet fyzikálních vlastností oceli byl založen na experimentálních údajích BISRA. K určení povrchových podmínek úlohy v krystalizátoru se využilo výsledků měření na provozních licích strojích. Popis odvodu tepla v sekundární oblasti vycházel z laboratorního výzkumu na teplém modelu. K zajištění numerické stability výpočtu v cylindrických souřadnicích byly odvozeny nové vztahy, udávající závislost mezi časovým krokem a prostorovým dělením. Funkce programu byla ukázána na odlévání kruhového předlitku ∅ 320 mm. Vypočtené povrchové teploty se dobře shodují s teplotami, zjištěnými experimentálně. Vytvořené programy jsou efektivním nástrojem k posouzení vlivu změn technologických parametrů na kinetiku teplotního pole předlitku.
LITERATURA [1] Rédr, M. – Příhoda, M. – Molínek, J.: Vliv rychlosti odvodu tepla přehřátí tekuté oceli na tuhnutí slitku při plynulém odlévání oceli. Hutnické listy XXXII, 1977, č. 1, s. 13–16. [2] Esser, F. – Kruse, H.: Beitrag zur Berechnung der thermischen Erstarrungsvorgänge durch Rechneranwendung Grundlangen. Neue Hütte, 17 Jg, Heft 11, November 1972. [3] Příhoda, M. aj.: Problematika odvodu tepla v krystalizátoru při lití kruhových bloků. In.: Sborník 7. mezinárodního metalurgického symposia METAL 98 – 2. díl. Ostrava, květen 1998, s. 31–37. ISBN 80–8612214–X. [4] Příhoda, M. aj.: Teplotní poměry u krystalizátoru bramového typu In.: Sborník 8. mezinárodního metalurgického veletrhu a symposia METAL 99 – II. díl. Ostrava, květen 1999, s. 114–121. ISBN 80–85988–36–4. [5] Příhoda, M. aj.: Stanovení součinitele přestupu tepla v sekundární oblasti chlazení při plynulém odlévání oceli. Hutnické listy LIV, 1999, č. 7/8, s. 29–32. ISSN 0018–8069. [6] Šmrha, L.: Tuhnutí a krystalizace ocelových ingotů. 1. vyd. Praha: SNTL, 1983. 308 s.
Práce byly řešeny v rámci komplexního projektu technologické inovace plynulého odlévání oceli v ČR - evid. číslo 106/96/K032 za finanční podpory Grantové agentury ČR.
