Modeliranje strjevanja pri kontinuiranem ulivanju z dvojno recipročno robno integralsko metodo Modelling of continuous casting solidification by boundary element method with dual reciprocity A.Košir, B.Sarler, Institut "Jožef Štefan", Univerza v Ljubljani, Jamova 39, Ljubljana . V članku pokažemo, kako uporabimo novo numerično metodo robnih elementov z dvojno recipročnostjo pri simulaciji strjevanja kovine ob pogojih, kakršni nastopajo pri kontinuiranem ulivanju. Kontinuirano ulivanje fizikalno opišemo s prenosom energije v nestisljivi snovi s fazno spremembo. Toplotni tok v snovi definira Fourierov zakon, pri čemer upoštevamo temperaturno odvisne lastnosti snovi. Parcialno energijsko transportno-difuzijsko enačbo za dvofazno zmes pomnožimo z Greenovo funkcijo in jo integriramo po območju. Zaradi temperaturno odvisne toplotne prevodnosti temperaturo transformiramo po Kirchhoffu. Z dvakratno uporabo recipročnosti transformiramo integrale po območju na njegov rob. Območje in njegov rob diskretiziramo. Polje na robu območja opišemo z odsekoma ravnimi in konstantnimi robnimi elementi. Polje na območju predstavimo s prostorskimi zlepki 1+r in z linearnimi časovnimi zlepki. Časovna integracija sledi novi shemi, ki sta jo predlagala Voller in Swaminathan. Analizirali smo občutljivost numeričnih rešitev za strjevanje neskončnega dvodimenzionalnega vogala in za enofazno enodimenzionalno klasično Štefanovo nalogo v odvisnosti od prostorske in časovne diskretizacije, Štefanovega števila in širine talilnega intervala. Ključne besede: kontinuirano ulivanje, metoda robnih elementov, metoda dvojne recipročnosti, Štefanova naloga, prenos toplote in snovi. This paper describes the application of the new dual reciprocity boundary element methodfor the energy transport in melting and solidification of metal in continuous casting process with realistic operational parameter s. The continuous casting is de ser i bed with energy transport in systems composed of incompressible phase-change material, with Fourier constitutive heat flux relation, and includes temperature dependent material properties. The diseretization of the Kirchhoffs transformed and heat source term reformulated governing equation is struetured by the Green's funetion of the Laplace equation, the 1+r space spline dual reciprocity boundary-only representation of the domain integrals, and straight line boundary elements with constant space and linear Ume splines. The timestep iterations follow the new Voller and Swaminathan scheme. The sensitivity of the results with respect to space-time diseretization, Štefan number, and melting interval was investigated on a two-phase analytičal solution for the solidification of an infinite reetangular corner and on the classical one-phase Štefan problem. Key words:continuous casting, boundary element method, dual reciprocity method, Štefan problem, heat and mass transfer. 1 Uvod enaka Raziskovanje sistemov s tekočo in trdno fazo običajno zahteva interdisciplinarno teoretsko, eksperimentalno in numerično modeliranje faznih prehodov, mehanike trdnin in trasportnih pojavov. Njegov vpliv na številne bazične in aplikativne raziskave je zelo širok, kar kaže tudi izčrpen zbornik1, opremljen z referencami in s ključnimi besedami. Odkar je Wrobel2 pokazal, daje numerična metoda robnih elementov (BEM) primerna tudi za diskretno aproksimativno reševanje nelinearnih transportnih enačb, se to metodo skuša uporabiti tudi pri problemih s taljenjem in strjevanjem. 2 Formulacija problema Vzemimo fiksirano povezano območje Q, ki ga omejuje rob T. Na tem območju naj bo snov s stalno gostoto po, temperaturno odvisnima toplotno prevodnostjo k? in specifično toploto c>, kjer T = S označuje trdno fazo in T=L tekočo fazo. Spremembo entalpije opisuje konvektivno-difuzijska enačba za kontinuumsko zmes F = -k V7" ■p rp d dt if, Po + fL Po HJ + + V-(/(PoF //J+/xPoK //J = -V.(/5F5+/£FJ+/5 qs+fLqL (D Konvektivno-difuzijsko enačbo (1) lahko zapišemo v obliki dT Poc-+Po dt (fsk-^. dT V, /A"— dT ■VT (4) / \ dfr dT -.v-{k Vr) + q- P0H Krajevno in časovno odvisni koeficient p, t) predstavlja volumski delež določene faze in je normiran z enačbo fs+fc,-1. Specifična entalpija I It za vsako izmed faz je enaka Hs=c5{TH)T„+\Tct{0)d6, (2) •"h H, = cl(Th)Th + fcAe) de+Hi,, (3) kjer je Th poljubna referenčna temperatura in Hh specifična talilna entalpija. Hitrosti delcev vsake izmed faz Vj> ( p, / ) sta v modelu čistega enoosneega stijevanja enaki ( p, t) = \l ( p, t )= V ( p, t) in različni v modelu čistega stebričastega stijevanja, pri čemer vzamemo, da trdna faza v opazovanem sistemu miruje ali se giblje s stalno hitrostjo, Vx ( p, t )= V ( p, t ), V5 ( P, t) = Vsys. Toplotna tokova v vsaki izmed faz sta odvisna od različnih toplotnih prevodnosti v vsaki izmed faz in sta po Fourieru dT dt pri čemer smo vpeljali oznaki za intenzivnost toplotnih izvorov in za razliko specifičnih entalpij med fazama q=fs;(P)=Ž|P_P"1 (27) T'(Pm^+A/) + Tr'-VT'(pmJ<0+A/) = FXpm>'o) + TVva-(pm,/0) + +qr'?(pm>'„+A/)+q;;9(pm>;0), (28) ki jo moramo pred reševanjem preurediti ustrezno robnim pogojem. Zgornji indeks i označuje vrednost spremenljivke pri /'-ti iteraciji. Matrični elementi so določeni z enačbami -T' Poco 1 t„+il*i + d T -T dl ), (29) At k. A t c'( Q, s)k0 At Učinkovitost transformacije (24,25) je močno odvisna od izbire prostorskih zlepkov <\ip„. Po predlogu Partridga et al.7 smo izbrali za prostorske zlepke funkcije s prvim redom razvoja, = 1. Numerični rešitveni postopek za nelinearno integralsko enačbo (14) sam po sebi zahteva iteracijo po časovnem koraku. Pri časovni iteraciji upoštevamo nelinearna člena A and T po nedavno razviti robustni in točni shemi Vollerja in Svvaminathana8. Bistvene prilagoditve te sheme za dvojno recipročno metodo vsebuje enačba (34). Osnovno enačbo (14) diskretiziramo z linearnimi časovnimi zlepki prek časovnega intervala [to, to + Ar]. Rob diskretiziramo z A'r robnimi odsekoma ravnimi elementi Fk z odsekoma konstantnimi vrednostmi. Prvih Nt točk polja pn, ki jih upoštevamo v zlepkih (27), sovpada z geometrijskimi središči robnih elementov, naslednjih Nq točk pa je lahko poljubno posejanih po notranjosti območja Q. Enačbo (14) rešimo, tako da konstruiramo sistem /=1,2,...^ algebraičnih enačb. Te enačbe določa diskretizirana enačba (14), kjer točka izvora s sovpada s točko polja pn. Izvedeni sistem algebraičnih enačb zapišemo v simbolični obliki p0c0-Y[r(p„,f0)]+-{ c=<&=v.fs^v«] -1 Ar (31) (32) Vrednost Ht0 + At m, kjer E označuje bodisi A bodisi T, pri (j +l)-iti iteraciji izračunamo iz prejšnjih iteracij i in (/' -1) z enačbo ^/rt+A/m (33) / \i+l dE -T ydT J V -//0+A/m dT (34) 7'fp»,'0 + Ar)-T,J(pjji,/0 +A/)j. časovno iteracijo ustavimo, ko absolutna vrednost razlike KirchhofTove spremenljivke v dveh zaporednih korakih v vseh mrežnih točkah pm ne presega vnaprej določenega pozitivnega števila Tg. 4 Numerični primeri Da bi preverili obnašanje rešitvenega postopka, smo si izbrali dva različna fizikalna testna primera. Prvi je standardni dvodimenzionalni testni pnmer strjevanja polneskončnega pravokotnega vogala, drugi je klasični enodimenzionalni enofazni Štefanov problem. Rešitev prvega problema smo primeijali s polanalitično rešitvijo9 in rešitev drugega problema z analitično rešitvijo10. Računsko območje predstavlja kvadrat 0 7 Numerično izračunane temperature Tc&\ = H/n(p)vt/»(?) smo primerjali z Noom = 10201 temperaturo, izračunano s primerjalno metodo 9- 10 na enakomerni mreži pcom , kjer so točke postavljene v sečišča premic v = konst. in x = konst. , kot kaže slika 1. Največja absolutna napaka 7max m povprečna absolutna napaka Tave numerične rešitve sta definirani z enačbama = max|rjptom„>0-7;,a(pc, n = 1,2,..., JV , j Ncom = -f\T ,(p ,t)-T fp ,n; / calVrcomn' ' aiia v* comn' /|7 (36) (37) N -i com Mreža Nr Nq N lM 12 16 28 2M 32 64 96 48 144 192 Korak M 'A t 0.01 2A t 0.005 3A t 0.001 Mejna razlika temperatur 7'a v dveh zaporednih iteracijah je enaka 0.001. Tabela 1: Vrednosti diskretizacijskih parametrov. Slika t: Razporeditev mrežnih točk v mreži M Robne točke so označene s krožnicami 0 in notranje točke s j*lnimi krogi • . Mreži *M in 3A/ sta sestavljeni iz podobno razporejenih, vendar različno gosto posejanih točk. Figure 1: Position of mesh points of 2M mesh. Boundary nodes are rep-resented by empty circles 0 and inner nodes by full circles •. Meshes 'M in 3M are similar, but vvith different number of points. 5 Zaključek V članku je predstavljen eden izmed prvih poskusov, kako numerično rešiti večdimenzionalni problem s premičnim robom in s prenosom energije, kjer se integracija skrči le na konstantne fizikalne količine z roba območja. Pokažemo, da je uporabljena metoda primerna pri Štefanovih nalogah obravnavanega tipa. Glavna prednost metode je preprostost uporabe različnih tipov robnih pogojev, preprosta generacija mreže, zmanjšanje števila spremenljivk in možnost uporabe pri problemih z geometrijskim premičnim robom. Glavna slabost metode je velikost končnega algebraičnega sistema enačb, ker moramo upoštevati tudi notranje točke. To slabost lahko odpravimo z učinkovito tehniko podstrukturiranja v kombinaciji z adaptivno strategijo, od katerih sta obe v razvoju. 6 Oznake Latinske črke c, co, C T specifična toplota, konstantni, spremenjivi del f volumski delež faze F toplotni tok h toplotna prestopnost Ii entalpija Hls razlika entalpij specifična talilna toplota toplotna prevodnost, konstantni, spremenjivi del število robnih, notranjih točk točka polja, izvora čas, začetni čas časovni korak temperatura Greenova funkcija začetni pogoj temperatura na robu temperatura solidus temperatura taljenja tempratura likvidus talilni interval hitrost moč toplotnega izvora skalama, vektorska funkcija matrični element konstantna gostota rob A-ti robni element območje konvektivni člen Kirchoffova spremenljivka izvorni člen A ali Y prostorski zlepek časovni zlepek integralski izraz Spodnji indeksi faza trdna, tekoča faza iteracijski indeks sumacijski indeks povprečje izračunano primerjalno maksimum sistem Zgornji indeksi Tabela 2: Testni primer II. Absolutna napaka numerično izračunanih temperatur na mreži 2M,3ATob času t = 0.1 za različne latentne toplote in različne širine talilnega intervala. Specifične talilne toplote so zapovrstjo ' Hlt>2Hh in enake °-25> °-50 oziroma 1.00. Prvi dve vrstici sta bili izračunani pri širini talilnega intervala A7W = 0.01 in zadnji dve pri A7W= 0.001. specifična talilna toplota napaka 1 ZJO "