Projekt 9: Dyfuzja ciepła - metoda Cranck-Nicloson.

(1)

Projekt 9: Dyfuzja ciepła - metoda Cranck-Nicloson.

Tomasz Chwiej 10 stycznia 2019

1 Wstęp

i= 0 1

n_x 2

0 1 2 n_y

j=

ρ

₁

Rysunek 1: Siatka węzłów użyta w obliczeniach z zaznaczonymi warunkami brzegowymi: Dirichleta (czerwony) i von Neumanna (niebieski).

Na zajęciach znajdziemy rozwiązanie zależnego od czasu równania dyfuzji (T = T (x, y, t))

∇

²

T = ∂T

∂t (1)

Równanie to zdyskretyzujemy na siatce, zapiszemy w postaci macierzowej i znajdziemy jego rozwią- zanie w kolejnych chwilach czasowych przy użyciu metody Cranck-Nicolson (do rozwiązania układu równań liniowych w tej metodzie wykorzystamy rozkład LU). Geometria układu wraz z siatką węzłów, w których wyznaczona zostanie temperatura pokazane jest na rys. 1.

1.1 Dyskretyzacja równania + metoda CN

Najperw definiujemy siatkę węzłów, dla położenia (x,y) i czasu (t)

∆x = ∆y = ∆ (2)

x

i

= i · ∆, i = 0, 1, 2, . . . , n

x

(3) y

_j

= j · ∆, j = 0, 1, 2, . . . , n

y

(4)

t

_n

= n · ∆t, n = 0, 1, 2, . . . (5)

T (x, y, t) → T (x

i

, y

j

, t

n

) → T

i,jⁿ

(6) (dolny indeks - położenie, górny - czas).

Zapisujemy równanie (1) korzystając ze schematu C-N 1

2 ( ∇

²

T

ⁿ

+ ∇

²

T

ⁿ⁺¹

)

= T

ⁿ⁺¹

− T

ⁿ

∆t (7)

(2)

zastępujemy pochodne ilorazami różnicowymi ^[ d

²

f /dx

²

= (f

_i+1

− 2f

i

+ f

_i₋₁

)/∆

²

^] i grupujemy wyrazy względem chwil czasowych

∆t 2∆

²

(

T

_i+1,jⁿ⁺¹

+ T

_iⁿ⁺¹_−1,j

− 4T

_i,jⁿ⁺¹

+ T

_i,j+1ⁿ⁺¹

+ T

_i,jⁿ⁺¹₋₁

) − T

_i,jⁿ⁺¹

= − ∆t 2∆

²

(

T

_i+1,jⁿ

+ T

_iⁿ_−1,j

− 4T

_i,jⁿ

+ T

_i,j+1ⁿ

+ T

_i,jⁿ₋₁

) − T

_i,jⁿ

(8)

Dokonujemy teraz reindeksacji węzłów (parę wskaźników (i, j) zastępujemy jednym l)

l = i + j · (n

x

+ 1), l = 0, 1, 2, . . . , N (9)

N = (n

x

+ 1) · (n

y

+ 1) (10)

j = f loor ( l

n

x

+ 1 )

(11)

i = l − j · (n

x

+ 1) (12)

Równanie 8 możemy zapisać w postaci macierzowej

A · ⃗T

ⁿ⁺¹

= B · ⃗T

ⁿ

+ ⃗ c (13)

[uwaga: pojawienie się wektora ⃗ c nie wynika bezpośrednio z równania (8), ale dodaliśmy je aby zapewnić spełnienie warunków brzegowych von Neumanna, dla WB Dirichleta przyjmujemy ⃗ c = 0. (góra/dół na rysunku 1],

oraz z uwzględnieniem tylko niezerowych elementów macierzowych

a

_l,l_−nx−1

T

_lⁿ⁺¹_−nx−1

+ a

_l,l₋₁

T

_lⁿ⁺¹₋₁

+ a

_l,l

T

_lⁿ⁺¹

+ a

_l,l+1

T

_l+1ⁿ⁺¹

+ a

_l,l+n_x₊₁

T

_l+nⁿ⁺¹

x+1

= b

l,l−nx−1

T

_lⁿ_−nx−1

+ b

l,l−1

T

_lⁿ₋₁

+ b

l,l

T

_lⁿ

+ b

l,l+1

T

_l+1ⁿ

+ b

l,l+nx+1

T

_l+nⁿ _x₊₁

+ c

l

(14) 1.1.1 Niezerowe elementy macierzowe z uwzględnieniem WB

• Wnętrze obszaru (szary obszar na rysunku)

i = 1, 2, . . . , n

_x

− 1 (15)

j = 1, 2, . . . , n

_y

− 1 (16)

a

_l,l_−n_x₋₁

= a

_l,l₋₁

= a

_l,l+1

= a

_l,l+n_x₊₁

= ∆t

2∆

²

(17)

a

l,l

= − 2∆t

∆

²

− 1 (18)

b

_l,l_−n_x₋₁

= b

_l,l₋₁

= b

_l,l+1

= b

_l,l+n_x₊₁

= − ∆t

2∆

²

(19)

b

l,l

= 2∆t

∆

²

− 1 (20)

• WB Dirichleta (lewy i prawy brzeg)

i = 0, n

_x

(21)

j = 0, 1, 2, . . . , n

_y

(22)

a

_l,l

= 1 (23)

b

l,l

= 1 (24)

c

l

= 0 (25)

• WB von Neumanna na górnym brzegu dla chwili n + 1

∂T

ⁿ⁺¹

∂y = −k

B

(T

ⁿ⁺¹

− T

B

) (26)

(3)

po zastąpieniu pochodnych ilorazami jest następujący T

_i,nⁿ⁺¹

y

− T

_i,nⁿ⁺¹_y₋₁

∆ = −k

B

(T

_i,nⁿ⁺¹_y

− T

B

) (27)

a po reindeksacji węzłów i po pogrupowaniu wyrazów a

_l,l−nx−1

T

_lⁿ⁺¹_−n

x−1

+ a

_l,l

T

_lⁿ⁺¹

= c

_l

(28)

(29) gdzie

i = 1, 2, . . . , n

_x

− 1 (30)

j = n

_y

(31)

a

_l,l_−n_x₋₁

= − 1

k

B

∆ (32)

a

_l,l

= 1 + 1

k

_B

∆ (33)

c

_l

= T

B

(wektor c) (34)

b

l,∗

= 0, (cały wiersz) (35)

• WB von Neumanna na dolnym brzegu dla chwili n + 1

∂T

ⁿ⁺¹

∂y = −k

D

(T

ⁿ⁺¹

− T

D

) (36)

po zastąpieniu pochodnych ilorazami jest następujący T

_i,0ⁿ⁺¹

− T

_i,1ⁿ⁺¹

∆ = −k

D

(T

_i,0ⁿ⁺¹

− T

B

) (37)

a po reindeksacji węzłów i po pogrupowaniu wyrazów a

_l,l

T

_lⁿ⁺¹

+ a

_l,l+nx+1

T

_l+nⁿ⁺¹

x+1

= c

_l

(38)

(39) gdzie

i = 1, 2, . . . , n

_x

− 1 (40)

j = 0 (41)

a

_l,l

= 1 + 1

k

_D

∆ (42)

a

_l,l+n_x₊₁

= − 1

k

_D

∆ (43)

c

_l

= T

_D

(wektor c) (44)

b

_l,_∗

= 0, (cały wiersz) (45)

1.1.2 Warunki Początkowe WP narzucamy na wektor startowy ⃗ T

⁰

T

_l⁰

= T

_A

, i = 0, j = 0, 1, 2 . . . , n

_y

, (lewy brzeg) (46) T

_l⁰

= T

_C

, i = n

x

, j = 0, 1, 2 . . . , n

y

, (prawy brzeg) (47)

T

_l⁰

= 0, w pozostałym obszarze (48)

(4)

2 Algorytm CN

Algorytm CN dla równania dyfuzji w wersji macierzowej i n i c j a l i z a c j a : A , B , ⃗ c, ⃗ T = ⃗ T

⁰

o b l i c z r o z k ł a d LU : A

FOR it =0 TO I T _ M A X S T E P 1 DO d = B ⃗ · ⃗T + ⃗c //T = T

ⁿ

r o z w i ą ż ukł . r ó w n ( LU ): A · ⃗T = ⃗d //T = T

ⁿ⁺¹

END DO

3 Zadania do wykonania

W obliczeniach do znalezienia rozkładu LU, rozwiązania układu równań liniowych, mnożenia macierz- wektor czy dodawania wektorów należy użyć biblioteki numerycznej np. GSL (lub innej dostępnej na Taurusie).

1. W obliczeniach należy użyć wartości parametrów: n

_x

= 40, n

_y

= 40, N = (n

_x

+ 1) · (n

y

+ 1) ,

∆ = 1, ∆t = 1, T

_A

= 40, T

_B

= 0, T

_C

= 30, T

_D

= 0, k

_B

= 0.1, k

_D

= 0.6, IT M AX = 2000.

2. Utworzyć macierze A[N ][N ], B[N ][N ] oraz wektor c[N ] i wypełnić je zgodnie z wzorami zamiesz- czonymi w sekcji 1.1.1

3. Utworzyć wektor startowy T [N ] i narzucić warunki początkowe (sekcja 1.1.2) 4. Znaleźć rozkład LU macierzy A (GSL)

5. Zaimplementować algorytm CN (sekcja 2) 6. Wykonać IT M AX kroków.

7. Dla it = 100, 200, 500, 1000, 2000 sporządzić mapy rozkładu temperatury w pomieszczeniu T (x, y) (50 pkt.)

8. W stanie ustalonym równanie dyfuzji redukuje się do ∇

²

T = 0 (brak zmian w czasie - znika po- chodna czasowa). Dla it = 100, 200, 500, 1000, 2000 proszę sporządzić mapy rozkładu ∇

²

T (x, y) i sprawdzić czy tak jest. (50 pkt.)

4 Przydatne procedury z biblioteki GSL

• rozkład LU

int g s l _ l i n a l g _ L U _ d e c o m p ( g s l _ m a t r i x * A , g s l _ p e r m u t a t i o n * p , int * s i g n u m ) Po wywołaniu procedury znaleziony rozkład LU jest wpisany do macierzy A

• rozwiązanie układu równań z użyciem LU

int g s l _ l i n a l g _ L U _ s o l v e ( c o n s t g s l _ m a t r i x * LU , c o n s t g s l _ p e r m u t a t i o n * p , c o n s t g s l _ v e c t o r * b , g s l _ v e c t o r * x )

• mnożenie macierz-wektor (BLAS 2)

(5)

y = α op(A) · ⃗x + β ⃗y op(A) = A, A

^T

, A

^H

int g s l _ b l a s _ d g e m v ( C B L A S _ T R A N S P O S E _ t TransA , d o u b l e alpha , c o n s t g s l _ m a t r i x * A , c o n s t g s l _ v e c t o r * x , d o u b l e beta , g s l _ v e c t o r * y )

( T r a n s A = C b l a s N o T r a n s , C b l a s T r a n s , C b l a s C o n j T r a n s )

• dodawanie dwóch wektorów (BLAS 1)

⃗

y = α⃗ x + ⃗ y

int g s l _ b l a s _ d a x p y ( d o u b l e alpha , c o n s t g s l _ v e c t o r * x , g s l _ v e c t o r * y )

5 Przykładowe wyniki

T(x,y)

’100.dat’ u 1:2:3

0 5 10 15 20 25 30 35 40 x

0 5 10 15 20 25 30 35 40

y

0 5 10 15 20 25 30 35 40

∇² T(x,y)

’100.dat’ u 1:2:4

0 5 10 15 20 25 30 35 40 x

0 5 10 15 20 25 30 35 40

y

0 0.02 0.04 0.06 0.08 0.1 0.12

Rysunek 2: Wyniki dla it = 100

Projekt 9: Dyfuzja ciepła - metoda Cranck-Nicloson.