Metody numeryczne

Materiał przygotował: Grzegorz Suder

1. Metoda różnic skończonych

Metoda różnic skończonych jest najstarszą z metod numerycznych, była już ona rozpatrywana
prze Gaussa. Jednakże dopiero wiek dwudziesty pozwolił na pełne wykorzystanie tkwiących
w niej zalet. Pomimo bowiem prostoty niezbędnych obliczeń jakie należy wykonać, ich liczba
nawet w najprostszym wypadku, jest często zbyt duża by można wykonać ją ręcznie.
Wykorzystując do tego celu komputer, można w stosunkowo krótkim czasie dokonać
obliczeń i przedstawić wynik graficznie.
Metoda różnic skończonych bazuje na przekształceniu pochodnych cząstkowych, w równania
różnic cząstkowych, opisujące proces fizyczny przy pomocy przybliżeń algebraicznych.

Rys. 1.1

Rozważmy dwuwymiarowe równanie Laplace’a:

0

)

,

(

)

,

(

2

2

2

2

=

∂

∂

+

∂

∂

y

y

x

f

x

y

x

f

[1.1]

dla

dowolnego obszaru z określonymi warunkami granicznymi (rys. 1.1). Cały obszar zostaje

podzielony siatką kwadratową o rozmiarach

y

x

∆

∆

,

(równie dobrze mogłaby to być siatka

prostokątna, ale kwadratowa znacznie ułatwia obliczenia). Wybieramy dowolny punkt na
siatce o współrzędnych (i, j) wraz z jego otoczeniem i określamy potencjał tego punktu.
Korzystając z zapisu funkcji w szereg Taylora, można otrzymać:

...

!

3

)

(

!

2

)

(

,

3

3

3

,

2

2

2

,

,

,

1

+

∂

∂

∆

+

∂

∂

∆

+

∂

∂

∆

+

=

+

j

i

j

i

j

i

j

i

j

i

x

x

x

x

x

x

φ

φ

φ

φ

φ

[1.2]

oraz

...

!

3

)

(

!

2

)

(

,

3

3

3

,

2

2

2

,

,

,

1

+

∂

∂

∆

−

∂

∂

∆

+

∂

∂

∆

−

=

−

j

i

j

i

j

i

j

i

j

i

x

x

x

x

x

x

φ

φ

φ

φ

φ

[1.3]

Po dodaniu stronami, otrzymujemy:

)

)

((

)

(

2

4

,

2

2

2

,

,

1

,

1

x

O

x

x

j

i

j

i

j

i

j

i

∆

+

∂

∂

∆

+

=

+

−

+

φ

φ

φ

φ

[1.4]

Zaniedbując (ze względu na czwartą i wyższe potęgi), resztę

)

)

((

4

x

O

∆

, otrzymujemy:

2

,

1

,

,

1

,

2

2

)

(

2

x

x

j

i

j

i

j

i

j

i

∆

+

−

=

∂

∂

−

+

φ

φ

φ

φ

[1.5]

Postępując analogicznie dla zmiennej y, otrzymujemy:

2

,

1

,

,

1

,

2

2

)

(

2

y

y

j

i

j

i

j

i

j

i

∆

+

−

=

∂

∂

−

+

φ

φ

φ

φ

[1.6]

Podstawiając wyniki do równania Laplace’a, otrzymujemy:

0

)

2

(

)

(

1

)

2

(

)

(

1

,

1

,

,

1

2

,

1

,

,

1

2

=

+

−

∆

+

+

−

∆

−

+

−

+

j

i

j

i

j

i

j

i

j

i

j

i

y

x

φ

φ

φ

φ

φ

φ

[1.7]

W myśl wcześniejszego założenia, że siatka jest kwadratowa, podstawiamy

h

y

x

=

∆

=

∆

, i

dzielimy przez h. Otrzymujemy w ten sposób równanie:

0

,

1

,

1

,

,

1

,

1

=

−

+

+

+

+

−

+

−

j

i

j

i

j

i

j

i

j

i

φ

φ

φ

φ

φ

[1.8]

bądź

4

1

,

1

,

,

1

,

1

,

+

−

+

−

+

+

+

=

j

i

j

i

j

i

j

i

j

i

φ

φ

φ

φ

φ

[1.9]

Podstawiając

V

=

φ

, otrzymujemy wzory na potencjał w punkcie (i, j).

Głównym problemem metody różnic skończonych jest fakt, że w przypadku nieregularnych
powierzchni, mogą wystąpić problemy z dopasowaniem siatki. Dodatkowo, rodzaj materiału
warunki symetrii oraz charakterystyki nieliniowe wymagają specjalnego potraktowania.
Zwiększenie dokładności obliczeń (przez uwzględnienie potęg wyższych od 4 w szeregu
Taylora) zwiększa zarazem czas jak i ilość obliczeń potrzebnych do uzyskania rozwiązania.

Rozwiązując problemy, z wykorzystaniem metody różnic skończonych, można wykorzystać
następujące sposoby:

Metoda macierzowa

-

Pierwszym krokiem, jest ustalenie dziedziny, w której znajduje się
rozwiązanie.

-

Następnie zostają nałożone warunki graniczne (ponieważ bez ich
nałożenia nie można przekształcić pochodnych w równania
różnicowe). To założenie jest istotne zarówno dla metody
analitycznej, jak i numerycznej.

-

Cała dziedzina zostaje podzielona siatką; dla wygody wybieramy
siatkę kwadratową (o wielkości oczka

y

x

∆

=

∆

).

-

Dla nieznanych punktów wewnętrznych, stosowane jest przybliżenie.
Do tego celu, wykorzystywanych jest pięć punktów, łącznie z
punktem, dla którego szukamy potencjału.

-

Po zapisaniu równań w macierzy, można wyliczyć wartości
szukanych potencjałów w danym punkcie.

-

Teraz można dokonać prezentacji wyniku, np. w formie graficznej.

Przykład

Rys. 1.2a

Rys. 1.2b

Rys. 1.2c

Niech będzie dany cylindryczny kondensator, o przekroju jak na rys. 1.2a. Ze względu na
jego symetrię, wystarczy rozważyć jedną ćwiartkę układu (rys. 1.2b). Po zastosowaniu
wyprowadzonego wzoru (

0

4

,

1

,

1

,

,

1

,

1

=

−

+

+

+

+

−

+

−

j

i

j

i

j

i

j

i

j

i

V

V

V

V

V

) dla punktów 1, 2, 3, 4, 5

(rys. 1.2c) o nieznanych potencjałach

5

4

3

2

1

i

,

,

,

V

V

V

V

V

, otrzymujemy następujące równania:

















=

−

+

+

+

=

−

+

+

+

=

−

+

+

+

=

−

+

+

+

=

−

+

+

+

0

4

100

0

0

4

100

100

0

4

0

100

0

4

0

100

0

4

0

100

5

4

4

4

5

3

3

4

2

2

3

1

1

2

2

V

V

V

V

V

V

V

V

V

V

V

V

V

V

V

,

[1.10]

co z kolei można zapisać w macierzy:

































=

































































−

−

−

−

−

−

−

−

100

200

100

100

100

4

2

0

0

0

1

4

1

0

0

0

1

4

1

0

0

0

1

4

1

0

0

0

2

4

5

4

3

2

1

V

V

V

V

V

.

[1.11]

Aby istniało rozwiązanie, musi istnieć macierz odwrotna, co oznacza, że liczba równań jest
równa liczbie niewiadomych. Rozwiązaniem dla tego przykładu, jest macierz odpowiadająca
badanej powierzchni (znajdującej się w pierwszej ćwiartce układu współrzędnych):





















100

0

0

0

100

100

100

100

100

100

65.477

80.953

58.333

52.381

51.191

[1.12]

.

Natomiast dla całej powierzchni, macierz będzie wyglądała następująco:

































100

100

100

100

100

100

100

100

100

100

100

100

0

0

0

0

0

100

100

100

100

100

100

100

100

100

100

100

100

80.953

58.333

52.381

51.191

52.381

58.333

80.953

65.477

65.477

80.953

58.333

52.381

51.191

52.381

58.333

80.953

[1.13]

Metoda iteracyjna

-

Pierwsze trzy punkty (dziedzina, warunek z granicą oraz podział na
siatkę) są takie same jak w przypadku poprzedniej metody. Różnica
polega na sposobie dokonania przybliżenia. O ile bowiem, metoda
macierzowa zakładała, że potencjały wewnętrzne są nieznane, o tyle
metoda jawna zakłada, że potencjały wewnętrzne mają wartość
równą zero.

-

Liczymy potencjał w pierwszym punkcie uwzględniając podane
założenia (wartości nieznanych potencjałów równe zeru); w każdym
następnym punkcie z uwzględnieniem wyniku otrzymanego dla
poprzedniego punktu.

-

Wszystkie obliczenia zostają powtórzone, przy czym potencjały
sąsiednie mają wartość różną od zera. Obliczenia przeprowadzamy n
razy, tj. dotąd, dopóki wynik nie zmieści się w granicy założonego
błędu.

-

Przybliżenie można wyrazić, korzystając ze wzoru:

∑

=

−

≤

−

N

i

n

i

n

i

e

V

V

N

1

1

1

, gdzie

i

V - potencjał w danym punkcie, n –

ilość iteracji, e – założony błąd.

Przykład

Rys. 1.3

Niech będzie dany obszar przedstawiony jak na rys. 1.3, o potencjałach V na brzegach,
równych 0 V, 20 V i 30 V. Na początku, w punktach 1, 2, 3, 4, 5, 6, 7 i 8 przyjmuje się
wartości potencjału równe zero;

0

8

7

6

5

4

3

2

1

=

=

=

=

=

=

=

=

V

V

V

V

V

V

V

V

. Po zastosowaniu

wzoru

4

1

,

1

,

,

1

,

1

,

+

−

+

−

+

+

+

=

j

i

j

i

j

i

j

i

j

i

V

V

V

V

V

dla kolejnych szukanych potencjałów, otrzymujemy

następujące równania:





























=

+

+

+

=

=

+

+

+

=

=

+

+

+

=

=

+

+

+

=



875

.

1

)

0

0

25

.

6

25

.

1

(

4

1

25

.

6

)

0

0

20

5

(

4

1

25

.

1

)

0

0

0

5

(

4

1

5

)

0

0

20

0

(

4

1

4

3

2

1

V

V

V

V

.

[1.14]

Po pierwszej iteracji, potencjały

8

1

,...,V

V

otrzymują nowe wartości. Uwzględniając otrzymane

wyniki w pierwszej iteracji, obliczenia zostają powtórzone w drugiej iteracji:



















=

+

+

+

=

=

+

+

+

=



187

.

2

)

875

.

1

0

0

875

.

6

(

4

1

875

.

6

)

25

.

6

25

.

1

20

0

(

4

1

2

1

V

V

.

[1.15]

Cykl ten, jest powtarzany, aż do osiągnięcia żądanej dokładności, zgodnie ze wzorem:

∑

=

−

≤

−

N

i

n

i

n

i

e

V

V

N

1

1

1

. Dla podanego przykładu, po przeprowadzeniu dziesięciu iteracji,

wartości potencjałów przedstawiają się następująco:

26

.

11

,

06

.

15

,

97

.

18

,

05

.

21

786

.

9

,

22

.

15

,

956

.

4

,

04

.

10

8

7

6

5

4

3

2

1

=

=

=

=

=

=

=

=

V

V

V

V

V

V

V

V

.

[1.16]

Na rys. 1.4 zostały przedstawione wartości potencjałów otrzymane dla poszczególnych
iteracji:

Rys. 1.4

Porównanie metod:

Metoda macierzowa

Metoda iteracyjna

•

wynik dla wewnętrznych
punktów jest dokładny

•

wolniejsze działanie

•

można otrzymać tylko wynik
przybliżony, zależny od ilości
iteracji

•

szybsze działanie i mniejsze
zużycie pamięci (zwłaszcza dla
bardzo dużych siatek)

Równanie Laplace’a zależne od czasu

Metodę różnic skończonych, można stosować również i dla innych równań spotykanych w
elektromagnetyzmie. Niech będzie dane równanie dyfuzji:

t

x

∂

∂

=

∂

∂

φ

α

φ

2

2

[1.17]

oraz równanie falowe:

2

2

2

2

t

x

∂

∂

=

∂

∂

φ

β

φ

.

[1.18]

Pochodną po czasie, można rozpisać w szereg Taylora, otrzymując:

t

t

k

i

k

i

∆

−

=

∂

∂

+

,

1

,

φ

φ

φ

,

2

1

,

,

1

,

2

2

)

(

2

t

t

k

i

k

i

k

i

∆

+

−

=

∂

∂

−

+

φ

φ

φ

φ

.

[1.19]

gdzie

t

∆

jest interwałem czasowym pomiędzy kolejnymi wartościami

i

φ

pojawiającymi się w

punkcie i, a przyrost k oznacza zmienną czasową. Postępując podobnie dla pochodnej po x,
otrzymujemy:

t

x

k

i

k

i

k

i

k

i

k

i

∆

−

=

∆

+

−

+

−

+

,

1

,

2

,

1

,

,

1

)

(

2

φ

φ

α

φ

φ

φ

[1.20]

dla równania dyfuzji, oraz

2

1

,

,

1

,

2

,

1

,

,

1

)

(

2

)

(

2

t

x

k

i

k

i

k

i

k

i

k

i

k

i

∆

+

−

=

∆

+

−

−

+

−

+

φ

φ

φ

β

φ

φ

φ

[1.21]

dla równania falowego, gdzie i jest zmienną przestrzeni, a

x

∆

odległością pomiędzy

sąsiednimi punktami. Zakładamy, że warunki graniczne są następujące:

)

,

( t

x

φ

jest dane dla

każdego t w x = 0 i x = L, oraz początkowe warunki dla

)

0

,

(x

φ

są wyszczególnione dla

każdego x. Do obliczenia, może być wykorzystana oczywiście metoda macierzowa jak i
iteracyjna. Niezależnie od wybranej metody, musi być ona kompatybilna z równaniem
różnicowym, oraz stabilna. Kompatybilność zapewnia, że rozwiązanie jest podobne do
rozwiązania oryginalnego równania. Stabilność można utracić, jeżeli suma błędów podczas
obliczeń będzie zbyt duża. Przykładowo, dla równania [1.20] stabilność jest zachowana w
przypadku spełnienia równania:

2

1

)

(

2

≤

∆

∆

x

t

α

[1.22]

a dla równania [1.21]:

1

1

2

≤













∆

∆

x

t

β

[1.23]

Rozwiązanie trójwymiarowego równania Laplace’a

Rozwiązanie problemu, jest trudniejszym zadaniem niż w przypadku dwuwymiarowym,
jednakże, tok postępowania jest podobny jak dla przypadku opisanego powyżej. Zasadniczą
różnicą jest fakt, że siatka jest tym razem trójwymiarowa. Uwzględniając fakt istnienia trzech
zmiennych, trójwymiarowe równanie Laplace’a można przedstawić jako:

0

6

1

,

,

1

,

,

,

1

,

,

1

,

,

,

1

,

,

1

,

,

=

−

−

−

−

−

−

+

−

+

−

+

−

k

j

i

k

j

i

k

j

i

k

j

i

k

j

i

k

j

i

k

j

i

V

V

V

V

V

V

V

[1.24]

dla metody macierzowej, oraz

6

1

,

,

1

,

,

,

1

,

,

1

,

,

,

1

,

,

1

,

,

+

−

+

−

+

−

+

+

+

+

+

=

k

j

i

k

j

i

k

j

i

k

j

i

k

j

i

k

j

i

k

j

i

V

V

V

V

V

V

V

[1.25]

dla metody iteracyjnej.

Przykład
Niech będzie dany obszar o potencjałach jak na rys. 1.5:

Rys. 1.5

Na obszar zostaje nałożona siatka o wymiarach 10x10x10. Potencjał

k

j

i

V

,

,

oblicza się

analogicznie jak dla przypadku dwuwymiarowego, biorąc pod uwagę sześć sąsiednich
punktów, jak na rys. 1.6:

Rys. 1.6

Jednym z możliwych rozwiązań graficznych dla równania Laplace’a w trzech wymiarach,
może być przedstawienie przekroju w punkcie x (bądź y).Inny sposób pokazany jest na (rys.
1.7):

Rys. 1.7

2. Metoda elementów skończonych

Metoda elementów skończonych, jest stosunkowo młoda; jej początki datuje się na 1950 r. W
ostatnich latach, zyskuje ona coraz większą popularność zarówno w obliczaniu pól
elektromagnetycznych, jak i innych dziedzinach wymagających rozwiązania problemów
inżynieryjnych. W przeciwieństwie do innych metod, metoda elementów skończonych
wymaga obliczeń z wykorzystaniem programów komputerowych, głównie ze względu na
dosyć skomplikowany algorytm.
Zasadniczo, działanie opiera się na podziale rozważanego obszaru, na podobszary, tzw.
elementy skończone. Każdy taki element, zdefiniowany jest poprzez ilość krawędzi; każde
przecięcie dwóch krawędzi, określa punkt. Podstawowym warunkiem rozwiązalności, jest
założenie, że każdy taki element ma skończony rozmiar; dodatkowo elementy muszą być ze
sobą kompatybilne.

Rys. 2.1a

Rys. 2.1b

przykładowe elementy skończone: trójkąt (2.1a) i prostokąt

(2.1b)

Metoda elementów skończonych, polega właściwie na odpowiednim dopasowaniu siatki do
powierzchni, a następnie znalezieniu minimum funkcji energii dla każdego elementu
skończonego. Dzięki zapisowi macierzowemu, można (poprzez dokonanie odpowiednich
przekształceń) znaleźć poszukiwane wartości potencjałów w odpowiednich punktach
(węzłach) siatki.
Zasadniczo, można wyróżnić następujące etapy postępowania:

1. dyskretyzacja, czyli podział płaszczyzny na elementy skończone (o

dowolnej ilości punktów, a więc i kształcie; dla najprostszego przypadku
jest to trójkąt). Im większe zmiany funkcji V w danym miejscu, tym
wielkość elementów jest mniejsza. Ważne jest, by nie było przerw w siatce
elementów, oraz by elementy nie nachodziły na siebie,

2. dobór funkcji aproksymującej V, by zachować ciągłość między

elementami,

3. utworzenie globalnego układu równań zapisanego w macierzy, z

uwzględnieniem warunków brzegowych,

4. rozwiązanie globalnego układu równań

Rozważmy dwuwymiarowe równanie Laplace’a:

0

2

2

2

2

2

=

∂

∂

+

∂

∂

=

∇

y

V

x

V

V

,

[2.1]

gdzie natężenie pola jest określone jako:









∂

∂

+

∂

∂

−

=

−

=

y

y

V

x

x

V

V

E

ˆ

grad



.

[2.2]

Teraz, należy zastosować zasadę ekstremum energii. Zasada ekstremum energii (zasada
minimum działania) mówi, że rzeczywisty rozkład pola elektrostatycznego zapewnia

minimum zgromadzonej energii. W rozważaniach dwuwymiarowych, energię można
przedstawić jako:

∫

∫

∫

∫

Ω

Ω

Ω

Ω

Ω

∇

=

Ω

























∂

∂

+













∂

∂

=

Ω

=

Ω

⋅

=

d

V

d

y

V

x

V

d

E

d

W

2

2

2

2

2

1

2

1

2

1

2

1

ε

ε

ε

D

E

.

[2.3]

By zachować ciągłość pomiędzy poszczególnymi elementami, należy odpowiednio dobrać
funkcję aproksymującą. Dla przypadku ogólnego, jest to:

...

2

2

+

+

+

+

+

+

=

fy

ex

dxy

cy

bx

a

V

.

[2.4]

Pod uwagę branych jest tyle wyrazów, ile jest punktów w elemencie skończonym. Dla
prostokąta, jest to:

dxy

cy

bx

a

V

+

+

+

=

, natomiast dla trójkąta:

cy

bx

a

V

+

+

=

. Zakładając,

że potencjał V w elemencie skończonym będzie przybliżony równaniem

cy

bx

a

V

+

+

=

,

można dla trójkąta (rys. 2.2)

Rys. 2.2

zapisać równania dla trzech punktów:











+

+

=

+

+

=

+

+

=

3

3

3

2

2

2

1

1

1

cy

bx

a

V

cy

bx

a

V

cy

bx

a

V

,

[2.5]

bądź w postaci macierzy:









































=





















c

b

a

y

x

y

x

y

x

V

V

V

3

3

2

2

1

1

3

2

1

1

1

1

[2.6]

gdzie (x

1

, y

1

), (x

2

, y

2

) i (x

3

, y

3

) są współrzędnymi wierzchołków trójkąta. Jeżeli teraz

przekształcimy macierz do postaci:









































=





















−

3

2

1

1

3

3

2

2

1

1

1

1

1

V

V

V

y

x

y

x

y

x

c

b

a

,

[2.7]

i podstawimy do równania

cy

bx

a

V

+

+

=

zapisanego w macierzy:

[

]





















=

c

b

a

y

x

V

1

,

[2.8]

to otrzymamy w ten sposób wartość potencjału w obszarze elementu skończonego:

[

]









































=

−

3

2

1

1

3

3

2

2

1

1

1

1

1

1

V

V

V

y

x

y

x

y

x

y

x

V

,

[2.9]

co z kolei można zapisać jako:

∑

=

=

3

1

)

,

(

i

i

i

y

x

N

V

V

, gdzie:

{

}

{

}

{

}

y

x

x

x

y

y

y

x

y

x

N

y

x

x

x

y

y

y

x

y

x

N

y

x

x

x

y

y

y

x

y

x

N

)

(

)

(

)

(

2

1

)

(

)

(

)

(

2

1

)

(

)

(

)

(

2

1

1

2

2

1

1

2

2

1

3

3

1

1

3

3

1

1

3

2

2

3

3

2

2

3

3

2

1

−

+

−

+

−

∆

=

−

+

−

+

−

∆

=

−

+

−

+

−

∆

=

.

[2.10]

natomiast

)

(

)

(

)

(

1

1

1

det

2

3

1

1

3

1

2

2

1

2

3

3

2

3

3

2

2

1

1

y

x

y

x

y

x

y

x

y

x

y

x

y

x

y

x

y

x

−

+

−

+

−

=





















=

∆

.

Funkcje N są nazywane funkcjami kształtu, i nie zależą od wartości potencjału V(x, y). Są one
równe 1 we własnym węźle, natomiast równe 0 w pozostałej części obszaru, np.

1

1

=

N

w

węźle 1 oraz

0

1

=

N

w pozostałych węzłach. Suma funkcji N jest równa 1 dla każdego

punktu wewnątrz elementu i na jego brzegach:

1

3

,

2

,

1

=

∑

=

i

i

N

.

W ten sposób, wyznaczenie funkcji V(x, y) w całym obszarze, zostało sprowadzone do
znalezienia wartości V w węzłach siatki.
Jeżeli przyjmiemy teraz, że dla każdego elementu energia wynosi:

∫

∇

=

e

e

dS

V

W

2

)

(

2

1

ε

,

[2.11]

gdzie

∑

=

∇

=

∇

3

1

)

,

(

i

i

i

y

x

N

V

V

, to energia każdego elementu można wyrazić wzorem:

(

)

j

i

j

e

j

i

i

e

V

dS

N

N

V

W

∑∑ ∫

=

=

∇

⋅

∇

=

3

1

3

1

)

(

2

1

ε

.

[2.12]

Przyjmując oznaczenia:

dS

N

N

N

e

j

i

e

j

i

∫

∇

⋅

∇

=

)

(

,

- macierz kwadratowa 3x3, oraz

T

]

[V -

transponowana macierz wartości potencjałów, otrzymujemy równanie zapisane w postaci
macierzowej:

]

[

]

[

]

[

2

1

)

(

T

)

(

V

N

V

W

e

e

ε

=

.

[2.13]

Teraz należy znaleźć gradienty funkcji kształtu:

{

}

{

}

{

}

y

x

x

x

y

y

N

y

x

x

x

y

y

N

y

x

x

x

y

y

y

y

N

x

x

N

N

ˆ

)

(

ˆ

)

(

2

1

ˆ

)

(

ˆ

)

(

2

1

ˆ

)

(

ˆ

)

(

2

1

ˆ

ˆ

1

2

2

1

3

1

1

3

2

2

3

3

2

1

1

1

−

+

−

∆

=

∇

−

+

−

∆

=

∇

−

+

−

∆

=

∂

∂

+

∂

∂

=

∇

.

[2.14]

Teraz należy znaleźć iloczyn skalarny gradientów funkcji, wykorzystując fakt, że dla
współrzędnych kartezjańskich

y

y

x

x

b

a

b

a

+

=

⋅

b

a

. Iloczyn zostaje scałkowany po

powierzchni; w ten sposób otrzymujemy następujące wartości N:

{

}

{

}



















−

−

+

−

−

∆

=

−

+

−

∆

=



)

)(

(

)

)(

(

4

1

)

(

)

(

4

1

3

1

2

3

1

3

3

2

)

(

12

2

2

3

2

3

2

)

(

11

x

x

x

x

y

y

y

y

N

x

x

y

y

N

e

e

,

[2.15]

które można wykorzystać w dalszych rozważaniach.

Przykład

Rys. 2.3a

Rys. 2.3b

W celu pokazania wykorzystania metody elementów skończonych, zostaną obliczone
potencjały dla obszaru jak na rys. 2.3a. Łatwo jest przewidzieć, że w punktach o potencjałach
równych

2

1

,V

V

wynik powinien wynosić

50

2

1

=

=

V

V

. Pierwszym krokiem, jest podzielenie

powierzchni na podobszary (czyli na trójkąty, spełniające podane wyżej założenia);
oczywiście zamiast podziału jak na rys. 2.3b, dozwolone są i inne kombinacje. Globalne
numeracja poszczególnych punktów jest rosnąca, poczynając od punktów, w których
potencjał nie jest znany, a kończąc na punktach z zadanym potencjałem. Numeracja lokalna
(cyfry w kółkach) jest przeciwna do kierunku ruchu wskazówek zegara, poczynając od punktu
leżącego w lewym, dolnym rogu:

Rys. 2.4

Dane zostały zapisane w poniższych tabelach:

Tabela 1

Punkt nr Współrzędne (x, y) Potencjał w (x, y)

1

1, 0

?

2

1, 1

?

3

2,1

100

4

0, 1

0

5

2, 0

100

6

0, 0

0

Tabela 2

Współrzędne

punktów lokalne

Element x

1

x

2

x

3

y

1

y

2

y

3

A

0

1

1

0

0

1

B

0

0

1

0

1

1

C

1

1

2

0

1

1

D

1

2

2

0

0

1

Korzystając z podanych wyprowadzeń, otrzymujemy, że

5

.

0

det

=

∆

. Obliczmy teraz różnice

współrzędnych dla elementu A:











=

−

−

=

−

=

−

1

1

0

1

2

3

1

2

3

x

x

x

x

x

x











=

−

=

−

−

=

−

0

1

1

2

1

1

3

3

2

y

y

y

y

y

y

.

[2.16]

Zatem poszczególne wartości N dla elementu A (ozn. przez

)

( A

ij

N

) będą przedstawiały się

następująco:

{

} {

}

{

}



















=

−

=

−

−

+

−

−

=

=

+

−

=

−

+

−

=

5

.

0

5

.

0

)

)(

(

)

)(

(

4

1

5

.

0

)

0

(

)

1

(

2

1

)

(

)

(

4

1

)

(

33

3

2

2

3

1

3

3

2

)

(

12

2

2

2

2

3

2

3

2

)

(

11

A

A

A

N

x

x

x

x

y

y

y

y

A

N

x

x

y

y

A

N



[2.17]

Dokonując analogicznych obliczeń dla elementów skończonych B, C i D, otrzymujemy
wyniki, które w naszym przypadku są identyczne z wynikami otrzymanymi dla elementu A:

[ ]

[ ]

[ ]

[ ]





















−

−

−

−

=

=

=

=

5

.

0

5

.

0

0

5

.

0

1

5

.

0

0

5

.

0

5

.

0

)

(

)

(

)

(

)

(

D

C

B

A

N

N

N

N

.

[2.18]

Dla przypadku ogólnego, macierze oczywiście nie będą sobie równe.
Ogółem mamy danych sześć punktów, a zatem macierz globalna powinna być rozmiarów
6x6. Patrząc na element A:

Rys. 2.5

można zauważyć, że macierz

)

(

]

[

A

N

można zapisać ze względu na punkty (6, 1, 2) jako:

[ ]





















=

22

21

26

12

11

16

62

61

66

)

(

N

N

N

N

N

N

N

N

N

N

A

,

[2.19]

podobnie i macierze B, C i D. Symbol

66

N oznacza teraz pozycję w macierzy globalnej

(szósta kolumna, szósty wiersz). Należy zauważyć jednak, że

66

N występuje dla dwóch

macierzy: dla elementów A i B (ponieważ punkt 6 jest punktem wspólnym dla elementów A i
B), więc będzie on sumą wartości macierzy

)

(

]

[

A

N

i

)

(

]

[

B

N

w tych punktach i wyniesie:

1

5

.

0

5

.

0

)

(

66

)

(

66

=

+

=

+

B

A

N

N

. Innymi słowy: wartości N o jednakowych indeksach w

macierzach A, B, C i D są sumowane, a następnie umieszczane w macierzy globalnej na
pozycji, którą określają te indeksy. Dla naszego przykładu, macierz globalna wyniesie więc:









































−

−

−

−

−

−

−

−

−

−

−

−

−

−

=

1

0

5

.

0

0

0

5

.

0

0

1

0

5

.

0

0

5

.

0

5

.

0

0

1

0

5

.

0

0

0

5

.

0

0

1

5

.

0

0

0

0

5

.

0

5

.

0

2

1

5

.

0

5

.

0

0

0

1

2

]

[N

.

[2.20]

Należy zaznaczyć, że zapisana w ten sposób macierz, niesie również informacje o tym, w
których punktach mamy dany potencjał, a w których nieznany; a przedstawić ją można jako:













−

−

=

2

1

1

2

]

[

ff

N

,













−

−

−

−

=

0

0

5

.

0

5

.

0

5

.

0

5

.

0

0

0

]

[

fp

N

,

























−

−

−

−

=

0

5

.

0

0

5

.

0

5

.

0

0

5

.

0

0

]

[

pf

N

,

























−

−

−

−

=

1

0

5

.

0

0

0

1

0

5

.

0

5

.

0

0

1

0

0

5

.

0

0

1

]

[

pp

N

[2.21]

gdzie:

f – punkty w których potencjał jest nieznany

p – punkty w których potencjał jest zadany

Teraz wykorzystana zostaje zasada ekstremum energii:

0

=

∂

∂

k

V

W

(gdzie k – numery punktów

dla globalnych oznaczeń). Równanie, możemy zapisać w następującej postaci:

0

]

[

]

[

]

[

]

[

]

[

]

[

]

]

[

]

[[

]

[

T

T

=





























∂

∂

=

∂

∂

p

f

pp

pf

fp

ff

p

f

k

f

k

V

V

N

N

N

N

V

V

V

V

W

.

[2.22]

Równanie jest prawdziwe tylko wtedy, gdy:

0

]

[

]

[

]]

][

[[

=

















p

f

fp

ff

V

V

N

N

, skąd po przekształceniu,

otrzymujemy:

]

][

[

]

[

]

[

1

p

fp

ff

f

V

N

N

V

−

−

=

[2.23]

co jak widać, jest wzorem na poszukiwane wartości potencjałów. Po podstawieniu danych,
otrzymujemy prosty układ równań:







=

+

−

=

−

50

2

50

2

2

1

2

1

V

V

V

V

,

[2.24]

którego rozwiązaniem jest

50

2

1

=

=

V

V

co jest zgodne z oczekiwaniami.

Metoda elementów skończonych, przezwycięża podstawowe trudności występujących w
metodzie różnic skończonych; można bowiem właściwie dopasować siatkę do nieregularnych
kształtów, a ponadto oferuje większą elastyczność i dokładność wyniku. Dla podanego
prostego przykładu (wykorzystującego równanie Laplace’a), obliczenia były stosunkowo
proste. Jednakże, w większości wypadków, zagadnienie staje się o wiele bardziej
skomplikowane, ze względu na bardziej złożony kształt powierzchni, ilość elementów
skończonych, czy też konieczność dokonania obliczeń dla równania Poissone’a.