7. Badania numeryczne CFD
7.1. Opracowanie modelu numerycznego CFD dla geometrii i warunków pracy wentylatora referencyjnego
7.1.1. Opis budowanego modelu numerycznego
Przestrzenna geometria (3D) modelowanego wentylatora przygotowana została w środowisku CAD. Wyróżniono trzy podstawowe bryły obejmujące obszar płynu w maszynie, tj. pojedynczy kanał łopatkowy, obudowę spiralną i lej wlotowy z odcinkiem rurociągu, co pokazano na rysunku 7.1.
Geometria CAD została zaimportowana do środowiska ANSYS Design Modeler, gdzie zidentyfikowano poszczególne bryły i opisano powierzchnie. W kolejnym kroku prze-prowadzono dyskretyzację modelu z wykorzystaniem narzędzi do generowania siatki. Siatkę obudowy spiralnej i kanału dolotowego wygenerowano w oprogramowaniu ANSYS Meshing [7,96,100,104]. W obszarach stref przyściennych wprowadzono zagęszczanie siatki poprzez tzw. warstwę inflacji, o zdefiniowanej grubości pierwszego elementu i zadanym współczyn-niku rozrostu równym 1,2 dla kolejnych 10 warstw. Szczególnie ważne z punktu widzenia
S t r o n a | 109
rozwiązania modelu jest wygenerowanie wysokiej jakości siatki w obrębie wirnika. Siatkę tę opracowano z wykorzystaniem narzędzia TurboGrid, umożliwiającego generowanie siatek strukturalnych o najwyższych parametrach jakościowych dla rotorów i statorów maszyn roto-dynamicznych [8].
Rys.7.1. Uproszczenie modelu rzeczywistego wentylatora na potrzeby symulacji numeryczne
(A-obudowa spiralna, B-kanał łopatkowy, C-odcinek rurociągu wlotowego)
Przepływ płynu przez wentylator opisany został dwoma podstawowymi równaniami bilanso-wymi mechaniki płynów [6,17,62,91,143]. Równaniem ciągłości (7.1):
𝑑𝜌
𝑑𝑡+ 𝑑𝑖𝑣(𝜌𝑣⃗) = 0 (7.1)
Równaniem zachowania pędu (równie Naviera-Stokesa)[6,17,62,91] wyrażono przez (7.2): 𝜌𝑑𝑣⃗ 𝑑𝑡 = 𝜌𝐹⃗𝑚+ 𝑑𝑖𝑣 [(−𝑝 − 2 3𝜇 𝑑𝑖𝑣(𝑣⃗)) 𝐼 + 2𝜇𝑆] (7.2) gdzie: 𝜌 ˗ gęstość płynu, 𝑣⃗ ˗ wektor prędkości,
𝐹𝑚 ˗ jednostkowa siła masowa, 𝑝 ˗ ciśnienie,
𝜇 ˗ lepkość dynamiczna płynu, I ˗ tensor jednostkowy,
S ˗ tensor prędkości deformacji
W modelu, celem przyspieszenia obliczeń założono, że przepływ jest izotermiczny, a czynnikiem jest powietrze o temperaturze T=25⁰C. Z badań bilansowych wiadomo, że przy-rost ciśnienia w badanym wentylatorze jest rzędu 1000 Pa. Co przy założeniu izotermicznego przepływu daje zmianę gęstości zgodnie z zależnością (7.3) na poziomie 1%.
𝑝 𝜌= 𝑅𝑖𝑇 = 𝑖𝑑𝑒𝑚 (7.3) 𝜌2 = 𝜌1∙𝑝2 𝑝1 = 𝜌1∙ 101000 100000= 1,01𝜌1
Badani a nu mer yc zne C FD
S t r o n a | 110
Założenie stałej gęstości ρ=const. (1,185 kg/m3) jak i lepkość dynamicznej μ=const. (1,831·10-5 Pa·s) skutkuje uproszczeniem równań transportu, odpowiednio do postaci (7.1a) i (7.2b): 𝑑𝑖𝑣(𝑣⃗) = 0 (7.1a) 𝑑𝑣⃗ 𝑑𝑡 = 𝐹⃗𝑚− 1 𝜌𝑔𝑟𝑎𝑑 𝑝 + 𝜗 ∇2𝑣⃗ (7.2a)
W równaniu (7.2a) wprowadzono lepkość kinematyczną 𝜗 = 𝜇/𝜌 oraz operator Laplaca ∇2. Rozwiązanie równań ruchu Naviera-Stokesa w sposób bezpośredni, czyli przeprowa-dzenie obliczeń typu DNS (Direct Numerical Simulations), wymaga bardzo gęstych siatek obliczeniowych i długiego czasu rozwiązania (obliczenia niestacjonarne). Dla przedmiotowe-go zagadnienia symulacja DNS jest bardzo trudna do przeprowadzenia – duże wymiary obiektu, wysokie prędkości przepływu determinują użycie ogromnych zasobów komputero-wych. W obliczeniach praktycznych (projektowych) wykorzystuje się najczęściej metodę uśredniania zmiennych względem czasu, prowadzącą do równań RANS (Reynolds-Averaged Navier–Stokes equations). Model RANS do obliczeń przepływów wentylatorowych zweryfi-kowano w licznych publikacjach naukowych [4,129,27,54,55,98] wskazanych w przeglądzie literatury. Przy założeniu średnio-ustalonego ruchu turbulentnego płynu, każdy z parametrów przepływowych w danej chwili czasu można zapisać jako sumę wartości uśrednionej i jej fluktuacji, jak w (7.4):
𝑝 = 𝑝 + 𝑝′
𝑣 = 𝑣 + 𝑣′ (7.4)
Kolejnym uproszczeniem, które przyspiesza rozwiązanie bez wpływu na jego jakość, w przypadku gazu o małej gęstości i przepływie izotermicznym, jest pominięcie w modelu numerycznym sił grawitacji. Ostatecznie po wyeliminowaniu sił ciężkości, równania (7.1a) i (7.2a) w zapisie wskaźnikowym dla przepływu średnio-ustalonego przyjmują postać (7.1c) i (7.2c): 𝜕𝑣̅𝑖 𝜕𝑥𝑖 = 0 (7.1c) 𝜕 𝜕𝑥𝑗(𝜌𝑣̅𝑖𝑣̅𝑗) = − 𝜕𝑝 𝜕𝑥𝑖 + 𝜕 𝜕𝑥𝑗[𝜇 ( 𝜕𝑣̅𝑖 𝜕𝑥𝑗+ 𝜕𝑣̅𝑗 𝜕𝑥𝑖)] − 𝜕𝑅𝑖𝑗 𝜕𝑥𝑗 (7.2c)
gdzie Rij jest tensorem naprężeń Reynoldsa, związanych z uśrednioną prędkością zmiany pędu na skutek pulsacji turbulentnych. W kartezjańskim układzie współrzędnych tensor ten ma postać zgodną z równaniem (7.5):
𝑅𝑖𝑗 = −𝜌𝑣𝑗′𝑣𝑖′ = [ −𝜌(𝑣𝑥′)2 −𝜌𝑣𝑥′𝑣𝑦′ −𝜌𝑣𝑥′𝑣𝑧′ −𝜌𝑣𝑦′𝑣𝑥′ −𝜌(𝑣𝑦′)2 −𝜌𝑣𝑦′𝑣𝑧′ −𝜌𝑣𝑧′𝑣𝑥′ −𝜌𝑣𝑗′𝑣𝑖′ −𝜌(𝑣𝑧′)2] (7.5)
S t r o n a | 111
Układ równań (7.1.c) i (7.2.c) nie jest układem zamkniętym, w tensorze Reynoldsa brak jest sześciu równań określających poszczególne składowe. W niniejszej pracy problem ułożenia dodatkowych równań rozwiązano poprzez obliczanie wpływu turbulencji według podejścia Boussinesqa [73]. Według tej hipotezy, istnieje związek między składowymi tensora Rey-noldsa, a tensorem prędkości deformacji. W zapisie wskaźnikowym, związek ten wyraża się równaniem (7.6).
𝑅𝑖𝑗 = −𝜌𝑣𝑗′𝑣𝑖′ = 𝜇𝑇(𝜕𝑣̅𝑖 𝜕𝑥𝑗 +
𝜕𝑣̅𝑗
𝜕𝑥𝑖) (7.6)
Współczynnikiem wiążącym naprężenia z polem prędkości jest lepkość turbulentna 𝜇𝑇. Jest to funkcja skalarna (często nieliniowa) wielu zmiennych, takich jak własności fizyko-chemicznych płynu czy charakteru przepływu. Istnieje wiele empirycznych/pół-empirycznych modeli matematycznych, rozwiązujących zagadnienie lepkości turbulentnej. Wiele z nich zostało zaadoptowanych do komercyjnych kodów komputerowych. Jako uzupełnienie równań RANS, w przedmiotowym zagadnieniu przyjęto dwurównaniowy model turbulencji SST k-ω opracowany przez Mentera [73,74]. Model SST (shear stress transport) jest modelem hybry-dowym, łączący najlepsze cechy modelu k-ε i k-ω. Model k-ε dobrze sprawdza się w przepływie swobodnym, natomiast model k-ω znacznie lepiej odwzorowuje przepływ w warstwie przyściennej. Dodatkowo SST wprowadza człon limitujący „nadprodukcję” ener-gii kinetycznej turbulencji w obszarach silnych dodatnich gradientów ciśnienia – cecha pożą-dana szczególnie w przepływach o nagłej zmianie kierunku, takich jak rozpatrywany prze-pływ przez wirnik promieniowy [32,47,86]. Model matematyczny uzupełniono warunkami jednoznaczności (brzegowymi i początkowymi) zgodnie z rysunkiem 7.2.
Rys.7.2. Definicja domen i wybranych warunków brzegowych
Zdefiniowano trzy domeny, dwie domeny stacjonarne, odpowiednio rurociąg wlotowy i obudowę spiralną oraz domenę obrotową ( ang. Rotating Domain). Warunki brzegowe zdefiniowano zgodnie z rysunkiem 7.2, przyjmując na wylocie z obudowy ciśnienie barome-tryczne równe 105 Pa, prędkość obrotową wirnika 2880 obr./min i zerową prędkość płynu na ściankach. W celu wyznaczenia charakterystyki przepływowej maszyny, na wlocie wprowa-dzano jako zmienny parametr strumień masy. Kontrolę eksperymentu prowadzono z użyciem ANSYS Design Exploration w procedurze „what-if” [6,7,103,104,116], gdzie zmienną
wej-Badani a nu mer yc zne C FD
S t r o n a | 112
ściową był strumień masy, a zmiennymi wyjściowymi przyrost ciśnienia statycznego i sprawność wewnętrzna wentylatora.
Ważnym elementem dla poprawności i stabilności obliczeń jest wygenerowania siatki o właściwych wymiarach i parametrach geometrycznych. Szczególną uwagę należy zwrócić na utrzymanie właściwej wartości Y+ (bezwymiarowej odległości od ścianki) [10,93], dla przyjętego w symulacji modelu turbulencji. W celu poprawnego modelowania podwarstwy laminarnej Y+, dla modelu turbulencji SST k-ω powinien przyjmować wartość równą 1 (w niektórych opracowaniach autorzy sugerują utrzymanie Y+ < 5).
W celu oceny jakości siatki i jej wpływu na rozwiązanie autor przeprowadził bania wrażliwości modelu pod tym kątem. Zdefiniowano dwa parametry siatki, które poddano weryfikacji obliczeniowej tj. wysokość pierwszego elementu siatki w warstwie przyściennej i wielkości elementów w przepływie rozwiniętym. Badania prowadzono dla stałej wartości strumienia masy, obserwując zmianę parametrów obliczanych wraz ze zmianą generowanej siatki numerycznej. Na rysunku 7.3. zestawiono wybrane wyniki przeprowadzonego ekspe-rymentu. Rysunek 7.3a pokazuje zmianę wartości Y+ na powierzchni łopatek w funkcji wysokości pierwszego elementu siatki w warstwie inflacji. Jako krytyczny parametr oblicze-niowy dla maszyny wybrano sprawność wewnętrzną, którą poddano badaniu wrażliwości na jakość generowanej siatki. Sprawność wewnętrzna jest wynikiem innych wielkości oblicza-nych tj. przyrostu ciśnienia i momentów sił na ścianach wirnika, które mocno zależą od modelowanego charakteru przepływu wokół łopatek (detekcja i modelowanie oderwań). Na rysunku 7.3b przedstawiono wpływ liczby elementów siatki (wymiaru elementów) na stabil-ność i poprawstabil-ność obserwowanego parametru tj. numerycznego względnego błędu obliczania sprawności wewnętrznej, zdefiniowanego zależnością (7.7):
𝛿𝜂 =𝜂𝑑𝑝 − 𝜂𝑌=1
𝜂𝑌=1 ∙ 100% (7.7)
gdzie:
𝜂𝑑𝑝 ˗ sprawność określona dla danej próby (dp = design point), 𝜂𝑌=1 ˗ sprawność dla siatki rekomendowanej przy Y+ = 1,
a) b)
Rys.7.3. Wpływ siatki na rozwiązanie, a) zależność Y+ od wysokości pierwszego elementu siatki,
S t r o n a | 113
Dodatkowo na rysunku 7.4. zestawiono wyniki obliczeń przepływu w wirniku dla dwóch skrajnych wartości Y+. Na rysunku 7.4a średni Y+ na łopatce wynosi 40, na rysunku 7.4b przyjmuje wartość równą 1. W przypadku zagęszczonej siatki, odpowiedniej dla danego modelu turbulencji możliwe jest właściwe zamodelowanie oderwań co widać na rysunku 7.4b.
a) b)
Rys.7.4. Prędkość w kanale międzyłopatkowym dla różnych wartości Y+, a) Y+=40, b) Y+=1
Ostatecznie przyjęto założenie generowania siatki przy wysokości pierwszego elemen-tu 0,01 mm, co daje około 0,6 miliona elementów na jeden kanał międzyłopatkowym mode-lowanego wentylatora. Siatkę dla wirnika przedstawiono na rysunku 7.5.
a) b)
Badani a nu mer yc zne C FD
S t r o n a | 114
Z punktu widzenia zbieżności rozwiązania nie tylko siatka odgrywa znaczącą rolę. Ważne jest również właściwe określenie skali czasowej. Z doświadczeń autora wynika, że automatyczna skala czasowa zadawana w solverze CFX często prowadzi do rozbieżności w rozwiązaniu. Dla modelowanych zagadnień przyjęto fizyczną skalę czasową, zmienną wraz z propagacją rozwiązania - im więcej przeliczonych iteracji, tym bardziej „skomplikowane” struktury przepływowe. Skalę czasową zakładano zgodnie z wyrażeniem (7.8).
𝛥𝑡 = 𝑘
𝜔 (7.8)
gdzie:
𝛥𝑡 ˗ fizyczna skala czasowa [s], 𝜔 ˗ prędkość kątowa wirnika [rad/s],
k ˗ współczynnik przyjmowany zgodnie z wykresem na rysunku 7.6.
liczba iteracji
Rys.7.6. Zmiana wartości współczynnika skali czasowej
7.1.2. Walidacja opracowanego modelu numerycznego w oparciu o badanie bilansowe