Rozkład temperatury płytki (1)
Cienka płytka aluminiowa o wymiarach 40�40 cm jest izolowana z wyjątkiem
jej boków, które są utrzymywane w stałej temperaturze (patrz następny slajd).
Przepływ ciepła w kierunku prostopadłym do płytki nie występuje. Gęstość
aluminium wynosi � = 2.7 g/cm3, zaś ciepło właściwe cp = 0.2174 cal/(g�"�C).
Natomiast współczynnik przewodzenia
ciepła k = 0.49 cal/(s�"cm�"�C). Jaki
jest ustalony rozkład temperatury płytki?
q  gęstość strumienia ciepła, cal/(m2�"s)
"z  grubość płytki, m
"x"z  powierzchnia w kierunku y, m2
"y"z  powierzchnia w kierunku x, m2
q(x)"y"z"t + q(y)"x"z"t = q(x + "x)"y"z"t + q(y + "y)"x"z"t
Strumień ciepła Strumień ciepła
wpływającego wypływającego
Rozkład temperatury płytki (2)
q(x)"y"z"t + q(y)"x"z"t = q(x + "x)"y"z"t + q(y + "y)"x"z"t |: "x"y"z"t
q(x + "x) - q(x,t) q(y + "y) - q(y,t)
- =
Warunki brzegowe Dirichleta
"x "y
T (0, y) = 0.625y2
"x 0, "y 0
T (L, y) =100
" q " q
- =
T (x, 0) = 0.625x2
"x " y
T (x, L) =100
Korzystając z prawa Fouriera
"T "T
qx = -k qy = -k
"x " y
Dostajemy
"2T "2T
+ = 0
"x2 " y2
równanie
eliptyczne
Różnice skończone
1
"2T "2T
2 2
f (x) H" [ f (y + "y) - 2 f (y) + f (y - "y)]
+ = 0
("y)2
"x2 " y2
1
2 2
f (x) H" [ f (x + "x) - 2 f (x) + f (x - "x)]
("x)2
Przyjmijmy h = "x = "y
1 1
[T (x + h, y) - 2T (x, y) + T (x - h, y)] = [T (x, y + h) - 2T (x, y) + T (x, y - h)]
h2 h2
Oznaczenie: Ti, j = T (xi, yj)
1 1
[Ti+1, j - 2Ti, j + Ti-1, j ] = [Ti, j+1 - 2Ti, j + Ti, j-1]
h2 h2
4Ti, j - Ti-1, j - Ti+1, j - Ti, j-1 - Ti, j+1 = 0, dla 1d" i, j d" n
Aby można było zastosować metodę Gaussa-Seidela
przedstawiamy je w postaci
1
Ti, j = [Ti-1, j + Ti+1, j + Ti, j-1 + Ti, j+1] = 0
4
Różnice skończone. Algorytm (1)
input n, M
for i = 0 to n+1 do
Ti0 �! g (xi, 0); Ti,n+1 �! g (xi, L)
Obliczenie warunków
T0i �! g (0, yi); Tn+1,i �! g (L, yi)
brzegowych
end do
for i = 1 to n do
Ustalenie początkowych
for j = 1 to n do
przybliżeń Ti j w punktach
Ti j �! 0
wewnętrznych siatki
end do
end do
for k = 1 to M do
for j = 1 to n do
for i = 1 to n do
Metoda Gaussa-Seidela
Ti j �! (Ti-1, j + Ti+1, j + Ti, j-1 + Ti, j+1 )/4
end do
end do
end do
1
Ti, j = [Ti-1, j + Ti+1, j + Ti, j-1 + Ti, j+1] = 0, dla 1d" i, j d" n
Układ równań
4
Różnice skończone. Algorytm (2)
input n, M
for i = 0 to n+1 do
Ti0 �! g (xi, 0); Ti,n+1 �! g (xi, L)
Obliczenie warunków
T0i �! g (0, yi); Tn+1,i �! g (L, yi)
brzegowych
end do
Ustalenie początkowych przybliżeń Ti j
w punktach wewnętrznych siatki
for k = 1 to M do
Zamiast M można powtarzać
for j = 1 to n do
iteracje tak długo aż spełniony
for i = 1 to n do
będzie warunek
Ti j �! (Ti-1, j + Ti+1, j + Ti, j-1 + Ti, j+1 )/4
end do
Tinew - Tiold
j j
< � , dla 1d" i, j d" n
end do
Tinew
j
end do
Do przyspieszenia zbieżności można zastosować po każdej
iteracji nadrelaksację (1 <  < 2):
Tinew = Tinew + (1- )Tiold
j j j
Rozkład temperatury płytki (3)
12
100 100 100
10
100
80
8
60
40
6
20
4
0
40
40
30
2
20
20
10
0 2 4 6 8 10 12
0
T (0, y) = 0.625y2
T (L, y) =100
12 punktów na każdej osi
Warunki brzegowe
T (x, 0) = 0.625x2
50 iteracji
T (x, L) =100
100
T [C]
100
1
9
0
00
9
0
8
0
7
0
8
9
0
0
7
5
0
0
8
6
4
0
0
0
3
0
9
5
0
7
0
0
6
4
0
0
2
0
1
3
0
8
0
0
5
0
4
7
0
9
2
0
0
0
6
1
0
3
0
0
Metoda elementów skończonych (1)
"2u "2u
+ = f (x, y)
 równanie Poissona, eliptyczne
"x2 " y2
Etapy:
" dyskretyzacja
" budowa równań opisujących
przybliżone rozwiązanie dla
każdego elementu
" łączenie elementów
" uwzględnienie warunków
brzegowych
" rozwiązanie układu równań
Metoda elementów skończonych (2)
Funkcja interpolująca (płaszczyzna):
u(x, y) = c0 + c1x + c2 y
Z układu równań wyznaczamy
współczynniki c0, c1 i c2.
u1 = c0 + c1x1 + c2 y1
u2 = c0 + c2x2 + c2 y2
u3 = c0 + c1x3 + c2 y3
c0, c1 i c2 podstawiamy
Otrzymujemy
u(x, y) = u1n1(x, y) + u2n2(x, y) + u3n3(x, y)
gdzie n1, n2 i n3  funkcje kształtu (są określone tylko na trójkącie leżącym
na płaszczyz
nie xy)
Metoda elementów skończonych (3)
Funkcje kształtu:
1
n1(x, y) = [(x2 y3 - x3 y2 ) + (y2 - y3)x + (x3 - x2 )y]
2"
1
n2(x, y) = [(x3 y1 - x1y3) + (y3 - y1)x + (x1 - x4)y]
2"
1
n3(x, y) = [(x1y2 - x2 y1) + (y1 - y2)x + (x2 - x1)y]
2"
gdzie (" - pole trójkąta)
" = (x1 - x3)(y2 - y3) - (x2 - x3)(y1 - y3)
Metoda elementów skończonych (4)
Funkcja interpolująca Nk jest określona
tylko na trójkątach otaczających węzeł k.
Poza tym obszarem jest równa 0.
Nk (x, y) = n1k (x, y) + n2k (x, y) +L�+ nm k (x, y)
k
gdzie mk  liczba elementów otaczających węzeł k.
Rozwiązane r.r.cz. ma postać
u(x, y) = u1N1(x, y) + u2N2(x, y) +K�+ uK NK (x, y)
gdzie K  liczba węzłów. Trzeba znalez
ć współczynniki u1, u2, ..., uK
które są rozwiązaniem układu równań
K K
�ł łł
k
0 =
�ł
k
"u "Nk "Ni +""N "Ni + f (x, y)śłdxdy dla i =1,2,L�, K
+"+"
&!
"x "x "y "y
�ł k =1 k =1 �ł
Metoda objętości kontrolnej (1)
Bilans ciepła w stanie ustalonym
Szukamy równania dla punktu (2, 4)
0 = Q1 - Q2 + Q3 + Q4 - Q5
h h h h h
0 = q1 "z - q2 "z + q3 "z + q4 "z - q5 "z
2 2 2 4 4
Q  strumień ciepła, cal/s
0 = 2q1 - 2q2 + q3 + q4 - q5
q  gęstość strumienia ciepła, cal/(cm2s)
Metoda objętości kontrolnej (2)
Prawo Fouriera wyrażone za
pomocą różnic skończonych
T24 - T14
q1 = -kA
h
T34 - T24
q2 = -kB
h / 2
T24 - T23
q3 = -kA
h
T24 - T23
q4 = -kB
h
0 = 2q1 - 2q2 + q3 + q4 - q5
Gęstość strumienia ciepła
spowodowana konwekcją
q5 = -kc (Ta - T24)
(3kA + 5kB + kc )T24 - 2kAT14 - 4kBT34 - (kA + kB )T23 = kcTa
Równanie dla punktu (2, 4)
Metody rozwiązywania r.r.cz.
metoda
Różnic Elementów Objętości
skończonych skończonych kontrolnej
np. np.
Multiphysics Fluent
Projekt najtańszego zbiornika (1)
Zaprojektuj najtańszy zbiornik w kształcie walca o pojemności V = 0.8 m3 do
przewozu toksycznych odpadów. Ponieważ zbiornik ma być zamontowany na
samochodzie jego średnica wewnętrzna nie może być większa niż Dmax = 1 m,
a jego długość nie większa niż Lmax = 2 m. Grubość ścian zbiornika specjalne
przepisy określają na w = 0.03 m. Gęstość materiału z którego ma być
wykonany zbiornik wynosi � = 8000 kg/m3, zaś jego cena cm = 12.6 zł/kg.
Koszt spawania wynosi cs = 56 zł/m.
C  koszt całkowity
ls  całkowita długość spawu
m  masa materiału zbiornika
2 2
�ł�ł D łł
D
�ł �ł
Vcylindra = LĄ + w�ł -
�ł�ł �ł �ł �ł śł
�ł�ł 2 łł �ł 2 łł śł
�ł �ł
2
D
D
�ł
lwew = 2Ą
Vplyty = Ą + w�ł w
�ł �ł
2
2
�ł łł
m = �(Vcylindra + 2Vplyty)
D
�ł
lzew = 2Ą + w�ł
�ł �ł
C = cmm + csls
ls = 2(lzew + lwew)
2
�ł łł
Projekt najtańszego zbiornika (2)
2
�ł łł
D
C(D, L) = cm�Ą w L(D + w) - 2�ł + w�ł + 4csĄ (D + w)  funkcja celu
�ł �ł �ł śł
2
�ł łł
�ł śł
�ł �ł
minC
D, L  zmienne decyzyjne
D,L
2
D
�ł �ł
 ograniczenie równościowe, Vzbiornika = 0.8 m3
Ą L = 0.8
�ł �ł
2
�ł łł
D d"1
 ograniczenia nierównościowe, Dmax = 1 m
L d" 2
Lmax = 2 m
Optymalizacja
2
1.5
1
0.5
0
1
1
0.5
1
0
0
-0.5
-1
0.5 -1
y
x
0
-0.5
-1
1
1
0.5
0
0
-0.5
-1
-1
y
x
f(x,y)
f(x,y)
Metoda złotego podziału odcinka (1)
Zakładamy, że na odcinku [a, b] funkcja f (x) posiada tylko jedno min  tzn.
jest unimodalna.
l1 + l2 l1
l0 = l1 + l2
=
l1 l2
l1 l2
=
Oznaczamy r = l2/l1
l0 l1
1
1+ r =
r
r2 + r -1 = 0
-1+ 1- 4(-1)
5 -1
r = = = 0.61803K�
2 2
Metoda złotego podziału odcinka (2)
input a, b, �, �
2
x �! a + r(b - a); y �! a + r (b - a)
u �! f (x); v �! f (y);
while |x - y| > � and |u - v| > � do
if u > v then
b �! x; x �! y; u �! v
2
y �! a + r (b - a) ; v �! f (y)
else
a �! y; y �! x; v �! u
2
x �! a + r (b - a) ; u �! f (x)
end if
end do
end do
5 -1
r =
2
Metoda złotego podziału odcinka. Przykład
w przedziale [0, 3]
Minimum f (x) = (x2 - 4)2 /8 -1
f(u) <= f(v) [v, b] f(u) > f(v) [a, u]
3 3
2 2
u = 1.8541 u = 2.2918
1 1
v = 1.1459 v = 1.8541
a v u b a v u b
0 0
-1 -1
0 1 2 3 0 1 2 3
f(u) <= f(v) [v, b] f(u) <= f(v) [v, b]
3 3
2 2
u = 1.8541 u = 2.0213
1 1
v = 1.5836 v = 1.8541
a v u b a v u b
0 0
-1 -1
0 1 2 3 0 1 2 3
xopt = 2.0000, f (xopt) = -1.0000
Funkcja fminbnd
x = fminbnd (fun, x1, x2)  znajduje minimum funkcji jednej zmiennej fun
w przedziale [x1, x2] metodą złotego podziału
x = fminbnd (fun, x1, x2, options)  uwzględnia opcje optymalizacyjne podane
w strukturze options
[x, fval] = fminbnd (...);
[x, fval, exitflag] = fminbnd (...);
min f (x)
fun  uchwyt funkcji albo m-plik
x
options  struktura tworzona przez funkcję optimset
x1 d" x d" x2
fval  wartość funkcji f w znalezionym punkcie x
exitflag  stan wyjścia z funkcji fmincon
f = @ (x) (x^2 - 4)^2/8 - 1;
[x, y] = fminbnd (f, 0 ,3)
x = 2.0000
Niektóre parametry struktury options:
y = -1.0000
TolX  dokładność końcowa rozwiązania
MaxFunEvals  maksymalna ilość obliczeń (wywołań) funkcji
MaxIter  maksymalna ilość iteracji
Metoda Neldera-Meada (1)
Sympleks w R2
Sympleks w R3
Algorytm Neldera-Meada służy do wyznaczania minimum funkcji f : RnR
bez badania tej funkcji na prostych i bez liczenia jej pochodnych.
Na początku trzeba przyjąć wartości parametrów ą > 0, � " (0, 1), ł > 0.
Można przyjąć ą = 1, � = 1/2 , ł = 1.
W każdym kroku dany jest układ n+1 punktów z Rn: {x(0), x(1), ... , x(n)}.
Metoda Neldera-Meada (2)
tworzymy
3 punkty
input x(0), x(1), ... , x(n), �, M
porządkujemy punkty x(i) tak aby f (x(0)) e" f (x(1)) e" ... e" f (x(n))
for k = 1 to M do
n
1
(i)
u �!
"x , v �! (1+ ą)u -ąx(0), w �! (1+ ł )v - łu
n
i=1
środek ściany
if f (v) < f (x(n)) then
sympleksu
1
if f (w) < f (x(n)) then x(0) �! w else x(0) �! v end
położonej
end
naprzeciw
2
if f (v) e" f (x(n)) '" f (v) d" f (x(1)) then x(0) �! v end
najgorszego
if f (v) > f (x(1)) then
punku x(0)
if f (v) d" f (x(0)) then x(0) �! v; w �! � x(0) + (1-� )u end
if f (w) d" f (x(0)) then
x(0) �! w
3
else
ściągnięcie
x(i) �! (x(i) +x(n))/2 dla i = 0, 1, ... , n-1
simpleksu
end
end
porządkujemy punkty x(i) tak aby f (x(0)) e" f (x(1)) e" ... e" f (x(n))
if | f (x(0)) - f (x(n))| / ( | f (x(0)) | + | f (x(n)) | ) < � then koniec end
end do
Metoda Neldera-Meada. Przykład (1)
Rys. 1
3
f (x, y) = x *(x - 4 - y) + y(y -1)
2.5
f (x(0)) e" f (x(1)) e" f (x(2))
2
Przypadek 1: x(0) �! w
w
1.5
x(1)
v
Rys. 2
1
3
u
0.5
2.5
0
0 1 2 3 4
2
x(2)
x(0)
x(2)
1.5
x(0) 1
u
0.5
v
w
Przypadek 1: x(0) �! v
0
0 1 2 3 4
x(1)
Metoda Neldera-Meada. Przykład (2)
Rys. 3
w
3
Przypadek 1: x(0) �! w
2.5
2
v v
x(2)
1.5
u
1
Rys. 4
x(2)
3
0.5
x(1)
2.5
u
0
0 1 2 3 4
x(0)
2
x(1)
1.5
1
0.5
Przypadek 3: ściągnięcie simpleksu
x(0)
0
0 1 2 3 4
Metoda Neldera-Meada. Przykład (3)
Rys. 5
x(0)
3
Przypadek 3: ściągnięcie simpleksu
2.5
x(2)
u
2
x(1)
1.5
v
Rys. 6
1
3
v
0.5
2.5
x(1)
0
0 1 2 3 4
x(2)
2 u
1.5
x(0)
1
0.5
Przypadek 3: ściągnięcie simpleksu
0
0 1 2 3 4
Metoda Neldera-Meada. Przykład (4)
Rys. 7
3
2.5
2
xopt = [3, 2]
f (xopt ) = 7
1.5
1
0.5
0
0 0.5 1 1.5 2 2.5 3 3.5 4
Funkcja fminsearch
x = fminsearch (fun, x0)  startując z punktu x0 znajduje minimum funkcji wielu
zmiennych fun metodą sympleksów Lagariasa
x = fminsearch (fun, x0, options)  uwzględnia opcje optymalizacyjne podane
w strukturze options
[x, fval] = fminsearch (...);
[x, fval, exitflag] = fminsearch (...);
fun  uchwyt funkcji albo m-plik
min f (x)
x
options  struktura tworzona przez funkcję optimset
fval  wartość funkcji f w znalezionym punkcie x
exitflag  stan wyjścia z funkcji fmincon
f = @ (x) x(1)*(x(1)-4-x(2))+x(2)*(x(2)-1);
[x, y] = fminsearch (f, [0 0])
Niektóre parametry struktury options:
x = 3.0000 2.0000
TolX  dokładność końcowa rozwiązania
y = -7.0000
TolFun  dokładność końcowa funkcji
MaxFunEvals  maksymalna ilość obliczeń (wywołań) funkcji
MaxIter  maksymalna ilość iteracji
Funkcja fmincon (1)
x = fmincon (fun, x0, A, b)  startując z punktu x0 znajduje minimum funkcji fun
przy ograniczeniach nierównościowych Ax d" b
x = fmincon (fun, x0, A, b, Aeq, beq)  znajduje minimum funkcji fun przy
ograniczeniach nierównościowych Ax d" b
i równościowych Aeq"x = beq. Jeśli
ograniczenia nierównościowe nie
występują podstaw A = [ ] i b = [ ].
fun  funkcja o argumencie wektorowym, przyjmująca
min f (x)
x
wartości skalarne, której minimum poszukujemy
(uchwyt funkcji albo m-plik). c(x) d" 0
ceq (x) = 0
x0  punkt startowy
Ax d" b
A, Aeq  macierze
Aeq x = beq
B, beq  wektory
lb d" x d" ub
x  rozwiązanie
Funkcja fmincon (2)
x = fmincon (fun, x0, A, b, Aeq, beq, lb, ub)  nakłada ograniczenia na zmienne
decyzyjne tak, że rozwiązanie
spełnia warunek
lb d" x d" ub
x = fmincon (fun, x0, A, b, Aeq, beq , lb, ub, nonlcon)  znajduje minimum przy
nieliniowych ograniczeniach
nierównościowych c(x) d" b oraz
równościowych ceq(x) = 0. Jeśli
ograniczenia lb i ub nie występują
podstaw lb = [ ] i ub = [ ].
lb  wektor, jeżeli przyjmiemy lb(i) =  Inf to zmienna deczyjna x(i) nie jest
ograniczona od dołu
ub  wektor, jeżeli przyjmiemy bb(i) = Inf to zmienna deczyjna x(i) nie jest
ograniczona od góry
nonlcon  funkcja definiująca nieliniowe oraniczenia c(x) d" b oraz ceq(x) = 0.
Funkcja fmincon (3)
x = fmincon (fun, x0, A, b, Aeq, beq, lb, ub , nonlcon, options)  uwzględnia opcje
optymalizacyjne podane w strukturze
options. Jeśli funkcja definiująca
ograniczenia nieliniowe nie
występuje podstaw nonlcon = [ ].
[x, fval] = fmincon (...);
[x, fval, exitflag] = fmincon (...);
[x, fval, exitflag, output] = fmincon (...);
options  struktura definiująca opcje optymalizacyjne tworzona przy pomocy
funkcji optimset
x  rozwiązanie
fval  wartość funkcji f w znalezionym punkcie x
exitflag  stan wyjścia z funkcji fmincon
output  struktura zawierająca informacje o obliczeniach
Funkcja optimset
options = optimset ('name1', value1, 'name2', value2, ...)  tworzy strukturę
zawierającą parametry
całkowania równań
różniczkowych.
'name1'  nazwa parametru
value1  wartość parametru
Parametr
TolFun  dokładność końcowa wyznaczenia funkcji celu (1e-3)
TolCon  dokładność końcowa przekroczenia ograniczeń
TolX  dokładność końcowa rozwiązania
MaxFunEvals  maksymalna ilość obliczeń (wywołań) funkcji
MaxIter  maksymalna ilość iteracji
Projekt najtańszego zbiornika (3)
function koszt = ZbiorOpt (x, w, ro, cm, cs)
D= x(1);
L = x(2);
function [c, ceq] = ogran (x)
R = D/2;
c = -0.1; % ogr. nieistotne
R2 = R^2;
ceq = pi*x(1)^2*x(2)/4 - 0.8;
Rw = R + w;
Rw2 = Rw^2;
m = ro*pi*(L*(Rw2-R2) + 2*Rw2*w);
Vzbiornika = 0.8 m3
ls=4*pi*(D + w);
koszt=cm*m + cs*ls;
f = @ (x) ZbiorOpt (x, 3, 8000, 12.6, 56);
x = fmincon (f, [1.0 2.0], [ ], [ ], [ ], [ ], [0; 0], [1; 2], @ogran);
L0
D0
x =
D d"1
D e" 0
1.0000 1.0186
L d" 2
L e" 0
D L