Zastosowanie metod identyfikacji w wybranych zagadnieniach przepływu biociepła

POLITECHNIKA ŚLĄSKA W GLIWICACH WYDZIAŁ MECHANICZNY TECHNOLOGICZNY

Katedra Wytrzymałości Materiałów i Metod Komputerowych Mechaniki

mgr inż. Marek Paruch

Zastosowanie metod identyfikacji w wybranych zagadnieniach

przepływu biociepła

Rozprawa doktorska

Promotor:

prof, dr hab. inż. Ewa Majchrzak

Gliwice 2005

, - - i w - ’,.

‘ I ' - . u - ; .

O f .

íZ-kAOO

W D a ^ )/ O Q '

Spis treści

SPIS TREŚCI

1 Ce! i zakres pracy...

4

2 Zadania odwrotne w przepływie biociepła - przegląd literatury

7

2.1. Równanie Pennesa... 7

2.2. Analiza wrażliwości... 8

2.3. Zadania odwrotne... 10

2.3.1. Metody gradientowe... 11

2.2.2. Metody ewolucyjne... 13

2.2.3. Metody hybrydowe...20

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

22

3.1. Wstęp... 22

3.2. Metoda elementów brzegowych ... 22

3.3. MEB z zastosowaniem wielokrotnej zasady wzajemności...24

3.4. Model numeryczny... 29

3.4.1. Elementy stałe... 29

3.4.2. Elementy liniowe - zadanie 2 D ... 30

3.5. MEB z zastosowaniem wielokrotnej zasady wzajemności - zadanie dodatkowe... 33

3.6. Model numeryczny - zadanie dodatkowe...37

3.6.1. Elementy stałe... 37

3.7. Obszary niejednorodne...38

4 Wyznaczanie rozkładu temperatury w tkance zaatakowanej nowotworem

42

4.1. Wstęp... 42

4.2. Opis matematyczny... 42

4.3. Analiza wrażliwości... 44

4.4. Wyniki obliczeń - zadanie 2D... 48

4.5. Wyniki obliczeń - zadanie 3D... 52

4.6. Analiza wrażliwości kształtu - zadanie 2 D ... 57

5 Identyfikacja parametrów termicznych oraz położenia i wiellcości

podobszaru nowotworowego

.

. 63

5.1. Wstęp... 63

5.2. Identyfikacja parametrów termicznych...64

5.2.1. Algorytm gradientowy...65

5.2.2. Algorytm ewolucyjny...66

5.2.3. Wyniki obliczeń... 67

5.3. Identyfikacja położenia i wielkości nowotworu... 78

5.4. Równoczesna identyfikacja parametrów termianych i geometrycznych nowotworu...84

6 Uwagi i wnioski końcowe...

.

...87

Literatura...

.

.88

Streszczenie.

.

.

..95

Summary.

.

.

96

_____________________________________________Spis treści

1. Cel i zakres pracy

ROZDZIAŁ 1

Ce! i zakres pracy

Pracę doktorską pt. „Zastosowanie metod identyfikacji w wybranych zagadnieniach przepływu biociepła” wykonano w Katedrze Wytrzymałości Materiałów i Metod Komputerowych Mechaniki w Politechnice Śląskiej w latach 2002-2005 pod kierunkiem prof. dr hab. inż. Ewy Majchrzak w ramach grantu promotorskiego nr 4 T07A 00626.

Zagadnienia związane z modelowaniem procesów cieplnych można podzielić na dwie grupy. Pierwsza z nich dotyczy tzw. zadań bezpośrednich, w których znane są zarówno równania i warunki brzegowo-początkowe (lub brzegowe) opisujące analizowany proces jak i wartości parametrów, które w tym modelu występują. Zadania, w których opis matematyczny jest niekompletny, np. brak jest pełnej informacji o warunkach brzegowych (warunku początkowym) lub też nieznane są wartości liczbowe pewnych parametrów, nazywane si^zadaniami odwrotnymi (zadaniami inwersyjnymi). Zadania odwrotne należą do tzw. problemów źle uwarunkowanych i niejednoznacznych, co powoduje, że metody ich rozwiązywania są znacznie trudniejsze niż w przypadku zadań bezpośrednich, a uzyskanie efektywnych rozwią^ń jest możliwe pod warunkiem, że dysponujemy dodatkowymi informacjami związanymi z przebiegiem modelowanego procesu, np. znamy wartości temperatury w wybranych punktach wewnętrznych obszaru.

Niniejsza praca doktorska dotyczy wybranych zadań odwrotnych występujących w problemach przepływu biociepła, a w szczególności tzw. zadań parametrycznych, które polegają na oszacowaniu wartości niektórych parametrów występujących w opisie matematycznym, zadań geometrycznych dotyczących identyfikacji położenia i kształtu podobszaru nowotworowego oraz tzw. zadań mieszanych związanych z równoczesną identyfikacją wartości parametrów termicznych i geometrycznych.

Rozpatrywano zdrową tkankę biologiczną oraz tkankę zaatakowaną nowotworem. Na podstawie wartości temperatury na powierzchni tkanki skórnej dokonano identyfikacji parametrów termicznych podobszaru nowotworowego oraz położenia i wielkości tego podobszaru. Do rozwiązania tak sformułowanego zadania zastosowano metody gradientowe, algorytmy ewolucyjne i hybrydowe. Zadania bezpośrednie i dodatkowe związane z analizą wrażliwości rozwiązywano wykorzystując metodę elementów brzegowych z zastosowaniem

/. Cel i zakres pracy

wielokrotnej zasady wzajemności. Rozważania teoretyczne zilustrowano licznymi przykładami obliczeń potwierdzającymi poprawność i efektywność prezentowanych algorytmów. Programy komputerowe opracowane w ramach rozprawy doktorskiej są programami autorskimi.

Tezę pracy można sformułować następująco:

Na podstawie rozkładu temperatury na powierzchni tkanki skórnej można wnioskować 0 obecności, parametrach termicznych i geometrycznych podobszaru nowotworowego.

Efektywnym narzędziem rozwiązania tak sformułowanego zadania odwrotnego są algorytmy gradientowe, ewolucyjne i hybrydowe połączone z metodą elementów brzego­

wych z zastosowaniem wielokrotnej zasady wzajemności

Praca składa się z sześciu rozdziałów, przy czym Avyniki badań własnych zamieszczono przede wszystkim w rozdziałach czwartym i piątym.

W rozdziale drugim dokonano przeglądu literatury związanego z zadaniami odwrotnymi przepływu ciepła w organizmach żywych. Omówiono równanie Pennesa opisujące rozkład temperatury w tkance biologicznej. Przedstawiono pokrótce metody analizy wrażliwości w ujęciu bezpośrednim i sprzężonym oraz przykłady zadań odwrotnych przepływu biociepła wraz z metodami ich rozwiązywania tzn. algorytmami gradientowymi 1 ewolucyjnymi.

Rozdział trzeci zawiera opis metody elementów brzegowych z zastosowaniem wielokrotnej zasady wzajemności. Metodę tę adaptowano na przypadek równania Pennesa, w którym składnik źródłowy jest zależny od temperatury. Pokazano również sposób zastosowania tego wariantu MEB do rozwiązywania zadań dodatkowych związanych z analizą wrażliwości. Rozważania przeprowadzono zarówno dla zadań płaskich (2D) jak i przestrzennych (3D), dla jednego oraz dwóch obszarów, na styku których przyjęto warunek brzegowy czwartego rodzaju.

W rozdziale czwartym (badania własne) wyznaczono rozkłady temperatury w obszarze tkanki zdrowej oraz tkanki zmienionej chorobowo. Przeprowadzono analizę wrażliwości otrzymanych rozwiązań ze względu na parametry termiczne i geometryczne podobszaru nowotworowego i sformułowano wnioski z przeprowadzonych badań.

Najważniejszą część rozprawy doktorskiej stanowi rozdział piąty. W rozdziale tym, na podstawie znajomości rozkładu temperatury na powierzchni tkanki skórnej dokonano identyfikacji parametrów termicznych podobszaru nowotworowego (odwrotne zadanie parametryczne), identyfikacji położenia i wielkości tego podobszaru (odwrotne zadanie geometryczne) oraz równoczesnej identyfikacji położenia, wielkości i parametrów

/. Cel i zakres pracy

termicznych podobszaru nowotworowego (odwrotne zadanie mieszane). Tak sformułowane zagadnienia rozwiązywano stosując algorytmy gradientowe, ewolucyjne i hybrydowe.

Obliczenia przeprowadzono zarówno dla dokładnych jak i zaburzonych pomiarów temperatury na powierzchni tkanki skórnej. Dokonano porównania efektywności i dokładności zastosowanych algorytmów identyfikacji.

Rozprawę zamyka rozdział szósty, w którym sformułowano uwagi i wnioski końcowe, spis literatury oraz streszczenia w języku polskim i angielskim.

2. Zadania odwrotne przepływie biocieph - przegląd literatury

ROZDZIAŁ 2

Zadania odwrotne w przepływie biociepła

przegiąd literatury

2.1. Równanie Pennesa

Podstawowym równaniem opisującym przepływ ciepła w obszarze tkanki biologicznej jest równanie Pennesa [44, 89, 99, 100, 101]

^ V [ X ( T ) W T ( x , / )] + (j:, /) + e _ , (2.1)

gdzie X [W/mK] jest współczynnikiem przewodzenia tkanki, c [J/m^K] - jej ciepłem właściwym odniesionym do jednostki objętości (iloczynem ciepła właściwego i gęstości masy). Równanie to zawiera w sobie dwa składniki źródłowe związane z perfuzją krwi Qp^rf i przemianami metabolicznymi Q„ei- Składnik perfiizyjny określa następująca zależność

QfKrf ~ ~ 0 . (2-2)

gdzie Gb [m^ krwi/s/ m^ tkanki] jest współczynnikiem perfuzji, cb - objętościowym ciepłem właściwym krwi, T^b - temperaturą krwi.

Składnik Q„ei traktuje się z reguły jako wartość stałą, ale zależną od warunków, w jakich w danej chwili znajduje się organizm (odpoczynek, lekki lub cięższy wysiłek fizyczny).

Wartość gpei/zmienia się w granicach 245-24500 [W/m^].

Dla stanu ustalonego równanie (2.1) przyjmuje postać

v [ ? t ( r ) v r ( x ) ] + 0 ^ ^ ( x ) + a „ = 0 (2.3) Równanie Pennesa dobrze przybliża rzeczywiste warunki przepływu ciepła w tkance zasilanej dużą liczbą równomiernie rozłożonych włoskowatych naczyń krwionośnych [97].

Jeżeli rozpatrujemy problem, w którym należy uwzględnić oddziaływanie cieplne pojedynczych dużych naczyń występujących z reguły parami (żyła - tętnica), to równanie Pennesa uzupełnia się dodatkowymi równaniami opisującymi przepływ ciepła w naczyniach krwionośnych [20, 30, 31, 58, 59, 97, 102]. Równania te wynikają z bilansów energii dla pojedynczych naczyń i przypominają w pewnym stopniu równania znane z teorii wymienników ciepła.

2. Zadania odwrotne w przepływie biociepła - przegląd literatury

Równanie Pennesa poddawano pewnym modyfikacjom [97, 98], ale obecnie jego postać bazowa jest powszechnie akceptowana, natomiast istotnym problemem jest właściwy dobór współczynników w tym równaniu, a w szczególności Gb i Q„,e, dla różnych typów tkanek zarówno zdrowych, jak i chorych [42, 93].

W teorii przepływu biociepła rozważa się również problemy, w których tkankę traktuje się jako obszar niejednorodny (wielowarstwowy). Takie podejście zaprezentowano między innymi w [35, 49-56] gdzie w tkance skórnej wyróżniono naskórek, skórę właściwą i obszar podskórny (parametry tych podobszarów były oczywiście zróżnicowane). Model matematyczny procesów cieplnych w obszarach niejednorodnych opisuje układ równań Pennesa uzupełniony warunkami ciągłości na granicach podobszarów oraz pozostałymi warunkami brzegowymi (brzegowo-początkowymi).

2.2. Analiza wrażliwości

Parametry termofizyczne występujące w równaniu Pennesa takie jak współczynnik przewodzenia ciepła, współczynnik perfuzji czy też składnik związany z metabolizmem są trudne do określenia i cytowane w literaturze wartości tych parametrów różnią się w sposób istotny (zależą one z całą pewnością od cech osobniczych takich jak płeć, wiek, a nawet wykonywany zawód). Składnik związany z metabolizmem jest ponadto zależny od warun­

ków, w jakich znajduje się organizm (spoczynek, ruch, zimno). Istotne różnice pojawiają się również w geometrii analizowanych obszarów, np. grubości naskórka, skóry właściwej.

W modelowaniu procesów przepływu biociepła ważnym problemem staje się więc określenie wpływu zmienności parametrów występujących w opisie matematycznym na otrzymywane wyniki obliczeń. Można tutaj wykorzystać metody analizy wrażliwości [16, 22, 23, 24, 25,35, 36, 39] zarówno w ujęciu bezpośrednim jak i sprzężonym.

Prace [21, 35, 36, 49-56, 77] dotyczą aplikacji analizy wrażliwości, a w szczególności jej wykorzystania do oceny wpływu zaburzeń parametrów wewnętrznych i zewnętrznych procesu na przebiegi pól temperatury w tkance biologicznej. Można więc badać skutki zmian takich parametrów jak ciepło właściwe tkanki, jej współczynnik przewodzenia, współczynnik perfuzji, wydajność źródła metabolicznego itp. na rozkłady temperatury w tkance. Przedmiotem analizy wrażliwości również może być wpływ zmian warunków brzegowych i początkowych na przebieg procesu [36, 76].

Jak wspomniano wcześniej, parametryczną analizę wrażliwości rozpatruje się bądź w ujęciu bezpośrednim bądź w ujęciu sprzężonym [22, 23, 35, 39]. Podejście bezpośrednie sprowadza się do konstrukcji zadania dodatkowego wynikającego z różniczkowania równań i warunków opisujących proces względem wyróżnionego parametru.

2. Zadania odwrotne przepływie biociepła - przegląd literatury

Przykładowo, równanie Pennesa dla stałych wartości c i X oraz obszaru ID jest postaci

0 < , < i : = + (2.4)

Równanie to uzupełnimy warunkami

x = 0: T { x , t ) = T^

x ^ L : q{x, /) = - X— ( 2. 5)

ox / = 0: T { x , t ) = T,

gdzie Th jest temperaturą brzegową, Tq - temperaturą początkową tkanki. Z fizycznego punktu widzenia zadanie (2.4), (2.5) może dotyczyć symulacji procesu nagrzewania tkanki przez źródło zewnętrzne.

Załóżmy, że będziemy badać wrażliwość tego procesu ze względu na zmianę parametru Qmei- Po zróżniczkowaniu równań i warunków po tym parametrze otrzymujemy następujący model wrażliwości

d U i x , t ) d ^ U ( x , t ) , .

0 < ^ < i : c — ^ = 1— A ^ - G , c , i / ( x , / ) + l (2.6)

oraz

>r = 0; t/(x,/) = 0 , d U ( x , t)

x = L : - X ---- ^— ¿ = 0 (2.7)

ox t = 0: U { x , t ) = 0

gdzie U = dT W omawianym przypadku model bazowy i model wrażliwości nie są ze sobą sprzężone i można rozpatrywać je odrębnie. Sytuacja taka nie jest typowa, z reguły obydwa modele są sprzężone i uzyskanie rozwią^nia zadania dodatkowego wymaga znajomości rozwiązania zadania podstawowego. Gdyby np. rozważanym parametrem było cie{>ło właściwe c, to odpowiedni model wrażliwości względem tego parametru byłby postaci

d U ( x , t ) d ^ U ( x , t ) , , d T ( x , t )

0 < x < L : c — ^ = X— ^ ^ - G , c , U { x j ) - c — ! ^ (2.8)

oraz

;c = 0; t/(;c,i) = 0

(2.9)

ÔX

t = 0: U { x , t ) = T,

gdzie U = d T l d c . Jak widać, w równaniu (2.8) pojawia się dodatkowy składnik źródłowy wynikający z zadania podstawowego.

2. Zadania odwrotne w przepływie biocieph - przegląd literatury

Do rozważań wprowadzimy następujący funkcjonał 5F

5c q{0, t) di (2.10)

który określa ilość ciepła doprowadzoną do tkanki przez powierzchnię kontaktu z zewnętrznym źródłem w czasie /.

W metodzie bezpośredniej wariację tego funkcjonału względem np. parametru c oblicza się z zależności

5c dc dx dt

(

.

)

W metodzie sprzężonej definiuje się jedno zadanie brzegowo-początkowe, którego postać jest zdeterminowana postacią funkcjonału [22, 23]. Nie wnikając w szczegóły, dla zadania podstawowego (2.4), (2.5) i funkcjonału (2.10) zadanie dodatkowe definiuje się następująco

' dr ~ dx^ GsC, W{x,t)

^ ( 0 , x ) = - l (2.12)

V { L , x ) = 0 W{ x, 0) = 0 0 < x < L : c-

x = 0:

x = L : x = 0:

przy czym , czyli tutaj czas „biegnie” w drugą stronę. Należy w tym miejscu zwrócić uwagę na fakt, że w podejściu sprzężonym bez względu na liczbę rozpatrywanych parametrów definiujemy tylko jedno zadanie dodatkowe, które nie ma żadnej interpretacji fizycznej. Rozwiązanie zadania bezpośredniego i dodatkowego pozwala na wyznaczenie wariacji funkcjonału (2.10) względem każdego z parametrów. Stosując odpowiednie wzory [22, 23, 35] np. wrażliwość funkcjonału (2.10) względem parametru c oblicza się na podstawie zależności

5^ r frr./ .

— = - W [ x , t - i ) — i;— ¿dxd/

6c J J ^ > dt (2.13)

2.3. Zadania odwrotne

Jeżeli znane są wartości wszystkich parametrów występujących w modelu matematycznym procesu przepływu biociepła i zadanie polega na wyznaczeniu rozkładów temperatury w analizowanym obszarze, to mówimy o tzw. zadaniu bezpośrednim {direct problem). Zadania, w których opis matematyczny jest niekompletny, np. brak jest pehiej informacji o warunkach brzegowych (warunku początkowym) lub też nieznane są wartości

2. Zadania odwrotne w przepływie biociepła - przegląd literatury

2, 6, 7, 8, 40, 79, 87, 94]. Zadania odwrotne należą do tzw. problemów źle uwarunkowanych i niejednoznacznych [32, 78, 95, 96] co powoduje, że metody ich rozwiązywania są znacznie trudniejsze niż w przypadku zadań bezpośrednich, a uzyskanie efektywnych rozwią^ń jest możliwe pod warunkiem, że dysponujemy dodatkowymi informacjami związanymi z przebiegiem modelowanego procesu, np. znamy wartości temperatury w wybranych punktach wewnętrznych obszaru.

Liczba prac poświęconych zagadnieniom odwrotnym przepływu biociepła jest stosunkowo niewielka [17, 18, 19, 33, 41, 48, 55]. Dotyczą one między innymi identyfikacji parametrów termofizycznych tkanki skórnej przy założeniu, że znany jest rozkład temperatury na jej powierzchni lub odtwarzania brzegowego strumienia ciepła na powierzchni tkanki skórnej, jeżeli znane są krzywe nagrzewania w punktach wewnętrznych tkanki. Zadania odwrotne rozwiązuje się wykorzystując najczęściej tzw. metody gradientowe [27, 92] oraz algorytmy ewolucyjne [3, 28, 75, 91]. W niniejszej pracy doktorskiej stosowano obydwa podejścia.

2.3.1. Metody gradientowe

Dla przykładu, rozpatrywać będziemy równanie Pennesa (2.4) opisujące rozkład temperatury w tkance skórnej traktowanej jako obszar jednowarstwowy (zadanie ID), którego parametiy termofizyczne są średnimi wartościami parametrów naskórka, skóry właściwej i obszaru podskórnego.

Założymy, że na powierzchni skóiy działa strumień ciepła czyli

x ^ 0 : g (0 ,/ ) = ę, (2.14)

natomiast na wewnętrznej powierzchni x = L : ę (x , t) = 0.

Temperatura początkowa tkanki skórnej jest również znana

/ = 0: T (x , 0 ) = T„ (2.15)

Rozwiązanie tak sformułowanego zadania bezpośredniego pod warunkiem, że znane są wartości wszystkich parametrów termofizycznych, pozwala między innymi wyznaczyć temperaturę na powierzchni skóry, a mianowicie

T j = T ( 0 , t f ) , / = 0 ,1 ,2 ,...,F (2.16) Załóżmy teraz, że na podstawie temperatury powierzchni skóry należy odtworzyć stałą wartość składnika metabolicznego Q„ei- W tym celu wprowadzamy kryterium najmniejszych kwadratów (funkcję celu)

S ( Q ^ ) ‘ t ( T ’ - T j ) " (2.17)

/-I

2. Zadania odwrotne w przepływie biociepła - przegląd literatury

W którym T/ oznacza znane wartości temperatury (por. wzór (2.16)) na powierzchni skóry (np. otrzymane z pomiarów), natomiast - wartości temperatury otrzymane z obliczeń dla arbitralnie założonej wartości składnika źródłowego Q„ei- Rozwiązanie tak sformułowanego zadania odwrotnego sprowadza się do znalezienia minimum funkcji S.

W przypadku zastosowania algorytmu gradientowego, funkcję S różniczkujemy względem nieznanej wartości Q„,ei, a następnie stosujemy warunek konieczny istnienia ekstremum, czyli

/-I

dT^ = 0 (2.18)

gdzie Q°„}est arbitralnie założoną wartością parametru natomiast dla A: > O wynikać będzie z poprzedniego kroku iteracji.

Funkcję 7(0, /^) rozwijamy w szereg Taylora z dokładnością do dwóch składników

dT^ ■ (2.19)

+ - y ^met ^mei}

gdzie oznacza wyznaczone z rozwiązania zadania bezpośredniego wartości

temperatury przy założeniu, że .

Wprowadzając (2.19) do (2.17), po prostych przekształceniach otrzymujemy

E

(

')‘

e . " = e i F r

(

.

)

gdzie

{ u > ) ‘ - dT^

(2.21)

q^,=qL

Wyznaczenie współczynników wrażliwości (2.21) wymaga zróżniczkowania równań tworzących opis matematyczny zadania bezpośredniego względem Q„e,- Różniczkowanie równania (2.4) względem tego parametru prowadzi do zależności (2.6), a różniczkowanie warunków brzegowo-początkowych daje

jc = 0: - 1 dx

dx

( ² . ²² )

/ = 0: U { x , ( ) = Q

Podsumowując, w każdym kroku iteracji należy rozwiązać zadanie bezpośrednie

2. Zadania odwrotne w przepływie biociepła - przegląd literatury

wyznaczyć nową wartość Q„,e,. Proces iteracyjny prowadzi do osiągnięcia żądanej dokładności lub do zadanej z góry liczby iteracji K.

Nawet tak prosty, przedstawiony wyżej przykład iteracyjnej metody gradientowej poszukiwania minimum uwidacznia kilka charakterystycznych cech tej grupy algorytmów [27, 91], Po pierwsze, rozwiązanie uzyskuje się poprzez wykonanie ciągu kroków (iteracji) wzdłuż odpowiednio wyznaczonych kierunków nazywanych kierunkami poszukiwań.

Kierunki te są tworzone na podstawie pochodnej (gradientu) funkcji celu, bądź drugiej pochodnej (hesjanu), jeśli uwzględnić większą liczbę składników rozwinięcia funkcji T w szereg Taylora (por. wzór (2.19)). Po drugie, poszukiwanie minimum funkcji celu rozpoczyna się od jednego punktu, który w pierwszej iteracji jest punktem startowym (w rozważanym przykładzie jest to punkt ). I po trzecie, punkt startowy musi leżeć w bliskim położeniu optimum, w przeciwnym razie algorytm może być rozbieżny lub też może zbiegać się do innego rozwiązania (minimum lokalnego).

2.3.2. Metody ewolucyjne

w chwili obecnej do rozwiązywania zadań identyfikacji bardzo często wykorzystywane są metody ewolucyjne, które naśladują procesy zachodzące w świecie organizmów żywych [3, 28, 75]. Metody te polegają na symulowaniu mechanizmów doboru naturalnego oraz dziedziczenia (teoria ewolucji) przez pokolenia potomne pewnych cech wspólnych po swoich rodzicach. Zgodnie z tą teorią największe szanse na przeżycie mają te osobniki, które są najlepiej przystosowane do środowiska. Szybki rozwój metod opartych na podejściu ewolucyjnym doprowadził do powstania nowej dziedziny informatyki, znanej obecnie jako obliczenia ewolucyjne.

W literaturze [75] można spotkać kilka różnych metod ewolucyjnych, jak np. algorytmy genetyczne, algorytmy ewolucyjne, strategie ewolucyjne, programowanie genetyczne i programowanie ewolucyjne. Choć różne metody ewolucyjne powstały niezależnie, obecnie różnice między poszczególnymi metodami wyraźnie się zacierają. Metody te nazywa się więc ogólnie technikami obliczeń ewolucyjnych, algorytmami ewolucyjnymi lub po prostu metodami ewolucyjnymi [91]. Aby podkreślić podobieństwo wszystkich metod czy algorytmów opartych na zasadach ewolucji, używa się wspólnego terminu programy ewolucyjne [75].

Można powiedzieć, że metody ewolucyjne charakteryzują się kilkoma cechami wspólnymi. W każdej z metod w sposób losowy generowana jest startowa populacja osobników, która podczas procesu ewolucji, podlega różnym operatorom ewolucyjnym, najczęściej krzyżowaniu i mutacji oraz procesowi selekcji. Każdy osobnik w populacji

2. Zadania odnrotne w przepływie biociepła- przegląd literatury

określany jest mianem chromosomu i złożony jest z pewnej liczby zmiennych projektowych (genów). Pojedynczy chromosom w populacji jest reprezentantem jednego rozwiązania zadania. Krzyżowanie polega na dowolnym wyborze par chromosomów z populacji oraz tzw. punktu krzyżowania, w którym dokonuje się wymiany materiału genetycznego (genów) pomiędzy wylosowanymi osobnikami. Mutacja polega na dowolnym wyborze chromosomów z całej populacji oraz zmianie wartości pojedynczego lub kilku genów w tych chromosomach. W procesie selekcji największe szanse „na przeżycie” mają osobniki (chromosomy) najlepiej przystosowane. W wyniku modyfikacji populacji osobników oraz procesu selekcji, w kolejnych krokach działania algorytmów ewolucyjnych otrzymuje się zwykle coraz to lepiej przystosowane osobniki, czyli lepsze rozwiązania rozważanego zadania.

Różnice pomiędzy poszczególnymi metodami ewolucyjnymi dotyczą m.in. sposobu reprezentacji i kodowania osobników, sposobu przeprowadzania procesu selekcji, kolejności występowania selekcji i operatorów genetycznych, parametrów algorytmu ewolucyjnego, które mogą być stałe lub zmienne w procesie ewolucji, sposobu uwzględniania ograniczeń w rozwiązywanych zadaniach identyfikacji itp.

W klasycznych algorytmach genetycznych, zwanych w literaturze prostymi algorytmami genetycznymi, do reprezentacji osobników w populacji używa się kodowania binarnego, natomiast jako operatory genetyczne stosuje się proste krzyżowanie i prostą mutację. Najczęściej stosowanym typem selekcji w algorytmach genetycznych jest metoda ruletki. W metodzie tej prawdopodobieństwo znalezienia się w populacji potomnej danego osobnika jest tym większe, im „lepsza” jest wartość jego funkcji przystosowania, przy czym wartość tej funkcji jest wartością optymalizowanej funkcji celu.

Programy ewolucyjne lub algorytmy ewolucyjne są uogólnieniem klasycznych algorytmów genetycznych. W algorytmach tych wykorzystuje się zwykle bardziej złożone struktury danych do reprezentacji populacji osobników. Nie mówi się tutaj o kodowaniu parametrów zadania, gdyż używana jest zazwyczaj reprezentacja zmiennoprzecinkowa lub całkowitoliczbowa, w której zmienne projektowe przetwarzane są w bezpośredni sposób. Stosuje się też nowe operatory genetyczne, tworzone np. na potrzeby rozwiązywanego zadania, które nazywa się operatorami ewolucyjnymi. Stosowanym typem selekcji w algorytmach ewolucyjnych jest wspomniana metoda ruletki lub inne rodzaje selekcji, jak np. selekcja turniejowa.

Algorytm ewolucyjny to inaczej program komputerowy, wykorzystujący pewną strukturę danych do reprezentacji chromosomów w populacji i przetwarzający tę populację w celu uzyskania lepszych rozwią^ń w czasie procesu optymalizacji.

W każdej iteracji takiego algorytmu przeprowadza się działania na populacji osobników

2. Zadania odwrotne w przepływie biociepła — przegląd literatury

P = cl\ ch2 ... ch^ , gdzie chi, dla i = 1, 2, m oznacza /-ty chromosom w populacji P o liczebności osobników m. Dowolny chromosom c/j, składa się z n genów geny, dla j = 1, 2, n, czyli chf = \genn gen^^ ... N a wartości poszczególnych genów geriij nałożone są dolne i górne ograniczenia genP < gen.j < g e n ° . Schemat blokowy algorytmu ewolucyjnego przedstawiono na rys. 2.1.

Rys. 2.1. Schemat blokowy algorytmu ewolucyjnego

W pierwszym etapie obliczeń tworzona jest populacja startowa chromosomów (zazwyczaj losowo), przy czym każdy chromosom chi reprezentuje jedno z możliwych rozwiązań problemu. Liczba osobników w populacji może być dowolna, natomiast liczba genów w chromosomie jest równa liczbie zmiennych projektowych.

Następnym etapem jest ocena wszystkich chromosomów w populacji, poprzez obliczenie wartości funkcji przystosowania dla każdego osobnika. Wartość tej funkcji zwana jest też przystosowaniem osobnika. W przypadku zadania identyfikacji, im większe przystosowanie osobnika, tym ma on większe szanse na „przetrwanie” w procesie selekcji, w którym eliminowane są zwykle osobniki najsłabsze. W zadaniach identyfikacji etap obliczania funkcji przystosowania jest etapem najdłuższym, gdyż dla każdego chromosomu w populacji należy rozwiązać jedno zadanie bezpośrednie.

Kolejny etap, po obliczeniu wartości funkcji przystosowania dla wszy'stkich osobników w populacji, to sprawdzenie warunku zatrzymania obliczeń. Warunek ten może mieć różną postać, w zależności od rozwiązywanego problemu. Zazwyczaj zatrzymanie algorytmu

2. Zadania odwrotne przepływie biociepła - przegląd literatury

następuje po zadanej z góry liczbie generacji, otrzymaniu satysfakcjonującego wyniku lub, gdy jego dalsze działanie nie poprawia wartości funkcji przystosowania.

W przypadku, gdy warunek zatrzymania obliczeń jest spełniony, następuje analiza parametrów najlepszego chromosomu z ostatniej populacji, którego geny zawierają dane dotyczące np. kształtu, bądź parametrów termofizycznych. W wyniku działania operatorów ewolucyjnych, które modyfikują poszczególne chromosomy i ich geny, najlepszy chromosom z ostatniej populacji nie musi być najlepszym osobnikiem w całym procesie ewolucji.

W przypadku nie spełnienia warunku zatrzymania, kolejnym krokiem jest stworzenie nowej populacji w wyniku selekcji chromosomów i działania operatorów ewolucyjnych.

Selekcja polega na wybraniu do następnego pokolenia chromosomów na podstawie obliczonych wartości funkcji przystosowania, przy czym największe szanse na utworzenie nowego pokolenia mają osobniki najsilniejsze, gdyż prawdopodobieństwo ich wyboru jest największe. Operatory ewolucyjne w postaci różnych odmian krzyżowań i mutacji modyfikują niektóre chromosomy oraz pojedyncze geny. O wyborze chromosomów do krzyżowania oraz genów do mutacji decydują wartości prawdopodobieństw odpowiednio krzyżowania i mutacji. Po utworzeniu nowej populacji następuje ocena wszystkich osobników tej populacji i cały proces ewolucji się powtarza, dopóki nie zostanie spełniony warunek zatrzymania obliczeń.

Poniżej przedstawiono i omówiono stosowane w pracy doktorskiej operatory ewolucyjne, a mianowicie:

• krzyżowanie proste,

• krzyżowanie arytmetyczne,

• mutacja równomierna,

• mutacja z rozkładem Gaussa,

• klonowanie.

Krzyżowanie proste polega na losowej zmianie dowolnie wybranej pary chromosomów (2.23). W tym celu losowany jest punkt krzyżowania z przedziału [1, w-l], gdzie n jest liczbą genów, w którym następuje zamiana genów pomiędzy wybranymi rodzicami

c\ = genu sen,,.

c h , = ^<^«22 ^««23 gen^ (2.23)

W wyniku krzyżowania prostego, dla punktu krzyżowania jak na lys. 2.2, chromosomy przyjmą postać

cĄ' = [g e « „ ge«,3 gen^, gen^

cH,=[gen^, gen^ gen^^ gen,, ge«,,. (2.24)

2. Zadania odwrotne w przepływie biociepła - przegląd literatury

CHROM OSOM 2

“I---r

J _________________

punkt krzyżowania i

C H R O M O S O M /

Rys. 2.2. Krzyżowanie proste

Krzytowanie arytmetycuie polega na utworzeniu dwóch osobników potomnych ch[

oraz c!^ (2.25), których geny są liniową kombinacją wartości genów dwóch chromosomów rodzicielskich (rys. 2.3)

ch^ = a*ci\ + { \ - a ) * ch^

= a * c } \

+(l-a)*c/^

gdzie ch\ oraz chi to osobniki rodzicielskie, natomiast a jest parametrem krzyżowania arytmetycznego, dobieranym w sposób losowy z przedziahi (O, 1). Operator ten ma charakter eksploracyjny i umożliwia szerokie przeszukiwanie przestrzeni rozwiązań.

(2.25)

C H R O M O SO M i

CH RO M O SO M 2

1

tt=0,33

C H R O M O S O M i'

CH RO M O SO M 2'

Rys. 2.3. Krzyżowanie arytmetyczne

M utacja równomierna polega na zmianie wartości jednego lub kilku genów w chromosomie (rys. 2.4). Każdy gen gen^ w chromosomie chi ma równe szanse na to, by ulec procesowi mutacji, zgodnie z prawdopodobieństwem jej wystąpienia. W celu dokonania mutacji dla każdego genu gentj losuje się liczbę z przedziahi (O, 1) i jeśli wylosowana liczba jest mniejsza lub rówTtia prawdopodobieństwu mutacji, to wartość danego genu

2. Zadania odwrotne w przepływie biociepła - przegląd literatury

zmieniana jest na wartość losową gen'^, spełniającą ograniczenia nałożone na gen [gen^ < geriy < g e n f) , czyli

ch^Ągen,, gen,^ gen,, gen,^

c/?; = [ge«i, geni2 gen,3 gen,, gen[^ (2.26)

C H R O M O S O M1

C H R O M O S O M i '

mutaqi podlega gen 2 oraz 5

j

Rys. 2.4. Mutacja równomierna

Mutacja z rozhtadem Gaussa polega na zmianie wartości jednego lub kilku genów w chromosomie, podobnie jak w mutacji równomiernej, z tym, że wartości genów zmieniane są z wykorzystaniem rozkładu Gaussa.

Klonowanie polega na migracji najlepszego osobnika z poprzedniej populacji, z określonym prawdopodobieństwem, do następnej, bez konieczności uczestnictwa w procesie selekcji. Operator klonowania może powodować pewne niebezpieczeństwo, które wiąże się z tym, że program ewolucyjny utknie w minimum lokalnym.

Selekcja turniejowa działa dwuetapowo, np. wykorzystując cieka\vy i efektowny pomysł uporządkowania [75], W metodzie tej wybiera się pewną liczbę k osobników i selekcjonuje najlepszego z tego ¿-elementowego zbioru do następnego pokolenia na podstawie funkcji przystosowania każdego z osobników. Proces ten powtarza się n razy, aż do uzyskania nowego pokolenia, gdzie n jest liczebnością populacji. Duża wartość k, tzw.

rozmiar turnieju, zwiększa napór selekcyjny tej metody.

Przykładową selekcję turniejową przedstawiono na rys. 2.5, dla ^ = 4.

© © 0

Rys. 2.5. Selekcja turniejowa

Dobór wszystkich parametrów algorytmu ewolucyjnego, jak liczba cłuomosomów

2. Zadania odwrotne w przepływie biociepła przegląd literatury

liczba osobników w turnieju), spoczywa na użytkowniku. Dobór optymalnych wartości tych parametrów nie jest łatwy i zależy zazwyczaj od konkretnego zadania oraz postaci optymalizowanej funkcji przystosowania.

Ideę algorytmu ewolucyjnego przedstawimy na przykładzie omówionym w poprzednim podrozdziale. Funkcję przystosowania można zdefiniować tak samo, jak w algorytmie gradientowym, czyli wzorem (2.17). Genem będzie tutaj składnik źródłowy spełniający ograniczenia <Q „„ ■ Ponieważ mamy tylko jedną niewiadomą do odtworzenia, więc w tym przypadku każdy chromosom ch, będzie składał się tylko z jednego genu, czyli

ch, = [gew,.], przy czym gen, =

Populację startową chromosomów można utworzyć losowo, czyli wylosować np. 100 różnych wartości parametru Qmeu, / = 1, 2, ...,100 (spełniających narzucone ograniczenia), co oznacza, że założymy liczbę osobników w populacji równą 100. Dla każdej wartości należy rozwiązać zadanie bezpośrednie i wyznaczyć wartość funkcji przystosowania (2.17), innymi słowy ocenić przystosowanie osobnika Q„eti- Najlepszym osobnikiem w populacji startowej jest ten, dla którego wartość funkcji (2.17) jest najmniejsza.

Następną populację chromosomów tworzymy stosując selekcję oraz operatory ewolucyjne. Oczywiście w rozważanym zadaniu, gdy chromosom składa się tylko z jednego genu krzyżowanie proste nie znajduje zastosowania. Załóżmy, że wprowadzimy mutację równomierną z prawdopodobieństwem 0.2. Dla każdego genu (chromosomu) losujemy liczbę z przedziahi (O, 1). Jeżeli wylosowana liczba jest mniejsza lub równa 0.2, wówczas zmieniamy wartość genu na nowo wylosowaną wartość spełniającą zadane ograniczenia. Jeśli wprowadzimy operator klonowania, to po wylosowaniu liczby z przedziału (O, 1), przy założeniu, że prawdopodobieństwo klonowania wynosi 0.3, naj­

lepszy osobnik (o najmniejszej wartości funkcji przystosowania) zostanie wprowadzony do nowej populacji pod warunkiem, że wartość wylosowanej liczby jest mniejsza lub równa 0.3.

Do stworzenia nowej populacji możemy zastosować np. selekcję turniejową.

Wybieramy w sposób losowy 4 osobników, a następnie na podstawie wartości funkcji przystosowania każdego z nich selekcjonujemy najlepszego osobnika do następnego pokolenia. W rozważanym zadaniu proces ten należy powtórzyć 100 razy po to, aby w kolejnym pokoleniu zachować liczebność populacji.

Opisany proces iteracyjny powtarzamy dla nowej populacji, tworząc w ten sposób kolejną. Jeśli spełniony zostanie warunek zatrzymania obliczeń, najlepszy chromosom końcowej populacji {Q„et) stanowi rozwiązanie omawianego zadania odwrotnego.

W przypadku, gdy warunkiem zatrzymania obliczeń jest z góry założona liczba pokoleń, wówczas najlepszy chromosom nie musi należeć do ostatniej populacji.

2. Zadania odwrotne w przepływie biociepła - przegląd literatury

Niewątpliwą zaletą algorytmów ewolucyjnych jest równoczesne przeszukiwanie całej dziedziny rozwiązań oraz duże prawdopodobieństwo osiągnięcia minimum globalnego.

Wadą natomiast - długi czas obliczeń związany między iimymi z wielokrotnym rozwiązywaniem zadania bezpośredniego. Ponadto, zastosowanie metod ewolucyjnych w wielu zadaniach (zwłaszcza 3D) jest znacznie prostsze niż wykorzystanie metod gradientowych.

2.3.3. Metody hybrydowe

Opisane powyżej metody identyfikacji posiadają pewne ograniczenia i wady, które mogą się nasilać zwłaszcza wraz ze wzrostem poziomu trudności zadania. Metody bazujące na znajomości gradientu mogą zmierzać do optimów lokalnych, trudniejsze może okazać się także wyznaczenie gradientu. Metody identyfikacji ewolucyjnej oka2oiją się natomiast bardzo czasochłonne [69, 70, 71, 72, 73], Alternatywą mogą być metody hybrydowe rozumiane jako algorytmy stanowiące różnego rodzaju połączenia między powyższymi metodami.

Celem połączenia obydwu algorytmów było stworzenie metody, która łączyłaby zalety obydwu algorytmów. Powinna ona zapewnić globalny charakter poszukiwań, co cechuje metody ewolucyjne, oraz powiima dość szybko i precyzyjnie dążyć do optimum, co ma miejsce w metodach gradientowych.

W niniejszej pracy doktorskiej zastosowano tzw. model strategii dwufazowej [86] - rys. 2.6

Rys. 2.6. Model strategii dwufazowej

Model strategii dwufazowej jest najprostszym przypadkiem metod hybrydowych.

Algorytm ewolucyjny stanowiący etap pierwszy ma za zadanie wygenerować „mocny”

punkt w pobliżu optimum globalnego, następnie algorytm gradientowy, stanowiący etap drugi, w kolejnych krokach szybko „prowadzi” do dokładnego rozwiązania. Jeżeli algorytm

2. Zadania odwrotne przepływie biociepła - przegląd Uteratury

ewolucyjny jest „odporny” na wpadanie w minima lokalne, to jest w stanie wygenerować punkt w pobliżu optimum globalnego. Algorytm gradientowy wykorzystuje swą zaletę szybkiego dojścia do optimum.

Problemom zastosowania algorytmów gradientowych, ewolucyjnych oraz metod hybrydowych do identyfikacji parametrów termicznych i geometrycznych podobszaru nowotworowego tkanki poświęcono rozdział 5 pracy doktorskiej.

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

ROZDZIAŁ 3

Metoda elementów brzegowych z zastosowa­

niem wielokrotnej zasady wzajemności

3.1. Wstęp

Metoda elementów brzegowych, obok metody elementów skończonych i metody różnic skończonych, jest coraz powszechniej stosowaną metodą rozwiązywania zadań brzegowych i brzegowo-początkowych [4, 5, 9, 10, 11, 12, 14, 15, 34, 45, 57]. Do najważniejszych jej zalet należy przede wszystkim duża dokładność aproksymacji warunków brzegowych, możliwość dokładnego odtworzenia rzeczywistej geometrii obszaru oraz stosunkowo niewielka liczba niewiadomych w tcw. układzie rozwiązującym (końcowym układzie równań), które są związane jedynie z brzegiem obszaru. W wielu zagadnieniach, np. przy wyznaczaniu bezźródłowych pól temperatury, dyskretyzacji podlega jedynie brzeg rozpatrywanego obszaru i dlatego często mówi się o tej metodzie, że zmniejsza wymiar zagadnienia o 1. Istnieje jednak duża grupa problemów brzegowych i brzegowo- początkowych, rozwiązanie któiych za pomocą klasycznego algorytmu MEB wiąże się z koniecznością dyskretyzacji zarówno brzegu jak i wnętrza obszaru. Aby „ocalić”

najcenniejszą zaletę MEB, czyli ograniczyć się jedynie do dyskretyzacji brzegu, rozwijane są takie jej warianty które, ogólnie rzecz biorąc, całkę po obszarze sprowadzają do całki brzegowej [13, 43, 80, 81, 82, 83, 84, 85]. Jedną z takich odmian jest MEB z zastosowaniem wielokrotnej zasady wzajemności, która zostanie omówiona w niniejszym rozdziale.

3.2. Metoda elementów brzegowych

Algorytm metody elementów brzegowych z zastosowaniem wielokrotnej zasady wzajemności przedstawimy dla następującego równania

X S / ^ T { x ) - k T { x ) ^ Q = ^ (3.1) gdzie X [W/mK] jest współczynnikiem przewodzenia ciepła, T oznacza temperaturę, k [W/m^K] jest stałym współczynnikiem, a Q [W/m^] stałym składnikiem źródłowym.

3. Metoda elementów brzegoiyych z zastosowaniem wielokrotnej zasady wzajemności

X - (jci, ^2) dla zadania płaskiego, x = (ati, X2, x%) dla zadania przestrzennego, natomiast

= (3-2)

i-i CWj gdzie w jest wymiarem zagadnienia.

Z matematycznego punktu widzenia, równanie (3.1) jest równaniem Poissona, w którym składnik źródłowy zależy od temperatury.

Porównując równanie (3.1) z równaniem Pennesa (wzór (2.1)) można zauważyć, że k^G^Ca, natomiast Q = kTg+Q^^i. Tak więc, równanie Pennesa jest szczególnym przypadkiem równania Poissona.

Równanie (3.1) uzupełniają następujące warunki brzegowe [36, 37]

* s r , : r ( * ) = r .

(3.3)

Y,-. q { x ) = - X ^ ^ = a [ T { x ) - r on

gdzie r = r , u T j u T j jest brzegiem rozpatrywanego obszaru Q, 7^ - znaną temperaturą na fragmencie brzegu T ,, qb - znanym strumieniem ciepła na fragmencie brzegu T2, a [W/m^K] - współczynnikiem wymiany ciepła, T “ - temperaturą otoczenia, dTIdn - pocłiodną w kierunku normalnym do brzegu

d T ( x ) ^ d T ( x )

= 2^ . - cosa. (3.4)

dn dx^

przy czym cos«^ sącosinusami kierunkowymi wersora normalnego do brzegu, zorientowa­

nego na zewnątrz obszaru.

W klasycznym wariancie metody elementów brzegowych [11, 14, 15, 45] równanie całkowe będące odpowiednikiem równania (3.1) jest następujące

B ( ^ ) T ( ^ ) + ¡ q ( x ) r ( ^ , ^ )d r =

(3-5) r ( x ) q ‘

(^, ;c)dr + J[-^r(:c) + e ] r (^, ;c)dQ

r n

gdzie ^ jest punktem obserwacji (punktem, w którym przyłożono skupione źródło ciepła), dla ^ G r : e (0 ,l) jest współczynnikiem zależnym od kształtu brzegu w pobliżu punktu dla 4 ^ Q : 5 ( ^ ) = 1, z kolei jest rozwiązaniem fundamentalnym (podstawowym), natomiast

(3.6)

Znajomość rozwiązania fundamentalnego jest niezbędna przy konstrukcji algorytmu MEB i dla omawianego zadania funkcja T* { i , x ) musi być rozwiązaniem równania

XV^7’- (ą ,x ) = -8(ą,;c) (3.7)

gdzie 8 {^ ,x) jest funkcją delta Diraca.

Postulat ten spełnia następująca funkcja [10, 13, 14, 44]

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

^ In -, zadanie 2D

(3.8) zadanie 3D

^%kr

gdzie r jest odległością między punktem obserwacji ^ a rozpatrywanym punktem x

e=i

Strumień ciepła (por. wzór (3.6)) wynikający z rozwiązania fundamentalnego można obliczyć analitycznie

gdzie

zadanie 2D

(3.10) zadanie 3D

47cr'

i^ = Z ( ^ e - ^ . ) c o s a , (3.11)

e=l

Jak można zauważyć, w brzegowym równaniu całkowym (3.5) występuje całka związana z wnętrzem obszaru co oznacza, że dla źródłowych pól temperatury, na etapie obliczeń numerycznych należy dyskretyzować zarówno brzeg F jak i wnętrze Q.

3.3. MEB z zastosowaniem wielokrotnej zasady wzajemności

w metodzie elementów brzegowych z zastosowaniem wielokrotnej zasady wzajemności występującą po prawej stronie równania (3.5) całkę po wnętrzu obszaru zastępuje się sumą całek po jego brzegu [80, 81, 82, 83, 84, 85].

Sposób postępowania jest następujący [26, 46, 47]. Oznaczymy przez / całkę po rozpatrywanym obszarze Q

I = \ [ Q ~ k { x ) ] r { ^ , x ) ń ę i (3.12)

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

i zapiszemy ją jako

I = \ [ Q - k T { x ) Y v : { ^ ^

czyli zakładamy znajomość funkcji takiej, że

r { ^ , x ) = v X { ^ , x ) Stosując drugą fonnułę Greena [45, 90] otrzymujemy

I = ^ V ^ [ Q - k T { x ) \ v ' { ^ , x ) d Q +

dr

lub

I = - k \ v ^ T { x ) v ; { ^ , x ) d n +

^ d v : { ą , x ) , . , J T { x ) [ Q - k T { x ) Y - ^ ^ - V : { ^ , x ) k dn

dr

Po wprowadzeniu zależności

oraz wyznaczeniu z równania (3.1)

mamy

' Q - k T { x ) ] v ; { i „ x ) d a -

I i[[■^^ W - 2 ] 2 ,'( i x ) - V ; { ^ , x )A g'(jc )]d r

(3.13)

(3.14)

(3.15)

(3.16)

(3.17)

(3.18)

(3.19)

czyli

/ z z _ Q - k T { x ) ] v ; { i „ x ) d Q-

^ \ A { ^ ^ x ) T { x ) d Y - ^ \ z \ { ^ , x ) ń Y - ^ \ v ; ( ą , x ) q ( x

)dP

W podobny sposób przekształcamy pierwszą z całek występującą po prawej stronie zależności (3.20). Tak więc

/ , = ^ i [ Q - k T ( x ) ] v ; ( ą , x ) d n = ^ j [ Q - k T ( x ) ] v % ' ( ą , x ) d n (3.21)

gdzie

v ; ( ą , x ) = v % ' ( ą , x ) (3.22)

Z drugiej formuły Greena otrzymujemy

/ , = 7 f v ^ [ e - A : r ( ; c ) ] K ; ( ą , x ) d Q +

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

czyli

dr

/,= -^ |v^ r(^ )F;(ą,^ )dQ +

^ n

^ J z; (ą, x)r(x)dr - e A I z; (ą, ^)dr - ^ J v; (^, (x)dr

(3.23)

(3.24)

przy czym

on Wprowadzając wzór (3.24) do (3.20) mamy

I = ^ \ [ Q - k T { x ) \ v ; { i „ x ) A ę i - ^ \ z ; { ^ , x ) ń T -

A, Q A p

k \ ^ f „ . /„ s„/ \

J z; (ą, x ) d r J z; (ą, x ) r ( x ) d r | z; (ą, x)r(x)dr -

A P A, j.

(3.25)

(3.26)

- r j K (^> W dr J F; (ą, (;c)dr

A ^P A ^P

Podobnie jak poprzednio, należy przekształcić pierwszą całkę po prawej stronie zależności (3.26). Postępując w ten sposób nieskończenie wiele razy uzyskuje się

' k ' '

Jz;( 4 ,i)dr+

^ /-I \ ^J r fk^

/»I p ; „ i V^A.y p

Po wykorzystaniu wzoru Q.27), równanie (3.5) przyjmuje postać

S ( y T-(ę) + Jq (x )T - (?, * ) d r + ¿ i i i v; (ą, *)dF =

r /=! \^/ r

r ( x ) ^ ' ( ą , ^ ) d r + £ ^ ¡ T { x ) z ; { s „ x ) d T - /-i V Ay p

A /-i lub po wprowadzeniu dodatkowych oznaczeń

fk^l-l u > r

z;{ ^ ,x )d T

(3.27)

(3.28)

3. Metoda elementów brzegowych z zasiosowaniem wielokrotnej zasady wzajemności

¡

U ; _V

/-O 'k'^/

r

T ( x ) z : ( % . x ) i T - i ±

^ /-I k

(3.29) z ; { ą , x ) d r

gdzie

oraz

orí Funiccje V,' (^, ;c) są rozwią^niami następujących równań

V X i ^ , x ) = T ' { ^ , x )

v X A ą , x ) = v ; { ą , x ) , 1=1,2,...

natomiast

C77

(3.30)

(3.31)

(3.32)

(3.33)

Oczywiście, otrzymanie efektywnego algorytmu jest możliwe jedynie pod warunkiem, że znana jest postać r o z w i^ ń F/ (^, ;c).

W przypadku zadań 2D są to funkcje [26, 80, 81, 82]

F / (ą ,x ) = : ^ r = ' U l n - + Ą 27U v r gdzie

Bo = 0, B ,=

_ A-\

4/

1

AV / + Ą -.

, / = 0,1,2,...

/ = 1,2,3,...

, / = 1,2,3,...

(3.34)

(3.35)

Odpowiedniki strumieni ciepła dla tych rozwiązań wyznacza się analitycznie i są one równe

(3.36)

d 2 1 - 2

z ; ( 4 - ) = - ^ r 1 '' A , - 21 A¡\n— + B¡

\ r

W zadaniach przestrzennych wykorzystujemy zależność [26, 84]

gdzie

v ; { % , x ) = - ^ r ‘- ' c „ / = 0,1,2,...

r = \ , C, = - , C, = —

° 1 2 ’ " 2 4

(3.37)

C , = 1

(2 / - l)(2 / - 3 ) / = 3,4,5,...

(3.38)

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

oraz

z ; (ą, x ) = - — {21 - \ y - % , I = 0,1,2,...

471 (3.39)

Dodatkowym problemem związanym z wykorzystaniem równania (3.28) jest zbieżność występujących w nim szeregów oraz oszacowanie popełnianych błędów przy zastąpieniu ich sumami skończonymi, co jest niezbędne w realizacji komputerowej metody.

Stosując kryterium d’ Alamberta zbieżności szeregów otrzymujemy

UJ fk'^

Uy

<1 (3.40)

oraz

r/tV UJ

fk

U . z;{ą,x) /

<1 (3.41)

czyli dla zadania 3D (por. wzory (3.37), (3.39)):

—r .2

X C,

<1

* , ( 2 M S

—r /+!

(3.42)

<1

skąd

r <

W

X (2 / - l)C , k { 2 l + \)C,7+1

Jak łatwo sprawdzić, dla / = 1 składnik po prawej stronie jest najmniejszy, więc

r < 2 l

(3.43)

(3.44)

Warunek zbieżności (3.44) ogranicza wymiaiy obszaru, dla którego można stosować ten wariant metody elementów brzegowych. W przypadku „za dużych obszarów” zaleca się podział rozpatrywanego obszaru Q na podobszary [81].

Dla zadań 2D (por. wzory (3.34), (3.36)) kryterium d’Alamberta prowadzi do zależności k^.2 - Am +

X - A, Inr + B, <1

X A i - 2 l { - Ą \ n r + B,)

(3.45)

<1

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

czyli

Ą+i

\nr A , A I n r - ^

4 l + (2/ + 2

A 1 + 2/

< — X k

\ n r - ^

\ n r - ^ A )

(3.46)

< —

3.4. Model numeryczny

Jeżeli szeregi występujące w zależności (3.29) są zbieżne, to równanie (3.29) można zapisać w postaci

■k"

« ( 4 ) n o + j ? w Z 7

r /»oV^y .V

i ' W I z ; { ^ , x ) d r

p /«oV^y

(3.47)

Innymi słowy, zamiast sumy całek brzegowycłi wygodniej jest rozpatrywać całki z sum odpowiednich szeregów.

W celu rozwiązania równania całkowego (3.47) brzeg obszaru dzielimy na N elementów brzegowych i całki występujące w zależności (3.47) zastępujemy sumą całek po tych elementach. Tak więc

N 00

£ J

f I (4'. - f i J ¿ i i ] ' " z ;(4'. .)dr,

j - 1 J- t = 0 \ ^ J A y - l r t - t \ ^ J

(3.48)

3.4.1. Elementy stałe

x e T . :

Dla stałych elementów brzegowych zakładamy, że

r W = r ( y ) = r ,

gdzie jest punktem centralnym elementu brzegowego F j . Równania (3.48) dla i = 1,2, przyjmują wówczas postać

(3.49)

3. Metoda ele/nentów brzegowych z zastosowaniem wielokrotnej zasady wzajemności

y-i r,. I-o

JM r,. /“O

Wprowadzamy następujące oznaczenia

Ą = J S ( T k / ( i ' . ^ K /=oV^y

oraz

X

f

• I

K

natomiast

A. J.i J- /.I \ KJ Uiiład równań (3.50) można zapisać jako

lub

gdzie

Z g a = Z " / . + z,

;-i ;-i

(3.50)

(3.51)

(3.52)

(3.53)

(3.54)

(3.55)

Często stosuje się również postać macierzową tego układu, a mianowicie G q - H T + Z

Uwzględniając zadane warunki brzegowe układ ten sprowadzamy do postaci A Y = jego rozwiązaniu temperatury w węzłach wewnętrznych obliczamy z zależności

>=i j.\

(3.56)

(3.57) : F , a po

(3.58)

3.4.2. Elementy liniowe - zadanie 2D

Dla brzegowych elementów liniowych (zarówno w sensie aproksymacji funkcji na elemencie, jak i kształtu elementu będącego odcinkiem) zakłada się liniową zależność temperatury i strumienia ciepła w obrębie elementu, czyli [12,45]

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

T{^) = N Ą x ’’ y N Ą x ^ )

q{T\) = N pq[ x' ’ ) + N , q [ x ^ )

(3.59)

gdzie:

_ 1- T1 1 + 11 N = ---Ł, N , =

'’ 2 * 2 ^(3.60)

przy czym ti6 [-1 , l], natomiast są współrzędnymi początku i końca liniowego elementu Ty - rys. 3.1.

Rys. 3.1. Liniowy element brzegowy Ty Całki występujące w zależności (3.48) przyjmą wówczas postać

oraz

'k - ''

?WZ

(3.61)

(3.62)

gdzie

oraz

;-o

(3.63)

(3.64)

(3.65)

(3.66)

3. Metoda elementów brzego-wych z zastosowaniem wielokrotnej zasady wzajemności

O l - ' i

1 /=0

/■ \/“l k _'i />iv^y

- ' ^"

’ -I /-I

Wykorzystano tutaj równanie parametryczne odcinka linii prostej

(3.68)

(3.69)

^ = (^ i.^ 2 )e r ^ :

przy czym (jak łatwo sprawdzić)

' 2 ^‘ 2 ^'

2 2

(3.70)

dr,=

j dii=

n /'l V L

U + ²

{ 2 J dri=-^dii (3.71)

gdzie /j jest długością elementu brzegowego.

Kolejnym etapem algorytmu jest utworzenie układu rozwiązującego, w którym niewia­

dome związane są z węzłami brzegowymi. Należy w tym miejscu podkreślić, że w węzłach brzegowych obszaru (np. narożach), w których występuje zmiana warunku brzegowego wprowadza się tzw. węzły podwójne, przy czym jeden z nich traktuje się jako koniec elementu Ty, natomiast drugi - początek elementu Ty+i (lys. 3.2).

© © © ©

Rys. 3.2. Węzeł pojedynczy i podwójny - elementy liniowe

Przyjmując numerację węzłów brzegowych r=\,2, układ równań (3.48) zapisuje się następująco

(3.72)

r « l r-1

przy czym dla węzłów pojedynczych

G , = G 5 + G ^ , ,

(3.73)

oraz dla węzłów podwójnych

Wzory (3.73) i (3.74) pokazują więc sposób „zszywania” elementów brzegowych Tj i Tj+i w węźle r. Należy przy tym podkreślić, że węźle brzegowym r=\ „zszywamy” element Fi z elementem F;v.

Układ równań (3.72) można również przedstawić w postaci

/ = 1,2...R (3.75)

________________ J. Metoda elementów brzegmyych z zastosowaniem wielokrotnej zasady wzajemności_______________

r=\

gdzie

przy czym [45]

i ^ r

H , - B „ i = r (3.76)

r - l r*i

oraz

(3.78)

A. r= l

Po rozwiązaniu układu (3.75) znane są wartości temperatur i strumieni ciepła we wszystkich węzłach brzegowych. Temperatury w dowolnych punktach należących do wnętrza obszaru wyznacza się na podstawie równania

H R

(3.79)

r « l r« l

W powyższym wzorze daszek nad literą/i został opuszczony zgodnie z (3.76), ponieważ dla węzłów wewnętrznych i^r.

3.5. MEB z zastosowaniem wielokrotnej zasady wzajemności - zadanie dodatkowe

w dalszych rozdziałach niniejszej pracy doktorskiej oprócz problemu podstawowego opisanego równaniem (3.1) z warunkami (3.3) rozwiązano również zadania dodatkowe związane z analizą wrażliwości.

Ogólnie rzecz biorąc, zadanie dodatkowe można przedstawić w postaci następującego równania [26,46,47]

x e C l : X V - U { x ) - k U { x ) + w T { x ) + b = 0 (3.80) gdzie U jest funkcją wrażliwości, T - temperaturą oraz w i ^ są stałymi. Równanie (3.80) uzupełnia się odpowiednimi warunkami brzegowymi.

Zastosowanie klasycznego podejścia MEB opisanego w podrozdziale 3.2 prowadzi do następującego równania (por. wzór (3.5))

B { t ) U { £ ) ^ \ w { x ) r { ^ ;c)dr =

r r r ^ . (3- S l )

f / (x )i* (£ , jr)d T + - k U { x ) + w T { x ) + b 7’*(£, .r)dQ

r Q

z rozwiązaniem podstawowym T ‘ { t, x ) zdefiniowanym zależnością (3.8), funkcją q ' { ę , x ) określoną wzorem (3.10) oraz flinkcją W { x ) = - X d U [ x ) l d n .

Zajmiemy się całką po obszarze Q, obszarze mianowicie

/= J[-A:i7(x) + wr(.r) + 6]r*(£,.T)di2 (3.82)

________________ i. Metoda eleTnentmy brzegowych : wtosowamem wielokrotne/ zasady wzcgemności_____________

Ponieważ znana jest funkcja V ' (£, x ) spełniająca zależność (3.14), więc

/ = \ \ - k U { x ) + w T { x ) ^ b ' \ v - v ; { ^ x ) d n n

Wykorzystujemy drugą formułę Greena

/ = J V ‘ - k U { x ) + w T [ x ) + b F;'{£, x )d O +

- k U { x ) + wT [ x ) + ¿>1- - V ' (ę, x ) — [~ k U ( x ) + w f ( x ) +

on ^

(3.83)

dn

(3.84) dr

czvli

1= - k V - U { x ) + w V ' T { x ) ] V; {Ł, x ) d n 4

fi [-k u {x )+ wr (x)+¿1.^^ v; (ć, x)

dn dn

(3.85)

dr

Wprowadzamy zależności (3.17), (3.18) oraz z równania (3.80) wyznaczamy

(3.86)

I wówczas

' - i

J.J k U [ x ) - w T [ x ) - b Z,’ x) - (ę, x ) łfF (x )-H - ę (x )'| dT

(3.87)

lub

3. Metoda elementów brzegowych z zaslosowaniem wielokrotnej zasady wzajemności

1 = a

v; {^, x ) ó n + k

X^ U { x ) z : { ^ , x ) A Y - ^ \ T { x ) z : { ^ , x ) d T - \ b \ z : { ^ , x ) ń T -

A. P A p (3.88)

A. p

W podobny sposób przekształcamy pierwszą z całek występującą po prawej stronie zależności (3.88), czyli

(3.89)

Ponieważ znamy funkcję F j(^ , jc) spełniającą równanie (3.22), więc

A n

v X ( ą , x ) d n (3.90)

i z drugiej formuły Greena otrzymujemy

/.=

n

- Ł v‘ U (x) + 2 ^ V ‘T (x)

j : L u ( x ) . 2 ^ T ( x ) Ą , - l a

dn (3.91)

e d u { x ) ^ ^^kwd T {x )

X dn X dn

dr

a po wykorzystaniu wzorów (3.18), (3.86), (3.25) mamy

/ .=

n L

- — U { x ) + ^ - ^ T { x ) + - ^ b - 2 — Q F ;(ą ,x )d Q +

^ j t / ( x ) z ; ( ą , x ) d r - 2^ J r ( x ) z ; ( ą , x ) d r - ^ z , J z ; ( ą , x ) d r + (3.92)

r A. p P

Podobnie jak poprzednio, należy przekształcić pierwszą całkę po prawej stronie zależności (3.92), czyli

- — U { x ) + 3— T { x ) + ^ b - l ^ Q

vv;(ą, x)dn =

- i l v > i ; w + 3 ^ v > r w v;{ą, ^ )d Q +

- - U { x ) + 2 ~ T { x ) + — b - 2 — Q 5 F ;(ą ,x ) (3.93)

k^ d U { x ) ^ ^k^wdT{x) dn

dr

dr-

3. Metoda elementów brzegowych zastosowaniem ■wielokrotnej sjsady w:sgemnaici

lub po zastosowaniu wzorów (3.18). (3.86) i (3.33) dla 1 = 3:

Q L

^ a ( . , ) + 4 ^ r ( . ) + ^ 4 - 3 ^ o l r , - ( Ę . , ) d Q +

— r

f/(.r)Ą(ć, j) d T - 3 ^ Jr(.t)z:(ę, .t)dT-|:jj 2 ;(i, .tjdTH.

A. [. A, j,

2 $ ^ e J z ; ( Ę , . t ) d r - t j r ; ( 5 , , : ) ( F ( x ) d r + 3 ^ J ( - ; { £ , . t ) , ( . v ) d r

A. r A- r r

Wprowadzając otrzj-mane zależności do wzoru (3.88) mamy

r =

l i J c ^ W z ; ( s , . ) d r - t f

IM V ^ / r

‘ ¿ I/«s ^ ' i N/-1

r z , - ( ć , j ) d r + o X ^ ( '- 0 /=*2 ^ u >

r;(ą,x)da-f

r ( . r ) z ; ( ć , - c ) d T -

z ; ( ć , . r ) d r -

3 3 / . nM

/«I

^

r

i. i

Postępując w ten sposób nieskończenie wiele razy uzyskuje się

' i ' '

^ = S t J

/-iV^y r ^

/ ^r \i^

\^J r

X ,

/»I ^

~ A.' ' ' V. > r

t-i

T { x ) z ; { ^ x ) ć r -

K a - t ) » ^ ( - T ) d r + t ^ / i ^ l f ^ '( ę , . r ) ę ( x ) d r

f i“i ^ V''‘/ f

Uwzględniając wzór (3.96), równanie (3.81) przyjmuje postać

r /»I f

fF(.t)i^'{ę.T)dr=

J

U(x)q

V ( Ł x ) r {. r ) ] d T

¿ , y ±I [ 2 ; ( ę , ^ ) d r ^ 0 y ™ ( / - i ) ' ’^ T ‘ f z ; { ę , . r ) d r

r X s, A. / j-

lub (por. podstawienia (330), (3.3i))

(3.94)

(3.95)

(3.96)

(3.97)

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

m m ( i ^ ) d r = ż i i T j t / w z ; (ą, . ) d r +

/-o r i ^ o \ X J p

x / - l

i £_{A, /«i p}J ' ^{i -}

u

(3.98) /~2

r

Jak można zauważyć, równanie dodatkowe (3.98) jest sprzężone z zadaniem podstawowym (3.1), (3.3), ponieważ wymaga ono znajomości zarówno funkcji T { x ) jak i funkcji q [ x ) na brzegu obszaru, które stanowią rozwiązanie problemu podstawowego.

3.6. Model numeryczny - zadanie dodatkowe

Równanie (3.98) zapiszemy w następującej postaci

w

1 r

z ; { ą , x ) d r -

l=> \'^y /»I \ X J

/ ' I \ l - 2

k

UJ

w X-

(3.99)

r 1 - 2 X )

3.6.1. Elementy stałe

w realizacji numerycznej, do dyskretyzacji brzegu obszaru zastosujemy stałe elementy brzegowe Y j , J = 1, 2, . . . , N( p o r . wzór (3.49)) i wówczas

lub

j-\ ^ r / - I

j7\X p /=iU> + Z ^ s J E ( / - i ) ( 7 ) z;{^ ,x)dTj

>-i r, /=2

U J

'>■

/ >

k

z c . » ; = z + z ( Ą i . - ) + K.

>1 y-l J»1 ’

gdzie G^j, H,j zdefiniowanowzorami(3.51), (3.56), natomiast

(3.101)