Przemysław Błasiak1
Streszczenie: W artykule przedstawiono problematyk˛e matematycznego modelowania zjawiska konwekcji. Przeprowadzono symulacj˛e numeryczn ˛a wymiany ciepła w warunkach konwekcji wymuszonej. Do oblicze´n wykorzystano standardowy model turbulencji SST oraz podstawowy model du˙zych wirów Large Eddy Simulation (LES) – model Smagorinskiego.
Słowa kluczowe: ANSYS CFX, konwekcja wymuszona, Large Eddy Simulation
1. Wprowadzenie
Przez konwekcj˛e rozumie si˛e proces transportu ciepła spowodowany makroskopowym ruchem płynu. Ruch ten mo˙ze by´c wymuszony, np. przez pomp˛e lub urz ˛adzenie mieszaj ˛ace, lub powsta´c w wyniku istnienia ró˙znic g˛esto´sci w obszarze płynu. W pierwszym przypadku ruch płynu wywołany jest siłami zewn˛etrznymi i nazywany jest konwekcj ˛a wymuszon ˛a. Przypadek drugi nazywany jest z kolei konwekcj ˛a naturaln ˛a lub swobodn ˛a. Siłami odpowiedzialnymi za ruch płynu s ˛a tutaj siły masowe, które wyst˛epuj ˛a jako rezultat ró˙znic g˛esto´sci płynu, powstałych w wyniku istnienia gradientów temperatury w polu przepływu. Siły masowe s ˛a w rzeczywisto´sci generowane przez gradient ci´snienia w płynie i w literaturze okre´slane s ˛a jako siły wyporu (ang. buoyancy forces). Termin konwekcja wymuszona odnosi si˛e do sytuacji, gdzie siły wyporu s ˛a zaniedbywalnie małe. W przypadku zagadnie´n, w których działaj ˛a siły zewn˛etrzne, a efektów działania sił masowych nie mo˙zna pomin ˛a´c, mówi si˛e o tzw. konwekcji mieszanej.
Przepływy konwekcyjne dzieli si˛e ponadto na zewn˛etrzne oraz wewn˛etrzne. Przepływ zewn˛etrzny dotyczy zwykle opływów ciał. Przypadek przepływu wewn˛etrznego to np. przepływ w rurze. Przy opływach ciał wyst˛epuj ˛a tzw. wielko´sci przepływu niezaburzonego np. temperatura T∞. Jako ich warto´sci przyjmuje si˛e zwykle warto´sci wyst˛epuj ˛ace w obszarze przed opływanym ciałem. Dla przepływów wewn˛etrznych bierze si˛e odpowiedni ˛a ´sredni ˛a warto´s´c temperatury Tb. Definicja ´sredniej temperatury jest wa˙zna przy obliczaniu współczynnika wnikania ciepła α.
Dla przypadku grzania mo˙zna go obliczy´c ze wzoru α = − λf
Tb− Tw
∂T
∂n
x=0
(1) gdzie λf to współczynnik przewodzenia ciepła płynu, Tw temperatura ´scianki, ∂T∂n
x=0 gradient temperatury płynu na powierzchni ´scianki.
W prezentowanej pracy podj˛eto prób˛e numerycznego modelowania wymiany ciepła w warunkach konwekcji wymuszonej w przypadku przepływu wewn˛etrznego. Szczegółowy opis postawionego zadania znajduje si˛e w punkcie 3. Do przeprowadzenia symulacji numerycznych wykorzystano komercyjne oprogramowanie CFD (Computational Fluid Dynamic) ANSYS CFX 13.0.
1Instytut Techniki Cieplnej i Mechaniki Płynów, Politechnika Wrocławska, Wybrze˙ze Wyspia´nskiego 27 50-370 Wrocław, przemyslaw.blasiak@pwr.wroc.pl
194
2. Modelowanie zjawiska konwekcji
2.1. Równania ruchu płynu i transportu energii cieplnej
Przy rozwi ˛azywaniu zagadnie´n przepływów z jednoczesn ˛a konwekcyjn ˛a wymian ˛a ciepła nale˙zy rozwi ˛aza´c układ równa´n ró˙zniczkowych cz ˛astkowych rz ˛adz ˛acych transportem masy, p˛edu oraz energii. Równania te wyra˙zaj ˛a tzw. zasady zachowania. W skład układu wchodz ˛a [1]:
1) równanie ci ˛agło´sci wyra˙zaj ˛ace prawo zachowania masy
∂ρ
∂t + ∇ · (ρu) = 0 (2)
gdzie ρ to g˛esto´s´c płynu, a u = (u, v, w) wektor pr˛edko´sci. Pierwszy wyraz oznacza pr˛edko´s´c zmiany g˛esto´sci w obj˛eto´sci kontrolnej na jednostk˛e obj˛eto´sci. Drugi natomiast pr˛edko´s´c strumienia masy przechodz ˛acego przez powierzchni˛e obj˛eto´sci kontrolnej na jednostk˛e obj˛eto´sci.
2) równanie Naviera-Stokes’a wyra˙zaj ˛ace prawo zachowania p˛edu ρ ∂u gdzie f jest jednostkowym wektorem sił masowych, µ współczynnikiem lepko´sci dynamicznej płynu, δij delt ˛a Kroneckera (i, j, k = 1, 2, 3). Składniki znajduj ˛ace si˛e w nawiasie kwadratowym reprezentuj ˛a tzw. lepki tensor napr˛e˙ze´n oznaczany w literaturze jako τij [2].
3) równanie energii, które stanowi zasad˛e zachowania energii ρ ∂e
∂t + u · ∇e
+ p (∇ · u) = ∂Q
∂t − ∇ · q + Φ (4)
Pierwszy wyraz po prawej stronie oznacza tempo wydzielania si˛e ciepła na jednostk˛e obj˛eto´sci z wewn˛etrznych ´zródeł ciepła, natomiast drugi wyraz to strumie´n ciepła przekazywanego przez przewodzenie. Składnik Φ oznacza tzw. funkcj˛e dyssypacji. Reprezentuje ona ciepło jakie zostaje wydzielone w procesie zamiany energii mechanicznej na energi˛e ciepln ˛a podczas deformacji płynu, spowodowanej siłami lepko´sci.
W celu domkni˛ecia powy˙zszego układu równa´n, nale˙zy go uzupełni´c o równania ustalaj ˛ace relacje pomi˛edzy poszczególnymi wła´sciwo´sciami płynu takimi jak: ρ, T, p, e. Przykładem jest równanie stanu, które wi ˛a˙ze ze sob ˛a zmienne ci´snienia p, g˛esto´sci ρ oraz temperatury T . Przy zało˙zeniu, ˙ze płyn mo˙zna traktowa´c jak gaz doskonały, obowi ˛azuje równanie stanu Clapeyron’a (5). Nale˙zy wtedy doda´c jeszcze równanie na energi˛e wewn˛etrzn ˛a gazu doskonałego (6):
p = ρ R T (5)
e = cυT (6)
gdzie R jest indywidualn ˛a stał ˛a gazow ˛a, a cυ ciepłem wła´sciwym przy stałej obj˛eto´sci.
2.2. Problem turbulencji
Ruch turbulentny jest wirowy i odznacza si˛e bardzo skomplikowan ˛a struktur ˛a. Bezpo´srednie rozwi ˛azanie układu równa´n rz ˛adz ˛acych przepływem mo˙zliwe jest jedynie na drodze numerycznej, ale pod warunkiem wykorzystania bardzo dokładnej siatki numerycznej. Jest to tzw. metoda DNS (Direct Numerical Simulation). Zwykle obliczenia t ˛a metod ˛a wymagaj ˛a tak du˙zej liczby oczek siatki numerycznej, ˙ze przy obecnych mo˙zliwo´sciach komputerów, czas symulacji jest zbyt długi i ich prowadzenie staje si˛e niepraktyczne. W celu rozwi ˛azania tego problemu przeprowadza si˛e
195
u´srednianie równa´n rz ˛adz ˛acych. S ˛a dwa główne podej´scia: u´srednianie czasowe oraz u´srednianie przestrzenne.
Pierwsza metoda prowadzi do tzw. równa´n Reynoldsa dla przepływu turbulentnego.
W literaturze okre´slana jest jako RANS (Reynolds Averaged Navier-Stokes equations). Procedura u´sredniania czasowego prowadzi do wygenerowania tzw. tensora napr˛e˙ze´n Reynoldsa. Pojawiaj ˛a si˛e dodatkowe wyrazy w równaniach i układ równa´n staje si˛e niedomkni˛ety (tzw. clossure problem).
Te dodatkowe wyrazy modelowane s ˛a za pomoc ˛a klasycznych modeli turbulencji, z których najlepiej znane s ˛a: k − ε oraz model napr˛e˙ze´n Reynoldsa [2]. Moce obliczeniowe wymagane w tym podej´sciu, pozwalaj ˛ace uzyska´c wyniki na zadowalaj ˛acym poziomie dokładno´sci, s ˛a umiarkowane.
Metoda u´sredniania przestrzennego, znana w literaturze jako LES (Large Eddy Simulation) jest form ˛a po´sredni ˛a oblicze´n turbulencji pomi˛edzy metod ˛a RANS i DNS. Daje ona wyniki dokładniejsze od metody RANS, przy wymaganej wi˛ekszej mocy obliczeniowej. Obecnie metoda ta jest coraz cz˛e´sciej wykorzystywana w obliczeniach nie tylko naukowych, ale i przemysłowych.
Zostanie ona bardziej szczegółowo opisana w kolejnym podpunkcie.
2.3. Large Eddy Simulation - LES
W podej´sciu LES wi˛eksze, trójwymiarowe, nieustalone wiry turbulentne s ˛a bezpo´srednio obliczane, podczas gdy mniejsze struktury wirowe s ˛a modelowane. Mo˙zna wyró˙zni´c cztery etapy [3]:
1. Filtrowanie - dekompozycja pr˛edko´sci u(u, v, w, t) na sum˛e dwóch składowych: filtrowanej
¯
u oraz resztkowej u0. Składowa filtrowana jest rozwi ˛azywana, natomiast resztkowa modelowana. Składowa u0 nazywana jest w literaturze jako subgrid − scale (SGS).
Przefiltrowane pole pr˛edko´sci ¯u, które jest trójwymiarowe i zale˙zne od czasu, reprezentuje ruch du˙zych struktur wirowych. Filtrowanie jest operacj ˛a matematyczn ˛a, która usuwa z całego pola pr˛edko´sci u składowe mniejsze od wielko´sci u˙zytego filtru - usuni˛ecie flaktuacji i pozostawienie jedynie wi˛ekszych wirów.
2. Utworzenie tzw. przefiltrowanych równa´n rz ˛adz ˛acych. W równaniach p˛edu pojawia si˛e tensor napr˛e˙ze´n resztkowych τrij (SGS stress tensor).
3. Wybranie modelu dla tensora SGS.
4. Rozwi ˛azanie układu równa´n numerycznie.
W celu domkni˛ecia układu równa´n przefiltrowanej pr˛edko´sci nale˙zy zamodelowa´c tensor napr˛e˙ze´n resztkowych τij. Pierwszym i najprostszym modelem jest tzw. model Smagorinskiego.
W modelu tym zakłada si˛e, ˙ze napr˛e˙zenia resztkowe zale˙z ˛a liniowo od przefiltrowanego tensora deformacji pr˛edko´sci ¯Sij Współczynnik proporcjonalno´sci νr(x, t) to tzw. lepko´s´c wirowa ruchu resztkowego. Przez analogi˛e do drogi mieszania Prandtla współczynnik ten okre´sla si˛e jako
νr = lS2 S = (C¯ S∆)2S¯ (9)
gdzie ¯S to charakterystyczny przefiltrowany tensor deformacji
196
S = 2 ¯¯ SijS¯ij12
(10) lS skala długo´sci Smagorinskiego, CS współczynnik Smagorinskiego i ∆ szeroko´s´c u˙zytego filtru.
3. Symulacja numeryczna
Rys. 1. a) Schemat zadania oraz warunki brzegowe, b) domena obliczeniowa, c), d) szczegóły siatki numerycznej
Schemat wraz z warunkami brzegowymi postawionego problemu wymiany ciepła w warunkach konwekcji wymuszonej oraz szczegóły siatki numerycznej, pokazane s ˛a na Rys. 1.
W cylindrze o ´srednicy wewn˛etrznej 140 mm i wysoko´sci 80 mm, wykonanym z aluminium, znajdowało si˛e powietrze, które zostało zamodelowane jako gaz doskonały. Grubo´s´c ´scianki cylindra wynosiła 5 mm. Wła´sciwo´sci materiałowe powietrza i aluminium zostały przedstawione w Tab. 1.
W osi cylindra umieszczony był obrotowy wał o ´srednicy 25 mm, na którym zamontowane były dwie łopatki o grubo´sci 3 mm, poruszaj ˛ace si˛e z pr˛edko´sci ˛a 10 obr/min. Ich zadaniem było wymuszenie przepływu wirowego gazu oraz zdzieranie warstwy przy´sciennej (boundary layer).
Warstwa ta stanowi główny opór cieplny w procesie wymiany ciepła. Mechaniczne usuni˛ecie warstwy granicznej powinno skutkowa´c zwi˛ekszeniem warto´sci współczynnika wnikania ciepła α, czyli polepszeniem transferu ciepła od gazu do ´scianki cylindra. Na wzrost wnikania ciepła ma tak˙ze wpływ szczelina pomi˛edzy ko´ncem łopatki a ´sciank ˛a, która wynosiła 0,5 mm.
Na powierzchni wału został nało˙zony warunek brzegowy Neumanna w postaci stałej g˛esto´sci strumienia ciepła ˙q = 7500 W/m2. Na ´sciance zewn˛etrznej cylindra nało˙zono warunek brzegowy Dirichleta stałej temperatury Tw = 12◦C. Powierzchnie górna i dolna cylindra były adiabatyczne.
Na powierzchniach wszystkich ´scianek nało˙zono warunek zerowej pr˛edko´sci, tzw. no-slip condition.
W chwili pocz ˛atkowej płyn był nieruchomy i miał temperatur˛e T = 12◦C.
197
Czas symulacji wynosił 25 s. Krok czasowy ∆t = 0.017 s, co odpowiada obrotowi łopatek o 1◦. Siatka numeryczna składała si˛e z około 6.5 miliona elementów typu heksagonalnego.
Do oblicze´n wykorzystano dwa modele turbulencji. Klasyczny model RANS, tzw. SST [2] oraz model LES Smagorinskiego. Stał ˛a Smagorinskiego ustawiono jako CS = 0, 1. W równaniach pomini˛eto człony od sił grawitacyjnych oraz funkcj˛e dyssypacji. Powietrze traktowano jako gaz nie´sci´sliwy, co jest uzasadnione małymi pr˛edko´sciami przepływu dochodz ˛acymi maksymalnie do 0,083 m/s.
Tab. 1. Wła´sciwo´sci fizyczne powietrza i aluminium u˙zyte w obliczeniach
Wła´sciwo´s´c fizyczna Powietrze Aluminium Przewodno´s´c cieplna λ, W/(m K) 0,0261 237
G˛esto´s´c ρ, kg/m3 1,185 2702
Ciepło wła´sciwe cp, J/(kg K) 1004,4 903 Lepko´s´c dynamiczna η, kg/(m s) 1, 831 · 10−5 –
4. Wyniki
Rys. 2. Rozkład temperatury w pobli˙zu ostrza łopatki w płaszczy´znie z = 40 mm po czasie t = 25 s: model SST (z lewej) oraz model LES Smagorinskiego (z prawej)
Rys. 3. Rozkład składowej stycznej pr˛edko´sci w płaszczy´znie z = 40 mm po czasie t = 25 s: model SST (z lewej) oraz model LES Smagorinskiego (z prawej)
198
Na Rys. 2. i Rys. 3. przedstawiono odpowiednio rozkłady temperatury i pr˛edko´sci w płaszczy´znie z = 40 mm, po czasie t = 25 s, obliczone za pomoc ˛a modelu SST oraz Smagorinskiego.
Mo˙zna zauwa˙zy´c, ˙ze najwi˛eksze gradienty temperatury oraz pr˛edko´sci wyst˛epuj ˛a w pobli˙zu ko´nca poruszaj ˛acej si˛e łopatki. Wyniki ró˙zni ˛a si˛e nieznacznie. Ponadto obliczono warto´sci współczynnika wnikania ciepła α ze wzoru (1). Jego warto´sci wyniosły αSST = 6, 97 W/(m2K) oraz αLES = 6, 62 W/(m2K). S ˛a one zani˙zone, poniewa˙z dla konwekcji naturalnej powinny mie´sci´c si˛e w zakresie 1 − 20 W/(m2K), a dla wymuszonej 20 − 100 W/(m2K).
5. Podsumowanie
W artykule przedstawiono problematyk˛e numerycznego modelowania konwekcji wymuszonej. Opisano podstawowy model Large Eddy Simulation - model Smagorinskiego.
Ponadto przeprowadzono symulacj˛e numeryczn ˛a konwekcyjnej wymiany ciepła. Obliczenia przeprowadzono z wykorzystaniem modelu turbulencji SST oraz Smagorinskiego. Obliczenia wykonano w trybie równoległym na 30 procesorach i trwały one trzy tygodnie. Zaprezentowano wyniki dla czasu symulacji t = 25 s. Niestety wyniki s ˛a jedynie jako´sciowe. Wida´c z nich, ˙ze najwi˛eksze zmiany temperatury i pr˛edko´sci wyst˛epuj ˛a w pobli˙zu ostrza łopatki. Obliczone warto´sci współczynnika wnikania ciepła α dla modelu SST i LES s ˛a zani˙zone. Konieczne jest udoskonalenie modelu. Najwi˛ekszy wpływ na bł˛edy mogła mie´c u˙zyta siatka oraz nieuwzgl˛ednienie zmian wła´sciwo´sci fizycznych powietrza z temperatur ˛a. Konsekwencj ˛a było niedokładne rozwi ˛azanie przepływu w obszarze warstwy przy´sciennej, co w procesie wymiany ciepła stanowi kluczowe zagadnienie i ma najwi˛ekszy wpływ na wyniki.
Literatura
[1] D.A. Anderson, J.C. Tannehill, R.H. Pletcher: Computational Fluid Mechanics and Heat Transfer, McGraw-Hill 1984.
[2] H.K. Versteeg, W. Malalasekera: An Intorduction to Computational Fluid Dynamics, 2nd ed., Pearson Education Limited 2007.
[3] S.B. Pope: Turbulent Flows, Cambridge University Press 2000.
„Obliczenia wykonano na komputerach Wrocławskiego Centrum Sieciowo-Superkomputerowego (http://www.wcss.wroc.pl)”
NUMERICAL MODELING OF FORCED TURBULENT HEAT CONVECTION
The article deal with numerical modeling of heat transfer in forced convection conditions.
Computational domain was consisted of a cylinder filled with a gas. Inside the cylinder there were two rotating scraper blades. Scrapers were mounted on a cylindrical shaft. Their task was stripping the boundary layer. Boundary layer is the main thermal resistance of the heat exchange process.
Mechanical removal of the fluid adhering to the walls is a method of heat transfer intensification.
The heat supplied to the fluid was modeled as a constant heat flux boundary condition on the shaft surface. On the outer surface of the cylinder the cooling jacket was modeled with the use of the constant temperature boundary condition. In the simulation classical turbulence model SST was compared with more complex model of Large Eddy Simulation (LES), so-called Smagorinsky model. Calculations were perormed on the supercomputers on 30 processors. Commercial CFD (Computational Fluid Dynamics) software ANSYS CFX 13.0 was employed.
199