Transcript Document

Metoda elementów
skończonych
Ludwik Antal - Numeryczna analiza pól elektromagnetycznych –W6
FEM - MES
Metoda elementów skończonych (MES) [Finite element
method (FEM)] ma swoje źródło w analizie strukturalnej.
Chociaż jej podstawy stworzył Courant już w 1943r.,
metoda jest stosowana w rozwiązywaniu zagadnień
elektromagnetycznych od roku 1968.
Metoda różnic skończonych (finite difference method - FDM) i
metoda momentów the (method of moments - MOM) to
metody koncepcyjnie prostsze i łatwiejsze w programowaniu
niż metoda elementów skończonych. Jednak MES jest
metodą zdecydowanie silniejszą i bardziej uniwersalną w
rozwiązywaniu problemów o złożonej geometrii i
niejednorodnych środowiskach. Aplikacja MES stworzona dla
określonej dyscypliny łatwo może być zastosowana do innej.
2
4 kroki
Zasadnicze etapy analizy MES
• dyskretyzacja regionu rozwiązania skończoną ilością
subregionów lub elementów,
• wyprowadzenie równań dla typowego elementu,
• złożenie wszystkich elementów w regionie rozwiązania,
• rozwiązanie uzyskanego układu równań.
3
Rozwiązanie problemu
Rozwiązanie problemu za pomocą metody elementów
skończonych wymaga następujących czynności:
1. Analizowany obszar dzieli się na pewną skończoną liczbę
geometrycznie
prostych
elementów,
tzw. elementów
skończonych.
2. Zakłada się, że są one połączone ze sobą w skończonej
liczbie punktów znajdujących się na obwodach. Najczęściej są
to punkty narożne. Noszą one nazwę węzłów.
3. Obiera się pewne funkcje jednoznacznie określające rozkład
analizowanej wielkości fizycznej wewnątrz elementów
skończonych, zależne od wartości tych wielkości fizycznych w
węzłach. Funkcje te noszą nazwę funkcji węzłowych lub funkcji
kształtu.
4. Równania
różniczkowe
opisujące
badane
zjawisko
przekształca się, przy pomocy tzw. funkcji wagowych, do
równań metody elementów skończonych. Są to równania
4
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ątkowe. 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ęźle.
6. Do układu równań wprowadza się warunki brzegowe przez
wykonanie
odpowiednich
modyfikacji
macierzy
współczynników układu równań oraz wektora prawych stron.
7. Rozwiązuje się układ równań otrzymując wartości
poszukiwanych wielkości fizycznych w węzłach.
5
8. W zależności od typu rozwiązywanego problemu, lub
potrzeb, oblicza się dodatkowe wielkości (energię, siły,
impedancje itp.).
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
Definicja:
Definicja
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.
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.
7
Funkcje interpolacyjne
Funkcje interpolacyjne nazywa się funkcjami węzłowymi, bądź
funkcjami kształtu.
Rząd elementu jest zawsze równy rzędowi funkcji
interpolacyjnych (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 w pozostałych
węzłach przyjmowały wartość zero.
8
Typowe elementy jedno-, dwu- i trójwymiarowe.
Typowe elementy
9
Rozwiązanie
równania
Laplacea
Rozwiązanie równania Laplace’a
 2V  0
Dyskretyzacja
Aby znaleźć rozkład potencjału
V (x, y) w dwuwymiarowym
regionie rozwiązania, dzielimy
region na elementy skończone.
Elementy 6, 8, i 9 to 4-węzłowe
czworokąty, a pozostałe 3-węzłowe
trójkąty.
Dla ułatwienia obliczeń w praktyce
stosuje się elementy tego samego
typu. Dlatego elementy
czworokątne należy podzielić na
trójkąty.
10
Szukamy aproksymacji dla potencjału Ve wewnątrz elementu e i
następnie zależności rozkładu potencjału od różnych elementów,
jako że potencjał jest ciągły poprzez granice międzyelementowe.
Aproksymacja rozwiązania dla całego regionu:
Szukamy aproksymacji
V(x, y) 
N
V (x, y)
e 1
e
gdzie N jest liczbą elementów trójkątnych
Najpopularniejszą formą aproksymacji Ve wewnątrz elementu
jest aproksymacja wielomianowa:
dla trójkąta
dla czworokąta
Ve(x, y)  a  bx  cy
Ve(x, y)  a  bx  cy  dxy
11
Pole elektryczne wewnątrz
Potencjał Ve wewnątrz elementu jest niezerowy, ale
zerowy na zewnątrz.
Zauważmy, że założenie liniowej kombinacji potencjału
wewnątrz elementu trójkątnego jest równoważne
przyjęciu, że pole elektryczne wewnątrz elementu jest
jednorodne tzn. ,
 Ve
Ve 
Ee  Ve  
ax 
a y   ba x  ca y 
y
 x

ax i ay – wektory jednostkowe
Ex  b
Ey  c
12
Wyprowadzenie równań
Wyprowadzenie równań dla typowego elementu
Potencjały Ve1, Ve2, i Ve3
opisane są równaniem:
Ve(x, y)  a  bx  cy
13
Ve
Współczynniki a, b i c są określone równaniem poprzednim,
jako
Podstawiając to do równanie na Ve otrzymamy
lub
3
Ve   i x, y Vei
i 1
14
Alfa i A
i – liniowe funkcje interpolacyjne (funkcje kształtu)
15
1 lub 0
Funkcje kształtu maja następujące własności:
1
x2 y3  x3 y2    y2  y3 x  x3  x2  y 
1 
2A
Dla x = x1 i y = y1
1 
1
x2 y3  x3 y2   x1 y2  x1 y3   x3 y1  x2 y1   2 A  1
2A
2A
Dla x = x2 i y = y2
2 
1
x2 y3  x3 y2   x2 y2  x2 y3   x3 y2  x2 y2   0
2A
16
Suma = 1
Bo sumowanie daje
↓
2A
↓
0
17
„A” jest powierzchnią elementu e (pole powierzchni trójkąta)
„A” jest dodatnie jeżeli
numeracja węzłów jest
przeciwna do ruchu
wskazówek zegara
18
Funkcje kształtu
Ilustracja funkcji kształtu
19
Funkcjonał
Funkcjonał korespondujący z równaniem Laplace’a,
∆V = 0, jest dany przez energię pola
Uwaga.
Podstawowym zagadnieniem rachunku wariacyjnego jest wyznaczenie takich
niewiadomych funkcji ui(x, y) (ekstremali), żeby całka I pewnej kombinacji tych
funkcji i ich pochodnych przybierała wartość ekstremalną. Całka I nazywana jest
funkcjonałem.
Największe znaczenie ma funkcjonał typu:
I  F u, u , u , x, y  dxdy
Typowe zagadnienie wariacyjne
 
x
y

I   F u, ux , u y , x, y  dxdy  min

Poszukuje się kombinacji liniowych ciągu funkcji, które od razu
ekstremalizują funkcjonał.
20
Fizycznie
Fizycznie funkcjonał We jest energią na jednostkę długości
związaną z elementem e.
3
Z równania
Ve   i x, y Vei
i 1
wynika
3
Ve  Vei i
i 1
Podstawiamy to do funkcjonału
21
Macierze
Jeżeli oznaczymy
to funkcjonał możemy zapisać w formie macierzowej jako
Gdzie t oznacza macierz transponowaną
⇗
macierz współczynników elementu
22
Cij może być uważany za sprzężenie
Cij
między węzłami i, j.
Na przykład
⇦A
Podobnie
A także
23
Asemblacja
Asemblacja wszystkich elementów
Energia związana z kompletem elementów
Gdzie
n jest liczbą węzłów, N jest liczbą elementów,
a [C] to tak zwana globalna macierz
współczynników która jest zestawem macierzy
współczynników indywidualnych elementów.
24
Niejednorodne
Równanie energii uzyskano dla założonej jednorodności regionu
(=const). Jeżeli region jest niejednorodny to dzieli się go na
elementy tak (rysunek) by każdy element skończony był
jednorodny.
25
epsilon
W takiej sytuacji, równanie
jest nadal obowiązujące, ale równanie
nie może być stosowane ponieważ =0r zmienia
się od elementu do elementu. By pokonać tę
trudność trzeba zastąpić  przez 0, i pomnożyć
funkcję podcałkową przez r .
26
Proces w którym macierze współczynników poszczególnych
elementów są składane w celu uzyskania globalnej macierzy
współczynników ilustruje przykład.
Globalna
Załóżmy, że siatka elementów skończonych składa się z trzech
elementów.
Numeracja: zewnętrzne 1 ,2, 3, 4 i 5
wewnętrzne 1, 2 i 3
– globalna
– lokalna ↺
27
Globalna 1
Globalna macierz współczynników
5x5 ponieważ siatka ma 5
węzłów (n = 5).
Cij jest sprzężeniem między węzłami i - j . Wyznaczamy Cij
wykorzystując fakt, że rozkład potencjału musi być ciągły na
granicy międzyelementowej. Wkład w i, j wyraz macierzy [C]
pochodzi od wszystkich elementów zawierających węzły i oraz j.
Na przykład: elementy 1 i 2 mają wspólny węzeł 1, stąd
28
Globalna 2
Węzeł 2 należy tylko do elementu 1, stąd
Węzeł 4 należy do elementów 1, 2, i 3; więc
Węzły 1 i 4 należą jednocześnie do elementów 1 i 2; stąd
Ponieważ nie ma sprzężenia (bezpośredniego połączenia) między węzłami 2 i 3,
Kontynuując postępowanie można wyznaczyć wszystkie wyrazy
globalnej macierzy współczynników przez dalszą analizę siatki
29
Globalna 3
∥
Odnotujmy, że mamy 27 wyrazów (9 dla każdego z 3 elementów)
w globalnej macierzy współczynników [C].
30
Własności macierzy [C]
Własności macierzy [C]:
1. Jest symetryczna (Cij = Cji ) tak jak macierz współczynników
elementu.
2. Ponieważ Cij = 0 jeśli nie ma sprzężenia między węzłami i i j
, oczekuje się, że dla dużej liczby elementów [C] stanie się
rzadka. Macierz [C] jest także pasmowa jeśli węzły
numerowane są właściwie. Można dowieść, że
3. Jest osobliwa. Chociaż nie jest to oczywiste można to
udowodnić przy użyciu macierzy współczynników elementu.
31
Rozwiązanie równań
Rozwiązanie równań
Z metod wariacyjnych wynika, że równanie Laplace’a jest
spełnione kiedy całkowita energia w regionie rozwiązania
jest minimalna. Tak więc żądamy by pochodne cząstkowe
W po każdej węzłowej wartości potencjału były zero tzn.
W rozpatrywanym przykładzie do równania W
podstawiamy [C] i różniczkujemy W po V1.
32
Otrzymujemy
czyli
Ogólnie
prowadzi do
gdzie: n - liczba węzłów w siatce.
33
Metody
Z ostatniego równania dla wszystkich węzłów k = 1, 2, ...n
otrzymujemy układ równań z którego znajdujemy rozwiązanie
Rozwiązanie może być uzyskane metodą iteracyjną lub metodą
pasmową.
34
Metoda iteracyjna
Jeżeli V1 jest węzłem swobodnym (o szukanym potencjale)
Ogólniej w węźle k siatki o n węzłach
gdzie k jest węzłem swobodnym
Ponieważ Cki = 0 jeśli węzeł k nie jest bezpośrednio połączony z
węzłem i, tylko węzły połączone bezpośrednio z k mają wkład w Vk .
Równanie może być użyte iteratywnie do wszystkich węzłów
swobodnych.
35
Proces iteracji zaczyna się przypisaniem znanych wartości
potencjałów do węzłów ustalonych i wartości zerowych lub średnich
potencjałów węzłom swobodnym ( o nieznanej wartości).
gdzie Vmin i Vmax są minimalną i maksymalną wartością V węzłów ustalonych.
Z takimi początkowymi potencjałami liczone są potencjały węzłów
swobodnych. W końcu pierwszej iteracji znane są nowe wartości, które
w następnej iteracji zastąpią stare.
Procedura jest powtarzana dopóki różnica między iteracjami nie osiągnie
żądanego poziomu.
36
Metoda pasmowa
Gdzie indeksy f and p, odpowiednio odnoszą się do węzłów
o szukanym (free) i zadanym (fixed) potencjale.
Ponieważ Vp jest stały, różniczkujemy po Vf i otrzymujemy
lub
37
Warunek brzegowy
Neuman’a
Warunek brzegowy Neumann’a
Warunek Neumann’a (∂V/∂n = 0) może
być warunkiem brzegowym lub określać
symetrię problemu. Załóżmy, że region
jest symetryczny wzdłuż osi y.
Wprowadzamy warunek ∂V/∂x = 0
wzdłuż osi y przyjmując:
38
Bezpośrednie i iteracyjne metody rozwiązywania układów równań
liniowych a macierze rzadkie
Bezpośrednie i iteracyjne
metody
Rozwiązywanie równań różniczkowych cząstkowych metodami MRS i MES
wymaga rozwiązywania układów równań liniowych z macierzami rzadkimi.
Metody rozwiązywania tych układów dzielą się na:
•metody bezpośrednie - dają rozwiązanie po skończonej liczbie kroków;
wykorzystują dekompozycję Gaussa, Choleskiego itd. Podstawowa
niedogodność stosowania metod bezpośrednich dla macierzy rzadkich to
pojawianie się nowych niezerowych elementów w macierzy w trakcie
obliczeń (ang. fill-in);
•metody iteracyjne - polegają na iteracyjnym ulepszaniu przybliżonego
rozwiązania do momentu osiągnięcia zadawalającej dokładności.
W każdym kroku wykorzystywana jest stosunkowo prosta procedura bazującej
na iloczynie macierzy przez wektor; procedura ta nie zmienia macierzy A
układu. Pozwala to zmniejszyć wymagania w stosunku do pamięci operacyjnej
w porównaniu z metodami bezpośrednimi oraz lepiej wykorzystać specyficzną
strukturę macierzy A. Dzięki zastosowaniu uwarunkowania wstępnego (ang.
preconditioning) udaje się też zmniejszyć wymaganą liczbę iteracji. Duże układy
równań, które zwykle charakteryzują się rzadkimi macierzami współczynników,
39
rozwiązuje się przeważnie przy użyciu metod iteracyjnych .