pobieranie: Modelowanie i Symulacja, Modelowanie obiektów o stałych rozłożonych
Pobierz dokument pdfPrzeglądaj wersję html pliku:
Laboratorium Modelowania i Symulacji Zespół Automatyki ZTMAiPC Wydział Elektryczny
LABORATORIUM MODELOWANIA I SYMULACJI
Ćwiczenie 8 Modelowanie obiektów o stałych rozłożonych
1. Cel ćwiczenia
Celem ćwiczenia jest zapoznanie się ze sposobami modelowania układów o stałych rozłożonych. W instrukcji zawarto podstawowe informacje na temat metody różnic skończonych wykorzystywanej do rozwiązywania jednowymiarowego zagadnienia przepływu ciepła. Dodatkowo przedstawiono metodę nadrelaksacji na przykładzie równania Laplace'a.
2. Obiekty o stałych rozłożonych 2.1 Jednowymiarowy, nieustalony przepływ ciepła w ciele stałym. Warunki początkowe i brzegowe
Rys. 1 przedstawia model obiektu o stałych rozłożonych. Modelowanym zjawiskiem jest ruch ciepła w pręcie o długości L = 1, przez który przepływa prąd elektryczny. Przepływ prądu powoduje generowanie ciepła wewnątrz pręta. Oba końce pręta są zanurzone w temperaturze T = 0. W tym przypadku pomimo tego, że obiekt modelowany ma trzy wymiary geometryczne zagadnienie określa się mianem jednowymiarowego. Jest to związane z faktem, że przepływ ciepła odbywa się tylko w jednym kierunku (od wnętrza walca na zewnątrz). Warunki brzegowe T(0,t) oraz T(1,T) określają wartość temperatury na końcach pręta (brzegach obszaru). Oprócz warunków brzegowych dla rozwiązania zagadnienia opisywanego cząstkowym równaniem różniczkowym niezbędna jest znajomość warunku początkowego. Podobnie jak w zagadnieniach opisywanych równaniami zwyczajnymi, warunki początkowe określają, jakie wartości musi przyjmować rozwiązanie w chwili t = 0.
Warunki brzegowe
T 0, t = 0
T 1, t = 0
Warunek początkowy
T x ,0 = sin x x 1 − x
Rysunek 1. Geometria, warunki brzegowe i warunek początkowy
Do opisu matematycznego opisywanego zjawiska stosuje się jednowymiarowe równanie nieustalonego przepływu ciepła:
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
1
Laboratorium Modelowania i Symulacji ∂ T x , t 1 ∂2 T x , t − = f x ,t, ∂t k ∂ x2
(1)
podawane wraz z warunkami jednoznaczności rozwiązania (brzegowymi i początkowym). W równaniu (1): T(x,t) – wartość temperatury w punkcie x, w chwili czasu równej t; k – współczynnik dyfuzyjności cieplnej (zależny od rodzaju materiału); f(x,t) – funkcja wydajności cieplnej objętościowych źródeł mocy. Można wykazać, że przy podanych warunkach brzegowych równanie (1) posiada następujące rozwiązanie:
T x , t =e t sin x x 1−x 2.2 Dwuwymiarowy, ustalony przepływ ciepła w ciele stałym
(2)
Na Rys. 2 przedstawiono dwuwymiarowe zagadnienie ustalonego przepływu ciepła wraz z warunkami brzegowymi.
T x ,2
Warunek brzegowy
T x ,2 = sin x sinh 2
x
Brzeg obszaru
T , y = 0
T 0, y =0
Element siatki przestrzennej
T x ,0 = 0
Rysunek 2. Dwuwymiarowy, ustalony przepływ ciepła
Zagadnienie opisuje równanie Laplace'a:
∂ T x , y ∂ T x , y =0. ∂ x2 ∂ y2
(3)
Problem brzegowy, ustalony nie wymaga podawania warunku początkowego. Jeżeli w równaniu (3) prawa strona zostanie zastąpiona niezerową funkcją, równanie nazywamy równaniem Poissona. Rozwiązanie równania (3) jest uzależnione od warunków brzegowych panujących na granicach (brzegach) rozpatrywanego obszaru. Przy uwzględnieniu warunków brzegowych jak na rys. 2 ma postać:
T x , y=sin xsinh y
(4)
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
2
Laboratorium Modelowania i Symulacji 3. Podstawy metody różnic skończonych 3.1 Wprowadzenie
Do rozwiązywania zagadnień brzegowych jedno i wielowymiarowych często stosuje się metodę różnic skończonych opartą na aproksymacji pochodnej funkcji ilorazem różnicowym „w przód”. Dla uproszczenia rozpatrzmy pochodną funkcji jednej zmiennej:
f ' x 0 ≃
f x 0 h− f x 0 h
(5)
gdzie: h oznacza niewielki rozmiar kroku. Aby wprowadzić pojęcie rzędu aproksymacji rozwińmy f x 0 h w szereg Taylora wokół punktu x0:
f x 0h− f x 0 h h2 = f ' x 0 f ' ' x 0 f ' ' ' x 0 RN h 2 6
(6)
gdzie RN oznacza składniki rozwinięcia wyższych rzędów. Przyjmując, że h jest bardzo małe, błąd popełniany przez stosowanie równania (5) jest w przybliżeniu równy
h f ' ' x 0 . Mówimy, że błąd aproksymacji w równaniu (5) jest rzędu h i oznaczamy go 2
O(h). Tego rodzaju aproksymacja nie jest satysfakcjonująca. Dla przykładu, dokładność 0.001, będzie wymagała bardzo małego kroku h. Z tego powodu wymagana jest aproksymacja rzędu h2 lub wyższego. Poniżej przedstawiono aproksymację rzędu h2, opartą na centralnym ilorazie różnicowym:
f ' x 0 ≃
f x 0h− f x 0 −h . 2h
(7)
W późniejszych rozważaniach niezbędna będzie również znajomość przybliżenia f''(x0):
f ' ' x 0 ≃
f x 0 h−2 f x 0 f x 0 −h h2
.
(8)
3.2 Aproksymacja rozwiązania jednowymiarowego równania przewodzenia ciepła
Poszukujemy rozwiązania równania (1) dla zadanych warunków brzegowych. Wprowadźmy oznaczenia: h=x i1− x i oraz t=t j1−t j jako odległość pomiędzy wartościami rozwiązania w punktach (xi, yj) (współrzędne siatki czasowo-przestrzennej). Oznaczając: T i , j =T x i , t j można równanie (1) zapisać jako:
T i , j1−T i , j 1 T i1, j −2T i , j T i− j , j = f i, j. k h2 t
Przekształcając (9) ze względu na Ti,j otrzymuje się:
(9)
T i , j1=1−2 T i , j T i1, j T i−1, j tf i , j .
gdzie: =
(10)
t , i = 1, 2, ... n kh 2
Wygodnie jest zapisać równanie (10) w postaci macierzowej. Przyjmując, że Tj oznacza wektor (T1j, T2j, ..., Tnj): Tj+1 = ATj + fj, gdzie:
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
3
Laboratorium Modelowania i Symulacji 1−2 1−2 A= 0 ⋱ ⋮ ⋱ 0 ⋯ 3.3 Metoda Cranka – Nicolsona.
0 ⋱ ⋱ 0
⋯ 0 ⋯ 0 ⋱ ⋮ ⋱ 1−2
(11)
W celu poprawy stabilności numerycznej oraz zwiększenia dokładności obliczeń, często stosuje się metodę Cranka – Nicolsona (C-N). W metodzie tej wprowadza się macierz B:
2 −1 B= 0 ⋮ 0
Oznaczając przez I – macierz jednostkową, można wprowadzić następujące macierze:
−1 2 ⋱ ⋱ ⋯
0 −1 ⋱ ⋱ 0
⋯ ⋯ ⋱ ⋱ −1
0 0 ⋮ −1 2
(12)
C =2 I B , D=2 I − B
zagadnienie brzegowe rozwiązuje się zgodnie z równaniem macierzowym:
(13) (14)
CT j1 =DT j 2 tf
j
3.4 Aproksymacja rozwiązania dwuwymiarowego zagadnienia ustalonego przepływu ciepła – metoda nadrelaksacji
Zagadnienie przedstawione w punkcie 2.2, zostanie rozwiązane metodą nadrelaksacji. W metodzie tej wyróżnić można następujące etapy: 1. Podział obszaru zagadnienia (utworzenie siatki) o odległościach pomiędzy węzłami odpowiednio: Δx w poziomie oraz Δy w pionie. 2. Przypisanie warunków brzegowych węzłom brzegowym, a węzłom wewnętrznym dowolnych wartości początkowych. 3. Wyznaczenie nowych wartości rozwiązania dla węzłów wewnętrznych: Równanie (3), uwzględniając (7) i (8) można zapisać jako:
2 RT i , j −T i−1, j −T i1, j −rT i , j−1−rT i , j1=0 x2 gdzie: r= , R = (1+r). Jeżeli Δx = Δy, to: r = 1, R = 2, oraz: y2 T i , j= T i1, j T i−1, j T i , j−1T i , j1 4
(15)
(16)
4. Kolejne iteracje równania (16) prowadzą do rozwiązania zagadnienia.
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
4
Laboratorium Modelowania i Symulacji 4. Program ćwiczenia.
UWAGA: Wszystkie wykresy powstające po uruchomieniu skryptów wraz z odpowiednimi opisami należy kopiować do programu Wordpad.
4.1 Wykreślanie rozwiązania dokładnego jednowymiarowego, nieustalonego przepływu ciepła
Zapoznać się z zawartością pliku skryptowego d1_dokladne.m Z linii poleceń programu MATLAB wprowadzić polecenie:
>> d1_dokladne
Wyedytować plik d1_dokladne.m. Zmienić chwile czasowe (np.: czas = 0.1, czas = 0.4, czas = 0.7). Zmienić zawartość legendy (polecenie legend()). Powtórnie uruchomić skrypt.
4.2 Wyznaczanie rozwiązania numerycznego jednowymiarowego, nieustalonego zagadnienia przepływu ciepła.
Zagadnienie brzegowe opisano w punkcie 2.1. Wyznaczyć rozwiązanie numeryczne zagadnienia dla następujących danych liczbowych:
J g =4 o cm sec C cm3 J J c=0.012 q ' ' ' =0.08 o 3 gm C sec cm k =2
W tym celu należy wprowadzić następujący skrypt w języku MATLAB:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Skrypt wyznacza rozwiązanie jednowymiarowego, nieustalonego % równania przewodnictwa cieplnego %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % parametry fizykalne (stałe) pii = 3.1415927; k = pii^2; rho = 4; c = 0.01*pii^2; q = 0.08; f = 2/(pii^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % parametry związane z algorytmem obliczeń numerycznych n = 9; % ilość punktów podziału % przestrzennego % (bez elementów brzegowych) h = 1/(n+1) deltaT = 0.001 x = 0:h:1; % krok przestrzenny % krok czasowy % siatka przestrzenna % współczynnik przenikania temperatury % gęstość materiału % ciepło właściwe materiału % wydajność cieplna źródeł wewnętrznych % ciepło Joule'a)
lambda = deltaT/(k*h^2) % parametr równania różnicowego
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
5
Laboratorium Modelowania i Symulacji
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % warunek brzegowy T = zeros(n+2,1000); % warunek początkowy for i=1:n+2 T(i,1) = sin(pi*x(i))+x(i)*(1-x(i)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Iteracje kolejnych kroków czasowych (wskaźnik: j) i % przestrzennych (wskaźnik: i) for j=2:1000 for i=2:n+1 T(i,j)=(1-2*lambda)*T(i,j-1)+lambda*(T(i+1,j-1)+ ... T(i-1,j-1))+deltaT*f; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Jeżeli zeros(..) to warunek zerowy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zapisać skrypt pod nazwą d1_iter.m. Uruchomić skrypt wpisując w oknie poleceń:
>>d1_iter
4.3 Wykreślanie rozwiązania jako wykresu dwuwymiarowego
W oknie poleceń programu MATLAB wprowadzić polecenie:
>>d1_wykres2d(x, [T(:,1), T(:,501), T(:,1000)]);
Wywołanie skryptu powoduje wykreślenie rozkładu temperatur wzdłuż obiektu dla trzech chwil czasowych: t = 0s, t = 0.5s, t = 0.999s
4.4 Wykreślanie rozwiązania jako wykresu trójwymiarowego
W oknie poleceń programu MATLAB wprowadzić polecenia:
[X, Y]=meshgrid(0:deltaT:0.999,x); d1_wykres3d(X,Y,T);
4.5 Zmiana wydajności wewnętrznego źródła ciepła
Wyedytować plik skryptowy d1_iter.m. Zmienić wydajność wewnętrznego źródła ciepła np.:
f = 4/(pii^2); % wydajność cieplna źródeł wewnętrznych % (ciepło Joule'a)
Uruchomić skrypt. Powtórzyć punkty 4.3, 4.4.
4.6 Zmiana warunków brzegowych
Wyedytować plik skryptowy d1_iter.m. Zgodnie ze wskazówkami prowadzącego zmienić warunki brzegowe (temperatury) panujące na końcach rozpatrywanego obiektu (np.:
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
6
Laboratorium Modelowania i Symulacji
T = zeros(n+2,1000)+0.2; % temperatura na brzegach równa 0.2)
Uruchomić plik d1_iter.m. Powtórzyć punkty 4.3, 4.4.
4.7 Zmiana warunku początkowego
Wyedytować plik skryptowy d1_iter.m. Zgodnie ze wskazówkami prowadzącego zmienić warunek początkowy:
% warunek początkowy for i=1:n+2; T(i,1) = 10*(sin(pi*x(i))+x(i)*(1-x(i))); end
Uruchomić skrypt. Powtórzyć punkty 4.3, 4.4.
4.8 Zmiana rozdzielczości siatki przestrzennej
Wyedytować plik skryptowy d1_iter.m. Zgodnie ze wskazówkami prowadzącego zmienić rozdzielczość siatki przestrzennej np.:
n = 3; % ilość punktów podziału przestrzennego
Uruchomić plik d1_iter.m. Powtórzyć punkty 4.3, 4.4
4.9 Wyznaczanie dokładnego rozwiązania dwuwymiarowego, ustalonego zagadnienia brzegowego
Problem dwuwymiarowego, ustalonego przepływu ciepła został przedstawiony w p. 2.2. W celu wizualizacji rozwiązania wprowadzić w oknie poleceń programu MATLAB wprowdzić polecenia:
>>d2_dokladne
4.10 Zmiana rozdzielczości siatki przestrzennej
Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_dokladne.m:
N = 30 % rozdzielczość siatki przestrzennej (siatka kwadratowa)
Uruchomić skrypt.
4.11 Rozwiązanie numeryczne – metoda nadrelaksacji
W celu rozwiązania zagadnienia w oknie poleceń programu MATLAB wprowadzić:
>>d2_iter
Skrypt rozwiązuje rozpatrywane zagadnienie przy zadanych warunkach brzegowych.
4.12 Zmiana liczby iteracji i ocena zbieżności metody
Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_iter.m, np..:
for k = 1:50 % liczba iteracji
Uruchomić skrypt dla k = 50, 200, 500. W celu oceny zbieżności, dla każdej wartości k wprowadzić następujące polecenie:
>>blad_srd = mean(mean(blad)) >>dokl_srd = mean(mean(T_num-T_dokl))
Pierwsze polecenie wyznacza średnią wartość różnicy pomiędzy wartością rozwiązania w kroku numer k i k-1. Jeśli metoda jest zbieżna wartość ta powinna maleć wraz ze wzrostem ilości iteracji. Drugie polecenie wyznacza średnią różnicę pomiędzy wartością rozwiązań: dokładnego i numerycznego. Zanotować wartości zmiennych blad_srd oraz dokl_srd dla poszczególnych k.
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
7
Laboratorium Modelowania i Symulacji 4.13 Zmiana rozdzielczości siatki przestrzennej
UWAGA: w każdym następnym punkcie uruchomienie skryptu d2_iter, powinno być poprzedzone wyczyszczeniem przestrzeni roboczej (polecenie clear). Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_iter.m, np..:
% liczba elementów siatki nx = 30; ny = 30;
Uruchomić skrypt.
4.14 Modyfikacje warunków brzegowych
Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_iter.m, :
for i = 1:nx T(ny,i) = 20*sin(x(i))*sinh(2) % krawędź górna obszaru end
Uruchomić skrypt. Zmienić warunek brzegowy:
for i = 1:nx T(ny,i) = 20 % krawędź górna obszaru end
Uruchomić skrypt. Zmienić warunki brzegowe:
for i = 1:nx T(ny,i) = 0 % krawędź górna obszaru end for i = 1:ny T(i,1) = 20; % krawędź lewa obszaru end
Uruchomić skrypt.
4.15 Wyznaczanie rozwiązania równania Poissona (równanie przewodnictwa przy istnieniu źródeł ciepła wewnątrz obszaru zagadnienia). Metoda nadrelaksacji.
Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_iter.m: 1. Ustawić zerowe warunki początkowe. 2. Ustawić parametry siatki przestrzennej (nx = ny = 30). 3. Ustawić liczbę iteracji k = 1000 4. Włączyć (usunąć komentarz) w skrypcie linie:
% Wymuszenie f = zeros(nx,ny); f(14:16,14:16) = 2000
Powyższe linie ustalają parametry wewnętrznego źródła ciepła (wartość oraz współrzędne przestrzenne). 5. Uruchomić skrypt. 6. Zmodyfikować współrzędne źródła np.:
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
8
Laboratorium Modelowania i Symulacji
f(5:7,4:3) = 2000
7. Uruchomić skrypt. 8. Ustawić warunek początkowy, np.:
for i = 1:ny T(i,1) = 20; % krawędź lewa obszaru end
9. Uruchomić skrypt. Zaobserwować zmiany pola temperatury wynikające z jednoczesnego występowania źródła ciepła i warunku brzegowego.
5. Opracowanie sprawozdania.
W sprawozdaniu należy umieścić wszystkie zarejestrowane wykresy wraz z opisami. Dodatkowo sprawozdanie powinno zawierać wnioski i spostrzeżenia. ZADANIE: Wykorzystując skrypt d2_iter.m, zasymulować model dwuwymiarowego ustalonego przepływu ciepła z wewnętrznymi źródłami ciepła (równanie Poissona) dla czterech źródeł rozmieszczonych w rogach obszaru. Rozmiar źródeł: 2x2 elementy siatki przestrzennej, wydajność źródeł 500 W/m2 (przewodność cieplna λ = 1). Wyniki symulacji dla dwóch siatek (nx=ny=15 i nx=ny=30) wraz z komentarzem zamieścić w sprawozdaniu.
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
9
Modelowanie i Symulacja, Modelowanie obiektów o stałych rozłożonych
Laboratorium Modelowania i Symulacji Zespół Automatyki ZTMAiPC Wydział Elektryczny
LABORATORIUM MODELOWANIA I SYMULACJI
Ćwiczenie 8 Modelowanie obiektów o stałych rozłożonych
1. Cel ćwiczenia
Celem ćwiczenia jest zapoznanie się ze sposobami modelowania układów o stałych rozłożonych. W instrukcji zawarto podstawowe informacje na temat metody różnic skończonych wykorzystywanej do rozwiązywania jednowymiarowego zagadnienia przepływu ciepła. Dodatkowo przedstawiono metodę nadrelaksacji na przykładzie równania Laplace'a.
2. Obiekty o stałych rozłożonych 2.1 Jednowymiarowy, nieustalony przepływ ciepła w ciele stałym. Warunki początkowe i brzegowe
Rys. 1 przedstawia model obiektu o stałych rozłożonych. Modelowanym zjawiskiem jest ruch ciepła w pręcie o długości L = 1, przez który przepływa prąd elektryczny. Przepływ prądu powoduje generowanie ciepła wewnątrz pręta. Oba końce pręta są zanurzone w temperaturze T = 0. W tym przypadku pomimo tego, że obiekt modelowany ma trzy wymiary geometryczne zagadnienie określa się mianem jednowymiarowego. Jest to związane z faktem, że przepływ ciepła odbywa się tylko w jednym kierunku (od wnętrza walca na zewnątrz). Warunki brzegowe T(0,t) oraz T(1,T) określają wartość temperatury na końcach pręta (brzegach obszaru). Oprócz warunków brzegowych dla rozwiązania zagadnienia opisywanego cząstkowym równaniem różniczkowym niezbędna jest znajomość warunku początkowego. Podobnie jak w zagadnieniach opisywanych równaniami zwyczajnymi, warunki początkowe określają, jakie wartości musi przyjmować rozwiązanie w chwili t = 0.
Warunki brzegowe
T 0, t = 0
T 1, t = 0
Warunek początkowy
T x ,0 = sin x x 1 − x
Rysunek 1. Geometria, warunki brzegowe i warunek początkowy
Do opisu matematycznego opisywanego zjawiska stosuje się jednowymiarowe równanie nieustalonego przepływu ciepła:
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
1
Laboratorium Modelowania i Symulacji ∂ T x , t 1 ∂2 T x , t − = f x ,t, ∂t k ∂ x2
(1)
podawane wraz z warunkami jednoznaczności rozwiązania (brzegowymi i początkowym). W równaniu (1): T(x,t) – wartość temperatury w punkcie x, w chwili czasu równej t; k – współczynnik dyfuzyjności cieplnej (zależny od rodzaju materiału); f(x,t) – funkcja wydajności cieplnej objętościowych źródeł mocy. Można wykazać, że przy podanych warunkach brzegowych równanie (1) posiada następujące rozwiązanie:
T x , t =e t sin x x 1−x 2.2 Dwuwymiarowy, ustalony przepływ ciepła w ciele stałym
(2)
Na Rys. 2 przedstawiono dwuwymiarowe zagadnienie ustalonego przepływu ciepła wraz z warunkami brzegowymi.
T x ,2
Warunek brzegowy
T x ,2 = sin x sinh 2
x
Brzeg obszaru
T , y = 0
T 0, y =0
Element siatki przestrzennej
T x ,0 = 0
Rysunek 2. Dwuwymiarowy, ustalony przepływ ciepła
Zagadnienie opisuje równanie Laplace'a:
∂ T x , y ∂ T x , y =0. ∂ x2 ∂ y2
(3)
Problem brzegowy, ustalony nie wymaga podawania warunku początkowego. Jeżeli w równaniu (3) prawa strona zostanie zastąpiona niezerową funkcją, równanie nazywamy równaniem Poissona. Rozwiązanie równania (3) jest uzależnione od warunków brzegowych panujących na granicach (brzegach) rozpatrywanego obszaru. Przy uwzględnieniu warunków brzegowych jak na rys. 2 ma postać:
T x , y=sin xsinh y
(4)
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
2
Laboratorium Modelowania i Symulacji 3. Podstawy metody różnic skończonych 3.1 Wprowadzenie
Do rozwiązywania zagadnień brzegowych jedno i wielowymiarowych często stosuje się metodę różnic skończonych opartą na aproksymacji pochodnej funkcji ilorazem różnicowym „w przód”. Dla uproszczenia rozpatrzmy pochodną funkcji jednej zmiennej:
f ' x 0 ≃
f x 0 h− f x 0 h
(5)
gdzie: h oznacza niewielki rozmiar kroku. Aby wprowadzić pojęcie rzędu aproksymacji rozwińmy f x 0 h w szereg Taylora wokół punktu x0:
f x 0h− f x 0 h h2 = f ' x 0 f ' ' x 0 f ' ' ' x 0 RN h 2 6
(6)
gdzie RN oznacza składniki rozwinięcia wyższych rzędów. Przyjmując, że h jest bardzo małe, błąd popełniany przez stosowanie równania (5) jest w przybliżeniu równy
h f ' ' x 0 . Mówimy, że błąd aproksymacji w równaniu (5) jest rzędu h i oznaczamy go 2
O(h). Tego rodzaju aproksymacja nie jest satysfakcjonująca. Dla przykładu, dokładność 0.001, będzie wymagała bardzo małego kroku h. Z tego powodu wymagana jest aproksymacja rzędu h2 lub wyższego. Poniżej przedstawiono aproksymację rzędu h2, opartą na centralnym ilorazie różnicowym:
f ' x 0 ≃
f x 0h− f x 0 −h . 2h
(7)
W późniejszych rozważaniach niezbędna będzie również znajomość przybliżenia f''(x0):
f ' ' x 0 ≃
f x 0 h−2 f x 0 f x 0 −h h2
.
(8)
3.2 Aproksymacja rozwiązania jednowymiarowego równania przewodzenia ciepła
Poszukujemy rozwiązania równania (1) dla zadanych warunków brzegowych. Wprowadźmy oznaczenia: h=x i1− x i oraz t=t j1−t j jako odległość pomiędzy wartościami rozwiązania w punktach (xi, yj) (współrzędne siatki czasowo-przestrzennej). Oznaczając: T i , j =T x i , t j można równanie (1) zapisać jako:
T i , j1−T i , j 1 T i1, j −2T i , j T i− j , j = f i, j. k h2 t
Przekształcając (9) ze względu na Ti,j otrzymuje się:
(9)
T i , j1=1−2 T i , j T i1, j T i−1, j tf i , j .
gdzie: =
(10)
t , i = 1, 2, ... n kh 2
Wygodnie jest zapisać równanie (10) w postaci macierzowej. Przyjmując, że Tj oznacza wektor (T1j, T2j, ..., Tnj): Tj+1 = ATj + fj, gdzie:
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
3
Laboratorium Modelowania i Symulacji 1−2 1−2 A= 0 ⋱ ⋮ ⋱ 0 ⋯ 3.3 Metoda Cranka – Nicolsona.
0 ⋱ ⋱ 0
⋯ 0 ⋯ 0 ⋱ ⋮ ⋱ 1−2
(11)
W celu poprawy stabilności numerycznej oraz zwiększenia dokładności obliczeń, często stosuje się metodę Cranka – Nicolsona (C-N). W metodzie tej wprowadza się macierz B:
2 −1 B= 0 ⋮ 0
Oznaczając przez I – macierz jednostkową, można wprowadzić następujące macierze:
−1 2 ⋱ ⋱ ⋯
0 −1 ⋱ ⋱ 0
⋯ ⋯ ⋱ ⋱ −1
0 0 ⋮ −1 2
(12)
C =2 I B , D=2 I − B
zagadnienie brzegowe rozwiązuje się zgodnie z równaniem macierzowym:
(13) (14)
CT j1 =DT j 2 tf
j
3.4 Aproksymacja rozwiązania dwuwymiarowego zagadnienia ustalonego przepływu ciepła – metoda nadrelaksacji
Zagadnienie przedstawione w punkcie 2.2, zostanie rozwiązane metodą nadrelaksacji. W metodzie tej wyróżnić można następujące etapy: 1. Podział obszaru zagadnienia (utworzenie siatki) o odległościach pomiędzy węzłami odpowiednio: Δx w poziomie oraz Δy w pionie. 2. Przypisanie warunków brzegowych węzłom brzegowym, a węzłom wewnętrznym dowolnych wartości początkowych. 3. Wyznaczenie nowych wartości rozwiązania dla węzłów wewnętrznych: Równanie (3), uwzględniając (7) i (8) można zapisać jako:
2 RT i , j −T i−1, j −T i1, j −rT i , j−1−rT i , j1=0 x2 gdzie: r= , R = (1+r). Jeżeli Δx = Δy, to: r = 1, R = 2, oraz: y2 T i , j= T i1, j T i−1, j T i , j−1T i , j1 4
(15)
(16)
4. Kolejne iteracje równania (16) prowadzą do rozwiązania zagadnienia.
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
4
Laboratorium Modelowania i Symulacji 4. Program ćwiczenia.
UWAGA: Wszystkie wykresy powstające po uruchomieniu skryptów wraz z odpowiednimi opisami należy kopiować do programu Wordpad.
4.1 Wykreślanie rozwiązania dokładnego jednowymiarowego, nieustalonego przepływu ciepła
Zapoznać się z zawartością pliku skryptowego d1_dokladne.m Z linii poleceń programu MATLAB wprowadzić polecenie:
>> d1_dokladne
Wyedytować plik d1_dokladne.m. Zmienić chwile czasowe (np.: czas = 0.1, czas = 0.4, czas = 0.7). Zmienić zawartość legendy (polecenie legend()). Powtórnie uruchomić skrypt.
4.2 Wyznaczanie rozwiązania numerycznego jednowymiarowego, nieustalonego zagadnienia przepływu ciepła.
Zagadnienie brzegowe opisano w punkcie 2.1. Wyznaczyć rozwiązanie numeryczne zagadnienia dla następujących danych liczbowych:
J g =4 o cm sec C cm3 J J c=0.012 q ' ' ' =0.08 o 3 gm C sec cm k =2
W tym celu należy wprowadzić następujący skrypt w języku MATLAB:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Skrypt wyznacza rozwiązanie jednowymiarowego, nieustalonego % równania przewodnictwa cieplnego %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % parametry fizykalne (stałe) pii = 3.1415927; k = pii^2; rho = 4; c = 0.01*pii^2; q = 0.08; f = 2/(pii^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % parametry związane z algorytmem obliczeń numerycznych n = 9; % ilość punktów podziału % przestrzennego % (bez elementów brzegowych) h = 1/(n+1) deltaT = 0.001 x = 0:h:1; % krok przestrzenny % krok czasowy % siatka przestrzenna % współczynnik przenikania temperatury % gęstość materiału % ciepło właściwe materiału % wydajność cieplna źródeł wewnętrznych % ciepło Joule'a)
lambda = deltaT/(k*h^2) % parametr równania różnicowego
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
5
Laboratorium Modelowania i Symulacji
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % warunek brzegowy T = zeros(n+2,1000); % warunek początkowy for i=1:n+2 T(i,1) = sin(pi*x(i))+x(i)*(1-x(i)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Iteracje kolejnych kroków czasowych (wskaźnik: j) i % przestrzennych (wskaźnik: i) for j=2:1000 for i=2:n+1 T(i,j)=(1-2*lambda)*T(i,j-1)+lambda*(T(i+1,j-1)+ ... T(i-1,j-1))+deltaT*f; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Jeżeli zeros(..) to warunek zerowy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zapisać skrypt pod nazwą d1_iter.m. Uruchomić skrypt wpisując w oknie poleceń:
>>d1_iter
4.3 Wykreślanie rozwiązania jako wykresu dwuwymiarowego
W oknie poleceń programu MATLAB wprowadzić polecenie:
>>d1_wykres2d(x, [T(:,1), T(:,501), T(:,1000)]);
Wywołanie skryptu powoduje wykreślenie rozkładu temperatur wzdłuż obiektu dla trzech chwil czasowych: t = 0s, t = 0.5s, t = 0.999s
4.4 Wykreślanie rozwiązania jako wykresu trójwymiarowego
W oknie poleceń programu MATLAB wprowadzić polecenia:
[X, Y]=meshgrid(0:deltaT:0.999,x); d1_wykres3d(X,Y,T);
4.5 Zmiana wydajności wewnętrznego źródła ciepła
Wyedytować plik skryptowy d1_iter.m. Zmienić wydajność wewnętrznego źródła ciepła np.:
f = 4/(pii^2); % wydajność cieplna źródeł wewnętrznych % (ciepło Joule'a)
Uruchomić skrypt. Powtórzyć punkty 4.3, 4.4.
4.6 Zmiana warunków brzegowych
Wyedytować plik skryptowy d1_iter.m. Zgodnie ze wskazówkami prowadzącego zmienić warunki brzegowe (temperatury) panujące na końcach rozpatrywanego obiektu (np.:
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
6
Laboratorium Modelowania i Symulacji
T = zeros(n+2,1000)+0.2; % temperatura na brzegach równa 0.2)
Uruchomić plik d1_iter.m. Powtórzyć punkty 4.3, 4.4.
4.7 Zmiana warunku początkowego
Wyedytować plik skryptowy d1_iter.m. Zgodnie ze wskazówkami prowadzącego zmienić warunek początkowy:
% warunek początkowy for i=1:n+2; T(i,1) = 10*(sin(pi*x(i))+x(i)*(1-x(i))); end
Uruchomić skrypt. Powtórzyć punkty 4.3, 4.4.
4.8 Zmiana rozdzielczości siatki przestrzennej
Wyedytować plik skryptowy d1_iter.m. Zgodnie ze wskazówkami prowadzącego zmienić rozdzielczość siatki przestrzennej np.:
n = 3; % ilość punktów podziału przestrzennego
Uruchomić plik d1_iter.m. Powtórzyć punkty 4.3, 4.4
4.9 Wyznaczanie dokładnego rozwiązania dwuwymiarowego, ustalonego zagadnienia brzegowego
Problem dwuwymiarowego, ustalonego przepływu ciepła został przedstawiony w p. 2.2. W celu wizualizacji rozwiązania wprowadzić w oknie poleceń programu MATLAB wprowdzić polecenia:
>>d2_dokladne
4.10 Zmiana rozdzielczości siatki przestrzennej
Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_dokladne.m:
N = 30 % rozdzielczość siatki przestrzennej (siatka kwadratowa)
Uruchomić skrypt.
4.11 Rozwiązanie numeryczne – metoda nadrelaksacji
W celu rozwiązania zagadnienia w oknie poleceń programu MATLAB wprowadzić:
>>d2_iter
Skrypt rozwiązuje rozpatrywane zagadnienie przy zadanych warunkach brzegowych.
4.12 Zmiana liczby iteracji i ocena zbieżności metody
Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_iter.m, np..:
for k = 1:50 % liczba iteracji
Uruchomić skrypt dla k = 50, 200, 500. W celu oceny zbieżności, dla każdej wartości k wprowadzić następujące polecenie:
>>blad_srd = mean(mean(blad)) >>dokl_srd = mean(mean(T_num-T_dokl))
Pierwsze polecenie wyznacza średnią wartość różnicy pomiędzy wartością rozwiązania w kroku numer k i k-1. Jeśli metoda jest zbieżna wartość ta powinna maleć wraz ze wzrostem ilości iteracji. Drugie polecenie wyznacza średnią różnicę pomiędzy wartością rozwiązań: dokładnego i numerycznego. Zanotować wartości zmiennych blad_srd oraz dokl_srd dla poszczególnych k.
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
7
Laboratorium Modelowania i Symulacji 4.13 Zmiana rozdzielczości siatki przestrzennej
UWAGA: w każdym następnym punkcie uruchomienie skryptu d2_iter, powinno być poprzedzone wyczyszczeniem przestrzeni roboczej (polecenie clear). Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_iter.m, np..:
% liczba elementów siatki nx = 30; ny = 30;
Uruchomić skrypt.
4.14 Modyfikacje warunków brzegowych
Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_iter.m, :
for i = 1:nx T(ny,i) = 20*sin(x(i))*sinh(2) % krawędź górna obszaru end
Uruchomić skrypt. Zmienić warunek brzegowy:
for i = 1:nx T(ny,i) = 20 % krawędź górna obszaru end
Uruchomić skrypt. Zmienić warunki brzegowe:
for i = 1:nx T(ny,i) = 0 % krawędź górna obszaru end for i = 1:ny T(i,1) = 20; % krawędź lewa obszaru end
Uruchomić skrypt.
4.15 Wyznaczanie rozwiązania równania Poissona (równanie przewodnictwa przy istnieniu źródeł ciepła wewnątrz obszaru zagadnienia). Metoda nadrelaksacji.
Zgodnie ze wskazówkami prowadzącego zmodyfikować skrypt d2_iter.m: 1. Ustawić zerowe warunki początkowe. 2. Ustawić parametry siatki przestrzennej (nx = ny = 30). 3. Ustawić liczbę iteracji k = 1000 4. Włączyć (usunąć komentarz) w skrypcie linie:
% Wymuszenie f = zeros(nx,ny); f(14:16,14:16) = 2000
Powyższe linie ustalają parametry wewnętrznego źródła ciepła (wartość oraz współrzędne przestrzenne). 5. Uruchomić skrypt. 6. Zmodyfikować współrzędne źródła np.:
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
8
Laboratorium Modelowania i Symulacji
f(5:7,4:3) = 2000
7. Uruchomić skrypt. 8. Ustawić warunek początkowy, np.:
for i = 1:ny T(i,1) = 20; % krawędź lewa obszaru end
9. Uruchomić skrypt. Zaobserwować zmiany pola temperatury wynikające z jednoczesnego występowania źródła ciepła i warunku brzegowego.
5. Opracowanie sprawozdania.
W sprawozdaniu należy umieścić wszystkie zarejestrowane wykresy wraz z opisami. Dodatkowo sprawozdanie powinno zawierać wnioski i spostrzeżenia. ZADANIE: Wykorzystując skrypt d2_iter.m, zasymulować model dwuwymiarowego ustalonego przepływu ciepła z wewnętrznymi źródłami ciepła (równanie Poissona) dla czterech źródeł rozmieszczonych w rogach obszaru. Rozmiar źródeł: 2x2 elementy siatki przestrzennej, wydajność źródeł 500 W/m2 (przewodność cieplna λ = 1). Wyniki symulacji dla dwóch siatek (nx=ny=15 i nx=ny=30) wraz z komentarzem zamieścić w sprawozdaniu.
Ćwiczenie 8 – Modelowanie obiektów o stałych rozłożonych
9
