Ćwiczenie 6a:Równanie dyfuzji w dwóch wymiarach ∗

(1)

Ćwiczenie 6a:Równanie dyfuzji w dwóch wymiarach

^∗

j_min j₁ j₂ j_max

i_min i₁ i₂ i₃ i₄ i_max

T₀

T₁ T₂

A

B

C

D E F

G H I J

indeks i

indeks j

Rysunek 1: Siatka obliczeniowa do opisu rozkładu temperatury wewnątrz korytarza (szary obszar), który wymienia ciepło z pomieszczeniami sąsiednimi (czerwony, niebieski) oraz z otoczeniem (przez zaznaczone na zielono ściany).

Sąsiednie pomieszczenie o temperaturze T1 graniczy z korytarzem wzdłuż czer- wonej linii, natomiast niebieską linią oznaczono sąsiednie pomieszczenie o tem- peraturze T₂. Na zewnątrz panuje temperatura T₀.

Zajmiemy się problemem rozkładu temperatury i bilansu cieplnego koryta- rza, prowadzącego z pomieszczenia o temperaturze T1 do pomieszczenia o tem- peraturze T2 (rys. 1). Kratka odpowiada siatce numerycznej. W rachunkach zajmiemy się czystą dyfuzją ciepła z całkowitym zaniedbaniem innych mecha- nizmów transportu (tzn. wymiana z otoczeniem będzie przez konwekcję, ale ta

∗Laboratorium z inżynierskich metod numerycznych, Wydział Fizyki i Informatyki Sto- sowanej AGH 2014/2015. Bartłomiej Szafran (bszafran@agh.edu.pl), Krzysztof Kolasiński (kolasinski@fis.agh.edu.pl), Elżbieta Wach (Elzbieta.Wach@fis.agh.edu.pl), Dariusz Żebrowski (Dariusz.Zebrowski@fis.agh.edu.pl)

1

(2)

znajdzie się tylko w warunkach brzegowych). Na korytarzu nie ma źródeł ciepła.

Temperatura powietrza T (x, y, t) opisana jest równaniem przewodnictwa cieplnego (dyfuzji ciepła)

∂T

∂t = k

ρc∇²T, (1)

gdzie ρ, c oraz k są odpowiednio gęstością powietrza, jego ciepłem właściwym i współczynnikiem przewodności termicznej.

Siatka i schemat różnicowy

Ustalamy imin = 0, i1 = 20, i2 = 40, i3 = 60, i4 = 80, imax = 100 oraz jmin = 0, j1 = 20, j2 = 30, jmax = 50. Pracujemy na siatce z krokiem przestrzennym ∆x = ∆y = 1. Przyjmujemy również k = 1, ρ = 1, c = 1.

Równanie (1) rozwiązywać będziemy iteracyjnym schematem Cranka-Nicolson.

Rozkład temperatury w n + 1 chwili czasowej (T_ijⁿ⁺¹) wyliczamy na podstawie rozkładu temperatury w kroku n-tym (potrzebujemy więc dwie tablice):

T_i,jⁿ⁺¹− Ti,jⁿ

∆t = k∇²T_i,jⁿ⁺¹

2ρc +k∇²T_i,jⁿ

2ρc . (2)

Rozpisując:

T_i,jⁿ⁺¹= T_i,jⁿ + k∆t

2ρc∆x²[T_i,jⁿ₋₁+ T_i,j+1ⁿ + T_iⁿ_−1,j+ T_i+1,jⁿ − 4Ti,jⁿ

+T_i,jⁿ⁺¹₋₁+ T_i,j+1ⁿ⁺¹ + T_iⁿ⁺¹_−1,j+ T_i+1,jⁿ⁺¹ − 4T_i,jⁿ⁺¹]. (3) Po przekształceniu otrzymujemy ostatecznie przepis iteracyjny dla kroku w chwili n+1:

T_i,jⁿ⁺¹= 1 1 + _2ρc∆x^4k∆t2

(T_i,jⁿ + k∆t

2ρc∆x²[T_i,jⁿ₋₁+ T_i,j+1ⁿ + T_iⁿ_−1,j+ T_i+1,jⁿ − 4Ti,jⁿ

+T_i,jⁿ⁺¹₋₁+ T_i,j+1ⁿ⁺¹ + T_iⁿ⁺¹_−1,j+ T_i+1,jⁿ⁺¹]). (4) Iteracje prowadzimy aż do uzyskania zbieżności (około 50 iteracji dla po- czątkowych kroków czasowych - liczba zależna od sposobu badania zbieżno- ści). W dalszych chwilach czasowych zbieżność będzie uzyskiwana bardzo szyb- ko, ponieważ rozkład temperatury przestanie się zmieniać (rozpatrywany układ znajdzie się w równowadze termicznej). Zbieżność można zbadać np. obliczając całkę z wartości bezwzględnej temperatury po całym pudle obliczeniowym - je- śli przestanie się znacząco zmieniać, możemy stwierdzić zbieżność dla danego t.

Na podstawie tak obliczonego rozkładu temperatury przechodzimy do kolejnego t := t + dt.

Warunki brzegowe

Korytarz wymienia ciepło z pomieszczeniem niebieskim i czerwonym (patrz rysunek 1). W czerwonym pokoju panuje T1= 21, w niebieskim T2= T1, na ze- wnątrz T0= 0. Dla punktów granicznych położonych w sąsiedztwie czerwonego lub niebieskiego wstawiamy na brzegu stałą temperaturę Tij = T1lub Tij = T2.

2

(3)

Na zewnętrznych ścianach korytarza (zaznaczone kolorem zielonym na ry- sunku 1) zadajemy konwekcyjne warunki brzegowe

h(T − T 0) = −k∂T

∂n, (5)

gdzie pochodna po prawej stronie oznacza normalną do powierzchni budynku, ale liczoną wewnątrz korytarza, a k jest współczynnikiem transmisji ciepła. Wa- runek (5) dla pionowych krawędziI (i = i2) iD (i = i4) ma postać

T_i,jⁿ⁺¹=

k

∆xT_iⁿ⁺¹_−1,j+ hT0

h +_∆x^k . (6)

Podobnie dla pozostałych pionowych krawędzi -B(i = i₁),G(i = i₃):

T_i,jⁿ⁺¹=

k

∆xT_i+1,jⁿ⁺¹ + hT0

h +_∆x^k . (7)

Na poziomych dolnych -A, E (j = j₂) orazC (j = j_min):

T_i,jⁿ⁺¹=

k

∆yT_i,j+1ⁿ⁺¹ + hT0

h +_∆y^k , (8)

na górnych -H (j = j₁) orazJ, F (j = j_max):

T_i,jⁿ⁺¹=

k

∆yT_i,jⁿ⁺¹₋₁+ hT0

h +_∆y^k . (9)

Na kantach wyżej wymienionych odcinków pochodne normalne policzymy po antydiagonali. Dlatego na narożnikachJ/I orazI/H, tj. (i₂, j_max) i (i₂, j₁):

T_i,jⁿ⁺¹=

√k

2∆xT_iⁿ⁺¹_−1,j−1+ hT0

h +√^k 2∆x

, (10)

...dlaF/G,G/H, tj. (i₃, j_max), (i₃, j₁):

T_i,jⁿ⁺¹=

√k

2∆xT_i+1,jⁿ⁺¹₋₁+ hT0

h +√^k 2∆x

, (11)

...dlaE/D,D/C, tj. (i₄, j₂), (i₄, j_min):

T_i,jⁿ⁺¹=

√k

2∆xT_iⁿ⁺¹_−1,j+1+ hT0

h +√^k 2∆x

, (12)

orazA/B iB/C, tj. dla punktów (i1, j2) i (i1, jmin):

T_i,jⁿ⁺¹=

√k

2∆xT_i+1,j+1ⁿ⁺¹ + hT0

h +√^k 2∆x

. (13)

Warunki brzegowe (6) - (13) są zmienne - należy je zadawać na początku każdej iteracji.

3

(4)

Zadanie 1: Idealna izolacja

Przyjmiemy ∆t = 2.5. Wstawić h = 0 (idealna izolacja termiczna ścian).

W chwili początkowej na korytarzu T = T₀. Prześledzić, jak korytarz ogrzewa się od sąsiadów: sprawdzić, ile ciepła wpływa od czerwonego i niebieskiego sąsiada (w funkcji czasu). Definiujemy strumień ciepła:

⃗

q =−k⃗∇T =−k





∂T

∂x

∂T

∂y



 . (14)

Dla czerwonego sąsiada: liczymy składową poziomą (zaznaczoną na czerwono we wzorze (14)) i całkujemy ją wzdłuż linii i = imin+ 1 [dla j∈ (j2, jmax)]. Dla niebieskiego sąsiada postępujemy identycznie, zmieniamy tylko znak wartości strumienia na przeciwny oraz sumujemy z odpowiednimi indeksami: i = i_max−1 dla j ∈ (j2, j_max). W stanie ustalonym bilans transferów powinien wyjść na zero. Narysować zależny od czasu strumień ciepła od sąsiadów oraz ich sumę (30 pkt). Narysować rozkład temperatury na korytarzu dla stanu ustalonego – gdy bilans wyjdzie na zero (30 pkt).

Zadanie 2: Nieszczelność

Ciepło ucieka przez ściany: wstawić h = 0.01. Powtórzyć algorytm. Naryso- wać pięć map rozkładu temperatury T aż do stanu ustalonego (40 pkt). Czy izotermy pozostaną prostopadłe do ścian?

4