1. Metoda elementu skończonego. Metoda Ritza.
1.1. Sformułowanie problemu. Dwu-wymiarowe przewodzenie ciepła w długim pręcie..
Rozważmy przekrój kwadratowy z brzegiem długiego pręta pokazany na Rys. 1 z nałożonymi warunkami brzegowymi,
Typu Dirichleta na części :
T = 50 dla y = 0 (1)
T = 100 dla y = 2
i typu Neumanna na części :
= 0 dla x = 0
= 0 dla x = 2, gdzie . (2)
Przewodzenie ciepła w przekroju dwuwymiarowym ze stałym współczynnikiem przewodzenia jest opisane równaniem Laplace’a
na . (3)
Jest to szczególny przypadek równania zachowania energii
gdzie
jest pochodną wędrowną temperatury. Q jest mocą cieplną na jednostkę objętości wynikającą ze źródeł zewnętrznych
Rys.1.Przekrój poprzeczny pręta z wymiarami i warunkami brzegowymi.
Chcemy rozwiązać ten problem numerycznie w sposób przybliżony przy użyciu tzw. metody wariacyjnej Ritza. W metodzie tej zastępujemy problem sformułowany za pomocą równań (1-3) przez równoważny problem minimalizacji odpowiednio zdefiniowanej całki, zwanej funkcjonałem.
Można udowodnić, że funkcja spełniająca warunki (1-3) minimalizuje następujący funkcjonał
(4)
w zbiorze (klasie) F tzw. funkcji próbnych zdefiniowanych na i spełniających warunek Dirichleta (1). Warunki Neumana (2) będą już automatycznie spełnione (dzięki odpowiedniej formie (4) funkcjonału).
Możemy otrzymać przybliżenie rozwiązania oryginalnego problemu opisanego równaniami (1-3) jeżeli ograniczymy klasę F , która jest zbyt obszerna, do klasy F’ (o skończonym wymiarze) funkcji ciągłych, które są tylko kawałkami różniczkowalne, ale zdolne do aproksymacji każdej funkcji ze zbioru F zdefiniowanych na i spełniających warunek Dirichleta (1). Oznacza to, że optymalizacja (minimalizacja) będzie przeprowadzona na węższym zbiorze funkcji. Warunki Neumana (2) będą tak jak w szerszym zbiorze automatycznie spełnione poprzez formę (4) funkcjonału.
Klasę F’ tworzymy następująco:
1) Dziedzinę dzielimy na skończoną liczbę figur geometrycznych zwanych elementami (mogą to być trójkąty, czworokąty dla 2D lub czworościany, sześcio-ściany dla 3D) tworzącymi w ten sposób siatkę. Wyróżnione punkty każdej z tych figur są nazywane punktami węzłowymi i powinno ich być tyle aby w jednoznaczny sposób definiowały interpolację wielomianową wewnątrz tego elementu.
2) Złożenie (sklejenie) w sposób ciągły tych interpolacji na elementach daje interpolację na całej dziedzinie . Zbiór wszystkich takich interpolacji daje zbiór na którym dokonana jest minimalizacja. (Minimum na jest oczywiście nie mniejsze niż minimum na zbiorze wszystkich funkcji próbnych .)
Rys. 2. Przykład siatki obliczeniowej.
Tak więc wartość funkcjonału dla dowolnej funkcji ze zbioru może być rozłożona na sumę udziałów elementowych
(5)
1.2. Tworzenie funkcji interpolacyjnych.
Rozważmy dowolny element trójkątny ‘e’ jak z Rys. 2. Dla trzech punktów węzłowych będących jego wierzchołkami jednoznaczną interpolacją jest interpolacja liniowa
(6)
gdzie , l=1,2,3 , są stałymi współczynnikami, innymi dla każdego elementu.
W celu zapewnienia ciągłości po sklejeniu tych interpolacji liniowych wyrazimy te stałe współczynniki poprzez wartości funkcji interpolacyjnej w punktach węzłowych , , gdzie i , j , k są numerami tych punktów. W tym celu wystarczy rozwiązać układ równań otrzymany przez podstawienie tych wartości trzykrotnie do równania (6) oraz odpowiadających par współrzędnych (xi , yi) , (xj , yj) , (xk , yk) tych punktów węzłowych
Układ ma jednoznaczne rozwiązanie jeśli wyznacznik macierzy współczynników jest różny od zera
W wyniku zastosowania wzorów Cramera otrzymujemy
(7)
, , (8)
i cykliczna zmiana indeksów powinna być zastosowana dla pozostałych.
Po podstawieniu wyrażeń na współczynniki do równania (6) otrzymujemy
(9)
(10)
Powyższe wzory definiują tzw. funkcje bazowe (czasem zwane funkcjami kształtu), które są wielomianami, charakterystycznymi dla obranej interpolacji (tu liniowymi). Posiadają własność
(11)
Używając tej własności, można rozszerzyć interpolację elementową na interpolację globalną na całym zbiorze według analogicznego wzoru jak w równaniu (9)
(12)
i sumowanie ma być wykonane po wszystkich punktach węzłowych siatki. Oczywiście, dla ustalonego punktu suma (12) redukuje się do wzoru (9) odpowiadającego elementowi zawierającemu ten punkt.
1.3. Minimalizacja funkcjonału.
Wymagane pochodne występujące pod całkami we wzorze (5) mogą być teraz łatwo obliczone
(13)
i po podstawieniu do całki reprezentującej udział elementowy otrzymamy
(14)
ze stałą funkcją podcałkową, która może być wyłączona poza znak całki. Po uwzględnieniu oczywistego wzoru
otrzymamy
(15)
Powyższe wyrażenie jest słuszne dla dowolnego elementu w siatce z tą różnicą, że inne trójki wskaźników a zatem i inne trójki wartości węzłowych , , występują dla różnych elementów.
Wyrażenie (15) jest kwadratowe względem wartości węzłowych , , (nie znanych jeszcze na tym etapie). Potrzebujemy układu równań dających jednoznaczne rozwiązanie dla tych niewiadomych.
Po zsumowania wszystkich udziałów elementowych zgodnie ze wzorze (5) otrzymujemy funkcję kwadratową wszystkich wartości węzłowych, o postaci ogólnej
(16)
Dla funkcji (w tym przypadku wielu zmiennych) możemy skorzystać z warunku koniecznego istnienia minimum.
Warunek wymaga by wszystkie pochodne cząstkowe funkcji względem wszystkich zmiennych, czyli względem wszystkich wartości węzłowych, były równe zeru
m = 1, ... ,Np (17)
Warto zauważyć, że pochodne cząstkowe wartości elementowych pod znakiem sumy są różne od zera jedynie względem tych wartości węzłowych, które należą do tego elementu.
Tak więc pochodne wartości elementowych można otrzymać przez różniczkowanie zależności (15)
(18)
co daje wyrażenia liniowe po prawej stronie względem wartości węzłowych.
Dzięki temu, to samo można przedstawić w postaci macierzowej
(19)
(20)
dla dowolnych indeksów i , j = 1, 2, 3.
Dzięki przedstawieniu macierzowemu liniowej postaci pochodnych cząstkowych, sumę w równaniach (17) można zrealizować równoważnie przez sumowanie macierzy będących udziałami elementowymi tego warunku minimalizacji.
Sumowanie, zwane procesem asemblacji (składania), polega na dodawaniu małych macierzy elementowych, takich jak w równaniu (19), do jednej wspólnej, dużej macierzy układu (17). W macierzy tej każdy wiersz i każda kolumna odpowiada jednemu numerowi węzła w siatce.
Dodawanie udziałów elementowych odbywa się więc tylko w określone miejsca odpowiadające numerom węzłów należących do tego asemblowanego elementu. Ten algorytm składania przedstawiony jest poniżej, w wyniku czego otrzymuje się układ równań liniowych o niewiadomych
(21)
1.4. Przykład numeryczny.
W każdym problemie numerycznym związanym z metoda elementu skończonego muszą być zdefiniowane następujące zbiory danych (tutaj odpowiadają one problemowi przewodnictwa ciepła w pręcie z siatką przedstawioną na Rysunku 2).
Tabela 1. Tabela połączeń elementów
Numer
elementu
Globalne numbery węzłów
1 2 3
i j k
1
4
5
2
3
6
...
viktus