Ćwiczenie 3a: Przepływ potencjalny
∗11 kwietnia 2016
Nieściśliwa i nielepka ciecz płynie przez rurę z przeszkodą w kształcie od- wróconej litery “L”.
1 50 100
1 50 100 150 200
y
x A B
C D
Rysunek 1: Przekrój rury z przewężeniem - schemat pudła obliczeniowego do opi- su cieczy przepływającej przez tę rurę. Rozpatrywany fragment znajduje się w obszarze [xmin, xmax]× [ymin, ymax] = [1, 200]× [1, 100]. Współrzędne (x, y) na- rożników przeszkody: A = (85, 85), B = (100, 85), C = (100, 70), D = (115, 70).
Interesuje nas rozkład prędkości. Rozwiązania poszukujemy na dwuwymia- rowej siatce o wymiarach 200× 100 punktów. Przyjmujemy krok przestrzenny
∆x = ∆y = 1.0.
Ze względu na nielepkość cieczy przepływ jest potencjalny (bezwirowy), tzn.
istnieje funkcja ϕ(x, y) (nazywana potencjałem przepływu) taka, że wektor pręd- kości cieczy (u, v) dany jest przez
u = ∂ϕ(x, y)
∂x , v = ∂ϕ(x, y)
∂y (1)
∗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
(u to prędkość w kierunku poziomym, v - w pionowym). Potencjał przepływu spełnia równanie Laplace’a
∇2ϕ(x, y) = 0. (2)
Zadanie 1. Linie strumienia (50 punktów)
Problem przepływu potencjalnego wygodnie rozwiązać używając funkcji stru- mienia ψ(x, y). Funkcja ta również spełnia równanie Laplace’a
∇2ψ(x, y) = 0. (3)
Rozwiązujemy równanie (3). Zastosujemy metodę różnic skończonych. Równanie Laplace’a dyskretyzujemy przy pomocy relaksacji (wzór z laboratorium nr 2):
ψi,j:= ψ(i−1),j+ ψi,(j−1)+ ψ(i+1),j+ ψi,(j+1)
4 . (4)
Iterację prowadzimy tylko na punktach spoza brzegu (omijamy brzegi siatki, brzegi przeszkody, oraz - wyłączone z siatki obliczeniowej - wnętrze przeszkody).
Zbieżność (z tolerancją np. tol = 10−6) badamy na podstawie całki działania:
a = 1
2
∫ xmax
xmin
∫ ymax
ymin
[(∂ψ(x, y)
∂x )2
+
(∂ψ(x, y)
∂y )2]
dxdy ≈
≈ 1
2
∑
i
∑
j
[(ψ(i+1),j− ψ(i−1),j
2∆x
)2
+
(ψi,(j+1)− ψi,(j−1)
2∆y
)2]
∆x∆y =
= 1
8
∑
i
∑
j
[(ψ(i+1),j− ψ(i−1),j)2
+(
ψi,(j+1)− ψi,(j−1))2]
. (5)
Warunki brzegowe (pozostałe punkty uzupełniamy zerami):
• na lewym i prawym brzegu siatki dajemy funkcję strumienia taką, jak dla przepływu swobodnego ψ0(x, y) = u0y - odpowiada to rozkładowi prędkości cieczy (u0, 0); przyjąć u0= 1.
• na całym górnym brzegu (górna krawędź rury i brzegi przeszkody) usta- lamy warunek brzegowy ψ(x, y) = ψ0(x = xmin, y = ymax). Dzięki temu górny brzeg będzie linią strumienia: ψ(x, y) = const. Prędkości cieczy są równoległe do linii strumienia, co daje nam odpowiednie warunki brzego- we na prędkość cieczy: znikanie v na osi oraz składowych prędkości cieczy normalnych do przeszkody.
• analogicznie postępujemy z dolną krawędzią rury. W tych punktach siatki ustalamy ψ(x, y) = ψ0(xmin, ymin).
Wyniki do uzyskania: narysować linie strumienia cieczy ψ(x, y) (w Gnu- plocie: konturami) oraz całkę działania dla kolejnych iteracji w skali logaryt- micznej na osi OX.
2
Zadanie 2. Linie potencjału (50 punktów)
Rozwiązać równanie (2). Stosujemy metodę relaksacji według dyskretyzacji (4) z badaniem zbieżności przy pomocy całki (5).
Warunki brzegowe:
• daleko od przeszkody ciecz nie odczuwa jej obecności (u = u0, v = 0) i potencjał dany jest przez ϕ(x, y) = ϕ0(x, y) = u0x. Potencjał przepływu swobodnego przyjąć na lewym i prawym brzegu siatki.
• metoda w skończonej liczbie iteracji powinna poradzić sobie z każdymi warunkami startowymi; jednak kładąc początkowo potencjał przepływu swobodnego (ϕ(x, y) = ϕ0(x, y)) na całej siatce obliczeniowej możemy przyspieszyć działanie metody dzięki wygodnemu warunkowi początkowe- mu.
Na dolnym i górnym brzegu oraz na krawędziach przewężenia zastosujemy warunki typu Neumanna, które trzeba ustalać na początku każdej iteracji.
Ciecz nie wnika w krawędzie - w związku z tym znika składowa normalna pręd- kości do brzegów rury i przeszkody (czyli pochodna ϕ po x na odcinkach piono- wych oraz po y na odcinkach poziomych).
• Przed każdą następną iteracją należy przepisać ϕi,jmax = ϕi,(jmax−1) dla i∈ [imin, iA−1]∪[iD+1, imax] oraz ϕi,jmin= ϕi,(jmin+1)dla i∈ [imin, imax] (górna i dolna krawędź rury).
• Na krawędziach przeszkody:
a) ϕiA,j = ϕ(iA−1),j dla j∈ (jA, jmax];
b) ϕiC,j = ϕ(iC−1),j dla j∈ (jC, jB);
c) ϕiD,j = ϕ(iD+1),j dla j∈ (jD, jmax];
d) ϕi,jA= ϕi,(jA−1) dla i∈ (iA, iB);
e) ϕi,jC = ϕi,(jC−1) dla i∈ (iC, iD).
• W narożnikach przeszkody - w punktach A, C, D z rys. 1 rozsądnie jest zastosować średnie arytmetyczne wartości potencjału z sąsiednich dwóch punktów, należących do wnętrza obszaru całkowania, tj:
– punkt A: ϕiA,jA= ϕ(iA−1),jA+ϕ2 iA,(jA−1), – punkt C: ϕiC,jC = ϕ(iC −1),jC+ϕiC ,(jC −1)
2 ,
– punkt D: ϕiD,jD =ϕ(iD +1),jD+ϕiD ,(jD −1)
2 .
• W narożniku B można przyjąć np. ϕ(iB,jB) = ϕ(iB−1),(jB−1) (ten punkt i tak nie bierze udziału w obliczeniach schematu relaksacji, więc nie ma wpływu na wyniki).
Wyniki do uzyskania: narysować linie stałego potencjału. Powinny wyjść lokalnie prostopadłe do linii strumienia. Drugi rysunek: wykres całki działania dla kolejnych iteracji, skala logarytmiczna.
3
