9. Metoda różnic skończonych (MRS)
Idea metody różnic skończonych została zaproponowana już przez Carla
Fridricha Gaussa, lecz dopiero w erze komputerowej znalazła ona szersze
zastosowanie. Jest to ogólna metoda rozwiązywania równań różniczkowych.
Polega ona na zastąpieniu równania różniczkowego układem równań różnicowych
i sprowadzeniu problemu do rozwiązywania układu równań algebraicznych.
9.1. Ilorazy różnicowe jako przybliżenia pochodnych
9.1.1. Ilorazy dla pierwszej pochodnej
Pochodna funkcji f(x) w punkcie x0 jest zdefiniowana jako (rys. 9.1):
d f (x0 ) f (x0 + "x) - f (x0 )
= lim . (9.1)
"x0
dx "x
f(x)
sieczna
f(x0 + "x)
styczna
f(x0)
x0 x0 + "x x
Rys. 9.1. Geometryczna interpretacja pochodnej: dla "x 0 sieczna (czerwona linia) dąży do
stycznej (zielona linia)
Geometrycznie, jest to tangens kÄ…ta nachylenia stycznej do krzywej y = f(x)
w punkcie x0. Stąd dla  małych h mamy w przybliżeniu
d f (x0 ) f (x0 + h) - f (x0 )
H" , (9.2)
dx h
1
9. Metoda różnic skończonych (MRS)
tzn. postępujemy następująco (rys. 9.2a):
" obliczamy wartość funkcji w punkcie x0,
" obieramy pewien bliski punkt x0 + h i obliczamy w nim wartość funkcji,
" odejmujemy od niej wartość funkcji w punkcie x0,
" dzielimy przez przyrost h zmiennej x.
Tak obliczona wartość jest tym bliższa dokładnej wartości pochodnej, im mniejsza
wartość h, gdyż wtedy sieczna jest  bliższa stycznej. Wzór (9.2) nazywamy
ilorazem różnicowym dla pochodnej. Dla h > 0 jest to iloraz różnicowy wprzód.
Istnieje też wersja zwana ilorazem różnicowym wstecz (rys. 9.2b):
d f (x0 ) f (x0 ) - f (x0 - h)
H" , (9.3)
dx h
Powyższe ilorazy biorą się z założenia o liniowej zmienności funkcji f(x)
w pobliżu punktu x0 i dlatego są obarczone dość dużym błędem, rzędu O(h2).
Dokładniejsze wyrażenia na postać różnicową pochodnej można uzyskać,
rozpatrując interpolację kwadratową, która prowadzi do tzw. symetrycznego
ilorazu różnicowego dla pierwszej pochodnej (rys. 9.2c):
d f (x0 ) f (x0 + h) - f (x0 - h)
H" , (9.4)
dx 2h
który jest obarczony mniejszym błędem, rzędu O(h3). Jest to widoczne na rysunku
9.2  sieczna dla ilorazu różnicowego symetrycznego ma kąt nachylenia znacznie
bliższy kątowi nachylenia stycznej niż dla ilorazów różnicowych wprzód i wstecz.
a) b) c)
f(x) f(x) f(x)
f(x0 + h) f(x0 + h)
f(x0) f(x0) f(x0)
f(x0  h) f(x0  h)
h h h h
x0 x0 + h x x0  h x0 x x0  h x0 x0 + h x
Rys. 9.2. Różnice: a) różnica wprzód b) różnica wstecz, c) różnica symetryczna
Przykład 9.1. Obliczyć przybliżoną wartość pochodnej funkcji f(x) = x3 w punkcie x0 = 2,
przyjmujÄ…c h = 0,1.
Dokładna wartość pochodnej wynosi
dx3
.
= 3x2 = 3Å" 22 = 12
x=2
dx
x=2
Wartość obliczona ilorazem różnicowym wprzód:
dx3 (2 + 0,1)3 - 23
H" =12,61.
dx 0,1
x=2
2
9.1. Ilorazy różnicowe jako przybliżenia pochodnych
Iloraz różnicowy wstecz daje
dx3 23 - (2 - 0,1)3
H" =11,41,
dx 0,1
x=2
natomiast iloraz symetryczny prowadzi do wyniku
dx3 (2 + 0,1)3 - (2 - 0,1)3
H" = 12,01.
dx 2Å" 0,1
x=2
9.1.2. Iloraz różnicowy dla drugiej pochodnej
W przypadku drugiej pochodnej najwygodniej jest rozpocząć od rozwinięcia
funkcji f(x) w szereg Taylora w otoczeniu punktu x = x0:
1 1
2 2
f (x) = f (x0 ) + f '(x0 )(x - x0 ) + f (x0 )(x - x0 )2 + ... (9.5)
1! 2!
PrzyjmujÄ…c w nim x = x0 + h, gdzie h > 0, dostajemy
1 1 1
(3)
2 2
f (x0 + h) = f (x0 ) + f '(x0 )h + f (x0 )h2 + f (x0 )h3 + O(h4 ) ,
1! 2! 3!
a dla x = x0  h mamy
1 1 1
(3)
2 2
f (x0 - h) = f (x0 ) - f '(x0 )h + f (x0 )h2 - f (x0 )h3 + O(h4) .
1! 2! 3!
Po dodaniu stronami powyższych równań i przekształceniu dostajemy
d2 f (x) f (x0 - h) - 2 f (x0 ) + f (x0 + h)
2 2
f (x0 ) = H" . (9.6)
dx2 x=x0 h2
Dokładność tego wyrażenia jest rzędu O(h4), a więc dość dobra.
Przykład 9.2. Obliczyć przybliżoną wartość drugiej pochodnej funkcji f(x) = x3 w punkcie
x0 = 2 przyjmujÄ…c h = 0,1.
Dokładna wartość drugiej pochodnej wynosi
d2 x3
,
= 6x = 6 Å" 2 =12
dx2 x=2 x=2
natomiast wzór różnicowy daje
d2 x3 (2 - 0,1)3 - 2 Å" 23 + (2 + 0,1)3 .
H" =12
dx2 x=2 0,12
Wynik jest dokładny, gdyż czwarta i wyższe pochodne f(x) równe są zeru.
3
9. Metoda różnic skończonych (MRS)
9.2. Ogólna idea metody różnic skończonych
Rozwiązywanie równania różniczkowego metodą różnic skończonych można
podzielić na cztery etapy:
1. Zapisanie równania różniczkowego w postaci różnicowej, tj. na zastąpienie
pochodnych odpowiednimi ilorazami różnicowymi.
2. Nałożenie na obszar działania siatki węzłów, w których poszukiwane będą
wartości funkcji pola.
3. Ułożenie równania różnicowego dla każdego węzła, z uwzględnieniem
warunków brzegowych.
4. Rozwiązanie powstałego układu liniowych równań algebraicznych ze względu
na nieznane wartości węzłowe funkcji pola.
Postępowanie to zobrazowano na rys. 9.3.
START
du ui+1 - ui-1

dx 2h
Znajdz
postać różnicową
d2 u ui+1 - 2ui + ui-1
rozwiązywanego równania

dx2 h2
Nałóż na obszar działania siatkę
węzłów
Ułóż równanie różnicowe dla
każdego węzła
Wprowadz
warunki brzegowe
Rozwiąż powstały układ równań
STOP
Rys. 9.3. Schemat rozwiÄ…zywania zagadnienia polowego za pomocÄ… MRS
4
9.3. Zagadnienia jednowymiarowe
9.3. Zagadnienia jednowymiarowe
9.3.1. Równanie różnicowe
W przypadku zagadnień jednowymiarowych ograniczymy się do równania
różniczkowego o postaci
d2Åš(x) dÅš(x)
+ g(x) - º2Åš(x) = - f (x) , (9.7)
dx2 dx
gdzie Åš(x) jest poszukiwanÄ… funkcjÄ…, º  staÅ‚Ä…, g(x) i f(x)  znanymi funkcjami
zmiennej x. Wykorzystując wzory (9.6) oraz (9.4), możemy zapisać postać
różnicową równania
Åš(x + h) - 2Åš(x) + Åš(x - h) Åš(x + h) - Åš(x - h)
+ g(x) - º2Åš(x) = - f (x) ,
h2 2h
co po przekształceniu daje
1 1
[1- g(x)h]Åš(x - h) - (2 + º2h2)Åš(x) + [1+ g(x)h]Åš(x + h) = - f (x)h2 . (9.8)
2 2
Powyższe równanie jest spełnione w przybliżeniu dla każdego x i  niezbyt dużego
kroku h.
W dalszym ciągu załóżmy, że rozpatrywane równanie (9.7) określone jest na
odcinku (obszarze działania pola) xmin d" x d" xmax. Podzielmy ten odcinek na N
części i wprowadz
my na nim N + 1 węzłów (dla uproszczenia równomiernie
rozłożonych)  rys. 9.4.
Wartości węzłowe
Åš(x)
Åši 1
Åši Åši+1
Węzły
h h h h h
x
x0 = xmin x1 x2 xi 1 xi xi+1 xN 1 xN = xmax
Obszar działania pola
Rys. 9.4. Węzły nałożone na obszar działania pola
Mamy zatem krok siatki węzłów
xmax - xmin
h = (9.9a)
N
oraz
xi = xmin + ih, i = 0,1, 2,...N. (9.9b)
5
9. Metoda różnic skończonych (MRS)
Równanie różnicowe (9.8) zapisane dla dowolnego węzła x = xi przyjmuje postać
1 1
(1- gih)Åši-1 - (2 + º2h2)Åši + (1+ gih)Åši+1 = - fih2 , (9.9c)
2 2
gdzie wprowadzono oznaczenia Śi = Ś(xi), fi = f(xi), gi = g(xi). Takie równania
układa się dla wszystkich węzłów, przy czym dla węzłów brzegowych (x0 oraz xN)
należy uwzględnić dodatkowo warunki brzegowe.
Przykład 9.3. Znalez
ć rozkład potencjału elektrycznego w rozległej płycie dielektrycznej
o gruboÅ›ci a = 1 cm i przenikalnoÅ›ci wzglÄ™dnej µr = 5 (rys. 9.5), naÅ‚adowanej Å‚adunkiem
elektrycznym o gÄ™stoÅ›ci objÄ™toÅ›ciowej Á = Á0 (przyjąć Á0/µ0 = 100 V·cm 2). Jedna
powierzchnia płyty ma potencjał 1 V, na drugiej zadano zerowy warunek Neumanna.
V(x) = ?
dV
V = 1
= 0
dx
0 a x
Á, µr
Rys. 9.5. Naładowana płyta (przekrój)
Z równań elektrostatyki
(a)
" Å" D = Á, D = µrµ0E, E = -"V
wynika dla potencjału elektrycznego V równanie Poissona
Á
. (b)
"2V = -
µrµ0
W rozpatrywanym przypadku przyjmujemy układ współrzędnych z osią Ox prostopadłą do
powierzchni płyty. Z dala od krawędzi płyty pole zależy tylko od x, a stąd otrzymujemy
równanie
d2V (x) Á0 .
(c)
= -
dx2 µrµ0
Równanie to rozwiążemy za pomocą MRS z N + 1 węzłami rozlokowanymi równomiernie
na odcinku x = <0, 1>. Na podstawie wzoru (9.9c) mamy postać różnicową powyższego
równania (º = 0, g(x) = 0):
, (d)
Vi-1 - 2Vi +Vi+1 = -h2 fi
gdzie oznaczono fi = Á(xi)/(µ0µr) oraz h = a/N. W rozpatrywanym przypadku fi = 20 V·cm 2
niezależnie od i oraz h = 1/N cm. Równania takie układamy dla wszystkich wewnętrznych
węzłów, tj. dla i = 1, 2, 3, & , N  1. Dla węzła i = 0 postępujemy inaczej, gdyż jest to
węzeł brzegowy. Ponieważ zadano w nim warunek Dirichleta V = 1 V, mamy równanie
. (e)
V0 = 1
Jeszcze inaczej postępujemy dla węzła i = N, w którym zadano zerowy warunek brzegowy
Neumanna, tzn.
6
9.3. Zagadnienia jednowymiarowe
dV
. (f)
= 0
dx
x=1
Zapisujemy chwilowo równanie (d) dla tego węzła
, (g)
VN -1 - 2VN +VN +1 = -h2 fN
przy czym węzeł N + 1 jest węzłem fikcyjnym (brak go w podziale obszaru działania pola).
Należy więc wyeliminować go z tego równania. W tym celu zapisujemy warunek
Neumanna (f) w postaci różnicowej
VN +1 -VN -1 ,
(h)
= 0
2h
z której wynika, że VN+1 = VN 1. Uwzględniając to we wzorze (g), dostajemy
. (i)
2VN -1 - 2VN = -h2 fN
Równania (d), (e) oraz (i) dają w sumie N + 1 równań na N + 1 nieznanych wartości Vi:
V0 = 1,
Å„Å‚
ôÅ‚
- 2V1 +V2 = -h2 f1,
ôÅ‚V0
ôÅ‚M
ôÅ‚
ôÅ‚ (j)
òÅ‚V - 2Vi +Vi+1 = -h2 fi ,
i-1
ôÅ‚M
ôÅ‚
ôÅ‚VN -2 - 2VN -1 +VN = -h2 fN -1,
ôÅ‚
ôÅ‚2VN -1 - 2VN = -h2 fN .
ół
Powstały układ równań można rozwiązać używając komputera. Jeżeli przyjrzymy się
postaci macierzowej tego układu równań,
1 0 0 L 0 0 0 V0 1
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
ïÅ‚1 - 2 1 L 0 0 0 śł ïÅ‚ śł ïÅ‚
V1 - h2 f1 śł
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
ïÅ‚ - 2 L 0 0 0 V2 - h2 f2
śł ïÅ‚ śł ïÅ‚ śł
0 1
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
. (k)
M = M
ïÅ‚M M M O M M M śł ïÅ‚ śł ïÅ‚ śł
ïÅ‚0 0 0 L - 2 1 0 śł ïÅ‚VN -2 śł ïÅ‚- h2 fN -2 śł
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
śł ïÅ‚- h2 fN -1
śł
ïÅ‚0 0 0 L 1 - 2 1 śł ïÅ‚V -1
N
ïÅ‚0 0 0 L 0 2 - 2śł ïÅ‚ śł ïÅ‚
VN - h2 fN śł
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚
to zauważymy, że macierz główna jest rzadka (w jednym wierszu występują prawie same
zera), pasmowa (wszystkie niezerowe elementy znajdują się w paśmie biegnącym wzdłuż
przekątnej) oraz symetryczna (prawie  odchyłki pojawiają się w równaniach dla węzłów
brzegowych). Są to ogólne cechy macierzy głównych układów równań generowanych przez
MRS.
Przyjmując dalej N = 5, dostajemy następujące wyniki:
V0 = 1, V1 = 4,6, V2 = 7,4, V3 = 9,4, V4 = 10,6, V5 = 11, (l)
które zobrazowano na rys. 9.6. Rozpatrywane zagadnienie ma rozwiązanie dokładne
Vd(x) =  10x2 + 20x + 1. Można sprawdzić, że jego wartości węzłowe są dokładnie równe
7
9. Metoda różnic skończonych (MRS)
podanym wyżej. Oznacza to, że w tym przypadku MRS generuje dokładne wartości
węzłowe, ale nie należy sądzić, że jest tak zawsze. Raczej przeciwnie  na ogół wartości
węzłowe są tylko zbliżone do wartości dokładnych.
Rys. 9.6. Rozwiązanie zagadnienia: MRS (punkty), rozwiązanie dokładne (linia ciągła)
9.3.2. Wprowadzanie warunków brzegowych
Z przykładu 9.3 wynika sposób wprowadzania warunków brzegowych do
równań MRS.
Warunek brzegowy Dirichleta
W przypadku warunków brzegowych typu Dirichleta jest to wyjątkowo
nieskomplikowane. Istotnie, dla warunku brzegowego Dirichleta w k-tym węz
le
(zwykle k = 0 lub N + 1) o postaci
Åš(xk ) = A , (9.10)
gdzie A znaną funkcją (tutaj stałą), zamiast równania różnicowego układa się
równanie o postaci
Åšk = A , (9.11)
W realizacji macierzowej odpowiada to wyzerowaniu k-tego wiersza macierzy
głównej układu równań, wpisaniu na przekątną 1 oraz wpisaniu na k-tą pozycję
wektora wyrazów wolnych wartości A.
Warunek brzegowy Neumanna
W przypadku warunku brzegowego Neumanna o postaci
dÅš
= B (9.12)
dn
x= x0
sytuacja jest również nieskomplikowana. Najpierw układamy równanie różnicowe
dla danego węzła brzegowego. Dla węzła 0 będzie ono miało postać (9.9c) z i = 0:
2
1 1
(1- g0h)Åš-1 - (2 + º h2 )Åš0 + (1+ g0h)Åš1 = - f0h2 .
2 2
8
9.3. Zagadnienia jednowymiarowe
Występuje w nim fikcyjny węzeł  1, który wyeliminujemy za pomocą warunku
Neumanna. Warunek ten w postaci różnicowej przyjmuje postać (rys. 9.7)
Åš-1 -Åš1
= B Ò! Åš-1 = Åš1 + 2hB .
2h
Uwzględniając to w poprzednim równaniu, dostajemy
2
1
- (2 + º h2 )Åš0 + 2Åš1 = - f0h2 - (1- g0h)2hB . (9.13a)
2
W równaniu tym nie występuje już fikcyjny węzeł i =  1. Analogicznie, dla
warunku Neumanna na brzegu x = xN dostaniemy
2
1
2ÅšN -1 - (2 + º h2 )ÅšN = - fN h2 - (1+ gN h)2hB . (9.13b)
2
Åš(x)
Åš 1
Åš0
Åš1
Węzeł
fikcyjny
Åš
h h
x 1 x0 = xmin x1 x
Rys. 9.7. Wprowadzanie warunku Neumanna lub Hankela
Warunek brzegowy Hankela
Również wprowadzenie warunku brzegowego
dÅš
= Ä…(Åš - ²) (9.14)
dn
x=x0
nie nastręcza trudności. Zapisany w postaci różnicowej
Åš-1 -Åš1
= Ä…(Åš0 - ²) Ò! Åš-1 = Åš1 + 2hÄ…(Åš0 - ²) ,
2h
pozwala na wyeliminowanie z równania różnicowego (9.9c) fikcyjnego węzła  1.
Ostatecznie otrzymuje siÄ™
2
1 1
[2hÄ…(1- g0h) - (2 + º h2)]Åš0 + 2Åš1 = - f0h2 + 2hÄ…²(1- g0h) . (9.15a)
2 2
Podobnie, dla węzła na brzegu x = xN dostaniemy
1 1
2ÅšN -1 +[2hÄ…(1+ gN h) - (2 + º2h2 )]ÅšN = - fN h2 + 2hÄ…²(1+ gN h) . (9.15b)
2 2
9
9. Metoda różnic skończonych (MRS)
Przykład 9.4. Wyznaczyć rozkład potencjału w obszarze rezystora cylindrycznego
o wysokości H = 1 cm i promieniach R1 = 0,2 mm i R2 = 2 mm, wypełnionego materiałem
o konduktywności ł = 1 S/m, jeżeli przez wewnętrzną okładkę wpływa prąd I = 1 A
(założyć równomierną gęstość), a zewnętrzna okładka ma potencjał równy zeru  rys. 9.8.
Rozwiązanie MRS porównać z dokładnym.
z
V(r) = ?
à = Ã0
V = 0
R1 R2 r
Á, µr
I
Rys. 9.8. Rezystor cylindryczny
Wychodząc z równań pola przepływowego
" Å" J = 0, J = Å‚E, E = -"V , (a)
dostajemy dla potencjału elektrycznego V równanie Laplace a:
"2V = 0, (b)
które w rozpatrywanych warunkach zależności V jedynie od odległości r od osi rezystora
redukuje siÄ™ do postaci
d2V (r) 1 dV (r)
(c)
+ = 0.
dr2 r dr
Wprowadzamy siatkę N węzłów rozłożonych równomiernie wzdłuż odcinka r = [R1, R2]:
R2 - R1
(d)
ri = R1 + ih, h = .
N
Równanie (c) zapisujemy w postaci różnicowej
Vi-1 - 2Vi +Vi+1 1 Vi+1 -Vi-1
(e)
+ = 0,
h2 ri 2h
a przemnożeniu przez h2 dostajemy
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚
h h
(f)
ìÅ‚ ÷Å‚ - 2Vi +
ìÅ‚1- 2ri ÷Å‚Vi-1 ìÅ‚ 2ri ÷Å‚ = 0.
ìÅ‚1+ ÷Å‚Vi+1
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
10
9.4. Zagadnienia dwuwymiarowe w układzie kartezjańskim
Równania takie układamy dla węzłów i = 1, 2, & , N  1. W węz
le i = 0 (tj. r = R1) musimy
uwzględnić warunek brzegowy
dV I dV I
(g)
Å‚ = Ò! = B, B = .
dn 2Ä„R1H dn 2Ä„Å‚HR1
r =R1 r =R1
Jest to warunek Neumanna; na podstawie wzoru (9.13a) otrzymujemy zatem
ëÅ‚ öÅ‚
h
(h)
ìÅ‚
- 2V0 + 2V1 = -2BhìÅ‚1- ÷Å‚
2R1 ÷Å‚.
íÅ‚ Å‚Å‚
Z kolei w węz
le i = N (tj. r = R2) mamy warunek Dirichleta V = 0, a stąd równanie
przyjmuje postać
(i)
VN = 0.
Rozwiązanie uzyskane za pomocą MRS dla różnej liczby węzłów przedstawiono na rys.
9.9. Naniesiono na nim również rozwiązanie dokładne, równe
I R2
(j)
V (r) = ln .
2Ä„Å‚H r
Jak wynika z otrzymanych wyników, w tym przypadku dla zbyt małej liczby węzłów
otrzymuje się wyniki bardzo odległe od rozwiązania dokładnego. Można to co prawda
nieco poprawić, stosując pewne modyfikacje, ale pokazuje to, że dokładność metod
numerycznych, w szczególności  MRS, może pozostawiać wiele do życzenia.
a) b) c)
Rys. 9.9. Wyniki obliczeń: a) N = 5, b) N = 20, c) N = 40.
9.4. Zagadnienia dwuwymiarowe w układzie kartezjańskim
9.4.1. Schemat pięciopunktowy
Rozważmy dwuwymiarowe niejednorodne równanie Helmholtza w układzie
współrzędnych kartezjańskich:
"2Åš(x, y) "2Åš(x, y)
2
+ - “ Åš(x, y) = - f (x, y) . (9.16)
"x2 "y2
gdzie “ oraz f sÄ… znanymi staÅ‚ymi lub funkcjami. Rozpatrywany przypadek
obejmuje także równanie Poissona (gdy “2 = 0) oraz Laplace a (gdy “2 = 0 i f = 0).
11
9. Metoda różnic skończonych (MRS)
Wprowadz
my prostokątną siatkę węzłów jak na rys. 9.10. Węzły znajdują się na
przecięciu wierszy i kolumn, przy czym mamy I + 1 kolumn numerowanych
indeksem i = 0, 1, & , I oraz J + 1 wierszy numerowanych indeksem j = 0, 1, & , J.
Węzeł o współrzędnych dyskretnych (i, j) ma współrzędne kartezjańskie (xi, yj),
przy czym
xi = xmin + ihx , y = ymin + jhy , (9.17a)
j
xmax - xmin ymax - ymin
hx = , hy = . (9.17b)
I J
Postać różnicową równania (9.16) dla węzła znajdującego się na przecięciu i-tej
kolumny i j-tego wiersza zapiszemy jako
Åši-1, j - 2Åši, j + Åši+1, j Åši, j-1 - 2Åši, j + Åši, j+1 2
+ - “ Åši, j = - fi, j , (9.17a)
2 2
hx hy
gdzie Śi,j = Ś(xi, yj), fi,j = f(xi, yj). Jest to równanie różnicowe ułożone wg. tzw.
schematu pięciopunktowego. Oznaczając ś = hx/hy, można je zapisać jako
2 2 2 2 2 2
Åši-1, j + Åši+1, j + Å› Åši, j-1 + Å› Åši, j+1 - (2 + 2Å› + “ hx )Åši, j = -hx fi, j . (9.18b)
Schemat taki można stosować w obszarach płaskich o różnych kształtach, jednak
najlepsze efekty daje dla prostokąta lub obszaru dającego się rozłożyć na sumę
prostokątów (np. kątownik, prostokąt z prostokątnym otworem  rys. 9.11). Chodzi
o to, aby węzły siatki wypadły na brzegu obszaru. Jeśli tak nie jest, pojawiają się
kłopoty z uwzględnieniem geometrii brzegu oraz warunków brzegowych.
y j
ymax
j = J
j (i, j)
yj
hy
hx
i
j = 0
ymin
i = I
i = 0 i
xmin xi xmax x
Rys. 9.10. Siatka dwuwymiarowa w układzie kartezjańskim nałożona na prostokąt
12
9.4. Zagadnienia dwuwymiarowe w układzie kartezjańskim
a) b) c)
Rys. 9.11. Obszary dwuwymiarowe z nałożoną siatką węzłów: a) kątownik, b) prostokąt z otworem
prostokątnym, c) kształt nieregularny  trudny do uwzględnienia w MRS
9.4.2. Wprowadzanie warunków brzegowych
W przypadku obszarów o brzegach wypadających na liniach węzłów warunki
brzegowe do równań MRS wprowadza się zasadniczo tak samo jak to omówiono w
punkcie 9.3.2. Jeżeli obszar ma nieregularny kształt, wprowadzanie warunków jest
kłopotliwe i nie jest tutaj omawiane.
Warunek Dirichleta
Jeżeli w węz
le (i, j) znana jest wartość funkcji pola
Åš = A , (9.19)
to zamiast równania (9.18b) należy zapisać równanie
Åši, j = A . (9.20)
Jest to więc dokładnie tak, jak w przypadku jednowymiarowym.
Warunek Neumanna
Jeżeli w węz
le brzegowym (i, j) zadany jest warunek Neumanna
"Åš
= B , (9.21)
"n
i, j
to mogą zajść przypadki pokazane na rys. 9.12.
a) b) c) d)
(i, j+1)
(i, j+1) (i, j+1) (i, j+1)
n n
hy
hy
hy hy
(i 1, j) hx (i, j)
hx hx hx hx (i+1, j)
(i, j) (i+1, j) (i 1, j) (i, j) (i, j) (i 1, j) (i+1, j)
(i 1, j) (i+1, j)
hy hy
n
n
(i, j 1) (i, j 1) hx (i, j 1)
(i, j 1)
Rys. 9.12. Możliwe przypadki przy wprowadzaniu warunku Neumanna lub Hankela w siatce XY
13
9. Metoda różnic skończonych (MRS)
Wszystkie one są łatwe do uwzględnienia, dlatego tylko pierwszy zostanie
omówiony. Postać różnicowa warunku (9.21) to wtedy
Åši-1, j -Åši+1, j
= B Ò! Åši-1, j = Åši+1, j + 2hxB . (9.22)
2hx
Uzyskana zależność pozwala na wyeliminowanie fikcyjnego węzła (i 1, j)
z równania różnicowego (9.18b), co prowadzi do zależności
2 2 2 2 2 2
2Åši+1, j + Å› Åši, j-1 + Å› Åši, j+1 - (2 + 2Å› Åši, j + “ hx )Åši, j = -hx fi, j - 2hxB . (9.23)
Warunek Hankela
Jeżeli w węz
le brzegowym (i, j) zadany jest warunek Hankela
"Åš
= Ä…(Åš - ²) , (9.24)
"n
i, j
to mogą zajść przypadki pokazane na rys. 9.12. Wszystkie je uwzględnia się
analogicznie, dlatego dalej rozpatrujemy tylko pierwszy z nich. Z postaci
różnicowej warunku brzegowego mamy
Åši-1, j -Åši+1, j
= Ä…(Åši, j - ²) Ò! Åši-1, j = Åši+1, j + 2hxÄ…(Åši, j - ²) , (9.25)
2hx
co pozwala wyeliminować z równania (9.18b) fikcyjny węzeł (i 1, j) i daje
2 2 2 2 2 2
2Åši+1, j + Å› Åši, j-1 + Å› Åši, j+1 - (2 + 2Å› + “ hx - 2Ä…hx )Åši, j = -hx fi, j + 2Ä…²hx . (9.26)
Przykład 9.5. W środowisku nieprzewodzącym umieszczono prostokątną płytkę
o przewodności ł, grubości h, długości 2a i szerokości 2b. Do przeciwległych boków płytki
przyłożono napięcie 2U, przy czym elektrody umieszczone są centralnie i mają długość 2d
(rys. 9.13a). Wyznaczyć rozkład potencjału, przepływający prąd, straty mocy oraz
rezystancję płytki. W obliczeniach numerycznych przyjąć ł = 1000 S/m, h = 1 mm,
a = 1 cm, b = 1 cm, d = 0,5 cm, U = 100 V.
a) b)
y
y
b
2b Å‚
"V
 U U
2d = 0
"n
V = 0
x
I I
V = U
"2V = 0 d
x
2a a
Rys. 9.13. Prostokątna płytka z prądem: a) obszar oryginalny, b) obszar przyjęty do analizy ze
względu na symetrię
14
9.4. Zagadnienia dwuwymiarowe w układzie kartezjańskim
Ze względu na symetrię obszaru działania pola i warunków brzegowych analizę
można ograniczyć do ćwiartki obszaru (rys. 9.13b). Z równań pola przepływowego
" Å" J = 0, J = Å‚E, E = -"V (a)
wynika dla potencjału V równanie Laplace a
"2V "2V
(b)
+ = 0.
"x2 "y2
Warunki brzegowe:
- elektroda: warunek Dirichleta V = U,
- ścianki płytki (y = b, x = a oprócz elektrody): zerowy warunek Neumanna "V/"n = 0
wynikający z braku możliwości przepływania prądu przez te ścianki,
- linia x = 0: zerowy warunek Dirichleta V = 0 wynikający z symetrii przyłożenia elektrod
względem osi Oy,
- linia y = 0: zerowy warunek Neumanna "V/"n wynikający z symetrii przyłożenia
elektrod względem osi Ox.
Na rozpatrywany obszar nakładamy siatkę węzłów wg rys. 9.10. Równanie różnicowe
odpowiadajÄ…ce równaniu (b) ma postać (9.18b) z “ = f = 0. Po prostym przeksztaÅ‚ceniu
zapiszemy je jako
2 2
Vi-1, j +Vi+1, j + Å› Vi, j -1 + Å› Vi, j +1
(c)
Vi, j = .
2
2 + 2Å›
Jest ono poprawne w każdym węz
le, z tym, że należy uwzględnić warunki brzegowe.
Ostatecznie mamy:
- w węzłach wewnętrznych (i, j), tj. dla i = 1, 2, & , I  1 oraz j = 1, 2, & , J  1 spełnione
jest równanie (c),
- na linii x = 0, tj. w węzłach (0, j), dla j = 0, 1, & , J zamiast równania (c) mamy
(d)
V0, j = 0,
- na linii y = 0 oprócz naroży, tj. w węzłach (i, 0) dla i = 1, 2, & , I  1 należy uwzględnić
zerowy warunek Neumanna, w związku z czym równanie (c) przyjmuje postać
2
Vi-1,0 +Vi+1,0 + 2Å› Vi,1
(e)
Vi,0 = ,
2
2 + 2Å›
- na linii y = b oprócz naroży, tj. w węzłach (i, J) dla i = 1, 2, & , I  1 należy uwzględnić
zerowy warunek Neumanna, w związku z czym równanie (c) przyjmuje postać
2
Vi-1,J + Vi+1,J + 2Å› Vi,J -1
(f)
Vi,J = ,
2
2 + 2Å›
- na elektrodzie V = U, tj. w węzłach (I, j) dla j = 0, 1, & , J', gdzie J' = [d/b] jest numerem
wiersza węzłów do którego sięga elektroda, mamy warunek Dirichleta
(g)
VI , j = U,
- na pozostałej części linii x = a (poza elektrodą) oprócz naroża, tj. w węzłach (I, j) dla
j = J' + 1, J' + 2, & , J  1 mamy zerowy warunek Neumanna, stąd równanie (c)
przyjmuje postać
15
9. Metoda różnic skończonych (MRS)
2 2
2VI -1, j + Å› VI , j -1 + Å› VI , j +1
(h)
VI , j = ,
2
2 + 2Å›
- w narożu (a, b), tj. w węz
le (I, J), który należy zarówno do brzegu x = a, jak i y = b,
należy uwzględnić zerowe warunki Neumanna w obydwu kierunkach; równanie (c) daje
2
2VI -1,J + 2Å› VI ,J -1
(i)
VI ,J = ,
2
2 + 2Å›
przy czym jeżeli J' = J, to pierwszeństwo ma równanie (g).
W ten sposób ułożyliśmy dla każdego węzła równanie. Wszystkich równań jest
(I + 1)(J + 1). Ze względu na to, że macierz główna jest rzadka i ma normę mniejszą od
jedności można zastosować iteracyjną metodę rozwiązywania takiego układu równań.
Polega to na przyjęciu dowolnych wartości potencjału w poszczególnych węzłach
i wyznaczeniu nowych wartoÅ›ci ze wzorów (c)÷(i). Proces ten należy powtarzać dopóki
nowe wartości węzłowe nie będą dostatecznie bliskie poprzednim.
W celu przyspieszenia zbieżności można zastosować metodę relaksacji (zob. punkt
3.2.5), która polega na tym, że nowe wartości węzłowe obliczane są jako
(j)
Vi(nowe) = Vi(stare) + "Vi, j ,
, j , j
gdzie
(k)
"Vi, j = É[(wzór (c) - (i)) -Vi(stare)],
, j
przy czym É jest parametrem relaksacji z zakresu (0, 2). Dla É = 1 uzyskuje siÄ™ zwyczajnÄ…
metodÄ™ iteracji, która jest dość wolno zbieżna (okoÅ‚o 300 iteracji dla ´max = 10 6 przy siatce
10×10, tzn. 11×11 = 121 wÄ™złów). W rozpatrywanym przypadku wartość É równa okoÅ‚o
1,6 zapewnia znacznie szybszą zbieżność (80 iteracji przy tej samej dokładności).
Jako kryterium zakończenia obliczeń można przyjąć warunek
I J I J
2 (nowe)
(l)
d" ´max )2 ,
"""Vi, j ""(Vi, j
i=0 j =0 i=0 j=0
gdzie ´max jest maksymalnym bÅ‚Ä™dem wzglÄ™dnym (np. 10 5). Powyższy warunek oznacza,
że poprawki są bardzo małe w porównaniu z wartościami węzłowymi, więc proces
rozwiązywania układu równań można uznać za zakończony.
Wyznaczanie prądu. Z dwóch ostatnich równań (a) wynika, że składowa normalna
gęstości prądu na brzegu płytki wynosi
"V
(m)
Jn = Å‚En = -Å‚ .
"n
W węz
le (0, j) na brzegu x = 0 postać różnicowa powyższego wyrażenia to
V0, j -V1, j
(n)
Jn 0, j = -Å‚ .
hx
Całkując po powierzchni brzegowej, dostaniemy całkowity prąd
b
(o)
Iwy = Å" dS = hd y
n
+"J +"J .
x=0 0
16
9.4. Zagadnienia dwuwymiarowe w układzie kartezjańskim
Całkowanie należy wykonać numerycznie, np. metodą trapezów:
J
Jn j-1 + Jn j
(p)
Iwy = h hy.
"
2
j=1
Jest to prąd płynący w górnej połowie płytki. Całkowity prąd jest równy I = 2Iwy.
Rezystancja płytki wynosi R = (2U)/I = U/Iwy, a straty mocy P = RI2. Przykładowe
wyniki obliczeń zaprezentowano na rys. 9.14.
a) b)
5×5 wÄ™złów, Iwy = 87,7 A, R = 87,7 &! 11×11 wÄ™złów, Iwy = 84,3 A, R = 84,3 &!
c)
Rys. 9.14. Obraz pola w płytce kwadratowej
z elektrodą zajmującą połowę ścianki
(czerwona linia); rysunki przedsta-
wiają ćwiartkę obszaru; zaznaczono
węzły (kropki), wartości węzłowe
potencjału (drobna czcionka), linie
ekwipotencjalne (od 0 do 100 V co
5 V), linie prądu (niebieskie strzałki);
kolor wypełnienia wskazuje na
gęstość wydzielanej mocy: jasne
obszary  duża moc, ciemne  mała:
a) dyskretyzacja zgrubna, ewidentnie
niedokładna, b) dyskretyzacja śred-
nia, widoczne załamania linii ekwi-
potencjalnych wskazują na zbyt małą
dokładność, c) dyskretyzacja drobna,
wydaje się dość dokładna
21×21 wÄ™złów, Iwy = 83,1 A, R = 83,1 &!
17
9. Metoda różnic skończonych (MRS)
9.5. Inne typowe siatki węzłów
9.5.1. Schemat pięciopunktowy w układzie biegunowym
Dla obszarów działania pola w kształcie koła, pierścienia, sektora koła lub
wycinka pierścienia wygodniejszy od kartezjańskiego jest układ współrzędnych
biegunowych, w którym laplasjan ma nieco inną postać:
"2Ś(r,Ć) 1 "Ś(r,Ć) 1 "2Ś(r,Ć)
"2Åš = + + . (9.27)
2
"r2 r "r r "Ć2
Na obszar działania pola nakładamy siatkę biegunową ( pajęczynę ) zbudowaną z
I + 1 łuków okręgów numerowanych indeksem i = 0, 1, & , I oraz z pęku J + 1
półprostych wychodzących ze środka układu współrzędnych i numerowanych
indeksem j = 0, 1, & , J  rys. 9.15. Dzięki temu węzły wypadają na granicy
obszaru. Dla ustalenia uwagi przypuśćmy, że obszarem działania jest pola jest
wycinek pierścienia o promieniach wewnętrznym rmin i zewnętrznym rmax
i rozpiętości kątowej od Ćmin do Ćmax. Wtedy skoki siatki wynoszą odpowiednio
rmax - rmin Ćmax - Ćmin
hr = , hĆ = , (9.28a)
I J
zaś węzeł o współrzędnych dyskretnych (i, j) ma współrzędne biegunowe
ri = rmin + ihr , Ćj = Ćmin + jhĆ . (9.28b)
Wtedy postać różnicowa laplasjanu dla węzła (i, j) to
Åši-1, j - 2Åši, j + Åši+1, j Åši+1, j -Åši-1, j Åši, j -1 - 2Åši, j + Åši, j+1
"2Åš H" + + . (9.29)
2 2
i, j
hr 2hrri ri2hĆ
y
hr
hĆ
(i, j)
j = J
j i
j
i = I
Ćmax
i
j = 0
Ćj
i = 0
Ćmin
rmin ri rmax x
Rys. 9.15. Siatka dwuwymiarowa w układzie biegunowym nałożona na wycinek pierścienia
18
9.5. Inne typowe siatki węzłów
Jeżeli mamy do czynienia z pełnym zakresem zmienności kąta Ć, tzn. Ćmin = 0
i Ćmax = 2Ą, to linie j = 0 i j = J należy utożsamić. Jeżeli z kolei rmin = 0, to okrąg
i = 0 degeneruje się do punktu. Degeneracja ta uniemożliwia skorzystanie ze
pięciopunktowego schematu różnicowego, gdyż węzeł środkowy sąsiaduje wtedy
nie z czterema innymi węzłami w kierunkach zmian r i Ć, lecz z wieloma węzłami
położonymi na okręgu i = 1. To wymusza szczególne traktowanie węzła
środkowego. W przypadku koła można pokazać, że
J -1
1, j
"Åš - JÅši=0
j=0
"2Åš H" , (9.30)
r=0
2
i=0 hr
natomiast dla sektora sytuacja zależy od rodzaju warunku brzegowego zadanego
w tym węz
le i kąta rozwarcia, lecz nie będziemy tego szczegółowo rozpatrywać.
a) b) c)
j = 1
j = 0
j = J i = 0
j = J 1
Rys. 9.16. Przypadku siatki biegunowej wymagające szczególnego traktowania: a) pierścień  pełny
zakres kąta powoduje utożsamienie linii j = 0 i j = J, b) sektor koła  zerowy promień
wewnętrzny rmin = 0  degeneracja łuku okręgu i = 0, w węz
le środkowym konieczne jest
specyficzne postępowanie zależne od warunków brzegowych, c) pełne koło  następuje
zarówno degeneracja środkowego okręgu do punktu oraz nałożenie linii j = 0 i j = J,
w punkcie środkowym można stosować wzór (9.30)
Przykład 9.6. Obliczyć rozkład gęstości prądu w przewodzie maszyny elektrycznej
ułożonym żłobku wyciętym w magnetycznym rdzeniu (rys. 9.17). Parametry pręta:
konduktywność Å‚, przenikalność magnetyczna µ, promieÅ„ przekroju poprzecznego a.
Parametry wymuszenia: prÄ…d sinusoidalny o wartoÅ›ci skutecznej natężenia I i pulsacji É.
Obliczyć rezystancjÄ™ prÄ™ta. W obliczeniach numerycznych przyjąć ɵła2/2 = 1.
Fe
(µr >> 1)
Cu
a
Å‚, µ
I, É
Rys. 9.17. Przewody w ułożone w rdzeniu; po prawej powiększenie otoczenia jednego przewodu
19
9. Metoda różnic skończonych (MRS)
Z uwagi na sinusoidalną zależność od czasu można stosować zapis zespolony.
Z powodu geometrii stosujemy biegunowy układ współrzędnych ze środkiem
pokrywającym się ze środkiem przekroju poprzecznego przewodu. Ponieważ długość
przewodu jest znacznie większa w porównaniu z jego wymiarami poprzecznymi, wektor
gęstości prądu ma tylko składową wzdłużną (z-ową). Z równań pola
(a)
" × H = J , "× E = - jɵH, J = Å‚E, " Å" J = 0
otrzymuje się równanie dla zespolonego wektora gęstości prądu
(b)
"2J - jɵłJ = 0.
Wprowadzając względną zespoloną gęstość prądu J' = J/J0, gdzie J0 = I/(Ąa2) jest średnią
gęstością prądu, dostajemy równanie dla składowej z-owej J':
2 2 2
"2 J (r,Ć) 1 "J (r,Ć) 1 "2 J (r,Ć) º2
(c)
2
+ + - j2 J = 0,
"r2 r "r r2 "Ć2 a2
gdzie
ɵł
(d)
º = a
2
jest tzw. względnym współczynnikiem wypierania. Symetria obszaru działania pola
pozwala na rozpatrywanie tylko górnej połowy przekroju pręta (rys. 9.18a).
a) b)
"J 2
y
y
= 0
2
"n "J Ä„
= jº2
j = N1
"n Ä…a
Ä… j = N
j = 0
x i = 0 i = M x
Rys. 9.18. Obszar przyjęty do analizy (a) oraz jego siatka biegunowa (b), w której pokazano
dodatkowo sposób traktowania węzła i = 0
Warunki brzegowe:
- na styku miedz
-żelazo zadany jest zerowy warunek Neumanna, wynikający z ciągłości
składowej stycznej natężenia pola magnetycznego
2
"J
(e)
= 0,
"n
- na styku miedz
-powietrze dokładnego warunku brzegowego nie znamy, lecz przyjmuje
się, że jest to stały warunek Neumanna; wtedy z prawa przepływu wyznaczamy
2
"J Ä„
(f)
= jº2 ,
"n Ä…a
gdzie 2Ä… jest kÄ…tem rozwarcia na powietrze (w rozpatrywanym przypadku Ä… = Ä„/2).
- na linii symetrii (y = 0) mamy zerowy warunek Neumanna (e), który wynika z tego, że
obraz pola w dolnej połowie jest symetrycznym odbiciem obrazu w górnej połowie.
20
9.5. Inne typowe siatki węzłów
Na analizowany obszar nakładamy siatkę biegunową (rys. 9.18b). Postać różnicową
równania (c) dla węzła (i, j) uzyskamy poprzez zależność (9.29). Po prostych
przekształceniach otrzymujemy
ëÅ‚
ëÅ‚1- 1 1 1 1 1 hr ÷Å‚J2 i,
öÅ‚J2 + ëÅ‚1+ öÅ‚J2 + J2 i, + J2 i, - 2ìÅ‚1+ + j 2 º2 öÅ‚ = 0. (g)
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚
i-1, j i+1, j
2 2 2
2i 2i i2hĆ j -1 i2hĆ j+1 ìÅ‚ i2hĆ a2 ÷Å‚ j
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
íÅ‚ Å‚Å‚
Równanie takie jest układane dla każdego węzła (i, j) siatki obszaru, i = 1, 2, & M, j = 0, 1,
& , M, przy czym w węzłach brzegowych należy uwzględnić warunki brzegowe. Wszystkie
one są typu Neumanna, co powoduje, że:
- jeżeli j = 0, to w miejsce J'i,j 1 należy wstawić J'i,1,
- jeżeli j = N, to w miejsce J'i,j+1 należy wstawić J'i,N 1,
- jeżeli i = M, to w miejsce J'i+1,j należy wstawić J'M 1,j; jeżeli dodatkowo j d" N1 (rys.
9.18b), to należy uwzględnić niezerowy warunek Neumanna (f), co powoduje, że po
prawej stronie równania (g) w miejsce zera należy wstawić
Ä„ 1
ëÅ‚ öÅ‚
(h)
- jº2 2hr 1+ .
ìÅ‚ ÷Å‚
Ä…a 2i
íÅ‚ Å‚Å‚
Punkt i = 0 wymaga szczególnego traktowania (nie można użyć równania (g)). Zauważmy,
że możemy w nim zastosować schemat pięciopunktowy o postaci (9.18b)  rys. 9.18b.
Uwzględniając dodatkowo zerwy warunek Neumanna, dostajemy równanie
2
ëÅ‚
hr öÅ‚
(i)
2 1,0 2 1,N 2 1,N / 2 ìÅ‚ 2 0,
J + J + 2J - 2ìÅ‚2 + j º2 ÷Å‚J j = 0.
a2 ÷Å‚
íÅ‚ Å‚Å‚
Tak powstały układ równań rozwiązujemy ze względu na wartości węzłowe J'i,j.
Obliczanie rezystancji. Rezystancję można obliczyć jako R = P/I2, gdzie P jest mocą
czynną wydzielaną w obszarze pręta. Moc tę można obliczyć na dwa sposoby: albo
sumujÄ…c elementarne moce wydzielane w obszarze przewodu, albo obliczajÄ…c moc
wnikającą do przewodu przez jego brzeg. Ten drugi sposób jest łatwiejszy
w implementacji. Moc pozorna wnikająca do przewodu wyraża się wzorem
"
(j)
S = - × H Å" dS.
+"+"E
pow.
boczna
Po przekształceniach dochodzi się do zależności na impedancję wewnętrzną przewodu
o przekroju kołowym odniesioną do rezystancji dla prądu stałego:
2 "
1 S j "J
(k)
2
Z = =
2
+"J2 dl,
RDC I ɵłS "n
C
która jest poprawna dla dowolnego przewodu o przekroju ograniczonym krzywą C i polu
tego przekroju równym S. W rozpatrywanym przypadku dostaniemy
Ä… Ä… " Ä…
2 "
j "J j Ä„ 1
ëÅ‚ öÅ‚
(l)
2
Z =
+"J2 a dĆ = 2+"J2 ìÅ‚ jº2 Ä…a ÷Å‚ adĆ = Ä… +"J2 dĆ,
2º2Ä„ "n 2º2Ä„
íÅ‚ Å‚Å‚
-Ä… 0 0
gdzie skorzystano z warunku Neumanna (f). Okazuje się więc, że impedancja jest średnią
wartością zespolonej względnej gęstości prądu na powierzchni przewodu stykającą się
z powietrzem. Całkowanie wykonujemy metodą trapezów, dostając zależność
21
9. Metoda różnic skończonych (MRS)
N1
2 M , j -1 2 M , j
J + J
1
(m)
2
Z H" .
"
N1 j =1 2
Wyniki obliczeÅ„ dla º = 1 i różnych dyskretyzacji pokazano na rys. 9.19. Widoczny jest
efekt wypierania prądu w kierunku powietrza (zjawisko naskórkowości). Nierównomierna
gęstość prądu zwiększa straty mocy, a przez to i rezystancję. Istotnie, względna wartość
rezystancji R' = ReZ' > 1, co oznacza, że rezystancja dla prądu zmiennego jest większa niż
rezystancja tego samego przewodu dla prądu stałego.
a) b)
5 półokrÄ™gów × 10 sektorów, R' = 1,43 10 półokrÄ™gów × 20 sektorów, R' = 1,38
Rys. 9.19. Obraz pola w poÅ‚owie przekroju poprzecznego przewodu dla º = 1 dla zgrubnej (a) oraz
drobnej (b) dyskretyzacji; zaznaczono linie jednakowej względnej gęstości prądu (jego
modułu), położenie węzłów (w przypadku (a)  także wartości modułu względnej gęstości
prądu); kolor symbolizuje gęstość mocy czynnej (jasne kolory  duże wartości, ciemne
kolory  małe wartości); dla zgrubnej dyskretyzacji widoczne załamania i braki
w wypełnieniu kolorem
9.5.2. Schemat pięciopunktowy w przypadku osiowosymetrycznym
W przypadku zagadnień o symetrii osiowej stosuje się układ cylindryczny,
a laplasjan przyjmuje postać
"2Åš(r, z) 1 "Åš(r, z) "2Åš(r, z)
"2Åš = + + , (9.31)
"r2 r "r "z2
przy czym z założenia funkcja Ś nie zależy od współrzędnej kątowej Ć. Obszar
działania pola jest wtedy torusopodobnym kształtem o stałym przekroju, a analizę
pola można ograniczyć właśnie do tego przekroju. Najlepsze wyniki uzyskuje się
dla przekroju prostokątnego lub dającego się rozłożyć na sumę prostokątów. Na
przekrój nakładamy siatkę węzłów analogiczną do rozpatrywanej w punkcie 9.4.1,
przy czym rolę x przejmuje współrzędna r, a rolę y  współrzędna z (rys. 9.20):
rmax - rmin zmax - zmin
hr = , hz = , (9.32a)
I J
ri = rmin + ihr , z = zmin + jhz . (9.32b)
j
22
9.5. Inne typowe siatki węzłów
z
j
zmax
j = J
j (i, j)
zj
hz
hr
i
j = 0
zmin
i = I
i = 0 i
rmin ri rmax r
Rys. 9.20. Siatka dwuwymiarowa w układzie osiowosymetrycznym nałożona na przekrój
Wtedy dla węzła (i, j) laplasjan przyjmuje postać
Åši-1, j - 2Åši, j + Åši+1, j Åši+1, j -Åši-1, j Åši, j-1 - 2Åši, j + Åši, j+1
"2Åš H" + + . (9.33)
2 2
i, j
hr 2hrri hz
Jeżeli rmin = 0, to wyrażenie r 1"Ś/"r zawiera dzielenie przez 0. Jego granica przy
r 0 wynosi
1 "Åš "2Åš
lim = . (9.34)
r0
r "r "r2
Co więcej, uwzględnienie symetrii obrotowej prowadzi do wniosku, że wartości
funkcji pola na linii i =  1 są takie jak na linii i = 1. Tak więc w przypadku rmin = 0
laplasjan dla i = 0 przyjmuje postać
Åš1, j -Åš0, j Åš0, j -1 - 2Åš0, j + Åš0, j+1
"2Åš H" 4 + . (9.35)
2 2
0, j
hr hz
Przykład 9.7. Walec dielektryczny o promieniu a i wysokości 2h oraz przenikalności
elektrycznej µ naÅ‚adowano Å‚adunkiem o staÅ‚ej gÄ™stoÅ›ci objÄ™toÅ›ciowej Á = Á0 (rys. 9.21a).
Zakładając, że potencjał powierzchni walca wynosi 0, obliczyć rozkład potencjału
wewnÄ…trz walca. W obliczeniach przyjąć: Á0/µ = 100 V/cm2, a = 1 cm, h = 1 cm.
23
9. Metoda różnic skończonych (MRS)
a) b)
z
a
V = 0
h
h
2h
"2V = 0
Á0
"V
µ
V = 0
= 0
"n
0
a r
Rys. 9.21. Naładowany walec: a) szkic sytuacyjny, b) obszar poddany analizie
Z równań pola elektrostatycznego
" Å" D = Á, D = µE, E = -"V (a)
dostajemy równanie Laplace a dla potencjału elektrycznego:
Á
(b)
"2V = - .
µ
Z uwagi na geometrię przyjmujemy układ współrzędnych cylindrycznych z osią Oz
pokrywająca się z osią symetrii walca. Symetria obszaru i warunków brzegowych wskazuje
też, że potencjał V nie zależy od współrzędnej kątowej Ć, więc analizę można ograniczyć
do przekroju osiowego Ć = const. Ponadto dolna połowa walca jest odbiciem lustrzanym
górnej, więc ostatecznie analizowany obszar można ograniczyć do górnej połowy przekroju
osiowego walca (rys. 9.21b).
W rozpatrywanej sytuacji równanie Poissona (b) upraszcza się do postaci
"2V 1 "V "2V Á0
(c)
+ + = -
"r2 r "r "z2 µ
z następującymi warunkami brzegowymi:
- V = 0 na powierzchni bocznej walca (r = a) oraz na górnej podstawie (z = h),
- "V/"n = 0 na osi walca (r = 0) oraz w płaszczyz
nie symetrii z = 0.
Na analizowany obszar nakładamy siatkę węzłów rozmieszczonych w I + 1
kolumnach oraz J + 1 wierszach. Na podstawie (9.33) równanie (c) w wersji różnicowej dla
węzła (i, j) przyjmuje postać:
1 1 Á0 (d)
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚
2 2 2
1- Vi-1, j + 1+ Vi+1, j + Å› Vi, j-1 + Å› Vi, j+1 - (2 + 2Å› )Vi, j = - hr2,
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚
2i 2i µ
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
gdzie ś = hr/hz. Dla węzłów brzegowych należy je zmodyfikować:
- jeżeli i = I (czyli r = a) lub j = J (czyli z = h), wprowadzamy zerowy warunek Dirichleta:
(e)
Vi, j = 0,
- jeżeli i = 0 (czyli r = 0), należy skorzystać z równania wynikającego z (9.35)
Á0 2
2 2 2
(f)
4V1, j + Å› V0, j-1 + Å› V0, j+1 - (4 + 2Å› )V0, j = - hr ,
µ
- jeżeli j = 0 (czyli z = 0), zamiast Vi,j 1 należy wziąć Vi,1 (dotyczy równań (d) jak i (f)).
Wyniki obliczeń dla wybranych dyskretyzacji przedstawiono na rys. 9.22.
24
9.6. Podsumowanie
a) b)
Rys. 9.22. Wyniki obliczeÅ„: a) podziaÅ‚ 10×10 (121 wÄ™złów), b) podziaÅ‚ 21×21 (441 wÄ™złów)
9.6. Podsumowanie
9.6.1. Ogólna idea metody różnic skończonych
Metoda różnic skończonych rozwiązywania równań różniczkowych
zwyczajnych i cząstkowych polega na sprowadzeniu tego równania do postaci
różnicowej. W tym celu wszystkie pochodne zastępuje się odpowiednimi ilorazami
różnicowymi. Następnie na obszar działania pola nakłada się siatkę węzłów,
w których poszukuje się wartości funkcji pola. Występuje zatem tyle
niewiadomych, ile jest węzłów. Dla każdego węzła układa się równanie różnicowe
wynikające z oryginalnego równania. Dla węzłów brzegowych uwzględnia się
ponadto warunki brzegowe. W efekcie powstaje układ liniowych równań
algebraicznych, w którym niewiadomymi są wartości funkcji pola w węzłach
siatki. Po rozwiązaniu układu równań dostajemy przybliżone rozwiązanie. Należy
podkreślić, że macierz układu równań jest symetryczna, pasmowa, a przede
wszystkim rzadka, co pozwala na stosowanie specyficznych, szybkich metod
rozwiązywania takich układów równań, w szczególności metody iteracyjnej.
9.6.2. Zalety i wady metody
Najważniejsze zalety metody różnic skończonych to:
" teoretycznie nadaje się do rozwiązania dowolnego rodzaju równania
różniczkowego (zwyczajnego, cząstkowego, liniowego, nieliniowego),
" pozwala na uwzględnianie niejednorodności materiałowych,
" jest łatwa pojęciowo i nieskomplikowana do zaimplementowania na
komputerze,
" generuje układ równań o rzadkiej macierzy głównej, co umożliwia stosowanie
szybkich iteracyjnych metod rozwiązywania takiego układu.
25
9. Metoda różnic skończonych (MRS)
Podstawowymi wadami sÄ…:
" trudności w uwzględnieniu geometrii w przypadku nieregularnych kształtów,
" dla zbyt grubej siatki może prowadzić do dużych błędów (zob. przykład 9.4),
a zbyt drobna siatka prowadzi do ogromnej liczby węzłów i niewiadomych,
" istnieje co prawda możliwość lokalnego zagęszczenia siatki węzłów, ale
komplikuje to znacznie układanie samych równań,
" raczej nie nadaje się do rozwiązywania równań pola w obszarach
nieograniczonych, gdyż wymaga nałożenia na obszar siatki skończonej liczby
węzłów,
" jako metoda numeryczna, nie pozwala na uzyskanie zależności ogólnych.
Zakres zastosowań jest dość spory: zagadnienia pola elektromagnetycznego,
pola temperatury, mechaniki, hydromechaniki. Jej odmiana FDTD (Finite-
-Difference-Time-Domain  MRS w dziedzinie czasu) stosowana jest do m.in.
rozwiązywania równań falowych, a w szczególności rozchodzenia się fal
elektromagnetycznych, np. generowanych przez telefony komórkowe lub
propagujÄ…cych w falowodach (rys. 9.23).
a) b)
Rys. 9.23. Obrazy fali TE10 elektromagnetycznej (a) oraz impulsu elektromagnetycznego (b) rozchodzÄ…cego siÄ™
w falowodzie (za http://www.cemtach.com/reference/software/toyFDTD/ToyFDTD1)
26