KOMÓRCE” 29
• Pod-krok 2: Dyfuzja
dωp
dt = ν∆ω. (3.98)
3.8. Algorytm numeryczny metody wirowej typu „Wir w
komórce”
Poniżej zostanie przedstawiony algorytm numeryczny, wedle którego realizowane były obliczenia metodą typu ”Wir w komórce”
Obszar obliczeniowy został pokryty siatką obliczeniową o równomiernym kro-ku h w każdym kierunkro-ku. Na każdym węźle siatki inicjowana była cząstka wirowa o intensywności
αp =
Z
dxdydy
ω(x, y, z, t0)dxdydz ≈ ω(xn, yn, zn, t0) · h3, (3.99)
gdzie xn, yn, zn były współrzędnymi węzła siatki.
W niniejszej pracy początkowe lokalizacje cząstek wybierane były na węzłach siatki kartezjańskiej. Alternatywą może być inicjalizacja pseudo-losowa [17] w której cząstki są inicjalizowane w losowych lokalizacjach wewnątrz komórek siatki.
W każdym kroku czasowym wykonywane są następujące obliczenia:
1. Wyznaczany był potencjał wektorowy A rozwiązując równanie Poissona (3.21) dla trzech składowych wirowości.
∆Ai = −ωi, i = 1, 2, 3. (3.100)
Na podstawie potencjału wektorowego wyznaczane było pole prędkości ze wzoru
u = ∇ × A. (3.101)
2. Cząstki z węzłów siatki przemieszcza się zgodnie z lokalnym polem prędko-ści. Do rozwiązania równania (3.96) wykorzystana została metoda Rungego-Kutty czwartego rzędu. Prędkość z węzłów na aktualne położenie cząstki była interpolowana za pomocą jądra M40 (3.80).
3. Zmiana intensywności cząstek (na skutek rozciągania linii wirowych) by-ła realizowana zgodnie ze wzorem (3.97) również metodą Rungego-Kutty czwartego rzędu. Wartości pochodnych prędkości pod operatorem gradien-tu były obliczane za pomocą metody różnic skończonych drugiego rzędu. Prędkość z węzłów siatki na aktualne położenia cząstek była interpolowana
za pomocą jądra M40 (3.80). W tym wypadku lepszym sposobem obliczania prawej strony równania (3.97) jest użycie transpozycji gradientu prędkości
dαp
dt = [∇u(xp, t)]
T · αp. (3.102)
Schemat ten dawał najlepsze wyniki jeżeli chodzi o zachowanie energii ki-netycznej. W opracowaniach [17, 81] autorzy wskazują, że obydwie posta-ci członów źródłowych dają te same wyniki jeżeli spełniony jest warunek ∇ · ω ≡ 0. Niestety obliczenia numeryczne wprowadzają pewną niewielką niezerową źródłowość pola wirowości. W tym wypadku w pracach [81, 85] postulowane jest, że najlepszy jest schemat transponowany.
4. Po zakończonym kroku przemieszczania cząstek ich intensywność była re-dystrybuowana z powrotem na węzły siatki. W niniejszej pracy było to wy-konywane w każdym kroku czasowym a nie po arbitralnie przyjętej liczbie kroków czasowych. Zastosowanie takiego rozwiązania pozwoliło na rozwią-zywanie równania dyfuzji na siatce, a co za tym idzie wykorzystanie dokład-nych schematów numeryczdokład-nych dopasowadokład-nych do tego równania. Dodatko-wo pozDodatko-woliło to na zredukowanie algorytmu o jedną operację interpolacji (nowej intensywności na stare położenia cząstek).
5. Ostatnim krokiem było rozwiązanie równania dyfuzji na siatce numerycznej. Równanie dyfuzji w niniejszej pracy było rozwiązywane za pomocą metody różnic skończonych. To podejście częściowo ogranicza zyski związane z wy-korzystaniem metod Lagrange’owskich, ponieważ wprowadza do obliczeń siatkę numeryczną oraz dyfuzję numeryczną. Dzięki temu jednak możliwe jest dokładne rozwiązanie równania dyfuzji oraz przyspieszenie obliczeń.
Rozdział 4
Badania numeryczne metodą
typu Wir-w-komórce na
pojedynczym procesorze
Przedstawiony w poprzednim rozdziale algorytm realizujący obliczenia metodą cząstek wirowych typu „Wir w komórce” został zaimplementowany w języku programowania Fortran 90.
Obszar obliczeniowy został opisany strukturalną siatką numeryczną. Na każ-dym węźle siatki została umieszczona cząstka wirowa. Początkowa intensywność cząstki była obliczana zgodnie ze wzorem (3.18). Następnie rozwiązywany był układ równań ruchu cząstki (3.96), (3.97). Prędkość na węzłach siatki była wy-znaczona z równania Poissona (3.21). Do jego rozwiązania została użyta biblioteka Fishpack [95,96], która wykorzystuje metodę szybkiego rozwiązywania równania Poissona (ang. Fast Poisson Solver - FPS) drugiego rzędu. Równania ruchu zo-stały rozwiązane metodą Rungego-Kutty czwartego rzędu. Wartości prędkości cząstki były interpolowane z węzłów siatki za pomocą jądra (3.80). Po każdym kroku czasowym wykonywana była redystrybucja wirowości na węzły siatki przy pomocy jądra interpolacyjnego (3.80). To kończyło krok czasowy przy symula-cji przepływu nielepkiego. Przedstawione poniżej wyniki odnoszą się do płynu nielepkiego.
4.1. Badanie numeryczne prędkości przemieszczania się
pier-ścienia wirowego
Aby sprawdzić poprawność implementacji metody został wykonany test, w któ-rym badana była prędkość przemieszczania się pierścienia wirowego w czasie.
Rysunek 4.1: Geometria pierścienia wirowego
Pierścień wirowy wraz z oznaczeniami przedstawiony jest na rys.4.1. Dla pier-ścienia wirowego o cienkim rdzeniu (r0/R0 << 1) i stałym rozkładzie wirowości wewnątrz rdzenia prędkość translacji można wyznaczyć ze wzoru Kelvina [29]:
U = Γ 4πR0 ln 8R0 r0 −1 4+ O r0 R0 . (4.1)
Istnieje także wzór na prędkość pierścienia wirowego z nieruchomym płynem wewnątrz (wirowość znajduje się jedynie na obwodzie rdzenia) nazywany wzorem Hicksa [88]: U = Γ 4πR0 ln 8R0 r0 −1 2+ O r0 R0 . (4.2)
W celu porównania wyników symulacji z wartościami teoretycznymi oraz sprawdzenia poprawności kodu, przeprowadzone zostały badania numeryczne, w których promień pierścienia i promień w chwili początkowej wynosiły R0 = 1.5, r0 = 0.3, a zmieniano początkową wartość cyrkulacji w przedziale Γ ∈ [1.06, 4.23]. Obszar obliczeniowy miał wymiary 2π ×2π ×2π i był podzielony siatką numerycz-ną o 128 × 128 × 128 węzłach. We wszystkich kierunkach były zadane okresowe warunki brzegowe. Warunek początkowy był realizowany poprzez nadanie odpo-wiedniej wartości wektora intensywności wszystkim węzłom siatki numerycznej, które spełniały równanie
rp=
q
x2+ y2− R0
2
+ z2< r02, (4.3) czyli leżały wewnątrz pierścienia wirowego.
Równania (3.96) i (3.97) były rozwiązywane w każdym kroku czasowym me-todą Rungego-Kutty czwartego rzędu z krokiem czasowym 4t = 0.01. Wyniki symulacji widoczne są na rys. 4.2. Przemieszczenie pierścienia wirowego dla naj-większej wartości cyrkulacji (Γ = 4.23) przedstawiono na rys.4.3 i rys.4.4.
Na rys.4.2widać dobrą zbieżność pomiędzy wynikami teoretycznymi i nume-rycznymi. Pewna rozbieżność pomiędzy wynikami numerycznymi a teoretyczny-mi może wynikać z faktu, że formuła (4.1) została wyprowadzona dla pierścienia