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.
