Rozkład temperatury płytki (1)
Cienka płytka aluminiowa o wymiarach 4040 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łaszczyznie 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
Wyszukiwarka
Podobne podstrony:
Wyklad studport 8wykład 13 i 14 stacjonarneBOwL wyklad Beamer 14Electronic Commerce wyklad ie? 14 prawoWyklad studport 9SS wyklad nr 14 pptWyklad 05 14 15 GWWyklad studport 5PSM Wyklad 13 14wyklad 10 14 12 20100214 13 10 2009, wykład nr 14 , Układ pokarmowy, cześć II Paul EszWyklad studport 11BOwL wyklad Beamer 14Wyklad studport 7Wyklad studport 12Wyklad 01 14 15 GWBOwL wyklad Beamer 14Wyklad 02 14 15 GWwięcej podobnych podstron