ZESZY TY N A U K O W E PO L IT EC H N IK I ŚLĄ SK IEJ Seria: E N ER G ET Y K A z. 133
2001 N r kol. 1486
Stanisław K U C Y PER A P olitechnika Śląska, G liw ice
Z A S T O S O W A N I E M E T O D Y B I L A N S Ó W E L E M E N T A R N Y C H I R A C H U N K U W Y R Ó W N A W C Z E G O D O R O Z W I Ą Z Y W A N I A O D W R O T N Y C H Z A G A D N I E Ń W S P Ó Ł C Z Y N N I K O W Y C H P R Z E W O D Z E N I A C I E P Ł A
Streszczenie. W pracy zastosowano metodę bilansów elem entarnych (M BE) i rachunek w yrów nawczy, nazyw any często uogólnioną m etodą najmniejszych kw adratów (UM NK) do rów noczesnego w yznaczania współczynnika przewodzenia ciepła X i ciepła w łaściwego c ciał stałych. Podano krótko charakterystykę (UM NK), założenia m odelu m atem atycz
nego oraz algorytm obliczeń num erycznych. Rozw iązanie zagadnienia w ym agało po
m iaru tem peratury w w ybranych punktach próbki. Dlatego do tego celu w ykorzystano skom puteryzow ane stanow isko badaw cze, na którym była m ożliw ość sterow ania cało
ś c ią pom iarów . P rzedstaw iono przykładow e w yniki badań. W yniki porów nano z dany
m i literaturow ym i.
THE A PPL IC A T IO N OF THE CONTROL VO LUM E M ETH OD A N D THE G E N E R A L SE D L E A ST SQU A R ES M ETHOD TO SOLVE INVERSE PA R A M E TE R PRO BLEM S
S um m ary. In this paper the control volum e m ethod and generalised least squares m ethod for estim ating sim ultaneously the therm al conductivity X and specific heat c have been applied. A short characteristic o f this m ethod, an assum ption o f a m athem ati
cal m odel and num erical algorithm has been given. To solve this problem the tem pera
ture had to be m easured in selected points o f the samples. For this purpose a com puter
ized test stand w as used w ith possibility o f controlling the m easurem ents. The sample results o f the these researches have been presented. The results obtained by this m ethod have been com pared w ith previously published ones.
1. W ST Ę P
D okładna znajom ość w łaściw ości cieplnych m ateriałów , tzn. w spółczynnika przew odzenia ciepła A, ciepła w łaściw ego c je st bardzo istotna z punktu widzenia:
- w ykorzystania m ateriałów o w ym aganych w łaściw ościach cieplnych w urzą
dzeniach, w których jed n y m z w ażnych kryteriów je st oszczędne w ykorzystanie energii cieplnej,
produkcji now ych m ateriałów o w ym aganych w łaściw ościach cieplnych, m atem atycznego m odelow ania różnego rodzaju procesów cieplnych.
D latego obserw uje się szeroki rozwój i stosowanie coraz to bardziej dokładnych metod w yznaczania w łaściw ości cieplnych m ateriałów. W klasycznym podejściu w yznaczanie każ
dego param etru cieplnego w ykonyw ano inną m etodą eksperym entalną. N a przykład ciepło w łaściw e zw ykle w yznaczano m etodam i kalorym etrycznym i, podczas gdy do w yznaczania w spółczynnika przew odzenia ciepła stosowano m etody ustalonego przew odzenia ciepła.
Ostatnio do w yznaczania w łaściw ości cieplnych m ateriałów szeroko stosow ane są m etody nieustalone przew odzenia ciepła. M etody te polegają na połączeniu rozw iązyw ania odw rot
nych zagadnień przew odzenia ciepła z danym i otrzym anym i z pom iarów dla stanu nieustalo
nego. C hociaż m etody te w ym agają bardziej skom plikow anych i droższych urządzeń to w sum ie są one tańsze, poniew aż um ożliw iają w yznaczanie rów nocześnie w ielu param etrów cieplnych. Jednak odw rotne zagadnienia przewodzenia ciepła należą do klasy problem ów źle uw arunkow anych. O znacza to, że ich rozw iązanie je st zw ykle bardzo czułe na niedokładności danych w ejściow ych, tzn. błędów pom iarow ych. W celu ustabilizow ania rozw iązania i zw iększenia dokładności obliczeń m uszą być stosowane specjalne techniki m atem atyczne.
D otychczas najbardziej rozpow szechnionym i m etodam i stabilizacji rozw iązań odw rotnych zagadnień w technice cieplnej była m etoda regularyzacji [6] i m etoda rachunku w yrów naw czego, nazyw ana często uogólnioną m etoda najm niejszych kwadratów. [1,3,4]. W zagad
nieniach cieplnych m etoda rachunku w yrów naw czego była szeroko stosow ana do uzgadnia
nia bilansów m asy i energii [3] oraz do rozw iązyw ania odw rotnych zagadnień brzegow ych i początkow ych przew odzenia ciepła [1,5]. W obecnej pracy zastosow ano j ą do w yznaczenia w spółczynnika przew odzenia ciepła i ciepła w łaściwego m ateriałów stałych. W rzeczyw isto
ści m etoda ta polega na poszukiw aniu m inim um odpowiednio zdefiniowanej form y kw adra
towej oraz uzgadnianiu w ielkości m ierzonych, obarczonych błędam i, które połączone są z w ielkościam i niew iadom ym i za pom ocą tzw. równań warunków. W pracy tej rów nania wa
runków zapisano w ykorzystując odpow iednio sform ułowany model m atem atyczny przepływ u ciepła w próbce. M odel ten sporządzono stosując metodę bilansów elem entarnych. M etoda ta należy do klasy m etod różnicow ych i ma przejrzystą interpretację fizykalną. D latego je st ona bardzo przydatna do m atem atycznego m odelow ania różnych procesów cieplnych. O trzym ane tą m etodą rów nania opisu ją pole tem peratury w próbce, którą nagrzew ano na stanow isku po
Zastosow anie m etody bilansów elem entarnych.. 75
m iarow ym .P oniew aż uogólniona m etoda najm niejszych kw adratów w ym aga większej liczby rów nań niż niew iadom ych, w ięc aby spełnić te w ym agania, rozpatryw ano nieustalony proces nagrzew ania próbki i rów nania w arunku sform ułow ano dla kilku kroków czasowych. N astęp
nie w ykorzystując ww. m etodę rozw iązano odw rotne zagadnienie w spółczynnikow e przew o
dzenia ciepła i w yznaczono param etry cieplne badanej próbki. Przedstaw iono przykładow e w yniki badań.
2. C H A R A K T E R Y ST Y K A Z A ST O SO W A N E J M E TO D Y O B L IC Z E N IO W E J
2.1. Istota m etody rachunku w yrów naw czego i założenia m odelu m atem atycznego
Isto ta m etody p olega n a tym [3], że m ając dany ciąg wielkości x\ (i= l,...,M ), otrzym any z pom iaru badanego procesu, oraz w stępnie oszacow ane w artości ciągu niew iadom ych w ystę
pujących w tym procesie/¡j (j= l,...,N P ), rów nania opisujące proces (tzw. rów nania w arun
ków) m ożna zapisać w następującej postaci:
P oniew aż pow yższe rów nania nie s ą dokładnie spełnione z uw agi n a błędy pom iaru w iel
kości m ierzonych, w ięc celem m etody je st znalezienie takich popraw ek w ielkości m ierzonych 5xi oraz popraw ionych w artości niew iadom ych /?*, aby dokładnie (w sensie przyjętego kryte
rium ) spełniony był następujący układ równań:
Najbardziej w iarygodne wartości poprawek wielkości mierzonych wyznaczane są z wa
runku m inim um ważonej sumy kwadratów tych poprawek, tzn.:
f k ( x \ •■■■• x m • P \ '■■■• Pnp ) ~ 0. dla k r \ , .. ., N Q . (1)
fk (iX
I "*■
SX\)»•••»
(XMd" )*
P\ (2 )(3)
R ów nania te m ożna zapisać w następującej postaci m acierzow ej:
G , - Y' 1 G x =>mi n, (4)
gdzie:
nij - średni błąd bezw zględny i-tej wielkości mierzonej, Ga— m acierz popraw ek w ielkości mierzonych,
V - m acierz zaw ierająca w artości m ] .
Aby m ożna było przeprow adzić proces uzgadniania, m uszą być spełnione dodatkow o nastę
pujące w arunki:
N P < N Q < N P + M oraz N Q < M . (5)
W klasycznej m etodzie rachunku w yrów naw czego problem m inim alizacji rozw iązyw any jest m etodą nieoznaczonych m nożników L agrange’a. [3]. W związku z tym rów nania w arunków m usiały być w prow adzane do algorytm u obliczeniow ego w postaci liniowej.
A X + BY + C = 0 , (6)
gdzie:
A ,B - odpow iednio m acierze w spółczynników przy w ielkościach m ierzonych i niewia- mych,
C - w ektor uw zględniający w arunki brzegowe, X - w ektor zaw ierający w ielkości mierzone, Y - w ektor zaw ierający wielkości w yznaczane.
W prezentow anej pracy do rozw iązania zagadnienia odw rotnego zastosow ano iteracyjną procedurę Levenberga-M arquardta [2], Pozwoliło to uniknąć procesu linearyzacji równań w arunków , a tym sam ym m ożliw ości wprow adzania dodatkow ych błędów oraz dzięki temu algorytm uogólnionej m etody najm niejszych kw adratów stal się bardziej uniw ersalny. W celu osiągnięcia odpow iedniej dokładności obliczeń przy form ułowaniu m odelu m atem atycznego przyjęto następujące założenia:
1. R ozpatruje się proces nagrzew ania nieskończenie rozleglej płyty, tzn. takiej, dla której rozm iary w zdłużne są znacznie w iększe od grubości.
2. M ateriał płyty posiada izotropowe w łaściw ości cieplne.
3. Z akłada się, że przed rozpoczęciem procesu nagrzew ania płyty rozkład tem peratury w całej jej objętości je s t w yrów nany i równy tem peraturze początkowej.
4. W chw ili r = 0 w osi płyty zaczyna działać stały strum ień ciepła q = idem W /m 2.
5. Z akłada się, że obie pow ierzchnie zewnętrzne płyty są dobrze chłodzone.
W ykorzystując pow yższe założenia sform ułow ano model m atem atyczny oraz algorytm obliczeń num erycznych.
Zastosow anie m etody bilansów elem entarnych.. 77
2.2. S form ułow anie m odelu m atem atycznego oraz algorytm u obliczeń num erycznych
D la jednow ym iarow ego przepływ u ciepła zagadnienie początkow o-brzegow e m oże być opisane znanym rów naniem różniczkow ym :
d T d .. d T . CP ~^— = ~z— ( k —— )■
O T O X O X
(7)
W arunki brzegow e:
- dla powierzchni grzanej:
- dla pow ierzchni p rzeciw ległej:
± k 0 7 d x
■ warunek początkowy:
T\ = T = T
l * = ± i w P
T ( x , r = 0) = Tp .
(8)
(9)
( 10)
D la potrzeb uogólnionej m etody najm niejszych kw adratów rów nania przew odzenia ciepła zapisano w postaci różnicow ej przeprow adzając dyskretyzacje obszaru rozw iązania rys. 1.
R y s. 1. G e o m e tria z a g a d n ie n ia g ra n ic z n e g o F ig . 1. G e o m e try o f th e b o u n d a ry p ro b le m
N astępnie do budow y postaci różnicowej m odelu wykorzystano metodę bilansów elem entarnych (M BE) i w rezultacie uzyskano następujący nieliniowy układ równań algebraicznych:
7j*+l = ( l + 6u )-7j* +6,2 -T2 + y,
T**' = b n ■ 7;* + (1 + b22) ■ r 2* + b22 ■ T}
T}*' = bn ■ T* + (1 + ¿>33) • r 3* + bu ■ Ti TtM = b„ - T i + ( l + bt t ) - T i> + biS T s , T r = b„ T : + ( i + b sj - T; ,
gdzie w yraz w olny y, uw zględnia w arunek źródła ciepła na pow ierzchni próbki:
Y\ 2 ■ c ■ p ■ zlx |
N atom iast w spółczynniki b\j określone są zależnościami:
- dla w ęzła „ 1 ” b\y m ożna zapisać:
bu = --- (gdy i = j ) ’ c . p . A x ] - ( A x l + ^ )
bij = --- -t— (gdy i * j ) . c p - A x r ( A x t + - ^ ~ )
A nalogicznie dla w ęzła „5” bij je st określone:
bo = --- k- ^ 1 — ( g d y ' = 7 ) , c - p - A x 5 ■( A xi + - ^ - )
i k - A r
b,j = ---(gdy i*j) . c p - A x 5 - ( A x s + ~ )
D la p ozostałych w ęzłów b,j m a postać:
2 k ■At 1 1 .
bj = — ■ ( - --- — + - ---— ) (gdy i = j ) , c- p - A x: Ax,_t + zlx, A xl + A xi+,
2 - k - A r .
j = --- — — --- — : (gdy I * j ) , c ■ p ■ A x: ■(Axi + Ax¿J
gdzie:
7’*,7’*+1 - tem peratura w i-tym w ęźle w k i k+ 1 kroku czasowym , Ax: - grubość i-tej w arstw y w podziale różnicow ym płyty, z lr - długość kroku czasowego.
( 11)
( 12)
(13)
(14)
( 15)
Zastosow anie m etody bilansów elem entarnych. 79
3. O PIS M E T O D Y P O M IA R O W E J
Schem at stanow iska pom iarow ego przedstaw iono na rys. 2. U kład pom iarow y składał się z cienkiego oporow ego grzejnika elektrycznego, dw óch sym etrycznie rozm ieszczonych w zglę
dem grzejnika w alcow ych próbek oraz czterech płytek m iedzianych z term oparam i. W alcowe
R y s. 2. O g ó ln y s c h e m a t u k ła d u p o m ia ro w e g o Fig. 2. S c h e m e o f th e te s t sta n d
pow ierzchnie próbek zaizolow ane były rozdrobnionym korkiem . Z ew nętrzne pow ierzchnie próbek chłodzone były (utrzym yw ano praw ie stałą ich tem peraturę) w o d ą z ultraterm ostatu.
Do pow ierzchni grzanych i chłodzonych próbek przylegały cienkie płytki m iedziane z przy
m ocow anym i term oparam i. W szystkie term opary były typu N i-NiCr. Do zasilania grzejnika użyto zasilacza sterow anego z m ikrokom putera. N apięcie U i prąd I z grzejnika m ierzone były poprzez kartę p om iarow ą na m ikrokom puterze. Strum ień ciepła dopływ ający do próbek:
Q = P = U ■ I . ( ’ 6)
G ęstość strum ienia ciepła dopływ ającego do każdej z próbek określa zależność:
2
P W
4 ~
na m
.2 ’ 2 (17)Term opary podłączone były poprzez w zm acniacz (ze w zględu na niską wartość sygnałów z term opar) i m ultiplekser (konieczność przełączania odczytyw anych w każdym kroku cza
sow ym sygnałów z term opar) oraz kartę pom iarow ą do m ikrokom putera. Term opary m ierzące tem peraturę pow ierzchni próbek od strony grzejnika oraz pow ierzchni chłodzonych były po
łączone szeregow o. D latego w celu dalszej obróbki danych w zbiorze w ynikow ym zapisyw a
no średnią ary tm ety czn ą w artość w skazań odpow iednich term opar. Do sterow ania całością pom iarów i rejestracji w yników napisano program kom puterowy. Program ten napisano w ję zy k u Pascal ze w zględu na posiadane przez ten języ k wew nętrzne oprogram ow anie bardzo dobrze nadające się do ww. celów. Program ten posiada w iele m ożliw ości operacyjnych.
4. P R Z Y K Ł A D O W E W Y N IK I BA DAN
P rzedm iotem badań było szkło organiczne (plexi). N a rys. 3 przedstaw iono otrzym any z m ikrokom putera przykładow y przebieg zm ian tem peratury powierzchni próbek w funkcji czasu.
Czas X, s-0.5
R y s. 3. P rz y k ła d o w y p rz e b ie g z m ia n te m p e ra tu ry p o w ie rz c h n i p ró b k i w fu n k c ji c zasu F ig. 3. T e m p e ra tu re o f th e sa m p le as a fu n c tio n o f th e tim e
Z astosow anie m etody bilansów elem entarnych. 81
Próbki w ykonane były w kształcie cylindrów o średnicy d = 40 m m i w ysokości 5 = 1 0 mm.
W celu w yznaczenia dyfuzyjności cieplnej w ykonano szereg pom iarów dla różnych wartości generow anej m ocy cieplnej. P rzykładow a w artość mocy cieplnej grzejnika w ynosiła P = 2.12 W, co daw ało rów now ażną gęstość strum ienia ciepła na je d n ą próbkę q = 844 W /m 2. W yko
rzystując dane pom iarow e oraz opracowany algorytm obliczeń num erycznych i program ko- puterow y, obliczono w spółczynnik przew odzenia ciepła i ciepło właściw e. W artości te w yno
szą X = 0.1842 W /(m K ), c = 1490 J/(kgK ). W literaturze w artości te dla plexi podaw ane są w granicach: w spółczynnik przew odzenia ciepła A = (0.174-0.20) W /(m K) i ciepło w łaściw e c = (1420-2090) J/(kgK ). Stąd w idać, że otrzym ane w yniki m ieszczą się w podanych przedziale.
5. W N IO SK I K O Ń C O W E
U ogólniona m etoda najm niejszych kwadratów należy do metod stochastycznie optymalnych.
Zastosow anie jej do rozw iązyw ania odwrotnych współczynnikowych zagadnień przewodzenia ciepła umożliwia:
• jednoznacznie określać najbardziej prawdopodobne wartości niewiadomych,
• uściślenie niedokładności wielkości mierzonych,
• na bieżąco kontrolow ać dokładność obliczeń.
Stąd m oże być z pow odzeniem stosowana do wyznaczania właściwości cieplnych różnych mate
riałów.
Pracę w ykonano w ram ach projektu badaw czego nr 8 T10B 012 18
LITER A TU R A
1. G uzik A., Styrylska T., Zastosow anie filtra c ji opartej na metodzie najm niejszych kw a
dratów do w yznaczania nieustalonych p ó l temperatury, M ateriały X V II Zjazdu Term o
dynam ików , Z akopane-K raków , 1999, str. 507
2. IM SL Library, U se r’s M anual V ersion 1.0, Houston, Texas 1987.
3. Szargut J i inni, R achunek w yrów naw czy w technice cieplnej, PAN W arszawa, 1984
4. Szargut J., S tyrylska T., K olenda Z., Num eryczne modelow anie procesów w ym iany ciepła i m asy uogólniona m etodą najm niejszych kwadratów, M ateriały IX Sym pozjum w ym iany
ciepła i m asy, A ugustów 1995, str. 301
5. Skorek J., Z astosow anie m etod stochastycznych i spektralnych do rozw iązyw ania g ra nicznych zadań odw rotnych przew odzenia ciepła, ZN. Pol. Śl. s. Energetyka z. 119, G li
w ice 1994
6. Tichonov A .N ., A rsenin V ., Solutions o f Ill-Posed Problems. W inston & Sons, W as
hington D.C. 1977
Recenzent: Prof. dr hab. inż. E ugeniusz Kalinowski
A bstract
In this paper a com bined experim ental and num erical m ethod o f determ ining the therm al properties, including the therm al conductivity and specific heat o f the different solid bodies, has been presented. The know ledge o f the true values o f the therm al properties o f m aterials is very im portant for designers and users o f m achines and equipm ent subjected to unsteady therm al influences. It is im portant w hile m athem atical m odelling various therm al processes, as well. The determ ination o f the therm al properties o f a solid body is based on the solution o f the inverse heat conduction problem in an investigated sam ple w ith given boundary condi
tions and a given geom etry. This problem belongs to the class o f ill-posed problem s. It means that the solution o f this problem is usually very sensitive to inaccuracies o f input data, e.g.
m easurem ents errors. To stabilize the solution, a special m athem atical m ethod has to be ap
plied. U p to the present day the m ost w idespread stabilization m ethods o f solving inverse problem s in the therm al techniques are the régularisation m ethod and the least squares ad
justm ent m ethod, often called generalised least squares m ethod. In this paper the generalised least squares m ethod has been applied. This m ethod is based on the search for the minim um o f a suitably defined quadratic function and the adjustm ent o f m easured values w ith errors in connective w ith unknow n values by the so-called constraint equation. In the standard gener
alised least squares m ethod the constraint equations have to be linear. But in the case o f the inverse heat conduction coefficient problem the constraint equations are nonlinear. Therefore,
Zastosow anie m etody bilansów elem entarnych.. 83
in the presented paper the m odified version o f the generalised least squares m ethod was ap
plied. This m ethod enables us to solve nonlinear problem s, w hereas in order to determ ine m easurem ent values in selected points o f the investigated m aterial the com puterized m eas
urem ent stand has been used. The stand consists o f tw o sym m etrical probes heated up by an electrical heater w hich is placed betw een them. The tem perature was m easured by four ther
m ocouples. To control the m easurem ents and to record the m easurem ent values, an own com puter program w as w orked out. Constraint equations w ere form ulated using a suitable m athem atical m odel o f heat conduction in the sample. This model was form ulated using the energy elem ental balance m ethod, the so-called control volum e m ethod. This m ethod belongs to the class o f differential m ethods and has an essential clear physical interpretation. There
fore it is especially effective in m athem atical m odelling o f various therm al problem s. The obtained equations described the tem perature distribution w ithin the sam ples, w hich a heated up on the test stand. B ecause the generalised least squares m ethod requires m ore equations than unknow n values, so to com ply w ith the requirem ents, the unsteady heating process o f the sam ple w as investigated and constraint equations for some tim e steps w ere form ulated. Next, using this m ethod, the inverse heat conduction problem was solved and the therm al conduc
tivity and specific heat o f the sam ple w ere determ ined. A Exem plary research results have been presented.
