MODELOWANIE INŻYNIERSKIE ISSN 1896-771X 44, s. 285-292, Gliwice 2012
METODA PURC W ANALIZIE NIEUSTALONEGO POLA TEMPERATURY W OBSZARACH PŁASKICH
EUGENIUSZ ZIENIUKa,DOMINIK SAWICKIb
aInstytut Informatyki, Uniwersytet w Białymstoku e-mail: ezieniuk@ii.uwb.edu.pl,
bWydział Mechaniczny, Politechnika Białostocka e-mail: sawicki.dominik1@gmail.com
Streszczenie. Obecność dyskretyzacji w klasycznej MES oraz MEB jest dość istotnym mankamentem. Alternatywą pozwalającą na uniknięcie wspomnianego problemu są parametryczne układy równań całkowych (PURC) niewymagające klasycznej dyskretyzacji podczas ich numerycznego rozwiązywania. Celem niniejszej pracy jest uogólnienie metody PURC i przedstawienie możliwości jej zastosowania do modelowania i symulacji zagadnień brzegowo-początkowych na przykładach dotyczących problemów temperaturowych.
1. WSTĘP
Modelowanie i symulacja zagadnień brzegowych i brzegowo-początkowych niezmiennie przyciąga zainteresowanie zarówno z naukowego, jak również i praktycznego punktu widzenia. Modelowanie tych zagadnień za pomocą równań różniczkowych czy też całkowych jest znane w literaturze od dawna. Znane są też różne metody ich analitycznego, jak również numerycznego rozwiązywania [1,2]. Metody analityczne są stosowane raczej do zagadnień zdefiniowanych w obszarach o elementarnych kształtach oraz z zadanymi elementarnymi warunkami brzegowymi. Przy takich założeniach są one z powodzeniem stosowane raczej tylko do zagadnień brzegowych.
W zagadnieniach brzegowo-początkowych dodatkowym problemem (w porównaniu z zagadnieniami potencjalnymi) jest to, że rozwiązanie zagadnienia w obszarze czy też tylko na jego brzegu jest zależne od dodatkowej zmiennej, którą jest czas. Analityczne rozwiązywanie takich zagadnień dla niektórych zagadnień brzegowo-początkowych jest również możliwe, ale najczęściej jest ono otrzymywane w postaci wyrażeń przybliżających, będących najczęściej szeregami możliwymi do otrzymania jedynie dla elementarnych kształtów obszaru [6,7]. W przypadku bardziej złożonych zagadnień nieustalonych są stosowane metody komputerowe [1].
Do najczęściej wykorzystywanych metod komputerowych, podobnie jak w zagadnieniach ustalonych, należą: metoda elementów skończonych (MES) [3] i metoda elementów brzegowych (MEB) [1,2]. Metody te w rozwiązywaniu zagadnień nieustalonych wymagają stosowania procesów iteracyjnych. O efektywności rozwiązywania takich zagadnień decyduje efektywność metody w poszczególnych krokach czasowych. Występowanie klasycznej dyskretyzacji we wspomnianych metodach jest dość istotnym mankamentem. Alternatywą
pozwalającą na uniknięcie wspomnianego problemu są parametryczne układy równań całkowych (PURC) eliminujące klasyczną dyskretyzację brzegu. PURC z dużą efektywnością były stosowane dla zagadnień potencjalnych definiowanych różnymi równaniami różniczkowymi [4].
Celem niniejszej pracy jest uogólnienie metody PURC i przedstawienie możliwości jej zastosowania do modelowania i symulacji zagadnień brzegowo-początkowych na przykładach dotyczących problemów temperaturowych mających rozwiązania dokładne.
2. PURC DLA ZAGADNIEŃ NIEUSTALONYCH
Równanie różniczkowe Fouriera opisujące nieustalone pole temperatury bez wewnętrznych źródeł jest przedstawiane w postaci [1]
), ) (
( 2
,t u t k
,t
c u x x
xΩ (1)
gdzie k [W/mK] jest to współczynnik przewodzenia ciepła, u(x,t) pole temperatury zależne od czasu, c[J/m3K] jest to ciepło właściwe odniesione do jednostki objętości, natomiast t jest to czas, zaś
). ( )
) ( ,
( 2
2 2
2 1 2 2
x ,t u x
,t t u
u
x x
x
Klasyczne brzegowe równanie całkowe (BRC) dla równania różniczkowego (1) jest przedstawiane w postaci [1]
), ( ) ( ) (
) ( ) 1 (
) ( ) 1 (
) , ( 5 0
0 0
*
0 0
x x
x ξ
x x
ξ x
x ξ ξ
dΩ ,t u ,t ,t , U
dt d ,t u ,t ,t , c P
dt dΓ ,t p ,t ,t , c U
t u .
Ω
F
t
t Γ
F t
t Γ
F F
F F
(2)
gdzie
) , ( exp 4 ) ( 4 ) 1 (
2
*
t t a
r t
t ,t πa
,t ,
U ξx F F F r2 (x11)2 (x22)2,
) , ( exp 4 ) ( ) 8
(
2 2
2
*
t t a
r t
t πa ,t kd ,t ,
P ξx F F F d
x1ξ1
n1
x2ξ2
n2.Po zastosowaniu analitycznej modyfikacji klasycznych BRC analogicznej, jaka była stosowana w przypadku zagadnień ustalonych [4], otrzymano PURC dla zagadnień nieustalonych, który jest przedstawiony za pomocą wyrażenia
).
( ) ( ) (
) ( ) (
) ( ) 1 (
) ( 5 0
0 0
1
* 1
1 1
1
0 1
y y
y,t ,t u ,t dΩ ,
s U
dt ds s J s,t u ,t ,s,t s P s,t p ,t ,s,t s c U
,t s u .
Ω
F l
t
t
j n
j s
s
j F lj
j F lj
F l
F j
j
(3)
Funkcje podcałkowe w (3) są przedstawiane w następującej postaci ) ,
( exp 4 ) ( 4 ) 1 , (
2 1
*
t t a t
t ,t πa
s,t s
Ulj F F F
) , ( )
) ( (
5 . ) 0 2 ( )
1 (
s s S s
s s S
Jj j j
) , ( exp 4 ) ( ) 8
, (
2 2
1 2
*
t t a t
t πa
d ,t k
s,t s
Plj F F F
2,
2 1
1n n
d
gdzie
2,
2 2 1
2
1 Sl(1)(s1)S(j1)(s),2 Sl(2)(s1)S(j2)(s), Si(s) [Si(1),Si(2)]T. )
i(s
S jest krzywą opisującą ity segment brzegu. Funkcja podcałkowa w drugiej całce (3) jest przedstawiana w postaci
) , (
exp 4 ) (
4 ) 1 ,
( 0
2 2
0 0
1
*
t t a t
t ,t πa
,t s
Ul F F F
y
gdzie 2 1222, 1 Sl(1)(s1) y1, 2 Sl(2)(s1)y2.
Na podstawie (3) można otrzymać tylko rozwiązania na brzegu rozpatrywanego zagadnienia brzegowo-początkowego. Podobnie jak w przypadku klasycznych BRC, mając rozwiązania na brzegu, można otrzymać rozwiązania w obszarze na podstawie tożsamości całkowej. Po otrzymaniu rozwiązań na brzegu na bazie równ. (3), rozwiązania w obszarze otrzyma się na podstawie zmodyfikowanej klasycznej tożsamości całkowej dla rozwiązań w obszarze. Przy podobnym postępowaniu jak w [4] zmodyfikowana tożsamość całkowa dla zagadnień nieustalonych przyjmuje następującą postać
).
( ) ( ) ˆ (
) ( ) ˆ (
) ( ) , ˆ ( ) 1
(
0 0
* 1
0 1
y y
y x
x x
x
dΩ ,t u ,t ,t , U
dt ds s J s,t u ,t ,s,t P s,t p ,t t ,s c U
,t u
Ω
F t
t
j n
j s
s
j F j
j F j
F
F j
j
(4)
Funkcje podcałkowe w pierwszej całce są przedstawiane w postaci ) ,
( exp 4 ) ( 4 ) 1
ˆ *( 2
t t a
r t
t ,t πa
,s,t
Uj F F F
x r2 r12 r22,
) , ( exp 4 ) ( ) 8
,
ˆ ( 2
2 2
*
t t a
r t
t πa
d ,t k
s,t
Pj F F F
x d r1n1 r2n2,
gdzie r1 x1S(j1)(s),
).
2 (
2
2 x S s
r (j )
Natomiast funkcja podcałkowa w drugiej całce (4) za pomocą
) (
4 ˆ )exp
( 4 ) 1 ,
ˆ (
0 2 0
0
*
t t a
r t
t ,t πa
,t
U F F F
y
x ˆ ˆ ˆ2,
2 2 1
2 r r
r
gdzie ˆ ,
1 1
1 x y
r
ˆ .
2 2
2 x y
r
3. ROZWIĄZYWANIE PURC DLA ZAGADNIEŃ NIEUSTALONYCH
Po zastosowaniu dyskretyzacji po zmiennej czasowej t , przy założeniu stałego kroku czasowego ttf tf1 PURC przedstawiony wzorem (3) sprowadza się do następującej postaci
), ( ) ( ) (
) ( ) ( ) (
) ( ) 1 (
) ( 5 0
1 1
1
* 1
1 1
1
1 1
y y
y,t ,t u ,t dΩ ,
s U
dt ds s J s,t u ,t ,s,t s P s,t p ,t ,s,t s c U
,t s u .
Ω
f f
f l
t
t
j n
j s
s
j f lj
j f lj
f l
f
f j
j
(5)
po scałkowaniu funkcji podcałkowych po zmiennej t , równania (5) przyjmują następującą postać
. d t , u , s U
ds s J t s, u s , s P t s, p s , s U ,t
s 0.5u
1 - f Ω
1 l
j n
1 j
s
s
f j 1 lj f j 1 lj f
l
j
1 j
) ( ) ( ) (
) ( ) ( ) ( ) ( ) ( )
(
* 1
y y
y
(6)
Funkcje podcałkowe Ulj,Plj,Ul zależą teraz od t i przedstawiane są w następującej postaci
,) ( 1
* i
lj E
4kπ s 1
, s
U ( ) ,
2
1 2
* 4a t
lj e
η 2 s d , s
P
l*( 1 ) 4a t,
2
te a 4 s 1 , s
U
gdzie k jest współczynnikiem przewodzenia ciepła, natomiast Ei
eksponencjalną funkcją całkową, przybliżoną następującym wzorem
n!
n 1 μ μ
ln γ E
n
1 n
1 n
i
gdzie ,
t 4a μ η
2
oraz 0,5772156649 jest stałą Eulera, a=k/c jest współczynnikiem dyfuzji, c ciepłem właściwym a t jest krokiem czasowym.
Tożsamość całkowa (4) dla rozwiązań w obszarze po scałkowaniu ma postać
U ,s p s,t P ,s u s,t J s ds U ,y u y,t d .) t ,
u( f-1
Ω j
n
1 j
s
s
f j j
f j j
f
j
1 j
) ( ) ( ) ( )
( ) ( ) ( ) ( )
(x x * x y
x
(7)
Funkcje podcałkowe Uj,Pj,U są przedstawiane w postaci następującej
, )*( i
j E
4kπ s 1 , x
U ,
2
t 4a
r
oraz
, )
(
2
2
* 4a t
r
j e
r 2 s d , x
P
( ) .
ˆ
* 4at
r2
te a 4 s 1 , x
U
4. ITERACYJNE ROZWIĄZYWANIE ZAGADNIEŃ
Rozwiązanie równania przedstawionego wzorem (5) wymaga procesu iteracyjnego.
Ogólny algorytm postępowania w przypadku tego procesu iteracyjnego przedstawiono schematycznie na rys.1
Rys. 1. Schemat procesu iteracyjnego.
Otrzymanie rozwiązania na brzegu za pomocą równania (5) oraz rozwiązań w obszarze za pomocą tożsamości całkowej (7) wymaga obliczania całek po obszarze . W klasycznej MEB całkowanie sprowadza się do podzielenia obszaru całkowania na komórki.
W PURC całkowanie po obszarze charakteryzuje się tym, że obszar jest traktowany w sposób globalny bez dzielenia na komórki. Zamiast stosowania kwadratur niższego rzędu na poszczególnych komórkach, jak jest stosowane w klasycznej MEB, zastosowano kwadraturę wyższego rzędu z dużą liczbą współczynników dla całego obszaru [5]. Do obliczenia całki w obszarze zastosowano technikę opartą na globalnym traktowaniu obszaru, bez konieczności dzielenia go na komórki. Do modelowania obszaru w sposób globalny zastosowano płaty powierzchniowe, natomiast do całkowania numerycznego zastosowano kwadratury wyższych rzędów.
Do zamodelowania obszaru całkowania w PURC oraz tożsamości całkowej wykorzystano płaty powierzchni Béziera stopnia 3. Przykładowe płaty zostały przedstawione na rys. 2a,b.
a) modelowanie b) modyfikacja
Rys. 2. Modelowanie i modyfikowanie obszarów do całkowania.
Jak pokazano na rys. 2b, kształt płata z rys. 2a można łatwo zmodyfikować poprzez przesunięcie jego niektórych punktów. Do numerycznego rozwiązania PURC zastosowano wcześniej testowaną na zagadnieniach brzegowych metodę kolokacji [4]. Zbadano wpływ liczby punktów kolokacji oraz liczby współczynników kwadratury całkowania numerycznego po obszarze na dokładność i stabilność otrzymywanych rozwiązań.
5. ANALIZA ROZWIĄZAŃ
Rozpatrywano pole temperatury w prostokątnym obszarze z kx1 kx2 1 o wymiarach
x1
L (0 x11) oraz
x2
L (0 x21). Początkowa wartość temperatury w obszarze wynosi
x ,x ,0
0Φ 1 2 . Warunki brzegowe są pokazane na rys. 3a i są przedstawiane w następującej postaci
1, 2,
1, 2,
1.0, Lx x t x Lx t
, , 0 0 , ,
,
0 2 1
dn t x d dn
t x
d
x1,x2,0
0.Rozwiązanie analityczne dla tak zdefiniowanego zagadnienia jest przedstawione w [6].
Problem ten był rozwiązywany numerycznie za pomocą omówionego w punkcie 3. algorytmu iteracyjnego z wykorzystaniem PURC w każdym kroku iteracyjnym. Do zamodelowania obszaru pokazanego na rys. 3a wykorzystano pojedynczy płat Béziera przedstawiony na rys.
2a. W celu jego zdefiniowania zadano tylko punkty brzegowe. Do numerycznego rozwiązywania PURC zastosowano metodę kolokacji, testowano wpływ liczby punktów kolokacji na otrzymywane wyniki jak również wpływ liczby współczynników kwadratury całkowania po obszarze.
a) kwadratowy b) okrągły
Rys. 3. Kształty obszarów z warunkami brzegowymi.
W tabl. 1 przedstawiono wpływ liczby współczynników kwadratury na średni błąd względny rozwiązań. Na podstawie zamieszczonej tabeli okazało się, że najmniejszy błąd względny uzyskano z zastosowaniem 40 współczynników w kwadraturze. Zwiększenie liczby współczynników z racji ich koncentracji głównie w bliskiej odległości od brzegu powoduje wzrost błędu rozwiązań. W następnej kolejności przeprowadzono testy dotyczące wpływu liczby punktów kolokacji na dokładność wyników przy 40 współczynnikach w kwadraturze zastosowanej do całkowania po obszarze. Okazało się, że rozwiązania ustabilizowały się przy 6 punktach kolokacji na poszczególnych segmentach boków.
Tabl. 1. Wpływ liczby współczynników (W) kwadratury przy 6 punktach kolokacji t Średni błąd względny [%]
20
W W 40 W 60
0,1 12,24727 12,24727 12,29615 0,2 2,114546 1,767501 1,632069 0,3 6,421704 1,408876 2,035915 0,4 8,385672 1,590537 2,49654 0,5 9,238565 1,283754 2,358713 0,6 9,621399 0,845153 2,04271 0,7 9,79113 0,444005 1,725105 0,8 9,861582 0,164877 1,461469 0,9 9,886179 0,235484 1,260587 1,0 9,890573 0,392009 1,114903
W kolejnym przykładzie rozpatrywano podgrzewanie nieskończonego walca (rys. 3b) o promieniu podstawy r=r0. Rozwiązanie analityczne dla takiego problemu jest przedstawiane w pracy [7]. Początkowa wartość temperatury w obszarze wynosi r
,0 0,0. Na powierzchni walca przyłożona jest stała temperatura o wartości
r0,t
1,0. Do zamodelowania obszaru wykorzystano płat Béziera z poprzedniego przykładu. Kształt płata został odpowiednio zmodyfikowany poprzez rozciągnięcie jego krawędzi, jak pokazano na rys. 2b. Zbadano wpływ liczby punktów kolokacji oraz liczby współczynników kwadratury całkowania numerycznego na otrzymywane wyniki. Wyniki przedstawiono w tabeli 2.Tab. 2. Wpływ liczby współczynników (W) kwadratury przy 6 punktach kolokacji.
t Średni błąd względny [%]
20
W W 40 W 60
0,1 12,72546 12,72546 12,72546 0,2 2,961319 4,665196 4,39615 0,3 0,629746 2,332742 1,982242 0,4 0,835457 1,37029 0,987728 0,5 1,355849 0,89664 0,498934 0,6 1,631038 0,642592 0,237129 0,7 1,783508 0,50039 0,090796 0,8 1,869999 0,419064 0,0158 0,9 1,919649 0,372048 0,041021 1,0 1,948321 0,344728 0,069021
Na podstawie przeprowadzonych testów okazuje się, że rozwiązania ustabilizowały się przy 6 punktach kolokacji i 60 współczynnikach w kwadraturze. Zwiększona liczba współczynników w porównaniu do wcześniejszego przykładu podyktowana jest modyfikacją obszaru i zwiększeniem jego pola powierzchni. Ważną zaletą tej strategii jest to, że modyfikowanie obszarów jest efektywne i nie wymaga jego dzielenia na komórki.
6. WNIOSKI
W pracy przedstawiono PURC dla zagadnień brzegowo-początkowych (3D) modelowanych równaniem różniczkowym Fouriera opisującym nieustalone pole temperatury.
Przedstawiono i zweryfikowano strategię procesu iteracyjnego z wykorzystaniem PURC na poszczególnych krokach tego procesu. Pokazano efektywność strategii globalnego całkowania w procesie iteracyjnym stosowanym w rozwiązywaniu zagadnień nieustalonych.
Testowano wpływ liczby współczynników w kwadraturze wyższego rzędu, zastosowanej do numerycznego całkowania oraz liczby punktów kolokacji na dokładność uzyskanych rozwiązań za pomocą PURC. Otrzymane wyniki dla zamieszczonych przykładów potwierdzają efektywność i wiarygodność stosowanej metody.
Praca naukowa finansowana ze środków na naukę w latach 2010-2013 jako projekt badawczy Nr N N519579538.
LITERATURA
[1] Majchrzak E.: Metoda elementów brzegowych w przepływie ciepła. Częstochowa: Wyd.
Pol. Częst. , 2001.
[2] Brebbia C.A., Telles J.C, Wrobel L.C: Boundary element techniques, theory and applications in engineering. New York: Springer, 1984.
[3] Zienkiewicz O.C., Taylor R.L.: The finite element method, Vol. 1-3. Oxford: Butterworth, 2000.
[4] Zieniuk E.: Bézier curves in the modification of boundary integral equations (BIE) for potential boundary-values problems. “International Journal of Solids and Structures”
2003, 9(40), p. 2301-2320.
[5] Bołtuć A., Zieniuk E.: Modeling domains using Bézier surfaces in plane boundary problems defined by the Navier-Lame equation with body forces. “Engineering Analysis with Boundary Elements” 2011, 35, p. 1116-1122.
[6] Agnieszka Fraska: Wyznaczanie niestacjonarnych pól temperatury – porównanie metod numerycznych w obszarach 2D. Zesz. Nauk. Pol. Pozn. nr 2 „Budowa Maszyn i Zarządzanie Produkcją” 2005, s. 5-16.
[7] Alok Sutradhar, Glaucio H. Paulino, L. J. Gray: Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method. “Engineering Analysis with Boundary Elements” 2002, 26, p.
119-132.
PIES METHOD IN ANALYSIS OF TRANSIENT TEMPERATURE DISTRIBUTION IN FLAT AREAS
Summary. The occurrence of discretization in classical FEM and BEM is a quite essential disadvantage. An alternative to avoid the problem are parametrical integral equations systems (PIES) that do not require the classical discretization while solving them numerically. The purpose of this paper is to generalize the PIES method and present its capabilities in application to modelling and simulation of initial-boundary value problems for transient heat conduction.