ZE SZ Y T Y N A U K O W E PO L IT EC H N IK I ŚLĄ SK IEJ Seria: EN ER G ET Y K A z. 133
2001 N r kol. 1486
Stanisław K U C Y PER A Politechnika Śląska, G liw ice
W Y Z N A C Z A N I E P R Z E W O D N O Ś C I C I E P L N E J I C I E P Ł A W Ł A Ś C I W E G O M A T E R I A Ł Ó W S T A Ł Y C H Z W Y K O R Z Y S T A N I E M S Z E R E G Ó W F O U R I E R A I M E T O D O P T Y M A L I Z A C J I
Streszczenie. W pracy w celu rozwiązania odwrotnego zagadnienia przewodzenia cie
pła i w yznaczenia właściwości cieplnych ciał stałych połączono metodę szeregów Fouriera z metodami optymalizacji. Podano założenia m odelu m atem atycznego oraz algorytm rozw iązania prostego zagadnienia przew odzenia ciepła. Zagadnienie odw rotne sform u
łow ano ja k o problem optym alizacyjny i w ykorzystując odpow iednie dane pom iarow e rozw iązano go m eto d ą bezgradientow ą H ooka-Jeevsa oraz zm odyfikow aną m etodą gra
d ientow ą N ew tona. Podano kró tk ą charakterystykę zastosow anych m etod optym alizacji.
O pisano skom puteryzow ane stanow isko badawcze. Przedstaw iono przykładow e w yniki badań. W yniki porów nano z danym i literaturow ym i.
D ET ER M IN A T IO N OF THE THERM AL CO N D U C TIV ITY A N D SPECIFIC HEAT OF SO LID M A TERIALS B Y APPLYING THE FO URIER SERIES A N D O PTIM IZATIO N M ETHODS
Sum m ary. In order to solve inverse heat conduction problem and to determ ine the therm al properties o f solids the Fourier series and optim isation m ethods have been com bined. The fundam entals o f the m athem atical model and algorithm o f solving the direct heat conduction problem have been given. The inverse problem was form ulated as an optim isation problem and, using proper m easurem ent data, it has been solved by applying the H ook-Jeeves non gradient m ethod and N ew ton m odified gradient method.
Short characteristic o f both m ethods has been given. The com puterized test stand has been described. The sam ple results have been presented. The results received by apply
ing the present m ethod have been com pared w ith date from literature.
1. W ST Ę P
W yznaczanie w łaściw ości cieplnych ciał stałych z zastosow aniem m etod optym alizacyj
nych i danych pom iarow ych polega na rozw iązaniu odw rotnego zagadnienia w spółczynniko
62 S. K ucypera
w ego przew odzenia ciepła i poszukiw aniu w artości jednego lub kilku w spółczynników , które o p isu ją w łaściw ości cieplne ciała. R ozw iązyw anie zagadnienia odw rotnego przew odzenia ciepła przeprow adzane je st na ogół w dw óch etapach:
W pierw szym etapie, na podstaw ie odpow iednio sform ułowanego m odelu m atem atyczne
go, rozw iązyw ane je st proste zagadnienie przew odzenia ciepła. W ówczas w yznaczane są w ielkości pom ocnicze lub ich charakterystyki (np. pole tem peratury ). W ielkości te w ykorzy
styw ane s ą w drugim etapie rozw iązyw anego problem u.
W drugim etapie, korzystając z danych pom iarow ych i w yznaczonych wcześniej wielkości, rozw iązyw ane je s t zagadnienie odw rotne i w yznaczane są wielkości docelowe.
W przypadku rozw iązyw ania zagadnienia prostego przew odzenia ciepła korzystnie je st w m iarę m ożliw ości tak sform ułow ać m odel m atem atyczny i przeprow adzić eksperym ent, aby w w yniku rozw iązania prostego zagadnienia przew odzenia ciepła w yznaczane w ielkości, np.
pole tem peratury, otrzym ać w postaci form uł analitycznych.
U w zględniając pow yższe w ym ogi w pracy tej tak sform ułow ano m odel m atem atyczny i przeprow adzono eksperym ent, że rozw iązując zagadnienie proste, rozkład tem peratury w próbkach otrzym ano w postaci szeregów Fouriera. N atom iast efektem końcow ym rozw iąza
nia odw rotnego zagadnienia w spółczynnikow ego było wyznaczenie w spółczynnika przew o
dzenia ciepła k i ciepła w łaściw ego c badanego m ateriału. Problem ten sform ułow ano jako zagadnienie optym alizacyjne i do jego rozw iązania zaproponow ano dw ie metody optym aliza
cyjne: bezgradientow ą m etodę H ooka-Jeevsa i gradientow ą m etodę zm iennej m etryki Davi- dona-F letchera-P ow ella (zm odyfikow aną m etodę Newtona).
Z aproponow ana m etoda rozw iązania odw rotnego zagadnienia w spółczynnikow ego prze
w odzenia ciepła p olega na połączeniu szeregów Fouriera z m etodam i optym alizacyjnym i.
Szeregi Fouriera otrzym ano ja k o rozw iązanie prostego zagadnienia przew odzenia ciepła, na
tom iast m etody optym alizacyjne w ykorzystano do sterow ania rozw iązaniem zagadnienia pro
stego i je g o m odyfikacji w celu w yznaczenia wartości w spółczynnika przew odzenia ciepła i ciepła w łaściw ego badanego m ateriału. Sprow adza się to do wielokrotnego rozw iązyw ania zagadnienia prostego z m odyfikow anym i wartościam i wielkości poszukiw anych. W w yniku połączenia eksperym entu i zaproponow anej m etody rozw iązania otrzym ano efektyw ny spo
sób w yznaczania param etrów cieplnych ciał stałych. Przedstawiono przykładow e wyniki.
W yznaczanie przew odności cieplnej i ciepła. 63
2. ZAŁ O ŻE N IA M O DELU M A T EM ATYC ZNEG O I ALG O RYTM O BLICZEŃ PR O STEG O ZA G A D N IEN IA PR ZEW O DZENIA CIEPŁA
2.1. Założenia m odelu m atem atycznego
W pracy podczas analizy teoretycznej rozpatrywano proces nagrzew ania nieskończenie roz
ległej płyty o grubości 26. W celu ułatw ienia rozwiązania zagadnienia i zapew nienia możliwie najlepszej zgodności m odelu m atem atycznego z rzeczywistymi warunkam i pom iarów przyjęto następujące założenia:
1. Rozpatruje się proces sym etrycznego nagrzew ania nieskończenie rozległej płyty, tzn.
płyty, której grubość je st znacznie m niejsza od rozm iarów w zdłużnych, czyli pole tem pe
ratury w płycie je s t polem jednow ym iarow ym . 2. M ateriał płyty m a izotropow e w łaściw ości cieplne.
3. Z akłada się, że w łaściw ości cieplne m ateriału nie zm ieniają się podczas badań.
4. 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ów ny tem peraturze początkowej.
5. W chw ili t = 0 na obu pow ierzchniach zew nętrznych płyty zaczyna działać stały strum ień ciepła q = idem , W /m 2.
6. Z akłada się, że całe w ygenerow ane ciepło w nika do w nętrza płyty, czyli w spółczynnik w nikania ciepła a do otoczenia płyty je st równy zero, co w pełni spełnione je s t w w arun
kach pom iarow ych, gdyż grzejnik oddaje ciepło do próbek tylko przez przewodzenie.
D la tak p rzyjętych założeń sform ułow ano m odel m atem atyczny i rozw iązano zagadnienie proste w ykorzystując m etodę szeregów Fouriera.
2.2. M odel m atem atyczny i algorytm obliczeń zagadnienia prostego
U w zględniając pow yższe założenia rów nanie opisujące jednow ym iarow y przepływ ciepła w płycie ( dla w spółrzędnych bezw ym iarow ych) m ożna opisać następującym rów naniem róż- niczkow ym [2]:
64 S. Kucypera
gdzie:
0 = T -Tp
, (2 )(3)
Fo = - k r
F o > 0 .
c-p-S
D la przyjętych założeń warunki brzegowe w yrażają następujące równania:
- dla powierzchni płyty:
(4)
80_ = +
q.S
(5a)
i=±;
- dla osi płyty:
80_ = o (5b)
P oczątkow a tem peratura w próbkach je st w yrów nana, w ięc w arunek początkowy m a postać:
0
L =o.
(5c)Rozwiązanie zagadnienia początkowo-brzegowego (l-5 c ) przyjęto jako sumę dwóch funkcji:
przy czym:
0 = 0 t + 0 2
0
.* I I +0,1
=0 .F o—0 2 | f-o = 0
Funkcja (9, jest rozw iązaniem następującego zagadnienia brzegowego:
80x d20t dFo dĘ2
( 6)
(7)
(8)
8 0,
W
= +
< 1-8
f-±/
(9)
W yznaczanie przew odności cieplnej i ciepła. 65
N atom iast funkcja 0 2 spełnia równania:
5(9, d 29 2
—- = --- f , (10)
dFo d ą 1
34
= o
.
Funkcję 6?, zakłada się w następującej postaci:
ą = A , - F o + Ą + Ą ? .
Rozw iązanie rów nań (10-11) przyjęto w postaci iloczynu funkcji:
9 2 = X ( t ) T( Fo) .
(11)
(12)
(13) Stąd rozwiązanie zagadnienia (1) z warunkami brzegowym i (5) we współrzędnych bezw ym ia
row ych m a następującą postać:
0( ą. Fo) =q - 8
F o Ą - f l - y ą *
n ¡-1 i
(14)
O statecznie rozkład tem peratury w próbce T ( x , z ) opisany je st następującym równaniem:
T(x, t) = Tp + ^~j~ k r 1 2 * r - i /
c p - 8 6 ] ---y l - y — c o s (i n — )e (15)
O trzym any z rozw iązania zagadnienia prostego rozkład tem peratury w próbce wykorzystano po uw zględnieniu wielkości pom iarowych do rozwiązania odwrotnego zagadnienia współczyn
nikow ego przew odzenia ciepła.
3. R O Z W IĄ Z A N IE O D W R O T N E G O ZA G A D N IE N IA PR Z E W O D Z E N IA CIEPŁA
G łów nym celem analizow anego problem u je st w yznaczenie ciepła w łaściw ego i w spół
czynnika przew odzenia ciepła badanego m ateriału. Do rozw iązania zagadnienia zapropono
w ano m etody optym alizacyjne, polegające na poszukiw aniu nieznanych składow ych w ektora
66 S. Kucypera
x poprzez m inim alizację sum kw adratów różnic pom iędzy w artościam i tem peratury zm ierzo
nej i obliczonej w tych sam ych punktach próbki. Stąd w arunek ten m ożna zapisać:
MF
f ( * ) = ' L ( TpM - Tpm ) => m in , (16)
p=\
gdzie:
M P - liczba punktów , w których m ierzona była tem peratura (obie zew nętrzne pow ierzchnie próbki),
Tp.obi, Tp.zm — odpow iednio tem peratury: obliczona i zm ierzona na stanow isku pom ia
row ym w p - tym punkcie próbki,
x - w ektor zaw ierający poszukiw anie wielkości: przew odność cieplną k i ciepło w łaściw e c m ateriału, czyli:
x = M r . (1?)
D odatkow ym celem pracy było znalezienie efektyw nego algorytm u rozw iązania odw rot
nego zagadnienia w spółczynnikow ego przew odzenia ciepła sform ułow anego zależnością (16). D latego przeanalizow ano m ożliw ość zastosow ania dw óch m etod optym alizacyjnych:
m etody bezgradientow ej H ooka-Jeevsa,
gradientow ej zm odyfikow anej m etody N ew tona tzw. m etody zm iennej m etryki Davidona- Fletchera-P ow ella (DFP).
Obie w ym ienione m etody należą do klasy m etod iteracyjnych.
Istota m etody H ooka-Jeevsa [1] polega na tym , że przystępując do obliczeń określa się ba
zę n liniow o niezależnych w ektorów d(l), d (n) (które są najczęściej w ersoram i karte- zjańskiego układu w spółrzędnych) oraz w spółczynnik kroku r>0. K ażda iteracja m etody składa się z dw óch etapów: próbnego i roboczego. Etap próbny służy do zbadania lokalnego zachow ania się funkcji celu F (xk) w niew ielkim otoczeniu punktu x k poprzez w ykonanie kro
ków próbnych o długości t we w szystkich kierunkach ortogonalnej bazy. Etap roboczy polega na przejściu do następnego punktu, w okół którego realizowany będzie następny etap próbny, lecz tylko w tedy, gdy przynajm niej jed en z w ykonyw anych kroków daje zm niejszenie w arto
ści funkcji celu. W przeciw nym razie etap roboczy je st pom ijany, a etap próbny w ykonyw any je st ponow nie ze zm niejszoną w artością kroku x. Podstaw ow ym kryterium zakończenia obli
czeń je s t x<£, gdzie e je st dokładnością obliczeń.
Z kolei podstaw ow ą id e ą m etody zm iennej m etryki [1] je st realizacja kolejnych iteracji i poszukiw ania m inim um funkcji celu w yznaczając kierunek poszukiw ań z następującej zależ
ności:
W yznaczanie przew odności cieplnej i ciepła.. 67
d (k) = - A - ' ( \ (k))■ V F ( x , k) ) , (18)
gdzie A (x(k)) je s t m a cie rz ą hesjanu w punkcie x(k). Stąd w m etodzie N ew tona w ym agana jest znajom ość m acierzy drugich pochodnych A(x) =<52F{x)ldx\dxy Ponadto m usi istnieć m acierz odw rotna, a jej elem enty są ciągłe i dodatnio określonym i elem entam i m acierzy hesjanu A(x).
następnie, w w yniku m inim alizacji w tym kierunku, otrzym uje się x*+l. Z a w arunek zakoń
czenia obliczeń przyjm uje się:
gdzie e je s t dokład n o ścią obliczeń.
4. O PIS ST A N O W ISK A PO M IA R O W E G O
U kład pom iarow y składał się z cienkiego (w postaci folii) oporow ego grzejnika elektrycz
nego, dw óch sym etrycznie rozm ieszczonych w zględem grzejnika w alcow ych (płytek) próbek oraz czterech cienkich płytek m iedzianych z term oparam i. C ałość zaizolow ana była rozdrob
nionym korkiem . C ienka w arstw a oporow a grzejnika zapew niała w yrów naną tem peraturę na jej pow ierzchni. Do zew nętrznych pow ierzchni grzejnika i próbek przylegały cienkie płytki
m iedziane z przym ocow anym i term oparam i.
W szystkie term opary były typu N i-NiCr. Do zasilania grzejnika użyto zasilacza stabilizow a
nego sterow anego z m ikrokom putera. N apięcie U i prąd I przepływ ający przez grzejnik m ie
rzone były poprzez kartę pom iarow ą na m ikrokom puterze. Z e w zględu na bardzo m a łą gru
bość grzejnika przyjęto, że strum ień ciepła dopływ ający do próbek rów na się m ocy grzejnika:
K ażda iteracja m etody polega na tym , ze m ając dany punkt x k generuje się kierunek d A (18), a
(19)
Q = P = U - I . (20)
G ęstość strum ienia ciepła dopływ ającego do każdej z próbek określa zależność:
2 P W
(21)
Schem at stanow iska pom iarow ego przedstaw iono na rys. 1.
68 S. K ucypera
X i cr
Rys. 1. S chem at układu pom iarow ego: ł-g rz e jn ik , 2—płytki z term oparam i, 3 -próbki, 4 -izo lacja cieplna Fig. 1. S chem e o f m easurem ent system : 1- heater, 2- plates w ith therm ocouples, 3- sam ples, 4- therm al insula
tion
Stała term oelektryczna użytych term opar s = 25 K/mV. Term opary podłączone były po
przez w zm acniacz (ze w zględu na n iską w artość sygnałów z term opar) i m ultiplekser ( ko
nieczność przełączania odczytyw anych w każdym kroku czasow ym sygnałów z term opar) oraz kartę p o m iaro w ą do m ikrokom putera. Tem peratura pow ierzchni próbek od strony grzej
nika i pow ierzchni zaizolow anych m ierzona by ła term oparam i połączonym i szeregow o. D la
tego w celu dalszej obróbki danych w zbiorze w ynikow ym zapisyw ano średnią arytm etyczną w artość w skazań odpow iednich term opar. Do sterow ania całością pom iarów i rejestracji w y
nikó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 w ew nętrzne oprogram ow anie (funkcje w ew nętrzne) bardzo dobrze nadające się do ww. celów. Program ten je st dosyć uniw ersalny i um ożliw ia między innymi:
1. Sterow anie z klaw iatury w artościam i mocy grzejnika i odłączanie grzejnika w przypadku korzystania z zew nętrznego źródła ciepła.
2. O bserwowanie na monitorze wartości prądu i napięcia pochodzącego z grzejnika.
3. Zadawanie z klawiatury liczby używ anych term opar (max 24).
4. O bserw ow anie na m onitorze zm ian tem peratury na w szystkich używ anych term oparach lub tylko w ybranych.
W yznaczanie przew odności cieplnej i ciepła. 69
5. Zadaw anie częstości odczytu i rejestracji wartości tem peratury z term opar oraz prądu i napięcia z grzejnika.
6. Zadaw anie z klawiatury m aksym alnych wartości obserwowanych tem peratur (w przypadku gdy układ nie dochodzi do stanu ustalonego.
7. Zadawanie z klawiatury czasu trw ania eksperymentu.
8. Zapisywanie w yników pom iarów do zbioru w celu dalszej ich obróbki.
Dokładniejszy opis stanow iska pom iarowego i sposobu pom iarów podano w [2,3],
5. P R Z Y K Ł A D O W E W Y N IK I BA D A Ń
Przedm iotem badań było szkło organiczne (plexi). N a rys. 2 przedstawiono przykładowe przebiegi zm ian tem peratury pow ierzchni próbek w funkcji czasu.
60 — ---
2 0
---
0 400 80
Czas, s
R y s.2 . 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 e k w fu n k c ji c za su F ig . 2. T e m p e ra tu rę d is trib u tio n as a fu n c tio n o f th e tim e
Próbki w ykonane były w kształcie płytek w alcow ych o średnicy d = 71,9 m m i wysokości 8 = 12.45 mm. W celu w yznaczenia w spółczynnika przewodzenia ciepła k i ciepła właściwego c wykonano szereg pom iarów dla różnych wartości generowanej mocy cieplnej. Przykładowa wartość mocy cieplnej grzejnika w ynosiła P = 4.8 W, co dawało równow ażną gęstość strumienia ciepła n a je d n ą p ró b k ę q = 591.4 W /m2.
70 S. Kucypera
W ykorzystując dane pom iarowe oraz opracowany na bazie przedstawionego algorytmu obli
czeń program kom puterowy, obliczono obiema metodami optymalizacyjnymi średnie wartości w spół-czynnika przewodzenia ciepła i ciepła właściwego. Wartości otrzymane obiema metoda
mi są praw ie takie same i w ynoszą Ar= 0.1839 W /(m K), c = 1437.6 J/(kgK).
6. W N IO SK I I U W AG I K OŃCOW E
O trzym ane w yniki badań porównano z danym i literaturow ym i [4], które dla badanego m ateriału podaw ane są w granicach: w spółczynnik przew odzenia ciepła k = (0.174-0.20) W /(m K ) i ciepło w łaściw e c = (1420-2090) J/(kgK ). Jak widać, przedstaw ione wyniki m ieszczą się w podanych granicach. Stąd podany algorytm m oże być zastosowany do w yznaczania param etrów term ofizycznych różnych m ateriałów stałych. Podstaw ow ą zaletą zaproponow anej m etody, podobnie ja k innych m etod odw rotnych przew odzenia ciepła je st to, że m oże być stosow ana do w yznaczania kilku param etrów term ofizycznych: np. a, c, k w czasie w ykonyw ania jednej serii pom iarów. D odatkow ą zaletą jej je st rów nież dość duża szybkość otrzym ania w yników . Stąd w dalszych badaniach zam ierza się zm odyfikować algorytm obliczeniow y w ten sposób, aby m ożna było w yznaczyć bezpośrednio ww.
w ielkości w funkcji tem peratury, podając do algorytm u obliczeń analityczne zależności na k i c w funkcji tem peratury. Jeżeli natom iast chodzi o porów nanie metod optym alizacyjnych, to m etoda gradientow a je s t szybsza pod w arunkiem przyjęcia startow ych wartości niew iadom ych bliskich rzeczyw istym , w przeciw nym przypadku algorytm staje się rozbieżny.
N ie zauw ażono tego dla m etody H ooka-Jeevsa. W przyszłości, w celu zwiększenia efektyw ności algorytm u obliczeniow ego, zam ierza się połączyć m etodę bezgradientow ą z gradientow ą.
Pracę w ykonano w ram ach projektu badaw czego nr 8 T l OB 012 18.
L ITER A TU R A
1. K ręglewski T i inni, M etody optymalizacji w ję zyk u FORTRAN, Warszawa 1984.
2. K ucypera S., N adziakiew icz J., W yznaczanie współczynnika przewodzenia ciepła i ciepła właściwego materiałów na podstaw ie charakterystyk tem peraturow ych próbek, M ateriały.
X X X V I Sym pozjonu M odelow anie w M echanice, t. 2, ZN KM T, G liw ice 1997.
W yznaczanie przew odności cieplnej i ciepła. 71
3. K ucypera S., Skorek., W yznaczanie właściw ości cieplnych ciał stałych z wykorzystaniem m etody filtra c ji dynam icznej, M ateriały X V II Z jazdu T erm odynam ików , Zakopane-
K raków , 1999, str. 689
4. R aźniević K ., Tablice cieplne z wykresam i, W NT, W arszaw a 1966
Recenzent: Dr inż. Antoni Guzik
Abstract
The determ ination o f the therm al properties o f solid bodies by m eans o f optim isation m ethod and m easurem ents data is based on the solution o f the inverse heat conduction coeffi
cient problem and on searching for the value o f one or m ore coefficients describing the ther
m al properties o f the body. The solution o f the inverse heat conduction problem is generally perform ed in tw o stages:
In the first stage, basing on a suitably form ulated m athem atical m odel, the direct heat con
duction problem is solved. The auxiliary m easures and their characteristics (for exam ple the tem perature field) are also determ ined. These w ill be used in the second stage o f the problem o f solving the algorithm .
In the second stage m aking use o f m easurem ent data and previously determ ined quantities, the inverse problem is solved and the final quantities are determ ined.
In the present paper th e m athem atical m odel has been form ulated and the experim ent car
ried out in the w ay that the tem perature distribution w as received in the Fourier series. The final effect o f solving the inverse heat conduction coefficient problem are the values o f the heat conduction coefficient k and specific heat c o f the investigated m aterial. This problem w as form ulated as an optim isation problem and solved by m eans o f optim isation m ethod e.g.
H ook-Jeeves m ethod and the m odified N ew ton m ethod.
In the presented paper an original m ethod o f the solving inverse heat conduction coeffi
cient problem w as proposed. This m ethod is based on com bining Fourier series w ith optim i
sation m ethod. F ourier series are obtained by solving the direct heat conduction problem , w hereas optim isation m ethods w ere used to control the calculations and to m odify the prob
lem in order to determ ine the values o f the heat conduction coefficient and specific heat o f the
72 S. Kucypera
material. The proposed m ethod is based on m ultiple solutions o f the direct problem . As a re
sult o f com bining the experim ent and the proposed way o f the solving problem an effective m ethod o f determ ination the therm al param eters o f the solid bodies has been found. The re
search results have been presented. The results received by m eans o f present m ethod has been com pared w ith the results by m eans o f previously used m ethods and data from literature.
