METODA RÓŻNIC SKOC
CZONYCH 2
(MRS)
Finite Difference Methods (FDM)
Różnice centralne
1
" 1
E"
0
"y 2hy
" 1
E"
-1 1
0
"x 2hx
-1
"2 1
E"
2 1 -2 1
"x2 hx
1
"2 1
E"
E"
2 -2
-2
"y2 h2
"y2 hy
1
Operator Laplace a (hx=hy=h)
1
"2 "2 1
+ E"
1 -4 1
"x2 "y2 h2
1
Klasyfikacja równań różniczkowych cząstkowych II rzędu
n n
"2u "u
k = 1,2,K, n
"a (xk )"xi"x + "b (xk )"xi + c(xk )u = f (xk ),
ij i
i, j=1 i=1
j
" eliptyczne: aij dodatnio określona
" hiperboliczne: jedna z wartości własnych aij jest ujemna, pozostałe
sÄ… dodatnie
" paraboliczne: det(aij)=0
Równanie przewodnictwa ciepła
(równanie paraboliczne)
"
ëÅ‚" - cÁ öÅ‚T + qv = 0
2
ìÅ‚ ÷Å‚
"t
íÅ‚ Å‚Å‚
T  temperatura
  współczynnik przewodzenia ciepła
c  ciepło właściwe
c  ciepło właściwe
Á  gÄ™stość materiaÅ‚u
qv  wewnętrzne z
ródła ciepła
Warunki brzegowe na brzegu obszaru:
1. rozkład temperatury
2. gęstość strumienia ciepła
3. określony sposób wymiany ciepła z otoczeniem (np. konwekcja)
Równanie falowe
(równanie hiperboliczne)
&&
u - c2"2u = 0
Warunki poczÄ…tkowe:
u(x,0) = u0(x)
&
&
u(x,0) = v (x)
u(x,0) = v0(x)
Warunki brzegowe:
u(x,t) = u(x,t)
"&!
Równanie membrany
(równanie eliptyczne)
p
"2w = -
S
p  ciśnienie
S  siła naciągająca
Warunki brzegowe:
w = 0
"&!
Równanie płyty (równanie bi-harmoniczne)
1
p
2 2
-8
"4w = -
D
1 -8 20 -8 1
p  obciążenie poprzeczne
D  sztywność płyty
2 2
-8
-8
1
Warunki brzegowe:
utwierdzenie swobodne podparcie brzeg swobodny
w = 0
w = 0
M = 0
n
"w
"M
"2w = 0 nt
Õ = = 0
Qn - = 0
"n
"s
Przykład: przepływ ciepła w pręcie
x = 0 x = a
g1(t) g2(t)
"T "2T
"T "2T
= c ; 0 < x < a, 0 d" t d" Ä
"t " x2
T(x,0 ) = f(x), 0 < x < a
T( 0,t) = g1(t)
Å„Å‚
òÅ‚T(a,t) = g2(t) , 0 d" t d" Ä
ół
Dyskretyzacja obszaru rozwiÄ…zania
przestrzeń: h = " czas: k = "
"x, "t
" "
" "
10
t
9
8
7
6
5
5
4
3
2
Czas
1
(j index)
0
0 1 2 3 4 5 6 7 8 9 10
przestrzeń
x
(i index)
Warunki brzegowe
10
9
8
7
6
T(a, t)
T(a, t)
T(0, t)
T(0, t)
5
= g2(t)
4
= g1(t)
3
2
1
0
0 1 2 3 4 5 6 7 8 9 10
Warunki poczÄ…tkowe : T(x,0) = f(x)
Schemat różnicowy
t
tj+1
T(x,t) (i,j+1)
tj
t
(i-1,j) (i,j) (i+1,j)
x
xi-1 xi xi+1
x
1
Tt = (Ti, j+1 - Ti, j )
czas: explicit
k
przestrzeń: centralny
c
cTxx = (Ti-1, j - 2Ti, j + Ti+1, j )
h2
h = " x = a / n, xi = ih
Å„Å‚
òÅ‚k = " t = T / m, t = jk
j
ół
1 c
Tt = cTxx Ò! (Ti, j+1 - Ti, j ) = (Ti-1, j - 2Ti, j + Ti+1, j )
k h2
ck
ck
Ti, j+1 = Ti, j + (Ti-1, j - 2Ti, j + Ti+1, j )
h2
= rTi-1, j + (1 - 2r)Ti, j + rTi+1, j
ck c" t
r = =
h2 " x2
Warunek stabilności (Courant):
0 < r d" 0.5
Explicit Euler Method: Stability
r = 1
Unstable !!
Schemat różnicowy
t
tj+1 (i-1,j+1) (i,j+1) (i+1,j+1)
T(x,t)
tj
t
(i,j)
x
xi-1 xi xi+1
x
1
czas: implicite ut = ( ui , j+1 - ui , j )
= -
= -
= -
+
+
+
k
c
przestrzeń: centralny
cuxx = ( ui-1, j+1 - 2ui , j+1 + ui+1, j+1 )
= - +
= - +
= - +
- + + + +
- + + + +
- + + + +
h2
" Implicit Euler method
1 c
(Ti, j+1 - Ti, j ) = (Ti-1, j+1 - 2Ti, j+1 + Ti+1, j+1)
k h2
- rTi-1, j+1 + (1 + 2r)Ti, j+1 - rTi+1, j+1 = Ti, j
" Macierz trójdiagonalna
T T + rT
T1, j+1 T1, j + rT0, j+1
1 + 2r - r Å„Å‚ üÅ‚ Å„Å‚ üÅ‚
1 + 2r - r Å„Å‚ üÅ‚ Å„Å‚ üÅ‚
îÅ‚ Å‚Å‚
îÅ‚ Å‚Å‚
+ +
ôÅ‚ ôÅ‚
ïÅ‚ śł
T2, ôÅ‚ ôÅ‚ T2, j
- r 1 + 2r - r
ôÅ‚ ôÅ‚ ôÅ‚ ôÅ‚
ïÅ‚ śłôÅ‚ j+1 ôÅ‚ ôÅ‚
ôÅ‚
T3, T3, j
=
ïÅ‚ - r 1 + 2r - r śłòÅ‚ j+1 żł òÅ‚
żł
ïÅ‚ śłôÅ‚ M ôÅ‚ ôÅ‚
ôÅ‚
O O - r
ïÅ‚ śłôÅ‚
ôÅ‚ ôÅ‚ ôÅ‚
ïÅ‚ śł Tn-1, j+1þÅ‚ ółTn-1, j + rTn, j+1þÅ‚
- r 1 + 2rûłół ôÅ‚ ôÅ‚ ôÅ‚
ôÅ‚
ðÅ‚
" Bezwarunkowo stabilna !
Crank-Nicolson method
Crank-Nicolson method
Implicit Euler method : pierwszego rzędu
Crank-Nicolson : drugiego rzędu
10
9
8
7
u(a, t)
u(a, t)
u(0, t)
u(0, t)
6
= g2(t)
5
= g1(t)
4
3
2
1
0
0 1 2 3 4 5 6 7 8 9 10
Warunki poczÄ…tkowe : u(x,0) = f(x)
" Średnia między dwoma poziomami czasu
1 c
(Ti, j+1 - Ti, j ) = (Ti-1, j - 2Ti, j + Ti+1, j )
k 2h2
c
+ (Ti-1, j+1 - 2Ti, j+1 + Ti+1, j+1)
2h2
" Macierz trójdiagonalna
" Macierz trójdiagonalna
r r r r
- Ti-1, j+1 + (1 + r)Ti, j+1 - Ti+1, j+1 = Ti-1, j + (1 - r)Ti, j + Ti+1, j
2 2 2 2
" Bezwarunkowo stabilna
Operatory dla równania przepływu ciepła
 "t
Ä… =
cÁ "x2
1
-Ä… 1+2Ä… -Ä… 1+Ä… -Ä…/2
-Ä…/2
1+2Ä…
-Ä… -Ä…
-Ä…/2 -1+Ä… -Ä…/2
-1
standardowy standardowy implicit
explicit implicit Crank-Nicholson
Równanie falowe
&& 2 2
u - c2u = 0
ui, j+1 - 2ui, j + ui, j-1 ui+1, j - 2ui, j + ui-1, j
- c2
"t2 "x2
"t2
"t2
( ) ( )
ui, j+1 = (2 - 2Ä…)ui, j + Ä…(ui+1, j + ui-1, j)- ui, j-1, Ä… = c2
"x2
-1
Warunek stabilności: ą<=1
2-2Ä…
Ä… Ä…
-1
Brzeg krzywoliniowy: operator Laplace a
2
"u "2u ("x)
uO H" uM + "x +
N
"x "x2 2
M
M
2
x
U
"u "2u (µ"x)
"y
uV H" uM - µ"x +
x
V M
"x "x2 2 W O
M
M
sumujac
S
"2u µuO + uV - (1 + µ)uM
= 2
= 2
2
2
µ"x
µ"x
"x2
"x2
µ(1 + µ)("x)
µ(1 + µ)("x)
M
Podobnie rozpisujemy operator w kierunku y
Operator Laplace a
ëÅ‚ öÅ‚
µuO + uV - (1 + µ)uM uS + uU - (1 + )uM ÷Å‚
ìÅ‚
"u = 2ìÅ‚ +
2 2
M ÷Å‚
µ(1 + µ)("x) (1 + )("y)
íÅ‚ Å‚Å‚
Nieregularna siatka węzłów
Krzywoliniowa MRS (CFD  Curvilinear Finite Difference)
t
y
(-1,1) (1,1)
s
x
x
(-1,-1) (1,-1)
I J
ux = ussx + uttx
J
x(s, t) = X ÕiI (s)È (t) a" F(X, s, t)
"" ij j
i=1 j=1
uy = ussy + utty
I J
J 2 2
y(s, t) =
uxx = usssx + 2ustsxtx + utttx + ussxx + uttxx
""Y ÕiI (s)È (t) a" F(Y, s, t)
ij j
i=1 j=1
K
I J
J
u(s, t) =
-1
""U ÕiI (s)È (t) a" F(U, s, t)
ij j
sx = J xt
i=1 j=1
J = xs yt - xt ys
Dowolna nieregularna siatka węzłów
y
x
Rozwijamy funkcje w szereg Taylora i tworzymy układ równań
Rozwijamy funkcje w szereg Taylora i tworzymy układ równań
qD = v
2 2
îÅ‚ Å‚Å‚
h1 k1 h1 / 2 k1 / 2 h1k1
ïÅ‚h k2 h2 / 2 k2 / 2 h2k2 śł
2 2
2
ïÅ‚ śł
q[m×5] = , DT = [ux uy uxx uyy uxy]
ïÅ‚K K K K K śł
ïÅ‚ śł
2 2
km hm / 2 km / 2 hmkm ûÅ‚
ðÅ‚hm
(k
vT = [u1(k ) K um )]
hi = xi - xk , ki = yi - yk
2
m
~
Åš = (ui(k ) - ui(k )) wi2
" (k )
i=1
" Åš
= 0 Ò! D = Qv
"D
p
wi,(k ) = (hi2 + ki2)
Operator Laplace a
1 1
1 4
1 -4 1 4 4
6 x -20
1 1
1 4