Numeryczne obliczenia cieplno-przepływowe CFD (Computational Fluid Dynamics) oferują ogromne możliwości przewidywania zjawisk i optymalizacji urządzeń i procesów, co sprawia, że są bardzo użyteczne w procesie rozwoju produktu. Pozwalają optymalizować projekt, skracać czas projektowania i redukować liczbę testów, prototypów czy też badań laboratoryjnych, co ma w efekcie wpływ na obniżenie kosztów.
Jednocześnie, jak przy każdym wyrafinowanym narzędziu, brak wiedzy i znajomości podstaw i metod, na których bazują nowoczesne kody do obliczeń CFD, może często prowadzić do sytuacji, w której zaawansowane narzędzie do modelowania zjawisk fizycznych, staje się w rękach niedoświadczonego użytkownika drogą zabawką do generowania kolorowych obrazków, niewiele mających wspólnego z rzeczywistością.
Poprawnie zaimplementowany algorytm numeryczny to taki, który jest zgodny i zbieżny. Mówiąc praktycznie – zgodność to nic innego, jak stwierdzenie i sprawdzenie, czy metoda numeryczna rozwiązuje rzeczywiście to równanie różniczkowe, jakie chcemy rozwiązać. Przerysowując problem – trzeba sprawdzić, czy zaimplementowana metoda numeryczna rzeczywiście rozwiązuje równanie przepływu, a nie np. elektromagnetyzmu. To zmartwienie producenta oprogramowania i stosunkowo banalny warunek. Jako użytkownicy możemy być spokojni, że tak jest. Druga rzecz to zbieżność. I tu dla nas będą już bardzo ważne wnioski. Otóż poprawnie zaimplementowana metoda numeryczna to taka, która wraz ze zmniejszaniem kroku czasowego i wielkości elementu siatki obliczeniowej będzie zbiegać do poprawnego rozwiązania równania, które wyliczamy. Zauważmy i podkreślmy, że nigdy nie ma gwarancji, że sensowne, bliskie rzeczywistości rozwiązanie uzyskujemy na dowolnie wybranej siatce. Uznaje się, że dobrej jakości symulacja CFD powinna mieć wykonane tzw. studium niezależności od siatki obliczeniowej. Co to znaczy? Wykonujemy naszą symulację na trzech różnych gęstościach siatki obliczeniowej i monitorujemy wyniki obliczeń. Na zgrubnej siatce wynik ten może odbiegać dość wyraźnie od pozostałych. Jeśli wyniki uzyskane na siatce o średniej gęstości różnią się niewiele z wynikami uzyskanymi na siatce o dużej gęstości, to możemy uznać, że znaleźliśmy taką gęstość siatki, która jest wystarczająca do uzyskania poprawnego wyniku. Siatki, które są rzadsze, są wówczas niewystarczające do poprawnego reprezentowania fizyki wszystkich zjawisk dziejących się w procesie.
Rysunki 1 i 2 obrazują prostą symulację laminarnego przepływu chłodnego płynu przez fragment wymiennika ciepła o gorących ściankach. Obliczenia zostały wykonane w programie QuickerSim CFD Toolbox for MATLAB na czterech różnych gęstościach siatki. Wszystkie cztery siatki, wraz z wykresem wizualizującym wartość mocy cieplnej odbieranej z tego wymiennika, w zależności od gęstości siatki, pokazano na rysunku 1.
Rysunek 2 obrazuje pole temperatury uzyskane na drugiej i czwartej siatce (licząc od najbardziej zgrubnej). Widać, że wynik ze zgrubnej siatki wyraźnie odbiega od pozostałych, natomiast w miarę zagęszczania siatki zaczynamy otrzymywać coraz bardziej poprawną, wiarygodną wartość odbieranego strumienia ciepła.
W praktyce inżynierskiej częste jest modelowanie zagadnień związanych z wymianą ciepła, jak również przepływami turbulentnymi. W obu przypadkach przy opływie ciał tworzą się dosyć cienkie termiczne lub hydrodynamiczne warstwy przyścienne, czyli obszary, gdzie na bardzo niedużej długości występują duże różnice temperatury (np. spadek temperatury w przestrzeni od gorącej ścianki do chłodnego płynu na bardzo krótkim odcinku – w niewielkiej odległości od ściany). Podobnie jest z hydrodynamiczną warstwą przyścienną, gdzie prędkość wzrasta, od zerowej na ścianie do dużej prędkości wewnątrz przepływu. Z tej fizyki zjawisk płyną wnioski dotyczące tego, jak poprawnie modelować taki przepływ. Wiemy, że w pobliżu ścianki będą występować duże gradienty poszczególnych wielkości przepływowych. Jednocześnie każdy kod numeryczny reprezentuje wartości tego pola w sposób dyskretny, a nie ciągły (przechowuje niewiadome albo w węzłach siatki, jak kody metody elementów skończonych, albo w środkach komórek, jak większość kodów metody objętości skończonych). Chcąc w sposób wiarygodny reprezentować tak strome gradienty pola przepływu w pobliżu ścianki musimy zadbać o to, by nasza siatka obliczeniowa dawała solwerowi szansę poprawnego reprezentowania takiego rozwiązania. Ilustruje to rysunek 3, gdzie pokazano bardzo rzadką siatkę obliczeniową, z nieadekwatnie dużymi elementami wokół opływanego ciała, oraz siatkę ze specjalnie wygenerowanymi dodatkowymi warstwami elementów w warstwie przyściennej. Na rzadkiej siatce solwer, najlepsze, co jest w stanie zrobić, to zadać zerową prędkość na ściance i dużo większą na węzłach siatki w pewnej odległości od ściany. Niemniej, przy tak dużych elementach implikuje to liniowe przybliżenie profilu prędkości wewnątrz elementu, które jest zupełnie nieprawdziwe i widać to w sposób bardzo przerysowany na „poszarpanym” polu prędkości. W przypadku siatki zagęszczonej w pobliżu ścianki solwer jest w stanie reprezentować zmiany prędkości (z dobrym przybliżeniem) w pobliżu ściany, i na kolejnych warstwach elementów następują one w sposób stopniowy – tym samym solwer jest w stanie sensownie, z dobrą dokładnością oddać szczegóły pola przepływu, co widać w prawej dolnej części rysunku 3.
Wnioski płynące z tego przykładu są jednak o wiele dalej idące. Na rysunku 3 zobaczyliśmy „ładnie” lub „brzydko” wyglądające pole prędkości. Niemniej, najbardziej istotne jest to, co się dzieje z gradientami. Są one albo poprawnie albo zupełnie błędnie przewidywane (nawet kilkukrotnie błędnie co do wartości). W obliczeniach przepływowych bardzo wiele ważnych wielkości interesujących dla inżyniera określonych jest gradientem jakiegoś pola: gęstość strumienia ciepła określona jest jako iloczyn przewodności cieplnej i gradientu pola temperatury:
Naprężenia lepkie na ściance określone są gradientem pola prędkości, jako:
Z kolei dla wielu przepływów bez dużych stref oderwania (na przykład przepływu w kanale, rurze itp.) naprężenia lepkie to jedyny, bądź główny składnik strat w przepływie i tym samym błędne reprezentowanie gradientów pola prędkości w symulacji w sposób bezpośredni przekłada się na zupełnie błędnie przewidziane spadki ciśnienia w urządzeniu. Z tych implikacji między poszczególnymi polami należy zdawać sobie sprawę, chcąc wykonywać wiarygodne obliczenia numeryczne. Rysunek 4 pokazuje ideę tego, jak zagęszczenie siatki w warstwie przyściennej wpływa na poprawne oddanie gradientów, a co za tym idzie – przewidywanie naprężeń lepkich, strumieni ciepła i podobnych wielkości określonych przez gradient jakiegoś pola.
Przywołajmy kolejny przykład związany z obliczeniami wymiennika ciepła. Jak wiemy, w fizyce istnieje prawo bilansu energii. Dla wymiennika ciepła możemy więc całkowity strumień ciepła przekazywany do przepływu liczyć na dwa różne sposoby. Możemy sprawdzić, ile mocy cieplnej zostało wniesione z chłodnym płynem wchodzącym do wymiennika i ile zostało wyniesione z rozgrzanym medium na wylocie. Różnica to oczywiście dostarczona moc. Można również to samo policzyć w inny sposób: bezpośrednio licząc gęstość strumienia ciepła w przepływie (a w szczególności na powierzchni gorących ścianek) z dobrze znanej zależności:
i tę gęstość strumienia ciepła scałkować po powierzchni ścianki, która dostarcza to ciepło. Z oczywistych względów jeden i drugi sposób powinien dawać ten sam wynik. Taka jest fizyka. A w symulacji numerycznej? Czy któryś sposób jest lepszy? Owszem. Mówiliśmy już, że mamy gwarancję zbieżności metody przy zmniejszaniu wielkości elementów w siatce obliczeniowej. A więc oczywiście jeden i drugi sposób liczenia strumienia ciepła w symulacji numerycznej dla odpowiednio dokładnej siatki da bardzo zbliżone wartości (różniące się od siebie o niewielkie błędy numeryczne). Niemniej drugi sposób, który wymaga policzenia gradientów pola temperatury jest potencjalnie o wiele bardziej wrażliwy na błędy w ich reprezentowaniu na złej siatce. Tym samym, ten sposób jest o wiele bardziej podatny na błędy i wprawdzie w granicy daje on w pełni poprawny wynik i nie można nic zarzucić poprawności solwera, natomiast można mieć zastrzeżenia do użytkownika, jeśli nie zdaje sobie sprawy z tych niebezpieczeństw i nie zadba zawsze o sprawdzenie bilansów podstawowych wielkości, które powinny być zachowane. Rysunek 5 przedstawia, jak dla czterech różnych gęstości siatek obliczeniowych zmienia się wartość strumienia ciepła odbieranego z wymiennika, gdy liczona jest jednym bądź drugim sposobem. Jednocześnie widać jedyną gwarantowaną w numeryce cechę: wraz z poprawą siatki oba wyniki zgadzają się coraz dokładniej, niemniej jeden sposób sprawdzania osiągów wymiennika jest bardziej wiarygodny niż drugi (choć oba poprawne).
Każdy z zaprezentowanych tu przypadków Czytelnik może samodzielnie zasymulować i przetestować w używanym przez siebie programie CFD, sprawdzić wpływ każdego z omówionych parametrów i zobrazować dla własnej wiedzy, jak duże zmiany wyników symulacji mogą powodować zmiany w siatce obliczeniowej i poszczególnych ustawieniach solwera. Dużo więcej uwagi i wiedzy wymaga poprawne i wiarygodne modelowanie przepływów turbulentnych.
W praktyce przemysłowej, duża część modelowanych przepływów, dotyczy właśnie tych, w których występuje turbulencja. Do modelowania turbulencji służy obecnie kilka różnych podejść. Najpowszechniejszym z nich, szeroko stosowanym w praktyce przemysłowej jest technika RANS (Reynolds Averaged Navier-Stokes). Jest to technika bazująca na uśrednieniu w czasie równań Naviera-Stokesa opisujących dynamikę płynu. Zakłada, że możemy turbulencję (która z natury rzeczy jest procesem niestacjonarnym) reprezentować jako uśrednioną w czasie i wszystkie efekty turbulentnego mieszania modelować w tymże uśrednionym (wtedy już stacjonarnym) polu przepływu. Technika ta jest niemal jedyną szeroko stosowaną, gdyż dwie pozostałe (LES i DNS) są o wiele bardziej kosztowne obliczeniowo.
Każdy program do obliczeń CFD bazujący na technice RANS ma możliwość uruchomienia i wyboru jednego z istniejących modeli turbulencji. W poprzedniej części artykułu dużo uwagi poświęciliśmy zagadnieniom tworzenia siatek obliczeniowych mogących wiernie reprezentować przepływ. Oczywiste jest więc już dla nas to, jak należy zagęścić siatkę w pobliżu ścianek materialnych. Co jednak z innymi ustawieniami? Jak warunki brzegowe wybrane przez użytkownika mogą wpłynąć na wyniki obliczeń? Czasem, ze względu na brak dokładnych danych, bywa, że użytkownik zadaje wartości warunków brzegowych w sposób nie do końca przemyślany. Kolejne pytanie dotyczy samej fizyki zjawiska. Przyjrzyjmy się tu dokładniej przykładowi wziętemu z jednej z dużych firm – modelowano pozornie prosty przepływ płynu wraz z wymianą ciepła w rurze. Prędkość przepływu i wymiary liniowe instalacji były takie, że liczba Reynoldsa była na poziomie 2000. W firmie porównywano wyniki uzyskiwane dwoma różnymi programami CFD: Ansys Fluent i QuickerSim. W pierwszych testach wartości temperatury płynu na wylocie z instalacji różniły się o kilkadziesiąt stopni między obydwoma programami. Okazało się, że błąd wynikał z nieświadomości użytkownika. Liczba Reynoldsa na poziomie 2000 oznacza przepływ laminarny w rurze. Tam nie ma turbulencji. Użytkownik musi być świadomy, że włączenie modelowania turbulencji (jak zrobił we Fluencie) oznacza rozwiązywanie dodatkowych równań. Z nich wyznaczana jest lepkość turbulentna. Na podstawie lepkości turbulentnej wyznaczana jest turbulentna dyfuzyjność termiczna. Modele typu RANS modelują ją jako iloraz lepkości turbulentnej i turbulentnej liczby Prandtla. I tu użytkownik musi być świadomy tego, jakie cechy ma wybrany przez niego model turbulencji i co rozwiązuje użytkowany program. Warto podkreślić, że modele turbulencji nie są ogólne. Nie jest tak, że w reżimie przepływu turbulentnego modelują turbulencję, a dla małych liczb Reynoldsa są w stanie samoistnie przewidzieć, że turbulencja nie wystąpi. Nie, modele turbulencji nie są mądre i nie odzwierciedlają fizyki. One sensownie przybliżają rzeczywistość w przypadkach, dla których zostały skalibrowane. W szczególności nie upraszczają się do przepływu laminarnego wtedy, gdy liczba Reynoldsa jest dostatecznie mała. Taką cechę wykazywałyby może dużo bardziej zaawansowane modele przejścia turbulentnego, umiejące również przewidywać relaminaryzację. Z tych jednak mało kto korzysta, bo potrzebne są tylko do bardzo specyficznych analiz i w zdecydowanej większości przypadków nie ma potrzeby ich wykorzystania. Przyjrzyjmy się zatem, z jakim błędem musimy się liczyć w przypadku wykonania obliczeń przepływu z wymianą ciepła na rozgrzanej, płaskiej płytce, jeśli przepływ z niską liczbą Reynoldsa będziemy modelować z wykorzystaniem modelu turbulencji. Rysunek 1 pokazuje poprawny obraz pola temperatury w takim przepływie, zaś rysunek 2 ilustruje rozkład lokalnej liczby Nusselta wzdłuż płytki (liczba Nusselta wyraża w bezwymiarowej formie wartość współczynnika przejmowania ciepła od płytki). Zauważalna jest olbrzymia niezgodność między wynikami obliczeń i rozwiązaniem teoretycznym.
Widoczna duża rozbieżność wyniku symulacji i rozwiązania analitycznego wywołana nieuzasadnionym włączeniem modelu turbulencji. Wykonajmy teraz identyczne obliczenia bez włączania modelu turbulencji (z fizyki zjawiska wiemy, że przepływ powinien być laminarny). Porównanie dla tego przypadku przebiegu liczby Nusselta wyznaczonej teoretycznie i numerycznie ilustruje rysunek 3. Tym razem widać już bardzo dużą zgodność z błędem względnym na poziomie pojedynczych procentów. Dokładnie to samo zaniedbanie było przyczyną błędu w obliczeniach firmy, o której wspominałem wcześniej. Obliczenia w programie Ansys Fluent użytkownik wykonał z włączonym modelem turbulencji k–ε, który w niefizyczny sposób podniósł mocno efektywną dyfuzyjność termiczną i tym samym bardzo mocno obniżył temperatury w przepływie. Po wyłączeniu modelu turbulencji, Fluent przewidywał już poprawnie pole temperatury – najwyższe wartości w przepływie wzrosły o kilkadziesiąt stopni i z dokładnością do pojedynczych stopni zgadzały się z przewidywaniami programu QuickerSim. Dość oczywisty, ale ważny wniosek płynie stąd taki, że to użytkownik z całą swą świadomością i znajomością metod CFD musi zadbać o wiarygodność symulacji. Program jest jedynie narzędziem do rozwiązania równań zadanych przez użytkownika.
Badając turbulencję przy wlotach urządzeń należy pamiętać, że o ile prędkość płynu na wlocie jest na tyle duża, że przepływ jest turbulentny, to w pobliżu takiego wlotu tworzą się strefy ścinania. Załóżmy, że interesuje nas wyznaczenie intensywności turbulencji w przepływie. Intensywność turbulencji definiowana jest jako stosunek wartości fluktuacji turbulentnych pola prędkości do średniej prędkości przepływu. Wyrażana jest najczęściej w procentach. W symulacji numerycznej może zostać policzona ze wzoru:
gdzie ti oznacza intensywność turbulencji, k energię kinetyczną turbulencji, zaś U prędkość referencyjną. Chcąc dokładnie wyznaczyć intensywność turbulencji, trzeba zatem dokładnie rozwiązać pole energii kinetycznej turbulencji. Właściwe do tego równanie, to równanie transportu (konwekcji-dyfuzji) z członem produkcyjnym postaci:
gdzie νt to lepkość turbulentna, a S oznacza tensor prędkości odkształcenia – ponownie mamy do czynienia z liczeniem wyrażenia, które zależy od gradientów pola prędkości. Tym razem jednak już nie przy ściance, a w warstwie ścinanej, na granicy wpływającego strumienia i niemal stojącego płynu. Poprawne wyznaczenie tych gradientów wymaga zastosowania zagęszczeń siatki w takiej strefie. Rysunki 4 i 5 ilustrują pole intensywności turbulencji, wyznaczone na siatce mającej takie zagęszczenie (Rys. 4) i bez tego zagęszczenia (Rys. 5). Widać, że w drugim przypadku pole intensywności turbulencji jest delikatnie niedoszacowane.
Modelowanie CFD jest zagadnieniem o bardzo złożonej sieci zależności. Najważniejszy w tym wszystkim jest użytkownik, który musi wzajemne zależności rozumieć, by móc poprawnie reagować na błędy i nieścisłości, i doprowadzić do zgodności obliczeń numerycznych z rzeczywistością. W bardziej zaawansowanych przypadkach powinien on wiedzieć, dlaczego stosowane modele tej rzeczywistości nie są w stanie oddać. Używane oprogramowanie musi dawać użytkownikowi kontrolę nad bogactwem modeli i siatką obliczeniową. Nie powinno być tylko „czarną skrzynką”, zabawką dającą bardzo zgrubne oszacowania lub stanowiącą tylko źródło kolorowych obrazków do folderów reklamowych, gdyż czułość niektórych obliczeń na różne ustawienia może być tak duża, że trudno mówić nawet o zgrubnych oszacowaniach. Dokładność w CFD wymaga o wiele więcej pracy i poświęcenia uwagi szczegółom.
Bartosz Górecki
b.gorecki@quickersim.com
artykuł ukazał się w dwóch częściach w wydaniach 11 (110) listopad i 12 (111) grudzień 2016
