NUMERYCZNA WIZUALIZACJA LINII PRZEPŁYWU CIEPŁA I GĘSTOŚCI STRUMIENIA CIEPŁA W PŁASKIM PRZEWODZENIU CIEPŁA METODĄ ELEMENTÓW BRZEGOWYCH

NUMERYCZNA WIZUALIZACJA LINII PRZEPŁYWU CIEPŁA I GĘSTOŚCI

STRUMIENIA CIEPŁA

W PŁASKIM PRZEWODZENIU CIEPŁA METODĄ ELEMENTÓW BRZEGOWYCH

Tomasz Janusz Teleszewski

^1a

1Katedra Ciepłownictwa, Politechnika Białostocka

at.teleszewski@pb.edu.pl

Streszczenie

W publikacji przedstawiono algorytm metody elementów brzegowych (MEB) wyznaczania linii przepływu ciepła w płaskim przewodzeniu ciepła metodą elementów brzegowych. Prezentowany algorytm stanowi alternatywę wo- bec najczęściej stosowanych obszarowych metod numerycznych, takich jak metoda różnic skończonych, metoda objętości skończonych, metoda elementów skończonych. W metodach tych stosowane są pracochłonne i skompli- kowane siatki wewnątrz obszaru, natomiast metoda elementów brzegowych wymaga jedynie dyskretyzacji brzegu.

Weryfikacja algorytmu została wykonana na podstawie znanych rozwiązań teoretycznych. W pracy przedstawiono przykład obliczeniowy przepływu ciepła przez dwa przewody centralnego ogrzewania obudowane wspólną izolacją cieplną. W celu wykonania symulacji komputerowych został napisany autorski program obliczeniowy.

Słowa kluczowe: przewodzenie ciepła, strumień ciepła, linie przepływu ciepła, metoda elementów brzegowych

NUMERICAL VISUALIZATION OF HEATLINE AND HEAT FLUX DENSITY IN TWO-DIMENSIONAL STEADY STATE THERMAL CONDUCTION USING BOUNDARY ELEMENT METHOD

Summary

The paper presents the numerical application of boundary element method (BEM) to calculate the heatline and heat flux density in two-dimensional steady state thermal conduction. The efficiency and the credibility of pro- posed algorithm were verified by numerical tests and the BEM solution is compared with theoretical results of steady state conduction in a plane wall with no heat generation and constant thermal conductivity. Numerical ex- ample is presented to illustrate steady state conduction in thermal isolation with a two pipe system of central heating. The BEM algorithm is an alternative to mesh method: Finite Difference Method, Finite Element Method or Finite Volume Method. The computer program was written in Fortran programming languages.

Keywords: heat conduction, heat flux, heatlines, boundary element method

1. WSTĘP

W wielu symulacjach cieplnych wykorzystuje się wizu- alizację linii przepływu ciepła oraz symulacje strumieni ciepła. Linie przepływu ciepła są ortogonalne do izoterm

(rys.1). W przypadku adiabatycznego brzegu linie przepływu ciepła są do niego równoległe, natomiast w przypadku brzegu o jednakowej temperaturze są do

niego prostopadłe. Z rozkładu gęstości strumienia ciepła w obszarze można odczytać intensywność przepływu ciepła. Wzrost wartości gęstości strumienia ciepła wska- zuje na wzrost intensywności przepływu ciepła podobnie jak zagęszczenie linii przepływu ciepła. Znajomość składowych strumienia ciepła również pozwala w sposób numeryczny wyznaczyć linie przepływu ciepła np. meto- dą Rungego-Kutty [1]. Linie przepływu ciepła mogą być wyznaczane w sposób graficzny przy zastosowaniu współczynników kształtu [2,3,4,5] jak i numeryczny.

Przykładowe rozwiązania linii przepływu ciepła w zagadnieniach cieplnych opartych na metodach siatkowych znajdują się w pracach: metodzie różnic skończonych [6], metodzie elementów skończonych [7,8,9,10], metodzie objętości skończonych [11] i metodzie wielosiatkowej [12]. Wadą metod obszarowych jest budowa pracochłonnych siatek. W prezentowanej meto- dzie elementów brzegowych dyskretyzacji wymaga jedynie brzeg obszaru.

Pole temperatury w płaskim ustalonym przepływie ciepła, w którym dominującym mechanizmem jest przewodzenie ciepła, jest opisane równaniem Laplace’a względem temperatury T [5]:

2 2

2 2 0

T T

x y

∂ +∂ =

∂ ∂ (1)

Składowe strumienia ciepła qx i qy zdefiniowane są zgodnie z prawem Fouriera [5]:

x

q T λ^∂x

= − ∂ (2)

y

q T λ^∂y

= − ∂ (3)

gdzie: λ jest współczynnikiem przewodzenia ciepła.

Wypadkowa gęstości strumienia ciepła jest wyznaczana z zależności:

2 2

x y

q= q +q (4)

W celu wyznaczenia linii przepływu ciepła Ξ wprowa- dzono zależność ortogonalną do izotach [13,14]:

T y λ x

∂Ξ= − ∂

∂ ∂ (5)

T x λ y

∂Ξ ∂

− = −

∂ ∂ (6)

Wobec czego równanie linii przepływu ciepła opisane jest zależnością:

x y

dΞ =q dy−q dx (7)

Wszystkie obliczenia zostały wykonane w autorskim programie obliczeniowym HEAT_CONDUCTION_2D opartym na metodzie elementów brzegowych napisanym w języku Fortran. Funkcje podcałkowe linii przepływu ciepła i gęstości strumienia ciepła zostały wyznaczone zgodnie z zależnościami (2-7).

Rys. 1. Linie przepływu ciepła i izotachy w przewodzeniu ciepła

2. BRZEGOWE RÓWNANIA

CAŁKOWE OPISUJĄCE PŁASKIE USTALONE PRZEWODZENIE CIEPŁA W OŚRODKU

JEDNORODNYM

Zagadnienie brzegowe dla równania Laplace’a (1) formułuje się w postaci przyjętego warunku brzegowego Dirichleta i Neumanna zakładającego znane wartości temperatury ( )T q% na części brzegu L^q (q∈L^q) i znane wartości strumienia ciepła q q%( ) na części brzegu L^T (q∈LT) (rys. 2) [15,16]:

( ) ( )

( ) ( )

( ) ( ) ( ) ( , ) ( ) ( , )

( ) ( , ) ( ) ( , ) ,

,

q T

T q

q T

L L

T q

L L

T q K dL T E dL

q K dL T E dL

L χ

− + + =

= − −

∈

∫ ∫

∫ ∫

p p q p q q p q

q p q q p q

p q

% % (8)

gdzie dla brzegu gładkiego χ( )p =1/2 oraz:

1 1

( , ) ln

2

K πλ r

 

=  

 pq

p q (9)

2

2 2

( ) ( )

( , ) 1 2

( ) ( )

x y

x x n y y n

E r

r x x y y

π

− + −

=

= − + −

p q p q

pq

pq p q p q

p q

(10)

Po wyznaczeniu niewiadomych T(p) na brzegu LT i q(p) na brzegu Lq, temperaturę w dowolnym punkcie (p∈A) rozpatrywanego obszaru (A) wyznacza się ze związku całkowego:

( ) ( )

( ) ( ) ( , ) ( ) ( , )

( ) ( ) , ( ) (A)

L L

T T E dL q K dL

L

= +

∈ ∈

∫ ∫

p q p q q p q

q p

(11) T

qr= −grad

1 2 3 4

T < < <T T T T1

T2

T3

T4

Ξ1

Ξ2

Ξ6 Ξ₅ Ξ4

Ξ3

Ξ7

Rys. 2. Szkic obrazujący zagadnienia brzegowe w obszarze płaskim

Wobec zależności (2-3) strumienie ciepła w kierunku x i y w przyjętym układzie współrzędnych w punktach (p) rozpatrywanego pola temperatury w obszarze (A) ograniczonym brzegiem (L) otrzymuje się, różniczkując funkcje podcałkowe (9-10) w wyrażeniu (11) odpowied- nio względem x i y:

( ) ( )

( ) ( )

( ) ( , ) ( , )

( ) ( ) ( )

( ) ( , ) ( , )

( ) ( ) ( )

x

L L

y

L L

T E K

q T dL q dL

x x x

T E K

q T dL q dL

y y y

λ λ

− − −

− − −

∂ ∂ ∂

= =

∂ ∂ ∂

∂ ∂ ∂

= =

∂ ∂ ∂

∫ ∫

∫ ∫

p p p

p p p

p p q p q

p q q

p p q p q

p q q

(12)

gdzie:

( )

( )

2

2 2

4

( , ) 1

2

( ) ( ) 2( )( )

( , )

2

x y

x x K

x r

y y x x n x x y y n

E

x r

π λ

π

∂ = − −

∂

− − − − − −

∂ =

∂

p q

p pq

p q p q p q p q

p pq

p q

p q (13)

( )

( )

2

2 2

4

( , ) 1

2

( ) ( ) 2( )( )

( , )

2

y x

y y K

y r

x x y y n x x y y n

E

x r

π λ

π

∂ = − −

∂

− − − − − −

∂ =

∂

p q

p pq

p q p q p q p q

p pq

p q

p q (14)

Znajomość składowych strumienia ciepła qx i qy pozwala wyznaczyć linie przepływu ciepła metodą Rungego- Kutty [1]. W pracy opracowano algorytm wyznaczania linii przepływu ciepła, przeliczając funkcje podcałkowe (13-14) zgodnie z definicją (5-6).

Po scałkowaniu funkcji podcałkowych (13-14) według definicji (5-6) i uwzględnieniu temperatury T oraz gęstości strumienia ciepła q na brzegu L wyznaczono zależność opisującą linie przepływu ciepła:

( ) ( )

( ) ( ) ( , ) ( ) ( , )

( ) ( ) , ( ) (A)

L L

T W dL q S dL C

L

Ξ = + +

∈ ∈

∫ ∫

p q p q q p q

q p

(15)

gdzie:

( , ) 1 2 lub

( , ) 1 2

y y

W arctg

x x

x x

W arctg

y y πλ

πλ

 − 

= −  − 

 − 

= +  − 

p q

p q

p q

p q

p q

p q

(16)

2

2 2

( ) ( )

( , ) 1 2

( ) ( )

y x

x x n y y n

S r

r x x y y

π

− − −

=

= − + −

p q p q

pq

pq p q p q

p q

(17)

gdzie: C jest to stała całkowania.

Na rys. 3 przedstawiono schemat blokowy autorskiego programu obliczeniowego HEAT_CONDUCTION_2D napisanego w języku fortran. Po zastąpieniu linii brze- gowej (L) układem linii cząstkowych, przy założeniu, że gęstości strumienia ciepła i temperatur na każdej linii są stałe, równania całkowe (8) można sprowadzić do ukła- du algebraicznych równań liniowych. Po wyznaczeniu niewiadomych temperatur i gęstości strumieni ciepła na brzegu L, pole temperatury, gęstość strumienia ciepła, linie przepływu ciepła wyznacza się całkując numerycz- nie funkcje podcałkowe całki (9-10, 13-14, 16-17).

Rys. 3. Schemat blokowy programu obliczeniowego HEAT_CONDUCTION_2D

3. WERYFIKACJA ALGORYTMU

W celu weryfikacji algorytmu przyjęto jednokierun- kowe przenikanie ciepła w płaskiej ścianie jednowar- stwowej o grubości 0.42m, zbudowanej z bloczków z betonu komórkowego o współczynniku przewodzenia ciepła 0.105W/(mK). Warunki brzegowe zagadnienia testowego zostały przedstawione na rys. 4.

( ) T q%

dL

q

r

_pq n_q

X Y

(A )

L

q

LT

nx

ny

( ) q q%

T q

L=L ∪L

p

Program

HEAT_CONDUCTION_2D (HC2D)

Czytanie danych: λ, dyskretyzacja brzegu (L), warunek brzegowy: , . T q%( ) q q%( )

Generowanie macierzy współczynników:

Generowanie wektora warunku brzegowego:

Rozwiązanie układu algebraicznych równań liniowych:

Wektor rozwiązań:

[ ]X =T( ), ( )q qq

( ) ( )

[ ] ( ) ( , ) ( ) ( , )

q T

q T

L L

A=

∫

q^qK^{p q}dL+

∫

T^qE^{p q}dL

( ) ( )

[ ] ( ) ( , ) ( ) ( , )

T q

T q

L L

B= −

∫

q^%qK^{p q}dL−

∫

T^%qE^{p q}dL [ ][ ] [ ]A X = B

Wyznaczenie pola temperatury T, składowych strumienia ciepła qx, qy, linii przepływu ciepła Ξ w obszarze A.

q=0 W / m2

Rys. 4. Warunki brzegowe w jednowarstwowej przegrodzie płaskiej

Pole temperatury w przekroju ściany wyznaczono ze wzoru analitycznego [3]:

( ) ( z w)

T w

T T x T x T

L

= + − (18)

gdzie: Tw, Tz są temperaturą na powierzchni ściany, L oznacza grubość ściany.

Gęstość strumienia ciepła została wyznaczona ze wzoru:

( ) ( z w)

T

T T q x dT

dx L

λ λ ⁻

= − = − (19)

Linie przepływu ciepła zostały wyznaczone poprzez scałkowanie zależności (18) zgodnie (5):

( ) (Tz Tw)

x y C

L

Ξ = − − + (20)

Do obliczeń przyjęto zerową linię prądu na wysokości y=0 (C=0). Błąd rozwiązania MEB gęstości strumienia ciepła i linii przepływu ciepła wyznaczono z następują- cych zależności:

MEB ^{TEO MEB}100%

TEO

q q

q q

δ = ⁻ (21)

MEB ^{TEO MEB}100%

TEO

δΞ =^Ξ Ξ^{− Ξ} (22)

gdzie indeksami TEO oznaczono wielkości teoretyczne, natomiast indeksem MEB wielkości wyznaczone metodą elementów brzegowych.

W tabeli 1 zestawiono błąd linii przepływu ciepła metody elementów brzegowych dla brzegu składającego się z 100 i 200 liniowych elementów w punktach przekro- ju przechodzących przez środek ścianki (y=0,21 m), natomiast w tabeli 2 przedstawiono błąd metody MEB dla gęstości strumienia ciepła. Maksymalny błąd wyzna- czania linii przepływu ciepła metodą MEB dla brzegu składającego się z 100 liniowych elementów nie przekra- cza 0.02%, natomiast dla brzegu zbudowanego z 200 elementów błąd MEB nie przekracza 0.01%. W przy-

padku gęstości strumienia ciepła maksymalny błąd MEB nie przekracza 0,04% dla brzegu zbudowanego z 100 elementów i 0,01% dla brzegu złożonego z 200 elemen- tów.

Tab. 1. Linie przepływu ciepła w przewodzeniu ciepła przez ściankę płaską- błąd rozwiązania MEB Współrzędne

węzłów

Rozwiązanie Rozwiązanie Błąd met.

teoretyczne num. MEB MEB

100 el.

xw yw ΞTEO ΞMEB ΞMEB

m m K K %

0,210 0,000 0,00000 0,00000 -

0,210 0,200 19,04762 19,04503 0,01360 0,210 0,400 38,09524 38,09126 0,01045 0,210 0,600 57,14286 57,13729 0,00975 0,210 0,800 76,19048 76,18291 0,00993 0,210 1,000 95,23810 95,23205 0,00635

200 el.

0,210 0,000 0,00000 0,00000 -

0,210 0,200 19,04762 19,04655 0,00561 0,210 0,400 38,09524 38,09311 0,00558 0,210 0,600 57,14286 57,13972 0,00549 0,210 0,800 76,19048 76,18649 0,00524 0,210 1,000 95,23810 95,23397 0,00433

Tab. 2. Gęstość strumienia ciepła w przewodzeniu ciepła przez przez ściankę płaską - błąd rozwiązania MEB Współrzędne

węzłów

Rozwiązanie Rozwiązanie Błąd met.

teoretyczne num. MEB MEB 100 el.

xw yw qTEO qMEB δqMEB

m m W/m² W/m² %

0,210 0,000 10,00000 10,00174 - 0,210 0,200 10,00000 9,99845 0,01553 0,210 0,400 10,00000 9,99626 0,03737 0,210 0,600 10,00000 10,00044 0,00444 0,210 0,800 10,00000 10,00197 0,01972 0,210 1,000 10,00000 10,00074 0,00744

200 el.

0,210 0,000 10,00000 9,99944 -

0,210 0,200 10,00000 9,99944 0,00559 0,210 0,400 10,00000 9,99945 0,00547 0,210 0,600 10,00000 9,99949 0,00507 0,210 0,800 10,00000 9,99964 0,00355 0,210 1,000 10,00000 10,00044 0,00441 Zagęszczenie podziału linii brzegowej konturu przekroju przewodu prostoliniowego powoduje zmniejszenie błędu metody elementów brzegowych. Graficzne rezultaty linii przepływu ciepła porównania MEB z rozwiązaniem teoretycznym (20) w wybranych punktach przekroju x=0,21m, zostały przedstawione na rys. 5.

Rys. 5. Porównanie rezultatów obliczeń MEB linii przepływu ciepła z rozwiązaniem teoretycznym (20) (x=0.21m)

4. PRZYKŁADY OBLICZENIOWE

y=0.00 X

o

Tz= −20 C

L=0.42 m

o

Tw= +20 C Y

( )

0.105 W / mK λ =

q=0 W / m2

Poniżej przedstawiono przykład obliczeniowy prze- wodzenia ciepła w izolacji wykonanej z wełny mineralnej (λ=0.042 W/(mK)) dwóch przewodów instalacji cen- tralnego ogrzewania. Do symulacji przyjęto rzeczywiste temperatury ścianek: TI=20 ^oC na izolacji, Tz=55 ^oC na ściance przewodu zasilającego oraz T^p=50 ^oC na ściance przewodu powrotnego. Różnica temperatur w instalacji centralnego ogrzewania wyniosła T∆ =5 ^oC. Wszystkie obliczenia wykonano metodą elementów brzegowych. Na rys. 6 przedstawiono warunki brzegowe i wymiary dla przykładu obliczeniowego: a=0.12m, b=0.08m, e=0.0187m, φ¹=φ²=0.0213m (DN15).

Na rys. 7 wykreślono pole temperatury w przekroju porzecznym wspólnej izolacji dwóch przewodów instala- cji centralnego ogrzewania.

Na rysunkach 8-9 przedstawiono składowe gęstości strumienia ciepła qx i qy, natomiast na rys. 10 zaprezen- towano wypadkową strumienia ciepła q w przekroju porzecznym wspólnej izolacji dwóch przewodów instala- cji centralnego ogrzewania.

Rys. 11 przedstawia wykreślone linie przepływu cie- pła wraz z zaznaczonym zwrotem przepływu ciepła w analizowanym przykładzie obliczeniowym.

Rys. 6. Warunki brzegowe w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania.

Rys. 7. Pole temperatury w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania (∆T=5 ^oC)

– rozwiązanie numeryczne MEB

Rys. 8. Składowa qx strumienia ciepła w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania

(∆T=5 ^oC) – rozwiązanie numeryczne MEB

Rys. 9. Składowa qy strumienia ciepła w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania

(∆T=5 ^oC) – rozwiązanie numeryczne MEB

Rys. 10. Wypadkowa gęstości strumienia ciepła w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania

(∆T=5 ^oC) – rozwiązanie numeryczne MEB e

a

T_Z T_P

T_I

Rys. 11. Linie przepływu ciepła w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania (∆T=5 ^oC) –

rozwiązanie numeryczne MEB

W drugim przykładzie obliczeniowym wykonano sy- mulację przepływu ciepła, w której założono większą temperaturę zasilania T^z=70 ^oC, w wyniku czego różnica temperatur między zasilaniem i powrotem, wyniosła

∆T=20 ^oC. Na rysunkach 12-14 wykreślono rezultaty obliczeń: pole temperatury, wypadkową gęstości stru- mienia ciepła i linie przepływu ciepła. W drugim przy- kładzie obliczeniowym wymiana ciepła odbywa się również miedzy przewodem zasilania i przewodem po- wrotnym instalacji centralnego ogrzewania.

Rys. 12. Pole temperatury w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania (∆T=20 ^oC)

– rozwiązanie numeryczne MEB

Rys. 13. Wypadkowa gęstości strumienia ciepła w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania

(∆T=20 ^oC) – rozwiązanie numeryczne MEB

Rys. 14. Linie przepływu ciepła w przekroju poprzecznym izolacji dwóch przewodów centralnego ogrzewania (∆T=20 ^oC)

– rozwiązanie numeryczne MEB

5. WNIOSKI

Metoda elementów brzegowych jest obecnie inten- sywnie rozwijana i trwają prace nad jej udoskonalaniem.

Przedstawiony algorytm metody elementów brzegowych wyznaczania linii przepływu ciepła i gęstości ciepła w przewodzeniu ciepła jest uzupełnieniem metody MEB i pozwala w sposób efektywny wyznaczać wymienione wielkości. Podstawową zaletą prezentowanego algorytmu jest eliminacja siatek wewnątrz obszaru, które są nieod- łącznym elementem klasycznych metod obszarowych takich jak metoda różnic skończonych, metoda elemen- tów skończonych czy metoda objętości skończonych.

Wyprowadzony w pracy algorytm również charakteryzu- je się małym błędem obliczeniowym.

W pracy przedstawiono przykład obliczeniowy symu- lacji przepływu ciepła w wspólnej izolacji dwóch prze- wodów instalacji centralnego ogrzewania, w którym wzrost różnicy temperatury między przewodem zasilają- cym i powrotnym powoduje wymianę ciepła między tymi przewodami.

Opracowanie zrealizowano w ramach pracy statutowej Politechniki Białostockiej.

Literatura

1. Press W.H., Flannery B.P., Teukolsky S.A., Vetterling W.T.: Numerical recipes in FORTRAN: example book The Art of Scientific Computing, Cambridge University Press, 1992.

2. Taler J., Duda P.: Solving direct and inverse heat conduction problems. Berlin: Springer Verlag, 2006.

3. Rao Y.V.: Heat transfer. India: Universities Press Lim., 2001.

4. Nag P.K.: Heat & mass transfer.2th ed. New Delhi: Tata McGraw-Hill Publ. Comp., 2007.

5. Holman J.: Heat transfer. 10^th ed. Tokyo: McGraw-Hill, 2009.

6. Hooman K., Gurgenci H., Dincer I.: Heatline and energy-flux-vector visualization of natural convection in a porous cavity occupied by a fluid with temperature-dependent viscosity. “Journal of Porous Media” 2010, Vol.

12, Iss. 3, p. 265 – 275.

7. Natarajan E., Basak T., Roy S.: Heatline visualization of natural convection flows within trapezoidal enclosures.

In: Proc. of the 5th IASME / WSEAS International Conference on Fluid Mechanics and Aerodynamics, August 25-27, 2007, Athens, Greece, World Scientific and Engineering Academy and Society, 2007, p.55 - 66.

8. Basak T., Chamkha A.J.: Heatline analysis on natural convection for nanofluids confined within square cavities with various thermal boundary conditions."International Journal of Heat and Mass Transfer" 2012, Vol.55, p.5526 - 5543.

9. Singh A.K., Roy S., Basak T.: A comprehensive Bejan’s heatline approach for natural convection heat transfer within inclined square cavities. ASME Heat Transfer Summer Conference, Rio Grande 2012, p. 1067 - 1076.

10. Kaluri R.S, Basak T.: Numerical visualization of heat flow and thermal mixing in various differentially heated square cavities using Bejan’s heatlines. ASME/JSME 8th Thermal Engineering Joint Conference Hawaii 2011, p.

T10051-T10051-10.

11. Speetjens M.F.M., Steenhoven A.A. van: Visualisation of heat transfer in unsteady laminar flows. “Computa- tional Thermal Sciences” 2011, Vol. 3(1), p.31 - 47.

12. Mahmud S., Fraser R.A.: Visualizing energy flows through energy streamlines and pathlines. "International Journal of Heat and Mass Transfer" 2007, Vol. 50(19-20), p.3990 - 4002.

13. Bejan A.: Convection heat transfer. 4^th ed. New Jersey : John Wiley 2013.

14. Costa V.A.F.: Bejan’s heatlines and masslines for convection visualization and analysis. “Applied Mechanics Reviews” 2006, Vol. 59, Iss. 3, p. 126 – 145.

15. Majchrzak E.: Metoda elementów brzegowych w przepływie ciepła. Częstochowa: Wyd. Pol. Częstochowskiej, 2001.

16. Brebbia C.A., Telles J.F.C., Wrobel L.C.: Boundary element techniques: theory and applications in engineering.

New York: Springer-Verlag, 1984.