Efektywna metoda identyfikacji nieustalonego stanu cieplnego ciśnieniowych elementów kotłów w czasie rzeczywistym

ZESZYTY NAUKOWE POLITECHNIKI ŚLĄSKIEJ Seria: ENERGETYKA z. 122

1994 Nr kol. 1262

Jan TALER, W iesław ZIMA

Instytut A p aratu ry Przemysłowej i E nergetyki, P olitechnika K rakow ska

EFEKTYWNA METODA IDENTYFIKACJI NIEUSTALONEGO STANU CIEPLNEGO CIŚNIENIOWYCH ELEMENTÓW KOILÓW W CZASIE RZECZYWISTYM

S tr e sz c z e n ie . W referacie przedstaw iono num eryczno-analityczną metodę rozw iązyw ania zagadnień odwrotnych przew odzenia ciepła, wy­

stępujących przy w yznaczaniu n ap rężeń cieplnych w ciśnieniowych, grubościennych elem entach kotłów i tu rb in . M etoda n ad aje się zarówno do zagadnień je d n o - ja k i wielowymiarowych. J e s t ona szczególnie przy d atn a do zastosow ania w kom puterow ych układach kontroli n a p rę ­ żeń cieplnych, pracujących w czasie rzeczywistym. Odznacza się ona dużą dokładnością, n aw et przy dużych, przypadkow ych błędach pom ia­

ru tem p eratu ry . Wzory obliczeniowe m ają przy tym p ro stą formę. Za­

stosow anie m etody zilustrow ano trzem a przykładam i.

EFFECTIVE METHOD OF IDENTYFICATION OF NONSTEADY STATE OF PRESSURIZED BOILER ELEM ENTS IN REAL TIME

Sum m ary. An accurate and stable sem inum erical m ethod is developed for solving inverse h e a t conduction problem s encountered in m onitoring of th erm al stresses in pressurized th ick -w alled elem ents of steam boilers and turbin es. The m ethod is ap p ro p riate for one and m ultidim ensional problems. Due to its efficiency, it can be readily im plem ented in an o n -lin e fatigue m onitoring system . The num erical resu lts show the c u rre n t m ethod of analysis is insensitive to m easu rem en t errors. The m ethod is very sim ple and is m uch more com putationally efficient th a n other m ethods. The application of the m ethod is illu strate d by th re e examples.

3 0 0 E K T H B H W K M ETO fl HßEHTHOMKAJJHH HECTAIfHO HAPHO rO T E n jIO B O rO COCTOAHHA K O H C T P yK IfH O H H blX 3AEM EH T0B KOTJIOB PABO TAIO m H X nO A AABAEHHEM B PEAAbHOM BPEMEHH

Pe3K)Me. B cTaTbe upcącTaBjicu hhcjichho— aHajnrmHecKHii m c t o a oiipeąejieHHsi o ö p a T H H x npoßjieMOB TciuiooßvfeHa, K O T o p H e BHCTynaioT ripa onpepejreHHM T e M n e p a T y p H H X H a n p s a c e H H H b tojictoctchhhx 3jre\ieHTax KOTJIOB H T y p Ö H H paßOTaiOIUHX IlO/i; p a B a e H H C M . MctO/I npHrOpHhlil fljia OflHO- h M H o r o M e p H H X sapai, o c o 6 c h h o k n p u M e H e H H i i b K O M n b i o T e p H H X cHc re Ma x K O H T p o j w T e M n e p a T y p H u x Haupsi/Kennii padoTaroignx b peajibHOM BpeMeHH.

OTJIHHaeTCH O H ÓOJIblHOH TOHHOCTbK), flä*e n p H ÖOJIblHHX, CJiynaHHHX o m n Ö K a x H 3 M e p e H H a T e M n e p a T y p H . BbiBCAeiiHhie cfiopMyjibi BbnracjieHHa paciipe/rejicuna T e M n e p a T y p u h HanpsLaceHHÜ hccjiojkhmc. n p H M e H e H i i e M C T o p a n o K a 3 a H o H a T p e x np HM ep a x.

W STĘP

W grubościennych elem entach bloków energetycznych pracujących przy wysokich tem p eratu rach i ciśnieniach mogą powstawać w czasie rozruchu i w yłączania z ru ch u wysokie naprężen ia cieplne, wywołane nierównom iernym rozkładem tem peratury.

W celu określenia naprężeń n a wew nętrznej powierzchni omywanej czynni­

kiem pod wysokim ciśnieniem niezbędna je s t znajomość te m p e ra tu ry tej po­

wierzchni, n a której zamocowanie czujnika je s t utrudnione.

W referacie proponuje się nową m etodę w yznaczania czasow o-przestrzen- nego rozkładu tem p e ra tu ry oraz naprężeń cieplnych w ew nątrz elem entu i na jego powierzchni wew nętrznej n a podstaw ie tem p e ra tu ry mierzonej na

powierzchni zewnętrznej tego elem entu.

Ponadto m ożna wyznaczyć gęstość stru m ien ia oraz w spółczynnik w nikania ciepła n a w ew nętrznej powierzchni elem entu.

Tak więc, bez w iercenia otworów w ściance odtworzyć m ożna zjawiska cieplne zachodzące n a w ew nętrznej powierzchni, co je s t szczególnie ważne w przypadku elem entów konstrukcyjnych pracujących pod wysokim ciśnieniem.

Sform ułowany w ten sposób problem je st odwrotnym, źle uw arunkow anym zagadnieniem przew odzenia ciepła, bardzo tru d n ym do poprawnego rozwiąza­

nia. W ynika to z ogromnej wrażliwości wyznaczanej tem p e ra tu ry lub gęstości n a powierzchni w ew nętrznej n a przypadkowe błędy pom iaru zaw arte w prze­

biegach czasowych tem p e ra tu r m ierzonych n a powierzchni zew nętrznej. N a­

w et niew ielkie błędy pom iaru tem p e ra tu ry mogą wywołać bardzo duże oscyla­

cje w przebiegach wyznaczanych param etrów [1]. Do w yznaczania pola tem ­ p e ra tu ry i naprężeń cieplnych w elem entach cylindrycznych n a podstawie pom iaru tem p e ra tu ry zewnętrznej powierzchni stosow ane są wzory w ażne dla ścianki płaskiej [2, 3],

Efektywna metoda identyfikacji nieustalonego stanu. 139

W przypadku cylindrycznych elem entów grubościennych tak ie uproszcze­

nie może prowadzić do dużych błędów.

Przedstaw iona w referacie m etoda n u m eryczno-analityczna rozwiązywa­

nia odwrotnych zagadnień przew odzenia ciepła odznacza się dużą dokładnoś­

cią i stabilnością, n aw et przy dużych błędach tem p e ra tu r. Może być wykorzy­

stana przy wyznaczaniu je d n o - ja k i wielowymiarowych pól tem p eratu ry . Opracowano ponadto oryginalny sposób w ygładzania zm ierzonych przebie­

gów tem p e ra tu ry za pomocą filtru cyfrowego, elim inujący w znacznym stop­

niu wpływ błędów pom iaru tem p e ra tu ry n a uzyskiw ane wyniki.

Zastosow anie proponowanej w referacie m etody zilustrow ano n a trzech przykładach.

W pierwszym, na podstawie zmierzonego przebiegu tem peratury na izolowa­

nej, nie ogrzewanej cieplnie powierzchni ściany płaskiej, odtworzono przebieg zmian gęstości strum ienia ciepła na przeciwległej powierzchni ogrzewanej.

W pozostałych dwóch przykładach n a podstaw ie pom iaru tem p e ra tu ry na powierzchni zew nętrznej zbiornika cylindrycznego wyznaczono rozkład tem ­ p eratury n a grubości ścianki oraz w spółczynnik w n ik an ia ciepła n a w ew nętrz­

nej powierzchni zbiornika w alczaka kotła w funkcji czasu. Obliczono również przebieg osiowych n aprężeń cieplnych n a w ew nętrznej powierzchni elem entu w funkcji czasu. W przykładzie drugim te m p e ra tu ra czynnika zm ienia się skokowo; w trzecim n a to m ia st rośnie najpierw ze s ta łą szybkością, a n a s tę ­ pnie pozostaje stała.

OPIS METODY

W zagadnieniach odwrotnych przew odzenia ciepła, n a podstaw ie zmierzo­

nego przebiegu tem p e ra tu ry w w ybranych p u n k tach leżących w ew nątrz ciała, odtw arzana je s t te m p e ra tu ra brzegu ciała albo gęstość stru m ie n ia lub współ­

czynnik w nikania ciepła. Z uw agi n a fakt, że je s t to problem źle postawiony, duży wpływ n a poszukiw ane rozw iązanie m ają przypadkow e błędy pom iaru tem p eratu ry [1].

W pracy [4] przedstaw iono efektyw ną m etodę rozw iązyw ania odwrotnych zagadnień przew odzenia ciepła, wychodząc z dobrze znanej m etody bilansów elem entarnych [5].

W niniejszej pracy w podobny spasób wyprowadzone zo staną rów nania pozwalające wyznaczyć tem p e ra tu rę i gęstość stru m ie n ia ciepła n a powierzch­

ni ciała cylindrycznego n a podstaw ie znanego przebiegu te m p e ra tu ry i gęstoś­

ci stru m ien ia ciała w punkcie położonym w ew nątrz ciała (rys. 1).

Zakładając, że w punkcie r z znane są z pom iaru przebieg tem p e ra tu ry f3(t) i gęstość stru m ien ia ciepła qE(t), wyznaczona zostanie te m p e ra tu ra wewnę­

trznej powierzchni cylindra T ^ t) oraz gęstość stru m ie n ia n a wew nętrznej powierzchni ą^ft).

Rys. 1. S c h e m a t p o d ziału śc ia n k i cylindrycznej n a objętości k o n tro ln e ; f3 i qE - wielkości zn a n e z p o m ia ru , T i i qs - w ielkości p o szu k iw an e

Fig. 1. O n e-d im en sio n a l m odel for cylindrical geom etry show ing elem en ts a n d control volu­

m e boundaries; f3 a n d qE - q u a n titie s from m easu rem en ts, T i a n d qs - un k n o w n q u an tities

Wymienione wielkości wyznaczone zostaną z rów nań bilansu ciepła dla objętości kontrolnej przedstaw ionej n a rys. 1.

dfą dt

N - 1

a3 Mfą) + MT4)

m i - - r 3

r 3 + r 4 (T4 - f3) - Mf3) + ^(T2)

m )

l - r 4

d lą dt

N - l

a 2

r 3 + r 4

< J

k(T2) + k(f3)

h - t 2)-

A.(T2) + k(TQ k(T2)

m t2)

r3 r 2 + r 3

1 - ^{r 2}

r 2 + r 3 (f3 - T2) - (T2 - TO-

ll)

Efektywna metoda identyfikacji nieustalonego stanu. 141

dTj dt

N - 1

r z _ r w

v y

a i MT i) + X( T2)

X(Ti) 1 -

r l + r 2 (T2 - T ! ) - A,(Ti) + A,(T0)

M ^ i )

1 - r l + r 2 (Ti - T0y

Na zew nętrznej i w ew nętrznej powierzchni ścianki cylindrycznej zadane są warunki brzegowe II rodzaju:

Z(T)dT dr

d r

= -<łE

= - q s

(2)

(3)

w których pochodne przybliżone zostaną ilorazam i różnicowymi centralnym i:

* t f s )

T4 - T 2

2Ar ^{- “ ÓE.} (4)

T2 - T q

Z rów nań (4—5) otrzym uje się odpowiednio:

T - t 2ArqE 4 2 m ’

2Arqs

t0 = t2

Z(Ti) ’

(5)

(

)

(7)

r * - r „ gdzie:

A r N - 1 ’

T0 i T4 - te m p e ra tu ry w węzłach pozornych.

W dalszych rozw ażaniach przyjm uje się, że własności term ofizyczne m ate­

riału ścianki: c, X oraz p nie zależą od tem p eratu ry .

Z pierwszego z rów nań (1), po uw zględnieniu rów nania (6), otrzym uje się:

Z drugiego z rów nań (1), po uw zględnieniu rów nania (8), otrzym uje się wzór n a tem p eratu rę Ti powierzchni wewnętrznej:

/ \

r 2 + r 3 (Ar)2 df3 (Ar)4 d2f3

i 2r4 ]ArqE

i \

r 4 (Ar)3 dqE

r 2 a d t + 4a2 d t2 + r 3 + r 4 l k ■'r 3 + r 4 aż, dt Ti = f3 +

Z ostatniego z rów nań (1), uwzględniając (7 - 9), otrzym uje się:

qs = qE

(9)

i \

r 4 (r3 - Ti) (r2 + r 3)

+ ( r 2 + r 3 )

[df3 Ark^{f l} + n

r 3 + r 4

V y r l r 2 dt 2a r l r 2

J

d2f3 (Ar)3k dt2 4a2

ri + r 2 1 1 r ir2 2r, 2r2

d3f3 (Ar)5A, rj + r 2

d t3 16a3 r 3r 2 ⁽10) r 4(r 2 + r 3) dqE (Ar)2i r l + r 2 1 1^ , d2qE (Ar)4( r i + r2)

(r3 + r 4) dt 2a r l r 2

V r l r 2

y dt2 4 a2 r l r 2

Przyjm ując r 3 = r 2 = r 3 = r 4, otrzym uje się z wyżej wyprowadzonych wzorów rów nania n a Tj i qs dla ścianki płaskiej, przedstaw ione w p racy [4]:

Ti - f3 + 1E 2 df3 l E4 d2f3 qEE l E 3 dqE

+ - + -

qs = qE

2 a dt 32 a2 d t2 1 E 5 d3f3

+ TT8 l a d t ’ ⁽¹¹⁾

Ę d f j

a dt T 16 a2 d t2 " 128

l3 . 3 E 3 d2f3

+ - a3 d t3 + l E 2 ^ 1 E4 d2qE 2 a dt + 32 a2 d t2 gdzie E=2Ar.

We wzorach określających tem p e ra tu rę w ew nętrznej powierzchni elem en­

tu T 3 oraz gęstości strum ienia ciepła n a tej powierzchni qs w ystępują pochod­

ne zmierzonej tem p eratu ry f3 i gęstości stru m ien ia ciepła qE po czasie. Z uw a­

gi n a fakt, że zarówno f3(t), ja k i qE(t) obarczone są przypadkow ym i błędami pom iaru, szczególną uwagę należy zwrócić n a popraw ne wyznaczanie tych pochodnych. Szeroko stosowane w analizie num erycznej przybliżanie pochod­

nych ilorazam i różnicowymi je s t zawodne, gdyż prowadzi do bardzo dużych błędów, szczególnie przy pochodnych wyższych rzędów.

Pochodne te m ożna by również obliczyć aproksym ując f3(t) i qE(t) w całym analizow anym zakresie wielom ianem n -teg o stopnia. Taki sposób postępowa­

n ia napotyka jed n ak w praktyce n a ogromne trudności z uw agi n a złożony przebieg f3(t) i qE(t), tru d n y do aproksym acji jed n ą funkcją w całym analizo­

w anym zakresie.

Efektywna metoda identyfikacji nieustalonego stanu. 143

Z tego też względu przedstaw iona zostanie nowa m etoda aproksym acji przebiegów f3(t) i qE(t), pozwalająca z dużą dokładnością odtworzyć niezabu- rzony błędam i przebieg sygnałów oraz ich pochodnych.

FILTRACJA CYFROWA ZMIERZONYCH PRZEBIEGÓW TEMPERATURY Schem at aproksym acji przebiegów czasowych wielkości zmierzonych je st taki sam dla te m p e ra tu r ja k i gęstości stru m ien ia ciepła. Dlatego też omówio­

na zostanie jedynie aproksym acja zmierzonego przebiegu tem p eratu ry . W niniejszej pracy do w ygładzania danych pom iarowych zastosowano fil­

trację cyfrową z uw agi n a wysoką jej skuteczność w elim inow aniu jej przypad­

kowych błędów pom iaru i prostą postać wzorów obliczeniowych. Pozwala ona z dużą dokładnością wyznaczyć rzeczywiste w artości funkcji i jej pochodnych na podstaw ie zakłóconych błędam i przypadkowymi danych pomiarowych. Fil­

tracja cyfrowa je s t odm ianą kroczącej aproksym acji wielomianowej, w której wprowadzono skalow anie czasu i przyjęto założenie, że p u n k ty pom iarowe są równolegle. W ygładzanie przeprow adza się m etodą najm niejszych kw adratów za pomocą ortogonalnych wielomianów G ram a [6], k tóre ap roksum ują grupę danych pom iarowych złożoną z N punktów.

Dokładność średniokw adratow ej aproksym acji wielomianowej je s t najw yż­

sza w środku analizow anego przedziału. Z tego względu do aproksym acji bierze się n iep arzy stą liczbę danych pomiarowych.

W ygładzanie za pomocą filtru cyfrowego przeprow adza się w analogiczny sposób ja k w przypadku aproksym acji wielomianowej kroczącej. Je śli filtr wykorzystuje pięć punktów pomiarowych, to przy pierw szych pięciu p u n k ­ tach: fj, f2, f3, f4 i fs obliczana je s t w artość w ygładzona y(t3), następ n ie bierze się drug ą grupę punktów: f2, f3, f4, f5 i f6 i oblicza się y(t4). W te n sposób otrzymuje się wygładzone wartości dla wszystkich punktów od i = 3 do i = N - 2.

N astępnie wygładza się dwie pierwsze i dwie ostatnie w artości funkcji, korzy­

stając jed n a k z innych postaci wzorów wygładzających. F iltry cyfrowe są dolnoprzepustowe, tj. zatrzym yw ane są sygnały o wysokiej częstotliwości, podczas gdy sygnały o niskiej częstotliwości są przepuszczane.

Załóżmy, że m am y N równoległych punktów pom iarowych (rys. 2), przy czym N > 3 je s t liczbą nieparzystą.

Podstaw iając L = (N—1)/2 w prowadzając now ą zm ienną s:

gdzie At je s t odległością sąsiednich punktów pom iarowych (czasem próbkowa­

nia). W spółrzędną s = 0 otrzym uje się w środku przedziału oraz L punktów

(L=4)

Fig. 2. S m o o th in g of e x p e rim e n ta l d a ta u s in g d ig ita l filte r (L=4)

pomiarowych zarówno po praw ej, ja k i po lewej stronie zerowego punktu środkowego.

Jeżeli wartości tem peratury fs = f(s) są dane w punktach s = -L , -L + 1,..., L, to średniokw adratow e przybliżenie tych w artości za pomocą ortogonalnego wie­

lom ianu G ram a trzeciego stopnia prowadzi przy L = 5 do następujących wzorów wygładzających n a wartości funkcji i jej pochodnych:

y(0) = — ( - 3 6 f 5 + 9f_4 + 44f_3 + 69f_2 + 841! + 89f0 + 841 + 69f2 + 44f3 + 9f4 - 3 6 f 5) ,

(14)

y(0) = 5 1 4 8 A t(3 0 0 f - 5 - 2 9 4 f -4 - 5 3 2 f -3 - 5 0 3 f -2 - 296f_i + 2961 + 503f2 + + 532f3 + 294f4 - 300f5) ,

(15)

Efektywna metoda identyfikacji nieustalonego stanu. 145

y(0) =

143(At)

f X , 2 3 , 2 , 3 2 , J_

ł-5 + 5 ' 15 “ ' 5 ' 3 " 5*1 “ 5 ~ 15*3 + (16) 2„ ,

+ gf4 + fs

y”’(0) =

143(At)

, 1 , 11- 23f 7 „ 23- 11-

- U + 5 L 4 + 15L 3 + 3 0 L 2 + 15 i - 1 5 ti - 3 0 i 2 - 1 5 i3 -

| f 4 + f5

(17)

WYZNACZANIE NAPRĘŻEŃ CIEPLNYCH

Znając rozkład tem p e ra tu ry n a grubości ścianki elem entu cylindrycznego, osiowe naprężenia cieplne n a w ew nętrznej powierzchni m ogą być wyznaczone z następującego wzoru:

aw = 3 ^ ( T ś r - T | r = rw), (18) gdzie:

rz

T śr = o 2 2 I r T (r ’ t ) d r

r z - K ^

Z uw agi n a zastosow aną m etodę num eryczną te m p e ra tu ra ścianki znana jest tylko w dyskretnych punktach, w związku z czym te m p e ra tu ra średnia Tśr jest obliczana w sposób przybliżony. W ykorzystując do obliczania całki we wzorze (19) m etodę trapezów otrzym uje się:

N

T śr = 2 ^ 2 X ( L - l T i - 1 + r iT i)A r ( 2 ° ) r z r w ; _ 2

Uwzględniając, że w zaproponowanym rozw iązaniu odwrotnym tem p e ra tu ­ ra ścianki wyznaczona je s t tylko w trzech punktach: rw, r śr i r z, dokładność obliczania tem p e ra tu ry średniej (20) może być w niektórych przypadkach niew ystarczająca. W związku z tym przeanalizow ano również inny sposób wyznaczania te m p e ra tu ry średniej interpolując wyznaczone z rozw iązania zagadnienia odwrotnego tem peratury: T 1; T2 i T3 parabolą:

gdzie:

T(r, t) = a(t) + b(t)r + c(t)r2, r w(Ti - T2)(2rz - r w) a = Tx

(ń r rw) 2 rz(rsr r w) L _ 2rz(TŁ - T2)

( 4 - 4 ) - 2rz(rsr - r w) ’ T2 - Tj

c = --- (Tsr ~ r w) — 2rz(rsr — r w) Wzór (19) n a tem p e ra tu rę średn ią przyjm uje postać:

2 (r z - 4 ) L rw(Ti - T 2)(2rz - r w) 'i ( 4 - 4 ) 2 1 ( 4 - 4 ) - 2rz(rsr - r )

V 7

(Ti - T2)

(r"sr r w) 2rz(rsr r w) ^{2 r z( r 3z - r l )} ( 4 - 4 )

3 4

Zarówno wzór (20) ja k i (25) m ają c h a ra k te r przybliżony.

(21)

(2 2)

(23)

(24)

(25)

PRZYKŁADY ZASTOSOWANIA PROPONOWANEJ METODY

Przedstaw ione zostaną trzy przykłady ilustrujące dokładność metody.

W pierwszym przykładzie n a podstawie pom iaru te m p e ra tu ry tylnej, izolowa­

nej powierzchni płyty odtworzona zostanie te m p e ra tu ra i gęstość strum ienia ciepła n a przeciwległej, nagrzew anej powierzchni. W drugim i trzecim przy­

kładzie wyznaczona będzie tem p e ra tu ra i współczynnik w nikania ciepła na w ew nętrznej powierzchni w alczaka kotła n a podstaw ie zmierzonego, czasowe­

go przebiegu tem p e ra tu ry n a zew nętrznej, izolowanej cieplnie powierzchni zbiornika. Wyznaczono również przebieg n aprężeń cieplnych n a wewnętrznej powierzchni walczaka.

P rz y k ła d 1

W przykładzie zastosowany będzie ogólnie przyjęty [1] sposób testowania efektywności m etod rozw iązywania odwrotnych zagadnień przew odzenia cie­

pła. Dane pomiarowe generuje się z rozw iązania bezpośredniego zagadnienia przew odzenia ciepła przy założonych, znanych w arun kach brzegowych i wa­

ru n k u początkowym. W danym przypadku przyjm uje się, że początkowa tern-

Efektywna metoda identyfikacji nieustalonego stanu. 147

peratura płyty je s t rów nom ierna. Czołowa pow ierzchnia płyty nagrzew ana jest strum ieniem ciepła o gęstości zmieniającej się w czasie w formie trójk ąta, podczas gdy ty ln a powierzchnia płyty je s t izolowana cieplnie.

Jako dane pomiarowe przyjm uje się obliczoną te m p e ra tu rę powierzchni w odstępach AF0 = 0,06 (AF0 = a ■ At/L2), do której dodano przypadkow e błędy pomiaru o rozkładzie norm alnym . W artość śred n ia generow anych błędów była rów na zeru, a odchylenie średniokw adratow e o | = 0,002, gdzie symbol

° f = Or~:— je s t bezwymiarowym, średnim odchyleniem kw adratow ym . Dane pomiarowe przedstaw iono w tabl. 1. W yniki obliczeń przy danych pomiarowych dokładnych i zaburzonych błędam i przedstaw iono odpowiednio na rys. 3a i 3b.

Fo

Rys. 3a. P rz e b ie g z m ia n g ęsto ści s tru m ie n ia ciepła n a czołowej p o w ierzch n i płyty; d a n e pom iarow e d o k ład n e, 1 - p rzeb ieg p rz y ję ty do obliczeń (w ielkość d o k ła d n a ), 2 - p rzeb ieg obliczony w g m eto d y ścisłej S te f a n a -B u rg g ra f a -L a n g fo rd a , 3 - p rz e b ie g obliczony w g propo­

now anej m etody

Fig. 3a. C a lc u la te d su rfa c e h e a t flu x for t r ia n g u l a r h e a t flux; e rr o rle s s d a ta , 1 - a p p lie d su rfa c e h e a t flu x , 2 - e s tim a te d h e a t flu x u s in g S te f a n - B u r g g r a f - L a n g f o r d m e th o d ,

3 - p re s e n t m e th o d

Fo

Rys. 3b. P rz e b ie g z m ian gęstości s tru m ie n ia ciepła n a czołowej p o w ierzch n i płyty; dane pom iarow e zab u rzo n e b łęd am i, 1 - p rzeb ieg p rz y ję ty do obliczeń (w ielkość dokładna), 2 - przebieg obliczony w g m etody ścisłej S te fa n a -B u rg g ra fa -L a n g fo rd a , 3 - przebieg obliczony

w g proponow anej m etody

Fig. 3b. C a lc u la te d su rfa c e h e a t flux for tr ia n g u la r h e a t flux; d a ta w ith ra n d o m errors, 1 - a p p lie d s u rfa c e h e a t flux, 2 - e s tim a te d h e a t flu x u s in g s te f a n - B u r g g r a f - L a n g f o r d

m ethod, 3 - p re s e n t m eth o d

Obliczenia przeprowadzono przy następujących danych: L = 0,1 m, X = 52 W/mK, a = 14,4- 10~6 m 2/s, At = 41,7 s (AFo = 0,06), = 105 W /m2, T0 = 20°C.

Z analizy rys. 3a i 3b wynika, że dokładność proponowanej m etody przybli­

żonej je s t bardzo dobra.

Efektywna metoda identyfikacji nieustalonego stanu. 149

T a b lic a 1

„D an e p o m ia r o w e ” d o k ła d n e i z a b u r z o n e b łę d a m i d la p ły t y n a g r z e w a n e j s t r u ­ m ie n ie m c ie p ła , z m ie n ia j ą c e g o s i ę w c z a s i e w f o r m ie t r ó j k ą ta

D an e p om iarow e d o k ład n e D a n e p o m iaro w e z a b u rz o n e b łę d a m i

Fo e = T - T °

qNL / K Fo e = T - T °

qNL /K Fo e = T - T °

qNL / K Fo 0 = T - T °

qNL / K

-0 ,2 4 0,0 0,66 0,127200 -0 ,2 4 0,00034 0,66 0,127722

-0 ,1 8 0,0 0,72 0,157880 -0 ,1 8 0,00281 0,72 0,155529

-0 ,1 2 0,0 0,78 0,189293 -0 ,1 2 0,00135 0,78 0,193275

-0 ,0 6 0,0 0,84 0,219593 - 0 ,0 6 -0 ,0 0 0 9 0 0,84 0,223066

0,0 0,0 0,90 0,247680 0,0 -0 ,0 0 0 2 0 0,90 0,248541

0,06 0,000007 0,96 0,272931 0,06 -0 ,0 0 4 2 4 1 0,96 0,275499

0,12 0,000374 1,02 0,295006 0,12 -0 ,0 0 0 6 3 9 1,02 0,293646

0,18 0,002171 1,08 0,313714 0,18 0,001307 1,08 0,311805

0,24 0,006323 1,14 0,328954 0,24 0,007161 1,14 0,330811

0,30 0,013381 1,20 0,340666 0,30 0,012874 1,20 0,342215

0,36 0,023656 1,26 0,348823 0,36 0,021764 1,26 0,348997

0,42 0,037319 1,32 0,353762 0,42 0,038954 1,32 0,350227

0,48 0,054465 1,38 0,356545 0,48 0,054821 1,38 0,355722

0,54 0,075145 1,44 0,358089 0,54 0,072952 1,44 0,359318

0,60 0,099389 1,50 0,358941 0,60 0,098381 1,50 0,361609

P r z y k ła d 2

Obliczono rozkład te m p e ra tu ry w w alczaku oraz współczynnik w nikania ciepła n a jego w ew nętrznej powierzchni n a podstaw ie pom iaru tem p e ra tu ry na jego zew nętrznej, izolowanej cieplnie (qE = 0) powierzchni.

Rozkład te m p e ra tu ry w w alczaku o średnicy w ew nętrznej 2rw = 1,3 m, średnicy zew nętrznej 2rz = 1,48 m obliczono m etodą prostych przy n a stę p u ją ­ cych danych: X = 49,5 W/mK, a = 1,3 10“5 m 2/s, T 0 = 20°C, Tcz = 100°C, a = 1000 W/m2K.

Przyjęto, że zbiornik znajdujący się w rów nom iernej tem p e ra tu rz e począt­

kowej zalano nagle cieczą o tem p eratu rze Tcz = 100°C. Jak o „dane pom iarowe”

przyjęto w artości te m p e ra tu ry n a zew nętrznej powierzchni, obliczone w odstę­

pach czasowych At = 24 s, (AF0 = a ■ At/(rz - r w)2 = 0,0385, do których dodawano liczby pseudolosowe o rozkładzie norm alnym o średniej w artości równej zeru, średnim odchyleniu kw adratow ym równym oy = 0,1394 K. W te n sposób generowano „dane pomiarowe" zaburzone błędam i (rys. 4)

W yniki obliczeń uzyskane z w ykorzystaniem proponowanej m etody przed­

stawiono n a rys. 5 - 7. Z rys. 5 i 6 w ynika, że zarówno rozkład tem p eratu ry

oi— i

ł-

G3 u-

<—D ■

ai-

<o

CL

E<u 110

100 90 80 70 60 50 40 30 20 10

0

/ /

/

/ /

2 0 0 400 600 800 1000 1200 1400 1600 czas t[s ]

R ys. 4. „D ane p o m iaro w e” zab u rzo n e b łę d a m i p rz y a n a liz ie szo k u cieplnego Fig. 4. T e m p e ra tu rę d a ta w ith ra n d o m e rro rs for th e r m a l sh o ck case

Rys. 5. R ozkład te m p e r a tu r n a gru b o ści śc ia n k i z b io rn ik a d la w y b ra n y c h p u n k tó w czaso­

w ych; 1 - 72 s, 2 - 33 s, 3 - 600 s, 4 - 864 s, 5 - 1128 s, 6 - 1392 s; lin ia ciąg ła - rozw iąza­

n ie n u m ery czn e, p u n k ty — ro zw iązan ie o d w ro tn e d la te m p e r a tu r zab u rzo n y ch Fig. 5. C o m p ariso n of c a lc u la te d a n d m e a s u re d r a d ia l te m p e r a tu r e s d is trib u tio n for vario­

u s tim es; 1 - 72 s, 2 - 33 s, 3 - 600 s, 4 - 864 s, 5 - 1128 s, 6 - 1392 s; solid lin e - d ire c t n u ­ m e ric a l solution, filled s q u a re s - in v e rse so lu tio n for in p u t d a ta w ith ra n d o m e rro rs

Efektywna metoda identyfikacji nieustalonego stanu. 151

Rys. 6. W yznaczony w spółczynnik w n ik a n ia ciepła w e w n ą trz pow ierzchni zbiornika cylindry­

cznego; 1 - „dane pom iarow e” zaburzone, 2 - „dane pom iarow e” dokładne, 3 - a = 1000 W /m K - w a rto ś ć p rz y ję ta do obliczeń

Fig. 6. C alcu lated h e a t tr a n s fe r coefficient a t in sid e su rface of cylindrical vessel; 1 - d a ta w ith random erro rs, 2 - errorless d a ta , 3 - a = 1000 W /m2K - applied h e a t tr a n s fe r coefficient

Rys. 7. O siow e n a p rę ż e n ia ciep ln e n a w e w n ę trz n e j p o w ierzch n i z b io rn ik a cylindrycznego wywołane sk okow ą z m ia n ą te m p e r a tu r y czy n n ik a; 1 - n u m e ry c z n e ro z w ią z a n ie z a g a d n ie ­ nia bezpo śred n ieg o d la N = 3, 2 - n u m e ry c z n e ro z w ią z a n ie z a g a d n ie n ia b ezp o śred n ieg o d la N = 1 1 , 3 - ap ro k sy m a c ja ro z k ła d u te m p e r a tu r y p a ra b o lą , 4 - ro z w ią z a n ie z a g a d n ie n ia

odw rotnego d la N = 3

Fig. 7. A xial th e r m a l s tre s s d u e to th e r m a l sh o ck a t in s id e su rfa c e ; 1 - d ire c t n u m e ric a l solution fo r N = 3, 2 - d ire c t n u m e ric a l so lu tio n for N = 11, 3 - ap p ro x im a tio n of te m p e ra tu re

d is trib u tio n by a po lynom ial of 2 - n d d egree, 4 - so lu tio n o f in v e rse p ro b lem fo r N = 3

w ściance, ja k i współczynnik w nikania ciepła n a w ew nętrznej powierzchni zbiornika mogą być odtworzone z dużą dokładnością n a podstaw ie pomiaru tem p e ra tu ry n a zewnętrznej powierzchni zbiornika.

Dokładność proponowanej m etody je s t więc bardzo dobra, naw et w w arun­

kach analizowanego u d a ru cieplnego.

Błędy w yznaczania m aksym alnych n aprężeń cieplnych n a wewnętrznej powierzchni w alczaka (rys. 7) są nieco większe. W ynika to z faktu, że ścianka podzielona została tylko n a trzy w arstw y. Dokładność uzyskanych wyników m ożna znacznie poprawić dzieląc ściankę w alczaka n a w iększą ilość w arstw - np. cztery.

P rz y k ła d 3

Dane przyjęte do obliczeń są analogiczne ja k w przykładzie 2. Różnica polega jedynie na przyjęciu przebiegu czasowego te m p e ra tu ry czynnika (rys. 8). W ykorzystując „dane pomiarowe” zaburzone błędam i, przedstawione n a rys. 8, wyznaczono osiowe n aprężenia cieplne n a w ew nętrznej powierzchni w alczaka (rys. 9). W tym przypadku wyniki uzyskane z wykorzystaniem interpolacji kw adratow ej przebiegu tem p e ra tu ry n a grubości ścianki odzna­

czają się bardzo dobrą dokładnością. W ynika to z malej szybkości zmian tem p e ra tu ry w ściance walczaka.

Rys. 8. „D ane p o m iaro w e” zab u rzo n e b łęd am i p rz y liniow o w z ra sta ją c e j, a n a s tę p n ie stałej te m p e ra tu rz e czy n n ik a

Fig. 8. T e m p e ra tu rę d a ta w ith ra n d o m e rro rs for ra m p flu id te m p e ra tu rę

Efektywna metoda identyfikacji nieustalonego stanu. 153

-10 -30

es

Ph

-50

S

*

-70

t>

a

cCJ

-90

c *

Q.Cd

-110

C

-130 -150

czas t[s]

Rys. 9. O siow e n a p rę ż e n ia cieplne n a p o w ierzch n i w e w n ę trz n e j z b io rn ik a cylindrycznego, w ywołane z m ia n ą te m p e r a tu r y czy n n ik a , p rz e d s ta w io n ą n a ry s. 8; 1 - n u m e ry c z n e ro zw ią­

zanie z a g a d n ie n ia b ezp o śred n ieg o d la N = 3, 2 - n u m e ry c z n e ro z w ią z a n ie z a g a d n ie n ia bezpośredniego d la N = 11, 3 - a p ro k s y m a c ja ro z k ła d u te m p e r a tu r y p a ra b o lą , 4 - ro z w ią z a ­

nie z a g a d n ie n ia odw rotnego d la N = 3

Fig. 9. A xial th e r m a l s tre s s a t in s id e su rfa c e of cy lin d rical v e ss e l for flu id te m p e ra tu re changes sh o w n in fig. 8 ,1 - d ire c t n u m e ric a l so lu tio n for N = 3, 2 - d ire c t n u m e ric a l so lu tio n for N = 11, 3 - a p p ro x im a tio n of te m p e r a tu r e d is trib u tio n by a polynom ial o f 2 - n d d egree,

4 - so lu tio n of in v e rse p ro b lem for N = 3

UWAGI KOŃCOWE

Przedstaw iona w referacie num eryczno-analityczna m etoda w yznaczania rozkładu te m p e ra tu r i naprężeń cieplnych w ciśnieniowych elem entach kot­

łów i tu rb in n a podstaw ie pom iaru tem p e ra tu ry n a ich zew nętrznej powierz­

chni nadaje się zarówno do ciał o prostych, ja k i złożonych kształtach. Dzięki zastosowaniu filtru cyfrowego do w ygładzania zmierzonych przebiegów tem ­ peratury, m etoda je s t odporna n a przypadkow e błędy pom iaru. Z uw agi na prostą postać wzorów obliczeniowych m etoda je s t szczególnie p rzy d atn a do wyznaczania naprężeń cieplnych w czasie rzeczywistym , co je s t szczególnie ważne w kom puterowych układach kontroli nap rężeń cieplnych i oceny trw a ­ łości resztkowej elem entów ciśnieniowych.

Wykaz sym b oli

a; - współczynnik w yrów nania tem p eratu ry , m 2/s c - ciepło właściwe, J/kgK

E - moduł Younga, M Pa

f3 - tem p e ra tu ra zm ierzona n a powierzchni zew nętrznej, °C fi - tem p e ra tu ra w i-ty m punkcie pomiarowym, °C

L - ilość punktów pomiarowych po prawej lub lewej stronie p u n k tu zero­

wego N - ilość węzłów

qE - gęstość stru m ien ia ciepła n a powierzchni zew nętrznej, W/m2

qjN- - m aksym alna gęstość stru m ien ia ciepła przy jego zm ianie w czasie w formie trójk ąta, W/m2

qs - gęstość strum ien ia ciepła n a powierzchni w ew nętrznej, W/m2 r ; - prom ień powierzchni w ew nętrznej i-te j objętości kontrolnej, m r w - prom ień powierzchni w ew nętrznej, m

r z - prom ień powierzchni zew nętrznej, m s - w spółrzędna bezwymiarowa

t - czas, s

Tcz — te m p e ra tu ra czynnika, °C T; - tem p e ra tu ra w i-ty m węźle, °C T0 - tem p e ra tu ra początkowa, °C

Tśr - średnia tem p e ra tu ra n a grubości ścianki, °C y(tj) - w artość wygładzona w i-ty m punkcie, °C a - współczynnik w nikania ciepła, W/m2K

ß - liniowy współczynnik rozszerzalności tem peraturow ej, 1/K

X - współczynnik przewodzenia ciepła, W/mK v - liczba Poissona

© - tem p e ra tu ra bezwymiarowa p - gęstość, kg/m3

of - średnie odchylenie kw adratow e ow - osiowe n aprężenia cieplne, M Pa Ar - szerokość objętości kontrolnej, m F„ - liczba Fouriera

LITERATURA

[1] Beck J. V., Blackwell B., C lair Ch. R. St.: Inverse h e a t conduction, III - posed problems. Wiley - Interscience Publication, New York 1985.

[2] Speitkam p L.: Bestim m ung von T em peraturdifferenzen in dicken D ruckbehälterw änden aus der zeitlichen Folge von Tem peraturm ew er-

Efektywna metoda identyfikacji nieustalonego stanu. 155

te n an der isolierten W andaußenseite, VGB K raftw erkstechnik, Vol. 68, 1988, H eft 2, ss. 182-186.

[3] L eithn er R., Pich R., E rlm an K., Steege F., T rung Chi.: Vergleich ver­

schiedener V erfahren zur B estim m ung der T em peraturdifferenz in dick­

w andigen B auteilen für die Lebensdauerberechnung, VGB F achtagung

„Dampfkessel u nd D am pfkesselbetrieb”, E ssen 1989.

[4] T aler J.: N um eryczno-analityczna m etoda identyfikacji wielowymiaro­

wej pola tem p e ra tu ry w oparciu o pom iar w ew nątrz ciała, XV Zjazd Termodynamików, M ateriały konferencyjne, tom 2, Gliwice-Kokotek 1 9 9 3 ,ss. 653-659.

[5] P a ta n k a r S.: N um erical H e a t T ran sfer and F luid Flow, H em isphere P ublishing Corporation, New York 1980.

[6] R alston A.: W stęp do analizy num erycznej. PWN, W arszaw a 1975.

Recenzent: Prof. dr hab. inż. G erard KOSMAN

Wpłynęło do Redakcji 10.08.1994 r.

A b stract

The pap er p resen ts a new sem inum erical m ethod for th e solution of inverse h e a t conduction problem s. The technique developed is ap p ro p riate for one and m ultidim ensional problems.

Because th e m ethod is employed. The digital filter based on 11 m easurem ent points is used to d eterm ine tem p e ra tu re tim e curve and its derivatives w ith high accurancy.

Three exam ples dem onstrate effectivness of th e m ethod. One te s t case is for a h e a t flux th a t varies in tim e in a tria n g u la r fashion. The second case is th erm al shock in thick walled cylindrical tubes. In th e th ird exam ple fluid tem p eratu re changes are asum ed to have ram p form. T em perature and th erm al stress distributions in cylindrical vessel are determ ined based on tem p eratu re m easurem en t a t outside surface.

N um erical experim ents show th a t th e m ethod can acu rately estim ate surface h e a t flux based on tem p e ra tu re m easu rem en t on th e outside surface of th e construction elem ent.