Projekt 8: Równanie adwekcji-dyfuzji – symulacja transportu masy metodą Cranka-Nicolson.

(1)

Projekt 8: Równanie adwekcji-dyfuzji – symulacja transportu masy metodą Cranka-Nicolson.

Tomasz Chwiej,Alina Mreńca-Kolasińska, Elżbieta Strzałka 11 grudnia 2018

1 Wstęp

x y

Rysunek 1: Geometria układu, w którym rozkład gęstości u(x, y, t) zmienia się ze względu na adwekcję w polu prędkości i pod wpływem dyfuzji. Na lewym i prawym brzegu nałożone są periodyczne warunki brzegowe – płyn przepływający przez prawy brzeg pojawia się na lewym brzegu.

Na zajęciach wykonamy symulację transportu masy rozwiązując równanie adwekcji-dyfuzji

∂u

∂t = −⃗v · ∇u + D∇

²

u (1)

wykorzystując niejawny schemat Cranka-Nicolson.

2 Dyskretyzacja, metoda CN, periodyczne warunki brzegowe, warunki początkowe

2.1 Dyskretyzacja

Wprowadzamy siatkę równoodległych węzłów w przestrzeni (x,y) oraz na osi czasu (t)

t = t

n

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

∆x = ∆y = ∆ (3)

x = x

_i

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

x

(4)

y = y

_i

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

y

(5)

u(x, y, t) = u(x

i

, y

j

, t

n

) = u

ⁿ_i,j

(6)

⃗ v = (v

^x

(x

i

, y

j

), v

^y

(x

i

, y

j

)) = (v

^x_i,j

, v

_i,j^y

) (7)

(2)

2.2 Metoda Cranka-Nicolson

Równanie (1) dyskretyzujemy stosując dwupunktowy iloraz różnicowy ”do przodu” dla pochodnej czasowej oraz symetryczne ilorazy różnicowe dla pochodnych przestrzennych

u

ⁿ⁺¹_i,j

− u

ⁿi,j

∆t = − 1

2 v

_i,j^x

[( u

ⁿ_i+1,j

− u

ⁿ_i_−1,j

2∆

) +

( u

ⁿ⁺¹_i+1,j

− u

ⁿ⁺¹_i_−1,j

2∆

)]

− 1

2 v

_i,j^y

[( u

ⁿ_i,j+1

− u

ⁿ_i,j₋₁

2∆

) +

( u

ⁿ⁺¹_i,j+1

− u

ⁿ⁺¹_i,j₋₁

2∆

)]

+ 1 2 D

[ u

ⁿ_i+1,j

+ u

ⁿ_i_−1,j

+ u

ⁿ_i,j+1

+ u

ⁿ_i,j₋₁

− 4u

ⁿ_i,j

∆

²

+ u

ⁿ⁺¹_i+1,j

+ u

ⁿ⁺¹_i_−1,j

+ u

ⁿ⁺¹_i,j+1

+ u

ⁿ⁺¹_i,j₋₁

−4u

ⁿ⁺¹_i,j

∆

²

]

(8) Równanie jest niejawne; aby je rozwiązać przenosimy u

ⁿ⁺¹_i,j

na lewą stronę i w każdej iteracji (przejście t

_n

→ t

n+1

) relaksujemy poniższe równanie 20 razy dla każdego węzła (i, j)

(

u

ⁿ⁺¹_i,j

⁾

^µ+1

=

( 1 1 +

^2D∆t_∆2

) {

u

ⁿ_i,j

− ∆t 2 v

^x_i,j

[( u

ⁿ_i+1,j

− u

ⁿ_i_−1,j

2∆

) +

( u

ⁿ⁺¹_i+1,j

− u

ⁿ⁺¹_i_−1,j

2∆

)

µ

]

− ∆t

2 v

^y_i,j

[( u

ⁿ_i,j+1

− u

ⁿ_i,j₋₁

2∆

) +

( u

ⁿ⁺¹_i,j+1

− u

ⁿ⁺¹_i,j₋₁

2∆

)

µ

]

+ ∆t 2 D

[ u

ⁿ_i+1,j

+ u

ⁿ_i_−1,j

+ u

ⁿ_i,j+1

+ u

ⁿ_i,j₋₁

− 4u

ⁿ_i,j

∆

²

+

( u

ⁿ⁺¹_i+1,j

+ u

ⁿ⁺¹_i_−1,j

+ u

ⁿ⁺¹_i,j+1

+ u

ⁿ⁺¹_i,j₋₁

∆

²

)

µ

]}

(9) gdzie µ = 0, 1, 2, . . . , 20 jest numerem iteracji Picarda. Na starcie przyjmujemy dla wszystkich węzłów

( u

ⁿ⁺¹_i,j

)

_µ=0

= u

ⁿ_i,j

(10)

2.3 Warunki brzegowe

Na lewy i prawy brzeg nakładamy periodyczne warunki brzegowe. Oznacza to, że

• lewym sąsiadem węzła (0, j) jest (n

x

, j)

• prawym sąsiadem węzła (n

x

, j) jest (0, j)

Dolny i górny brzeg oraz obrys zastawki omijamy (rów. adwekcji spełnia je automatycznie, bo ⃗ v

brzeg

= 0, ale rów. dyfuzji już nie i musimy o tym pamiętać).

2.4 Pole prędkości

Pole prędkości wygenerujemy mając do dyspozycji funkcję strumienia (wczytywaną z pliku). Poza brzegami i zastawką prędkości wyliczamy z definicji, v

^x

=

^∂ψ_∂y

oraz v

^y

= −

^∂ψ_∂x

, zastępując pochodne symetrycznymi ilorazami różnicowymi

v

^x_i,j

= ψ

_i,j+1

− ψ

i,j−1

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

_x

− 1, j = 1, 2, . . . , n

_y

− 1 (11) v

^y_i,j

= − ψ

i+1,j

− ψ

i−1,j

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

_x

− 1, j = 1, 2, . . . , n

_y

− 1 (12) Na zastawce ustalamy

x y

(3)

na dolnym i górnym brzegu

V

_i,0^x

= V

_i,n^y _y

= 0, i = 1, . . . , n

x

− 1 (14) a na lewym i prawym brzegu przepisujemy prędkości z sąsiednich węzłów

V

_0,j^x

= V

_1,j^x

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

_y

(15)

V

_n^x_x_,j

= V

_n^x_x_−1,j

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

_y

(16) 2.5 Warunek początkowy

Jako warunek początkowy [u(x, y, t = 0)] przyjmujemy poniższy rozkład gęstości u(x, y, t = 0) = 1

2πσ

²

exp (

− (x − x

A

)

²

+ (y − y

A

)

²

2σ

²

)

(17)

gdzie: (x

A

, y

A

) to położenie środka gęstości, a σ to rozmycie.

2.6 Algorytm dla równania AD Pseudokod

i n i c j a l i z a c j a p a r a m e t r ó w : n

x

, n

y

, i

1

, i

2

, j

1

, D , ∆ u t w ó r z t a b l i c e : u

⁰

, u

¹

,ψ ,v

^x

,v

^y

w c z y t a j f u n k c j ę s t r u m i e n i a : ψ w y z n a c z p o l e p r ę d k o ś c i : v

^x

, v

^y

w y z n a c z : v

max

, ∆t

i n i c j a l i z a c j a g ę s t o ś c i : u

⁰

= wzór (17) FOR IT =1 TO I T M A X S T E P 1 DO

// s t a r t i t e r a c j i P i c a r d a

i n i c j a l i z a c j a k o l e j n e g o k r o k u : u

¹

= u

⁰

FOR K =1 TO 20 S T E P 1 DO

FOR i =0 TO n

_x

ST E P 1 DO FOR j =1 TO n

y

-1 S T E P 1 DO

IF (( i , j ) ∈ ( zastawka )) THEN C O N T I N U E

EL S E IF ( i ==0 || i ==n

x

) T H E N u

¹_i,j

= wzór(9) + periodyczne WB EL S E

u

¹_i,j

= wzór(9) END IF

END DO END DO END DO

// z a c h o w u j e m y r o z w i ą z a n i e do n a s t ę p n e g o w y w o ł a n i a u

⁰

= u

¹

w y z n a c z i z a p i s z do p l i k u : c, x

_sr

END DO

(4)

3 Zadania do wykonania

1. W obliczeniach użyć parametrów: n

x

= 400, n

y

= 90, i

1

= 200, i

2

= 210, j

1

= 50, ∆ = 0.01, σ = 10 · ∆, x

A

= 0.45, y

_A

= 0.45.

2. Wczytać funkcję strumienia z pliku ”psi.dat”, format:

int i , int j , d o u b l e ψ

_i,j

3. Wyznaczyć pole prędkości. Stworzyć mapy prędkości v

^x

(x, y) oraz v

^y

(x, y). (20 pkt)

Uwaga: V

^x

powinno być wszędzie dodatnie (przepływ w prawą stronę), a V

^y

dodatnie przed zastawką i ujemne za nią.

Wyznaczyć krok czasowy ∆t według wzoru:

∆t = ∆

4v

_max

, (18)

gdzie v

max

to maksymalny moduł wektora prędkości (należy najpierw znaleźć punkt i, j, dla którego moduł prędkości |⃗v

i,j

| = ^√( v

^x_i,j

)

₂

+

( v

^y_i,j

)

₂

jest największy).

4. Zaimplementować algorytm dla równania adwekcji-dyfuzji.

5. Wykonać symulację dla D = 0 (brak dyfuzji). Sporządzić wykresy:

• całki gęstości c(t)

c(t

_n

) = ^∑

i,j

u

ⁿ_i,j

· ∆

²

, (19)

Uwaga: wartość c(t) może odchylać się od 1 o ±0.005.

• średniego położenia pakietu w kierunku ’x’

x

sr

(t

n

) = ^∑

i,j

x

i

· u

ⁿ_i,j

· ∆

²

, (20)

• 5 map rozkładu u(x, y, T

k

) dla T

_k

= k · t

max

/5, k = 1, . . . , 5, t

max

= IT M AX · ∆. War- tość t

_max

(lub inaczej: IT M AX) dobrać (empirycznie) tak aby na wykresie x

_sr

(t) widoczne były 3 maksima.

Za komplet wyników (50 pkt)

6. Powtórzyć obliczenia z poprzedniego punktu dla D = 0.1, stworzyć wykresy c(t), x

_sr

(t) oraz mapy u(x, y, T

K

), k = 1, . . . , 5. (30 pkt)

Uwaga: wartość c(t) powinna maleć w czasie (pakiet znika) ze względu na dyfuzję.

4 Przykładowe wyniki

(5)

0 0.5 1 1.5 2 2.5 3 3.5 4

0 0.1 0.2 0.3 0.4 0.5 0.6 xsr, c(t)

t xsr(D=0)

c(D=0) x_sr(D=0.1) c(D=0.1)

Rysunek 2: c(t) oraz x

_sr

(t) z dyfuzją oraz bez. Przejście od maksimum do minimum oznacza że pakiet przechodzi przez prawy brzeg (x = x

_max

) i pojawia się na lewym brzegu (x = 0).

’u1000.dat’ u 1:2:3

0 50 100 150 200 250 300 350 400 x

0 10 20 30 40 50 60 70 80 90

y

0 1 2 3 4 5 6 7

8 ’u2000.dat’ u 1:2:3

0 50 100 150 200 250 300 350 400 x

0 10 20 30 40 50 60 70 80 90

y

0 0.5 1 1.5 2 2.5 3 3.5

4 ’u3000.dat’ u 1:2:3

0 50 100 150 200 250 300 350 400 x

0 10 20 30 40 50 60 70 80 90

y

0 0.5 1 1.5 2 2.5 3

Rysunek 3: Rozkład gęstości w trzech wybranych chwilach czasowych dla D = 0.1. Na trzecim rysunku

widać, jak pakiet przechodzi przez prawy brzeg i pojawia się na lewym brzegu.