Projekt 7: Równania Naviera-Stokesa - symulacja przepływu cieczy lepkiej.

Projekt 7: Równania Naviera-Stokesa - symulacja przepływu cieczy lepkiej.

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

1 Wstęp

A B

C

D E

F

x y

W ej cie W yj cie

Rysunek 1: Geometria układu przez który przepływa ciecz lepka. Z lewej strony Wejście (ścianka A), z prawej Wyjście (ścianka C). Szarym kolorem zaznaczono nieprzepuszczalne ścianki, które stanowią brzeg (ścianki: B, D, E, F)

Na zajęciach naszym zadaniem będzie rozwiązanie układu równań Naviera-Stokesa

∇ ² ψ = ζ (1)

µ ∇ ² ζ = ρ ( ∂ψ

∂y

∂ζ

∂x − ∂ψ

∂x

∂ζ

∂y )

(2) dla układu pokazanego na rysunku (ρ - gęstość, µ - lepkość, ζ - funkcja wirowości, ψ - funkcja stru- mienia).

1.1 Dyskretyzacja

Wprowadzamy siatkę równoodległych węzłów

∆x = ∆y = ∆ (3)

x = x i = ∆ · i, i = 0, 1, 2, . . . , n x (4) y = y _j = ∆ · j, j = 0, 1, 2, . . . , n y (5)

ψ(x, y) = ψ(x _i , y _j ) = ψ _i,j (6)

ζ(x, y) = ζ(x _i , y _j ) = ζ _i,j (7)

1

W równaniach (1) i (2) pochodne zastępujemy ilorazami różnicowymi centralnymi [df /dx = (f _i+1 − f i −1 )/(2∆), d ² f /dx ² = (f i+1 −2f i +f i −1 )/∆ ² ] i przegrupowujemy wyrazy tak aby ψ i,j oraz ζ i,j wyrazić za pomocą pozostałych wyrazów

ψ i,j = 1 4

(

ψ i+1,j + ψ i −1,j + ψ i,j+1 + ψ i,j −1 − ∆ ² ζ i,j

)

(8) ζ i,j = 1

4 (ζ i+1,j + ζ i −1,j + ζ i,j+1 + ζ i,j −1 )

− Ω ρ

16µ [(ψ i,j+1 − ψ i,j −1 )(ζ i+1,j − ζ i −1,j ) − (ψ i+1,j − ψ i −1,j )(ζ i,j+1 − ζ i,j −1 )] (9) Układ równań (8) i (9) rozwiążemy metodą relaksacji Gaussa-Seidla. Uwaga: parametr Ω wprowadza- my do równania (9), aby zapewnić stabilność relaksacji.

1.2 Warunki brzegowe

Zakładamy że prędkość cieczy [⃗ V = (u, v)] na wejściu/wyjściu ma tylko jedną niezerową składową V ⃗ _we/wy = (u _we/wy , 0). Korzystając z wykładu możemy zapisać

u we = Q we

2µ (y − y j

)(y − y n

) (10)

u _wy = Q _wy

2µ y(y − y n

) (11)

przy czym gradienty ciśnienia Q _we i Q _wy muszą się różnić, bo przekroje na we/wy są różne. Korzystając z zasady zachowania masy (strumień wpływający=strumień wypływający ) ^∫ _we u _we dy = ^∫ _wy u _wy dy dostajemy warunek

Q _wy = Q _we y _n ³

− y j ³

− 3y j

y ² _n

+ 3y _j ²

y _n

y ³ _n

(12)

1.2.1 WB dla ψ

• brzeg A (wejście) ψ _0,j = Q _we

2µ [ y ³

3 − y ²

2 (y _j

+ y _n

) + y · y j

· y n

]

, i = 0, j = j ₁ , . . . , n _y (13)

• brzeg C (wyjście) ψ _n

_,j = Q _wy

2µ [ y ³

3 − y ² 2 y _n

]

+ Q _we · y _j ²

( −y j

+ 3y _n

)

12µ , i = n _x , j = 0, 1, . . . , n _y (14) Uwaga: drugi wyraz w równaniu (14) zapewnia sklejanie się WB w węzłach (n x , 0) i (n x , n y ).

Otrzymujemy go z warunku ψ _we (0, y _j

) = ψ _wy (x _n

, 0).

• brzeg B

ψ _i,n

= ψ _0,n

, i = 1, 2, . . . , n _x − 1, j = n y (15)

• brzeg D

ψ i,0 = ψ 0,j

, i = i 1 , . . . , n x − 1, j = 0 (16)

• brzeg E

ψ _i

_,j = ψ _0,j

, i = i ₁ , j = 1, . . . , j ₁ (17)

• brzeg F

ψ i,j

= ψ 0,j

, i = 1, . . . , i 1 , j = j 1 (18)

2

1.2.2 WB dla ζ

• brzeg A (wejście)

ζ 0,j = Q we

2µ (2y − y j

− y n

), i = 0, j = j 1 , . . . , n y (19)

• brzeg C (wyjście)

ζ _n

_,j = Q _wy

2µ (2y − y n

), i = n _x , j = 0, . . . , n _y (20)

• brzeg B

ζ _i,n

= 2

∆ ² [ψ _i,n

₋₁ − ψ i,n

], i = 1, . . . , n _x − 1 (21)

• brzeg D

ζ _i,0 = 2

∆ ² [ψ _i,1 − ψ i,0 ], i = i ₁ + 1, . . . n _x−1 (22)

• brzeg E

ζ i

,j = 2

∆ ² [ψ i

+1,j − ψ i

,j ], j = 1, . . . j 1 − 1 (23)

• brzeg F

ζ _i,j

= 2

∆ ² [ψ _i,j

₊₁ − ψ i,j

], i = 1, . . . i ₁ (24)

• wierzchołek E/F

ζ _i

_,j

= 1

2 [ζ _i

_−1,j

+ ζ _i

_,j

₋₁ ] (25)

2 Algorytm relaksacji równań NS

Zgodnie z równaniem (1) możemy zapisać:

∇ ² ψ − ζ = δ (26)

Wartość parametru δ powinna maleć w trakcie relaksacji teoretycznie do zera, co możemy wykorzystać do kontroli zbieżności, np. całkując błąd w wybranych węzłach

Γ =

n ∑

−1

i=1

[ψ i+1,j

+ ψ i −1,j

+ ψ i,j

+1 + ψ i,j

−1 − 4ψ i,j

− ∆ ² ζ i,j

], j 2 = j 1 + 2 (27)

Relaksację można przeprowadzić opierając się na poniższym pseudokodzie (z kontrolą błędu Γ - wy- świetlaną na ekranie)

u s t a l a m y WB : ψ

FOR IT =1 TO I T _ M A X S T E P 1 DO IF ( IT < 2 0 0 0 ) T H E N

Ω = 0 EL S E

Ω = 1 END IF

FOR i =1 TO n _X − 1 STEP 1 DO FOR j =1 TO n Y − 1 STEP 1 DO

IF ( (i, j) ̸= BRZEG ) THEN ψ i,j = wzór (8)

3

ζ _i,j = wzór (9) END IF

END DO END DO

m o d y f i k a c j a WB : ζ k o n t r o l a b ł ę d u : Γ END DO

3 Zadania do wykonania

1. W obliczeniach należy użyć wartości parametrów: ∆ = 0.01, ρ = 1, µ = 1, n x = 200, n y = 90, i ₁ = 50, j ₁ = 55, IT M AX = 20000

2. Napisać dwie funkcje, w których zmieniane będą WB dla ψ i ζ.

3. Zaimplementować algorytm relaksacji równań NS.

Uwaga: dla it < 2000, zakładamy że Ω = 0 co oznacza, że wyłączamy chwilowo drugi wyraz w rów. (9), bo może on powodować niestabilność. Dla it > 2000, gdy ψ i ζ w pewnym stopniu się ustabilizują, włączamy ten wyraz (Ω = 1) do obliczeń i kontynuujemy relaksację.

4. Wykonać relaksację równań NS dla Q = −1000. Po jej zakończeniu sporządzić: wykres konturowy ψ, wykres konturowy ζ, mapę rozkładu składowej poziomej prędkości u(x, y) = ∂ψ/∂y, mapę rozkładu składowej pionowej prędkości v(x, y) = −∂ψ/∂x. ( 60 pkt)

5. Wykonać relaksację równań NS dla Q = −4000 i sporządzić rysunki jak w poprzednim punkcie.

Dobrać tak ilość konturów aby na wykresie ψ był wyraźnie widoczny wir. (30 pkt)

6. Wykonać relaksację równań NS dla Q = +4000 (odwracamy przepływ: ciecz płynie w lewo).

Sporządzić wykres konturowy ψ - czy wir utworzy się przed przeszkodą? (10 pkt)

4