NUMERICKÝ MODEL TUHNUTÍ KRUHOVÉHO PŘEDLITKU PRO ON-LINE MONITORING NUMERICAL MODEL OF ROUND BLANK SOLIDIFICATION FOR

METAL 2008 13. –15. 5. 2008, Hradec nad Moravicí ___________________________________________________________________________

NUMERICKÝ MODEL TUHNUTÍ KRUHOVÉHO PŘEDLITKU PRO ON-LINE MONITORING NUMERICAL MODEL OF ROUND BLANK SOLIDIFICATION FOR ON-LINE MONITORING David Dittela René Pyszkoa Pavel Fojtíka Miroslav Příhodaa Jiří Molíneka Michal Adamikb a

VŠB-TU Ostrava, katedra tepelné techniky, 17. listopadu 15, 708 33 Ostrava – Poruba, ČR b

TŘINECKÉ ŽELEZÁRNY, a. s., Průmyslová 1000, 739 70 Třinec – Staré Město, ČR

Abstrakt Znalost teplotního pole předlitku v průběhu jeho odlévání a chladnutí umožňuje řešit problematiku vnitřní struktury, povrchové jakosti, mechanických vlastností předlitku, metalurgickou délku, změnu tloušťky licí kůry pro různé licí rychlosti a přehřátí oceli. Z hlediska kvality odlévané produkce je nutno tyto parametry dále sledovat a optimalizovat. Na katedře tepelné techniky, Fakulty metalurgie a materiálového inženýrství, VŠB – TU Ostrava byl vytvořen program, který řeší teplotní pole předlitků explicitní diferenční metodou sítí. Podstata této metody spočívá v rozdělení předlitku na síť a v aproximaci základní diferenciální rovnice na odpovídající rovnici diferenční. Řešením je soustava algebraických rovnic, která je odvozena z příslušných okrajových podmínek. Pro řešení teplotního pole je tedy důležité stanovení podmínek jednoznačnosti řešení, konkrétně pak povrchových, počátečních, geometrických a fyzikálních podmínek. Vzhledem k charakteru úlohy jsou při výpočtu použity povrchové podmínky ve formě hustot tepelných toků, popřípadě součinitelů přestupu tepla, a to na základě toho, zda se jedná o primární, sekundární či terciární oblast chlazení. Pro vlastní výpočet 3D teplotního pole kruhového předlitku byla použita rovnice Fourierova – Kirchhoffova, která zohledňuje rychlost pohybujícího se předlitku. Program řeší teplotní pole ve všech zónách chlazení plynulého odlévání oceli. K vlastnímu výpočtu program vyžaduje parametry lití, chemické složení oceli, teploty solidu a likvidu, rozmístění a typy trysek. V rámci výzkumných prací je vyvíjen systém, který spojuje numerický model tuhnutí s funkcí predikce průvalu. Systém bude informačně propojen s technologickým procesem a bude dynamicky reagovat na změny provozních veličin. Systém umožní zvýšení výrobnosti, zlepšení kvality produkce a zamezení vzniku průvalu. Jeho výsledky mohou být dále využity pro optimalizaci řízení licího stroje. Abstract The knowledge of round blank temperature field during its casting and cooling makes possible to solve the problematic of inner structure, surface quality, blank mechanical properties, metallurgical length, casting crust thickness growth for different casting velocities 1

METAL 2008 13. –15. 5. 2008, Hradec nad Moravicí ___________________________________________________________________________ and steel overheating. From the view of cast production quality it is necessary to monitor and optimize these parameters. In the Department of Thermal Engineering, Faculty of Metallurgy and Materials Engineering, VSB – TU Ostrava was made program that solves blanks temperature field by using of explicit differential net method. Principle of this method consists in blank separation on net and in approximation of elementary differential equation on corresponding difference equation. The solution is system of algebraic equations that is deduced from relevant boundary conditions. For solution of temperature field it is also important definiteness solution conditions determination, concrete surface, initial, geometrical and physical conditions. In light of exercise character are, used by the solution surface condition in form of heat flow densities, eventually heat transfer coefficients on the basic, if it goes about primary, secondary or tertiary cooling zone. On heat flow density determination in primary zone and heat transfer coefficient determination in secondary zone is focused in the long term attention. There are given initial and geometrical conditions, because there is known not only steel temperature in tundish, but geometrical shape of CC blank. For the own 3D solution of round blank temperature field was used Fourier´s – Kirchhoff´s equation, which makes provision for velocity of moving blank. The given program solves temperature field in all cooling zones of continuous steel casting. For the own solution the program demands casting parameters, chemical constitution, liquid and solid temperature, location and type of nozzles. Within the researching project is developed system that connects numerical solidification model with function of rush prediction and will be connect with technological process and will dynamical answer to process signals and prevention of rush origin. Its conclusions can be continuously used for optimize of casting machine operating. 1. ÚVOD Metoda sítí je oblíbená metoda pro numerické řešení parciálních diferenciálních rovnic, které se v aplikacích transportů tepla a hmoty vyskytují. Podstata metody záleží v aproximaci základní diferenciální rovnice s příslušnými okrajovými podmínkami odpovídající rovnicí diferenční, jež má tvar soustavy algebraických rovnic. Aproximace je tím dokonalejší, čím přesnějšími výrazy nahrazuje derivace. Náhrada se provádí v diskrétních místech tvořených uzly sítě, které pokrývají zkoumanou oblast. Metoda sítí má značné potenciální možnosti uplatnění v podmínkách počítačového modelování, neboť je pro ni charakteristická opakovatelnost jednoduchých algebraických operací. Při praktickém řešení je třeba se zaměřit na konvergenci, přesnost a stabilitu řešení. 1.1 Aplikace metody sítí pro 3D řešení teplotního pole kruhového předlitku Při řešení 3D teplotního pole metodou sítí je síť na předlitku generována způsobem, uvedeným na obr. 1. U kruhového předlitku se vyskytují tři typy elementů, a to elementy vnitřní, obvodové (vnější) a osové. 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ří obvodové elementy, jejichž šířka je ∆r/2. Při výpočtu teplotního pole je podstatný správný výběr podmínek jednoznačnosti řešení, tedy podmínek geometrických, fyzikálních, počátečních a povrchových. Geometrická podmínka je známa, neboť plyne z konkrétního tvaru odlévaného předlitku, počáteční podmínka pak souvisí s teplotou oceli v mezipánvi. Fyzikální podmínky zahrnují součinitel tepelné vodivosti, měrnou tepelnou kapacitu a hustotu oceli. V programu se využily poznatky G. Wölka [2], kde je na základě experimentálních měření odvozena závislost výše uvedených veličin na teplotě a chemickém složení oceli ve tvaru polynomu třetího stupně. Pro konkrétní veličinu P platí rovnice

2

METAL 2008 13. –15. 5. 2008, Hradec nad Moravicí ___________________________________________________________________________ 3

3

P( A, t ) = ∑∑ bi , j ⋅ At ⋅ t j

(1)

i =0 j =0

kde

A je procento legujících prvků (%), bi , j - konstanty závislé na složení oceli,

t

-

teplota (°C).

Mezi fyzikální podmínky pak dále patří teplota likvidu a solidu, které definují oblast tekuté a tuhé fáze a také vnitřní tepelný objemový zdroj, který je v průběhu tuhnutí předlitku představován uvolňováním latentního tepla tuhnutí oceli. Proces přestupu tepla u vnějších a obvodových uzlů je dále ovlivněn i typem povrchové podmínky. Ta může být v případě simulace teplotního pole ∆r kruhového předlitku ZPO dvojího typu. Podmínka II. druhu (Neumannova), při které je známa hodnota hustoty ∆ r/2 tepelného toku q. Tato situace je charakteristická pro odvod tepla v krystalizátoru. Pro zjištění hustoty tepelného toku se vychází z měření teplotních gradientů v různých místech po výšce a obvodu pracovního povrchu ∆ z/2 krystalizátorové vložky, což dovoluje vypočítat místní hodnoty tepelných toků. ∆z Další možnost, jak stanovit hodnoty tepelného toku vychází z měření parametrů chladicí vody krystalizátoru. Nevýhodou je zjištění pouze střední hodnoty tepelného toku. Jako optimální postup se jeví kombinace obou výše uvedených postupů. Podmínka III. druhu (Fourierova) se aplikuje v případech kdy je známa hodnota součinitele přestupu tepla α – typicky pro sekundární a terciární zónu Obr. 1. Síť kruhového předlitku chlazení. V sekundární zóně chlazení se Fig. 1. Net of the round blank pro určení hodnoty součinitele přestupu tepla ostřikem využívá laboratorního výzkumu na teplém fyzikálním modelu na katedře tepelné techniky VŠB-TU Ostrava. 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 se řeší z kriteriální rovnice, udávající závislost mezi kritériem Nusseltovým, Prandtlovým a Grashofovým. Teplo je v této zóně odváděno z předlitku převážně sáláním do okolního vzduchu (kolem 90 %) zbytek tvoří přirozená konvekce (10 %). Pro výpočet 3D teplotního pole kruhového předlitku byla použita rovnice Fourierova – Kirchhoffova, která zohledňuje rychlost pohybujícího se předlitku. Při zanedbání ostatních složek rychlosti kromě rychlosti vz má rovnice tvar ∆φ

3

METAL 2008 13. –15. 5. 2008, Hradec nad Moravicí ___________________________________________________________________________

ρ⋅

∂i = − div (− λ ⋅ grad t + ρ ⋅ vz ⋅ i ) + qV ∂τ

kde

(W.m-3)

(2)

ρ je hustota (kg.m-3), i

τ λ t vz qV

-

měrná entalpie (J.kg-1), čas (s), součinitel tepelné vodivosti (W.m-1.K-1), teplota (°C), rychlost lití ve směru osy z (m.s-1), vydatnost vnitřního objemového tepelného zdroje (W.m-3).

Pokud jsou vnitřní tepelné objemové zdroje nulové a příslušné derivace v rovnici (2) se nahradí konečnými rozdíly, vznikne výraz  6  Pi ∑   ∆τ i=1 ⋅ + S ⋅ vz ⋅ (i6 − i0 ) i0+ − i0 = ∆V  ρ   

kde

(J.kg-1)

(3)

i0 je měrná entalpie v čase τ (J.kg-1), i0+ Pi ∆τ ∆V S

-

měrná entalpie v čase τ + ∆τ (J.kg-1), tepelné toky v jednotlivých směrech souřadného systému (W), časový krok (s), elementární objem (m3), elementární plocha (m2).

První člen na pravé straně rovnice (3) vyjadřuje konduktivní a druhý člen konvekční složku. Vzorec (3) lze považovat za univerzální pro všechny typy uzlů, pouze elementární objemy, elementární plochy a hodnoty tepelných toků budou nabývat pro každý typ elementu a povrchové podmínky odlišných hodnot. Pro vnitřní uzlové body lze elementární objem ∆V a elementární plochu S vyjádřit následujícími vztahy: ∆V = ∆ϕ ⋅ ∆r ⋅ ∆z ⋅ r

(m3)

(4)

kde ∆ϕ je elementární úhel (rad), ∆r - vzdálenost mezi dvěma síťovými uzly ve směru poloměru r (m), ∆z - vzdálenost mezi dvěma síťovými uzly ve směru osy z (m). S = ∆ϕ ⋅ ∆r ⋅ r

(m2)

(5)

Po dosazení jednotlivých tepelných toků reprezentující konduktivní složku a po rozepsání konvekční složky lze rovnici (3) s přihlédnutím k rovnicím (4) a (5) zapsat ve tvaru:

4

METAL 2008 13. –15. 5. 2008, Hradec nad Moravicí ___________________________________________________________________________ ∆τ ⋅ λ ⋅ ∆ϕ ⋅ ∆z ∆τ ⋅ λ ⋅ ∆r ⋅ ∆z ⋅ (i1 + i2 − 2 ⋅ i0 ) + ⋅ (i3 − i0 ) r + ∆r  ∆ϕ  2 ⋅ r ⋅ sin ln  r  2  i0+ − i0 = + c p ⋅ ρ ⋅ ∆z ⋅ ∆r ⋅ r ⋅ ∆ϕ

(J.kg-1)

(6)

∆τ ⋅ λ ⋅ ∆ϕ ⋅ ∆z ∆τ ⋅ λ ⋅ r ⋅ ∆ϕ ⋅ ∆r ⋅ (i4 − i0 ) + ⋅ (i5 + i6 − 2 ⋅ i0 ) r ∆z ln ∆ϕ ⋅ r ⋅ ∆r ⋅ v z ⋅ ∆τ r − ∆r + + ⋅ (i6 − i0 ) c p ⋅ ρ ⋅ ∆z ⋅ ∆r ⋅ r ⋅ ∆ϕ ∆z ⋅ ∆r ⋅ r ⋅ ∆ϕ

kde

i1 − i6 cp

jsou je

měrné entalpie v jednotlivých směrech souřadného systému (J.kg-1), měrná tepelná kapacita při konstantním tlaku (J.kg-1.K-1).

Poněvadž je obecně entalpie funkcí teploty, lze z rovnice (6) přímo zjišťovat hodnoty teplot v jednotlivých uzlových bodech. Rovnice (6) tedy udává celkem jednoduchý postup řešení teplotního pole předlitku. Aby byl numerický postup podle rovnice (6) prakticky použitelný, musí být konvergentní a numericky stabilní [3]. Je tedy nutno dodržet určitou závislost mezi časovým intervalem ∆τ a prostorovým dělením. Nedodržení podmínky stability v jediném bodě znamená divergenci příslušné teploty, která se postupně rozšíří i do všech ostatních uzlů. U praktických úloh se ovšem obvykle vyskytují největší gradienty v ose a na povrchu těles, takže odpovídající podmínky stability jsou pro osové elementy přísnější než pro vnější a vnitřní body. Pro zachování celkové stability řešení je třeba brát v úvahu tu nejpřísnější podmínku, která se v systému vyskytuje. Pro vnitřní element lze odvodit z rovnice (6) následující vztah pro numerickou stabilitu řešení vnitřního uzlu kruhového předlitku:

∆τ ≤

1     λ  1 1 2  vz 1 ⋅ + + + 2 + r + ∆r r ∆z  ∆z cp ⋅ ρ  2  ∆ϕ  ∆r ⋅ r ⋅ ln  ∆r ⋅ r ⋅ ln  r ⋅ ∆ϕ ⋅ sin   r r − ∆r  2   

(s)

(7)

Rozdíl mezi teplotou hodnotou vypočítanou z rovnice (6) a skutečnou teplotou, je dán přesností, s jakou se konečné rozdíly nahradí příslušnými derivacemi v rovnici (2). Je zřejmé, že čím je časový interval a prostorové dělení sítě menší, tím je i menší chyba řešení. Konkrétní velikost chyby lze stanovit rozložením dané funkce v nekonečnou řadu pomocí Taylorova rozvoje [1]. Obdobně se sestavují další typy těchto rovnic reprezentující další typy možných elementů. Pokud je na části předlitku zadána adiabatická hranice, pak příslušné tepelné toky odpadají a rovnice (3) se zjednoduší o tyto chybějící tepelné toky.

5

METAL 2008 13. –15. 5. 2008, Hradec nad Moravicí ___________________________________________________________________________ 2. POČÍTAČOVÁ SIMULACE Samozřejmou činností před případnou realizací numerického modelu v on-line režimu je jeho validace s ohledem na reálné ZPO. Jinými slovy, je nutné vytvořit simulační aplikaci, jež na základě reálných vstupních dat z provozu generuje výstupní veličiny modelu a ty jsou pak následně porovnávány s výstupními veličinami reálného ZPO. Základní algoritmy simulační aplikace jsou pak shodné s těmi, které budou implementovány do aplikace pro on-line režim. Simulační aplikace je navržena tak, aby poskytovala všechen nezbytný uživatelský komfort při konfiguraci modelu a režimu výpočtu, dále nabízí přehlednou prezentaci informací o průběhu simulace a zobrazení s následnou archivací získaných výsledků – viz obrázek 2. Konfiguraci modelu a režimu výpočtu je možné realizovat buď ručně pomocí aplikačního rozhraní, kde je uživateli umožněna modifikace všech konfiguračních parametrů, nebo přečtením jednotlivých parametrů z konfiguračního souboru. Před samotným zahájením simulačního výpočtu je nejprve na základě základních vstupních parametrů, jako je vzdálenost jednotlivých uzlů sítě, délka jednotlivých zón ZPO, rozmístění a typ chladicích trysek, vygenerována výpočetní síť. Následně, zpracováním dalších vstupních údajů, mezi které patří chemické složení oceli, výška hladiny v krystalizátoru, teplota lití, tlak a průtok chladícího média v příslušných zónách atd., jsou definovány parametry výpočetní sítě, okrajové podmínky, nezbytné funkční závislosti a je také vyhodnocena podmínka numerické stability. Jako poslední úkon v rámci preprocesingu je provedeno přidělení typu jednotlivým uzlům výpočetní sítě dle zvoleného režimu výpočtu (počátek lití, jen lití nebo konec lití).

Obr. 2. Hlavní obrazovka programu Fig. 2. Main screen of program

6

METAL 2008 13. –15. 5. 2008, Hradec nad Moravicí ___________________________________________________________________________ V průběhu simulace je samozřejmě možné ručně modifikovat některé vstupní parametry, např. rychlost lití, průtok chladícího média, výšku hladiny v krystalizátoru atd. Model je schopen reagovat na tyto změny v reálném čase. Dynamiku modelu pak uživatel může sledovat jak na jednotlivých grafech, které aplikace poskytuje, tak také na velmi názorném barevném zobrazení teplotního pole chladnoucího předlitku a zobrazení polohy solidu a likvidu. Po skončení simulace je samozřejmostí možnost archivace simulačních výsledků do souboru a nabídka přímého exportu vykreslených grafů. Důležitou vlastností aplikace s ohledem na validaci numerického modelu je možnost přímého načítání reálných provozních dat z datového souboru, který obsahuje hodnoty vstupních parametrů v daných časových okamžicích. Tato vlastnost umožňuje simulovat průběh již proběhlého procesu lití na konkrétním stroji a získat tak představu o reálnosti a přesnosti numerického modelu. 2.1 Vstupní parametry a výsledky simulace Odlévaný formát byl kruhového profilu průměru 530 mm, tedy r = 265 mm, délka primární zóny (krystalizátoru) 600 mm, délka sekundární zóny 4770 mm a délka terciární zóny 29140 mm. Vzdálenosti mezi dvěma sousedními elementy jsou: ∆r = 10 mm a ∆z = 10 mm. Za těchto podmínek je celkový počet bodů ve výpočtové síti 372 708. Latentní teplo 270000 J.kg-1, emisivita se v programu definuje jako funkce teploty povrchu předlitku. Rozmístění chladicích trysek v sekundární oblasti ZPO bylo nadefinováno dle reálných rozměrů uvedených v technických výkresech pro tuto oblast. Pro posouzení rovnoměrnosti ostřiku v závislosti na množství a tlaku chladicí vody existuje na katedře tepelné techniky v Ostravě modelové zařízení, které umožňuje stanovení ostřikových charakteristik trysek, které jsou používány na zařízeních plynulého odlévání. Výstupní parametry modelu, jako je průběh teploty po délce předlitku, teplotní pole kruhového předlitku, tloušťka licí kůry, likvidus, solidu, zachycují následující obrázky.

Obr. 3. Tloušťka licí kůry Fig. 3. Thickness of casting skin

Obr. 4. Teploty kruhového předlitku Fig. 4. Temperatures of round blank

Obr. 5. Likvidus, solidus Fig. 5. Liquidus, solidus 7

METAL 2008 13. –15. 5. 2008, Hradec nad Moravicí ___________________________________________________________________________

Obr. 6. Teplotní pole předlitku Fig. 6. Temperature field of blank Jak již bylo výše zmíněno, při řešení teplotního pole explicitní diferenční metodou je nezbytné zajistit numerickou stabilitu řešení. Během výpočtu se tedy generují jednotlivé „lokální“ stability řešení pro jednotlivé uzly kruhové sítě. Pro „globální“ stabilitu řešení se pak volí ta hodnota časového kroku, která je nejmenší ze všech „lokálních“ hodnot. Pro kruhové předlitky se jedná o středové elementy, následují vnější elementy, u kterých je jejich stabilita výrazně ovlivněna povrchovou podmínkou. Největší časový krok dávají vnitřní elementy sítě, které nejsou povrchovými podmínkami bezprostředně ovlivněny. Řešení teplotního pole je tím časově náročnější, čím se volí jemnější rastr sítě. Proto je vhodné nalézt určitý kompromis mezi numerickou stabilitou řešení a jemností sítě. 3. ZÁVĚR V příspěvku jsou uvedeny základní charakteristiky simulačního PC programu pro výpočet tuhnutí kruhového předlitku explicitní diferenční metodu sítí. Rozebrány jsou jednotlivé podmínky jednoznačnosti, které proces tuhnutí předlitku ovlivňují. Jsou zde diskutovány problémy numerické stability řešení a jemnosti sítě s ohledem na celkový výpočtový čas. Matematický on-line model umožňuje kromě simulace teplotních polí i posouzení dalších technologických parametrů jako např. metalurgická délka, tloušťka licí kůry atd. Programem lze rovněž simulovat vliv doby tuhnutí na tloušťku licí kůry a vliv změny licí rychlosti na metalurgickou délku. Program je schopen načítat externí data (teploty) z jednotlivých kampaní plynulého odlévání ZPO a graficky je porovnávat s vypočítanými hodnotami a proces tuhnutí předlitku tímto verifikovat. Dále umožňuje také grafické znázornění průběhu hustoty, měrné tepelné kapacity, tepelné vodivosti a měrné entalpie oceli v závislosti na teplotě.

LITERATURA [1] PŘÍHODA, M., et al. Nové poznatky z výzkumu plynulého odlévání oceli. 1. vyd. Ostrava: Ediční středisko VŠB-TU Ostrava, 2001. 175 s. ISBN 80-248-0037-3. [2] WÖLK, G. Stahl und Eisen 91, 1971, č. 11, s. 282-286. [3] PŘÍHODA, M., RÉDR, M. Základy tepelné techniky. 1. vyd. SNTL Praha, 1991, 680 s. ISBN 80-03-00366-0. [4] KUNEŠ, J. Modelování tepelných procesů. SNTL Praha, 1989, 424 s. ISBN 80-03-00134-X.

Dílčí část výzkumu probíhala v rámci řešení projektu evidenční číslo 106/07/0938 GA ČR a 102/06/1332 GA ČR.

8