» Blog » Schematy dyskretyzacji przestrzennej w ANSYS Fluent

Schematy dyskretyzacji przestrzennej w ANSYS Fluent

Spis Treści:

Wstęp

Rozmowy osób zajmujących się na co dzień symulacją komputerową bywają fascynujące, ale również enigmatyczne. Szczególnie dla osoby postronnej, nie znającej tematu. Takie konwersacje pełne są bowiem dziwacznych, specjalistycznych zwrotów jak: „ubolewam nad skośnością tych komórek”, „zaraza! Ale dywergencja!” albo klasyczne „to już piąty raz mam dzisiaj Segmentation Violation”. Jednym z takim terminów, który wywołują pełne zdziwienia twarze w towarzystwie jest „Schematy dyskretyzacji przestrzennej”. Temat, który często jest bagatelizowany przez osoby symulujące, a który naszym zdaniem warto dobrze poznać.

Dyskretyzacja, czyli podział

Na początek kilka słów wyjaśnienia czym dokładnie są schematy dyskretyzacji przestrzennej. Do symulacji komputerowych oprogramowanie ANSYS Fluent stosuje powszechnie wykorzystywaną w kodach przepływowych Metodę Objętości Skończonych (ang. Finite Volume Method). Technika ta, oparta na kontroli objętości przekształca równania rządzące w równania algebraiczne, które można rozwiązać numerycznie. Transformacja równań i modeli funkcji ciągłych na ich odpowiedniki w formie dyskretnej możliwa jest dzięki podziałowi objętości reprezentującej model na mniejsze, skończone fragmenty (dyskretyzacja przestrzenna). Można stwierdzić zatem, że im większa jest liczba elementów, tym większa jest dokładność przeprowadzonych obliczeń w programie. Zdyskretyzowane w ten sposób równania są rozwiązywane iteracyjnie, aż do osiągnięcia zbieżności.

Rysunek 1. Podział domeny obliczeniowej.

Domyślnie Fluent rozwiązuje wartości dyskretne skalara w środku komórki, jednak w warunkach zbieżności rozwiązania, potrzebne są wartości interpolowane między środkami komórek. Szacuje się je poprzez zastosowanie jednego z dostępnych schematów dyskretyzacji przestrzennej. W zależności od zastosowanych modelów i rozpatrywanych równań schematy te mogą być różne. Najczęściej jednak wymienia się trzy (i na nich skupimy się w tym artykule): First Order Upwind, Second Order Upwind i Quadratic Upwind Interpolation for Convective Kinetics (w skrócie QUICK). Słowo „upwind” w nazwach oznacza że są to tzw. schematy pod wiatr. Uzyskują one wartości nominalne z wartości branych z komórki wyżej w stosunku do kierunku normalnej prędkości przepływu. Poniżej przedstawiono w sposób graficzny, jak działają wspomniane trzy schematy dyskretyzacji.

Rysunek 2. Uproszczona siatka obliczeniowa.

First Order, najprostszy z tej trójki, kopiuje wartość komórki sąsiadującej pod prąd. Schemat pierwszego rzędu jest, jak się można domyśleć, mniej dokładnym algorytmem w porównaniu z schematem drugiego rzędu oraz QUICK. Z drugiej strony, jest szybki i bardziej stabilny. Często wybierany jako punkt wyjścia skomplikowanych rozwiązań. Oczywiście jest teoretyczna możliwość, aby schematem pierwszego rzędu uzyskać rozwiązanie zgodne z rzeczywistością, jednak wymaga to nieskończenie małych elementów. Jest to więc tak samo prawdopodobne jak bezawaryjna jazda starym Oplem nad morze.

Rysunek 3. Zasada działania schematu First Order Upwind.

Przy zastosowaniu schematu dyskretyzacji Second Order, wielkości na ścianach komórek są obliczane przy użyciu metody wielowymiarowej rekonstrukcji liniowej. W tym podejściu dokładność wyższego rzędu jest osiągana na powierzchniach komórek poprzez rozszerzenie szeregu Taylora rozwiązania skoncentrowanego na komórce wokół środka ciężkości komórki. Bardziej ludzkim językiem, oznacza to że na podstawie dwóch punktów pod wiatr, wyznaczana jest funkcja liniowa, którą interpoluje się wartość na interesująca na ścianę.

Rysunek 4. Zasada działania schematu Second Order Upwind.

Dlatego główną różnicą między schematem pierwszego i drugiego rzędu jest liczba punktów użytych do obliczeń, tj. jeden punkt „upstream” dla pierwszego rzędu i dwa dla drugiego rzędu. W praktyce zwykle korzysta się ze schematu drugiego rzędu. Charakterystyczną cechą Second Order jest to, ze lokalnie potrafi przekroczyć teoretyczne granice danej wielkości. Przykładowo, gdy jakaś wielkość ф w domenie powinna zawierać się między 1 a 0, możemy zauważyć wartości poza tym zakresem, np. 1.04.

Schemat QUICK może zapewnić lepszą dokładność niż schemat drugiego rzędu w przypadkach siatek prostokątnych lub hexahedralnych. Głównie za sprawą zastosowania 3 punktów do interpolacji – dwóch pod wiatr i jednej z wiatrem. Dzięki trzem punktom możemy wyznaczyć funkcję kwadratową wielkości, co jest dokładniejsze niż prosta funkcja liniowa. Dla siatek hybrydowych, QUICK działa tylko dla części hexahedralnej, a w elementach tetrahedralnych działa wtedy Second Order.

Rysunek 5. Zasada działania schematu QUICK.

Podsumowując:

First Order Upwind: łatwy do zbieżności i szybki kosztem mniejszej dokładności;

Second Order Upwind: bardziej wymagający i trudniejszy do osiągnięcia zbieżności ale za to dokładniejszy. Zdarza się, że przekracza teoretyczne min i max;

QUICK: schemat dedykowany elementom czworobocznym i sześciościennym. Najdokładniejszy z opisywanych.

Niefizyczna dyfuzja

Największym problemem schematów niższego rzędu jest występowanie zjawiska zwanego dyfuzją numeryczną. Jest to błąd numeryczny charakteryzujący się nadmiernym rozmywaniem się (dyfuzją) wielkości liczonych. Poniżej przedstawiono wynik prostego testu, mającego na celu ukazanie różnic w działaniu różnych schematów dyskretyzacji. W przedstawionym przykładzie analizowano zjawisko dyfuzji numerycznej w równaniu energii.

Jako model testowy wykorzystano symetryczny model kanału przelotowego powietrza, na który składają się dwa wloty powietrza o różnych temperaturach i prędkościach oraz jeden wylot.

Rysunek 6. Podgląd na analizowaną geometrię.

Wlot numer 1: Prędkość V= 0.2 m/s , Temperatura T= 10 °C

Wlot numer 2: Prędkość V= 0.4 m/s , Temperatura T= 90 °C

Do testów wykorzystano cztery siatki:

NR 1: Tetrahedralna o rozmiarze elementu 10 mm (ok 100 tysięcy elementów)

NR 2: Tetrahedralna o rozmiarze elementu 3.5 mm (ok 500 tysięcy elementów)

NR 3: Hybrydowa o rozmiarze elementu 7 mm (ok 100 tysięcy elementów)

NR 4: Hybrydowa o rozmiarze elementu 3.5 mm (ok 400 tysięcy elementów)

Rysunek 7. Podgląd na analizowane siatki w osi symetrii.

Na początek, porównajmy wyniki z dwóch schematów dyskretyzacji dla pierwszej zgrubnej siatki. Dlaczego nie trzech? Cóż, z definicji schemat QUICK będzie działał tylko dla siatek hexahedralnych. Dla siatki tetrahedralnej działa tak samo jak 2nd Order.

Rysunek 8. Podgląd na analizowane siatki. Widok na wylot.

To co powinniśmy zaobserwować to strugę gorącego powietrza, która będzie rozpraszać się i mieszać z zimnym powietrzem. Pytanie tylko, jak mocno struga ta będzie się rozpraszać. Przy niższych schematach dyskretyzacji powinniśmy zauważyć sztucznie zwiększone mieszanie, na skutek błędu obliczeniowego (dyfuzji numerycznej).

Rysunek 9. Rozpraszająca się struga gorącego powietrza.

Okazuje się, że schemat dyskretyzacji ma olbrzymie znaczenie na wynik. Poniższy rysunek przedstawia różnice między First Order Upwind a Second Order Upwind w wartościach temperatury na płaszczyźnie symetrii.

Rysunek 10. Kontury temperatury na symetrii modelu. Górny widok 1st Order, dolny 2nd Order.

2nd Order odznacza się dłuższym, widocznym śladem oddziaływania sięgającym wylotu kanału. Różnicę widać również na polu prędkości, choć różnice te są dużo mniejsze niż w przypadku pola temperatury.

Rysunek 11. Pola prędkości na symetrii modelu. Górny lewy widok – 1st Order, górny prawy – 2nd Order, dolny – różnica między nimi.

Rozbieżności uwydatniają się na wylocie z domeny. Wyraźnie widać na nim jak struga gorącego powietrza właściwie nie dociera do wylotu, całkowicie rozmywając się po powierzchni. Na konturach prędkości również widać spore różnice, szczególnie w obszarze rdzenia przepływu.

Rysunek 12. Rozkład temperatury na wylocie dla schematu 1st Order (lewa) i 2nd Order (prawa).

Rysunek 13. Rozkład prędkości na wylocie dla schematu 1st Order (lewa) i 2nd Order (prawa).

Co ciekawe, sytuacja niewiele poprawia się po zagęszczeniu siatki. Mimo 5-o krotnego zwiększenia liczby elementów, procentowy błąd w wartościach między schematem pierwszego rzędu i wyższego utrzymuje się na podobnym poziomie. Mniejszy rozmiar elementu powoduje, że kontur poprawia się dla obydwu schematów, jednak różnica między nimi pozostaje podobna. W 1st Order widoczny zaczyna być nawet niewielki skok temperatury w osi kanału, jednak w dalszym ciągu struga jest mocno rozmyta.

Rysunek 14. Kontur temperatury dla 2nd order dla dwóch siatek numerycznych.

Rysunek 15. Rozkład temperatury na wylocie dla schematu 1st Order (lewa) i 2nd Order (prawa). Siatka zagęszczona.

Rysunek 16. Rozkład prędkości na wylocie dla schematu 1st Order (lewa) i 2nd Order (prawa). Siatka zagęszczona.

Teoretycznie, dalsze zagęszczenie siatki ostatecznie powinno spowodować zbiegnięcie się tych dwóch wyników. Zgodnie bowiem z tym co napisano wcześniej, im mniejsze objętości użyte zostają do obliczeń, tym bardziej wszystkie metody powinny zbiegać się do wartości rzeczywistych.

Rysunek 17. Wpływ zagęszczenia siatki na wynik, w zależności od użytego schematu dyskretyzacji.

Podobny test wykonano na siatce hybrydowej, z przeważającą liczbą elementów hexahedralnych. Tutaj jednak, okazało się że wyniki są niemal identyczne, niezależnie od użytego schematu dyskretyzacji czy stopnia zagęszczenia siatki.

Rysunek 18. Rozkład temperatury na wyjściu dla siatki hexadralnej. Po lewej siatka zgrubna, po prawej gęsta.

Pokazuje to, dlaczego w świecie symulacji komputerowych siatki hexahdralne (a w szczególności strukturalne) są tak cenione. Siatka hexahedralna jest dużo mniej wrażliwa na dyfuzję numeryczną niż siatki tetrahedralne. Choć należy tu oczywiście zaznaczyć, że nie zawsze. Magiczny efekt uzyskamy jedynie wtedy, gdy dominujący kierunek przepływu będzie zgodny z ułożeniem siatki hexadralnej. Jeśli przepływ będzie realizowany na ukos elementów, efekt dyfuzji będzie spotęgowany.

Rysunek 19. Rozkład temperatury w przepływie nielepkim, realizowanym na ukos elementu hexahedralnego. Po prawej porównanie wyniku z siatką tetrahedralną.

Ogólny pogląd na występujące rozbieżności w konturze temperatury pomiędzy siatkami przedstawia rysunek poniżej. Przedstawiono na nim rozkład temperatury na wylocie przy użyciu schematu Second Order Upwind.

Rysunek 20. Rozkład temperatury na wylocie z domeny w zależności od siatki. Schemat 2nd order.

Wszystko da się podzielić

Jak widać, schematy dyskretyzacji jeśli użyte nieumiejętnie mogą znacznie wpłynąć na końcowy rezultat symulacji. Jest to wyjątkowo ważne przy siatkach tetrahedralnych, które z definicji nie są ułożone zgodnie z dominującym kierunkiem przepływu. W takim układzie dyfuzja numeryczna jest znaczna, więc istotne jest aby stosować schematy dyskretyzacji wyższego rzędu. W przypadku siatek hexahedralnych to zalecenie jest dalej w mocy, choć jak pokazało doświadczenie przy dobrym ułożeniu elementów w stosunku do przepływu problem się rozwiązuje. Dyfuzja numeryczna jest wtedy minimalna niezależnie od użytego schematu dyskretyzacji.

Na szczęście domyślne ustawienia Fluenta w kwestii dyskretyzacji są w większości przypadków wystarczające. Zwrócić należy jedynie uwagę na schematy użyte do modeli turbulencji, oraz większości modeli zaawansowanych (np. modele radiacyjne), które ustawione są na 1st Order Upwind.

Pozostaje oczywiście kwestia innych schematów, których tutaj nie omówiono, a które najczęściej dedykowane są dla konkretnego równania. Na przykład, schematy dyskretyzacji ciśnienia czy interfejsu międzyfazowego maja zupełnie inne działanie i właściwości. Moglibyśmy oczywiście poświęcić im kolejne strony tego artykułu, jednak obawiam się że surowa Pani Redaktor z naszego działu marketingu nie zezwoli na dzieło o objętości Władcy Pierścieni. Mam natomiast nadzieję, że przynajmniej te podstawowe informacje zawarte w artykule pozwolą Czytelnikowi na bardziej świadome i dokładniejsze ustawienia symulacji.

 AutorzyMaciej Kryś, MESco Sp. z o.o
  Natalia Zelek, Akademia Górniczo-Hutnicza