4. METODY APROKSYMACYJNEGO ROZWI
ZYWANIA
RÓWNAC
RÓŻNICZKOWYCH
4.1. Uwagi wstępne i koncepcje podstawowe
W wielu przypadkach, ważnych z punktu widzenia zastosowań technicznych, znalezienie
rozwiązania równania różniczkowego w postaci analitycznej jest trudne lub wręcz niemożliwe do osiągnięcia.
Taka sytuacja towarzyszy znacznej większości przypadków, szczególnie gdy w problemach mechaniki
ośrodka ciągłego mamy do czynienia z nieregularną geometrią układów, bądz
gdy właściwości materiałów
są funkcjami zmiennych pól. Takie problemy nazywamy zadaniami nieliniowymi. Metoda elementów
skończonych często w środowisku inżynierów bywa utożsamiana ze sposobem rozwiązywania problemów
mechaniki, przy czym zapomina się że wszystkie te zagadnienia modelowane są za pomocą równań
różniczkowych. Poniżej pragniemy pokazać, że metoda elementów skończonych jest przede wszystkim
metodą rozwiązywania dowolnych równań różniczkowych.
Spośród metod znajdowania rozwiązań przybliżonych, które są użyteczne dla lepszego zrozumienia
istoty MES, chcemy zwrócić uwagę na następujące:
- metoda Ritza,
- metoda wariacyjna Rayleigha - Ritza,
- metody ważonych reziduów.
Wymieniane powyżej metody znajdywania przybliżonych rozwiązań równań różniczkowych są
metodami całkowymi, gdyż opierają się na formach całkowych, nie zaś bezpośrednio na równaniach
różniczkowych.
Metoda Ritza jest na tyle prosta, że nie wymaga dodatkowych informacji poszerzających klasyczny
kurs matematyki. Metody wariacyjne, operujące pojęciami klasycznego rachunku wariacyjnego, wymagają
choćby znajomości podstawowych formuł tego rachunku. Nie będziemy ich tutaj przedstawiali, jedynie za-
proponujemy Czytelnikowi zapoznanie siÄ™ z fragmentami literatury z
ródłowej.
Zanim przedstawimy metody Ritza i Rayleigha-Ritza, spróbujmy na przykładzie prostego problemu
brzegowego prześledzić ideę koncepcji generalnej. Nasze rozważania zilustrujemy przykładem prostego
przypadku pojedynczego równania różniczkowego z jedną tylko zmienną niezależną. By przedstawić różnice
między kolejno omawianymi metodami będziemy analizować zawsze to samo równanie różniczkowe,
którego rozwiązanie analityczne jest proste do osiągnięcia. Tym samym łatwo nam będzie porównać
dokładność otrzymywanych przez aproksymację wyników z rozwiązaniami dokładnymi.
Przyjmijmy równanie różniczkowe zapisane symbolicznie:
f (T (x)) = 0 (4.1)
21 Metody aproksymacyjnego rozwiązywania równań różniczkowych
opisujące dowolne zagadnienie w przestrzeni &! , gdzie T opisuje funkcję zmian x (może to być np.
temperatura, funkcja ugięcia belki), zaś &! jest obszarem działania funkcji rządzonej przez prawo
zdefiniowane za pomocą operatora różniczkowego f .
Warunki brzegowe zapiszemy w postaci:
g1(T (x)) = 0 na obszarze “1 ,
g2 (T (x)) = 0 na obszarze “2 , (4.2)
gdzie “1 i “2 zawierajÄ… te części &! , które ograniczajÄ… brzeg ( rys. 4.1).
Rys. 4.1 Obszar jednowym iarowy działania funkcji f (T (x))
Określmy rozwiązanie równania (4.1) łącznie z warunkami (4.2) w postaci funkcji:
n
T '= T '(x,a1,a2,...,an ) = Å" Ni (x) , (4.3)
"ai
i=1
która ma jedną bądz
więcej niewiadomych ai , zaś funkcje Ni spełniają dokładnie warunki
brzegowe (4.2). Funkcje Ni nazywamy funkcjami próbnymi. Problem redukuje się więc do trafnego wyboru
funkcji próbnych Ni i znalezienia parametrów ai . Ogólnie należy się więc spodziewać, że ciąg
przyjętych aproksymacji
T '= T '(x,ai ) = ai Å" Ni dla i =1,2,...,n (4.4)
Rozdział 4 22
polepsza rozwiązanie wraz ze wzrostem liczby stosowanych funkcji próbnych Ni . Przyjmowane funkcje Ni
muszą być ciągłe i różniczkowalne do najwyższego rzędu występującego w całkowej formie równania.
Nie powinno nikogo zaskoczyć, że T ' , zastępujące w (4.1) funkcję T (x) , nie spełnia dokładnie
tegoż równania, to znaczy, że f (T ') `" 0 . Różnicę między rozwiązaniem dokładnym a przybliżonym
oznaczymy przez R(x,ai ) i zapiszemy:
f (T '(x, ai )) = R(x, ai ) (4.5)
Rezidua R zależą od x oraz ai . Z dokładnym rozwiązaniem mamy do czynienia wtedy, gdy R = 0 dla
wszystkich punktów obszaru &! . Dla rozwiązań przybliżonych R zasadniczo różni się od zera, jakkolwiek w
wybranych punktach obszaru &! warunek ten może być spełniony.
4.1.1. Metoda Ritza
Podstawowa metoda Ritza wymaga, by dla aproksymacji I rzędu był spełniony następujący warunek:
R(x,a1) = 0 . (4.6)
&!
Powyższe wyrażenie prowadzi do równania algebraicznego z niewiadomym parametrem a1 . Dla poprawnie
postawionego problemu z dobrze dobraną funkcją próbną rozwiązanie a1 zawsze istnieje. Za funkcje Ni (x)
przyjmowane mogą być wielomiany, funkcje cykliczne oraz inne funkcje ciągłe i różniczkowalne. Zilustrujmy
to, co dotąd powiedziano następującym przykładem.
Przykład. Rozwiążmy równanie różniczkowe zwyczajne o postaci:
2
d T
+1000Å" x2 = 0 dla 0 d" x d"1 (4.7)
dx2
przy zachowaniu następujących warunków brzegowych:
T (0) = 0 oraz T (1) = 0 (4.8)
Za funkcję próbną przyjmijmy funkcję nie zakłócającą problemu, spełniającą deklarowane warunki
brzegowe:
N1 = x Å" (1- x2) . (4.9)
23 Metody aproksymacyjnego rozwiązywania równań różniczkowych
Poszukiwane rozwiÄ…zanie przyjmujemy wiec w postaci:
T ' = a1 Å" N1(x) = a1 Å" x Å" (1- x2) . (4.10)
Podstawiając je do równania wyjściowego (4.7), rezidua wynoszą:
2 2
d T ' d T '
R = +1000Å" x2 , gdzie = -6Å" x Å" a1, (4.11)
dx2 dx2
tak więc w końcu
R(x,ai ) = -6Å" a1 Å" x +1000Å" x2 . (4.12)
Spełnienie wymagania (4.6) pociąga konieczność dobrania współczynnika a1 na podstawie równania:
1
(- 6 Å" a1 Å" x +1000 Å" x2)dx = 0 (4.13)
0
skąd otrzymujemy jedno równanie algebraiczne z niewiadomą a1 w postaci:
1
îÅ‚
ëÅ‚
x3
= 0 , (4.14)
ïÅ‚ìÅ‚- 3Å" a1 Å" x +1000 3
ìÅ‚
íÅ‚
ðÅ‚
0
1000
z którego wynika, że a1 = , wobec czego otrzymane rozwiązanie przybliżone ma następująca
9
formÄ™:
1000
T '= Å" x Å"(1- x2 ) . (4.15)
9
Warto w tym miejscu uzmysłowić sobie, że rozwiązanie równania różniczkowego przerodziło się
niepostrzeżenie z problemu analizy matematycznej w problem algebraiczny. Otrzymany wynik będzie
porównany z rozwiązaniem dokładnym. W tym miejscu prosi się dociekliwego Czytelnika, by zechciał
samodzielnie znalez
ć rozwiązanie nieskomplikowanego przecież równania (4.7).
Rozdział 4 24
Rys. 4.2 Interpretacja wydajności z
ródła ciepła i warunków brzegowych
Analizowane równanie Jest formą równania przewodnictwa cieplnego w izolowanym pręcie, z
wewnętrznym z
ródłem energii. Wydajność tego z
ródła jest proporcjonalna do x2 (rys.4.2).
Wyższych temperatur należy się spodziewać bliżej końca, dla x = 1, chociaż na obu brzegach powinien być
spełniony warunek T = 0 . Miejsce maksymalnej temperatury wypada dla x = 0.577 . Należy podkreślić, że
funkcje próbne, które w tym przypadku rozciągnięte są nad całym analizowanym obszarem, nie powinny
zaburzać opisywanego zjawiska fizycznego, ale powinny dokładnie spełniać warunki brzegowe. W
przyszłości technika MES wykaże, że można wyeliminować powyższe ograniczenia nakładane na funkcje
próbne, które to ograniczenia są trudne do spełnienia w większości praktycznych problemów technicznych.
4.1.2. Metoda wariacyjna Rayleigha-Ritza
Typowy problem jednowymiarowego rachunku wariacyjnego polega na znalezieniu takiej funkcji T(x) , by
zminimalizować bądz
zmaksymalizować całkę:
b
I = F(x,T(x),Tx (x))dx , (4.16)
a
dT
gdzie Tx = zaś F jest funkcjonałem (funkcją, której argumentami są funkcje). Można wykazać,
dx
że funkcjonał odpowiadający równaniu z poprzedniego przykładu ma postać:
25 Metody aproksymacyjnego rozwiązywania równań różniczkowych
2
1 dT
ëÅ‚
F = - ìÅ‚
+1000 Å" x2 Å"T (4.17)
2 dx
íÅ‚
Konsekwentne sformułowanie wariacyjne problemu z analizowanego przykładu wymaga więc ekstremalizacji
następującego wyrażenia:
2
1
ëÅ‚
ìÅ‚- 1 ëÅ‚ dT +1000 Å" x2 Å"T dx
I = (4.18)
ìÅ‚
ìÅ‚
2 dx
íÅ‚
0
íÅ‚
Idea polega więc na znalezieniu takiej funkcji T(x) , która minimalizuje I . Zanim pokażemy, w
jaki sposób powyższe równanie może być używane do znalezienia przybliżonego rozwiązania
analizowanego problemu, zauważmy, że:
- funkcja T(x) , która ekstremalizuje wyrażenie na I , jest funkcją spełniającą równanie różniczkowe i
zadane warunki brzegowe,
- wyjściowe równanie różniczkowe zawiera wyrażenie drugiego rzędu, a sformułowanie wariacyjne - tylko
pochodne pierwszego rzędu (jest to tzw. słabe sformułowanie).
Przykład (cd.). Przyjmijmy, podobnie jak poprzednio, funkcję próbną w postaci:
T' = a1 Å" N1(x) = a1 Å" x Å"(1- x2 ) ; (4.19)
podstawiając to wyrażenie do wzoru (4.18), otrzymujemy:
2
1
ëÅ‚
ìÅ‚- 1 ëÅ‚ dT' +1000 Å" x2 Å"T' dx .
I(a1) = (4.20)
ìÅ‚
ìÅ‚
2 dx
íÅ‚
0
íÅ‚
Zależność ta jest funkcją jedynej niewiadomej a1 . Stacjonarność będzie zapewniona przez spełnienie
warunku
dI
= 0 . (4.21)
da1
Wobec tego, że
dT'
= (-3Å" x2 +1) Å" a1 ,
da1
Rozdział 4 26
ostatecznie otrzymujemy (4.20) w jawnej postaci:
1
ëÅ‚- 1 2
I(a1) = Å" a1 Å"(1 - 3Å" x2) +1000 Å" x2 Å"(a1 Å" x Å"(1 - x2)) dx (4.20)
ìÅ‚
2
íÅ‚
0
a po wykonaniu przepisanych operacji
1000 2 dI 1000 4
I(a1) = Å" a1 - Å" a12 wiÄ™c także = - Å" a1 = 0
12 5 da1 12 5
5000
i w końcu poszukiwaną wartość a1 = .
48
Tym razem otrzymane rozwiązanie jest następujące:
5000
T'(x) = Å" x Å" (1 - x2 ) . (4.22)
48
Porównania dotąd otrzymanych rozwiązań dokonano na rysunku 4.3.
Rys. 4.3. Porównanie rozwiązań m etodą Ritza i m etodą wariacyjną
27 Metody aproksymacyjnego rozwiązywania równań różniczkowych
4.2. Metoda ważonych reziduów
Metoda ważonych reziduów jest silnym narzędziem znajdowania przybliżonych rozwiązań równań
różniczkowych powszechnie stosowanych w problemach inżynierskich. Poniżej prezentujemy cztery
najbardziej popularne wersje tej metody. Pewną zaletą ważonych reziduów jest możliwość ominięcia
sformułowań wariacyjnych, które częstokroć mogą sprawiać wiele trudności. Omówimy więc po krotce,
posiłkując się ciągle tym samym przykładem, dwie wersje metody kollokacyjnej, znaną metodę
najmniejszych kwadratów oraz tzw. metodę Galerkina i bezpośrednio z niej wywodzącą się MES. Na
koniec porównamy wyniki z poprzedniego przykładu uzyskane przy zastosowaniu różnych metod.
Przedstawiając koncepcję ogólną, ograniczymy się, tak jak to miało miejsce poprzednio,
do pojedynczego równania z jedną tylko zmienną niezależną
f (T(x)) = 0 , działającego w obszarze &! (4.23)
z określonymi warunkami brzegowymi:
g1(T(x)) = 0 na obszarze “1 ,
g2(T(x)) = 0 na obszarze “2 , (4.24)
Podobnie przybliżymy rozwiązanie za pomocą funkcji próbnych Ni (x) i nieznanych współczynników ai :
n
T' = T'(x,a1, a2,...,an ) = Å" Ni (x) (4.25)
"a
i
i=1
i otrzymamy spełnienie równania wyjściowego z dokładnością do reziduów R(x,ai ) :
f (T'(x, ai )) = R(x,ai ) (4.26)
Metody ważonych reziduów zakładają wyznaczenie parametrów ai przez spełnienie określonego warunku:
w(x) Å" Ri (x, a ) = 0 , dla i = 1,2,...,n , (4.27)
&!
gdzie funkcje w(x) są tzw. funkcjami wagowymi. Wybór tych funkcji różnicuje wersje metody ważonych
reziduów. Zauważmy, że przyjmując w zadaniu z jedną tylko niewiadomą stalą funkcję wagową w(x) = 1 ,
otrzymamy znaną już nam wersję metody Ritza.
Rozdział 4 28
4.2.2. Punkt kollokacji
W tej wersji metody ważonych reziduów za funkcję wagi przyjmuje się wyrażenie
wi (x) = ´ Å" (x - xi ) , (4.28)
gdzie ´ peÅ‚ni funkcjÄ™ delty Kroneckera, która jest zdefiniowana nastÄ™pujÄ…co:
b
´ Å"(xi - x)dx =1, dla xi "(a,b) ,
a
b
´ Å"(xi - x)dx = 0 , dla xi " (a,b) , (4.29)
a
gdzie współrzędna xi opisuje położenie punktu kollokacji i wybierana jest na etapie formułowania
ograniczenia. Przy stosowaniu n funkcji próbnych wymaga się - w celu wyznaczenia wszystkich
współczynników ai - zastosowania n punktów kollokacji. Ograniczenia wyglądają wówczas
następująco:
´ Å"(x - xi ) Å" R(x,ai ) = 0 , dla i =1,2,...,n , (4.30)
&!
Formułując je w n punktach xi , otrzymujemy układ równań algebraicznych, z którego
wyznaczymy niewiadome ai .
Przykład (cd.). Wracając do naszego przykładu, wielkości reziduów są zdefiniowane
jak poprzednio: R(x,ai ) = -6Å" a1 Å" x +1000Å" x2 . Przy przyjÄ™tej tylko jednej funkcji próbnej
1
wymaga się spełnienia ograniczenia typu (4.30) w jednym tylko punkcie, na przykład x1 = .
2
Wówczas otrzymujemy:
2
1 1 1000
- 6Å" a1 Å" +1000Å"ëÅ‚ = 0 , czyli a1 = (4.31)
ìÅ‚
2 2 12
íÅ‚
Rozwiązanie końcowe jest tym razem w postaci:
1000
T '(x) = Å" x Å"(1- x2 ) . (4.32)
12
29 Metody aproksymacyjnego rozwiązywania równań różniczkowych
4.2.2. Podobszar kollokacji
Funkcje wagi w tej metodzie dobiera się w taki sposób, że ich wartości są równe 1 w danej części obszaru
&!i , natomiast na pozostałej części obszaru &! są równe zeru. Podobszarów kollokacji definiuje się tyle ile
jest przyjętych funkcji próbnych. Sytuację ilustruje rysunek 4.4, który przykładowo dzieli cały obszar na trzy
podobszary. Matematycznie funkcje wag w poszczególnych podobszarach wyrażone są następująco:
1 x "&!1
Å„Å‚
w1(x) =òÅ‚
ół0 x "&!1,
1 x "&!2
Å„Å‚
w2 (x) =òÅ‚
ół0 x "&!2 ,
1 x "&!3
Å„Å‚
w3(x) =òÅ‚
ół0 x "&!3
Rys. 4.4. Przyjęte funkcje wagowe w podobszarze kollokacj
Podobszary &!i nie zawierają wspólnych elementów, lecz w sumie wypełniają cały obszar
&! .Te założenia w konsekwencji prowadzą do następującego układu równań, którego rozwiązanie
pozwala wyznaczyć niewiadome współczynniki ai :
Rozdział 4 30
R(x, ai ) = 0 , dla i = 1,2,...,n , (4.33)
&!i
Przykład (cd.). W naszym przypadku przyjęcie jednej funkcji próbnej wymaga potraktowania całego obszaru
&! jako jednego tylko podobszaru kollokacji &! = &!1 . Tak więc &!1 = 0,1 , a całość sprowadza się do
spełnienia warunku identycznego z warunkiem Ritza:
R(x, a1) = 0 , (4.34)
&!
skąd także
1
1000
(- 6 Å" a1 Å" x +1000 Å" x2)dx = 0 oraz a1 = (4.35)
9
0
i końcowy wynik jest identyczny z tym, który otrzymaliśmy poprzednio stosując propozycję Ritza:
1000
T'(x) = Å" x Å" (1 - x2 ) (4.36)
9
4.2.3. Metoda najm niejszych kwadratów
Metoda najmniejszych kwadratów jest powszechnie znana wszystkim tym, którzy stają przed problemem
opracowania wyników eksperymentu. Polega ona na zminimalizowaniu wyrażenia typu:
2
(R(x, ai )) dx = 0 , (4.37)
&!
co oznacza, że cząstkowe pochodne powyższej całki względem kolejnych współczynników ai mają być
równe zeru. Otrzymujemy więc układ równań algebraicznych o postaci:
"I " " "R
2 2
= (R) dx = (R(x, ai )) dx = R dx = 0,
"a1 "a1 &! "a1 "a1
&! &!
"I " " "R
2 2
= (R) dx = (R(x, ai )) dx = R dx = 0,
"a2 "a2 &! "a2 "a2
(4.38)
&! &!
"I " " "R
2 2
= (R) dx = (R(x,ai )) dx = R dx = 0.
"an "an &! "an "an
&! &!
31 Metody aproksymacyjnego rozwiązywania równań różniczkowych
Widzimy wiec, że role funkcji wagowych w powyższych równaniach pełnią pochodne reziduów, a więc
"R
wi (x) = dla i = 1,2,...,n .
"ai
Przykład (cd.). Wracając znów do analizowanego przykładu, mamy funkcję reziduów w znanej już nam
dobrze postaci (4.12), funkcja zaś wagowa tym razem określona jest jako
"R
= -6 Å" x , (4.39)
"ai
i w zwiÄ…zku z czym podstawiajÄ…c do (4.27), otrzymujemy
1
1
- 6 Å" x Å"(- 6 Å" a1 Å" x +1000 Å" x2)dx = [12 Å" x3 Å" a1 -1500 Å" x4] = 0 ,
0
0
skÄ…d
1000
a1 = (4.40)
48
RozwiÄ…zanie otrzymane tym sposobem wynosi:
1000
T'(x) = Å" x Å" (1 - x2 ) . (4.41)
8
4.2.4. Metoda Galerkina
Ostatnią z proponowanych tu metod opierających się na ważonych reziduach jest metoda
zaproponowana przez Galerkina. Polega ona na przyjęciu, że rolę funkcji wagowej pełni przyjmowana
funkcja próbna. Zakładamy więc, że wi (x) = Ni (x) i żądamy - w przypadku używania kilku funkcji
próbnych jednocześnie - spełnienia warunków:
Ni (x) Å" R(x, a1, a2,..., an )dx = 0 , dla i =1,2,...,n , (4.42)
&!
Przykład (cd.). Ponieważ, jak dotąd, zawsze przyjmowaliśmy funkcję próbną w postaci
N1(x) = x Å"(1- x2 ) , wymagamy wiÄ™c teraz speÅ‚nienia warunku.
Rozdział 4 32
1
x Å"(1 - x2)Å"(- 6 Å" a1 Å" x +1000 Å" x2)dx = 0 (4.43)
0
5000
który prowadzi do wielkości współczynnika a1 = . Otrzymana tak funkcja T '(x) ma
48
postać:
5000
T '(x) = Å" x Å"(1- x2 ) (4.44)
48
Tym razem wynik pokrywa siÄ™ z wynikiem otrzymanym w wariacyjnej metodzie Rayleigha-Ritza.
Spróbujmy na koniec tych rozważań porównać otrzymane wyniki z prostym do uzyskania
rozwiązaniem dokładnym. Postać tego rozwiązania jest następująca:
1000
T (x) = Å" x Å"(1- x2 ) (4.45)
12
Zwróćm y uwagę, że potęga zm iennej w nawiasie jest wyższa od tej, którą przyjmowaliśm y w funkcji
próbnej. M ożna by powiedzieć, że nasze poszukiwania skazane były na częściowe tylko powodzenie,
gdyż posługiwaliśm y się tylko jedną funkcją próbną, która odbiegała postacią od rozwiązania
dokładnego. Trzeba w tym m iejscu podkreślić raz jeszcze, że zwiększając liczbę wyrazów wielom ianu,
Rys. 4.5. Porównanie rozwiązań otrzym anych różnym i m etodam i
33 Metody aproksymacyjnego rozwiązywania równań różniczkowych
a tym samym zwiększając liczbę niewiadomych współczynników, nasze rozwiązania zbliżałyby się do
dokÅ‚adnego. GdybyÅ›my od razu przyjÄ™li funkcjÄ™ próbnÄ… w postaci N1(x) = x Å"(1- x3) , każda z
wymienionych tu metod doprowadziłaby nas do rozwiązania dokładnego. Celowo nie uczyniliśmy tego, gdyż
w przypadkach ważnych z punktu widzenia zastosowań praktycznych nie sposób jest przewidzieć postać
równania, bo nie są znane rozwiązania analityczne.
Rysunek 4.5 przedstawia krzywe naszego rozwiązania. Linią ciągłą zaznaczono rozwiązanie
dokładne. Widzimy, że w analizowanym przypadku metody najmniejszych kwadratów i podobszarów
kollokacji prowadziły do określenia temperatur znacznie wyższych niż dokładne.
Podkreślmy raz jeszcze, co jest niezwykle ważne przy analizowanej aproksymacji, że przyjmowane
funkcje próbne rozpięte są nad całym obszarem &! . Dalej, przyjmując koncepcję podziału całego obszaru
na części (elementy), będziemy omawiać aproksymację tylko na fragmencie obszaru, co w znacznym
stopniu ułatwi nam zadanie. Warto w tym miejscu wspomnieć, że w przypadku choć trochę skomplikowanych
warunków brzegowych nie można zaproponować funkcji próbnych spełniających wymieniane powyżej
wymagania co do dokładnego spełnienia warunków brzegowych i niezakłócania problemu fizycznego.
4.2.5. Przykład MES w aproksym acji Galerkina
Zasadnicze nowe elementy, które chcemy teraz wprowadzić, polegać będą na przyjęciu nawet
prostszej, bo liniowej funkcji próbnej, ale przykładanej tylko lokalnie na pewnym podobszarze, który nazywać
będziemy elementem skończonym.
Załóżmy więc, że obszar zawarty między punktami x = a i x = b , (a d" x d" b) został podzielony na
M podobszarów (elementów), tak jak to pokazano na rysunku 4.6. Proces ten nazywamy dyskretyzacją.
Przeanalizujmy zatem element o dwóch węzłach w punktach xk oraz x . Dowolną wielkość, która na tym
j
obszarze ma podlegać zmianom (może to być np. przemieszczenie lub jak w naszym
przykładzie temperatura) oznaczmy przez y(x) . Niech na obszarze elementu zmiana tej wielkości będzie
opisana za pomocą zależności
Rys. 4.6. Podział obszaru jednowym iarowego na elem enty.
Rozdział 4 34
ye(x) = m Å" x + b (4.46)
gdzie e oznacza, że aproksymacja dotyczy wyłącznie elementu, a wartości współczynników m
i b wyznaczymy z następujących warunków:
dla x = x ye(x ) = y ,
j j j
dla x = xk ye(xk ) = yk . (4.47)
Podstawiając te warunki do równania (4.46), po rozwiązaniu układu równań z
niewiadomymi m i b otrzymujemy:
ëÅ‚ - y
yk - y yk j
j
ìÅ‚
m = , b = y - Å" x (4.48)
j j
ìÅ‚
xk - x xk - x
j j
íÅ‚
i funkcję opisującą zmianę wielkości ye(x) na obszarze elementu w postaci:
ëÅ‚ ëÅ‚ - x
x
xk - x
j
ìÅ‚ ìÅ‚
ye(x) = Å" y + Å" yk , (4.49)
j
ìÅ‚ ìÅ‚
xk - x xk - x
j j
íÅ‚ íÅ‚
przy czym wielkości y i yk pełnią rolę podobną do tej, którą pełniły współczynniki ai w
j
poprzednich rozważaniach. Powyższą zależność można więc zapisać jako
ye(x) = N (x) Å" a + Nk (x) Å" ak . (4.50)
j j
Przyjmowane więc funkcje próbne (kształtu) są tym razem funkcjami liniowymi i wyrażają się
następującymi zależnościami:
x - x
xk - x
j
N (x) = , Nk (x) = (4.51)
j
xk - x xk - x
j j
35 Metody aproksymacyjnego rozwiązywania równań różniczkowych
Rys. 4.7. Liniowe funkcje próbne w obszarze elem entu
Zmienność tych funkcji pokazano na rysunku 4.7. Ogólnie można więc zapisać zmienność pola
ye(x) stosujÄ…c notacjÄ™ macierzowÄ…:
ye(x) = N Å" ae , (4.52)
gdzie dwuelementowa macierz N gromadzi funkcje kształtu, a wektor ae zawiera wartości zmiennych
pola punktach węzłowych.
Przykład (cd.). Powróćmy do naszego przykładu. Pamiętamy, że w metodzie Galerkina funkcjami wagi
e
były funkcje próbne N(x) . Przyjmijmy więc teraz równanie opisujące zmianę temperatury T (x) dla
pojedynczego elementu skończonego w przedziale od x do xk .
j
e
T (x) = N (x) Å"Tj + Nk (x) Å"Tk , (4.53)
j
gdzie Tj i Tk odpowiadają teraz zmiennym ai w metodzie Galerkina, a funkcje kształtu mają
znaną liniową postać. Równania Galerkina można zapisać w postaci:
xk
M
N (x)(Re(x,Tj ,Tk ))dx = 0 ,
" j
e=1
x
j
(4.54)
xk
M
Nk (x)(Re(x,Tj ,Tk ))dx = 0 ,
"
e=1
x
j
Rozdział 4 36
przy czym M odpowiada liczbie elementów przyjętych w dyskretyzacji. Układ (4.54) można
również zapisać macierzowo
xk
M
T
N (x)(Re(x,Tj ,Tk ))dx = 0 . (4.55)
"
e=1
x
j
W naszym przykładzie rezidua Re dla typowego elem entu wyrażają się znaną zależnością:
2 e
d T
Re = +1000 Å" x2 . (4.56)
dx2
Otrzymujemy więc:
xk
2 e
M
ëÅ‚ d T
T
ìÅ‚
N (x)ìÅ‚ +1000 Å" x2 dx = 0 . (4.57)
"
dx2
e=1 íÅ‚
x
j
Na tym etapie pomińmy sumowanie i poprzestańmy na analizie jednego tylko elementu. Przedstawmy
powyższe równanie w postaci dwóch całek:
xk xk
2 e
d T
T T
N (x) Å" dx + N (x) Å"(1000 Å" x2)dx = 0 ; (4.58)
dx2
x x
j j
całkując pierwsze wyrażenie przez części, otrzymujemy:
xk xk T e xk
e
îÅ‚ dN dT
dT
T T
(4.59)
ïÅ‚N Å" dx x - dx Å" dx dx + N Å"(1000 Å" x2)dx = 0 ,
ðÅ‚
x x
j j
j
e
ponieważ przyjÄ™to w aproksymacji, że T = N Å" ae , skÄ…d:
xk xk T
xk
e
îÅ‚ dN dN
dT
T T
ïÅ‚N Å" dx x - dx Å" dx Å" aedx + N Å"(1000 Å" x2)dx = 0 , (4.60)
ðÅ‚
x x
j j
j
xk T
dN dN
e
gdzie wyrażenie Å" dx nazywane jest macierzÄ… sztywnoÅ›ci elementu i oznaczane przez K . Ze
dx dx
x
j
względu na to ,że wektor ae , opisujący wartości temperatur w węzłach, nie zależy od zmiennej x ,
powyższe równanie możemy zapisać macierzowo w postaci:
37 Metody aproksymacyjnego rozwiązywania równań różniczkowych
e e
K Å" ae = f
. (4.61)
e
Macierz sztywności K o wymiarach 2x2 jest macierzą symetryczną i osobliwą. Pozostałe wektory ae i
e
f mają wymiary 2x1. Wyniki działań we wzorze (4.60) przedstawiają się następująco:
1
1 îÅ‚ -1
e
K = Å"ïÅ‚
xk - x
j ðÅ‚-1 1
xk
dT(x )
dT dT(xk )
îÅ‚ = j
T T T
ïÅ‚dx Å" N x dx Å" N (xk ) - dx Å" N (x j ) , (4.62)
ðÅ‚
j
0 1
îÅ‚ îÅ‚
T T
gdzie N (xk ) =ïÅ‚ i N (x ) =ïÅ‚
j
ðÅ‚1 ðÅ‚0
oraz
1000
xk 4 4
îÅ‚ - 4 Å" xk Å" x + 3Å" x
xk 4
T j j
12
N Å"(1000 Å" x2)dx = . (4.63)
ïÅ‚
4
xk - x xk 4 - 4 Å" x Å" xk 4 + x
ïÅ‚3Å"
j j j
x
ðÅ‚
j
Rys. 4.8. Dyskretyzacja MES przykładu w aproksym acji Galerkina (podział na 5 elem entów)
Dokonajmy teraz dyskretyzacji rozważanego przedziału na pięć równych części (elementów) o długości 0.2
każdy. W wyniku tego podziału otrzymamy obraz, jaki pokazano na rysunku 4.8. Niewiadomymi są teraz
temperatury w sześciu węzłach. Zgodnie z (4.62) macierze sztywności elementów są identyczne i wynoszą:
1
1 îÅ‚ -1 5 - 5
îÅ‚
1 îÅ‚ -1 1
e
K = Å"ïÅ‚ = Å"ïÅ‚ =
xk - x
ðÅ‚
j ðÅ‚-1 1 0.2 - 0.0 ðÅ‚-1 1 ïÅ‚- 5 5
2 3 4 5
K1 = K = K = K = K . (4.64)
Rozdział 4 38
Wektory prawych stron równania (4.61) dla elementów otrzymamy w postaci:
1000
0
îÅ‚ îÅ‚-1 îÅ‚ - 4 Å" 0.2 Å" 0.04 + 3Å" 0.04
dT(x2 ) 0.24
dT(x1)
1
12
f =ïÅ‚ Å" +ïÅ‚ Å" + ,
dx dx 0.2 - 0.0ïÅ‚ - 4 Å" 0.0 Å" 0.24 + 0.04
ðÅ‚1 ðÅ‚0
ðÅ‚3Å" 0.24
0
îÅ‚ îÅ‚-1 0.667
dT(x2 ) îÅ‚
dT(x1)
1
f =ïÅ‚ Å" +ïÅ‚ Å" +ïÅ‚ ,
dx dx
ðÅ‚1 ðÅ‚0 ðÅ‚2.000
0
îÅ‚ îÅ‚-1 7.300
dT(x3) îÅ‚
dT(x2 )
2
f =ïÅ‚ Å" +ïÅ‚ Å" +ïÅ‚ ,
dx dx
ðÅ‚1 ðÅ‚0 ðÅ‚11.30
0
îÅ‚ îÅ‚-1 22.00
dT(x3)
dT(x4 ) îÅ‚
3
f =ïÅ‚ Å" +ïÅ‚ Å" +ïÅ‚ , (4.65)
dx dx
ðÅ‚1 ðÅ‚0 ðÅ‚28.60
0
îÅ‚ îÅ‚-1 44.70
dT(x5) îÅ‚
dT(x4 )
4
f =ïÅ‚ Å" +ïÅ‚ Å" +ïÅ‚ ,
dx dx
ðÅ‚1 ðÅ‚0 ðÅ‚54.00
0
îÅ‚ îÅ‚-1 75.40
dT(x6 ) îÅ‚
dT(x5)
5
f =ïÅ‚ Å" +ïÅ‚ Å" +ïÅ‚ .
dx dx
ðÅ‚1 ðÅ‚0 ðÅ‚87.30
Podobnie jak to czyniliśmy poprzednio, dokonujemy agregacji macierzy sztywności całego
układu oraz wektora prawych stron zależności
K Å" a = f (4.66)
gdzie teraz macierz K oraz wektory a i f opisują stan układu odniesiony do
współrzędnych globalnych. Po scaleniu otrzymujemy następujący układ równań:
5
îÅ‚1 îÅ‚
îÅ‚ - 5 0 0 0 0 T 0.667 - dT(x1) / dx
ïÅ‚- 5 10 - 5 0 0 0 ïÅ‚T2 ïÅ‚ 9.33
ïÅ‚ ïÅ‚ ïÅ‚
ïÅ‚0 - 5 10 - 5 0 0 ïÅ‚T3 ïÅ‚ 33.33
Å" = (4.67)
ïÅ‚0 0 - 5 10 - 5 0 ïÅ‚T4 ïÅ‚ 73.33 .
ïÅ‚ ïÅ‚ ïÅ‚
ïÅ‚5 ïÅ‚
ïÅ‚
0 0 0 - 5 10 - 5 T 129.40
ïÅ‚ ïÅ‚
ïÅ‚
ðÅ‚0 0 0 0 - 5 10 ðÅ‚T6 ðÅ‚87.3 + dT(x6 ) / dx
Wprowadzając warunki brzegowe (w węzłach 1 i 6 temperatura jest równa zeru), czyli na przykład
wykreślając z układu równań pierwszy i ostatni wiersz oraz pierwszą i ostatnią kolumnę z macierzy sztywności,
otrzymujemy układ równań z czterema niewiadomymi. Rozwiązaniem tego układu są następujące wyniki:
T1 = 0,T2 = 16.5,T3 = 31.2,T4 = 39.2,T5 = 32.5,T6 = 0.
39 Metody aproksymacyjnego rozwiązywania równań różniczkowych
A
atwo sprawdzić, że otrzymane wyniki prawie w ogóle nie odbiegają od rozwiązania analitycznego
oczywiście co do wartości temperatur w wybranych przez dyskretyzację punktach. Warto zauważyć, że
zabieg dyskretyzacji mimo przyjęcia funkcji próbnych w prostej postaci funkcji liniowych prowadzi do
osiągnięcia bardzo dużej dokładności aproksymacji rozwiązania analitycznego.
Zadania
1. Przeanalizuj problem (4.7) z warunkami (4.8) wszystkimi przedstawionymi w tym
rozdziale m etodam i, przyjm ując funkcje próbne w postaci:
a) N1(x) = x Å"(1 - x4),
b) N1(x) = sin(Ä„x)
2. Znajdz
rozwiązanie problem u (4.7) m etodą Ritza przyjm ując funkcję próbną w
postaci:
N1(x) = x Å"(1 - x3).
W jaki sposób rozwiązanie to odpowiada rozwiązaniu dokładnemu? Uogólnij otrzymany
wynik.
3. Przyjmij równanie różniczkowe w postaci:
2
d y
+ 6 Å" y = 10 Å" x dla 0 d" x d" 2 ,
dx2
przy spełnieniu warunków brzegowych y(0) = 1 i y(2) = 0 . Sprawdz
, że funkcjonał, który
podlega ekstremalizacji w metodzie Rayleigha-Ritza jest wyrażony w postaci:
2
2
îÅ‚
1 dy
ëÅ‚
I = y -10 Å" x Å" y - ìÅ‚
dx .
ïÅ‚3Å"
2 dx
íÅ‚
ïÅ‚
0
ðÅ‚
4. Przeanalizuj problem (4.7) metodą punktów kollokacji, przyjmując dwie funkcje próbne
N1(x) = x Å"(1 - x2),
N1(x) = x Å"(1 - x4),
1 2
oraz dwa punkty kollokacji dla współrzędnych x = i x = .
3 3
Rozdział 4 40
5. Rozwiąż problem z zadania 4 metodą najmniejszych kwadratów.
6. Przeanalizuj problem (4.7) metodą punktu kollokacji, zakładając
x = 0.2, x = 0.4, x = 0,6, x = 0.8 . SporzÄ…dz
wykresy i porównaj otrzymane rozwiązania z
wynikiem dokładnym.
7. Oblicz ugięcie belki swobodnie podpartej o rozpiętości l obciążonej ciężarem
równomiernie rozłożonym o intensywności q przyjmując funkcje próbne:
N1(x) = x Å"(1 - x),
Ä„ Å" x
ëÅ‚1
N1(x) = x Å" - sin ,
ìÅ‚
2 Å" l
íÅ‚
a)metodą punktu kollokacji dla x = 0.5 i podobszaru kollokacji przyjętego dla całego
przedziału,
b) metodą najmniejszych kwadratów.
Porównaj otrzymane ugięcie w środku rozpiętości z wynikiem dokładnym .