Wykonano symulacje przepływu powietrza przez komorę, aby zbadać wpływ modelu turbulencji na wyniki obliczeń numerycznych. Analizowano zwłaszcza rozkład linii prądu w komorze i profile składowej poziomej wektorów prędkości wzdłuż wybranych prostych przechodzących przez komorę. Wyniki numeryczne poddano walidacji z danymi pomiarowymi.
Geometria układu
Dwuwymiarową geometrię analizowanego układu, zaproponowaną przez Nielsena [4, 5], obecnie wykorzystuje się do testowania modeli/procedur numerycznych, m.in. modeli turbulencji zaimplementowanych w programach CFD. Geometrię stanowi przekrój przez pomieszczenie (zwane dalej komorą), z wlotem powietrza w lewym górnym rogu i wylotem w prawym dolnym rogu. Model geometryczny, dla którego wygenerowano siatkę i przeprowadzono symulacje, zawierał dodatkowo kanał wlotowy o długości 1 m i kanał wylotowy o długości 3 m. Schemat geometrii pokazano na rysunku obok. Linią przerywaną zaznaczono odcinki, dla których zostały wykreślone profile składowej poziomej i składowej pionowej prędkości. Komora symulacyjna ma proporcje geometryczne opisane w raporcie [4]: L/H = 3,0; h/H = 0,056; t/H = 0,16. Analizy numeryczne wykonano dla komory o długości L = 9 m, dlatego pozostałe wymiary były następujące: H = 3 m, h = 0,168 m, t = 0,48 m. Grubość przestrzeni symulacyjnej (różnica współrzędnych z obu ścian o warunku brzegowym typu symetria) była równa 0,1 m.
Siatka numeryczna
Zastosowano siatkę strukturalną. Długość i/lub wysokość komórek ma rozkład sinusoidalny. Przy ścianach komory oraz w obszarze wlotu i wylotu komórki mają najmniejsze rozmiary, a wraz ze wzrostem odległości od ścian przestrzeni symulacyjnej ich długość (obszar głównej komory) oraz wysokość (obszar głównej komory, wlotu i wylotu) rosną sinusoidalnie. Taki rozkład wielkości komórek wynika ze spodziewanych dużych gradientów analizowanych wielkości. W takich obszarach wymagana jest szczególnie duża dokładność obliczeń. Widok siatki w kierunku osi z przedstawiono na rysunku po prawej stronie.
W programie ANSYS ICEM zdefiniowano nazwy powierzchni ograniczających rozpatrywany obszar. Te nazwy zastosowano także w programie ANSYS CFX. Ogólna liczba węzłów numerycznych wynosiła 22 190, a liczba komórek – 10 820.
| Nazwa powierchni w CFX | Typ warunku brzegowego | Wartość | Jednostka |
| Inlet | Inlet, Normal Speed (u) | 0,455 | m/s |
| Inlet, Turbulence intensity | 5 | % | |
| Outlet | Outlet, Average Static Pressure | 0 | Pa |
| Wall_up | Wall, No Slip Wall, Smooth Wall | – | – |
| Wall_down | Wall, No Slip Wall, Smooth Wall | – | – |
| Symmetry_1 | Symmetry | – | – |
| Symmetry_2 | Symmetry | – | – |
Obliczenia numeryczne
Obliczenia numeryczne dotyczyły izotermicznego przepływu ustalonego powietrza o temperaturze 25°C, traktowanego jako gaz doskonały, o współczynniku lepkości dynamicznej 1,831∙10–5 Pa∙s. Działanie siły grawitacji pominięto (moduł wyporności był wyłączony). Ciśnienie referencyjne w komorze wynosiło 1 atm. Warunki brzegowe na poszczególnych powierzchniach symulowanej domeny zestawiono w tabeli.
Jako kryterium zakończenia obliczeń iteracyjnych przyjęto średniokwadratową wartość znormalizownych residuów równań bilansu masy i pędu na poziomie 1 ∙ 10–6. Warunek ten został spełniony przy najmniejszej liczbie iteracji (600) w symulacji z użyciem modelu turbulencji k–ε. W przypadku pozostałych modeli wymagana była większa liczba iteracji: dla modelu SST – 1112, dla modelu k–ω – 1661, a dla modelu SSG Reynolds Stress – 1909.
Wyniki obliczeń
Na poniższych rysunkach pokazano linie prądu oraz pola modułu prędkości dla czterech badanych modeli turbulencji.
Charakterystyka wirów
Na podstawie uzyskanych wyników można stwierdzić, że wybór modelu turbulencji ma istotny wpływ na przewidywanie zjawisk zachodzących podczas przepływu gazu przez komorę testową. Dowiedziono [5, 6], że w komorze wytworzą się trzy wiry: jeden duży o zwrocie wektorów prędkości zgodnym z ruchem wskazówek zegara (dalej określany mianem prawoskrętnego) oraz dwa znacznie mniejsze, mające zwrot wektorów prędkości przeciwny do ruchu wskazówek zegara (dalej nazywane lewoskrętnymi). Dwa mniejsze wiry to tzw. strefy recyrkulacji, zlokalizowane w lewym dolnym i w prawym górnym rogu komory. Wyniki dla standardowego modelu k–ε odzwierciedlają istnienie wszystkich wirów stwierdzonych empirycznie (nr 1, 2 i 3), w tym wiru głównego z centrum o przybliżonych współrzędnych (x, y) ≈ (2/3L + 1, 1/2H) [m]. Obliczenia z wykorzystaniem pozostałych modeli turbulencji sugerują, że w komorze powstają dodatkowe wiry.
Wszystkie zastosowane modele wykazują zgodność co do położenia centrum największego wiru – w przybliżeniu na wysokości 1,5 m i w odległości ok. 6,5 m od lewej ściany głównej komory. W przypadku modelu k–ε linie prądu głównego wiru (nr 1) wypełniają prawie całą komorę główną, z wyjątkiem obszaru o wysokości w zakresie h÷3h przy górnej ściance, obszaru o szerokości równej ok. t przy prawej ściance komory oraz trójkątnego obszaru w lewej dolnej części komory (wir nr 3). Na rysunku obszar bezwirowy oddzielono od obszarów wirowych za pomocą czarnych kreskowanych krzywych. Defektem modelu k–ε jest niedoszacowanie wielkości wirów nr 2 i 3.
W przypadku modelu k–ω (rys. 4) wyniki symulacji sugerują występowanie dodatkowego prawoskrętnego wiru (nr 4) w lewym dolnym rogu komory. Zajmuje on obszar o szerokości ok. 40 cm. Rozmiar wiru lewoskrętnego (nr 3, rys. 4), będącego odpowiednikiem strefy recyrkulacji nr 3 przewidzianej przez model k–ε, jest przeszacowany, ponieważ na linii y = 0,084 m rozciąga się on na przedział 0,4÷4,2 m. Podobnie przeszacowany jest rozmiar wiru nr 2. Można stwierdzić, że model k–ω źle oddaje charakter przepływu powietrza przez lewą dolną część komory testowej.
Model SST daje w wyniku cztery wiry; lewoskrętna strefa recyrkulacji nr 3 zajmuje całą lewą dolną część komory, a jej rozmiar jest przeszacowany – podobnie jak w przypadku użycia modelu k–ω. Wir nr 4 jest zlokalizowany nad wirem nr 3, a jego orientacja jest taka jak wiru głównego nr 1. Wielkość i kształt obszaru zajmowanego przez wir nr 2 są zbliżone do parametrów odpowiadającego mu wiru przewidzianego przez model k–ω.
Model SSG Reynolds daje wyniki sugerujące występowanie największej liczby wirów. Na rysunku widać ich pięć lub sześć. Linie prądu wiru nr 1 z dala od centrum nie mają kształtu eliptycznego, ponieważ opływają obszar zajęty przez wiry 3, 4 i 5. Wielkość obszaru zajmowanego przez wir nr 2 w prawym górnym rogu komory jest nieco mniejsza niż w przypadku modeli SST i k–ω, ale większa niż wynikająca z modelu k–ε.
Na ostatnim rzsunku widać, że obszar występowania dość szybkiego przepływu (wartość prędkości powyżej 0,1 m/s) jest ograniczony do prawej połowy komory, a prędkość w obszarze wirów 3, 4 i 5 jest bardzo mała. Ponadto można stwierdzić, że model SSG Reynolds przewiduje szybszy ruch wirowy cieczy wokół głównego wiru (nr 1) niż pozostałe modele turbulencji.
Składowe wektora prędkości
Uzyskane profile składowej poziomej u prędkości płynu podzielonej przez wartość prędkości wlotowej u0, pochodzące z symulacji wykorzystujących różne modele turbulencji, porównano z danymi eksperymentalnymi. Okazało się, że standardowy model turbulencji k–ε najlepiej opisuje przepływ gazu przez analizowaną komorę – profile prędkości wzdłuż wybranych linii są najbardziej zbliżone do punktów eksperymentalnych. Wyjątkami są tu dwie strefy: w prawym górnym (x od 9,4 do 10,0 m) oraz w lewym dolnym rogu komory (x od 1,15 do 2,1 m). W pierwszej z nich wraz ze wzrostem współrzędnej x składowa pozioma prędkości u maleje po wartościach dodatnich aż do zera, co jest niezgodne z danymi doświadczalnymi, które wskazują na ujemną wartość składowej u w przedziale 9,4÷10,0 m. W drugiej ze stref model k–ε sugeruje, że składowa u jest ujemna, czyli że jest to obszar głównego prawoskrętnego wiru, z doświadczenia wynika zaś, że występuje tam strefa lewoskrętnej recyrkulacji. Przyczyną tych niezgodności jest błąd modelu k–ε, polegający na niedoszacowaniu wielkości wirów nr 2 i 3. W prawym górnym rogu komory testowej najlepsze wyniki daje model SSG Reynolds, gdyż szerokość wiru Δx jest bardzo bliska wartości eksperymentalnej. Wadą rozwiązania numerycznego jest niedoszacowanie wartości prędkości w obszarze małego wiru. Stosunek u/u0 nie spada tam poniżej –0,04, natomiast z eksperymentu wynika, że minimum wynosi –0,20 dla x = 9,59 m.
Centrum głównego prawoskrętnego wiru jest zlokalizowane na prawo od linii wyznaczającej 2/3 długości komory. Wynika to z tego, że dla x = 7 m składowa v jest dodatnia dla y = 1,5 m, czyli płyn przemieszcza się w górę komory, opływając od lewej strony centrum wiru. Wartość bezwzględna składowej pionowej prękości nie przekracza ok. 0,093 m/s (model SSG Reynolds, linia x = 7 m, y = 1,3 m), podczas gdy wartości bezwzględne składowej poziomej osiągają ok. 0,50 m/s (modele SST i k–ω, linia y = 2,916 m, x = 1,0 m). Wynika to z połączenia dwóch faktów: z prawa zachowania masy oraz z tego, że pole przekroju komory w płaszczyźnie poziomej jest trzykrotnie większe od pola przekroju w płaszczyźnie pionowej.
Wnioski
Wybór modelu turbulencji stosowanego przez solvery w programach CFD ma wpływ na wyniki symulacji. Nieodpowiedni dobór tego modelu może prowadzić do błędnych wyników, wyraźnie odbiegających od rzeczywistego charakteru przepływu, który można poznać tylko na drodze eksperymentu. Prawdopodobną przyczyną występowania największej liczby wirów w symulacji z użyciem modelu SSG Reynolds jest to, że korzysta ona bezpośrednio z tensora naprężeń Reynoldsa bez uproszczeń. Umożliwia on opis wirów o małej intensywności, np. zanikających, a jednocześnie jest podatny na lokalne, chwilowe zmiany wartości jego składowych. To może prowadzić do obliczenia cyrkulujących pól prędkości tam, gdzie w rzeczywistości one nie występują. W przypadku symulowanej geometrii (prostopadłościennej komory) zastosowanie modelu SSG Reynolds powoduje nadprodukcję wirowości.
Najlepszą zgodność wyników symulacji z danymi eksperymentalnymi daje standardowy model k–ε.
Literatura
1. www.cfd-online.com/Wiki/Turbulence_modeling (dostęp: 07.01.2019 r.).
2. Jaworski Z. „Numeryczna mechanika płynów w inżynierii chemicznej i procesowej”. Warszawa: AOW EXIT, 2005.
3. Wilcox D.C. “Turbulence modeling for CFD”. Wyd. 2. Anaheim, Kalifornia: DCW Industries, Inc., 1998.
4. Lemaire A.D. (red.). “Room Air and Contaminant Flow, Evaluation of Computational Methods. Subtask-1 Summary Report”. Wyd. 2. Delft (Holandia), grudzień 1993.
5. Nielsen P.V. “Specification of a two-dimensional test case”. Listopad 1990, ISSN 0902–7513 R9040. 6. http://homes.civil.aau.dk/pvn/cfd-benchmarks/two_d_benchmark_test (dostęp: 07.02.2019 r.).