Metoda elementów brzegowych z zastosowaniem wielokrotnej zasady wzajemności

(1)

Z E S Z Y T Y N A U K O W E

P O L I T E C H N I K I Ś L Ą S K I E J

ANDRZEJ J. N O W A K

METODA ELEMENTÓW BRZEGOW YCH Z ZASTO SO W AN IEM

WIELOKROTNEJ ZA SA D Y W ZAJEM NOŚCI

F Y K A

z. ne

G L I W I C I

1 9 9 3

(2)

P O L I T E C H N I K A Ś L Ą S K A

ZESZYTY NAUKOWE Nr 1202

ANDRZEJ J.

METODA ELEM ENTÓW BR ZEG O W YCH Z ZA STO SO W A N IEM

WIELOKROTNEJ ZASAD Y W ZAJEM NOŚCI

^ 9 ) 1 0 ) ^

N O W A K

G L I W I C E 1 9 9 3

(3)

O PIN IO D A W C Y

Prof. dr hab. inż. R y s z a r d P a rk iłn y Prof. zw. dr Inż. Jan S z a r g u t

K O L E G IU M R E D A K C Y J N E

R E D A K T O R N A C Z E L N Y — Prof. dr hab. inż. J a n B a n d ro w sk i R E D A K T O R D Z IA Ł U — D oc. dr hab. inż. Z b ig n ie w Rudnicki S E K R E T A R Z R E D A K C J I — M g r Elżbleła Leśko

R E D A K C J A M g r A n n a B łażK ie w icz

R E D A K C J A T E C H N IC Z N A Alicja N o w a c k a

W y d a n o za z g o d g Rektora Politechniki Ślgskie j

P L IS S N 0372-9796

W y d a w n ic tw o Politechniki Ślgskie j ul. K u ja w sk a 3, 44-100 G liw ice

N a k ła d 150+83 A r k . w y d . 8,5 A r k . d ru k . 7.375 P a p ie r o f f s e t , k l . I I I 70x100 80g O d d a n o d o d ru k u 1.03.93 P o d p is , d o d ru k u 1.13.93 D ru k u k o ń c z , w k w ie t n iu 1913

Z a m . 110|93 C en a z ł 11.900,—

Fotokopie, druk i o p ra w ę

w y k o n a n o w Z a k ła d z ie G raficzn ym Politechniki Ś lg sk ie j w G liw icach

(4)

Bez koleżeńskich dyskusji w In s t y t u c ie T e c h n ik i C ie p ln e j Politechniki Śląskiej w Gliwicach oraz w the W e s s e x I n s t it u t e o f T e c h n o lo g y University o f Portsmouth w Wielkiej Brytanii ta praca nigdy by nie powstała.

(5)

(6)

Spis treści

1. W s t ę p 17

1.1. Podstaw y m etody elementów b r z e g o w y c h ... 20

1.1.1. Równanie całkowe M E B ... 20

1.1.2. D yskretyzacja równania całkowego M E B ... 23

1.1.3. M acierze w pływ u . . . ... . . . . . . . . . 26

1.1.4. Uwzględnianie warunków b r z e g o w y c h ... 28

1.1.5. R ozw iązanie w punktach wewnętrznych ciała . . . 29

1.1.6. P ola ź r ó d ł o w e ... 29

1.2. M eto d y transform acji całek po obszarze w całki brzegowe . . 31

2. P o d s t a w y w ie lo k r o t n e j z a s a d y w z a je m n o ś c i 34 2.1. Zależności ogólne . . . ... .... . 34

2.2. D yskretyzacja brzegowego równania c a łk o w e g o ... 36

2.3. R ozw iązania podstawowe wyższych r z ę d ó w ... 39

2.4. Zagadnienia zbieżności szeregu W Z W ... 40

2.5. Im plem entacja komputerowa W Z W ... ... 43

3. Z a s to s o w a n ie W Z W d o m o d e lo w a n ia p ó l t e m p e r a t u r y 44 3.1. Ustalone przew odzenie ciepła ... 44

3.2. Nieustalone przew odzenie c i e p ł a ... 50

3.2.1. D yskretyzacja czasowa i p rz e s trz e n n a ... 54

3.2.2. K roczen ie w czasie . . ... 57

3.3. N ieliniow e zagadnienia, brzegowe przew odzenia ciepła . . . . 58

4. Z a s to s o w a n ie W Z W w z a g a d n ie n ia c h s p r ę ż y s t o ś c i i t e r m o - s p r ę ż y s t o ś c i 68 4.1. Podstawowe zależności teorii s p rę ż y s to ś c i... 68

4.1.1. Naprężenia te r m ic z n e ... 70

4.2. Równania całkowe sprężystości i terrnosprężystości . . . 70

4.2.1. D yskretyzacja równań całkowych sprężystości i ter- mosprężystości ... 73

4.3. Transform acja cajek po obszarze w całki b r z e g o w e ... 74

P o d s ta w o w e ozn aczen ia 11

5

(7)

4.3.ł. Zagadnk lia sprężystości z obciążeniem siłam i grawi

tacyjnym i i d o ś r o d k o w y m i... 76 4.3.2. Ustalone zagadnienia term osp rężystości... 79 4.3.3. Nieustalone zagadnienia te rm o s p rę ż y s to ś c i... 80

5. Z a s to s o w a n ie W Z W d o a n a liz y r ó w n a n ia H e lm h o lt z a 85 5.1. Zależności o g ó l n e ... 85 5.2. O bliczanie wartości własnych ... 88 6. Z a s to s o w a n ie W Z W w in ż y n ie r ii ją d r o w e j 91

6.1. Zagadnienia krytyczności reaktora jądrowego - przybliżenie dyfuzyjne je d n o g r u p o w e ... 91 6.1.1. Transform acja zagadnienia brzegowego w równanie

c a ł k o w e ... 94 6.1.2. Rozwiązania wyższych r z ę d ó w ... 95 6.1.3. Obliczanie parametru geom etrycznego reaktora . . . 95 7. Z a s to s o w a n ie W Z W w in n y c h d z ie d z in a c h m e c h a n ik i 97 7.1. Zagadnienia opływ u c i a ł ... 97 7.2. Drgania harmoniczne w cienkich płytach sprężystych . . . . 98

8. P o d s u m o w a n ie i w n io s k i k o ń c o w e 100

B ib lio g r a fia 102

A R o z w ią z a n ia p o d s t a w o w e w y ż s z y c h r z ę d ó w d la z a g a d n ie ń

p o t e n c ja ln y c h 111

B S to s u n e k w s p ó łc z y n n ik ó w lic z b o w y c h w r o z w ią z a n iu p o d

s t a w o w y m 114

Streszczenie 115

(8)

e q u a t io n ... . . ; 94 6.1.2. Higher order fundamental s o lu t io n s ... 95 6.1.3. Calculation o f reactor geom etric b u c k lin g ... 95 7. A p p lic a t io n o f M R M in o t h e r b ra n c h e s o f m e c h a n ic s 97 7.1. Fluid flow p ro b le m s ... 97 7.2. Harm onic vibration o f thin elastic p l a t e s ... 98

8. F in a l r e m a r k s an d c o n c lu s io n s 100

B ib lio g r a p h y 102

A H ig h e r o r d e r fu n d a m e n ta l s o lu tio n s fo r p o t e n t ia l p r o

b le m s 111

B C o e ffic ie n t r a t io in fu n d a m e n ta l s o lu tio n 114

S u m m a r y 125

(10)

Inhaltsverzeichnis

W i c h t i e F o r m e lz e ic h e n 11

1. E in fü h r u n g 17

1.1. Grundlagen der R an delem en tm eth ode... 20

1.1.1. Integralgleichung der R E M ... 20

1.1.2. Diskretisierung der R E M -In teg ra lg leich u n g ... 23

1.1.3. E in flu ß m a trizen ... 26

1.1.4. Berücksichtigung der R an d b ed in gu n gen ... 28

1.1.5. Lösung für In n e n p u n k te ... 29

1.1.6. Potentialfelder m it internen m it Q u e l l e n ... 29

1.2. M ethoden der Umwandlung von Gebietsintegrale in Randin tegrale ... ... . ... 31

2. G r u n d la g e n fü r d ie M e t h o d e d e r m e h r fa c h e n R e z i p r o z i t ä t 34 2.1. A llgem eine G ru n d la g e n ... 34

2.2. Diskretisierung der Randintegralgleichung ... 36

2.3. Fundamentallösungen höherer O r d n u n g ... 39

2.4. Problem e der Konvergenz von M M R - R e ih e n ... 40

2.5. Computerrealisierung der M M R ... 43

3. A n w e n d u n g d e r M M R z u r M o d e llie r u n g v o n T e m p e r a t u r fe ld e r n 44 3.1. Stationäre W ä r m e le itu n g ... 44

3.2. Instationäre W ä r m e le it in g ... 50

3.2.1. Zeit- und Raum diskretisierung... 54

3.2.2. Zeitschritt-Prozedur ... 57

3.3. Nichtlineare Randwertproblem e der W ä rm e le itu n g ... 58

4. A n w e n d u n g d e r M M R in d e r E la s t iz it ä t u n d T h e r m o e la s t i- z i t ä t 68 4.1. Grundlagen der E la s t iz it ä t ... 68

4.1.1. Tem peraturspannungen... 70

4.2. Integralgleichungen der Elastizität und Th erm oelastizität . . 70

4.2.1. Diskretisierung der Integralgleichungen der Elastizität und T h e r m o e la s t iz it ä t... 73

9

(11)

4.3. Umwandlung der Gebietsintcgrale in R a n d in t e g r a le ... 74 4.3.1. Aufgaben der Therm oelastizität unter Berücksichti

gung von G ravitations- und Zentripetalkraft Belas

tungen ... 76 4.3.2. Stationäre Problem e der T h e r m o e la s tiz itä t... 79 4.3.3. Instationäre Problem e der T h e r m o e la s tiz itä t... 80

5. A n w e n d u n g d e r M M R z u r A n a ly s e d e r H e lm h o lt z - G le i-

ch u n g 85

5.1. Allgem ein e G ru n d la g e n ... 8-5 5.2. Berechnung von E igen w erten ... 88

6. A n w e n d u n g d e r M M R in d e r K e r n e n e r g ie t e c h n ik 91 6.1. Anwendung der Diiïusionshâhcrung und des Einsgruppen-

energievchrfahrens zur Untersuchung der Kritikalit.ät des A t o m r e a k t o r s ... 91 6.1.1. Umwandlung des Randwertproblem s in eine Integral

gleichung ... 94 6.1.2. Lösungen höherer Ordnung ... 95 6.1.3. Berechnung des geometrischen Parameters des Reaktors 95 7. A n w e n d u n g d e r M M R in a n d e r e n G e b ie t e n d e r M e c h a n ik 97 7.1. U m ström u n gsp rob lem e... 97 7.2. Harmonische Schwingungen von dünnen elastischen St äben . 98

8. Z u s a m m e n fa s s u n g u n d S c h lu ß fo lg e ru n g e n 100

L it e r a t u r 102

A F u n d a m e n ta llö s u n g e n h ö h e r e r O r d n u n g d e r P o t e n t ia l

p r o b le m e 111

B V e r h ä lt n is d e r Z a h le n k o e ffiz ie n te n in d e r F u n d a in e n -

t a llö s u n g 114

K u r z fa s s u n g 125

(12)

Podstawowe oznaczenia

Skalarne wielkości fizyczne oznaczone są w pracy kursywą. W y tłu sz

czenia stosowano dla m acierzy powstałych w wyniku dyskretyzacji równań całkowych. N ad sym bolam i wielkości wektorowych umieszczono strzałki.

W szystkie występujące w pracy wielkości fizyczne w yrażone są w je d nostkach m iędzynarodowego układu SI lub sprowadzone zostały do postaci bezwymiarowej.

W części pracy dotyczącej sprężystości i term osprężystości stosowano konwencję sumacyjną Einsteina, w której dwukrotnie występujący indeks oznacza sumowanie. Dla zagadnień dwuwym iarowych, w których indeksy 1 i 2 odpowiadają kolejno osiom x i y. mamy zatem

G m a m — ^ j “ł" ^2 ®mrn ^11 “ł~ ®22

W stosowanej notacji różniczkowanie pola a według współrzędnej geo

metrycznej x m oznacza się przecinkiem D a U ,m —

0 x 1H

Sym bole łacińskie

ii współczynnik wyrównywania tem peratury

A macierz główna układu równań liniowych M E B . (1.29) ,lj współczynnik w rozwiązaniu podstawowym , (2.25), (2.26) /i*01 uogólniony człon źródłow y w równaniu różniczkowym Poisso-

na

//•'* ciąg laplasjanów funkcji źródła definiowany przez równanie (2-13)

/)*/* - średnia, wartość funkcji w rozważanym obszarze

b', macierz kolunuiowa zawierająca wartości funkcji />*■'* w punk

tach węzłowych należących do n-tego elementu brzegowego.

(2.19)

/>,„ składowa siły masowej w kierunku ni

\>'u jednostkowe obci;|żenie punktowe działające w kierunku ⁱⁱⁱ

(13)

ómL -■ maksymalna wartość funkcji 6(j) w rozważanym obszarze

Bg

- param etr geom etryczny reaktora, (6.4)

B j - współczynnik w rozwiązaniu podstawowym , (2.25), (2.27) B j - macierz kolumnowa zawierająca wartości funkcji 6,jl w punk

tach węzłowych

Bm

- param etr m ateriałow y reaktora, (6.2) c - pojemność cieplna właściwa (ciep ło właściwe)

c, - współczynnik zależny od kąta wewnętrznego brzegu f w punk

cie ^i. Dla brzegu gładkiego ^c, = 0.5 D współczynnik dyfuzji w równaniu (6.1) D j - całka po obszarze; (j

=

0,1,2,...)

cm składowa wektora jednostkowego w kierunku m E - moduł sprężystości Younga

f - macierz kolumnowa wyrazów wolnych w liniow ym układzie równań M E B , (1.29)

f j funkcja interpolacyjna w Podwójnej Zasadzie W zajem ności, (1.37)

y przyśpieszenie graw itacyjne

g 'n - macierz wierszowa zawierająca te współczynniki macierzy wpływ u G j, które odpow iadają całkowaniu z punktu obser

wacji iw zdłuż elementu brzegowego l 'n G j macierz w pływ u M E B dla strumienia ciepła h współczynnik wnikania ciepła

łij" - macierz wierszowa zawierająca te współczynniki macierzy w pływ u

H; ,

które odpow iadają całkowaniu z punktu obser

wacji i wzdłuż elementu brzegowego ]'„

H j - macierz wpływu M E B dla tem peratury

składowa normalna gęstości prądu neutronów, (6.14) macierz kolumnowa gęstości prądu neutronów, (6.15) współczynnik przewodzenia ciepła

efektyw ny współczynnik mnożenia, (6.1)

funkcja McDonalda. (zm odyfikowana funkcja Bcsscla) rzędu j normalna zewnętrzna do brzegu 1’

M

liczba wyrazów obciętego szeregu (2.31) ,V liczba elem entów brzegowych

prn składowa, (w kierunku m ) siły powierzchniowej

!>"m rozwiązanie podstawowe teorii sprężystości, (4 .1 1)

Plmk ~ rozwiązanie stowarzyszone z rozwiązaniem podstawowym te

orii sprężystości, (1.16)

P macierz sił powierzchniowych, (4.21)

r odległość dwóch punktów (punktu obserwacji i punktu całko

wania)

R - macierz reprezentująca wartości całki obszarowej w sformuło

waniach (2.18) i (4.21)

ii

12

(14)

Rm - reszta szeregu, (2.35)

q - składowa normalna gęstości strumienia ciepła, (1.4)

q‘ - analog strumienia ciepła dla rozwiązania podstawowego, (1.5) - pom ocnicza wielkość definiowana równaniami (2.12) i (2.28) q n - macierz kolumnowa zawierająca wartości strumienia ciepła w punktach węzłowych należących do n-tego elementu brzego

wego, (1.22)

qv - wydajność wewnętrznych źródeł ciepła, (1.6)

Q - macierz kolumnowa zawierająca wartości strumienia ciepła we wszystkich punktach węzłowych, (1.28)

t

- czas

T

- tem peratura

T - macierz kolumnowa zwierająca wartości tem peratury we wszystkich punktach węzłowych, (1.28)

T

- rozw iązanie szczególne równania Poissona, (1.34)

T0 - warunek początkow y dla zagadnień przew odnictw a cieplnego, (3-2)

T * - rozwiązanie podstawowe spełniające równanie (1.3)

T*ti) _ rozw iązanie podstawowe j- te g o rzędu spełniające równania (2.11) i (2.24)

T n - macierz kolumnowa zawierająca wartości tem peratury w punktach węzłowych należących do n-tego elementu brzego

wego, (1.21)

um - odkształcenie w kierunku m

M*m ~ rozwiązanie podstawowe teorii sprężystości, (4.10)

- rozwiązanie podstawowe pierwszego rzędu teorii sprężystości, (4.28)

U - macierz przemieszczeń, (4.21)

ttfk) - ciąg funkcji definiowanych równaniem (2.14)

- macierz kolumnowa zawierająca wartości funkcji w punk

tacji węzłowych należących do n-tego elementu brzegowego, (2 .20)

W , - m acierz kolumnowa zawierająca wartości funkcji w punk

tach węzłowych x ,y - współrzędne glpbalne

X - macierz kolumnowa niewiadomych w układzie równań linio

wych M E B , (1.29) Yt - punkt bieżący, rys. 1.1

Yz

- punkt działania, źródła, rys. 1.1

13

(15)

Sym bole greckie

a - współczynnik rozszerzalności termicznej liniowej, (4.6) otj - współczynnik liczbowy w aproksymacji stosowanej w P o

dwójnej Zasadzie W zajem ności, (1.37) 7 - współczynnik w równaniu Helm holtza (5.1) T - brzeg obszaru

fi

<5(yb,

Y.)

- funkcja Diraca działająca w punkcie

Yz

6im funkcja Kroneckera

At

- krok czasu

ćtm składowe tensora odkształcenia, (4.2)

0 współczynnik definiujący typ schematu różnicowego fi - moduł odkształcenia postaciowego

u współczynnik Poissona w równaniach sprężystości lub śre

dnia liczba neutronów rozszczepieniowych na jeden akt rozszczepienia

£, i] - współrzędne lokalne

o - gęstość

<r,m - składowe tensora naprężeń

<jJm składowe tensora naprężeń termicznych

a"m - składowe tensora stowarzyszonego z rozwiązaniem podsta

w ow ym teorii sprężystości, (4.20) - tensor < m i-te g o rzędu, (4.34)

makroskopowy przekrój czynny na absorpcję neutronów,

(

6

.

1

)

E/ makroskopowy przekrój czynny na rozpraszanie neutro

nów, (6.1)

<j> strumień neutronów, (6.1)

# 1 macierz wierszowa zawierająca lokalne funkcje kształtu M E B , (1.21) i (1.22) lub macierz strumienia neutronów, (6.15)

wm składowa (w kierunku m ) prędkości kątowej

fi obszar w przestrzeni Ił2 lub /?' ograniczony brzegiem V

Indeksy dolne

u oznacza numer elementu brzegowego i oznacza punkt obserwacji

oznacza różniczkowanie według współrzędnej geom etrycznej j

7 odnosi się do warunku brzegowego Ncumamia / odnosi się do warunku brzegowego Dirichleta w odnosi się do punktów wewnątrz obszaru

(16)

Indeksy górne

* oznacza rozwiązanie podstawowe pochodna względem czasu druga pochodna w zględem czasu oznacza rząd rozwiązania podstawowego oznacza numer kroku czasu

7 - oznacza macierz transponowaną oznacza, rozwiązanie szczególne

Inne sym bole

V V 2

m On

gradient

operator Laplace'a (laplasjan)

różniczkowanie wzdłuż zewnętrznej normalnej do brzegu

(17)

.

'

' I;,, i ' • * : • :

' ' ■• ijk

■

V

(18)

Rozdział 1 W stę p

llosnące wym agania odnośnie do jakości w ytw orów , ich trwałości, niezawo

dności działania itp. stawiają przed konstruktorami i producentami coraz poważniejsze zadania. Począwszy od etapu projektowania, poprzez prawie wszystkie stadia produkcji przeprowadza się dzisiaj powszechnie analizę większości cech konstrukcyjnych i eksploatacyjnych wytworu. Podejściu ta

kiemu sprzyja burzliw y rozwój techniki komputerowej i idący z nim w pa

rze im ponujący postęp w dziedzinie nowoczesnych technik obliczeniowych.

Jesteśmy świadkami powstawania nowych, rozległych dziedzin nauki pole

gających na zastosowaniu komputerowych m etod obliczeniowych do proje

ktowania, konstruowania i sterowania. Term iny takie jak C om p u ter Aided Design (K om puterow e Wspomaganie Projektow ania) czy też C om pu ter A i

ded Engineering (Kom puterow e Wspomaganie Obliczeń Inżynierskich) na trwałe weszły do języka naukowego.

Coraz istotniejszy staje się problem optym alizacji konstrukcji i procesów.

Rosnące wciąż ceny podstawowych nośników energii, a co za tym idzie - również je j bardziej przetworzonych form, w ym agają rozwiązań tanich, a j e dnocześnie spełniających określone warunki. Ze względów finansowych dąży się do ograniczenia liczby budowanych prototypów i instalacji doświadczał nycli. Zdobywanie informacji o obiekcie czy nowej konstrukcji coraz częściej zastępuje się modelowaniem m atem atycznym .

O bliczenia inżynierskie z dziedziny w ytrzym ałości m ateriałów, przepły

wu ciepła lub hydromechaniki prowadzi się dzisiaj prawie w yłącznie m eto

dami numerycznymi. M etod y analityczne w ym agają bowiem zwykle, l rud

nych w obecnej dobie do zaakceptowania, uproszczeń geom etrii analizowa

nego obiektu i innych założeń upraszczających. Ponadto, rozwiązania uzys

kiwane m etodam i analitycznym i zawierają zw yk le szeregi lub całki, których obliczenie również wymaga użycia komputera i m etod analizy numerycznej.

Spośród m etod numerycznych, które znalazły powszechne zastosowanie w obliczeniach inżynierskich,należy wymienić (w kolejności powstawania):

(19)

18 /. W stęp

• m etodę różnic skończonych (M U S )

• m etodę elementów skończonych (M E S )

• m etodę elementów brzegowych (M E R )

Każda z tych metod może być uważana za szczególną wersję ogólniej

szej m etody ważonej residualnej (M W R ) [18]. Różnice pom iędzy metodam i wynikają z przyjęcia różnych funkcji wagowych oraz wykonania całkowania przez części różną liczbę razy.

Popularność i skuteczność metod numerycznych wynika przede wszyst

kim z ich następujących cech:

1. M etod y numeryczne pozwalają modelować zagadnienia brzegowe w obszarach o praktycznie dowolnych kształtach. Trzeba oczyw iście pa

miętać o tym . że poszczególne m etody stosują różne reprezentacjo nu

meryczne kształtu. T ym samym, różny jest w poszczególnych m eto

dach stopień trudności przy tworzeniu modelu geom etrycznego roz

ważanego obiektu.

2. M etody numeryczne pozwalają uwzględniać praktycznie dowolne nie

liniowości zagadnień brzegowych. Mogą to być zarówno nieliniowości typu zależności właściwości materiałowych od poszukiwanego rozwią

zania. jak i nieliniowości warunków brzegowych, źródeł, czy nawet nieliniowości geometrycznych. O czywiście stopień komplikacji m ate

matycznych wynikających z uwzględniania nieliniowości jest różny w poszczególnych m etodach.

3. Program y komputerowe, stanowiące implementacje algorytm ów po

szczególnych metod, są programami uniwersalnymi. Oznacza to. że za pomocą jednego programu można analizować różne obiekty z różnymi warunkami brzegowymi i początkowym i. P rzy uwzględnieniu ponadto podobieństwa zjawisk wyrażającego się podobieństwem zagadnień brzegowo-początkowych ten sam kod komputerowy może służyć do analizy różnych zjawisk fizycznych.

4. W spółczesne pakiety M RS. M ES i M E R współpracują zw ykle z grafi

cznymi p i f - oraz postprocisortuni. T y m samym liczba danych wstęp

nie przygotowywanych przez użytkownika programu jest niewielka.

Trzeba również dodać, że /m- i poslpm n sory nie tylko autom atyzują proces przygotowania danych, ale również umożliwiają graficzną we

ryfikację zbioru danych oraz wizualizację rezultatów obliczeń.

Z punktu widzenia użytkownika programu istotne jest to. że m etody różnic skończonych oraz elementów skończonych są tzw . m ilo ilu w i obszaro

wymi. O znacza to, że ich stosowanie wym aga dyskrctyzacji całego obszaru.

(20)

1. W stęp 19

M etoda elem entów brzegowych sprowadza rozwiązanie problemu do nume

rycznego rozwiązania równania całkowego. W niektórych sytuacjach prak

tycznych występujące w tym równaniu całki są całkami jed yn ie po brzegu ciała. T y m samym, rozwiązując równanie całkowe, należy poddać dyskretyzacji tylko brzeg obszaru. Rozw iązanie w wybranych punktach wewnętrz

nych uzyskuje się w następnym etapie (tzw . postprocessingu) na podstawie znajomości wartości brzegowych. W praktyce oznacza to istotne uproszcze

nia na etapie przygotowania danych dla program ów komputerowych. Dla zagadnień dwuwym iarowych użytkownik zajmuje się jed yn ie obiektem je d nowym iarowym , jakim jest linia brzegowa. W przypadku zagadnień trójw y

miarowych dyskretyzuje się powierzchnię brzegową - obiekt dwuwym iarowy.

W tym sensie mówi się często o obniżeniu wymiarowości zagadnienia brze

gowego poprzez zastosowanie m etody elementów brzegowych.

N ależy jednak zaznaczyć, że stosując m etodę elementów brzegowych (M E B ) do m odelowania praktycznych problemów inżynierskich, nie zawsze uzyskuje się sformułowania zawierające tylko całki brzegowe. Częstokroć bezpośrednie zastosowanie M E B prowadzi do równań całkowych zawiera

jących również całki po wnętrzu obszaru. W konsekwencji, w istotny sposób komplikuje się etap rozwiązania równania całkowego. Dyskretyzacji należy bowiem poddać zarówno brzeg, jak i wnętrze obszaru. Nic więc dziwnego, że od wielu lat prowadzone są prace nad uzyskaniem takich sformułowań M E B , które pozwalają uniknąć podziału obszaru na elementy m im o wystąpienia całek po jego wnętrzu. W y m a g a to oczywiście zastąpienia danych o wnę

trzu obszaru dodatkowym i informacjami o zachowaniu się rozwiązania na jego brzegu. Zagadnienia omawiane w niniejszej pracy dotyczą tzw. w ie lo k r o t n e j z a s a d y w z a je m n o ś c i W Z W (ang. M ultiple Reciprocity M ethod).

M etoda pozwala właśnie bez dyskretyzacji w nętrza obszaru rozwiązywać za pomocą M E B te zagadnienia inżynierskie, które normalnie takiej dyskrety

zacji wym agają.

M etod a W Z W jest oryginalną m etodą opracowaną przez autora roz

prawy pod koniec lat 80. Polega ona na przekształceniu występujących w równaniach M E B całek po obszarze w szereg całek po brzegu tego obszaru.

W procesie transformacji wykonuje się całkowanie przez części nieskończoną liczbę razy otrzym ując w granicy całkowicie brzegowe sformułowanie me

tody elem entów brzegowych.

W dalszej części tego rozdziału zaprezentowano ogólny opis m etody elementów brzegowych stanowiący jednocześnie podstawę przeglądu lite

ratury dla problemów omawianych w pracy. Idea m etody została zaprezen

towana na przykładzie zagadnienia przewodzenia ciepła, przy czym podane zależności m ogą być również wykorzystane w sposób bezpośredni do m ode

lowania innych pól potencjalnych.

O gólny opis W Z W jest treścią rozdziału drugiego, podczas gdy następne rozdziały pracy przedstawiają zastosowanie tej techniki do analizy wybra

nych zagadnień brzegowych. R ozdziały 1 oraz -i. w których opisano zasto

(21)

20 1. W stęp

sowanie W Z W do modelowania pól temperatury, są indywidualnym i osiąg

nięciami autora niniejszej rozprawy. R ezu ltaty zaprezentowane w rozdziale 4 są w większości plonem wieloletniej współpracy z Wessex Institute o f T e

chnology w W ielk iej Brytanii. W Z W znalazła tam grono kontynuatorów, k tórzy wykorzystując doświadczenie autora, zastosowali tę technikę do roz

w iązyw ania zagadnień sprężystości. A n aliza stacjonarnych i niestacjonar

nych problem ów termosprężystości zaowocowała wspólnym i publikacjami w m iędzynarodowych czasopismach.

A u tor pracy jako pierwszy zastosował W Z W do rozw iązyw ania równania Helm holtza. Idea ta pozw oliła innym badaczom na opracowanie efektyw nych m etod obliczania wartości własnych zagadnień brzegowych. Problem y te przedstawiono w rozdziale 5 pracy. M ateriał stanowiący treść rozdziału 6 jest bezpośrednim efektem wym iany doświadczeń naukowych z Japan A to m ie Research Institute. R ozdział 7 opracowany został na podstaw ie lite

ratury, a jego celem było jed yn ie przedstawienie potencjalnych możliwości W Z W .

G łów ny nacisk położono w pracy na wyprowadzenie brzegowych sfor

mułowań M E B oraz na aspekty numeryczne wielokrotnej zasady w zajem ności ograniczając do niezbędnego m inimum zagadnienia takie jak badanie zbieżności szeregu, dowody jednoznaczności rozwiązań itp. Praktyczne za

stosowania W 'Z W zaprezentowano rozwiązując szereg zagadnień brzegowych mechaniki.

1.1. Podstaw y m etody elementów brzegowych

M ateriał zawarty w tym podrozdziale nie jest system atycznym wykładem m etody elementów brzegowych i zawiera jed yn ie te informacje, które są niezbędne do zrozum ienia istoty m etody i jej zastosowania do analizy za

gadnień potencjalnych. W szczególności przedm iotem rozważań są ustalone procesy przewodzenia ciepła. Pełny wykład M E B znaleźć można w licznych monografiach zagranicznych, np. [3], [9], [12]. W Polsce, jak do tej pory, nie ukazała się żadna monografia poświęcona w yłącznie M E B , jednakże ob

szerne inform acje na tem at zastosowania m etody do modelowania procesów cieplnych znaleźć można w artykule [5] oraz w monografiach [19] i. [75]. M E B poświęcone też były prace doktorskie i habilitacyjne, np. [14]. [37], [39].

1.1.1. Równanie całkowe M E B

Jak wspomniano wcześniej, M E B jest szczególną wersją m etody ważonej residualnej (M W R ), [18]. Równania M E B m ożna zatem w yprowadzić z rów

(22)

1.1. Podstaw y m etod y elem entów brzegowych 21

nań M W R , np. [9] i [12], przyjm ując jako funkcje wagowe tzw . rozw iąza nie podstawowe oraz jeg o pochodną wzdłuż normalnej zewnętrznej do brzegu ciała. W zagadnieniach cieplnych rozwiązanie podstawowe jest polem tem peratury w nieskończonym obszarze (całej przestrzeni) w yw ołanym działa

niem punktowego źródła ciepła.

A ltern atyw n ie równania M E B m ożna uzyskać poprzez w ykorzystanie tzw. zasady w za je m n oś ci [32], [38], czyli sym etrycznej tożsamości Greena.

Zasada ta, powszechnie stosowana np. przy analizie zagadnień sprężystości, nie zdobyła większej popularności w przepływ ie ciepła, praw dopodobnie z braku prostej interpretacji fizycznej w tej gałęzi nauki. Z m atem atycznego punktu widzenia zasada wzajemności jest form ułą dwukrotnego całkowania przez części i stanowi związek pom iędzy całką po obszarze i całką po jego brzegu. Zależnie od tego czy rozważany obszar jest dwu- czy trójw ym iarow y, całka po obszarze jest całką powierzchniową bądź całką objętościową. D la

tego też, celem zachowania ogólności zapisu, element różniczkowy obszaru oznaczym y sym bolem d ii, zaś element różniczkowy brzegu obszaru sym bo

lem dT.

Zasada wzajemności jest spełniona przez dwa dowolne pola, w szczegól

ności przez dwa pola tem peratury:

T

oraz

T"

J (TV2r - TV 2 T) dfi = J { t ~ -*

^dr ^(1-1)

gdzie: fi - rozważany obszar, T - brzeg obszaru fi,

^ - różniczkowanie wzdłuż normalnej zew nętrznej do brzegu r.

Poniew aż pola

T

i T * m ogą być w ybrane zupełnie dowolnie, najdogo

dniej jest utożsam ić pole

T

z poszukiwanym polem tem peratury, za pole zaś

T"

przyjąć rozwiązanie podstawowe. Rozw iązanie podstawowe spełnia następujące równanie różniczkowe [27], [63]

V2T * = 6 (Y b,Y z) (1.2) gdzie: Yz - punkt działania źródła, rys. 1.1,

Yb - punkt bieżący,

6 - dystrybucja Diraca (delta-funkcja).

W ielkość T * zależy od dwóch param etrów: współrzędnych punktu działa

nia źród ła Yz i współrzędnych punktu bieżącego % , w którym obserwujem y efekty działania źródła. Przem ieszczając początek układu współrzędnych do punktu Yz łatw o stwierdzić, że rozwiązanie podstawowe T * jest ostatecznie funkcją odległości r punktów Yz i Yb

ln r problem y dw uw ym iarow e T * =

2 tt

(1.3) 1 1 problem y trójw ym iarow e

47r r

(23)

22

1. W stęp

Różniczkując pola T i T " wzdłuż normalnej zewnętrznej do brzegu ciała, otrzym ujem y gęstości strumienia ciepła q i q*

= - k d T dn d T *

dn

(1.4)

(1.5)

gdzie: k - współczynnik przewodnictwa cieplnego.

Załóżm y chwilowo, że pole T jest ustalonym polem tem peratury i tym samym, że opisuje go równanie Poissona o postaci

k V 2T + qv = 0 w fi

gdzie: qv - wydajność wewnętrznego źródła ciepła.

(1.6)

Rys. 1.1. Globalny (x ,y ) i lokalny (£, Tj) układy współrzęd

nych dla fu n k cji T ' Fig. 1.1. Global (x ,y ) and local

( ( , r ) ) co-ordinate sys

tems f o r fu n ction T *

W ystępujące w równaniu (1.6) wewnętrzne źródła ciepła mogą być albo rzeczywistym i źródłam i ciepła działającym i w rozpatry

wanym obszarze fi, albo źródłam i fikcyjnym i. Źródła fikcyjne two

rzone są przez zgrupowanie w funkcji qv wszystkich członów o- pisujących efekty niestacjonarne, nieliniowe itp. Z punktu widzenia M E B źródła fikcyjne traktowane są dokładnie tak samo jak źródła rzeczywiste.

Podstawiając w miejsce laplasja- nów w równaniu (1.1) odpow ie

dnie wyrażenia z równań (1.2) i (1.6), otrzym ujem y

i L

T * (Y z,Y b) - q v(Y b) + T ( Y b) 6 (Y z,Y b)

d n (Y b) d n (Y b) J

d n (Y b) =

d T (Y b) (1.7)

Różniczkowanie i całkowanie w równaniu (1.7) przebiega według współ

rzędnych punktu bieżącego Yb, co zaznaczono symbolicznie, wprowadzając ten punkt jako argument operacji różniczkowania i całkowania.

(24)

1.1. Podstawy m etod y elementów brzegowych 23

W ykorzystując następnie własność filtracji funkcji delta Diraca i doko

nując prostych operacji podstawienia, otrzym ujem y ostatecznie

k c (Y z) T ( Y z) + j q ' ( Y z,Y b) T ( Y b)d T (Y b) =

= / T ' ( Y z,Y b) q ( Y b) d r ( Y b) - [ r ( Y z,Y b) q v(Y b) d n ( Y b) (1.8)

J r Jn

gdzie: c (Y z) - współczynnik zależny od położenia punktu działania źródła Yz.

Jeśli punkt Yz leży poza rozważanym obszarem fi, współczynnik c (Y z) = 0. Dla punktów leżących wewnątrz obszaru f i współczynnik c (Y z) = 1, w przypadku zaś kiedy punkt Yz należy do brzegu T, współczynnik ten jest równy kątowi wewnętrznemu brzegu w punkcie Yz [9], [12], [27]. Dla brzegu gładkiego otrzym u je się więc c (Y z) = 0.5;

W dalszej części pracy, celem uproszczenia zapisu, tem peraturę w punk

cie Yz oznaczono przez Ti, współczynnik zaś c (Y z) przez c,-. Jednocześnie pominięto argumenty zarówno funkcji podcałkowych, jak i samych operacji całkowania i różniczkowania. Równanie (1.8) można więc zapisać w postaci

k c i T i + [ q * T d Y =

f

T " q d T -

f

T * q v dft (1.9)

Jr Jr Jn

Równanie (1.9) jest równaniem całkowym ze względu na poszukiwane pole tem peratury T oraz gęstość strumienia ciepła q i stanowi całkową re

prezentację rozwiązania wyjściowego zagadnienia brzegowego (1.6).

1.1.2. Dyskretyzacja równania całkowego M E B

Analityczne rozwiązanie równania (1.9) m ożna uzyskać jed yn ie w obsza

rach o bardzo prostych kształtach. W przypadku ogólnym poddaje się je dyskretyzacji i rozw iązuje numerycznie.

Przystępując do dyskretyzacji równania (1.9), zauważmy, że dla bezźró- dłowych pól tem peratury równanie to upraszcza się do postaci zawierającej tylko całki po brzegu ciała

k c { T i +

J

^q*^T^{d T}⁼

J

^{T ’}^{q d T}

(1-10)

P om inięty w równaniu (1 1 0 ) wyraz poddano analizie w dalszej części pracy.

N a wstępie dokonano podziału brzegu F na N części. R zeczyw iste od

cinki brzegu ciała zastąpiono tzw . elem entam i brzegowymi. Zbiór w szy

stkich elementów brzegowych Fn stanowi m odel geom etryczny analizowa

nego ciała. Sytuację tę ilustruje dla zagadnień dwuwym iarowych rys.1.2.

(25)

24 1. W stęp

W najprostszym przypadku elementy brzegowe są odcinkami linii prostych.

M ożna jednak również używać elementów wyższego rzędu, jak krzywe sto

pnia drugiego, splajny itd. Elem enty specjalne, jak odcinki okręgu, elipsy, odcinki powierzchni kulistej itp., pozwalają na bardzo dokładne odwzoro

wanie kształtu ciał zbudowanych z wymienionych elementarnych obiektów geom et ry czny ch.

Z każdym z elementów brzegowych związany jest lokalny układ współ

rzędnych krzywoliniowych zaczepiony w środku elementu. Dla zagadnień płaskich jest to układ jednowym iarowy, w którym oś ( przebiega wzdłuż elementu, oraz -1 < £ < 1. D la zagadnień przestrzennych wym agany jest układ dwuwym iarowy ^ r).

Przebieg tem peratury i gęstości strumienia ciepła w obrębie elementu aproksymuje się, wykorzystując lokalne funkcje kształtu. Zależnie od przy

jętej aproksymacji rozróżnia się:

• elementy stałe, dla których zakłada się stałość tem peratury i strumie

nia ciepła w elemencie, rys. 1.3

T( o = Ti =

t2

(1.11)

9(0 = 9i = 92 (142)

Związek pom iędzy globalnym i lokalnym układem współrzędnych dla elementów stałych w yznaczają liniowe funkcje kształtu rozpięte na wę

złach 1 i 2

* (0 = |(1 - 01 + |(1 + 0 2 (1.13)

j/ ( 0 = ^ ( i - 0 » i . +5 ( 1 + Osfe (L 1 4 )

W arto w tym miejscu zauważyć, że elementy stałe należą do grupy elementów subparametrycznych, tzn. takich, dla których rząd aprok

symacji funkcji jest niższy niż rząd aproksymacji brzegu.

• elem enty liniowe, dla których tem peratura i gęstość strumienia ciepła aproksymowane są liniowo pom iędzy wartościami w węzłach 1 oraz 2, rys. 1.4

T( o = ! (1 - 0 I i + 1(1 + 0 1 * (1.15) 9(0 = J ( 1 - 0 «1 + 1 (1 + 0 9» (1-16)

Transformacja z układu lokalnego do globalnego przebiega zgodnie z równaniami (1.13) i (1.14). Ponieważ tym razem rząd aproksymacji funkcji jest taki sam jak rząd aproksymacji brzegu, elementy liniowe są elementami izoparametrycznymi.

(26)

1.1. Podstawy m etod y elementów brzegowych 25

A y

Rys. 1.2. M odel geometryczny ob- Rys. 1.3. Elem ent stały w global- szaru - elementy brze- nym i lokalnym układzie

gowe współrzędnych

Fig. 1.2. Geom etrical model o f Fig. 1.3. Constant element in glo- the region - boundary bal and local co-ordinate

elements system

• elementy kwadratowe, dla których przebieg tem peratury i strumienia ciepła w elemencie jest kwadratową funkcją lokalnej współrzędnej ( . Dla elementów kwadratowych wprowadza się więc trzeci węzeł, zw ykłe umieszczony w środku elementu, rys. 1.5

*(0

⁼ | ć ( ć - i ) * i + ( 1 - Ć 2)*3 + ^ ( £ + l)'*a (1-17)

y

(0 = 2 ^{^ ~} ^{+ 1)3/2} ⁽ 1 ^- 18 ⁾

T (0 = + ( 1 - ? ) T 3 + \ t ( t + l ) T 2

(1.19)

9 (0 = k ć - l ) 9 i + ( l - £ 2) 9 3 + k ( £ + l)92 (1.20)

W ięcej szczegółów na tem at różnych typów elementów brzegowych, w tym również na tem at elementów superparametrycznych i elementów nie

ciągłych, znaleźć można np. w monografiach [9] i [12].

W zapisie m acierzowym aproksymacje tem peratury i strumienia ciepła przyjmują dla poszczególnych typów elementów brzegowych następującą ogólną postać

(1.21) ( 1.22)

(27)

26 I. W stęp

gdzie; T „ - macierz kolumnowa zawierająca wartości tem peratury w punktach węzłowych należących do n-tego elementu brzegowego,

q n - macierz kolumnowa zawierająca wartości strumienia cie

pła w punktach węzłowych należących do n-tego elementu

, brzegowego,

- macierz wierszowa zawierająca lokalne funkcje interpola

cyjne.

Rozm iar poszczególnych m acierzy kolumnowych i wierszowych zależy od typu elementu brzegowego i wynosi: dla elementu stałego 1, dla liniowego 2, dla kwadratowego 3.

Rys. 1.4. Elem ent liniowy w glo- Rys. 1.5. Elem ent kwadratowy w balnym i lokalnym ukła- globalnym i lokalnym li

dzie wspólrzęinych kładzie współrzędnych Fig. 1.4. L in e a r element in global Fig. 1.5. Quadratic element in glo-

and local co-ordinate sy- bal and local co-ordinate

stem system

1.1.3. M acierze w p ły w u

D zieląc brzeg F na elementy brzegowe, zastąpiliśmy całki w równaniu (1.10) sumami całek po poszczególnych elementach Fn

N . N

k Ci Ti + L (f ' T d V - = E /, ( i .23)

n = l '/l'" n = i ■'l >•

Z kolei aproksymacje (1.21) i (1.22) pozwalają na przekształcenie rów

nania (1.23) do postaci

N . N .

kaTi + Y,

t „ / # T <?*

d\\ =

^{q» /}

*T T (lv«

^{( 1 -2'1)}

n=l n=l

(28)

1.1. Podstawy m etod y elem entów brzegowych 27

Całki wzdłuż poszczególnych elementów brzegowych T n oblicza się nu

merycznie w lokalnym układzie współrzędnych, po uwzględnieniu Jakobianu transformacji pom iędzy globalnym i lokalnym układem współrzędnych. Dla każdego elementu brzegowego uzyskuje się dw ie macierze wierszowe o roz

miarze jak macierz 4>T

(1.25)

(1.26)

Ponieważ rozwiązanie podstawowe zależy od odległości r punktu obser

wacji Yz i punktu bieżącego Yb, niektóre z całek są całkami osobliwym i.

Dzieje się tak wtedy, gdy punkt obserwacji należy do elementu brzego

wego, po którym całkujemy, por. rys.1.6 i rys.1.7. Podczas gdy całki re

gularne m ogą być z wysoką dokładnością obliczane zw ykłym i kwadratu

rami Gaussa [74], całki osobliwe wym agają osobnego traktowania. Stoso

wane są m iędzy innymi specjalne transformacje (np. transformacja Tełlesa [79]), które grupują w ęzły kwadratury w okolicach osobliwości. A ltern aty

wnie całki te można obliczać analitycznie lub pośrednio na podstawie sumy całek regularnych. Szczegóły znaleźć m ożna w monografiach M E B .

h" = / #VdT„

J

r„

g.n = [ $

t

t dTn*

JTn

Rys. 1.6. Całkowanie po elemen- Rys. 1.7. Całkowanie po elemen

cie regularnym cie osobliwym

Fig. 1.6. Integration along regu- Fig. 1.7. Integration along singu-

lar element lar element

Podstawiając równania (1.25) i (1.26) do równania (1.24), otrzym ujem y dla każdego punktu obserwacji

N N

k

c, T,

+ h in

t„ = £ Sm (1 -27)

n = 1 71=1

Punkt obserwacji i, zwany również punktem kolokacji, można w zasa

dzie wybrać dowolnie. Jednakże generując równania M E B , umieszcza się

(29)

28 1. W stęp

go kolejno w węzłach elementów brzegowych. O trzym uje się w' ten sposób tyle równań typu (1.27), ile jest punktów węzłowych. Zbierając odpowiednio tem peratury węzłowe w macierz kolumnową T , węzłowe strumienie ciepła w macierz kolumnową Q , macierze wierszowe zaś h ,n oraz g m w tzw. macierze wpływu, odpowiednio H i G , otrzym uje się ostatecznie następujący układ równań będący dyskretną form ą równania całkowego

H T = G Q (1.28)

U kład ten jest bezpośrednim związkiem m acierzowym pom iędzy tem pe

raturam i i gęstościami strumienia ciepła na brzegu ciała.

1.1.4. Uwzględnianie warunków brzegowych

Warunki brzegowe narzucają w każdym punkcie węzłowym wartość tem pe

ratury lub strumienia ciepła. Dodatkowo dla warunku brzegowego trzeciego rodzaju występuje liniowy związek pom iędzy tym i wielkościami. General

nie stwierdzić zatem należy, że tylko N wielkości w układzie równań (1.28) jest znanych z warunków brzegowych, pozostałych zaś N należy wyznaczyć przez rozwiązanie systemu (1.28).

Przegrupowując odpowiednio kolumny w układzie (1.28) i mnożąc prawą stronę przez wielkości znane z warunków brzegowych, otrzym ujem y ostate

cznie

A X = f (1.29)

gdzie: A - macierz główna układu,

X - macierz kolumnowa niewiadomych, f - macierz kolumnowa w yrazów wolnych.

Elem enty m acierzy kolumnowej niewiadomych oraz współczynniki ma

cierzy głównej zależą od rodzaju warunku brzegowego w danym węźle

!

^{— G i n} dla warunku brzegowego I rodzaju

H in dla warunku brzegowego II rodzaju (1.30)

H in — h G in dla warunku brzegowego II I rodzaju

{

^ąn dła warunku brzegowego I rodzaju

(1.31)

T n dla warunków brzegowych II i I I I rodzaju

gdzie: h - współczynnik wnikania ciepła.

M acierz główna układu równań (1.29) jest macierzą pełną i niesyme

tryczną. Rozwiązanie uzyskuje się klasycznymi metodam i numerycznymi, jak np. m etodą elim inacji Gaussa. P rzy większych rozmiarach macierzy A stosuje się podział m acierzy na bloki.

(30)

1.1. Podstawy m etod y elem entów brzegowych 29

1.1.5. Rozwiązanie w punktach wewnętrznych ciała

Po rozwiązaniu równania (1.29) i wyznaczeniu rozkładu tem peratury i stru

mienia ciepła w zdłuż brzegu ciała m ożna w yznaczyć tem peraturę w do

wolnie wybranym punkcie wewnętrznym . O peracja ta nie w ym aga żadnego dodatkowego podziału obszaru i polega na przemieszczeniu punktu i w rów

naniu (1.27) w nowe położenie. N ależy więc pow tórzyć całkowanie po brzegu ciała, tym razem z punktu wewnętrznego obszaru. W arto w tym miejscu zaznaczyć, że wszystkie obliczane w tym etapie całki są całkami regular

nymi i utworzenie wewnętrznych m acierzy wpływu H „ i G „ nie przedstawia większych trudności. Pam iętając o tym , że dla wszystkich punktów wewnę

trznych w spółczynniki c, są równe 1, tem peraturę T w w tych punktach w y  znacza się przez prostą operację mnożenia m acierzy w pływ u przez brzegowe wartości tem peratury i strumienia ciepła

T „ = - H „ T +

GWQ

(1.32)

K olejn e etapy rozw iązyw ania zagadnień brzegowych, opisanych równa

niami różniczkow ym i, za pom ocą M E B przedstawiono na rys.1.8.

1.1.6. P ola źródłowe

Dla zagadnień brzegowych opisywanych równaniem Laplace’a równanie (1.9) zawiera tylko całki po brzegu T obszaru 0 i przekształca się do postaci (1.10). D yskretyzacji podlega w ięc jed yn ie brzeg obszaru. W zagadnieniach ze źródłam i (rzeczyw istym i bądź fikcyjnym i) oprócz całek brzegowych w równaniu (1.9) pojaw ia się również całka p o całym obszarze fi. W pierw szych pracach traktujących o M E B , np. [8], [13], całki te obliczano wprost poprzez dyskretyzację obszaru. Jak ju ż jednak wspomniano, proces dys

kretyzacji obszaru jest kłop otliw y i pracochłonny, szczególnie w przypadku obszarów trójw ym iarow ych. Ponadto, całkowanie po wnętrzu obszaru należy wykonać tyle razy, ile jest brzegowych punktów kolokacji, co w zdecydowany sposób odbija się na efektywności m etody. Jeśli wewnętrzne źródła ciepła qv opisane są znaną funkcją położenia, całka po wnętrzu obszaru nie w pro

wadza żadnych nowych niewiadomych. M E B traci jednak swą największą zaletę.

Gipson w swej monografii [20] poświęconej zastosowaniu M E B do roz

wiązywania równania Poissona proponuje wykorzystanie m etody M on te Car- lo do obliczania całki po obszarze fi. O prócz licznych przykładów numery

cznych autor zam ieszcza wydruk programu kom puterowego realizującego zaproponowany algorytm .

(31)

I. W stęp

:{0

, i

zagadnienie brzegowe 1

zasada wzajemności | V

równanie calkpwe I_________

dyskretyzacja

dyskretna postać róionania calowego

warunki brzegowe j

w

ukfad równań algebraicznych

r o z w ią z a n ie

wielkości brzegowe

I

uńetkości wewnątrz obszaru

Uy s. I.S. // upij rn.:iri(t~i/ini ma ziif/ailini ii b r ;i (jo irijrli :n po m oru M i . l i l-ig . I.S. S u m s.sin shifii.s o f s o lriiifi b ouiiild ru p ro b liiiis by lili.M

(32)

1.2. M e to d y transform acji całek . 31

Całka po całym obszarze w równaniu (1.9) m oże być w wielu sytuacjach praktycznych przekształcona w całki brzegowe. T y m samym udaje się unik

nąć dyskretyzacji wnętrza obszaru. Zagadnienia związane z taką właśnie transformacją całki po obszarze w równoważne je j całki brzegowe są treścią niniejszej pracy.

1.2. M eto d y transformacji całek po obszarze w całki brzegowe

M etody transformacji omówione w niniejszym rozdziale dotyczą jedyn ie tych całek po obszarze, które występują w M E B , a precyzyjniej w równaniu (1.9). Całki te m ają następującą postać

D 0 = [ T ’ qv d n (1.33)

Jn

Stosowane dotychczas w M E B m etody transformacji całek (1.33) można podzielić na trzy główne grupy.

Pierw sza grupa m etod związana jest z wykorzystaniem rozwiązań szcze

gólnych równania Poissona, [2], [17]. Rozwiązanie szczególne, oznaczone przez T , spełnia oczywiście równanie (1.6). Zatem

9„ = -J fc V 2f (1.34)

Podstawienie zależności (1.34) do (1.33) um ożliwia wykonanie całkowa

nia przez części

D 0 = - k [ T ‘ V 2f d n =

Ja

= - k

J

^{T V}²T *< fn +

j

( r mq - q‘ T

)

^{d r} ^(1.35)

Po uwzględnieniu definicji rozwiązania podstawowego i skorzystaniu z własności filtracji delty Diraca powyższe równanie zapisuje się następująco

D 0 - - k Ci f i + J ( r * q - q* f ) d r (1.36)

Dla przypadku gdy rozwiązanie szczególne nie jest znane (a tak jest na ogół) Nardini i Brebbia zaproponowali na początku lat 80. [40], [41] globalną interpolację członu źródłowego qv

J

*• = £ / , « ; (1.37)

3

(33)

32 1. W stęp

gdzie: f j - funkcje interpolacyjne, a j - współczynniki aproksymacji,

J - liczba wyrazów w wyrażeniu aproksymacyjnym.

Funkcje interpolacyjne f j dobiera się tak, aby m ożna było łatw o znaleźć rozwiązanie szczególne równania

co ostatecznie prowadzi do następującego wyrażenia na całkę (1.33)

W arto w tym miejscu podkreślić, że zapewnienie odpowiedniej dokład- p .sci m etody wym aga, aby liczba wyrazów J w aproksymacji (1.37) była większa niż liczba węzłów brzegowych. Technika ta, znana obecnie jako podwójna zasada wzajemności P Z W (ang. D ual R eciprocity M ethod), w yko

rzystuje więc oprócz węzłów brzegowych również tzw . bieguny wewnętrzne (ang. internal p o lts). Liczba tych punktów oraz ich lokalizacja powinny oczy

wiście uwzględniać kształt funkcji qv. Innym i słowy, bieguny wewnętrzne powinny być głównie lokowane w obszarach o strom ym przebiegu funkcji źródła. Dokładniejsze informacje o P Z W znaleźć m ożna w licznych publika

cjach, np. [11], [64], [65], [67], [81], [82]. W szystkie one podkreślają ogólność m etody, jej efektywność numeryczną oraz wyjątkow o prostą implementa

cję komputerową. Jednocześnie przeprowadzone obliczenia testowe potw ier

dzają jej wysoką dokładność. N a początku 1992 roku ukazała się monografia [66] poświęcona m atem atycznym podstawom i praktycznym zastosowaniom podwójnej zasady wzajemności.

Druga grupa metod, zaproponowana przez Tanga [76], opiera się na za

stąpieniu aproksymacji (1.37) szeregiem Fouriera funkcji qv. Rozw inięcia w szereg dokonuje się w obszarze nie m niejszym niż fi. Ponieważ wyznaczenie współczynników tego rozwinięcia wym aga w zasadzie całkowania po wnętrzu obszaru, autor proponuje wprowadzić do rozważań obszar pom ocniczy fi o prostym kształcie, który zawiera w sobie obszar fi. Prosty kształt obszaru pomocniczego, np. prostokąt dla zagadnień dwuwymiarowych, czy prosto

padłościan dla zagadnień trójwym iarowych, pozwala na analityczne w yzna

czenie współczynników Fouriera. Rozw inięcia dokonuje się zatem faktycznie (1.38)

J

(1.39)

J

(34)

1.2. M e to d y transform acji całek . 33

w ’’ nieco za dużym ” obszarze f i . Następnie, wykorzystując ortogonalność funkcji bazowych i postępując podobnie jak w P Z W autor transformuje całkę (1.33) w szereg całek brzegowych.

M etod a opracowana w oryginale dla zagadnień potencjalnych oraz za

gadnień sprężystości była również stosowana przez Itagaki i Brebbię [25] do modelowania procesów dyfuzji neutronów. Jednakże ze względu na skompli

kowaną im plem entację komputerową nie zyskała ona dotąd większej popu

larności.

Trzecia grupa m etod związana jest z tzw . tensorem Galerkina [9], [12], [15], [16], [78]. Pozw ala ona na transformację całki (1.33) w całki brzegowe dla wybranych typów funkcji qv. Rozw inięciem tej koncepcji jest technika opracowana przez autora rozprawy pod koniec lat 80. [10], [50], [51], [56], [57], [60] i znana w literaturze anglosaskiej jako M ultiple R eciprocity M ethod (wielokrotna zasada wzajemności). Pozw ala ona modelować nie tylko pro

cesy cieplne, ale również wiele innych zagadnień inżynierskich opisywanych równaniami różniczkowym i innych typów. M iędzy innymi zastosowanie w ie

lokrotnej zasady wzajemności (W Z W ) do rozwiązywania równania Helmhol- tza om ówione jest w pracy [58], natomiast obliczanie wartości własnych tego równania w pozycjach [28], [29], [30] i [31]. Pow er H. i Power B.F. w ykorzy

stali W Z W do analizy zjawisk hydrodynam icznych [68], [69], wykazując, że uzyskane w ten sposób rozwiązanie jest całkową reprezentacją rozwiązania analitycznego. Neves i Brebbia rozszerzyli zastosowanie m etody na zagad

nienia sprężystości [43], zaś Neves, Brebbia, W róbel i Nowak na problem y term osprężystości dla zagadnień stacjonarnych [44], [45], [46] i niestacjonar

nych [47], [48]. Swego rodzaju podsumowaniem badań nad wykorzystaniem W Z W do modelowania zagadnień sprężystości i term osprężystości jest praca doktorska Nevesa [42]. Ostatnie, znane z literatury, zastosowania W Z W to problem y dyfu zji neutronów analizowane przez Itagaki i Brebbię [21], [22], [23], [24], [26] oraz zagadnienia drgań harmonicznych w cienkich sprężystych płytach analizowane przez V. Sladka, J. Sladka i M . Tanakę [72].

M etod a W Z W operuje tzw . rozwiązaniami podstawowymi wyższych rzę

dów i przekształca całkę (1.33) w sposób rekurencyjny (stąd nazwa m etod y).

W rezultacie otrzym uje się w granicy równanie całkowe, które jest dokład

nym i w pełni brzegow ym sformułowaniem problemu.

Podstaw y W Z W wyjaśnione są w następnym rozdziale, podczas gdy dal

sze części pracy dotyczą zastosowania m etody do rozw iązyw ania wybranych zagadnień inżynierskich.

Początkow o m etoda wykorzystywana była do rozw iązyw ania liniowych zagadnień brzegowych. Ostatnio [53], [54], [55] pole zastosowań W Z W zo

stało poszerzone o klasę nieliniowych problemów cieplnych. W ym a ga to je dnak jeszcze pewnego usystematyzowania i dalszych intensywnych badań.

(35)

R ozdział 2

P o d staw y wielokrotnej zasady wzajemności

2.1. Zależności ogólne

W niniejszym rozdziale przedstawiono ideę wielokrotnej zasady wzajem ności. Jak wspomniano ju ż wcześniej, technika ta pozwala przekształcić całkę po obszarze występującą w równaniu całkowym (1.9) w równo

ważne jej całki po brzegu ciała. Rozważania dotyczą całki z uogólnionym członem źródłow ym , który oznaczono przez 6(0*. Człon ten może być rze

czyw istym członem źródłow ym lub może zawierać funkcje reprezentujące źródła fikcyjne. W szczególności dla ustalonego przewodzenia ciepła mamy

bW = <h (2.1)

Górny indeks ’’ o ” został wprowadzony celem ujednolicenia notacji, co wyjaśniono w dalszej części pracy. Podobnie, indeksem ” o” opatrzono roz

wiązanie podstawowe zdefiniowane równaniem (1.3). T y m samym całkę (1.33) po obszarze ft można zapisać jako

D 0 = f rr m l/0 )dn (2.2)

Jn

Początkow o zakłada się również, że funkcja //0) jest znaną i gładką funkcją współrzędnych. Inne przypadki przeanalizowano w następnych roz

działach pracy.

R ozpocznijm y od wprowadzenia Iz w. rozwiązaniu podstawowego pierw

szego rzędu, czyli funkcji 7" (l) związanej z klasycznym rozwiązaniem pod

stawowym zależnością

y 2 y - ( l ) _ y..(0) ( 2 . 3 )

Jednocześnie analog strumienia ciepła, czyli funkcja </**'*, wynika z równania O T ' o

(36)

2.1. Zalcżn ości ogóln e 35

W ykorzystując sym etryczną tożsamość Grcena, można przedstawić całkę

Da

w sposób następujący

D 0 = [ T * (0)

6(0)

d n =

f

^V²^{r * ( ) )}⁶^{(0 )rfn =}

Jn J n

_ i f (7’*(1)U ,(0)_9*(1)6(°)) d r + f T ' {1)V 2bi0)d n =

k J r Jn

= 7 / (T*(1) w (0) - 9*(,) fc(0)) <*r + £>, -(2.5)

k J r

gdzie: u;*0' - analog strumienia ciepła będący pochodną norm alną uo

gólnionego członu źródłowego:

(

2

.

6

)

Całka D i w zależności (2.5) jest definiowana jako

D i = f T * (1) V26(0)<Zfi =

f

^(2.7)

prz y c z y m

j(i) = V 26(o> (2.8)

Ponieważ funkcja 6*°^ jest znaną funkcją współrzędnych, je j laplasjan może być obliczony analitycznie. D aje to nową, znaną funkcję położenia (/’ ). Z kolei różniczkując M1) wzdłuż normalnej zewnętrznej, otrzym ujem y funkcję stanowiącą odpowiednik strumienia ciepła

wW = ~ k i £ r (29)

Proces analitycznego różniczkowania znanej funkcji nie stanowi więk

szego problemu przy obecnym poziom ie techniki komputerowej. Szeroko do

stępne oprogramowanie, jak np. pakiety D erive [73] czy M athem atica [80], pozwalają zautom atyzować większość operacji ’ analitycznych’ , w tym ope

rację różniczkowania. Za ich pom ocą uzyskuje się również gotow e procedury komputerowe (w Fortranie, Pascalu itp .) obliczające wartości wyznaczonych pochodnych.

W arto w ty m miejscu zauważyć, że postać całki D\ jest form alnie identy

czna z postacią całki D „ (różn ią się one jed yn ie rzędem rozwiązania podsta

wowego i rzędem laplasjanu funkcji źródła). T y m samym całka m oże być przekształcona w sposób analogiczny jak całka D a . W rezultacie otrzym u

jem y

D , = i / ( T * (2) - q ' {2) 6( 1 )) dT + D 2 (2.10)

k Jr

(37)

36 2. Podstawy w ielokrotnej zasady wzajemności

Opisana procedura m oże być w oczyw isty sposób uogólniona. Definiu

jąc m ianowicie ciąg rozwiązań podstawowych wyższego rzędu określonych za

leżnością rekurencyjną

V2T .(;+D _ r .(i) j _ o, 1,2, . . . (2.11) ... d T '(l)

ciąg laplasjanów fu n k cji źródła

6(J) = V26(i_1) j = 1 ,2 , . . . (2.13) , > d b ^

" ('’ = - t a r (2'14)

oraz prowadząc proces transformacji zgodnie z rys.2.1, otrzym ujem y na

stępujący szereg całek brzegowych reprezentujący w yjściową całkę (2.2) i 00 r

D a = - / ( T * (J+1) wu ) - 9*<>+1>6< » ) d r (2.15) j =o •'r

N ależy pokreślić, że równanie (2.15) jest dokładną reprezentacją wyjścio

wej całki (2.2) - żadne uproszczenia nie zostały poczynione. Z m atem aty

cznego punktu widzenia szereg (2.15) jest rezultatem wykonania całkowania przez części nieskończoną liczbę razy.

Ostatecznie podstawiając (2.15) do równania (1.9) otrzym ujem y na

stępujące, w pełni brzegowe sformułowanie wyjściowego zagadnienia brze

gowego (1.6)

k Ci Ti +

J

q "{0) T d T =

J

T ’ (0)q d r -

jT O + i) w

U)

- qmV+ ' H u ) ) d r (2.16) j=o

Sformułowanie to, nie zawierając całek po obszarze, nie wym aga dyskretyzacji obszaru. W ten sposób zachowana została największa zaleta M E B , jaką jest prosta generacja siatki numerycznej. W ym agana jest natomiast znajomość na brzegu obszaru nie tylko samego członu źródłowego b^°\ ale również ciągu jeg o laplasjanów.

2.2. Dyskretyzacja brzegowego równania całkowego

Ogólne zasady dyskretyzacji równań całkowych M E B przedyskutowano już w podrozdziałach 1.1.2 oraz 1.1.3. Osobnego om ówienia w ym aga zatem je-

(38)

2.2. Dyskretyzacja brzegowego równania całkowego 37

[ 7 ^ 0 I

Rys. 2.1. Proces transform acji całek po obszarze w całki brzegowe - wielo

krotna zasada wzajemności

Fig. 2.1. Transform ation o f domain integrals into boundary ones - the M u ltip le R eciprocity M ethod

dynie dyskretyzacja członu źródłowego, w yrażonego w równaniu (2.16) sze

regiem całek brzegowych.

Jak pokazano wcześniej, podział brzegu T na N elem entów brzegowych pozwala na zastąpienie całki po brzegu ciała, sumą całek po poszczególnych elementach Fn