A
KADEMIAG
ÓRNICZO-H
UTNICZA IM. S
TANISŁAWAS
TASZICA WK
RAKOWIEWYDZIAŁ ELEKTROTECHNIKI, AUTOMATYKI, INFORMATYKI I IN ˙ZYNIERII BIOMEDYCZNEJ KATEDRA AUTOMATYKI I IN ˙ZYNIERII BIOMEDYCZNEJ
Rozprawa doktorska
Sterowanie procesów cieplnych z wykorzystaniem modeli
niecałkowitego rz ˛edu
mgr in˙z. Anna Obr ˛
aczka
Promotor: Prof. dr hab. in˙z. Wojciech Mitkowski
Praca powstała przy współudziale ´srodków pochodz ˛
acych z projektu
„Do-ctus - Małopolski fundusz stypendialny dla doktorantów“ współfinansowanego
ze ´srodków Unii Europejskiej w ramach Europejskiego Funduszu Społecznego.
Serdecznie dzi˛ekuj˛e promotorowi, Panu Prof. dr
hab. in˙z. Wojciechowi Mitkowskiemu, za cenne
uwagi i rady, a tak˙ze za wszelk ˛
a pomoc przy
po-wstawaniu tej pracy.
Dzi˛ekuj˛e równie˙z pracownikom Katedry
Automa-tyki i In˙zynierii Biomedycznej AGH za
meryto-ryczne dyskusje na temat pracy oraz pomoc.
Dzi˛ekuj˛e Rodzinie i Przyjaciołom za wsparcie i
cierpliwo´s´c.
Spis tre´sci
Spis rysunków ... 8 Spis tablic... 10 Wykaz oznacze´n... 11 1. Wst˛ep... 13 1.1. Cele pracy ... 14 1.2. Struktura pracy ... 15 2. Pochodne ułamkowe... 16 2.1. Definicja Grünwalda-Letnikowa ... 16 2.2. Definicja Riemanna-Liouville’a ... 17 2.3. Definicja Caputo ... 182.4. Zwi ˛azki mi˛edzy definicjami ... 18
2.5. Pochodna ułamkowa Riesza ... 19
2.6. Przykłady ... 19 2.6.1. Funkcja liniowa ... 19 2.6.2. Funkcja kwadratowa... 20 2.7. Transformata Laplace’a ... 22 2.7.1. Przykład ... 23 3. Przewodzenie ciepła ... 25
3.1. Ustalone przewodzenie ciepła ... 25
3.2. Nieustalone przewodzenie ciepła ... 26
3.3. Ułamkowe równanie przewodzenia ciepła ... 27
3.4. Współczynnik a ... 28
4. Obiekt eksperymentalny... 30
4.1. Opis eksperymentu laboratoryjnego ... 30
4.2. Modelowanie ... 32
SPIS TRE ´SCI 6 4.2.1. Pochodna Riesza... 33 4.2.2. Pochodna Caputo... 35 4.3. Identyfikacja parametryczna... 38 4.3.1. Równania ró˙zniczkowe... 38 4.3.2. Algorytmy identyfikacji... 40 4.3.3. Testy numeryczne ... 41 4.3.3.1 Test 1 . . . 41 4.3.3.2 Test 2 . . . 43 4.3.4. Wybór metody ... 43 4.4. Weryfikacja i porównanie ... 46 5. Porównanie modeli... 48 5.1. Model I ... 48 5.2. Model II ... 50 5.3. Model III... 50 5.4. Model IV ... 53 5.5. Model klasyczny... 55
5.6. Podsumowanie porównania modeli... 57
6. Model w dziedzinie Laplace’a... 59
7. Sterowanie... 63
7.1. Zamkni˛ety układ regulacji ... 63
7.2. Regulator PID ... 63
7.3. Regulator FOPID ... 65
7.3.1. Idealna posta´c Bodego transmitancji... 66
7.3.2. Projektowanie regulatora FOPID ... 67
7.4. Ocena regulatora... 68
7.5. Porównanie regulatorów dla modelu inercyjnego II rz˛edu z opó´znieniem ... 71
7.6. Porównanie regulatorów dla modelu ułamkowego z opó´znieniem ... 75
7.7. Wpływ parametrów na FOPID ... 80
8. Podsumowanie ... 82
Bibliografia... 84
Dodatek ... 90
SPIS TRE ´SCI 7
Abstract ... 94
Spis rysunków
2.1 Funkcja liniowa oraz jej pochodne: rz˛edu 1 oraz 0.5. . . 20
2.2 Funkcja kwadratowa oraz jej pochodne: rz˛edu 1 oraz 0.5. . . 21
2.3 Kolejne pochodne funkcji f (x) = x2, od 0 do 2 z krokiem 0.2. . . 22
2.4 Odpowied´z na skok jednostkowy układu o transmitancji s0.51+1 oraz 1 s+1. . . 24
4.1 Pojedyncze zdj˛ecie rozkładu temperatury wykonane kamer ˛a termowizyjn ˛a SC-660. . . 31
4.2 Sze´s´c zdj˛e´c rozkładu temperatury z ró˙znych chwil czasowych. . . 31
4.3 Schemat stanowiska laboratoryjnego. . . 32
4.4 Wykres przedstawiaj ˛acy zebrane pomiary. . . 32
4.5 Przykład RFDE z parametrami: β = 1.3, Kβ = 0.25, oraz warunkiem pocz ˛ at-kowym: g(x) = sin2(2x). . . 34
4.6 Porównanie rozwi ˛azania analitycznego i numerycznego metod ˛a SG równania RFDE z parametrami: β = 1.3, Kβ = 0.25, oraz warunkiem pocz ˛atkowym: g(x) = sin2(2x). . . 35
4.7 Przykład FDWE z parametrem: α = 1.7 oraz warunkiem pocz ˛atkowym: g(x) = sin(x). . . 36
4.8 Porównanie rozwi ˛azania analitycznego i numerycznego metod ˛a GMMP równa-nia FDWE z parametrem: α = 1.7 oraz warunkiem pocz ˛atkowym: g(x) = sin(x). 37 4.9 Przykład RFADE z parametrami: α = 1.8, Kα = 0.25, β = 0.2, Kβ = 0.25, oraz warunkiem pocz ˛atkowym: g(x) = x2(π − x). . . 39
4.10 Kryterium Q w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Gaussa-Newtona. . . 43
4.11 Czas oblicze´n w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Gaussa-Newtona. . . 44
4.12 Kryterium Q w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Levenberga-Marquardta. . . 44
SPIS RYSUNKÓW 9
4.13 Czas oblicze´n w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu
Levenberga-Marquardta. . . 45
4.14 Kryterium Q w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Sim-plex Neldera-Meada. . . 45
4.15 Czas oblicze´n w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Simplex Neldera-Meada. . . 46
5.1 Model I dla warto´sci parametrów: a = 2.28 · 10−7 oraz α = 1.48. . . 49
5.2 Wykres bł˛edu procentowego εndla ułamkowego modelu I. . . 49
5.3 Model II dla warto´sci parametrów: a = 9.61 · 10−8 oraz β = 1.00. . . 51
5.4 Wykres bł˛edu procentowego εndla ułamkowego modelu II. . . 51
5.5 Model III dla warto´sci parametrów: a = 6.57 · 10−9, α = 2.00 oraz β = 2.00 . 52 5.6 Wykres bł˛edu procentowego εndla ułamkowego modelu III. . . 53
5.7 Model IV dla warto´sci parametrów: a = 3.25 · 10−7, c = 2.89 · 10−3oraz α = 0.19 54 5.8 Wykres bł˛edu procentowego εndla ułamkowego modelu IV. . . 55
5.9 Model klasyczny dla warto´sci parametrów: a = 1.64 · 10−7. . . 56
5.10 Wykres bł˛edu procentowego εndla modelu klasycznego. . . 56
6.1 Wymuszenie, pomiary oraz odpowied´z modeli. . . 60
6.2 Pomiary oraz odpowied´z modeli. . . 61
6.3 Wykres bł˛edu procentowego εndla modelu GAP. . . 62
6.4 Wykres bł˛edu procentowego εndla modelu GBP. . . 62
7.1 Zamkni˛ety układ regulacji. . . 63
7.2 Odpowied´z na skok jednostkowy systemu z modelem A oraz regulatorem PID. 64 7.3 Odpowied´z na skok jednostkowy systemu z modelem B oraz regulatorem PID. . 65
7.4 Charakterystyka Bodego obiektu o transmitancji (7.5) dla γ ∈ (1, 2). . . 66
7.5 Odpowied´z na skok jednostkowy systemu z modelem A oraz regulatorem FO-PID rz˛edu γ = 0.8. . . 69
7.6 Odpowied´z na skok jednostkowy systemu z modelem B oraz regulatorem FO-PID rz˛edu γ = 0.8. . . 70
7.7 Graficzna interpretacja wska´zników: troraz em. . . 70
7.8 Graficzna interpretacja wska´zników: tn, emoraz e2. . . 71
7.9 Charakterystyka Bodego systemu z modelem A oraz regulatorem PID. . . 72
7.10 Charakterystyka Nyquista systemu z modelem A oraz regulatorem PID. . . 72
7.11 Charakterystyka Bodego systemu z modelem A oraz regulatorem FOPID. . . . 73
7.12 Charakterystyka Nyquista systemu z modelem A oraz regulatorem FOPID. . . . 74
7.13 Porównanie odpowiedzi systemu z modelem A oraz regulatorem PID i FOPID. 75 7.14 Porównanie sygnałów steruj ˛acych regulatorow PID i FOPID dla systemu z mo-delem A. . . 76
7.15 Charakterystyka Bodego systemu zmodelem B oraz regulatorem PID. . . 77
7.16 Charakterystyka Nyquista systemu z modelem B oraz regulatorem PID. . . 77
7.17 Charakterystyka Bodego systemu z modelem B oraz regulatorem FOPID. . . . 78
7.18 Charakterystyka Nyquista systemu z modelem B oraz regulatorem FOPID. . . . 78
7.19 Porównanie odpowiedzi systemu z regulatorem PID oraz FOPID. . . 79
7.20 Porównanie sygnałów steruj ˛acych regulatorow PID i FOPID dla systemu z mo-delem B. . . 80
7.21 Odpowiedzi systemu z regulatorem FOPID dla ró˙znych warto´sci γ. . . 81
7.22 Odpowiedzi systemu z regulatorem FOPID dla ró˙znych warto´sci ωc. . . 81
Spis tablic
3.1 Warto´sci dyfuzyjno´sci termicznej dla przykładowych materiałów . . . 294.1 Podsumowanie testów dla poziomu szumu δ = 0 - brak zakłóce´n. . . 42
4.2 Podsumowanie testów dla poziomu szumu δ = 0.01. . . 42
4.3 Podsumowanie testów dla poziomu szumu δ = 0.02. . . 42
4.4 Podsumowanie testu 2. . . 46
5.1 Porównanie modeli . . . 57
7.1 Porównanie wska´zników regulatorów dla systemu z modelem A . . . 73
Wykaz oznacze ´n
GLDα
a – pochodna w sensie Grünwalda-Letnikova rz˛edu α RLDα
a – pochodna w sensie Riemanna-Liouville’a rz˛edu α RLIα
a – całka w sensie Riemanna Liouville’a rz˛edu α CDα
a – pochodna w sensie Caputo rz˛edu α
dxe – cecha górna (funkcja sufit) x: dxe = min{k ∈ Z : k ≥ x}
n k – symbol Newtona n! – silnia [ ] – przedział domkni˛ety ( ) – przedział otwarty
Γ(x) – funkcja gamma (gamma Eulera)
Tn(f, x0) – rozwini˛ecie funkcji f w wielomian Taylora stopnia n w
punk-cie x0
L{f (t); s} – transformata Laplace’a funkcji f (t)
L−1{F (s); s} – odwrotna transformata Laplace’a funkcji F (s)
f (t) ∗ g(t) – splot funkcji f i g ρ – gesto´s´c [mkg3]
cp – ciepło wła´sciwe [kgKJ ]
λ – współczynnik przewodnictwa cieplnego materiału [mKW ] a – dyfuzyjno´s´c termiczna [ms2]
xj – j∆x
tm – m∆t
u(m)j – aproksymacja warto´sci dokładnej u(x, t) dla x = xj i t = tm
exp (x) – funkcja ekspotencjalna
εmax – maksymalny bł ˛ad procentowy
εn – bł ˛ad procentowy dla ka˙zdego punktu
M RE – ´sredni (procentowy) bł ˛ad wzgl˛edny (Mean Relative Error) G(s) – transmitancja
SPIS TABLIC 12
P ID – regulator proporcjonalno-całkowo-ró˙zniczkuj ˛acy
F OP ID – regulator proporcjonalno-całkowo-ró˙zniczkuj ˛acy ułamkowego rz˛edu
Am – zapas modułu (amplitudy)
Ωm – zapas fazy
ωg – pulsacja odci˛ecia modułu
ωp – pulsacja odci˛ecia fazy
|z| – moduł lizcby zespolonej z arg[z] – argument lizcby zespolonej z
tn – czas narastania odpowiedzi systemu na skok
tr – czas regulacji odpowiedzi systemu na skok
em – odchylenie maksymalne odpowiedzi systemu na skok
e2 – druga amplituda odpowiedzi systemu na skok
κ – przeregulowanie odpowiedzi systemu na skok J – całka z kwadratu uchybu regulacji
1. Wst˛ep
Idea rachunku ułamkowego, nazywanego te˙z rachunkiem niecałkowitego rz˛edu, jest równie stara jak rachunek całkowitego rz˛edu. Ju˙z w 1695 roku, kiedy Leibniz i Newton sformułowali standardowy rachunek ró˙zniczkowy, L’Hôspital napisał list do Leibniza, w którym rozwa˙zał po-chodn ˛a stopnia 12. Od tego czasu wielu wybitnych badaczy, takich jak: N. H. Abel, B. Riemann, P. S. Laplace, J. Liouville, J. Fourier czy F. Riesz, zajmowało si˛e tym zagadnieniem, kład ˛ac podstawy pod teori˛e ułamkowego rachunku ró˙zniczkowego. Przez około 300 lat, była to jednak gał ˛a´z rozwijana jako czysto teoretyczna dziedzina matematyki.
Zmieniło si˛e to dopiero w drugiej połowie XX wieku, kiedy znaleziono pewne zastosowa-nia pochodnych i całek niecałkowitego rz˛edu w wielu dziedzinach nauki np. ekonomii, biolo-gii, in˙zynierii czy fizyce [47]. Szczególnie w ci ˛agu ostatnich kilkunastu lat temat ten stał si˛e popularny, zarówno w Polsce jak i na ´swiecie, czego wynikiem s ˛a liczne publikacje w for-mie ksi ˛a˙zek oraz artykułów naukowych. Z tych pierwszym mo˙zna wymieni´c takie prace jak: [41, 51, 58], w których zawarto obszern ˛a teori˛e rachunku ułamkowego, [66, 68] - gdzie opisano metody numeryczne dla równa´n niecałkowitego rz˛edu, czy te˙z pozycje po´swi˛econe zastoso-waniom w praktyce teorii rachunku ułamkowego, np. [7, 21]. Z ciekawszych pozycji mo˙zna jeszcze wymieni´c na przykład: [27] - w której przedstawiono metody rozwi ˛azywania równa´n ró˙zniczkowych ułamkowego rz˛edu zarówno w sposób analityczny jak i numeryczny lub [69] - w której omówione s ˛a ułamkowe systemy o dynamice chaotycznej. Je˙zeli chodzi o artykuły, to mo˙zna zauwa˙zy´c, ˙ze pojawiaj ˛a si˛e one coraz liczniej. Zarówno w czasopismach o tematyce czysto matematycznej jak i tych z dziedzin systemów dynamicznych, modelowania kompute-rowego czy fizyki i chemii. Opisywane s ˛a w nich ró˙znego typu równania i ich rozwi ˛azania: np. [15, 22, 23, 32, 40], metody numeryczne: np. [12, 38, 42, 47, 50, 59], zagadnienia teorii stero-wanie: np. [2, 28] oraz coraz cz˛e´sciej liczne zastosowania: [24, 25, 26, 29, 54, 53, 48, 63, 70] - w automatyce, [14] - w robotyce, [16, 45] - w elektronice, [44] - elektrotechnice, [4, 20] - w medycynie, [9] - w chemii czy nawet [36] - w finansach.
Wszystkie te publikacje wskazuj ˛a na coraz powszechniejsze u˙zycie rachunku ró˙zniczko-wego ułamkoró˙zniczko-wego, zarówno do modelowania ró˙znych zjawisk fizycznych i procesów, jak do
1.1. Cele pracy 14
zagadnie´n teorii sterowania. Dlatego postanowiono w niniejszej pracy wykorzysta´c wła´sciwo´sci tego rachunku do zamodelowania procesu przepływu ciepła oraz do zaprojektowania ułamko-wego regulatora PID.
1.1. Cele pracy
Układy niecałkowitego rz˛edu znane s ˛a w matematyce od XVII wieku, jednak dopiero w ostatnich latach podj˛eto próby stosowania rachunku ró˙zniczkowo-całkowego ułamkowego rz˛edu do zagadnie´n klasycznej teorii sterowania. W literaturze specjalistycznej pojawiaj ˛a si˛e przesłanki, ˙ze rachunek tego typu mo˙ze by´c zastosowany do opisu przepływów cieplnych w ró˙znych o´srodkach [60]. Badania te nie s ˛a jednak zaawansowane na tyle aby umo˙zliwi´c ich praktyczne, przemysłowe wykorzystanie. Klasyczne równanie przepływu ciepła (równanie Fo-uriera), opisuje procesy wymiany ciepła w wystarczaj ˛aco dokładny sposób dla znacznej klasy problemów przemysłowych, nie uwzgl˛ednia jednak niektórych specyficznych zjawisk w nich wyst˛epuj ˛acych. W ramach doktoratu planowana jest budowa serii przykładowych modeli nie-całkowitego rz˛edu, jak naj´sci´slej opisuj ˛acych rzeczywiste zachowanie układu. Kolejnym eta-pem b˛edzie analiza porównawcza oraz weryfikacja zaproponowanych modeli. Zostan ˛a one tak˙ze porównane z klasycznym modelem przepływu ciepła. Mo˙zna zatem sformułowa´c pierw-szy zasadniczy cel pracy nast˛epuj ˛aco:
Cel 1: Znalezienie modelu niecałkowitego rz˛edu, lepiej opisuj ˛acego proces przepływu cie-pła oraz jego weryfikacja.
Regulatory proporcjonalno-całkowo-ró˙zniczkowe - PID, s ˛a cz˛esto wykorzystywane w za-stosowaniach przemysłowych, a badania naukowe poszukuj ˛ace nowych rozwi ˛aza´n i algoryt-mów opartych na PID s ˛a prowadzone w wielu placówkach na całym ´swiecie [6]. Wraz z roz-wojem rachunku ró˙zniczkowego ułamkowego oraz wzrostem liczby jego zastosowa´n w au-tomatyce, pojawiły si˛e liczne publikacje na temat ró˙znych regulatorów PID niecałkowitego rz˛edu. W ostatnich latach licznie pojawiaj ˛a si˛e artykuły o regulatorach PIλ - [34] oraz PIλDµ -[10, 39, 46, 48, 63]. Mo˙zna w nich znale´z´c ró˙zne podej´scia do doboru nastaw takich regulato-rów. Niektórzy autorzy postuluj ˛a tak˙ze, ˙ze tego typu regulatory w przyszło´sci zast ˛api ˛a klasyczne regulatory PID, ze wzgl˛edu na swoje lepsze wła´sciwo´sci [10]. W pracy planowane jest zapro-ponowanie regulatora typu PID niecałkowitego rz˛edu steruj ˛acego procesem przepływu ciepła. Policzony zostanie równie˙z zwykły regulator PID, w celu porównania i oceny obu rozwi ˛aza´n. Drugi zasadniczy cel pracy mo˙zna sformułowa´c nast˛epuj ˛aco:
Cel 2: Zaprojektowanie regulatora PID niecałkowitego rz˛edu dla procesu przepływu cie-pła oraz porównanie go z klasycznym regulatorem PID.
1.2. Struktura pracy 15
1.2. Struktura pracy
Praca składa si˛e z sze´sciu rozdziałów zasadniczych, wst˛epu oraz podsumowania. We wst˛e-pie opisano motywacj˛e bada´n, cele pracy oraz jej struktur˛e. Rozdział drugi, zatytułowany „Po-chodne ułamkowe“ zawiera podstawowe definicje oraz własno´sci z rachunku ró˙zniczkowego ułamkowego. Zamieszczono tam tak˙ze kilka prostych przykładów ilustruj ˛acych ró˙znice mi˛e-dzy pochodn ˛a klasyczn ˛a, a ułamkow ˛a. W nast˛epnym rozdziale - „Przewodzenie ciepła” wpro-wadzone wzory fizyczne opisuj ˛ace przepływ ciepła, na których oparte s ˛a pó´zniejsze modele. Czwarty rozdział opisuje metodyk˛e przeprowadzonego eksperymentu. Na pocz ˛atku opisany jest przeprowadzony eksperyment oraz metoda uzyskania z niego pomiarów. Nast˛epnie przed-stawiony jest schemat, według którego dalej b˛ed ˛a weryfikowane kolejne modele. Kolejnym rozdziałem jest „Porównanie modeli“ - wprowadzone s ˛a tam kolejno cztery modele niecałko-witego rz˛edu oraz dla porównania standardowy model przepływu ciepła. S ˛a one symulowane, identyfikowane oraz weryfikowane według metodyki opisanej w rozdziale czwartym. Na koniec wybrany zostaje najlepszy model i jest to realizacja pierwszego celu pracy. Rozdziała szósty za-wiera opis modelu w dziedzinie Lapalce’a - jest to tylko przygotowanie i wyja´snienie pewnych rzeczy przed kolejnym rozdziałem. W rozdziale siódmym przedstawiona jest metoda projekto-wania zaproponowanego regulatora niecałkowitego rz˛edu. Policzony zostaje tu tak˙ze regulator PID w celu pó´zniejszego porównania go z proponowanym. Nast˛epnie dokonana jest ocena i porównanie obu regulatorów. Przedstawiony jest wpływ pewnych parametrów na dynamik˛e re-gulatora niecałkowitego rz˛edu. Badania te przeprowadzone s ˛a na przykładzie dwóch modeli przepływu ciepła: inercyjnego drugiego rz˛edu z opó´znieniem oraz modelu ułamkowego rz˛edu z opó´znieniem.
2. Pochodne ułamkowe
Rozdział po´swi˛econy jest wprowadzeniu podstawowych poj˛e´c z rachunku ró˙zniczkowego niecałkowitego rz˛edu, od pochodnej do transformaty Laplace’a. Ma on tak˙ze pokaza´c przyj˛ete oznaczenia i symbole, które b˛ed ˛a u˙zywane pó´zniej w pracy, gdy˙z w literaturze nie s ˛a one jeszcze ujednoznacznione.
Pochodna rz˛edu całkowitego w danym punkcie zale˙zy od lokalnego zachowania si˛e funkcji w otoczeniu rozwa˙zanego punktu. Inaczej jest z pochodn ˛a rz˛edu ułamkowego - zale˙zy ona od zachowania funkcji w całym rozwa˙zanym przedziale. W szczególno´sci wi˛ec, jest ona w stanie uwzgl˛edni´c wpływ warunków brzegowych. Głównie ta cecha powoduje jej liczne zastosowania w modelowaniu zjawisk, gdzie wpływ warunków brzegowych jest istotny.
W kolejnych podrozdziałach zostan ˛a omówione podstawowe definicje pochodnych niecał-kowitego rz˛edu.
2.1. Definicja Grünwalda-Letnikowa
Definicja ta została wprowadzona przez Antona Karla Grünwalda (1838-1920) w roku 1867 w Pradze oraz równolegle przez Aleksieja Wasilewicza Letnikowa (1837-1888) w Moskwie w roku 1868. Jest ona uogólnieniem definicji klasycznej pochodnej funkcji.
Definicja 2.1.1 (Pochodna funkcji) Pochodn ˛a funkcji f (x) w punkcie x0 nazywamy granic˛e
(o ile istnieje):
df (x)
dx = lim∆x→0
f (x0+ ∆x) − f (x0)
∆x . (2.1)
Stosuj ˛ac wzór (2.1) dwa razy otrzymamy formuł˛e na pochodn ˛a drugiego rz˛edu:
d2f (x)
dx2 = lim∆x→0
f (x0+ ∆x) − 2f (x0) + f (x0− ∆x)
(∆x)2 , (2.2)
oraz n (n ∈ N) razy na pochodn ˛a n-tego rz˛edu: dnf (x) dxn = lim∆x→0 1 (∆x)n n X k=0 (−1)kn k f (x − k · ∆x). (2.3)
2.2. Definicja Riemanna-Liouville’a 17
Definicja pochodnej ułamkowej rz˛edu α (α ∈ R) jest uogólnieniem tego wła´snie wzoru. Definicja 2.1.2 (Pochodna w sensie Grünwalda-Letnikowa) Pochodn ˛a rz˛edu α (α ∈ R+),
na przedziale[a, b] (a 6 x 6 b), w sensie Grünwalda-Letnikowa nazywamy operator GLDα a [66]: GLDα af (x) = lim ∆x→0 1 (∆x)α x−a ∆x X k=0 (−1)kα k f (x − k · ∆x). (2.4)
Dodatkowo wprowadzamy oznaczenie:
(∆α∆xf )(x) = 1 (∆x)α x−a ∆x X k=0 (−1)kα k f (x − k · ∆x). (2.5)
Jest to ułamkowy odpowiednik ilorazu ró˙znicowego.
2.2. Definicja Riemanna-Liouville’a
W tym przypadku definicja pochodnej jest otrzymywana z definicji całki ułamkowej Riemanna-Liouville’a, która jest uogólnieniem wzoru na całk˛e n-tego rz˛edu (2.6), oraz z in-nego dobrze znain-nego w klasycznym rachunku ró˙zniczkowym wzoru (2.7).
Lemat 2.2.1 Funkcja f jest całkowalna w sensie Riemanna na przedziale [a, b]. Wtedy, dla ka˙zdegox ∈ [a, b] oraz n ∈ N zachodzi:
Ianf (x) = 1 (n − 1)!
Z x
a
(x − ξ)n−1f (ξ)dξ. (2.6)
Lemat 2.2.2 Funkcja f jest n razy ró˙zniczkowalna na przedziale [a, b] oraz m, n ∈ N i m > n. Wtedy zachodzi: dn dxnf = dm dxmI m−n a f. (2.7)
Uogólniaj ˛ac wzór (2.6) na rz˛edy niecałkowite, przy u˙zyciu funkcji Gamma (Γ) otrzymamy definicj˛e:
Definicja 2.2.3 (Całka Riemanna-Liouville’a) Całk ˛a Riemanna-Liouville’a rz˛edu α (α ∈ R+), zdefiniowan ˛a na przedziale[a, b], nazywamy operatorRLIaα [66]:
RLIα af (x) = 1 Γ(α) Z x a (x − ξ)α−1f (ξ)dξ. (2.8)
Korzystaj ˛ac z tej definicji i uogólniaj ˛ac wzór (2.7) otrzymujemy definicj˛e pochodnej ułamkowej w sensie Riemann-Liouville’a.
2.3. Definicja Caputo 18
Definicja 2.2.4 (Pochodna w sensie Riemann-Liouville’a) Pochodn ˛a rz˛eduα (α ∈ R+), na
przedziale[a, b] (a 6 x 6 b), w sensie Riemanna-Liouville’a, nazywamy operatorRLDα a [66]: RLDα af (x) = 1 Γ(n − α) d dx nZ x a (x − ξ)n−α−1f (ξ)dξ, n = dαe. (2.9)
2.3. Definicja Caputo
Wprowadzona przez Michaela Caputo w 1967 w pracy „Linear model of dissipation whose Q is almost frequency independent-II"(Geophys. J. R. Astr. Soc. 13: 529–539). Definicja ta tak˙ze wykorzystuje definicj˛e całki ułamkowej Riemanna-Liouville (2.8), ale w porównaniu z definicj ˛a pochodnej RL, odwrócona jest kolejno´s´c całkowania i ró˙zniczkowania: (2.10) zamiast (2.7), co wpływa na inn ˛a struktur˛e ko´ncow ˛a zdefiniowanej pochodnej.
Lemat 2.3.1 Funkcja f jest n razy ró˙zniczkowalna na przedziale [a, b] oraz m, n ∈ N i m > n. Wtedy zachodzi: dn dxnf = I m−n a dm dxmf. (2.10)
Definicja 2.3.2 (Pochodna w sensie Caputo) Pochodn ˛a rz˛eduα (α ∈ R+), na przedziale[a, b]
(a 6 x 6 b), w sensie Caputo, nazywamy operatorCDα a [66]: CDα af (x) = 1 Γ(n − α) Z x a (x − ξ)n−α−1 d dξ n f (ξ)dξ, n = dαe. (2.11)
2.4. Zwi ˛
azki mi˛edzy definicjami
Przedstawione definicje pochodnych ułamkowych, w ogólnym przypadku nie s ˛a sobie rów-nowa˙zne. Jednak przy pewnych zało˙zeniach zachodz ˛a mi˛edzy nimi pewne zwi ˛azki. Zostały one przedstawione poni˙zej. Pierwsz ˛a zale˙zno´sci ˛a jest zwi ˛azek mi˛edzy pochodn ˛a w sensie Caputo, a pochodn ˛a w sensie Riemanna-Liouville’a - opisuje j ˛a poni˙zsze twierdzenie.
Twierdzenie 2.4.1 Załó˙zmy, ˙ze α6 0 oraz n = dαe. Dodatkowo, przyjmijmy ˙zeRLDα
af istnieje
i jest(n − 1) razy ró˙zniczkowalna w a. Wtedy, prawie wsz˛edzie zachodzi [58]:
CDα af =
RL Dα
a(f − Tn−1(f, a)), (2.12)
gdzieTn−1(f, a) oznacza rozwini˛ecie funkcji f w wielomian Taylora stopnia n − 1, w punkcie
a.
Nast˛epne twierdzenie jest istotne do wyprowadzenie zwi ˛azku mi˛edzy definicj ˛a Grünwalda-Letnikowa, a pozostałymi dwoma.
2.5. Pochodna ułamkowa Riesza 19
Twierdzenie 2.4.2 Niech α 6 0, n = dαe oraz funkcja f jest n razy ró˙zniczkowalna na prze-dziale[a, b], wtedy [58]
GLDα af (x) = n−1 X k=0 f(k)(a)(x − a)k−α Γ(k + 1 − α) + 1 Γ(n − α) Z x a f(n)(ξ)dξ (x − ξ)1+α−n. (2.13)
Korzystaj ˛ac z powy˙zszego twierdzenia oraz zale˙zno´sci (2.12) mo˙zemy zapisa´c zwi ˛azek mi˛edzy wszystkimi trzema definicjami:
Twierdzenie 2.4.3 Niech α 6 0, n = dαe oraz funkcja f jest n razy ró˙zniczkowalna na prze-dziale[a, b], wtedy [58]
GLDα
af (x) = Tn−1(f, a)(x) +CDaαf (x) = RL Dα
af (x). (2.14)
2.5. Pochodna ułamkowa Riesza
W pracy u˙zywana b˛edzie pochodna w sensie Riesza. Jest to pochodna niecałkowitego rz˛edu, u˙zywana w równaniach z pochodnymi cz ˛astkowymi, do opisu pochodnej po zmiennej prze-strzennej. Pochodna w sensie Riesza na przedziale [a, b] (a6 x 6 b) funkcji dwóch zmiennych f (x, y) mo˙zna zdefiniowa´c nast˛epuj ˛aco [68]:
∂α ∂|x|αf (x, y) = − 1 2 cos (πα2 )( RL aD α x + RL x D α b)f (x, y), α 6= 1, (2.15) gdzie: RL aDαx = 1 Γ(n − α) ∂ ∂x nZ x a (x − ξ)n−α−1f (ξ, y)dξ, n = dαe (2.16)
jest pochodn ˛a lewostronn ˛a w sensie Riemanna-Liouville’a, natomiast:
RL xD α b = 1 Γ(n − α) ∂ ∂x nZ b x (x − ξ)n−α−1f (ξ, y)dξ, n = dαe (2.17)
jest pochodna prawostronn ˛a w sensie Riemanna-Liouville’a.
2.6. Przykłady
2.6.1. Funkcja liniowa
Rozwa˙zamy funkcj˛e liniow ˛a f (x) = x. Pierwsza pochodna:
d
dxf (x) = 1.
2.6. Przykłady 20
Nast˛epnie przyjmijmy α = 0.5 i policzmy pochodn ˛a funkcji f rz˛edu α na przedziale [0, 5]. Korzystaj ˛ac ze wzoru (2.11) otrzymujemy:
CDαf (x) = 1 Γ(0.5) Z x 0 (x − ξ)−0.5 d dξξdξ = 1 √ π Z x 0 1 √ x − ξdξ = 2√x √ π .
Funkcj˛e f oraz jej obliczone pochodne rz˛edu 1 i 0.5 przedstawiono na rysunku 2.1.
0 1 2 3 4 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x f(x) f(x) = x df(x) dx = 1 df0.5(x) dx0.5 = q 4x π
Rysunek 2.1. Funkcja liniowa oraz jej pochodne: rz˛edu 1 oraz 0.5.
2.6.2. Funkcja kwadratowa
Rozwa˙zamy funkcj˛e kwadratow ˛a f (x) = x2. Pierwsza pochodna:
d
dxf (x) = 2x.
Teraz obliczamy pochodn ˛a rz˛edu α = 0.5 na przedziale [0, 5]:
CDαf (x) = 1 Γ(0.5) Z x 0 (x − ξ)−0.5 d dξξ 2dξ = √1 π Z x 0 2ξ √ x − ξdξ = 8√x3 3√π.
2.6. Przykłady 21 0 1 2 3 4 5 0 5 10 15 20 25 x f(x) f(x) = x2 df(x) dx = 2x df0.5(x) dx0.5 = 83 q x3 π
Rysunek 2.2. Funkcja kwadratowa oraz jej pochodne: rz˛edu 1 oraz 0.5.
Ogólnie wzór na pochodn ˛a ułamkow ˛a funkcji f (x) = x2 (na [0, 5]) dla dowolnego α b˛edzie
wygl ˛adał nast˛epuj ˛aco:
CDαf (x) = 1 Γ(n − α) Z x 0 (x − ξ)n−α−1( d dξ) nξ2dξ, n = dαe.
Na rysunku 2.3 pokazano kolejne pochodne funkcji kwadratowej dla α od 0 do 2 z krokiem 0.2.
2.7. Transformata Laplace’a 22 0 1 2 3 4 5 0 5 10 15 20 25 x d α d x α f (x ) α = 0 α = 0.2 α = 0.4 α = 0.6 α = 0.8 α = 1 α = 1.2 α = 1.4 α = 1.6 α = 1.8 α = 2
Rysunek 2.3. Kolejne pochodne funkcji f (x) = x2, od 0 do 2 z krokiem 0.2.
2.7. Transformata Laplace’a
Podobnie, jak definicja samej pochodnej ułamkowej, tak i definicja transformaty Laplace’a, jest uogólnieniem poj˛e´c z rachunku klasycznego. Transformat˛e Laplace’a wyra˙za si˛e wzorem:
F (s) = L{f (t); s} = Z ∞
0
e−stf (t)dt, (2.18)
natomiast odwrotn ˛a transformat˛e wzorem:
f (t) = L−1{F (s); t} = 1 2πj
Z c+j∞
c−j∞
estF (s)ds, (2.19)
gdzie c jest stał ˛a, która jest wi˛eksza od cz˛e´sci rzeczywistych wszystkich punktów funkcji na płaszczy´znie s, w których funkcja F(s) nie istnieje. Korzystaj ˛ac teraz z zale˙zno´sci na transfor-mat˛e pochodnej n - tego rz˛edu:
L{fn(t); s} = snF (s) − n−1 X k=0 sn−k−1f(k)(0) = snF (s) − n−1 X k=0 skf(n−k−1)(0), (2.20)
oraz transformat˛e splotu funkcji:
2.7. Transformata Laplace’a 23 gdzie: f (t) ∗ g(t) = Z t 0 f (t − τ )g(τ )dτ = Z t 0 f (τ )g(t − τ )dτ (2.22)
mo˙zna uogólni´c wzór (2.18) dla pochodnych rz˛edu niecałkowitego w sensie Riemanna-Liouville’a: L{RLDα 0f (t); s} = s αF (s) − n−1 X k=0 sk[RLDα−k−1 0 f (t)]t=0, n − 1 ≤ α < n. (2.23)
Praktyczne zastosowanie tej transformaty jest ograniczone ze wzgl˛edu na brak fizycznej interpretacji warto´sci granicznych pochodnych ułamkowych dla dolnego ograniczenia: t = 0 [58].
Analogicznie mo˙zna otrzyma´c formuł˛e dla pochodnej ułamkowej w sensie Caputo:
L{CDα 0f (t); s} = s αF (s) − n−1 X k=0 sα−k−1f(k)(0), n − 1 ≤ α < n, (2.24)
której zastosowanie w praktyce jest du˙zo szersze, poniewa˙z wyst˛epuj ˛ace tutaj pochodne dla t = 0 maj ˛a fizyczne interpretacje, na przykład: f0(0) mo˙ze by´c pr˛edko´sci ˛a pocz ˛atkow ˛a, a f00(0) - pocz ˛atkowym przyspieszeniem [58].
2.7.1. Przykład
Rozwa˙zany jest układ rz˛edu niecałkowitego:
CDα
0x(t) = Ax(t) + Bu(t), y(t) = Cx(t), (2.25)
gdzie A jest macierz ˛a stanu, B macierz ˛a wej´scia, natomiast C macierz ˛a wyj´scia. Przyjmijmy teraz nast˛epuj ˛ace warto´sci: A = −1, B = 1, C = 1 oraz α = 0.5, nasz układ mo˙zna teraz zapisa´c:
CD0.5
0 x(t) = −x(t) + u(t), y(t) = x(t), (2.26)
Przechodz ˛ac w dziedzin˛e operatora Laplace’a otrzymamy:
s0.5X (s) = −X (s) + U (s), Y(s) = X (s), (2.27)
sk ˛ad nast˛epnie policzy´c mo˙zna transmitancj˛e rozwa˙zanego układu:
G(s) = Y(s) U (s) =
1
s0.5+ 1. (2.28)
Łatwo teraz otrzyma´c odpowied´z na skok jednostkowy: h(t) = 1(t) takiego układu:
H(s) = 1 sG(s) = 1 s 1 s0.5+ 1. (2.29)
Na rysunku 2.4 przedstawiono odpowied´z skokow ˛a rozwa˙zanego układu oraz, dla porównania, odpowied´z układu o transmitancji G1(s) = s+11 .
2.7. Transformata Laplace’a 24 0 5 10 15 20 25 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t h(t) 1 s0.5+1 1 s+1
Rysunek 2.4. Odpowied´z na skok jednostkowy układu o transmitancji s0.51+1 oraz
1 s+1.
3. Przewodzenie ciepła
W rozdziale tym przedstawiono krótkie wprowadzenie w tematyk˛e przewodzenia ciepła z punktu widzenia fizyki. Omówiono najwa˙zniejsze wzory dla poszczególnych przypadków oraz ich interpretacj˛e fizyczn ˛a. Rozdział ten nie ma na celu zagł˛ebiania si˛e w t ˛a tematyk˛e, a tylko krótkie przypomnienie najwa˙zniejszych poj˛e´c, które b˛ed ˛a wykorzystywane w dalszych cz˛e´sciach pracy.
Ciepło jest jedn ˛a z form energii, obok energii elektrycznej, mechanicznej, chemicznej czy j ˛adrowej. Energia jest miar ˛a ogólnie pojmowanego ruchu materii. Podej´scie energetyczne do ciepła zawdzi˛eczane jest badaniom Jamesa Joule’a z 1843 roku, w których pokazał on równo-wa˙zno´s´c energii cieplnej i mechanicznej. Mo˙zna zatem dla ciepła stosowa´c zasad˛e zachowania energii, która w tym przypadku wyra˙zona jest przez I zasad˛e termodynamiki. II zasada termo-dynamiki, mówi natomiast o przenoszeniu energii cieplnej, które jest wynikiem szczególnej własno´sci materii - temperatury [62].
Przenoszenie energii cieplnej nazywane jest równie˙z przepływem ciepła lub wymian ˛a ciepln ˛a i mo˙zna je podzieli´c na trzy podstawowe mechanizmy: przewodzenie, konwekcj˛e i promieniowanie cieplne [30]. Niniejsza praca skupia si˛e tylko na przewodzeniu ciepła, jako mechanizmie odgrywaj ˛acym zasadnicza rol˛e przy procesach cieplnych, dotycz ˛acych ciał sta-łych.
3.1. Ustalone przewodzenie ciepła
Przewodzenie ciepła wyst˛epuje mi˛edzy ciałami o ró˙znej temperaturze pozostaj ˛acymi ze sob ˛a w bezpo´srednim kontakcie. Mówimy o nim, gdy przekazywana jest energia bezładnego ruchu - poprzez drgania cz ˛astek, a nie w wyniku ich makroskopowego (uporz ˛adkowanego) ru-chu. Ten rodzaj ruchu jest charakterystyczny dla ciał stałych [62].
Przewodzenie ciepła w najbardziej ogólnym przypadku opisuje si˛e nieliniowym równaniem ró˙zniczkowym przewodnictwa cieplnego, zwanym równie˙z równaniem Fouriera-Kirchhoffa
3.2. Nieustalone przewodzenie ciepła 26 ((3.1)): 52T + 1 λq˙v = ρcp λ DT Dt , (3.1)
gdzie: T - temperatura, t - czas, λ - współczynnik przewodzenia ciepła, cp - ciepło wła´sciwe
(przy stałym ci´snieniu), ˙qv - g˛esto´s´c strumienia wewn˛etrznych ´zródeł ciepła, ρ - g˛esto´s´c,
nato-miast operatorDtD oznacza tutaj operator Stokesa (nazywany inaczej operatorem ró˙zniczkowania w˛edrownego lub pochodn ˛a substancjaln ˛a) i jest definiowany jako (3.2):
D Dt = ∂ ∂t + ∂x ∂t ∂ ∂x + ∂y ∂t ∂ ∂y + ∂z ∂t ∂ ∂z. (3.2)
To ogólne równanie Fouriera - Kirchhoffa redukuje si˛e w ró˙zny sposób przy odpowiednich zało˙zeniach [30]. Dla przypadku ustalonego przewodzenia ciepła równanie (3.1) sprowadza si˛e do równania ró˙zniczkowego Poissona (3.3):
52T = −1
λq˙v, (3.3)
gdy dodatkowo przyjmie si˛e, ˙ze to ustalone pole jest bez´zródłowe otrzymujemy jeszcze prostsze równanie ró˙zniczkowe Laplace’a (3.4):
52T = 0. (3.4)
W zastosowaniach najcz˛e´sciej wystarczy rozwa˙za´c zmiany temperatury tylko wzdłu˙z jednej zmiennej przestrzennej. W tej sytuacji i przy zało˙zeniach jak do równania (3.4), problem prze-pływu ciepła sprowadza si˛e do rozwa˙zania równa´n ró˙zniczkowych zwyczajnych, których roz-wi ˛azanie nie stanowi zwykle problemu [30]. Ciekawszym problemem jest przypadek nieusta-lonego pola temperatury.
3.2. Nieustalone przewodzenie ciepła
Z nieustalonym przewodzeniem ciepła mamy styczno´s´c bardzo cz˛esto - opisuje ono zmiany temperatury ´scian na skutek wahania temperatur na zewn ˛atrz lub wewn ˛atrz budynku. Innym przykładem mog ˛a by´c ró˙zne procesy technologiczne: nagrzewanie metali w piecach grzew-czych lub ochładzanie i nagrzewanie innych ciał [30].
W pracy tej podj˛eto prób˛e modelowania procesów cieplnych w ciałach stałych. W takim przypadku pochodna substancjalna upraszcza si˛e do pochodnej cz ˛astkowej po czasie, poniewa˙z nie ma przepływu czynnika. Dodatkowo, gdy zało˙zy si˛e brak wewn˛etrznych ´zródeł ciepła, to
˙
qv = 0 i równanie (3.1) przyjmuje posta´c (3.5):
52T = ρcp
λ ∂T
3.3. Ułamkowe równanie przewodzenia ciepła 27
Jest to tzw. równanie Fouriera i jest ono podstawowym modelem klasycznym przewodzenia cieplnego rozwa˙zanym w tej pracy. Jest to równanie ró˙zniczkowe cz ˛astkowe, które w zale˙z-no´sci od przyj˛etych warunków pocz ˛atkowych i brzegowych, mo˙ze przyjmowa´c ró˙zne stopnie skomplikowania. Rozwi ˛azania analityczne znane s ˛a dla pewnych przypadków szczególnych -dla jednorodnych o´srodków, stałych w czasie warunków pocz ˛atkowych czy podstawowych brył geometrycznych. Szczegółowe omówienie tych rozwi ˛aza´n mo˙zna znale´z´c m. in. w [30], [62] czy [56]. Stosowane s ˛a tak˙ze metody operatorowe, takie jak przekształcenie Laplace’a, Fo-uriera, Mellina i Hankela oraz funkcje Greena [30]. Wi˛ecej o tych metodach jest np. w [8], [49] oraz [61].
W praktyce rzadko jednak spotyka si˛e takie uproszczone przypadki, cz˛e´sciej wyst˛epuj ˛a zmienne warunki brzegowe, nieregularne geometryczne kształty czy zale˙zno´sci własno´sci fi-zycznych od temperatury [30]. Opracowane metody analityczne opisuj ˛a wtedy zjawisko prze-wodnictwa cieplnego z niewystarczaj ˛ac ˛a dokładno´sci ˛a lub te˙z w ogóle niemo˙zliwe jest ich u˙zy-cie. Stosowane s ˛a wtedy ró˙znego rodzaju metody numeryczne - zwłaszcza metoda elementów sko´nczonych (MES) - [71], metoda ró˙znic sko´nczonych (MRS) - [13] i metoda elementu brze-gowego (MEB) - [67]. Główn ˛a wad ˛a tych metod jest wra˙zliwo´s´c na dobór odpowiedniej siatki oraz długi czas oblicze´n.
3.3. Ułamkowe równanie przewodzenia ciepła
W pracach [19], [57] oraz [60] przedstawione zostały uogólnienia równania Fouriera (3.5). Mo˙zna z nich wyprowadzi´c ułamkowe równania ciepła w ró˙znych postaciach. Zostan ˛a one przedstawione poni˙zej i posłu˙z ˛a w tej pracy jako modele przewodzenia cieplnego niecałkowi-tego rz˛edu. Pierwsze z nich to czasowe ułamkowe równanie przewodnictwa cieplnego:
∂αT
∂tα = a4T, 0 < α < 2. (3.6)
Równanie (3.6) opisuje cały zakres, od lokalnego przewodzenia ciepła, poprzez standardowe przewodnictwo, a˙z do balistycznego przewodzenia ciepła [60]. Kolejne rozwa˙zane w tej pracy równanie to przestrzenne ułamkowe równanie przewodnictwa ciepła:
∂T ∂t = a
∂βT
∂|x|β, 0 < β 6 2, (3.7)
lub w wy˙zszych wymiarach: ∂T
∂t = −a(−4)
β
2T, 0 < β 6 2. (3.8)
Czasami, w oparciu o dodatkowe zało˙zenia fizyczne, w równaniach (3.7) i (3.8) rozwa˙za si˛e zakres 1 < β 6 2 [60]. Trzeci model niecałkowitego rz˛edu jest zbudowany na podstawie
3.4. Współczynnik a 28
ogólnego przestrzenno-czasowego, ułamkowego równania przewodnictwa cieplnego:
∂αT ∂tα = a ∂βT ∂|x|β, 0 < α 6 2, 0 < β 6 2, (3.9) lub ∂αT ∂tα = −a(−4) β 2T, 0 < α 6 2, 0 < β 6 2. (3.10)
Ułamkowe równania ciepła zostały wybrane do modelowania, poniewa˙z autorka uwa˙za, ˙ze szczególne własno´sci pochodnej rz˛edu niecałkowitego pozwol ˛a na uzyskanie wi˛ekszej szcze-gółowo´sci w opisie rozwa˙zanego zjawiska, jakim jest przewodnictwo cieplne.
3.4. Współczynnik a
Wyst˛epuj ˛acy w powy˙zszych równaniach współczynnik a nazywany jest w fizyce współ-czynnikiem wyrównania temperatur lub inaczej dyfuzyjno´sci ˛a termiczn ˛a. Pojawił si˛e on ju˙z wcze´sniej, nie wprost, w równaniu (3.1), dla przypomnienia:
a = λ ρcp , (3.11) gdzie: · λ - przewodno´s´c cieplna - [ W m·K], · ρ - g˛esto´s´c - [kg m3], · cp - ciepło wła´sciwe - [kg·KJ ].
Mianownik, czyli wyra˙zenie ρcp okre´sla obj˛eto´sciow ˛a pojemno´s´c ciepln ˛a. Odwrotno´s´c
współ-czynnika wyrównania temperatur, czyli 1a okre´sla bezwładno´s´c ciepln ˛a materiału. Jednostk ˛a współczynnika a jest [ms2]. Dla niektórych materiałów mo˙zna znale´z´c warto´sci dyfuzyjno´sci termicznej w tablicach fizycznych. W tabeli 3.1 znajduje si˛e kilka warto´sci dla przykładowych materiałów.
3.4. Współczynnik a 29
Tablica 3.1. Warto´sci dyfuzyjno´sci termicznej dla przykładowych materiałów
materiał a[ms2] złoto 1.27 · 10−4 silikon 8.8 · 10−5 aluminium 8.418 · 10−5 ˙zelazo 2.3 · 10−5 stal 304A w 27oC 4.2 · 10−6 kwarc 1.4 · 10−6 mur 5.2 · 10−7 szkło 3.4 · 10−7 woda w 25oC 1.43 · 10−7 alkohol 7 · 10−8
4. Obiekt eksperymentalny
W rozdziale tym zostanie opisana metodyka przeprowadzania eksperymentu. Według tego schematu w nast˛epnym rozdziale b˛ed ˛a sprawdzane kolejne rozwa˙zane modele, aby na ko´ncu wybra´c jeden najlepszy. Potrzebny jednak jest do tego schemat post˛epowania, którego przygo-towanie opisano wła´snie poni˙zej. Zaczyna si˛e on od wykonania eksperymentu na rzeczywistym obiekcie oraz uzyskaniu z niego danych pomiarowych. Nast˛epnie model jest implementowany w ´srodowisku MATLAB i w tym celu niezb˛edne było wybranie odpowiedniej metody nume-rycznej. Kolejnym krokiem jest identyfikacja parametryczna - podobnie jak poprzednio nale˙zało wybra´c algorytm identyfikacji. Na koniec obliczone zostaj ˛a wska´zniki statystyczne do oceny modelu.
4.1. Opis eksperymentu laboratoryjnego
W celu weryfikacji przyj˛etego modelu wykonany został eksperyment pomiarowy. Cegł˛e o wymiarach 120×250×65 mm ustawiono na grzałce elektrycznej i podgrzano z jednej strony od temperatury 16.5oC do 150oC. Zmiana temperatury zaj˛eła niecałe 7 minut. Proces ten został za-rejestrowany kamer ˛a termowizyjn ˛a SC-660 firmy FLIR SYSTEMS. Przy jej u˙zyciu wykonano seri˛e zdj˛e´c w kolejnych chwilach czasowych. Przykładowe zdj˛ecie przedstawiono na rysunku 4.1. Cało´s´c eksperymentu zarejestrowano na 214 zdj˛eciach zawieraj ˛acych informacje o rozkła-dzie temperatury. Na kolejnych uj˛eciach mo˙zna zaobserwowa´c post˛ep fali cieplnej. Na rysunku 4.2 przedstawiono 6 przykładowych zdj˛e´c z ró˙znych chwil czasowych. Odległo´s´c mi˛edzy ka-mer ˛a a podgrzewana cegł ˛a wynosi 50 cm (wzdłu˙z osi y z rysunku 4.3), a temperatura otoczenia 21.1oC. Rysunek 4.3 zawiera schemat stanowiska laboratoryjnego. Na podstawie zapisanych
zdj˛e´c otrzymano macierze z warto´sciami temperatury. Z tych macierzy wyodr˛ebniono wektory warto´sci temperatury wzdłu˙z osi szeroko´sci cegły, czyli osi x z rysunku 4.3. Z jednego zdj˛ecia otrzymano wektor temperatury od punktu x = 0 do punktu x = L (L = 120 mm). Wyod-r˛ebniaj ˛ac taki wektor z kolejnych zdj˛e´c otrzymano macierz rozkładu temparatury w czasie i przestrzeni. W pracy rozwa˙zany jest tylko rozkład temperatury wzdłu˙z osi x, jako
aproksyma-4.1. Opis eksperymentu laboratoryjnego 31
Rysunek 4.1. Pojedyncze zdj˛ecie rozkładu temperatury wykonane kamer ˛a termowizyjn ˛a SC-660.
Rysunek 4.2. Sze´s´c zdj˛e´c rozkładu temperatury z ró˙znych chwil czasowych.
cja rozkładu na powierzchni x×y. W przyszłych badaniach mo˙zna rozszerzy´c modele o rozkład tak˙ze w osi y. Otrzymane pomiary przedstawiono na wykresie 4.4.
4.2. Modelowanie 32
Rysunek 4.3. Schemat stanowiska laboratoryjnego.
Rysunek 4.4. Wykres przedstawiaj ˛acy zebrane pomiary.
4.2. Modelowanie
Przedstawione w rozdziale 5 modele matematyczne zostały zaimplementowane w MATLA-Bie, w celu przeprowadzenia symulacji numerycznych i identyfikacji parametrów. Konieczna
4.2. Modelowanie 33
była zamiana ich na posta´c dyskretn ˛a. W literaturze opisano kilka metod zamiany pochodnych ułamkowego rz˛edu na schematy ró˙znicowe dogodne do oblicze´n numerycznych. W artykule [47] porównano trzy metody aproksymuj ˛ace pochodn ˛a typu Caputo, na przykładzie ułamko-wego równania dyfuzji. W pracy [66] omówionych zostało kilka metod numerycznych: m.in. schematy ró˙znicowe oparte na definicji Grünwalda-Letnikova, metoda Lubicha dla wy˙zszych rz˛edów czy dwie metody przybli˙zania szeregiem: Taylora i dekompozycja Adomiana. Pocz ˛ at-kowo wi˛ekszo´s´c z nich dotyczyła głównie pochodnej typu Caputo i operowała na równaniach ró˙zniczkowych wzgl˛edem czasu. Ostatnio jednak w coraz wi˛ekszej ilo´sci prac mo˙zna spotka´c metody numeryczne dostosowane do pochodnych ułamkowych wzgl˛edem zmiennej przestrzen-nej. S ˛a to mi˛edzy innymi prace [12], [37] czy [38] gdzie rozwa˙zane s ˛a ró˙zne sposoby aprok-symacji pochodnej ułamkowej w sensie Riesza, czy [68] - praca doktorska, w której testowano zarówno znane metody jak i zaproponowano nowe, do ró˙znego typu równa´n ró˙zniczkowych cz ˛astkowych z pochodnymi niecałkowitego rz˛edu.
Opieraj ˛ac si˛e na testach i wnioskach przedstawionych w tych pracach, postanowiono do dalszych bada´n wybra´c metod˛e „Shifted Grünwald“ (SG) dla pochodnych Riesza oraz metod˛e GMMP dla pochodnych Caputo.
4.2.1. Pochodna Riesza
Metoda SG została po raz pierwszy przedstawiona w [38], gdzie były weryfikowane nu-meryczne schematy dla pewnej klasy ułamkowych cz ˛astkowych równa´n ró˙zniczkowych. Jest oparta na definicji lewo- i prawo- stronnej pochodnej ułamkowej Grünwalda-Letnikova, ale lekko zmodyfikowanej: (4.1) i (4.2). U˙zywana b˛edzie nast˛epuj ˛aca notacja: xj = j∆x, tm =
m∆t i u(xj, tm) ' u(m)j , gdzie u (m)
j oznacza numeryczn ˛a estymacj˛e warto´sci dokładnej u(x, t)
dla x = xj i t = tm. 0Dxβju(xj, tm) = 1 (∆x)β j+1 X r=0 ωrβu(m)j−r+1+ O(∆x), (4.1) LDβxju(xj, tm) = 1 (∆x)β J −j+1 X r=0 ωβru(m)j+r−1+ O(∆x). (4.2)
Aproksymacja pochodnej ułamkowej dla tej metody mo˙ze by´c zatem wyra˙zona nast˛epuj ˛aco:
∂βu(x, t) ∂|x|β ≈ −(∆x)−β 2 cosπβ2 "j+1 X r=0 ωβru(m)j−r+1+ J −j+2 X r=0 ωrβu(m)j+r−1 # , (4.3)
przy czym współczynniki ωrβliczone s ˛a według nast˛epuj ˛acej zasady (podstawiaj ˛ac odpowiednio β za γ oraz r za i):
ω0γ = 1, ωγi = (−1)iγ(γ − 1) . . . (γ − i + 1)
i! dla i = 1, 2, . . . (4.4)
4.2. Modelowanie 34
Rysunek 4.5. Przykład RFDE z parametrami: β = 1.3, Kβ = 0.25, oraz warunkiem pocz ˛
atko-wym: g(x) = sin2(2x).
Przykład Rozwa˙zmy ułamkowe równanie dyfuzji Riesza (Riesz fractional diffusion equation) - w skrócie RFDE:
∂
∂tu(x, t) = Kβ ∂β
∂|x|βu(x, t) , (4.5)
gdzie 0 ≤ t ≤ T, 0 < x < L, 1 < β ≤ 2, Kβ 6= 0, oraz warunki brzegowe (4.6) i pocz ˛atkowe
(4.7):
u(0, t) = u(L, t) = 0 , (4.6)
u(x, 0) = g(x) . (4.7)
Analityczne rozwi ˛azanie mozna zapisa´c (4.8) i (4.9).
u(x, t) = ∞ X n=1 bnsin nπx L exp −Kβ nπ L β · t , (4.8) gdzie: bn = 2 L Z L 0 g(ξ) sinnπξ L dξ (4.9)
Na rysunku 4.5 pokazano przykładowe rozwi ˛azanie RFDE, przy warunku pocz ˛atkowym g(x) = sin2(2x), warto´sciach parametrów β = 1.3, Kβ = 0.25 oraz dla L = π i T = 5.
4.2. Modelowanie 35
Rysunek 4.6. Porównanie rozwi ˛azania analitycznego i numerycznego metod ˛a SG równania RFDE z parametrami: β = 1.3, Kβ = 0.25, oraz warunkiem pocz ˛atkowym:
g(x) = sin2(2x).
Korzystaj ˛ac teraz z metody numerycznej SG mo˙zna równanie RFDE aproksymowa´c nast˛e-puj ˛aco: u(m+1)j = −Kβ(∆t)(∆x) −β 2 cosπβ2 "j+1 X r=0 ωβru(m)j−r+1+ J −j+2 X r=0 ωrβu(m)j+r−1 # + u(m)j , (4.10)
Na rysunku 4.6 pokazano porównanie przykładowego rozwi ˛azania RFDE z jego aproksymacj ˛a metod ˛a SG.
4.2.2. Pochodna Caputo
Pochodna ułamkowa po czasie typu Caputo b˛edzie aproksymowana według schematu za-proponowanego w pracy [18] i nazwanego tam GMMP (od nazwisk autorów pracy). U˙zywaj ˛ac tej samej notacji, co dla poprzedniego schematu, mo˙zna go zapisa´c nast˛epuj ˛aco:
∂αu(x, t) ∂tα ≈ 1 (∆t)α m X k=0 ωkαhu(m−k)j − u(0)j − u0(0)j tm i , (4.11)
przy czym współczynniki ωkαdefiniuje si˛e tak˙ze korzystaj ˛ac ze wzoru (4.4).
4.2. Modelowanie 36
Rysunek 4.7. Przykład FDWE z parametrem: α = 1.7 oraz warunkiem pocz ˛atkowym: g(x) = sin(x).
Przykład W tym przykładzie rozwa˙zane b˛edzie ułamkowe równanie dyfuzji fali (fractional diffusion-wave equation) - w skrócie FDWE:
∂α
∂tαu(x, t) =
∂2
∂x2u(x, t) , (4.12)
gdzie 0 ≤ t ≤ T, 0 < x < L, 1 < α ≤ 2, oraz warunki brzegowe (4.13) i pocz ˛atkowe (4.14):
u(0, t) = u(L, t) = 0 , (4.13) u(x, 0) = g(x), ∂u(x, t) ∂t t=0 = h(x). (4.14)
Je˙zeli przyjmie si˛e: L = π, g(x) = sin(x) oraz h(x) = 0, to analityczne rozwi ˛azanie mo˙zna zapisa´c w postaci równania (4.15).
u(x, t) = Eα(−tα)sin(x) , (4.15)
gdzie Eα to funkcja Mittaga-Lefflera [42]. Na rysunku 4.7 pokazano przykładowe rozwi ˛azanie
FDWE, przy warto´sciach α = 1.7 i T = 5.
Korzystaj ˛ac teraz z metody numerycznej GMMP mo˙zna równanie FDWE aproksymowa´c nast˛epuj ˛aco: u(m+1)j = (∆t) α (∆x)2 h u(m)j+1− 2u(m)j + u(m)j−1i− m X k=1 ωkαu(m+1−k)j + u(0)j m X k=0 ωαk + u0(0)j tmωkα, (4.16)
4.2. Modelowanie 37
Rysunek 4.8. Porównanie rozwi ˛azania analitycznego i numerycznego metod ˛a GMMP równa-nia FDWE z parametrem: α = 1.7 oraz warunkiem pocz ˛atkowym: g(x) = sin(x).
Na rysunku 4.8 pokazano porównanie przykładowego rozwi ˛azania FDWE z jego aproksymacj ˛a metod ˛a GMMP.
Podsumowuj ˛ac, w badanych modelach, pochodne ułamkowe typu Riesza, b˛ed ˛a aproksymo-wane metod ˛a Shifted-Grünwald według wzoru (4.3), pochodne ułamkowe typu Caputo metod ˛a GMMP (4.11), natomiast pochodne klasyczne, pierwszego i drugiego rz˛edu, b˛ed ˛a przybli˙zane przez schematy ró˙znicowe odpowiednich stopni (4.17) i (4.18):
∂ ∂tu(xj, tm) ' u(xj, tm+ 1) − u(xj, tm) ∆t , (4.17) ∂2 ∂x2u(xj, tm) '
u(xj+1, tm) − 2u(xj, tm) + u(xj−1, tm)
(∆x)2 . (4.18)
Aby, na podstawie modelu, otrzyma´c estymat˛e przepływu ciepła w testowanym materiale, ko-nieczna jest znajomo´s´c warunków pocz ˛atkowych i brzegowych. Wektor warto´sci temperatury wzdłu˙z osi x cegły (rysunek 4.3) z pocz ˛atku eksperymentu został u˙zyty jako warunek pocz ˛ at-kowy do symulacji. Warunek brzegowy otrzymano z serii wszystkich zdj˛e´c z punktu (x = 0) na brzegu cegły.
4.3. Identyfikacja parametryczna 38
4.3. Identyfikacja parametryczna
Nast˛epnym etapem była identyfikacja parametrów badanego modelu. Przeprowadzono se-ri˛e testów, aby wybra´c odpowiedni algorytm. Poni˙zej opisano cz˛e´s´c z nich. Najwa˙zniejsza, z punktu widzenie dalszych bada´n, jest odporno´s´c algorytmu na niepewno´sci pomiarowe obecne w wynikach eksperymentu. Głównie pod tym k ˛atem były konstruowane testy. Nale˙zy zwróci´c uwag˛e, ˙ze ich wynik ´swiadczy tylko o przydatno´sci metody dla konkretnego problemu, rozwa-˙zanego w niniejszej pracy, nie okre´sla natomiast ogólnych i uniwersalnych zasad, obowi ˛azuj ˛ a-cych dla dowolnego zagadnienia.
4.3.1. Równania ró˙zniczkowe
Testy przeprowadzono na dwóch ułamkowych równaniach ró˙zniczkowych cz ˛astkowych:
1. Równaniu adwekcji-dyspersji Riesza (Riesz fractional advection-dispertion equation) - w skrócie RFADE.
2. Równaniu dyfuzji Riesz (Riesz fractional diffusion equation) - w skrócie RFDE.
Oba te równania zostały wybrane do testów, poniewa˙z s ˛a one równaniami ró˙zniczkowymi cz ˛ ast-kowymi z pochodn ˛a w sensie Riesza, a zatem tak ˛a jaka wyst˛epuje w badanych modelach, a dodatkowo znane jest ich analityczne rozwi ˛azanie, co pozwoli na weryfikacje algorytmów identyfikacji. Poni˙zej podane zostan ˛a te równania wraz z ich warunkami granicznymi oraz roz-wi ˛azaniami analitycznymi. Równanie RFDE było ju˙z opisane w poprzednim podrozdziale 4.2, w przykładzie. Równanie RFADE mo˙zna zapisa´c nast˛epuj ˛aco:
∂ ∂tu(x, t) = Kα ∂α ∂|x|αu(x, t) + Kβ ∂β ∂|x|βu(x, t) , (4.19)
gdzie 0 ≤ t ≤ T, 0 < x < L, 1 < α ≤ 2, 0 < β < 1, Kα 6= 0, Kβ 6= 0, oraz warunki brzegowe
(4.20) i pocz ˛atkowe (4.21):
u(0, t) = u(L, t) = 0 , (4.20)
u(x, 0) = g(x) . (4.21)
Rozwi ˛azanie analityczne (4.22) i (4.23):
u(x, t) = ∞ X n=1 bnsin nπx L exp − Kα nπ L α + Kβ nπ L β t , (4.22) gdzie: bn = 2 L Z L 0 g(ξ) sinnπx L dξ (4.23)
4.3. Identyfikacja parametryczna 39
Równanie to jest u˙zywane do modelowania transportu znaczników pasywnych w postaci płynu w o´srodkach porowatych lub transportu w materiałach podpowierzchniowych [68]. Na rysunku 4.9 przedstawiono przykładowe rozwi ˛azanie równania RFADE. Jak mo˙zna zauwa˙zy´c, RFDE
Rysunek 4.9. Przykład RFADE z parametrami: α = 1.8, Kα = 0.25, β = 0.2, Kβ = 0.25, oraz
warunkiem pocz ˛atkowym: g(x) = x2(π − x).
mo˙ze by´c szczególnym przypadkiem prezentowanego powy˙zej RFADE, gdy przyjmiemy Kβ =
0.
4.3. Identyfikacja parametryczna 40
4.3.2. Algorytmy identyfikacji
Zadaniem algorytmu identyfikacji b˛edzie zatem znalezienie wła´sciwych warto´sci parame-trów: α, Kα, β, Kβ dla RFADE oraz α, Kα dla RFDE, tak aby zminimalizowa´c kwadratowe
kryterium jako´sci Q: Q = Z L 0 Z T 0 ε2(x, t)dtdx , (4.24)
gdzie ε(x, t) wyra˙za bł ˛ad miedzy stanem otrzymanym z pomiarów u(x, t) a estymat ˛a obliczon ˛a z modelu u(x, t):
ε(x, t) = u(x, t) − u(x, t) . (4.25)
Mo˙zna zauwa˙zy´c, ˙ze rozwa˙zany problem jest nieliniowy ze wzgl˛edu na α i β, dlatego ko-nieczne jest u˙zycie algorytmów nieliniowych. Porównane zostały trzy metody: algorytm Le-venberga–Marquardta (LMA), algorytm Gaussa–Newtona (GNA) oraz metoda Simplex oparta na algorytmie Neldera-Meada. Wszystkie one s ˛a powszechnie stosowane do rozwi ˛azywania nieliniowego zagadnienia najmniejszych kwadratów. Teoretycznie LMA jest bardziej odporny ni˙z GNA - gwarantuje znalezienie optymalnych parametrów przy starcie z dowolnego punktu pocz ˛atkowego, nawet bardzo daleko od poprawnego rozwi ˛azania. Z drugiej jednak strony GNA w niektórych przypadkach jest du˙zo szybszy. Poni˙zej skrótowo przedstawiono idee działania tych algorytmów. Przyjmijmy wektor parametrów:
θRF ADE = [α, Kα, β, Kβ]T dla RFADE oraz θRF DE= [α, Kα]T dla RFDE. (4.26)
Załó˙zmy tak˙ze: t = kτ , x = mξ, kryterium Q mo˙zna wtedy zapisa´c nast˛epuj ˛aco:
Q = K X k=0 M X m=0 ε2m,k , εm,k = um,k− um,k , (4.27)
gdzie: um,k = u(mξ, kτ ) i um,k = u(mξ, kτ ). Algorytmy minimalizuj ˛a kryterium (4.27),
esty-muj ˛ac w ka˙zdej kolejnej iteracji wektor parametrów (4.26), a˙z do uzyskania rozwi ˛azania speł-niaj ˛acego zadane warunki „STOP“. Mo˙ze to by´c maksymalna liczba iteracji, czas oblicze´n, mi-nimalna zmiana warto´sci funkcji celu lub rozwi ˛azania, a tak˙ze dowolne ich kombinacje. GNA został zaimplementowany w MATLABie przez autork˛e zgodnie ze wzorem (4.28):
θGNi+1 = θGNi − (Jθθ00)−1Jθ0θ=θGN i
, (4.28)
gdzie Jθ0 jest gradientem, natomiast Jθθ00 jest aproksymacj ˛a hesjanu funkcji σθi
m,k, nazywanej
funkcj ˛a wra˙zliwo´sci wyj´scia, a wyra˙zaj ˛acej si˛e przy u˙zyciu wzorów (4.29) - (4.31).
Jθ0 = −2 K X k=0 M X m=0 εm,kσθm,ki , (4.29)
4.3. Identyfikacja parametryczna 41 Jθθ00 ≈ 2 K X k=0 M X m=0 σθi m,k(σ θi m,k) T , (4.30) σθi m,k = ∂um,k ∂θi . (4.31)
Pozostałe dwa algorytmy, tj. LMA i Simplex były testowane przy u˙zyciu wbudowanych funkcji pakietu MATLAB: fminsearch i lsqnonlin. Szczegółowe ich opisy mo˙zna znale´z´c w literaturze, np. [43], [35] dla algorytmu Levenberga-Marquadta oraz [31] dla algorytmu Simplex Neldera-Meada.
4.3.3. Testy numeryczne
Do testów u˙zyto dwóch, ju˙z wcze´sniej wspomnianych, równa´n (4.19) i (4.5) z warunkiem pocz ˛atkowym:
u(x, 0) = g(x) = x2(π − x) , (4.32)
oraz warto´sciami brzegowymi: L = π, T = 5. W celu zasymulowania niepewno´sci pomiaro-wych, zostały wygenerowane zakłócenia, w postaci szumu białego, o odchyleniu standardowym δ, które dodano na wyj´scie modelu:
um,k = um,k+ Rm,k , (4.33)
gdzie R jest macierz ˛a zawieraj ˛ac ˛a losowe warto´sci o rozkładzie normalnym, z parametrami: µ = 0 - warto´s´c ´srednia oraz δ - odchylenie standardowe - zmieniane na potrzeby testów.
4.3.3.1. Test 1
W pierwszym eksperymencie sprawdzano odporno´s´c metody na zakłócenia na przykładzie RFADE. Jest to wa˙zne, poniewa˙z w pracy, identyfikacja b˛edzie przeprowadzana na podstawie pomiarów rzeczywistych, a zatem obci ˛a˙zanych bł˛edem pomiarowym. Symulacje wykonywano dla ró˙znych warto´sci odchylenia standardowego szumu δ. W tabelach 4.1 - 4.3 znajduje si˛e ze-stawienie otrzymanych wyników. Warunkiem „STOP“ zatrzymania algorytmu był brak zmiany warto´sci kryterium Q w kolejnych iteracjach o wi˛ecej ni˙z: Qtol = 10−6 i był on ustawiony
tak samo dla wszystkich algorytmów. W tabelach umieszczono cztery zwracane warto´sci: roz-wi ˛azanie - czyli zidentyfikowane parametry θ, czas obliczania, ilo´s´c wykonanych iteracji oraz warto´s´c kryterium Q.
Jak mo˙zna zauwa˙zy´c algorytm Simplex był to najszybsz ˛a i najdokładniejsz ˛a metod ˛a, nieza-le˙znie od poziomu szumu. Algorytm Gaussa-Newtona zwracał tak˙ze do´s´c dobre rozwi ˛azania, zwłaszcza dla małego poziomu szumu, jednak czas oblicze´n był wi˛ekszy. By´c mo˙ze mo˙zna by zmniejszy´c ten czas poprzez optymalizacj˛e kodu. Metoda LMA zawsze zwracała to samo złe
4.3. Identyfikacja parametryczna 42
Tablica 4.1. Podsumowanie testów dla poziomu szumu δ = 0 - brak zakłóce´n.
parametry LMA GNA Simplex
θ = α Kα β Kβ 1.4106 0.3921 0.1255 0.1079 1.7831 0.2623 0.1000 0.2377 1.7855 0.2576 0.1395 0.2424 czas[s] 77.72 93.93 33.74 iteracje 126 162 185
Q 2.2479E-4 2.4126E-5 2.6396E-7
Tablica 4.2. Podsumowanie testów dla poziomu szumu δ = 0.01.
parametry LMA GNA Simplex
θ = α Kα β Kβ 1.4104 0.3915 0.1255 0.1086 1.7659 0.2668 0.1000 0.2331 1.7878 0.2573 0.1548 0.2428 czas[s] 160.39 51.45 33.41 iteracje 261 90 173 Q 0.1817 0.1812 0.1643
Tablica 4.3. Podsumowanie testów dla poziomu szumu δ = 0.02.
parametry LMA GNA Simplex
θ = α Kα β Kβ 1.4121 0.3924 0.1254 0.1078 1.9961 0.2087 0.1000 0.2911 1.7959 0.2502 0.1540 0.2495 czas[s] 110.05 103.15 32.38 iteracje 179 181 174 Q 0.6452 0.6682 0.0676
rozwi ˛azanie - mo˙ze to by´c minimum lokalne lub algorytm ten jest wolniejszy od pozostałych. Zostanie to sprawdzone w kolejnym eksperymencie.
4.3. Identyfikacja parametryczna 43
4.3.3.2. Test 2
W drugim eksperymencie sprawdzano wpływ doboru punktu startowego na przykładzie RFDE. Ustalono warto´s´c odchylenia standardowego: δ = 0.1. Warunek zako´nczenia identy-fikacji przyj˛eto tak samo jak w poprzednim te´scie. Porównywane b˛ed ˛a teraz czas oblicze´n i warto´s´c kryterium Q. Identyfikacja startowała z ró˙znych punktów pocz ˛atkowych θ: α ∈ [1.1; 2] z krokiem = 0.1 oraz Kα ∈ [0.1; 1] z krokiem = 0.1. Na wykresach 4.10 - 4.15 przedstawiono
uzyskane rezultaty.
Rysunek 4.10. Kryterium Q w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Gaussa-Newtona.
Ponownie metoda Simplex wypadła najlepiej. Jest najszybsza i do´s´c dokładna. LMA jest wolniej zbie˙zna, wi˛ec aby otrzyma´c porównywalne wyniki, potrzebowała znacznie wi˛ekszej ilo´sci iteracji i czasu. Metoda GNA jest dobra, ale tylko dla punktów le˙z ˛acych blisko rozwi ˛ a-zania, dla bardziej odległych rozwi ˛aza´n zwraca złe wyniki. W tabeli 4.4 znajduje si˛e podsumo-wanie drugiego testu.
4.3.4. Wybór metody
Po przeprowadzeniu testów, zdecydowano si˛e wybra´c, do dalszych bada´n, metod˛e Simplex opart ˛a na algorytmie Neldera-Meada. Jest ona najbardziej odporna na zakłócenia i zły wybór
4.3. Identyfikacja parametryczna 44
Rysunek 4.11. Czas oblicze´n w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Gaussa-Newtona.
Rysunek 4.12. Kryterium Q w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Levenberga-Marquardta.
4.3. Identyfikacja parametryczna 45
Rysunek 4.13. Czas oblicze´n w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Levenberga-Marquardta.
Rysunek 4.14. Kryterium Q w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Sim-plex Neldera-Meada.
punktu startowego, co jest istotne w planowanych symulacjach. Testy wykazały najszybszy czas jej działania przy jednocze´snie poprawnym rozwi ˛azywaniu postawionego problemu.
4.4. Weryfikacja i porównanie 46
Rysunek 4.15. Czas oblicze´n w funkcji pocz ˛atkowych warto´sci parametrów dla algorytmu Simplex Neldera-Meada.
Tablica 4.4. Podsumowanie testu 2.
GNA LMA Simplex
´sredni czas oblicze´n 2 s 150 s 8 s ´srednia warto´s´c kryterium Q 1000 0.027 0.168
Wracaj ˛ac do metodyki, po przeprowadzeniu identyfikacji, otrzymujemy warto´sci parame-trów modelu. Nale˙zy teraz zweryfikowa´c poprawno´s´c otrzymanych rezultatów. W tym celu wyznaczone zostan ˛a wska´zniki statystyczne, opisane w kolejnym podrozdziale.
4.4. Weryfikacja i porównanie
W celu weryfikacji i oceny modelu obliczanych zostanie kilka wska´zników statystycznych. Poni˙zej wymieniono i krótko omówiono ka˙zdy z nich.
1. Maksymalny bł ˛ad procentowy - εmax:
εmax= max n∈N | bTn− Tn| b Tn · 100% ! , (4.34)
4.4. Weryfikacja i porównanie 47
2. ´Sredni (procentowy) bł ˛ad wzgl˛edny - MRE (Mean Relative Error):
M RE = 1 N N X n=1 εn, (4.35)
który operuje na całej dziedzinie. Lepiej oddaje ogólne dopasowanie modelu.
3. Bł ˛ad procentowy dla modelu. Jest on obliczony dla ka˙zdego punktu zgodnie ze wzorem:
εn=
| bTn− Tn|
b Tn
· 100%, (4.36)
i przedstawiany b˛edzie na wykresie dla wszystkich punktów dost˛epnych z pomiarów.
Przedstawione powy˙zej wska´zniki posłu˙z ˛a do porównania wszystkich badanych modeli i wy-brania najlepszego z nich.
Podsumowuj ˛ac metodyk˛e, przeprowadzana b˛edzie, dla ka˙zdego modelu, nast˛epuj ˛aca proce-dura: implementacja numerycznej aproksymacji modelu według schematów opisanych w roz-dziale 4.2, identyfikacja parametrów na podstawie pomiarów z eksperymentu, ocena za pomoc ˛a wska´zników statystycznych. Wyniki zostan ˛a zaprezentowane w nast˛epnym rozdziale 5. Po wy-konaniu tych operacji, modele zostan ˛a porównane i wybrany b˛edzie najlepszy, tj. taki, który naj-´sci´slej odpowiada rzeczywistemu, obserwowanemu zjawisku. Maj ˛ac ju˙z wybrany model dalsze badania b˛ed ˛a przeprowadzane tylko na nim.
5. Porównanie modeli
W rozdziale tym przedstawiono wyniki bada´n przeprowadzonych według schematu, opi-sanego w poprzednim rozdziale, dla ka˙zdego z pi˛eciu testowanych modeli: czterech modeli niecałkowitego rz˛edu oraz, dla porównania, modelu klasycznego. Oprócz parametrów otrzy-manych z identyfikacji modeli, podano równie˙z obliczone warto´sci wska´zników statystycznych oraz wykresy zarówno samych modeli jak i ich bł˛edów w stosunku do pomiarów. Na ko´ncu rozdziału dokonano krótkiego porównania otrzymanych wyników i wybrano najlepszy model.
5.1. Model I
Model ten to czasowe, ułamkowe równanie przewodnictwa cieplnego (3.6): ∂αT (x, t)
∂tα = a
∂2T (x, t)
∂x2 , 0 < α < 2, (5.1)
z pochodn ˛a ułamkow ˛a po czasie w sensie Caputo. Warunki pocz ˛atkowe i brzegowe zostały wzi˛ete z pomiarów, jak opisano w rozdziale 4.2. Na podstawie metod opisanych równie˙z w tym rozdziale oraz przy u˙zyciu tych samych oznacze´n, model został przybli˙zony schematem numerycznym: 1 (∆t)α m X r=0 ωkαhTj(m−k)− Tj(0)− Tj0(0)tm i = aT (m) j+1 − 2T (m) j + T (m) j−1 (∆x)2 , 0 < α < 2. (5.2)
Po wykonaniu symulacji i identyfikacji na podstawie pomiarów z eksperymentu, otrzymano nast˛epuj ˛ace warto´sci parametrów:
a = 2.28 · 10−7, α = 1.48, (5.3)
natomiast warto´s´c kryterium (4.24) dla tych parametrów wyniosła QI = 1.3073406 · 106. Na
rysunku 5.1 pokazano model I dla otrzymanych warto´sci parametrów. Zidentyfikowany model został porównany z pomiarami otrzymanymi z eksperymentu. Na rysunku 5.2 przedstawiono wykres bł˛edu procentowego dla tego modelu. Jest on obliczony dla ka˙zdego punktu zgodnie ze wzorem (4.36).
5.1. Model I 49
Rysunek 5.1. Model I dla warto´sci parametrów: a = 2.28 · 10−7oraz α = 1.48.
Rysunek 5.2. Wykres bł˛edu procentowego εndla ułamkowego modelu I.
5.2. Model II 50
Nast˛epnie zostały policzone pozostałe wska´zniki niedokładno´sci. Wyniosły one odpowied-nio:
εImax= 12.0648%, (5.4)
M REI = 2.9258%. (5.5)
5.2. Model II
Model II to przestrzenne ułamkowe równanie przewodnictwa cieplnego (3.7):
∂T (x, t) ∂t = a
∂βT (x, t)
∂|x|β , 0 < β 6 2, (5.6)
z pochodna ułamkow ˛a po zmiennej przestrzennej w sensie Riesza oraz warunkami granicznymi jak poprzednio.. Na podstawie metod opisanych w rozdziale 4.2 oraz przy u˙zyciu tych samych oznacze´n, model został przybli˙zony schematem numerycznym:
Tj(m+1) − Tj(m) ∆t = −a(∆x)−β 2 cosπβ2 "j+1 X r=0 ωrβTj−r+1(m) + J −j+2 X r=0 ωrβTj+r−1(m) # , 0 < β 6 2. (5.7) Po wykonaniu symulacji i identyfikacji na podstawie pomiarów z eksperymentu, otrzymano nast˛epuj ˛ace warto´sci parametrów:
a = 9.61 · 10−8, β = 1.00, (5.8)
natomiast warto´s´c kryterium 4.24 dla tych parametrów wyniosła QII = 1.1779685 · 106. Na
rysunku 5.3 pokazano model II dla otrzymanych warto´sci parametrów. Zidentyfikowany model został porównany z pomiarami otrzymanymi z eksperymentu. Na rysunku 5.4 przedstawiono wykres bł˛edu procentowego dla tego modelu. Jest on obliczony dla ka˙zdego punktu zgodnie ze wzorem (4.36).
Nast˛epnie zostały policzone pozostałe wska´zniki niedokładno´sci. Policzone bł˛edy wyniosły odpowiednio:
εIImax= 10.6120%, (5.9)
M REII = 2.1497%. (5.10)
5.3. Model III
Model III to czasowo-przestrzenne ułamkowe równanie przewodnictwa cieplnego (3.9):
∂αT (x, t)
∂tα = a
∂βT (x, t)
5.3. Model III 51
Rysunek 5.3. Model II dla warto´sci parametrów: a = 9.61 · 10−8oraz β = 1.00.
Rysunek 5.4. Wykres bł˛edu procentowego εndla ułamkowego modelu II.
5.3. Model III 52
z pochodna ułamkow ˛a po czasie w sensie Caputo oraz pochodn ˛a ułamkow ˛a po x w sensie Riesza. Warunki pocz ˛atkowe i brzegowe wzi˛eto z pomiarów, jak poprzednio. Analogicznie jak poprzednio model został przybli˙zony schematem numerycznym:
1 (∆t)α m X r=0 ωkα h Tj(m−k)− Tj(0)− Tj0(0)tm i = −a(∆x) −β 2 cosπβ2 "j+1 X r=0 ωβrTj−r+1(m) + + J −j+2 X r=0 ωrβTj+r−1(m) # , 0 < α 6 2 0 < β 6 2. (5.12) Po wykonaniu symulacji i identyfikacji na podstawie pomiarów z eksperymentu, otrzymano nast˛epuj ˛ace warto´sci parametrów:
a = 6.57 · 10−9, α = 2.00, β = 2.00, (5.13)
natomiast warto´s´c kryterium (4.24) dla tych parametrów wyniosła QIII = 4.3654473 · 105. Na
rysunku 5.5 pokazano model III dla otrzymanych warto´sci parametrów. Zidentyfikowany model
Rysunek 5.5. Model III dla warto´sci parametrów: a = 6.57 · 10−9, α = 2.00 oraz β = 2.00
został porównany z pomiarami otrzymanymi z eksperymentu. Na rysunku 5.6 przedstawiono wykres bł˛edu procentowego dla tego modelu. Jest on obliczony dla ka˙zdego punktu zgodnie ze wzorem (4.36).
5.4. Model IV 53
Rysunek 5.6. Wykres bł˛edu procentowego εndla ułamkowego modelu III.
Nast˛epnie zostały policzone pozostałe wska´zniki niedokładno´sci i wyniosły one odpowied-nio:
εIIImax = 6.5591%, (5.14)
M REIII = 1.4663%. (5.15)
5.4. Model IV
Model IV został zaproponowany przez autork˛e, na podstawie obserwacji wcze´sniej uzyski-wanych wyników. Jest to modyfikacja modelu klasycznego o czynnik zawieraj ˛acy pochodn ˛a ułamkow ˛a w sensie Caputo:
∂T (x, t) ∂t + c ∂αT (x, t) ∂tα = ∂2T (x, t) ∂x2 , 0 < α 6 2, a, c ∈ R, (5.16)
z pochodna ułamkow ˛a po czasie w sensie Caputo oraz pochodnymi klasycznymi całkowitego rz˛edu. Warunki graniczne sa takie same jak dla poprzednich modeli. Model został przybli˙zony
5.4. Model IV 54 schematem numerycznym: 1 ∆t h Tj(m+1)− Tj(m)i+ c (∆t)α m X r=0 ωkα h Tj(m−k)− Tj(0)− Tj0(0)tm i = = a (∆x)2 h Tj+1(m)− 2Tj(m)+ Tj−1(m) i , 0 < α 6 2 a, c ∈ R. (5.17) Po wykonaniu symulacji i identyfikacji na podstawie pomiarów z eksperymentu, otrzymano nast˛epuj ˛ace warto´sci parametrów:
a = 3.25 · 10−7, c = 2.89 · 10−3, α = 0.19, (5.18)
natomiast warto´s´c kryterium (4.24) dla tych parametrów wyniosła QIV = 8.6098516 · 104. Na
rysunku 5.7 pokazano model IV dla otrzymanych warto´sci parametrów. Zidentyfikowany model
Rysunek 5.7. Model IV dla warto´sci parametrów: a = 3.25·10−7, c = 2.89·10−3oraz α = 0.19
został porównany z pomiarami otrzymanymi z eksperymentu. Na rysunku 5.8 przedstawiono wykres bł˛edu procentowego dla tego modelu. Jest on, tak jak poprzednio, obliczony dla ka˙zdego punktu zgodnie ze wzorem (4.36).
Nast˛epnie zostały policzone pozostałe wska´zniki niedokładno´sci. Wyniosły one odpowied-nio:
εIVmax = 5.6331%, (5.19)
5.5. Model klasyczny 55
Rysunek 5.8. Wykres bł˛edu procentowego εndla ułamkowego modelu IV.
5.5. Model klasyczny
Dla porównania, identyfikacja została równie˙z wykonana dla modelu klasycznego, opartego na prawie Fouriera (3.5).
∂T (x, t) ∂t = a
∂2T (x, t)
∂x2 , (5.21)
Model został przybli˙zony schematem numerycznym z wykorzystaniem odpowiednich schema-tów ró˙znicowych: 1 ∆t h Tj(m+1)− Tj(m)i= a (∆x)2 h Tj+1(m)− 2Tj(m)+ Tj−1(m)i. (5.22)
Po wykonaniu symulacji i identyfikacji na podstawie pomiarów z eksperymentu, otrzymano nast˛epuj ˛ace warto´sci parametrów:
a = 1.64 · 10−7. (5.23)
Zidentyfikowany model został porównany z pomiarami otrzymanymi z eksperymentu. Na rysunku 5.10 przedstawiono wykres bł˛edu procentowego dla tego modelu. Jest on obliczony dla ka˙zdego punktu zgodnie ze wzorem (4.36). Natomiast warto´s´c kryterium (4.24) dla tego para-metru wyniosła QK = 1.1575619 · 106. Na rysunku 5.9 pokazano model klasyczny dla
otrzyma-nych warto´sci parametrów. Nast˛epnie zostały policzone pozostałe wska´zniki niedokładno´sci:
εKmax= 10.5502%, (5.24)
