mgr inż. Piotr Kamil Korman1, dr hab. inż. Sławomir Pietrowicz, prof. PWr2
1 Politechnika Wrocławska, Wydział Chemiczny, Zakład Inżynierii Chemicznej
2 Politechnika Wrocławska, Wydział Mechaniczno-Energetyczny, Katedra Termodynamiki, Teorii Maszyn i Urządzeń Cieplnych
Walidacja modeli turbulencji na przykładzie komory testowej Annex 20 Spis treści
1. Wstęp teoretyczny ... 2
1.1. Pojęcie turbulencji ... 2
1.2. Modele turbulencji ... 2
1.3. Wielkości charakteryzujące turbulencje ... 3
1.4. Charakterystyka metody LES ... 4
1.5. Charakterystyka metod RANS ... 5
1.6. Modele turbulencji wykorzystane w symulacji ... 6
1.7. Modele k-epsilon i k-omega ... 6
1.8. Model SST ... 6
1.9. Model SSG Reynolds ... 6
2. Cel symulacji ... 7
3. Parametry modelu ... 7
3.1. Geometria układu ... 7
3.2. Parametry geometryczne ... 8
3.3. Siatka numeryczna ... 8
3.4. Warunki symulacji ... 9
3.5. Kryteria zakończenia obliczeń numerycznych ... 10
4. Wyniki obliczeń i analiza ... 11
4.1. Charakterystyka wirów ... 14
4.2. Składowe wektora prędkości wzdłuż wybranych profili ... 16
4.3. Lepkość turbulentna ... 19
5. Wnioski ... 20
6. Bibliografia ... 20
1 Spis symboli
Symbol Jednostka* Znaczenie
D m∙s-3 gęstość spektralna dyssypacji
E m∙s-2 = J∙kg-1∙m-1 gęstość spektralna energii kinetycznej fluktuacji f s-1 częstotliwość fluktuacji prędkości
H m wysokość komory
h m wysokość wlotu do komory
k m2∙s-2 = J∙kg-1 właściwa energia kinetyczna burzliwości
l m charakterystyczna wielkość wiru
p Pa ciśnienie
S s-1 tensor szybkości deformacji
S’ s-1 tensor szybkości deformacji fluktuacyjnych
T Pa tensor naprężeń lepkościowych
t s czas
t m wysokość wylotu z komory
T’ Pa tensor naprężeń Reynoldsa
T’ m2∙s-2 właściwy tensor naprężeń Reynoldsa
u m∙s-1 wektor prędkości
ui m∙s-1 składowa wektora prędkości w kierunku i
x, y, z m współrzędne położenia
Δ m średnia geometryczna długości boków komórki obliczeniowej ε m2∙s-3 szybkość dyssypacji energii kinetycznej fluktuacji
κ m-1 liczba falowa wiru
µ Pa∙s współczynnik lepkości dynamicznej ν m2∙s-1 współczynnik lepkości kinematycznej Π m2∙s-3 tensor naprężeń ciśnieniowych
ρ kg∙m-3 gęstość
ϕ wielkość skalarna
τ s czas
ω s-1 właściwa szybkość dyssypacji energii kinetycznej fluktuacji Indeksy górne
¯ część odfiltrowana (LES), średnia czasowa (RANS)
’ część fluktuacyjna
Indeksy dolne
f wielkość po operacji filtrowania
i, j, k indeksy kierunków układu współrzędnych T odniesienie do warunków turbulentnych η odniesienie do mikroskali Kołmogorowa
* Dla wektorów i tensorów podano jednostkę elementów wektora/tensora
2 1. WSTĘP TEORETYCZNY
1.1. Pojęcie turbulencji
Turbulencję można zdefiniować jako losowe fluktuacje wielkości charakteryzujących przepływ, m.in.
prędkości i ciśnienia. Fluktuacje te są zmienne w czasie oraz w przestrzeni. Empirycznie stwierdzono, że przepływ płynu w przewodzie zamkniętym zmienia się z laminarnego na turbulentny, gdy liczba Reynoldsa (Re) przekracza wartość ok. 2100 – 2300. Wizualnym świadectwem turbulencji jest występowanie wirów (ang. eddy lub vortex).
Matematycznie takie obszary cechuje niezerowa długość wektora rotacji prędkości zdefiniowanego jako rot(𝒖) = ∇ × 𝒖 = [𝜕𝑢𝑧
𝜕𝑦 −𝜕𝑢𝑦
𝜕𝑧 ,𝜕𝑢𝑥
𝜕𝑧 −𝜕𝑢𝑧
𝜕𝑥 ,𝜕𝑢𝑦
𝜕𝑥 −𝜕𝑢𝑥
𝜕𝑦]. (1) 1.2. Modele turbulencji
Symulowanie przepływów turbulentnych wymaga uwzględnienia w rozwiązaniu numerycznym niestabil- ności w przestrzeni i w czasie. Rozwój numerycznej mechaniki płynów w XX wieku zaowocował wieloma mo- delami przepływów turbulentnych, które stanowią uproszczenie w stosunku do najbardziej precyzyjnej, ale jedno- cześnie najbardziej czasochłonnej metody symulowania przepływów, jaką jest DNS. DNS (ang. Direct Numerical Simulation) to bezpośrednia symulacja numeryczna, która rozwiązuje bezpośrednio równania transportu, w tym równanie Naviera-Stokes’a. Równanie to dla cieczy nieściśliwej (div(u) = 0), jednorodnej (ρ = const) i znajdującej się poza polem grawitacyjnym ma postać [3]:
𝜌𝜕𝑢𝑖
𝜕𝑡 + 𝜌𝑢𝑗∙𝜕𝑢𝑖
𝜕𝑥𝑗 = − 𝜕𝑝
𝜕𝑥𝑖 +𝜕𝑇𝑗𝑖
𝜕𝑥𝑗, (2) gdzie T oznacza tensor naprężeń lepkościowych zdefiniowany jako 2µS, zaś S to symetryczny tensor szybkości deformacji:
𝑺 = [
𝑆11 𝑆12 𝑆13 𝑆21 𝑆22 𝑆23 𝑆31 𝑆32 𝑆33
], (3) gdzie
𝑆𝑖𝑗 =1 2(𝜕𝑢𝑖
𝜕𝑥𝑗 +𝜕𝑢𝑗
𝜕𝑥𝑖). (4) Wyniki w metodzie DNS są reprezentowane bez uproszczeń, ale metoda ta wymaga bardzo dużej mocy obliczeniowej. Aby przyspieszyć obliczenia, zaproponowano różne modele turbulencji oparte na określonych uproszczeniach. Metody te można podzielić na kilka grup przedstawionych poniżej [1, 2]:
• LES (ang. Large Eddy Simulation) – metoda dużych wirów, w której równania transportu są rozwiązywane bezpośrednio, ale algorytmy filtrujące uśredniają parametry przepływu w obszarach występowania małych wirów; przez mały wir rozumie się tutaj wir o wielkości rzędu rozmiaru komórki obliczeniowej lub mniejszej;
obliczane numerycznie wielkości, np. prędkość płynu, są rozdzielane na dwa człony, mianowicie odfiltrowaną część rozwiązywalną (resolvable-scale part) i część fluktuacyjną o skali mniejszej niż rozmiar siatki nume- rycznej (subgrid-scale part);
• RANS (ang. Reynolds Averaging of Navier-Stokes Equation Simulation) – metoda uśredniania Reynoldsa;
w tej grupie metod obliczane wielkości rozdziela się na dwie składowe: wartość uśrednioną po czasie 𝜙 i od- chylenie (fluktuację) 𝜙′; część metod z tej grupy wykorzystuje pojęcie lepkości turbulentnej µT (ang. turbulent viscosity lub eddy viscosity);
• DES (ang. Detached Eddy Simulation) – model będący połączeniem metod LES i RANS, który wykonuje obliczenia metodą RANS w warstwie przyściennej, a poza nią posługuje się metodą LES.
3 1.3. Wielkości charakteryzujące turbulencje
Turbulencje występujące podczas przepływu mają różne konsekwencje, m.in. zamiana energii kinetycznej wirów w energię cieplną na skutek działania sił lepkości. Rozmiary wirów wahają się od skali mikro do skali makro, dlatego, aby opisać tak rozległe ich spektrum, stosuje się funkcje gęstości spektralnej, które charakteryzują wybrane parametry wirów dla różniczkowych przedziałów zmiennej niezależnej. Są one analogiczne do tych uży- wanych przy opisie promieniowania elektromagnetycznego. Jedną z przyjmowanych zmiennych niezależnych [2]
jest częstotliwość fluktuacji prędkości f i obliczana na jej podstawie liczba falowa wiru κi, która określana jest wzorem
𝜅𝑖 =2𝜋𝑓
𝑢𝑖 , (5) gdzie 𝑢𝑖 to średnia prędkość w kierunku i. Charakterystyczna wielkość wiru w kierunku i (li) jest określana jako iloraz średniej prędkości i częstotliwości wiru. Wprowadza się także średnią liczbę falową κ, która jest długością wektora κ = [κ1, κ2, κ3].
Istotną funkcją, która umożliwia wyznaczenie wielkości wykorzystywanych przez część modeli turbulen- cji, w tym przez badany model k-epsilon, jest gęstość spektralna energii kinetycznej fluktuacji Ei [2]. Jeżeli zosta- nie ona scałkowana, to otrzyma się wartość wariancji i-tej składowej prędkości, mianowicie:
∫ 𝐸𝑖(𝜅𝑖)𝑑𝜅𝑖
∞ 0
= 𝑢𝑖′2. (6) Całkowanie trójwymiarowej gęstości spektralnej daje w wyniku:
∫ 𝐸(𝜅)𝑑𝜅
∞ 0
= ∫ 2𝜋𝜅2(𝐸1(𝜅1) + 𝐸2(𝜅2) + 𝐸3(𝜅3))𝑑𝜅
∞ 0
=1
2(𝑢1′ 2+ 𝑢2′ 2+ 𝑢3′ 2) = 𝑘, (7) gdzie k jest kinetyczną energią właściwą burzliwości, zapisywaną też w postaci wektorowej z konwencją suma- cyjną Einsteina, mianowicie:
𝑘 =1
2𝒖𝑖′𝒖𝑖′. (8) Drugą istotną wielkością uwzględnianą w modelach turbulencji jest ε, czyli szybkość dyssypacji (rozpra- szania) energii kinetycznej fluktuacji. Jest ona równa całce z gęstości spektralnej dyssypacji
∫ 𝐷(𝜅)𝑑𝜅
∞ 0
= ∫ 2𝜅2𝜇
𝜌𝐸(𝜅)𝑑𝜅
∞ 0
= 𝜀. (9) Wielkość ε można także wyprowadzić z równania transportu kinetycznej energii burzliwości k [3]. Człon tego równania opisujący rozpraszanie energii wskutek działania sił lepkości ma postać:
𝜀 = 2𝜇
𝜌∙ 𝑆𝑖𝑗′ ∙ 𝑆𝑖𝑗′ = 2𝜈 ∙1 2(𝜕𝑢𝑖′
𝜕𝑥𝑗+𝜕𝑢𝑗′
𝜕𝑥𝑖) ∙1 2(𝜕𝑢𝑖′
𝜕𝑥𝑗 +𝜕𝑢𝑗′
𝜕𝑥𝑖) = 𝜈 ∙𝜕𝑢𝑖′
𝜕𝑥𝑗 ∙ (𝜕𝑢𝑖′
𝜕𝑥𝑗 +𝜕𝑢𝑗′
𝜕𝑥𝑖), (10) gdzie Sij’ to składowa symetrycznego tensora szybkości deformacji fluktuacyjnych. Dyssypacja energii kinetycz- nej wirów polega na jej zamianie w inne formy energii, np. w energię cieplną.
Ścisły matematycznie opis przepływów burzliwych wykorzystujący równania transportu bez żadnych uproszczeń byłby bardzo czasochłonny, dlatego badacze, m.in. Rosjanin Andriej Kołmogorow, podejmowali próby stworzenia teorii ułatwiających opis przepływów turbulentnych. Kołmogorow wprowadził pojęcia mikro- skali długości, prędkości i czasu, odpowiednio η, uη i τη oraz stworzył teorię opartą na trzech hipotezach [2]:
• przepływy turbulentne są lokalnie izotropowe w sensie statystycznym dla odpowiednio małej skali i odpo- wiednio dużych liczb Reynoldsa,
• w przepływie burzliwym dla odpowiednio dużych liczb Reynoldsa i odpowiednio dużych wartości liczb falowych (odpowiadających tzw. obszarowi uniwersalnej równowagi) statystyczne rozkłady
4
prawdopodobieństwa różnych małych wirów są do siebie podobne; nie zależą one od warunków przepływu i geometrii układu, ale tylko od lepkości ν i od szybkości dyssypacji energii ε,
• w przepływie turbulentnym dla wystarczająco wysokich liczb Reynoldsa w tzw. bezwładnościowym pod- zakresie obszaru uniwersalnej równowagi statystyczne rozkłady prawdopodobieństwa różnych małych wi- rów są do siebie podobne i zależą tylko od szybkości dyssypacji energii ε.
Opis teorii Kołmogorowa i jej modyfikacji można znaleźć w literaturze z zakresu modelowania turbulencji, m.in. w książce Wilcoxa [3]. Jednym z parametrów występujących w modelach burzliwości jest turbulentna skala długości l, którą można zinterpretować jako rozmiar największych wirów występujących w symulowanym ukła- dzie. Ta wielkość jest zatem zależna od rozpatrywanej geometrii, w tym od średnicy hydraulicznej.
1.4. Charakterystyka metody LES
Metoda LES charakteryzuje się tym, że obliczane numerycznie wielkości, np. prędkość płynu, są rozdzie- lane na dwa człony, mianowicie odfiltrowaną część rozwiązywalną (resolvable-scale part) i część fluktuacyjną oskali mniejszej niż rozmiar siatki numerycznej (subgrid-scale part). Operacja rozdziału, czyli filtrowanie dotyczy zarówno wielkości skalarnych, jak i wektorowych:
𝜙(𝒙, 𝑡) = 𝜙𝑓(𝒙, 𝑡) + 𝜙𝑓′(𝒙, 𝑡), (11) 𝒖(𝒙, 𝑡) = 𝒖𝑓(𝒙, 𝑡) + 𝒖𝑓′(𝒙, 𝑡). (12) Do uśredniania stosuje się różne filtry, np. klatkowy i gaussowski [2]. Działanie opiera się na odpowiednim zastosowaniu operacji splotu dwóch funkcji – funkcji filtrującej G i funkcji wielkości badanej. Przykładowo dla filtra klatkowego operację filtrowania można zapisać jako
𝜙𝑓(𝒙, 𝑡) = ∭ 𝐺(𝒙 − 𝝃, Δ) ∙ 𝜙(𝝃, 𝑡)𝑑3𝝃
𝑉𝑘
, (13) gdzie Vk to objętość prostopadłościennej komórki obliczeniowej, a Δ to średnia geometryczna długości boków komórki. Funkcja filtrująca ma postać [3]:
𝐺(𝒙 − 𝝃, Δ) = {1/Δ3 dla |𝑥𝑖− 𝜉𝑖| < Δ𝑥𝑖/2,
0 dla |𝑥𝑖 − 𝜉𝑖| ≥ Δ𝑥𝑖/2. (14) Po wykonaniu filtrowania można zapisać równanie Naviera-Stokes’a w postaci [3]:
𝜕𝑢𝑓𝑖
𝜕𝑡 + 𝜕
𝜕𝑥𝑗(𝑢𝑓𝑖𝑢𝑓𝑗) = −1 𝜌∙ 𝜕𝑝
𝜕𝑥𝑖 +𝜇 𝜌∙ 𝜕
𝜕𝑥𝑘(𝜕𝑢𝑓𝑖
𝜕𝑥𝑘), (15) przy czym pochodną konwekcyjną wyraża się w dość skomplikowany sposób przy pomocy tensorów naprężeń, mianowicie:
𝑢𝑓𝑖𝑢𝑓𝑗 = 𝑢𝑓𝑖𝑢𝑓𝑗+ 𝐿𝑖𝑗 + 𝐶𝑖𝑗+ 𝑅𝑖𝑗, (16) gdzie:
𝐿𝑖𝑗 = 𝑢𝑓𝑖𝑢𝑓𝑗− 𝑢𝑓𝑖𝑢𝑓𝑗, (17) 𝐶𝑖𝑗 = 𝑢𝑓𝑖𝑢𝑓𝑗′ + 𝑢𝑓𝑗𝑢𝑓𝑖′ , (18) 𝑅𝑖𝑗 = 𝑢𝑓𝑖′ 𝑢𝑓𝑗′ . (19) Tensory L, C i R są określane mianem tensorów naprężeń Leonarda, naprężeń mieszanych i naprężeń Reynoldsa SGS (subgrid-scale). Wprowadza się także nowy tensor LES, oznaczany jako Q, który jest sumą tensorów
5
naprężeń mieszanych i naprężeń Reynoldsa. Równanie Naviera-Stokes’a zapisuje się także w innej postaci, uwzględniającej wspomniane tensory
𝜕𝑢𝑓𝑖
𝜕𝑡 + 𝜕
𝜕𝑥𝑗(𝑢𝑓𝑖𝑢𝑓𝑗) = −1 𝜌∙ 𝜕𝑃
𝜕𝑥𝑖 + 𝜕
𝜕𝑥𝑗(𝜇 𝜌∙𝜕𝑢𝑓𝑖
𝜕𝑥𝑗 + 𝜏𝑖𝑗), (20) przy czym
𝜏𝑖𝑗 = − (𝑄𝑖𝑗 −1
3𝑄𝑘𝑘𝛿𝑖𝑗), (21) 𝑃 = 𝑝 +1
3𝜌𝑄𝑘𝑘𝛿𝑖𝑗, (22) Zagadnieniem, które pozostaje do rozpatrzenia, jest obliczenie tensora naprężeń SGS, co dokonywane jest z użyciem modeli podsiatkowych opisanych w literaturze [1, 3], takich jak model Smagorinskiego. Korzystanie z metody LES jest uzasadnione tym, że tylko średnie i duże wiry są zależne od środowiska i geometrii, w której występują, zaś małe wiry (o skali podsiatkowej) są niezależne od geometrii [2].
1.5. Charakterystyka metod RANS
Ta grupa metod wykorzystuje uśrednianie wielkości skalarnych i wektorowych w dziedzinie czasu. Lo- kalne wartości, podobnie jak w metodzie LES, są w niej rozdzielane na dwa składniki, przy czym w przypadku RANS są to średnia czasowa i odchylenie. Dla procesów ustalonych można zapisać równania [2]:
𝜙(𝒙, 𝑡) = 𝜙(𝒙) + 𝜙′(𝒙, 𝑡), (23) 𝒖(𝒙, 𝑡) = 𝒖(𝒙) + 𝒖′(𝒙, 𝑡). (24) Istotnym faktem jest to, że średnia fluktuacja dowolnej wielkości liniowej jest równa 0. Równanie Naviera-Sto- kes’a dla przepływu turbulentnego z rozdziałem wielkości na średnie i odchylenia przyjmuje postać [3]:
𝜌𝜕𝑢𝑖
𝜕𝑡 + 𝜌𝑢𝑗∙𝜕𝑢𝑖
𝜕𝑥𝑗 = −𝜕𝑝
𝜕𝑥𝑖+ 𝜕
𝜕𝑥𝑗(2𝜇𝑆𝑗𝑖− 𝜌𝑢𝑗′𝑢𝑖′). (25) To właśnie to równanie jest określane mianem RANS (Reynolds-averaged Navier-Stokes equation). Występujący w nim składnik −𝜌𝑢𝑖′𝑢𝑗′= 𝜌𝑇𝑖𝑗′ to pomnożony przez gęstość element symetrycznego właściwego tensora naprężeń Reynoldsa T’. Trzy równania (25) dla poszczególnych składowych i równanie zasady zachowania masy (div(u) = 0) stanowią układ czterech równań z dziesięcioma niewiadomymi (sześć niewidomych to elementy ten- sora T’). Z tego powodu opis przepływów turbulentnych metodami typu RANS wymaga zastosowania modeli do obliczania elementów tensora naprężeń Reynoldsa. Gama tych modeli jest obszerna. W raporcie wykorzystano cztery z nich. W literaturze symbolem T’ oznacza się także tensor naprężeń Reynoldsa, który jest iloczynem gę- stości i tensora właściwego.
Równania różniczkowe pozwalające na obliczenie elementów tensora T’ uzyskuje się poprzez wykorzy- stanie własności momentów równania Naviera-Stokes’a. Otrzymane równania mają postać [2]:
𝜕𝑇𝑖𝑗′
𝜕𝑡 + 𝑢𝑘𝜕𝑇𝑖𝑗′
𝜕𝑥𝑘 = −𝑇𝑖𝑘′ 𝜕𝑢𝑗
𝜕𝑥𝑘− 𝑇𝑖𝑘′ 𝜕𝑢𝑖
𝜕𝑥𝑘+ 𝜀𝑖𝑗 − 𝛱𝑖𝑗 + 𝜕
𝜕𝑥𝑘[𝜈𝜕𝑇𝑖𝑗′
𝜕𝑥𝑘 + 𝐶𝑖𝑗𝑘], (26) przy czym
𝛱𝑖𝑗 =𝑝′ 𝜌 ∙ (𝜕𝑢𝑖′
𝜕𝑥𝑗+𝜕𝑢𝑗′
𝜕𝑥𝑖), (27) 𝜀𝑖𝑗 = 2𝜈 ∙𝜕𝑢𝑖′
𝜕𝑥𝑘
𝜕𝑢𝑗′
𝜕𝑥𝑘, (28)
6 𝐶𝑖𝑗𝑘 = 𝑢𝑖′𝑢𝑗′𝑢𝑘′ +1
𝜌(𝑝′𝑢𝑖′𝛿𝑗𝑘+ 𝑝′𝑢𝑗′𝛿𝑖𝑘). (29) Π to tzw. tensor naprężeń ciśnieniowych, a δij to funkcja delta Kroneckera.
1.6. Modele turbulencji wykorzystane w symulacji
W symulacji wykorzystane zostały cztery różne modele turbulencji z palety modeli dostępnej w programie ANSYS CFX. Były to modele:
• k-epsilon (k-ε) w wersji standardowej,
• k-omega (k-ω) w wersji standardowej,
• SST,
• SSG Reynolds.
Poniżej zaprezentowany został zwięzły opis tych modeli turbulencji. Równania rozwiązywane przez solvery programu ANSYS CFX można znaleźć w literaturze, m.in. na stronie CFD-Online [1] i w podręczniku Wilcoxa [3, rozdz. 4.3 i 6.3].
1.7. Modele k-epsilon i k-omega
Model turbulencji k-epsilon to jeden z dwurównaniowych modeli z grupy RANS. Występują w nim dwa parametry, mianowicie k i ε, które zostały zdefiniowane w 1.3. Model k-epsilon należy do grupy metod opartych na założeniach Boussinesq’a oraz na pojęciu lepkości turbulentnej µT, która zdefiniowana jest w sposób uwikłany wzorem
𝑇𝑖𝑗′ = 2𝜇𝑇
𝜌 𝑆𝑖𝑗 −2
3𝑘𝛿𝑖𝑗 = 2𝜇𝑇 𝜌 (𝜕𝑢𝑖
𝜕𝑥𝑗 +𝜕𝑢𝑗
𝜕𝑥𝑖) −2
3𝑘𝛿𝑖𝑗, (30) gdzie 𝑺 to średni tensor szybkości deformacji.
W standardowym modelu k-epsilon lepkość turbulentna jest obliczana wg wzoru 𝜇𝑇 = 𝜌𝐶𝜇𝑘2
𝜀 , (31) przy czym Cµ = 0,09 = const.
W modelu k-omega jedno z równań różniczkowych odnosi się do wielkości ω, czyli tzw. właściwej szyb- kości dyssypacji energii. W przypadku standardowego modelu Wilcoxa jest ona zdefiniowana wzorem:
𝜔 = 𝜀
𝑘𝛽∗, (32) przy czym β* = 0,09.
Wadą modelu k-omega jest wrażliwość na zmiany warunków brzegowych na wlocie przy przepływach swobodnych [1].
1.8. Model SST
Model SST (Menter’s Shear Stress Transport) również jest dwurównaniowy i łączy w sobie oba wcześniej omówione modele. W uproszczeniu połączenie polega na tym, że w obszarze warstwy przyściennej obliczenia są wykonywane zgodnie z modelem k-omega, a w rdzeniu płynu zgodnie z modelem k-epsilon.
1.9. Model SSG Reynolds
Model ten zalicza się do podgrupy metod RANS określanej mianem RSM (Reynolds stress equation mo- del). SSG oznacza inicjały trzech autorów modelu: Spezialego, Sarkara i Gatski’ego. W odróżnieniu od trzech
7
poprzednich modeli nie korzysta z pojęcia lepkości turbulentnej µT, tylko oblicza bezpośrednio wartości elemen- tów właściwego tensora naprężeń Reynoldsa T’. To pozwala na dokładne symulowanie zjawisk zachodzących przy zanikających lub słabo rozwiniętych przepływach turbulentnych [1]. Uwzględnia także kierunkowe efekty naprężeń Reynoldsa. Model SSG Reynolds charakteryzuje się specyficznym sposobem obliczania tensora naprę- żeń ciśnieniowych Π. Opis tego sposobu można znaleźć w literaturze [3].
2. CEL SYMULACJI
Wykonano symulacje przepływu powietrza przez komorę dwuwymiarową o rozmiarach zbliżonych do rozmiarów rzeczywistych pomieszczeń. Celem było zbadanie wpływu wyboru modelu turbulencji na wyniki ob- liczeń numerycznych. W szczególności badano rozkład linii prądu w komorze oraz profile składowej poziomej wektorów prędkości wzdłuż wybranych prostych przechodzących przez komorę.
Podczas wykonywania analizy wykorzystywano program ANSYS ICEM służący do wykonania modelu geometrycznego oraz do wygenerowania siatki numerycznej, zaś programem ANSYS CFX posłużono się celem przygotowania i przeprowadzenia symulacji.
3. PARAMETRY MODELU
3.1. Geometria układu
Geometria symulowanego układu została zaproponowana przez P. V. Nielsena [4, 5] podczas prac nad aneksem 20 Międzynarodowej Agencji Energii (IEA) zatytułowanym Air Flow Patterns within Buidings. Jest ona obecnie wykorzystywana do testowania modeli numerycznych, m.in. modeli turbulencji. Geometrię stanowi prze- krój przez pomieszczenie, zwane dalej komorą, z wlotem powietrza w lewym górnym rogu i z jego wylotem w prawym dolnym rogu. Schemat komory jest ukazany na rys. 3.1, a rys. 3.2 pokazuje całą geometrię (wraz z od- cinkiem wlotowym i wylotowym) z wymiarowaniem. Linią przerywaną zaznaczono proste, dla których wykre- ślone zostały profile składowej poziomej prędkości.
Rys. 3.1. Komora główna będąca częścią dwuwymiarowej geometrii wykorzystywanej w symulacjach [4].
Oznaczenia: L – długość komory, H – wysokość komory, h – wysokość kanału dolotowego, t – wysokość kanału wylotowego
8
Rys. 3.2. Widok geometrii wykorzystywanej w symulacjach wraz z kanałem wlotowym i wylotowym [opracowanie własne]
3.2. Parametry geometryczne
Komora symulacyjna ma proporcje geometryczne opisane w raporcie [4]. Są one następujące:
• L/H = 3,0,
• h/H = 0,056,
• t/H = 0,16.
Oznaczenia w powyższych proporcjach, jak na rys. 3.1.
Symulacje wykonano dla komory o długości L równej 9,0 m, dlatego pozostałe rozmiary miały wartości:
H = 3,0 m, h = 0,168 m, t = 0,48 m. Model geometryczny, dla którego wygenerowano siatkę i przeprowadzono symulacje, zawierał dodatkowo kanał wlotowy o długości 1,0 m i kanał wylotowy o długości 3,0 m widoczne na rys. 3.2. Grubość przestrzeni symulacyjnej (różnica współrzędnych z-owych obu ścian o warunku brzegowym typu symetria) była równa 0,1 m, a liczba komórek obliczeniowych wzdłuż osi z wynosiła 1. Tym sposobem uzyskano symulację typu dwuwymiarowego.
3.3. Siatka numeryczna
Zastosowana została siatka strukturalna. Długość i/lub wysokość komórek ma rozkład sinusoidalny. Przy ścianach komory oraz wlotu i wylotu komórki są najmniejsze, 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 wraz z odległością od ścian. Taki rozkład wielkości komórek wynika z tego, że dobrą praktyką w stosowaniu metod numerycznych jest zagęszczanie siatki bądź zmniejszanie kroku czasowego w obszarach, względnie w przedziałach czasu, w których spodziewane są duże gradienty badanych wielkości. W takich obsza- rach wymagana jest szczególnie duża dokładność obliczeń. Inną metodą zagęszczania siatki w pobliżu ścian jest utworzenie numerycznego obszaru warstwy przyściennej. W programie ICEM nie korzystano jednak z tej opcji.
W tabelach 3.1 i 3.2 przedstawiono wybrane parametry siatki numerycznej, zaś widok siatki w kierunku osi z jest przedstawiony na rys. 3.3.
9
Tabela 3.1. Parametry siatki numerycznej Rodzaj Strukturalna (heksaedralna)
Liczba węzłów 22190
Liczba komórek 10820
Objętość 2,8603 m3
Współrzędne
xmin 0 m xmax 13,0 m ymin 0 m ymax 3,0 m zmin 0 m zmax 0,1 m
Tabela 3.2. Parametry fragmentów siatki numerycznej Nazwa powierzchni
w programie CFX Liczba węzłów (nodes) Liczba komórek oblicze- niowych (elements)
Maksymalna smukłość
Inlet 30 14 47,4821
Outlet 42 20 100
Wall_up 510 254 100
Wall_down 522 260 139,905
Symmetr_1 11095 10820 139,905
Symmetry_2 11095 10820 139,905
Rys. 3.3. Widok geometrii z boku z zaznaczonymi liniami siatki numerycznej
3.4. Warunki symulacji
Symulacja dotyczyła przepływu ustalonego powietrza o temperaturze 25°C traktowanego jako gaz dosko- nały. Przepływ był izotermiczny, a działanie siły grawitacji zostało pominięte (moduł wyporności ‘buoyancy’ był wyłączony). Stała wartość współczynnika lepkości dynamicznej równy 1,831 ∙ 10-5 Pa∙s. Ciśnienie referencyjne w komorze wynosiło 1 atm. Zastosowano inicjalizację automatyczną. Rodzaj wykorzystywanych warunków brze- gowych na poszczególnych powierzchniach symulowanej domeny pokazano w tabeli 3.3.
10
Tabela 3.3. Warunki brzegowe Nazwa powierzchni w pro-
gramie 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 - -
Symmetr_1 Symmetry - -
Symmetry_2 Symmetry - -
3.5. Kryteria zakończenia obliczeń numerycznych
Kryterium zakończenia obliczeń iteracyjnych była średniokwadratowa wartość znormalizownych residuów równań bilansu masy i pędu na poziomie 1∙10-6.
Na wykresach 3.1–3.4 przedstawiono wartości znormalizowanych residuów podczas kolejnych kroków iteracyjnych dla równań opisujących zasadę zachowania masy oraz pędu, rozwiązywanych iteracyjnie przez solver programu ANSYS CFX.
Na podstawie wykresów residuów można stwierdzić, że spełnienie warunku zakończenia obliczeń wyma- gało najmniejszej liczby iteracji, gdy wykorzystywany był model turbulencji k-ε (600 iteracji). Pozostałe badane modele wymagały przeprowadzenia większej liczb iteracji, mianowicie:
• model SST – 1112,
• model k-ω – 1661,
• model SSG Reynolds Stress – 1909.
Wykres 3.1. Zależność wartości średniokwadratowych znormali- zowanych residuów równania ciągłości i równań zachowania pędu od numeru iteracji; symulacja dla modelu k-ε
Wykres 3.2. Zależność wartości średniokwadratowych znormalizo- wanych residuów równania ciągłości i równań zachowania pędu od numeru iteracji; symulacja dla modelu k-ω
11
Wykres 3.3. Zależność wartości średniokwadratowych znormalizo- wanych residuów równania ciągłości i równań zachowania pędu od numeru iteracji; symulacja dla modelu SST
Wykres 3.4. Zależność wartości średniokwadratowych znormalizo- wanych residuów równania ciągłości i równań zachowania pędu od numeru iteracji; symulacja dla modelu SSG Reynolds
4. WYNIKI OBLICZEŃ I ANALIZA
Na rys. 4.1–4.4 przedstawiono pola wektorów prędkości, na rys. 4.5–4.8 linie prądu, a na rys. 4.9–4.12 pola modułu prędkości w symulowanym układzie dla czterech badanych modeli turbulencji.
Rys. 4.1. Pole prędkości – model k-ε
Rys. 4.2. Pole prędkości – model k-ω
12
Rys. 4.3. Pole prędkości – model SST
Rys. 4.4. Pole prędkości – model SSG Reynolds
Rys. 4.5. Linie prądu – model k-ε
13
Rys. 4.6. Linie prądu – model k-ω
Rys. 4.7. Linie prądu – model SST
Rys. 4.8. Linie prądu – model SSG Reynolds
14 4.1. Charakterystyka wirów
Rysunki 4.1–4.12 pozwalają stwierdzić, że wybór modelu turbulencji ma istotny wpływ na przewidywanie zjawisk zachodzących podczas przepływu gazu przez komorę testową. Podczas eksperymentów [5, 6] dowie- dziono, ż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ętny”) oraz dwa znacznie mniejsze, mające zwrot wekto- rów prędkości przeciwny do ruchu wskazówek zegara (dalej nazywane „lewoskrętnymi”). Dwa mniejsze wiry stanowią tzw. strefy recyrkulacji i zlokalizowane są w lewym dolnym i w prawym górnym rogu komory. Istnienie stref recyrkulacji potwierdzają wykresy 4.3 i 4.4. Na wykresie 4.3, obrazującym dolną część komory, punkty od- powiadające danym empirycznym w przedziale [1,0; 2,1] m wskazują na dodatnią wartość składowej poziomej u wektora prędkości, zaś w głównym prawoskrętnym wirze wartość u jest ujemna. Wykres 4.4, ukazujący profil u/u0 w górnej części komory, zawiera punkty danych empirycznych o ujemnej wartości składowej poziomej wek- tora prędkości w przedziale [9,4; 10,0] m, natomiast w górnej części głównego wiru składowa u jest dodatnia.
Wyniki uzyskane przy wykorzystaniu standardowego modelu k-ε (rys. 4.1 i 4.5) odzwierciedlają istnienie wszyst- kich wirów stwierdzonych empirycznie (nr 1, 2 i 3, rys. 4.5), 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 suge- rują, że w komorze powstają dodatkowe wiry.
Wszystkie zastosowane modele wykazują zgodność w kwestii położenia centrum największego wiru. Znaj- duje się ono w przybliżeniu na wysokości 1,5 m w odległości ok. 6,5 m od lewej ściany głównej komory. W przy- padku modelu k-ε linie prądu głównego wiru (nr 1, rys. 4.5) wypełniają prawie całą komorę główną, za 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, rys. 4.5). Dla lepszego zobrazowania obszarów bezwirowych, na rysunku linii prądu dla modelu k-ε zostały one oddzielone od obszarów wirowych za pomocą czarnych kreskowanych krzywych. Defektem modelu k-ε jest niedoszacowanie wielkości wirów nr 2 i 3.
W przypadku zastosowania modelu k-ω (rys. 4.2 i 4.6) wyniki symulacji sugerują występowanie dodatko- wego prawoskrętnego wiru (nr 4, rys. 4.6) w lewym dolnym rogu komory. Zajmuje on obszar o szerokości ok. 40 cm (wykres 4.3). Rozmiar wiru lewoskrętnego (nr 3, rys. 4.6), będącego odpowiednikiem strefy recyrkula- cji 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. Podobna uwaga dotyczy wiru nr 2, którego rozmiar obliczony przy pomocy modelu k-ω jest zbyt duży. Można stwierdzić, że model k-ω źle oddaje charakter przepływu powietrza przez lewą dolną część komory testowej.
Model SST zastosowany do symulacji daje w wyniku cztery wiry; wir oznaczony numerem 3 na rys. 4.7 jest lewoskrętny i zajmuje całą lewą dolną część komory. Rozmiar lewej dolnej strefy recyrkulacji obliczony z wykorzystaniem modelu SST jest przeszacowany, podobnie jak w przypadku użycia modelu k-ω. Wir oznaczony numerem 4 zlokalizowany jest nad wirem nr 3, a jego orientacja jest taka, jak wiru głównego nr 1. Można zauwa- żyć, że wielkość i kształt obszaru zajmowanego przez wir nr 2 są zbliżone do parametrów odpowiadającego mu wiru przewidzianego przez model k-ω.
Symulacja z wykorzystaniem modelu SSG Reynolds daje wyniki sugerujące występowanie największej liczby wirów. Na rysunku linii prądu (rys. 4.8) widać ich pięć lub sześć. Rozstrzygnięcie tego, czy obszar ozna- czony numerem 6 na rys. 4.8 jest osobnym wirem, czy tylko odgałęzieniem głównego wiru nr 1, wymagałoby dokładniejszych rysunków; tak szczegółowy opis nie jest jednak wymagany do stwierdzenia, który model turbu- lencji daje wyniki najbardziej zbliżone do eksperymentalnych. Można zauważyć, że linie prądu wiru nr 1 z dala od centrum nie mają kształtu eliptycznego, co wynika z tego, że opływają one obszar zajęty przez wiry 3, 4 i 5.
Charakterystyczne jest to, że wielkość obszaru zajmowanego przez wir nr 2 w prawym górnym rogu komory jest nieco mniejsza niż w przypadku zastosowania modeli SST i k-ω, ale większa niż w wynikach uzyskanych przy pomocy modelu k-ε.
15
Rys. 4.9. Pole modułu prędkości – model k-ε
Rys. 4.10. Pole modułu prędkości – model k-ω
Rys. 4.11. Pole modułu prędkości – model SST
16
Rys. 4.12. Pole modułu prędkości – model SSG Reynolds
Na rys. 4.4 i 4.12 można zauważyć, ż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, zaś 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łów- nego wiru nr 1 niż pozostałe modele turbulencji.
4.2. Składowe wektora prędkości wzdłuż wybranych profili
Na wykresach 4.1–4.4 przedstawione są profile składowej poziomej u prędkości płynu podzielonej przez wartość prędkości wlotowej u0, pochodzące z symulacji wykorzystujących ww. modele turbulencji oraz pocho- dzące z eksperymentu. Dane eksperymentalne są opatrzone słupkami błędu. Wartości błędu pochodzą z danych eksperymentalnych dostępnych na stronie internetowej [6]. Długości słupków błędu nad i pod punktami pomiaro- wymi są równe i wynoszą u’RMS/u0,gdzie u’RMS to średniokwadratowa wartość prędkości fluktuacyjnej w danym punkcie. Prędkość płynu na wlocie do komory (u0) była zadana w symulacji jako warunek brzegowy. Wykresy 4.5–4.8 ukazują z kolei profile składowej pionowej v prędkości płynu podzielonej przez wartość prędkości wloto- wej u0. Wykresy potwierdzają, że model turbulencji k-ε standard najlepiej opisuje przepływ gazu przez analizo- waną komorę, ponieważ w przypadku tego modelu profile prędkości wzdłuż wybranych linii są najbardziej zbli- żone do punktów eksperymentalnych. Widoczne jest to na wykresach 4.1–4.3. Wyjątkami są tu dwie strefy: w pra- wym górnym (wyk. 4.4, x od 9,4 do 10,0 m) oraz w lewym dolnym rogu komory (wyk. 4.3, 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 do- datnich 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 z ww. stref model k-ε sugeruje, że składowa u jest ujemna, czyli że jest to obszar głównego prawoskrętnego wiru, zaś z doświadczenia wynika, że występuje tam strefa lewoskrętnej recyr- kulacji. Przyczyną tych niezgodności jest wspomniany w punkcie 4.1 błąd modelu k-ε polegający na niedoszaco- waniu wielkości wirów nr 2 i 3.
W prawym górnym rogu komory testowej najlepsze wyniki daje 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.
17
Wykres 4.1. Profil ilorazu u/u0 wzdłuż linii pionowej x = 4,00 m (u0 = 0,455 m∙s-1)
Wykres 4.2. Profil ilorazu u/u0 wzdłuż linii pionowej x = 7,00 m (u0 = 0,455 m∙s-1)
Wykres 4.3. Profil ilorazu u/u0 wzdłuż linii poziomej y = 0,084 m (u0 = 0,455 m∙s-1)
Wykres 4.4. Profil ilorazu u/u0 wzdłuż linii poziomej y = 2,916 m (u0 = 0,455 m∙s-1)
-0,40 -0,20 0,00 0,20 0,40 0,60 0,80 1,00
0,00 0,50 1,00 1,50 2,00 2,50 3,00
u/u0
y [m]
x = 4 m, k-ε x = 4 m, k-ω x = 4 m, SST
x = 4 m, SSG Reynolds x = 4 m, empiryczne
-0,60 -0,40 -0,20 0,00 0,20 0,40 0,60 0,80 1,00
0,00 0,50 1,00 1,50 2,00 2,50 3,00
u/u0
y [m]
x = 7 m, k-ε x = 7 m, k-ω x = 7 m, SST
x = 7 m, SSG Reynolds x = 7 m, empiryczne
-0,60 -0,50 -0,40 -0,30 -0,20 -0,10 0,00 0,10 0,20 0,30 0,40
1,00 3,00 5,00 7,00 9,00
u/u0
x [m]
y = 0,084 m, k-ε y = 0,084 m, k-ω y = 0,084 m, SST
y = 0,084 m, SSG Reynolds y = 0,084 m, empiryczne
-0,60 -0,40 -0,20 0,00 0,20 0,40 0,60 0,80 1,00 1,20 1,40
1,00 3,00 5,00 7,00 9,00
u/u0
x [m]
y = 2,916 m, k-ε y = 2,916 m, k-ω y = 2,916 m, SST
y = 2,916 m, SSG Reynolds y = 2,916 m, empiryczne
18
Wykres 4.5. Profil ilorazu v/u0 wzdłuż linii pionowej x = 4,00 m (u0 = 0,455 m∙s-1)
Wykres 4.6. Profil ilorazu v/u0 wzdłuż linii pionowej x = 7,00 m (u0 = 0,455 m∙s-1)
Wykres 4.7. Profil ilorazu v/u0 wzdłuż linii poziomej y = 0,084 m (u0 = 0,455 m∙s-1)
Wykres 4.8. Profil ilorazu v/u0 wzdłuż linii poziomej y = 2,916 m (u0 = 0,455 m∙s-1)
Na wykresach 4.5–4.8 przedstawione są wartości składowej pionowej v prędkości podzielonej przez u0
wzdłuż wybranych prostych. Wykres 4.6 pozwala stwierdzić, że 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 składowej poziomej osiągają wartość bezwzględną 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: 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.
-0,04 -0,02 0,00 0,02 0,04 0,06 0,08 0,10 0,12 0,14
0,00 0,50 1,00 1,50 2,00 2,50 3,00
v/u0
y [m]
x = 4 m, k-ε x = 4 m, k-ω x = 4 m, SST
x = 4 m, SSG Reynolds
-0,05 0,00 0,05 0,10 0,15 0,20 0,25 0,30
0,00 0,50 1,00 1,50 2,00 2,50 3,00
v/u0
y [m]
x = 7 m, k-ε x = 7 m, k-ω x = 7 m, SST
x = 7 m, SSG Reynolds
-0,10 -0,08 -0,06 -0,04 -0,02 0,00 0,02 0,04 0,06 0,08
1,00 3,00 5,00 7,00 9,00
v/u0
x [m]
y = 0,084 m, k-ε y = 0,084 m, k-ω y = 0,084 m, SST
y = 0,084 m, SSG Reynolds
-0,10 -0,08 -0,06 -0,04 -0,02 0,00 0,02 0,04 0,06
1,00 3,00 5,00 7,00 9,00
v/u0
x [m]
y = 2,916 m, k-ε y = 2,916 m, k-ω y = 2,916 m, SST
y = 2,916 m, SSG Reynolds
19 4.3. Lepkość turbulentna
Rysunki 4.13–4.16 ukazują pole lepkości turbulentnej wewnątrz badanej domeny. Widać, że wartość maksymalna lepkości turbulentnej jest największa w przypadku modeli turbulencji k-ε i k-ω, dla których wynosi ok. 0,015 Pa∙s. Wartość maksymalna obliczona z użyciem modeli SST i SSG Reynolds jest ok. dwukrotnie mniejsza. Można zaobserwować, że obszary o stounkowo wysokich wartościach lepkości turbulentnej dla poszczególnych modeli pokrywają się w przybliżeniu z centrami głównych wirów.
Rys. 4.13. Pole lepkości turbulentnej – model k-ε
Rys. 4.14. Pole lepkości turbulentnej – model k-ω
Rys. 4.15. Pole lepkości turbulentnej – model SST
20
Rys. 4.16. Pole lepkości turbulentnej – model SSG Reynolds
5. WNIOSKI
• Wybór modelu turbulencji stosowanego przez solvery w programach CFD ma wpływ na wyniki przepro- wadzanych symulacji. Nieodpowiedni dobór modelu turbulencji do opisywania zjawisk przepływowych danego rodzaju może prowadzić do błędnych wyników znacząco odbiegających od rzeczywistego charak- teru zjawiska.
• Stwierdzenie, który model turbulencji najtrafniej opisuje przepływ danego typu, jest możliwe tylko wtedy, gdy dysponuje się danymi eksperymentalnymi. Pomiary doświadczalne powinny być wykonane w kilku różnych miejscach badanej geometrii, np. wzdłuż kilku prostych, a zagęszczenie punktów pomiarowych powinno być odpowiednio duże.
• 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ń. W efekcie, z jednej strony umożliwia on opis wirów o małej intensywności, np. takich zanikających, ale z drugiej strony jest on podatny na lokalne, chwilowe zmiany wartości składowych tensora. To może powodować obliczenie cyrkulujących pól prędkości tam, gdzie w rzeczywistości nie występują.
• W przypadku symulowanej geometrii, którą stanowi prostopadłościenna komora, stosowanie modelu SSG Reynolds prowadzi do nadprodukcji wirowości.
• Najlepszą zgodność wyników symulacji z danymi eksperymentalnymi daje wybór standardowego modelu k-epsilon.
• Wszystkie cztery symulacje zostały przeprowadzone z użyciem tej samej siatki numerycznej. Być może dobór siatki o większej liczbie komórek lub o innym rozkładzie ich wielkości spowodowałby, że wyniki symulacji w mniejszym stopniu obiegałyby od pomiarów doświadczalnych, w szczególności w przypadku modeli SST i SSG Reynolds.
6. BIBLIOGRAFIA
1. https://www.cfd-online.com/Wiki/Turbulence_modeling (dostęp 07.01.2019).
2. Jaworski Zdzisław, Numeryczna mechanika płynów w inżynierii chemicznej i procesowej. AOW EXIT, Warszawa, 2005.
3. Wilcox David C., Turbulence modeling for CFD. Wyd. 2. DCW Industries, Inc., Anaheim, Kalifornia, 1998.
4. Room Air and Contaminant Flow, Evaluation of Computational Methods. Subtask-1 Summary Report, edytor:
Lemaire A.D. Wyd. 2. Delft (Holandia), grudzień 1993.
5. Nielsen Peter 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).
