ZESZY TY 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. 130
1998 N r kol. 1397
N ehad A L .-K H A LID Y Janusz SK O R EK
Politechnika Śląska, G liw ice
I D E N T Y F I K A C J A R O Z K Ł A D U T E M P E R A T U R Y P O D C Z A S S P A W A N I A Z A P O M O C Ą M E T O D Y F I L T R A C J I D Y N A M I C Z N E J
Streszczenie.
P raca dotyczy identyfikacji pola tem peratury w m ateriale poddaw anym spawaniu w iązką plazm y lub lasera. Do rozw iązania zadania odw rotnego w ykorzystano metodę optym alnej filtracji dynam icznej (filtr K alm ana) w połączeniu z m etodą bilansów elem entarnych. Podstaw ow ym celem była tu identyfikacja pola tem peratury w obszarze ciała stałego. P odejście to um ożliw ia uw zględnienie zm ian w łasności term ofizycznych m ateriału w raz z tem peraturą oraz pozw ala uw zględnić w pływ wielkości błędu pom iarów.W celu przetestow ania dokładności zastosow anych m etod analizow ano w pływ następujących param etrów technologicznych i obliczeniow ych: w łaściw ości term ofizycznych m ateriału, liczby i położenia punktów pom iarow ych i w ielkości błędu pom iaru. Dokładność obliczeń określano przez porów nanie uzyskanych w yników z w ynikam i analogicznego zagadnienia bezpośredniego.
IDETIFICATION OF TEM PERATURE PROFILE D U R IN G W ELDING APPLYING STO CH ASTIC APPROACH
Summary.
A K alm an filter m ethod in conjuction with the nonuniform control volume m ethodology is adopted as a num erical tool for solving tw o-dim ensional inverse heat conduction problem s for a m oving body. This w ork focuses on estim ation o f the transient tem perature in the solid region o f the w orkpiece during w elding process that involves a m oving heat source. Influence o f several param eters as therm al properties o f m aterial, location o f m easurem ent points and m easurem ent errors on the accuracy o f calculations were exam ined.6 N. A l-Khalidy, J. Skorek
1. W STĘP
O dw rotne zagadnienia przepływ u ciepła to takie problem y, w których nie są znane w szystkie przyczyny w yw ołujące dane zjawisko (np. warunki brzegowe), są natom iast znane niektóre skutki w yw ołane tym i przyczynam i (np. wartości tem peratury w niektórych punktach ciała). D o rozw iązyw ania zagadnień odw rotnych stosuje się obecnie w iele różnych m etod [1]. W ażną grupę zagadnień odw rotnych stanow ią problem y identyfikacji. Jedną z najbardziej efektyw nych m etod rozw iązyw ania problem ów identyfikacji je st m etoda filtracji dynam icznej, zw ana m eto d ą filtru Kalm ana. M etoda ta znalazła zastosow anie do rozw iązyw ania w ielu problem ów technicznych [2-5].
Z asadniczym celem przedstaw ionej pracy je st opracow anie modeli m atem atycznych i algorytm ów obliczeniow ych rozw iązania liniowego i nieliniow ego zagadnienia odwrotnego dla procesów ze zm ianą fazy (ciecz - ciało stale) oraz z ruchom ym źródłem ciepła w celu identyfikacji pola tem peratury podczas plazm ow ego i laserowego spawania metali. W celu przetestow ania dokładności rozw iązania oraz w celu sym ulacji danych pom iarow ych dla rozw iązania zagadnienia bezpośredniego [6].
2. Z A Ł O Ż E N IA M O D ELU M A TEM A TY C ZN EG O
Przedstaw iono podstaw ow e założenia i sform ułow anie m odelu m atem atycznego procesu spaw ania za p o m o c ą luku plazm ow ego. W iązka plazm y lub w iązka laserow a o odpowiedniej m ocy je st kierow ana na połączenie dw óch części spaw anego materiału. W m odelu przyjęto, że układ w spółrzędnych (x, y) porusza się w raz ze źródłem energii (ciepła), tzn. wraz z w iązką plazm y (rys. 1). Z akłada się, że energia prom ieniow ania wiązki je st całkowicie absorbow ana przez m ateriał w tzw. „oczku” (w technologii spawalniczej używ a się również określenia „dziurka do klucza” - ang. K eyhole). Ponadto przyjm uje się (zgodnie z danymi eksperym entalnym i) kołow y kształt oczka.
R ów nania energii dla dw uw ym iarow ego obszaru ze zm ianą fazy R sform ułowano oddzielnie dla obszaru fazy ciekłej R] i fazy stałej R s. O bydw a obszary są oddzielone p ow ierzchnią rozdziału faz 91. R ów nania nieustalonego pola tem peratury m ają następującą postać:
Identyfikacja rozkładu tem peratury. 7
dla fazy stałej
c s
+ p C s(T )u d! ^ X' y ^ = VXS (T)VTS (x, y ,t) -
^ T (x , y ,t)
( 1)oraz dla fazy ciekłej
aT,(x,y,t ) + gr,(x^y,t) = VX(VT) (X;y> t) _ h £ ) T(Xi y> t) (x> y) e Ri
5t 5x 6
(2 )
gdzie: p - gęstość, c - ciepło w łaściw e, u - prędkość spawania, T - tem peratura liczona względem tem peratury otoczenia, X - przew odność cieplna, h - w spółczynnik wnikania ciepła, 5 - grubość płyty, t - czas. Indeksy s i 1 odnoszą się do fazy stałej i ciekłej. W rów naniu dla fazy ciekłej pom inięto człon zw iązany z konwekcją.
R ys. 1. S p a w a n ie p ły t m e ta lo w y c h lu k ie m p la z m o w y m Fig. 1. W e ld in g o f th e s te e l p la te s u s in g p la sm a arc
W spółczynnik w nikania ciepła na górnej i dolnej pow ierzchni płyty w ynika z zależności:
h(T) = hc + hr (3)
8 N. A l-Khalidy, J. Skorek
Zastosow anie techniki odwrotnej pozw ala pom inąć w rozw ażaniach w ym ienione ciepła w obszarze fazy ciekłej. Stąd rów nanie (2) w prow adzono do m odelu bezpośredniego tylko w celu sym ulacji danych w ejściow ych do testow ania m odelu zagadnienia odw rotnego. Wyniki rozw iązania zagadnienia bezpośredniego przedstaw iono w pracy [6].
W m odelu przyjęto następujący w arunek brzegowy na granicy fazy ciekłej i oczka
T = T b - T „ a t r arc= ^ x2+ y 2 , (4)
gdzie Tb oznacza tem peraturę parow ania m ateriału a rarc je st prom ieniem oczka.
W arunek sym etrii w zględem osi x m a postać:
y r per
— 1 = 0; -—*■ = 0 dla y = 0 (5)
dy dy
N a kraw ędziach płyty m etalu przyjęto warunki brzegowe drugiego rodzaju (powierzchnie adiabatyczne)
^ = 0 — = 0 (6)
dx dy
Do dyskretyzacji problem u zastosow ano m etodę bilansów elem entarnych [7]. Ważną zaletą tej m etody je st to, że całkow e równania zachow ania energii są spełnione autom atycznie zarówno w obrębie elem entu różnicow ego, ja k i w całym obszarze. W celu uniknięcia oscylacji rozw iązania (ze w zględu na obecność członu konw ekcyjnego w rów naniu energii) zastosow ano tzw. technikę p o d wiatr.
Jak ju ż w spom niano, rozw iązanie zagadnienia bezpośredniego Tex stanowiło podstaw ę do sym ulacji w yników pom iarów T mea, w ykorzystyw anych następnie do rozwiązania zagadnienia odw rotnego:
T(xm, ym, tk) = T ex(m) m = 1 ,2 ,..., M , (7)
gdzie M oznacza ca łkow itą liczbę punktów pom iarow ych, a k je st num erem kroku czasu.
Identyfikacja rozkładu tem peratury... 9
Do rozw iązania zagadnienia nie je st potrzebny w arunrk początkow y, co je st zaletą metody filtru K alm ana.
3. R O Z W IĄ Z A N IE ZA G A D N IE N IA O D W RO TN EGO Z W Y K O R ZY STA N IEM M ETO D Y O PT Y M A L N E J FILTRA CJI D Y NA M ICZN EJ
P rzedstaw iono zastosow anie m etody optym alnej filtracji dynam icznej do identyfikacji nieustalonego pola tem peratury w obszarze stałym spawanego m ateriału. Przyjęto, Ze w łaściw ości m ateriału są stałe lub są funkcjam i tem peratury. Do identyfikacji w ykorzystuje się w yniki pom iarów tem peratury w punktach usytuow anych w obszarze fazy stałej. Wyniki pom iarów są sym ulow ane poprzez rozw iązanie zagadnienia bezpośredniego [6], Analiza odw rotna dotyczy identyfikacji tem peratury w ew nątrz i na brzegu rozpatryw anego obszaru.
W m etodzie filtru K alm ana estym acji (rekursyw nej) podlegają zarówno wartości identyfikow anych w ielkości, ja k i ich wariancje. Jest to w ażna cecha m etody, um ożliw iająca ocenę w iarygodności uzyskanych rezultatów. M etoda składa się z dw óch etapów:
- predykcji, tzn. oceny pola tem peratury na bazie m odelu m atem atycznego, filtracji, tzn. korekcji rezultatów na podstaw ie w yników pom iarów.
W celu zastosow ania m etody filtracji m odel zjaw iska pow inien być zapisany w standardow ej postaci tzw. rów nania stanu (w ogólności nieliniowego):
T em peratury w ęzłow e (zgrupow ane w w ektorze stanu T) s ą traktow ane ja k wielkości stochastyczne. Z akłada się, że relacja pom iędzy w ielkościam i Zk+i zm ierzonym i w chwili czasu k+ 1 a tem peraturam i w ęzłow ym i Tk+i m a następującą postać (tzw. rów nanie wyjść):
gdzie H je s t m acierzą w yrażającą relacje pom iędzy w ielkościam i m ierzonym i a identyfikow anym i, a V|<+i oznacza w ektor błędów pom iarow ych (o zerowej wartości średniej i o m acierzy kow ariancji V).
Tk+i - Fk+i (T k ) (8)
Zk+i - H Tk+i + Vk+i, (9)
10 N. A l-Khalidy, J. Skorek
W skrócie algorytm m etody przedstaw ia się następująco:
- Szacuje się w artości tem peratury w w ęzłach
Tk+1k
w chwili czasu k+1 na podstaw ieoszacow ania tem peratury
Tk
i macierzy kow ariancji W k w poprzednim kroku czasu k. Do tego celu w ykorzystuje się równanie stanu (8) oraz prawo propagacji błędów.- W chw ili czasu k+1 dokonuje się pom iarów tem peratury w w ybranych punktach obszaru.
- O szacow anie w pierw szym kroku w artości koryguje się na podstaw ie rów nania filtru:
t kłl= f kł, k+Kkt,[zkłl- H f kłlJ ,
( 1 0 )gdzie m acierz K k+i oznacza tzw. m acierz w zm ocnienia filtru.
K k+, = W k ł l HTVk;1l; ( l l )
gdzie indeks T oznacza transpozycję macierzy.
M acierz kow ariancji a posteriori w yników estym acji W k+i m a następującą postać:
W
k+1
= W k+1. k - W k+1
. kH T( H W k+l. kH T + V k+i) 'lH W k+l. k (12)D la zagadnienia nieliniow ego zastosow ano procedurę iteracyjną do korygow ania zależnych od tem peratury w łaściw ości term ofizycznych materiału.
4. W Y N IK I O B L IC Z E Ń
Zaproponow ana m etoda została zastosow ana do identyfikacji pola tem peratury w ew nątrz i na brzegu spawanej płyty. Rozpatrzono problem dw uw ym iarow ego przew odzenia ciepła.
W łaściw ości term ofizyczne są funkcjam i tem peratury (zagadnienie nieliniow e). Prędkość źródła przyjęto ró w n ą 1 mm/s. W yniki pom iarów sym ulow ano w ykorzystując rozwiązanie zagadnienia bezpośredniego Tex. W artość tę zaburzano następnie losow ym błędem zgodnie z zależnością:
Tcr(i, k) = T mea(i, k) + (l - 2t)w, (13) gdzie w oznacza m aksym alną w artość błędu, a x oznacza liczbę losow ą z przedziału [0, 1],
Identyfikacja rozkładu tem peratury. 11
P rzyjęto następujące podstaw ow e dane do obliczeń:
N = 20 x 15 = 300 elem entów różnicow ych, m = 1 1,
X = 14,42 + 0,0169*T - 2,44E-06*T 2,
c = 484,6 + 0,159*T + 18,07E-06*T2.
At = 0,08 s.
Przyjęto, że punkty pom iarow e są usytuow ane na linii o współrzędnej y = 8 mm.
1200.00
1 0 0 0 .0 0
800.00
o c iIh 3 600.0003
t-<
<D O- i 1 0 0 .0 0 H
2 0 0 .0 0
0. 00
0.00 0.02 0. ¡34 0.(36 0.(38 0.'l0 O dległość, m
R y s. 2 . R o z k ła d te m p e ra tu ry w z d łu ż o si x (u = 1 m m /s, t = 12 s, b łąd p o m ia ru = 0 ) F ig. 2. T e m p e ra tu rę d is trib u tio n a lo n g (u = 1 m m /s, t = 12 s, m e a su re m e n t e rro r = 0 )
W yniki identyfikacji pola tem peratury przedstaw iono na rysunkach od 2 do 5. Pole tem peratury w zdłuż osi x dla przypadku braku błędów pom iarow ych przedstaw iono na rys. 2 i rys. 3. W idoczna je s t bardzo dobra zgodność z rozw iązaniem dokładnym . Jest ona szczególnie w yso k a w obszarze, gdzie w ystępują w ysokie gradienty tem peratury. D okładność z m niejsza się w obszarach o w zględnie m ałych gradientach tem peratury. W idoczne jest, że stan pseudoustalony osiągany je st najszybciej w pobliżu działania źródła (wiązki plazm y). Po osiągnięciu pseudoustalonego rozkładu tem peratury w m ateriale dokładność identyfikacji staje się bardzo w ysoka (rys.3).
*
X u
«
** 4 *
X
* * 4X 4 K 4
* *• łx *
* X
,K - x ‘ 4 4X *
* 4 X x ^ X
*X 4 X +
► *y X * 4 4 x X
*x* * X**■♦5 * 4
: 5 :
* * . * ¥ 5. l i i i x i s ł i i l l i
■* x l
5
i I 1
N12 N. A l-K halidy, J. Skorek
R y s. 3. R o z k ła d te m p e ra tu ry w z d łu ż osi x (u = 1 m m /s, t = 6 0 s, b łą d p o m ia ru = 0 ) Fig. 3. T e m p e ra tu re d is trib u tio n a lo n g x ax is (u = 1 m m /s, t = 6 0 s, m e a su re m e n t e rro r = 0)
4 0 .0 0 60 .0 0
Czas, s
R y s. 4. P rz e b ie g te m p e ra tu ry w w ę ź le n r 5 2 (i = 3, j = 12) Fig. 4. T e m p e ra tu re d is trib u tio n at th e n o d e 52 (i = 3, j = 12)
Identyfikacja rozkładu tem peratury. 13
A nalizow ano rów nież w pływ błędów pom iarów na dokładność estym acji tem peratury. Na rys. 4 przedstaw iono w yniki obliczeń dla błędu w zględnego, w ynoszącego 5%. Dokładność estymacji je st ciągle bardzo dobra (szczególnie po osiągnięciu stanu ustalonego). Błąd względny identyfikacji tem peratury (Tex - T inv)/Tcx w funkcji czasu dla w ybranych węzłów oraz dla różnych w artości błędu pom iaru przedstaw iono na rys. 4. Błąd pom iarów przyjęto 0% (linia ciągła) oraz 5% (linia przeryw ana).
R y s. 5. W z g lę d n y b łąd e s ty m a c ji te m p e ra tu ry w p ły c ie (1 - x = 0,031 m , y = 0 ,0 1 2 5 m ; 2 - x = 0,0 5 m , y “ 0 ,0 0 8 7 5 m )
Fig. 5. R e la tiv e e rro r o f te m p e ra tu re e stim a tio n in th e p la te (1 - x = 0,031 m , y = 0 ,0 1 2 5 m , 2 - x = 0,0 5 m , y = 0 ,0 0 8 7 5 m )
5. UW AGI K O Ń C O W E I W NIOSK I
1. Z astosow anie techniki odw rotnej je s t korzystnym narzędziem analizow ania procesów spaw ania, g dyż um ożliw ia opis bardzo złożonych procesów cieplnych bez konieczności rozpatryw ania trudnych do opisu m atem atycznego procesów zachodzących w obszarze fazy ciekłej.
14 N. A l-K halidy, J. Skorek
2. A nalizow ano w pływ szeregu param etrów (liczba punktów pom iarow ych, położenie punktów pom iarow ych, prędkość spawania, moc wiązki plazm owej, dokładność pom iarów ) na rozkład tem peratury w m ateriale i dokładność obliczeń. Uzyskano zadow alające rezultaty, szczególnie w przypadku, gdy punkty pom iarow e są położone stosunkow o blisko granicy faz.
3. O pracow any m odel, oparty na m etodzie optymalnej filtracji dynam icznej (filtr Kalmana).
um ożliw ia identyfikację dw uw ym iarow ego nieustalonego pola tem peratury w spawanym m ateriale. D okładność m odelu była zw eryfikow ana za pom ocą eksperym entów num erycznych. D okładność rozw iązania zależy w głównej m ierze od następujących param etrów : liczby i rozm ieszczenia punktów pom iarow ych, błędu pom iaru, prędkości spaw ania, w łaściw ości term ofizycznych m ateriału oraz rozm iarów elementów różnicow ych i kroku czasu.
4. Testy num eryczne potw ierdziły przydatność m odelu filtracji dynam icznej do analizy problem u spaw ania plazm ow ego. M odel m oże być rów nież w ykorzystany do analizy' innych problem ów ze zm ianą fazy. M ożliw a w ydaje się także m odyfikacja m odelu w celu identyfikacji położenia granicy faz.
LITER A TU R A
1. Beck J. B lackw ell B. and Clair St. C.R., Inverse H eat Conduction- III Posed problems:
W iley-Tnterscience, N ew York, 1985.
2. N orton J.: A n introduction to Identification, Academ ic Press: London, 1986.
3. M urakam i A. and H asegaw a T., A plication o f Kalm an filtering to inverse problems:
Theoretical and A pplied M echanics, 1993,42, pp. 3-14.
4. M urakam i A. and Ham aguchi T., Identification o f param eters o f groundw ater flow model by using eztended K alm an filter-FE: Proceeding o f the 2nd International Sym posium on Inverse P roblem -ISIP, Paris/France, 1994.
5. A l.-K halidy N. and Skorek J., Optymal Dynamic Filtration Approach for Inverse Heat C onduction P roblem s w ith M oving Body: Inverse Problem s in Engineering, 1997, 4, pp.
209-229, USA.
Identyfikacja rozkładu tem peratury. 15
6. A l.-K halidy N ., Enthalpy technique for solution o f Stefan problem s: A pplication to the keyhole plasm a arc w elding process involving m oving heat source: Int. Com m . Heat M ass transfer, 1995, 22, pp. 779-790.
7. P atankar S., N um erical H eat Transfer and Fluid Flow, Hem isphere, W ashington, 1980.
Recenzent: Prof.dr hab. Rom an Domański
W płynęło do R edakcji 03.09.1998 r.
A b stra c t
The inverse problem w ith is considered herein deals w ith estim ation o f the transient tem perature in teh solid region o f a w orkpiece during the w elding process based on the transit tem perature m easurem ents at a num ber o f therm ocouples located in the solid region. The fusion zone not analyzed. The problem belongs to the class o f boundary inverse heat conduction problem s. In practice, such a problem is o f great interest w hen only the m easurem ent o f the tem peratures at a num ber o f sensors enables us to identify the tem perature fiels w ithin the calculation dom ain. In such a case heat transfer and fluid flow equations in the fusion zone need not be considered. In this study the tem perature data are sim ulated basing on the solution o f the direct problem . The solution o f the inverse problem is based on the theory o f the optim aldynam ic filtration m ethod. Basing on the num erical results the follow ing conclusions are sum m arized.
1. The optim al dynam ic filtration m ethod gives very good and stable results o f calculation (The accuracy o f the num erical solution for inverse problem is verified by com paring the results w ith the reference direct solution).
2. The presented m ethod perm its treatm ent o f a tem perature dependent therm al properties o f m aterials and disturbances in the m easurem ent data.
3. One o f the advantages o f the algorithm is that the initial tem perature can be unknow n (the initial inverse problem is o f significant interesting in m any engineering application).
4. The resukts can be very im portant for determ ination the solid structure during w elding and heat treatm ent process.
