Primerjava lastnosti numeričnih metod kontrolnih prostornin in robnih elementov z dvojno recipročnostjo
Comparison of Control Volume and Dual Reciprocity Boundary Element Numerical Method
A. Košir1, B. Šarler, FS, Univerza v Ljubljani
Prejem rokopisa - received: 1995-10-04; sprejem za objavo - accepted for publication: 1995-12-22
V delu sta uporabljeni numerični metodi: metoda kontrolnih prostornin in metoda robnih elementov z dvojno recipročnostjo, s katerima z iterativno implicitno shemo rešimo nalogo konvekcijsko-difuzijskega prenosa toplote v sistemu s faznim prehodom v poljubno razsežnem prostoru. Temperaturno polje je izračunano iz tako diskretiziranih enačb, da te upoštevajo močno nelinearno zvezo med temperaturo in specifično entalpijo pri temperaturi faznega prehoda. Z Voller-Svvaminathanovim nastavkom pri vsaki iteraciji znotraj časovnega koraka napovemo novi prostorninski delež trdne faze. Točnost obeh numeričnih metod primerjamo z enorazsežnim konvektivno-difuzijskim problemom, za katerega smo našli analitično rešitev. Predstavljeni numerični shemi sta bili pri prikazanih primerih konvergentni.
Ključne besede: prenos toplote, fazni prehod, strjevanje, Štefanova naloga, gibajoči se rob, kontinuirano ulivanje, konvekcijsko-difuzijski pojavi, metoda robnih elementov, metoda dvojne recipročnosti, metoda kontrolnih prostornin
Two numerical fixed grid methods, the control volume method and the boundary element method with dual reciprocity, are presented with an updated iterative implicit scheme to solve arbitrary dimensional phase change problems with convection and diffusion heat transport. The temperature field is deduced from the resolution of the governing equations whose discretization takes into account the highly nonlinear relation between temperature and specific heat at the melting point. At each iteration within one tirne step an updated volume fraction of solid phase is found from the Voller-Swaminathan procedure. The accuracy of the proposed numerical methods has been checked on one dimensional convection-conduction test problems for vvhich i/ve have found an analytical solution. For sever al examples presented, the numerical schemes have demonstrated convergence properties.
Key vvords: heat transfer, phase change, solidification, Štefan problem, moving boundary, continuous casting, convective-diffusive phenomena, boundary element method, dual reciprocity method, control volume method
1 Uvod
Učinkovito reševanje konvekcijsko-difuzijske enačbe navadno zahteva posebne in dobro premišljene numerične sheme. Predmet predhodnih1"2 in pričujoče raziskave, v kateri nadaljujemo razvoj in preskus numerične metode robnih elementov z dvojno recipročnostjo in študij njenih lastnosti, so sistemi s faznim prehodom med kapljevinsko in trdno fazo.
V raziskavi smo na izbranih analitično izračunljivih primerih5 prenosa toplote v eni razsežnosti v snovi, ki se strjuje
' Aleš Košir, dipl.inž.fiz.
Fakulteta za strojništvo
Laboratorij za dinamiko fluidov in tcrmodinamiko
1000 Ljubljana, Aškerčeva 6
ali tali, primerjali numerični metodi kontrolnih prostornin3 in robnih elementov z dvojno recipročnostjo2. Zanimala nas je predvsem stabilnost, točnost, konvergenca s krajšanjem časovnega koraka in večanjem diskretizacije vsake od metod pri različnih začetnih in robnih pogojih ter pri različnih izvedbah numeričnih metod. Fizikalno preproste primere smo izbrali zato, da smo študirali predvsem lastnosti numeričnih metod.
2 Vodilne enačbe
Klasična numerična metoda v uporabi je standardna ental-pijska metoda, ki rešuje diskretizirano energijsko enačbo v en-talpijski obliki:
p ^ + pv • Vh = V • (KVT) ot
(1)
Metoda je preprosta za uporabo in učinkovita ter ohranja energijo sistema, vendar pri nelinearnosti enačbe pri neeksplicitni časovni diskretizaciji zahteva uporabo dodatnih iteracijskih shem znotraj enega časovnega koraka. Med temi se je za posebej učinkovito izkazala shema Vollerja in Swami-nathana4.
Navadno izberemo za spremenljivke v vozelnih točkah diskretiziranega območja vrednosti za entalpije in temperature. Z diskretizacijsko enačbo izračunamo na naslednjem koraku neznane entalpije in iz njih s konstitutivno zvezo T = T (h) vrednosti za temperaturno polje. Entalpija v vozelni točki predstavlja povprečno entalpijo v kontrolni prostornini. Če je ta entalpija med entalpijo trdne in kapljevinske faze, pomeni, da se snov v tej prostornini strjuje. Temperatura v taki kontrolni prostornini je med strjevanjem ves čas enaka temperaturi tališča. Te konstantne vrednosti povzročijo, da se temperature v kontrolnih prostorninah v bližini medfaznega roba s časom spreminjajo stopničasto, a se z oddaljevanjem tudi zgladijo.
Naslednji klasični numerični način je s temperaturno formulacijo energijske enačbe, pri kateri prenesemo latentno toploto strjevanja v člen z navideznimi toplotnimi izvori:
f)T	3f
p Cp— + pcpV • VT = V • (KVT) + Lp + Lpv ot	ot
Vfs (2)
Poleg te enačbe, ki opisuje fizikalni model, pa na medfaznem robu v točki določeni s temperaturo faznega prehoda T = T„„ velja Štefanov pogoj, ki opisuje ohranitev energije:
9T i „ 3T i
K T- Xm-E - K
I Xrn+E :
3xm	3x„
Lp (vm-v) Sgn(T | T | xm+e)	(3)
Tako zapisan Štefanov pogoj velja povsem splošno. Pri tem v,„ označuje hitrost medfaznega roba glede na izbrani opazovalni sistem in v hitrost gibanja snovi. Znano je, da standardna entalpijska metoda in temperaturna formulacija energijske enačbe v obliki (2) že vsebujeta Štefanov medfazni pogoj (3). Tega potrebujemo pri sestavljanju analitičnih rešitev.
3 Analitična rešitev
Omejimo se na model kontinuiranega ulivanja čistih aluminijastih plošč, tako širokih, da dovolj daleč od stranskih robov ni čutiti njihovega vpliva, in tako tankih, da smemo trditi, da ni prečnega temperaturnega gradienta. Toplotne lastnosti trdne in kapljevinske faze naj bodo konstantne, homogene in v obeh fazah enake. Prenos toplote v takem ingotu pri dani stalni hitrosti ulivanja opisuje enorazsežna linearna kon-vekcijsko-difuzijska enačba v brezdimenzijski obliki:
3h 8h 92T — + Pe — = Ste —-3t dx dx2
(4)
Vse snovne lastnosti smo vključili v dve brezdimenzijski števili, Pecletovo in Štefanovo. Štefanov pogoj se brezdimenzi-jsko zapiše v obliki:
I xm-£ " I Xm-£ J =
3T
Ste (—
oxm
dxrr
,dxm . _ i	_ i *
- Pe) Sgn(T I xm+E - T | xm-e;
Pri numerični analizi kontinuiranega ulivanja iščemo rešitev enačbe (4) po dovolj dolgem času. ko so izzveneli prehodni pojavi in se je stanje ustalilo. Cilj naloge je poiskati z izotermo faznega prehoda določeno lego medfaznega roba in določiti temperaturno polje obeh faz. Rešitve nalog s faznim prehodom so sestavljene iz rešitev nalog brez faznega prehoda. Rešitve so vedno monotone funkcije.
Pri Dirichletovih robnih pogojih na obeh krajiščih intervala pči, x2], f(xi) = Ti in f(x2) = T2 in pri procesnih parametrih Pe * 0 je rešitev v enem območju enaka:
T,(x) = Ci exp(Pe x) + D,	(6)
in v drugem
T2(x) = C2 exp(Pe x) + D2	(7)
Rešitvi morata ustrezati robnima pogojema:
T, =T1(x,) = C, exp(Pexi) + D,	(8)
T2 = T2(x2) = C2 exp(Pe x2) + D2	(9)
temperaturnemu medfaznemu pogoju:
Tm = T,(xm) = C[ exp(Pe xm) + D,	(10)
Tm = T2(xm) = C2 exp(Pe xm) + D2	(11)
in Štefanovemu medfaznemu pogoju:
Ste (C, - C2) exp(Pe xm) = - Sgn(T2-Ti)	(12)
Vpeljimo naslednje oznake:
Ei=exp(Pexi)	(13)
E2 = exp(Pe x2)	(14)
E = exp(Pe xm)	(15)
sgn = Sgn(T2-T,)	(16)
Iz temperaturnega pogoja na medfaznem robu izrazimo konstanti D/ in D2:
Di = Tm - C|E	(17)
D2 = Tm - C2E	(18)
Iz robnih pogojev izrazimo konstantni Ci in C2: ti-Tm
C2 =
Ei-E T2-Tm E2-E
(19)
(20)
Lega medfaznega roba je določena z rešitvijo kvadratne enačbe:
{ Ste(T 2-T i )+sgn }E2+{ Ste[(T2-T i )E2-E, ]-
-sgn(El+E2)}E+sgnEiE2 = 0	(21)
v kateri označimo koeficiente
aE2 + bE + c = 0	(22)
tako da se rešitvi izrazita po znanem obrazcu:

1 , -b±Vb2-4ac
2a
(23)
Korena izračunamo tako, da prvemu ustreza znak plus in drugemu znak minus. Nekatere značilne rešitve pri karakterističnih vrednostih parametrov so predstavljene v tabeli 1.
Fazni prehod poteka izotermno pri temperaturi Tm
-10-
Slika 1: Lega medfaznega roba v odvisnosti od Pe in Ste pri referenčnem primeru xi = 0, X2 = 1. Ti = 1, T2 = 0, Tm = 0,95 Figure 1: Interphase boundary position, dependent on Pe and Ste in reference test čase xi = 0, X2 = 1, Ti = 1. T2 = 0, Tm = 0,95
Tabela 1: Odvisnost lege medfaznega roba od parametrov. Pri temperaturnem polju, ki raste s koordinato x, je fizikalno smiselen prvi koren xm.i kvadratne enačbe (21), sicer pa drugi xm.2. V tabeli
sta navedena fizikalno nesmiselna korena samo pri prvih treh primerih. Zadnji primer kaže, da s to analitično rešitvijo ne moremo
zajeti primerov z ničelno hitrostjo Table X: Dependence of interphase boundary position on parameters.
In temperature fields which increase vvith x coordinate, the physically correct solution of the quadratic equation (21) is xm.i, otherwise it is Xm.2- The table shows physically impossible roots only in the first three places. The final exainple shovvs that with this solution we cannot take into account examples with zero velocity
XI	X2	T,	T2	Tm	Ste	Pe	Xra.l	Xm.2
0	1	1	0	0,95	1	2	-0,028234	0,681661
1	0	1	0	0,95	1	2	-0,334824	0.988250
0	1	0	1	0,95	1	2	0.988250	-0,334824
0	1	1	0	0,5	1	2	0,864614	
0	1	1	0	0,95	10	2	0,209455	
0	1	1	0	0,95	0,1	2	0,954876	
0	1	1	0	0,95	0	2	1	
0	1	1	0	0,95	1	10	0,935564	
0	1	1	0	0,95	1	1	0,398643	
0	1	1	0	0,95	1	0,1	0,058178	
0	1	1	0	0,95	1	0	/	/
4 Numerični shemi
4.1 Metoda kontrolnih prostornin
Pri numeričnem reševanju smo uporabili metodo kontrolnih prostornin na ekvidistantni mreži, predstavljeni na sliki 2, primer a. Za diskretizacijo konvektivnega člena smo uporabili privetrno shemo in za diskretizacijo v času delno ali popolnoma implicitno shemo. Znotraj enega časovnega koraka se z iteracijo približujemo vrednostim za temperaturno polje in za prostorninski delež trdne faze pri naslednjem časovnem koraku, dokler ni razlika med dvema zaporednima iteraci-jskima približkoma dovolj majhna. Za napovedovanje pros-torninskega deleža trdne faze v enačbi (2) smo uporabili Voller-Swaminathanovo shemo
^fl+ali-l
f+Ali — jrt+Ati-1 + _5_ ('p+At i __ 'pi+Al i-1^
3T
Natančnost rezultatov numerične metode ocenimo z napako, ki jo izračunamo kot razliko med numeričnimi in analitičnimi vrednostmi za temperaturno polje v vozelnih točkah. V tabeli 2 so predstavljene vrednosti za povprečno absolutno napako ea, za koren povprečnega kvadrata napake eR in za največjo napako eM-
4.2 Metoda robnih elementov z dvojno recipročnostjo
Kot pri metodi kontrolnih prostornin smo tudi v tem primeru uporabili temperaturno formulacijo, izbrali stalno mrežo, predstavljeno na sliki 2, primeri b, c in č, in diskretizi-rali navidezno enorazsežno območje z geometrijsko konstantnimi robnimi elementi s konstantnimi interpolacijskimi funkcijami na njih. Na stranskih robovih smo predpisali homogene Neumannove robne pogoje. Temperaturno polje in njegov gradient smo v času diskretizirali z delno implicitno metodo in kot fundamentalno rešitev uporabili rešitev Laplaceove enačbe v dveh razsežnostih. Za radialne interpolacijske funkcije metode dvojne recipročnosti smo izbrali funkcije 1 + r. Z njimi prenesemo vpliv notranjih točk na rob območja in se izognemo integraciji po notranjosti območja. Integrali po robnih elementih so bili izračunani z numerično integracijo z 32 Gaussovimi točkami. Kot pri metodi kontrolnih prostornin smo tudi tu za napovedovanje prostorninskega deleža trdne faze uporabili Vol-ler-Swaminathanovo shemo. Metoda je obširneje predstavljena v4
Stacionarni analitični primer je z obema numeričnima metodama rešen s približevanjem stacionarnemu stanju ob začetnem pogoju
T(x, 0) = 1
V tabeli 3 so predstavljene vrednosti za napako numerične metode.
5 Rezultati
Rezultati so predstavljeni v obliki tabel 2 in 3 z vrednostmi za povprečno absolutno napako ea, za koren povprečnega kvadrata napake eR in za največjo napako eM v odvisnosti od Pecletovega in Štefanovega števila.
Tabela 2: Napake med analitično in numerično rešitvijo z metodo kontrolnih prostornin. Črka a označuje mrežo in v oklepaju število
točk znotraj območja v vrsti Table 2: Error between analytical and numerical solutions using the control volume method. Letter a stands for the grid and the number in brackets for the number of points in a single row
Ste		ea	eR	eM
Pe = 1				
0,2 1	a(14)	0.582E-2	0.627E-2	0,712E-2
03	a(14)	0.567E-2	0.666E-2	0.702E-2
Pe = 2				
1	a(14)	0.587E-2	0.596E-2	0.633E-2
Pe = 5				
0,2	a(14)	0,333E-1	0.473E-1	0,832E-1
1	a(14)	0.288E-1	0.322E-1	0.542E-1
OO	a(14)	0.172E-1	0.288E-1	0.357E-1
6 Sklep
Mreža a(9)
0.1 0.05
o
0.1 0.05
Mreža b(9) 0
Mreža c(9)
				
■ s.	0	loto-	lO+OtO+O+O+O+OI	s -
				
OIOIOIOIOIOIOIOIOIO
OIOIOIOIOIOIOIOIOIO
0.5
Z opisanim modelom smo ugotovili, kakšen je vpliv treh brezdimenzijskih števil na lego medfaznega roba pri izbranih stalnih vrednostih za preostale parametre naloge. Brž ko Štefanovo število zraste preko 3, se njegov vpliv na lego medfaznega roba ustali. Pecletovo število najmočneje vpliva na lego medfaznega roba na intervalu | Pe | < 1. Če je Pe > 10 ali Pe < -1, je pri referenčnem primeru medfazni rob tik ob robu območja.
Preglednici napak numeričnih metod 2 in 3 kažeta, kako se napaka povečuje z večanjem Pecletovega števila in s približevanjem Štefanovega števila proti nič. Metoda robnih elementov kaže pri enakih razmerah bolj stabilno vedenje in daje bolj točne rezultate, posebej pri večjih Pecletovih številih.
Tabela 4: Pregled definicij za brezdimenzijska števila Table 4: Overview of definitions of nondimensional numbers
Mreža č(9)
0.2 0.1 o
OIOIOIOIOIOIOIOIOIO
OIOIOIOIOIOIOIOIOIO
Slika 2: Uporabljene diskretizacije območja z devetimi točkami v
vrsti znotraj območja Figure 2: Discretization of domains with nine internal points in one row
X
X = —
a
t = ■
Kt
pcPa~ T-min(Ti,T2)
h =
T2-T11 h
Ste = Pe
cP IT2-T11
L
pcPav
K
Tabela 3: Napake med analitično in numerično rešitvijo z BEM DRM. Črke b, c in č označujejo mrežo in v oklepaju število točk znotraj območja v vrsti Table 3: Error between analytical and numerical solutions using the boundary element method vvith dual reciprocity. Letters b, c and č stand for the grid and the number in brackets for the number of points in a single row
Ste		ea	eR	eM
Pe = 1				
0,2	c(14)	0,498E-2	0,45 5E-2	0,593E-2
0,5	c(9)	0,402E-2	0,432E-2	0,551E-2
1	c(14)	0,310E-2	0,340E-2	0,430E-2
OO	c(14)	0.286E-2	0,310E-2	0.407E-2
Pe = 2				
0,2	c(14)	0.548E-2	0,578E-2	0,669E-2
0,5	c(14)	0.327E-2	0,354E-2	0,460E-2
	b(9)	0.423E-2	0,459E-2	0,575E-2
	c(9)	0,413E-2	0,443E-2	0,563E-2
	č(9)	0,420E-2	0,451E-2	0,572E-2
1	c(14)	0,315E-2	0,341E-2	0,445E-2
OO	c(14)	0,294E-2	0.319E-2	0.418E-2
	_ c(9)	0,384E-2	0.413E-2	0,530E-2
Pe = 5				
0,2	c(14)	0,877E-2	0.103E-1	0,171E-1
	c(9)	0,630E-2	0.968E-2	0,234E-1
1	c(14)	0,582E-2	0,637E-2	0.727E-2
OO	c(14)	0.535E-2	0,566E-2	0.668E-2
Tabela S: Pregled uporabljenih oznak. Znak pomeni, daje količina pod njim brezdimenzijska Table 5: Symbols used. The symbo! ~ denotes a non-dimensional quantity
oznaka	enota_opis_
a	m	polovična debelina plošče
C,	C1.2 proste konstante cp	J/kg K specifična toplota
D,	D 1,2	proste konstante
E,	Ei,2	pomožne oznake ea, eR, eM napake
E	m	majhen odmik
fs	1	prostorninski delež trdne faze
h	J/kg	specifična entalpija
i	števec iteracij znotraj časovnega
koraka
K	W/m K koeficient toplotne prevodnosti
L	J/kg	specifična talilna toplota
A.i.2	korena karakterističnega polinoma
diferencialne enačbe
p	kg/m3 gostota
Sgn	1	funkcija predznaka
sgn	pomožna oznaka
T	K	temperatura
Ti, Ti	K	temperatura na levem oziroma
desnem robu
Ta	K	temperatura okolice
Tm	K	temperatura tališča
t	s	čas
At	s	velikost časovnega koraka
v	m/s	hitrost vlečenja ingota
Vm	m/s	hitrost medfaznega roba
X	m	koordinata
Xm	m	lega medfaznega roba
Xml.2	m	korena kvadratne enačbe za
medfaznega roba
Zahvala
Delo je nastalo v okviru projektov Dvoelementno modeliranje trdno-kapljevinskih sistemov, podprtega s strani MZT, in Modeling in material science and processing v okviru COST-512.
7 Literatura
1	Šarler, B., Košir, A.: Solution of melting and solidification problems by the dual reciprocity boundary element method, Lewis, R. W. (ed.): Numerical Methods in Thermal Problems, Vol VIII., Pineridge Press, Swansea 1993, 139-150
2	Šarler, B., Košir, A.: Solution of conduction-convection heat transport in systems with solid-liquid phase change by dual reciprocity boundary element method, in: Lewis, R. W. (ed.): Numerical Methods in Thermal Problems, Vol IX., Pineridge Press, Svvansea 1995, 737-748
3	Patankar, S. V.: Numerical heat transfer and fluidflow, Hemisphere, New York, 1980
4VoIler, V. R., Swaminathan, C. R.: General source-based method for solidification phase change, Num. Heat Transfer, 19B, 1991, 175-189
5Pardo, E., Weckman. D. C.: A fixed grid finite element technique for modelling phase change in steady-state conduction-advection problems, Int. J. Num. Met. Eng., 29, 1990, 969-984