METODA ELEMENTU SKOC
CZONEGO 
WPROWADZENIE
Finite Element Method (FEM)
Magdalena Gierszewska"
20 stycznia 2005
"
http://fatcat.ftj.agh.edu.pl/<" dziopa/
1
Spis treści
1 Wprowadzenie 3
2 Podstawy metody FEM 3
2.1 Metoda elementu skończonego . . . . . . . . . . . . . . . . . 3
2.2 Algorytm - Etapy rozwiÄ…zania problemu za pomocÄ… FEM . . 4
2.3 Definicja elementu skończonego . . . . . . . . . . . . . . . . . 7
3 Siatka elementów skończonych 8
3.1 Numeracja elementów i węzłów . . . . . . . . . . . . . . . . . 8
3.2 Tablica koneksji . . . . . . . . . . . . . . . . . . . . . . . . . 8
4 Zastosowania metody elementu skończonego 10
4.1 Rozwiązywanie równań różniczkowych . . . . . . . . . . . . . 10
4.2 Modelowanie równoległe z wykorzystaniem Metody Elemen-
tów Skończonych. . . . . . . . . . . . . . . . . . . . . . . . . . 10
2
1 Wprowadzenie
Metoda elementów skończonych ( z ang. Finite Element Method) jest to
obecnie jedna z najszerzej stosowanych metod rozwiązywania różnych pro-
blemów inżynierskich. Jej uniwersalność, polegająca na łatwości schematy-
zacji różnych obszarów o skomplikowanej geometrii, także niejednorodnych i
anizotropowych, kwalifikuje ją jako dobre narzędzie do modelowania proble-
mów fizycznych. Rozwój metody elementów skończonych przebiegał i nadal
przebiega, równolegle z rozwojem techniki komputerowej. Pierwsze prace sto-
sujące metodę elementów skończonych zostały opublikowane w latach czter-
dziestych ubiegłego wieku. Początkowo obliczenia przeprowadzane za po-
mocą metody elementów skończonych dotyczyły obiektów o bardzo prostych
geometriach (najczęściej modelowanych jako jednowymiarowe) i stałych wła-
snościach materiałowych oraz zjawisk opisanych liniowymi równaniami róż-
niczkowymi. Od lat siedemdziesiątych metodę elementów skończonych za-
częto stopniowo stosować do rozwiązywania problemów nieliniowych, ale da-
lej dla obiektów o stosunkowo prostych geometriach, modelowanych jako 1D
lub 2D. Gwałtowny rozwój techniki komputerowej w latach osiemdziesiątych,
związany z coraz większą mocą obliczeniową komputerów oraz możliwością
operowania i przechowywania bardzo dużych zbiorów informacji, umożliwił
zastosowanie metody elementów skończonych do obliczeń problemów nieli-
niowych dla obiektów o dowolnie złożonych geometriach, szczególnie 3D.
2 Podstawy metody FEM
2.1 Metoda elementu skończonego
Zajmijmy się przestrzenią liniową F, w której jest zdefiniowany iloczyn ska-
larny | . Elementami tej przestrzeni są funkcje z określonym działaniem
dodawania funkcji i mnożenia funkcji przez liczbÄ™. Niech U ‚" F bÄ™dzie li-
niową podprzestrzenią F, w której jest określony ten sam iloczyn skalarny
Ć
| co w F i niech L będzie liniowym operatorem różniczkowym określonym
na przestrzeni U o wartościach należących do F.
Ć
Załóżmy, że L jest operatorem symetrycznym i dodatnio określonym tzn.:
Ć Ć
1. Lu|v = u|Lv dla każdego u, v " U
Ć
2. Lu|u e" 0 dla wszystkich u " U z wyjÄ…tkiem u = 0.
Można więc wykazać, że wzór:
Ć
u|v L = Lu|v
Ć
3
dla każdego u, v " U określa nowy iloczyn skalarny w przestrzeni U.
Iloczyn skalarny u|v nazywać będziemy iloczynem energetycznym względem
Ć
operatora L, a przestrzeń U z tym iloczynem skalarnym  przestrzenią ener-
Ć
getycznÄ… operatora L.

Normę ||u||L = u|u L nazywać będziemy normą energetyczną. Niech prze-
Ć
Ć
strzenie U,F oraz operator różniczkowy spełniają w/w założenia i niech po-
nadto UN oznacza N - wymiarową podprzestrzeń U.
Rozpatrzmy równanie różniczkowe
Ć
Lu = f, (1)
gdzie f jest daną funkcją należącą do F przy założeniu, że równanie to ma w
Ć
podprzestrzeni U jednoznaczne rozwiązanie. Rozwiązanie równania Lu = f
można przybliżyć funkcjami z przestrzeni UN, przy czym jako przybliżone
Ć
rozwiązanie Lu = f przyjmuje się taki element uN " UN, że:
||u - uN||L = min ||u - k
||L (2)
Ć Ć
k
"UN
Można udowodnić, że istnieje dokładnie jeden taki element uN. Metodę
tę nazywamy metodą Rayleigha-Ritza. Jeżeli elementami przestrzeni UN są
wielomianowe funkcje sklejane, to takÄ… metodÄ™ nazywamy metodÄ… elementu
skończonego, a funkcje un - elementem skończonym aproksymującym rozwią-
zanie u z przestrzeni UN.
2.2 Algorytm - Etapy rozwiÄ…zania problemu za pomocÄ…
FEM
Rozpatrujemy ogólny przypadek, gdy elementami UN są funkcje dowolnego
typu, niekoniecznie funkcje sklejane. Niech funkcje Ći(i = 1, ..., N) stanowią
bazę przestrzeni UN można przedstawić w postaci:
k
= c1Ć1 + c2Ć2 + . . . + cNĆN, (3)
4
gdzie ci sÄ… liczbami rzeczywistymi.
Elementu uN poszukiwać będziemy w tej postaci. Mamy:
g(c1, . . . , cN) = ||u - uN||2 (4)
Ć
L
N

= ||u - ciĆi||2
Ć
L
i=1

N N


= u - ciĆi u - cjĆj

Ć
L
i=1 j=1
N N

= u|u L - 2 cj u|Ćj L + cicj Ći|Ćj L
Ć Ć Ć
j=1 i,j=1
Element z przestrzeni UN, którego współczynnikami w rozwinięciu (5) są
liczby c", . . . , c" takie, że
1 N
g(c", . . . , c" ) = min g(c1, . . . , cN) (5)
1 N
(c1,...,cN )
jest poszukiwanym elementem uN.
"g
W celu wyznaczenia minimum funkcji g(c1, ..., cN) obliczamy pochodne ,
"ck
k = 1, ..., N, i przyrównujemy je do zera
N

"g
= -2 u|Ćk L + 2 cj Ćk|Ćj L = 0. (6)
Ć Ć
"ck
j=1
Ć
Funkcja u jest rozwiązaniem równania Lu = f, zatem
Ć
u|Ćk L = Lu|Ćk = f|Ćk (7)
Ć
po podstawieniu (6) do równań (7) otrzymamy
N

Ćk,jcj =< f, Ćk > (8)
j=1
gdzie k = 1, . . . , N, a Ćk,j = f|Ćk L.
Ć
Równania (8) tworzą układ N równań z N niewiadomymi ck. Można udo-
wodnić, że macierz tego układu jest macierzą nieosobliwą. Układ (8) posiada
więc jednoznaczne rozwiązanie. Rozwiązanie to daje punkt w którym funkcja
g osiÄ…ga minimum.
Po ustaleniu więc bazy Ći, i = 1, . . . , N, w przestrzeni UN, należy rozwią-
zać układ (8). Można również pokazać, że macierz układu (8) jest macie-
rzą symetryczną i dodatnio określoną. Są to fakty ważne z numerycznego
5
punktu widzenia. Dla takich układów opracowanych jest wiele efektywnych
algorytmów ich rozwiązywania. Możemy zatem przedstawić kolejne etapy
rozwiÄ…zywania problemu za pomocÄ… FEM:
1. Analizowany obszar dzieli się myślowo na pewną skończoną liczbę geo-
metrycznie prostych elementów, tzw. elementów skończonych.
2. Zakłada się, że te połączone są ze sobą w skończonej liczbie punk-
tów znajdujących się na obwodach. Najczęściej są to punkty narożne.
Noszą one nazwę węzłów. Poszukiwane wartości wielkości fizycznych
stanowią podstawowy układ niewiadomych.
3. Obiera się pewne funkcje jednoznacznie określające rozkład analizowa-
nej wielkości fizycznej wewnątrz elementów skończonych, w zależności
od wartości tych wielkości fizycznych w węzłach. Funkcje te noszą na-
zwę funkcji węzłowych lub funkcji kształtu.
4. Równania różniczkowe opisujące badane zjawisko przekształca się, po-
przez zastosowanie tzw. funkcji wagowych, do równań metody elemen-
tów skończonych. Są to równania algebraiczne.
5. Na podstawie równań metody elementów skończonych przeprowadza się
asemblację układu równań, tzn. oblicza się wartości współczynników
stojących przy niewiadomych oraz odpowiadające im wartości prawych
stron. Jeżeli rozwiązywane zadanie jest niestacjonarne, to w obliczaniu
wartości prawych stron wykorzystuje się dodatkowo warunki począt-
kowe. Liczba równań w układzie jest równa liczbie węzłów przemno-
żonych przez liczbę stopni swobody węzłów, tzn. liczbę niewiadomych
występujących w pojedynczym węz
le.
6. Do tak utworzonego układu równań wprowadza się warunki brzegowe.
Wprowadzenie tych warunków następuje poprzez wykonanie odpowied-
nich modyfikacji macierzy współczynników układu równań oraz wek-
tora prawych stron.
7. Rozwiązuje się układ równań otrzymując wartości poszukiwanych wiel-
kości fizycznych w węzłach.
8. W zależności od typu rozwiązywanego problemu, lub potrzeb, oblicza
się dodatkowe wielkości.
9. Jeżeli zadanie jest niestacjonarne, to czynności opisane w pkt. 5, 6, 7 i 8
powtarza się aż do momentu spełnienia warunku zakończenia obliczeń.
Może to być np. określona wartość wielkości fizycznej w którymś z
węzłów, czas przebiegu zjawiska lub jakiś inny parametr.
6
2.3 Definicja elementu skończonego
Element skończony jest prostą figurą geometryczną (płaską lub przestrzenną),
dla której określone zostały wyróżnione punkty zwane węzłami, oraz pewne
funkcje interpolacyjne służące do opisu rozkładu analizowanej wielkości w
jego wnętrzu i na jego bokach. Funkcje te nazywa się funkcjami węzłowymi,
bÄ…dz
funkcjami kształtu. Węzły znajdują się w wierzchołkach elementu skoń-
czonego, ale mogą być również umieszczone na jego bokach i w jego wnętrzu.
Jeżeli węzły znajdują się tylko w wierzchołkach, to element skończony jest
nazywany elementem liniowym (ponieważ funkcje interpolacyjne są wtedy
liniowe). W pozostałych przypadkach mamy do czynienia z elementami wyż-
szych rzędów. Rząd elementu jest zawsze równy rzędowi funkcji interpola-
cyjnych (funkcji kształtu). Liczba funkcji kształtu w pojedynczym elemencie
skończonym jest równa liczbie jego węzłów. Funkcje kształtu są zawsze tak
zbudowane, aby w węzłach których dotyczą ich wartości wynosiły jeden, a
pozostałych węzłach przyjmowały wartość zero.
7
3 Siatka elementów skończonych
3.1 Numeracja elementów i węzłów
Wszystkie elementy skończone oraz węzły muszą być ponumerowane. Naj-
częściej dąży się do tego, aby numeracja węzłów zapewniała jak najmniejszą
szerokość pasma współczynników niezerowych macierzy układu równań. Po-
nieważ macierz ta najczęściej jest symetryczna, stąd w pamięci komputera
trzeba przechowywać tylko fragment jej niezerowego pasma, np. główną dia-
gonalę i część powyżej tej diagonali. Daje to dużą oszczędność miejsca w
pamięci komputera oraz możliwość szybszego rozwiązania układu równań.
Układ równań buduje się najczęściej w ten sposób, że odpowiednie warto-
ści współczynników macierzy i wektora prawych stron oblicza się dla kolej-
nych elementów skończonych i sumuje w odpowiednich miejscach macierzy
 i wektora b, które sÄ… zwiÄ…zane równaniem liniowym Âx = b. Ponieważ
własności materiałowe mogą być inne w różnych podobszarach, dlatego aby
skrócić czas obliczeń korzystnym jest, aby elementy skończone w tych pod-
obszarach posiadały kolejne numery. Wprowadzenie warunków brzegowych
następuje poprzez węzły elementów brzegowych.
3.2 Tablica koneksji
Numery węzłów tworzących poszczególne elementy skończone zapisuje się w
tzw. Tablicy koneksji. Węzły wszystkich elementów muszą być zapisane
w takim samym obiegu - lewoskrętnym bądz
prawoskrętnym. Najczęściej
stosuje się zapis w obiegu lewoskrętnym, ponieważ równania metody elemen-
tów skończonych są wyprowadzane w ten sposób, że wtedy otrzymuje się
dodatnią wartość jakobianu przekształcenia, tj. wielkości charakterystycz-
nej dla poszczególnych typów elementów skończonych. Przykład siatki dla
FEM przedstawiono na rys.1. Przyjętej numeracji węzłów odpowiada ma-
cierz układu równań pokazana na rys.2. Nietrudno zauważyć, że macierz ta
jest macierzÄ… rzadkÄ… o specjalnej strukturze znanej w literaturze jako tzw.
strzałka.
8
Rysunek 1: Siatka elementów skończonych
Rysunek 2: Macierz rzadka tzw. strzałka
9
4 Zastosowania metody elementu skończonego
4.1 Rozwiązywanie równań różniczkowych
Metodę elementu skończonego można stosować do rozwiązania zarówno rów-
nań różniczkowych zwyczajnych, jak i cząstkowych. Opis poniżej dotyczy
rozwiązywania równania różniczkowego zwyczajnego, równaniem różniczko-
wym cząstkowym zajęłam się w kolejnym podrozdziale na przykładzie rów-
nania dyfuzji. Rozpatrzmy następujące zagadnienie brzegowe.
Wyznaczyć funkcję k
spełniajacą:
Ć Ż
Lk
(x) = k
(x) + Ã(x)k
(x) = f(x), (9)
z warunkami brzegowymi w postaci: k
(0) = Ä… oraz k
(1) = ², gdzie Ã(x),
Å»
f(x) - ustalone funkcje ciÄ…gÅ‚e na przedziale < 0; 1 >, Ã(x) e" 0 i Ã(x) różna
od 0 na < 0; 1 >, a Ä…, ² - ustalone liczby rzeczywiste. Jak wiadomo, tak
postawione zagadnienie ma jednoznaczne rozwiÄ…zanie:
k
" C2(< 0; 1 >). Nasze zagadnienie sprowadzamy do zagadnienia z jedno-
rodnymi warunkami brzegowymi. Można to zrobić na przykład w następu-
jący sposób:
niech v(x) " C2(< 0; 1 >) będzie ustaloną funkcją spełniającą warunki
v(0) = Ä…, v(1) = ². Funkcja u(x) = k
(x) - v(x) spełnia więc równanie:
Ć
Lu(x) = u (x) + Ã(x)u(x) = f(x), (10)
wraz z warunkami brzegowymi w postaci: u(0) = 0, u(1) = 0, natomiast
Ż Ć
f(x) = f - Lv(x).
4.2 Modelowanie równoległe z wykorzystaniem Metody
Elementów Skończonych.
Jako przykład zastosowania metody elementu skończonego zajmę się niejed-
norodnym równaniem przewodnictwa cieplnego zależnego od czasu (równanie
dyfuzji), gdzie niejednorodność stanowi z
ródło ciepła:
"Q "T
lim = "x"("T ) = cÁ , (11)
"V 0"t0
"V "t "t
gdzie T oznacza temperaturę, t - czas, - współczynnik przewodnictwa ciepl-
nego, Á - gÄ™stość gazu, c - ciepÅ‚o wÅ‚aÅ›ciwe, Q -ciepÅ‚o wymienione przez objÄ™-
tość "V w ciągu czasu "t.
Równanie to możemy zapisać w postaci:
"T "2T
= D (12)
"t "2x
10
gdzie D nazywamy współczynnikiem dyfuzji. Wartość D dla strumienia cie-
pła zgodnego z osią X (a tak zakładamy dla uproszczenia obliczeń) wynosi
D = (1/3) < v >< l >.
Zastosowanie powyższych dwóch równań wraz z odpowiednimi warunkami
brzegowymi i początkowymi pozwala otrzymać temperaturowy opis dyfuzji.
Uwzględniając pojęcie entalpii (różniczkując wyrażenie na entalpie po tempe-
raturze) dochodzimy do uniwersalnego wzoru, opisującego równanie dyfuzji:
n+1
AnT = bn
Rozwiązanie poprzedza wyznaczenie macierzy współczynników An oraz wek-
tora prawych strony bn . Do implementacji równoległej tych obliczeń można
wykorzystać powszechnie stosowaną w obliczeniach równoległych koncepcję
dekompozycji geometrycznej , która sprowadza się tutaj do podziału dużego
i obliczeniowo czasochłonnego zagadnienia elementów skończonych na mniej-
sze i łatwiejsze do operowania pod-zagadnienia, które dają się efektywnie
rozwiązywać. Proces opisanego podziału zwany jest również dekompozycją
obszaru (ang. domain decomposition). Stosując powyższe podejście, siatka
elementów skończonych zostaje podzielona na skończoną liczbę podobszarów
(podsiatek) tak, aby obciążenie obliczeniowe przypadające na poszczególne
pod-siatki było w przybliżeniu takie samo, co umożliwia efektywną realiza-
cję obliczeń dla poszczególnych pod-siatek przez różne procesy, utożsamiane
często z różnymi procesorami. W rezultacie, opisane podejście umożliwia
implementację równoległą (rozproszoną) zagadnienia modelowania dla sia-
tek o bardzo dużej liczbie elementów poprzez wykorzystanie wielu proce-
sów (komputerów) połączonych wzajemnie za pomocą sieci komunikacyjnej.
W przypadku ogólnym, gdy mamy do czynienia z nieregularnymi siatkami
FEM, zadanie podziału siatki na podsiatki, które są zatem przetwarzane
w sposób równoległy na poszczególnych procesorach, nie jest łatwe, gdyż
trzeba kierować się nie tylko chęcią zrównoważenia obciążenie przypadające
na poszczególne procesory. Wymagana jest również minimalizacja operacji
komunikacyjnych, potrzebnych do wymiany danych pomiędzy procesorami.
Komunikacja ta jest wynikiem sprzężenia występującego pomiędzy węzłami
granicznymi pod-obszarów. Przykład opisanego podziału przedstawiony jest
na rys. 3, odpowiadającym siatce FEM zawierającej 13 węzłów i 16 ele-
mentów trójkątnych. Została ona podzielona na cztery pod-siatki, z których
każda jest odwzorowana na odpowiadający jej procesor. W każdej podsiatce
będziemy wyróżniać trzy typy węzłów:
1. Węzły wewnętrzne czyli takie, które są sprzężone tylko z węzłami na-
leżącymi do danej podsiatki.
11
2. Węzły graniczne to węzły, które są sprzężone z węzłami należącymi do
innych pod-siatek. Węzły wewnętrzne i graniczne nazywane są węzłami
lokalnymi.
3. Nazwę węzłów zewnętrznych nadamy tym węzłom lokalnym innych
podobszarów, które są sprzężone z węzłami lokalnymi rozpatrywanego
podobszaru.
Rysunek 3: Przykładowy podział siatki elementów skończonych na podsiatki
wraz z przydziałem węzłów do procesorów.
12
Literatura
[1] Z.Fortuna, B.Macukow, J.WÄ…sowski: Metody numeryczne Wyd. Na-
ukowo - Techniczne, Warszawa 1993.
[2] A.N.Matwiejew: Fizyka czÄ…steczkowa PWN, Warszawa 1989.
[3] Tao Pang:Metody obliczeniowe w fizyce: PWN, Warszawa 2001.
[4] Strona www: http://icis.pcz.czest.pl/ roman.
[5] Strona www: http://www.owlnet.rice.edu.
13