Ćwiczenie 3b: Przepływ stacjonarny cieczy nielepkiej w rurze z dwoma zagięciami.

(1)

Ćwiczenie 3b: Przepływ stacjonarny cieczy nielepkiej w rurze z dwoma zagięciami.

5 maja 2016

W ramach projektu należy znaleźć rozkład funkcji strumienia ψ(x, y) oraz po- tencjału ϕ(x, y) dla przepływu cieczy nielepkiej w układzie jak na rysunku 1.

Ciecz wpływa z lewej strony do rury, która zmienia następnie swój przekrój, a wypływa z prawej przez rurę o niewielkim przekroju. Kropki na rysunku ozna- czają brzeg, krzyżyki to obszar roboczy w którym poszukujemy rozwiązania.

Parametry wspólne dla obu zadań:

nx= 400, ny= 200, ja1 = 0, ja2 = 50, jb1= 150, jb2= 200, i1= 100, i2= 300,

∆x = ∆y = ∆ = 0.01, u0= 1.0.

x x x x x x x x x x x x x x

x x x x x x x x x x x x

WE

WY

0 1 2 i1 i2 n_x

j_a1=0 j_a2 n_y

j_b1 j_b2=n_y

C

D

Rysunek 1: Schemat siatki używanej przy poszukiwaniu ψ i ϕ.

1 Funkcja strumienia ψ(x, y)

Aby znaleźć funkcję strumienia, należy rozwiązać rówanie Poissona (a dokładniej Laplace’a), czyli∇²Ψ = 0. Rozwiązanie poszukujemy na siatce ψ(xi, yj) = ψi,j,

1

(2)

gdzie: xi= ∆· i, yi= ∆· j oraz i = 0, 1, . . . , nxi j = 0, 2, . . . , ny, z wyłączeniem punktów brzegowych (tam zadamy warunki Dirichleta) oraz obszarów bariero- wych (kolor szary na rysunku 1). Dla obszaru roboczego (krzyżyki na rysunku) stosujemy metodę nadrelaksacji:

ψ_i,j= (1− ω)ψi,j+ ωψ_i+1,j+ ψ_i_−1,j+ ψ_i,j+1+ ψ_i,j₋₁

4 (1)

W obliczeniach przyjąć parametr zbieżności ω = 1.9.

Warunki brzegowe Dirichleta ustalamy przed relaksacją:

1. Wejście: ψ_0,j= u₀· (yj− yja1) dla j = j_a1, . . . , j_a2 2. Wyjście: ψn_x,j = u0· (yj− yj_b1) dla j = jb1, . . . , jb2

3. Dolny brzeg (niebieska linia): ψi,j= ψ0,j_a1

4. Górny brzeg (czerwona linia): ψi,j= ψ0,j_a2

Zadania do wykonania:

1. Wyznaczyć rozkład funkcji strumienia ψ(x, y) dla obszaru zaznaczone- go krzyżykami na rysunku 1 stosując metodę nadrelaksacji (wzór 1). Na starcie, dla obszaru wewnętrznego (czyli poza brzegiem) przyjąć ψi,j= 0.

Wykonać wykres konturowy funkcji strumienia. (25 pkt)

Uwaga: Aby usprawnić obliczenia dla obszaru roboczego, którego kształt jest dość skomplikowany, najlepiej utworzyć dodatkową dwuwymiarową tablicę typu integer. Przed relaksacją do tablicy wpsujemy zera dla wszyst- kich węzłów brzegowych i tych położonych w barierach (kolor szary na rys.

1) i jedynki dla pozostałych węzłów (krzyżyki). W pętli iteracyjnej poten- cjał w węźle (i,j) modyfikujemy tylko wtedy, gdy w pomocniczej tablicy mamy wpisaną 1.

2. W trakcie relaksacji należy monitorować funkcjonał energii dla równania Poissona. Całkę liczymy wykorzystując do przybliżenia pierwszej pochodnej centralny dwupunktowy iloraz różnicowy:

S =

∫ ∫ dxdy1

2(∇ψ)²=∑

i

∑

j

∆² 2

[(ψi+1,j− ψi−1,j

2∆

)2

+

(ψi,j+1− ψi,j−1

2∆

)2]

(2) Uwaga: Całkę liczymy tylko w obszarze wewnętrznym (krzyżyki→ tablica pomocnicza). Proces relaksacji przerywamy jeśli spełniony jest warunek:

|^S^it+1_S_it^−S^it| < 10⁻¹⁰. Wykonać wykres zmian wartości funkcjonału S w funkcji numeru iteracji w podwójnej skali logarytmicznej. (25 pkt)

2

(3)

2 Potencjał przepływu ϕ(x, y)

Potencjał znajdujemy również z równania Poissona ∇²ϕ = 0, do rozwiązania którego również wykorzystamy metodę nadrelaksacji:

ϕi,j= (1− ω)ϕi,j+ ωϕ_i+1,j+ ϕ_i_−1,j+ ϕ_i,j+1+ ϕ_i,j₋₁

4 (3)

Dla parametru zbieżności przyjąć wartość ω = 1.9.

Warunki brzegowe: na wejściu i na wyjściu ustalamy warunki brzegowe typu Dirichleta, natomiast na pozostałych krawędziach (A,B,C,D,E,F) warunki typu von Nemanna (zerowanie się pochodnej normalnej do brzegu):

1. Wejście (Dirichlet): ϕ_0,j = u₀· x = 0 dla j = ja1, . . . , j_a2 2. Wyjście (Dirichlet): ϕn_x,j = u0· ∆ · nxdla j = jb1, . . . , jb2

3. Krawędź A (von Neuman):

∂ϕ/∂y = 0 =⇒ ϕi,j_a2 = ϕi,j_a2−1 dla i = 1, 2, . . . , i1− 1 4. Krawędź B (von Neuman):

∂ϕ/∂x = 0 =⇒ ϕi₁,j = ϕi₁+1,j dla j = ja2+ 1, . . . , ny− 1 5. Krawędź C (von Neuman):

∂ϕ/∂y = 0 =⇒ ϕi,n_y = ϕi,ny−1 dla i = i1+ 1, . . . , nx− 1 6. Krawędź D (von Neuman):

∂ϕ/∂y = 0 =⇒ ϕi,ja1 = ϕ_i,j_a1₊₁ dla i = 1, 2, . . . , i₂− 1 7. Krawędź E (von Neuman):

∂ϕ/∂x = 0 =⇒ ϕi2,j = ϕ_i₂_−1,j dla j = j_a1+ 1, . . . , j_b1− 1 8. Krawędź F (von Neuman):

∂ϕ/∂y = 0 =⇒ ϕi,j_b1 = ϕi,j_b1+1dla i = i2+ 1, . . . , nx− 1

9. Na połączeniu krawędzi A-B dajemy średnią z przyległych węzłów:

ϕi₁,j_a2 = (ϕi₁−1,ja2+ ϕi₁,j_a2+1)/2

10. Na połączeniu krawędzi E-F dajemy średnią z przyległych węzłów:

ϕi₂,j_b1 = (ϕi₂,j_b1−1+ ϕi₂+1,j_b1)/2

Uwaga: Warunki Dirichleta narzucamy tylko raz, przed relaksacją. Warunki typu von Neumanna musimy aktualizować po wykonaniu jednej iteracji, ponieważ zmienia się potencjał w obszarze obliczeniowym.

Zadania do wykonania:

1. Znaleźć rozkład potencjału ϕ(x, y) na siatce wykorzystując schemat re- laksacyjny (3). Na starcie w obszarze obliczeniowym (krzyżyki) przyjąć ϕ(x, y) = 0. Sporządzić wykres konturowy zrelaksowanego potencjału ϕ(x, y).

(25 pkt)

3

(4)

2. W trakcie relaksacji równania Poissona należy monitorować wartość funk- cjonału energii, który liczymy podobnie jak dla funkcji strumienia

S =

∫ ∫ dxdy1

2(∇ϕ)²=∑

i

∑

j

∆² 2

[(ϕi+1,j− ϕi−1,j

2∆

)2

+

(ϕi,j+1− ϕi,j−1

2∆

)2]

(4) Sumowanie we wzorze (4) odbywa się oczywiście tylko po węzłach z obszaru roboczego (krzyżyki =⇒ tablica pomocnicza). Proces relaksacji przery- wamy jeśli spełniony jest warunek:|^S^it+1_S_it^−S^it| < 10⁻¹⁰. Sporządzić wykres wartości funkcjonału S w funkcji numeru iteracji w podwójnej skali logarytmicznej. (25 pkt)

4