7. Badania numeryczne kolektorów dolotowych
7.1. Założenia teoretyczne do modelu
Zadaniem dynamiki płynów jest analiza ruchu dużej liczby indywidualnych cząstek, które mają na siebie wzajemny wpływ. W przypadku płynów są to atomy lub cząsteczki. Podstawowym założeniem jest ośrodek ciągły czyli gęstość płynu jest wystarczająco duża. Oznacza to, że najmniejsza badana objętość zawiera wystarczającą liczbę cząsteczek, dla których można obliczyć prędkość i średnią energię kinetyczną.
W tym przypadku możliwe jest zdefiniowanie prędkości, ciśnienia, temperatury, gęstości i innych parametrów dla każdej objętości kontrolnej. Podstawowe równania dynamiki płynów można wyprowadzić bazując na trzech zasadach zachowania (bilansach): masy, pędu i energii oraz wiedzy o lokalnych właściwościach przepływu. Wielkości bilansowe przyłożone są zawsze w środku objętości kontrolnej. Ich zmiana w czasie może być wywołana albo rozpływem strumieni przyłożonych na ściankach, albo źródłami.
Podejście lokalne powstało na celu opisania trójwymiarowego rozkładu pól prędkości, ciśnienia, temperatur, gęstości w przepływającym płynie. Obejmuje ono takie zjawiska jak np. kawitacja wody, przejścia laminarno-turbulentne, oderwania strumienia itp.
Bilans masy oznacza, że suma strumieni mas wpływająca do obiektu technicznego jest równa sumie strumieni mas wypływających. W ujęciu lokalnym równanie zachowania masy w postaci zachowawczej (numerycznej) wyraża się następująco [85, 88, 90, 93, 97]:
∂tρ + div(ρv⃗ ) = 0 (7.1)
gdzie: ∂t – zmiana w czasie, div – dywergencja strumienia masy, pędu i energii po wszystkich powierzchniach objętości cząstki płynu, ρ – gęstość masy [kg/m3], pv⃗ – gęstość pędu [(kgm/s)/m3], v⃗ – prędkość płynu, t – czas.
W układzie kartezjańskim równanie bilansu masy przyjmuje postać:
∂
∂t(ρ) + ∂
∂xi(ρvi) = 0; i, k = x, y, z (7.2) gdzie: x, y, z – współrzędne układu prostokątnego.
Bilans pędu określa, że suma sił dynamicznych i sił naporu płynu wpływającego do obiektu technicznego pomniejszona o siły reakcji ścianek i siły ciężkości jest równa sumie sił dynamicznych i naporu płynu wypływającego z obiektu. Bilans pędu jest
„kontynualną” wersją II zasady dynamiki Newtona. W postaci zachowawczej równanie zachowania pędu przedstawia się następująco [33, 35, 36, 37]:
∂t(ρv⃗ ) + div(ρv⃗ ⨂v⃗ + pI⃡) = div(τ⃡c) + ρb⃗ (7.3)
gdzie: ρv⃗ ⨂v⃗ – konwekcyjny strumień pędu [N/m2], pI⃡ – sprężysty strumień pędu [N/m2], τ⃡c – dyfuzyjny strumień pędu [N/m2], ρb⃗ – źródło pędu [N/m3], τ⃡c – całkowite lepkie naprężenie.
W układzie kartezjańskim bilans pędu jest równy:
∂
∂t(ρvi) + ∂
∂xk(ρvivk+ pδik) = ∂
∂xk(τikc ) + ρbi (7.4) Bilans energii określa, że strumień energii dysponowanej wpływający do obiektu technicznego powiększony przez strumień energii cieplnej (ciepło) i strumień energii mechanicznej (pracę) pobrany w obiekcie jest równy strumieniowi energii odpadowej wypływającej z niego. Bilans energii jest I zasadą termodynamiki. W ujęciu lokalnym równanie w postaci zachowawczej (numerycznej) wyraża się następująco:
∂t(ρe) + div (ρ [e +p
W układzie kartezjańskim bilans energii przyjmuje postać:
∂
∂t(ρe) + ∂
∂xi(ρ[e +p
ρ]vi) = ∂
∂xi(τikc vk+ qic) + ρbivi (7.6) Zachowanie pewnej wielkości przepływu oznacza, że całkowita zmiana zachodząca wewnątrz określonej objętości może być zdefiniowana jako efekt wypadkowy działania objętości płynu transportowanej przez granicę tej objętości, z wszelkich sił wewnętrznych i źródeł, oraz sił zewnętrznych działających na tą objętość. Wielkość objętości przemieszczającej się poza granicę objętości kontrolnej nazywa się strumieniem. Strumień może być rozbity na dwie różne składowe: wywołaną przez konwekcję i przez ruch cząstek w płynie podczas spoczynku. Drugi składnik ma naturę dyfuzyjną – jest proporcjonalny do gradientu rozważanej wielkości i stąd jest pomijalny przy rozkładzie homogenicznym. Dodatkowo przy rozważaniu bilansów masy, pędu i energii konieczne jest zdefiniowanie tzw. skończonej objętości kontrolnej.
Dowolny skończony obszar przepływu, zamknięty granicą dΩ i umiejscowiony na stałe w przestrzeni nazywa się objętością kontrolą Ω (rys. 7.1). Na powierzchni objętości zdefiniowane jest elementarną powierzchnię dS i związane z nią normalny wektor n⃗ .
Rys. 7.1. Definicja skończonej objętości kontrolnej [13]
Ogólne prawo zachowania stanowi, że dla objętości kontrolnej prawo zachowania stosowane do przykładowej wielkości skalarnej na jednostkę objętości – U określa, że jej zmienność w czasie jest równa sumie udziału przepływu konwekcyjnego – ilości U przechodzącej przez granicę objętości kontrolnej z prędkością v; następnie dyfuzji strumienia – zgodnie z prawem Ficka i ostatecznie ze źródeł objętości jak i powierzchni.
Prawo zachowania dla wielkości skalarnej U przedstawia się następująco:
∂
Jeżeli wielkość U byłaby wektorowa, równanie (7.7) przyjmuje postać:
∂ źródeł objętości,Q̿̿̿̿ – tensor źródeł powierzchniowych. S
Równanie (7.8) jest formą zapisu całkowego prawa zachowania i według niej wykonywane są obliczenia w oprogramowaniu CFD. Ta forma ma dwie ważne cechy [11, 15, 67, 78, 87]:
a) jeżeli nie ma źródeł objętości, to zmiana wielkości U zależy wyłącznie od strumienia przepływającego przez granicę objętości kontrolnej a nie zależy od strumienia wewnątrz objętości,
b) prawo zachowania jest obowiązujące także w przypadku pojawienia się nieciągłości w przepływie płynu takich jak fala uderzeniowa lub nieciągłość styku.
Z powyższych ogólnych równań zachowania można wyprowadzić odpowiednie równania trzech praw zachowania w dynamice płynów: masy, pędu i energii. Równanie zachowania masy dla płynów jednofazowych oznacza, że w rozpatrywanej objętości płynu masa nie może powstać ani nie może zniknąć. Nie ma także udziału dyfuzji strumienia ponieważ dla płynu w spoczynku, każda zmiana masy oznaczałaby przemieszczenie się cząstek płynu. Dla powierzchni kontrolnej dS, na którą działa prędkość przepływu v⃗ oraz wektor normalny n⃗ , własnością zachowaną jest gęstość ρ.
Dla całej objętości kontrolnej Ω zmiana całkowitej masy w czasie jest równa:
∂
∂t∫ ρdΩ
Ω
(7.9)
Przepływ masy płynu przez określoną powierzchnię w przestrzeni jest równa iloczynowi gęstości, pola tej powierzchni oraz składowej wektora prędkości prostopadłej do powierzchni. Stąd udział konwekcji strumienia w przepływie przez powierzchnię dS wynosi:
ρ(v⃗ ∙ n⃗ )dS (7.10)
Ostatecznie w zapisie całkowym równanie zachowania masy przyjmuje postać:
∂
Równanie zachowania pędu wywodzi się z drugiej zasady dynamiki która stanowi, że zmiana strumienia pędu płynu przepływającego przez objętość kontrolną jest równa wszystkim działającym na tą objętość siłom. Dla najmniejszej części objętości kontrolnej uzyskujemy:
ρv⃗ dΩ (7.12)
Zmiana pędu w objętości kontrolnej w czasie wynosi:
∂
∂t∫ ρv⃗ dΩ
Ω
(7.13)
W tym przypadku własnością zachowaną (stałą) jest iloczyn gęstości i prędkości:
ρv⃗ = [ρu, ρv, ρw]T (7.14) Tensor konwekcji strumienia, który określa przeniesienie pędu przez granicę objętości kontrolnej, w układzie kartezjańskim zawiera trzy składowe:
a) w osi x: ρuv⃗ , b) w osi y: ρvv⃗ , c) w osi z: ρwv⃗ .
Udział tensora konwekcji strumienia w zachowaniu pędu jest równy:
− ∮ ρv⃗ (v⃗ ∙ n⃗ )dS
∂Ω
(7.15)
Dyfuzja strumienia jest równa zero, ponieważ nie jest możliwa dyfuzja momentu dla płynu w spoczynku. Możemy wyróżnić dwie siły działające w objętości kontrolnej:
a) siły zewnętrzne lub siły masowe, działające bezpośrednio na masę objętości. Są to np. siła grawitacji, siła wyporu hydrostatycznego, siły Coriolisa lub siły odśrodkowe. W niektórych przypadkach mogą to też być siły elektromagnetyczne, b) siły powierzchniowe, działające na powierzchnię objętości kontrolnej, ich źródłem
są [9, 64]:
– rozkład ciśnienia wywołany przez płyn otaczający objętość kontrolną,
– ścinanie i naprężenia normalne wynikające z tarcia pomiędzy płynem i powierzchnią objętości kontrolnej
Z powyższego wynika, że siła masowa przypadająca na jednostkę objętości, oznaczona jako ρf e, odpowiada składnikowi źródeł objętości w równaniu (7.8). Udział sił masowych (zewnętrznych) w równaniu zachowania pędu jest równy:
∫ ρf edΩ
Ω
(7.16)
Źródła powierzchniowe sił przedstawione są na rys. 7.2 złożone są z dwóch części – izotropowego ciśnienia i tensora naprężeń lepkościowych τ̿:
QS
̿̿̿̿ = −pI̿ + τ̿ (7.17)
gdzie: I̿ – tensor jednostkowy.
Podsumowując wszystkie powyższe składowe zgodnie z wzorem (7.8) otrzymujemy równanie zachowania momentu dla objętości kontrolnej Ω:
∂
∂t∫ ρvΩ ⃗ dΩ+ ∮ ρv⃗ (v⃗ ∙ n∂Ω ⃗ )dS= ∫ ρf Ω edΩ− ∮ pn∂Ω ⃗ dS+∮ (τ̿ ∙ n∂Ω ⃗ )dS (7.18)
Rys. 7.2. Siły powierzchniowe działające na elementarnej powierzchni objętości kontrolnej [13]
Równanie zachowania energii wynika z pierwszej zasady termodynamiki, która zastosowana do objętości kontrolnej reprezentującej element płynu oznacza, że zmiana energii wewnątrz elementu jest równa sumie strumienia ciepła dopływającego do elementu i zmiany pracy wykonanej na elemencie przez siły masowe i powierzchniowe.
Całkowita energia E przypadająca na jednostkę masy płynu jest równa sumie wewnętrznej energii i energii kinetycznej:
E = e +|v⃗ |2
2 = e +u2+ v2+ w2
2 (7.19)
Własnością zachowaną jest w tym przypadku energia przypadająca na jednostkę objętości - 𝜌𝐸, i jej zmiana w objętości kontrolnej w czasie jest określona:
𝜕
𝜕𝑡∫ 𝜌𝐸𝑑Ω
Ω
(7.20)
Udział konwekcji strumienia w równaniu zachowania energii wyniesie:
− ∮ ρE(v⃗ ∙ n⃗ )dS
∂Ω
(7.21)
W odróżnieniu od równania zachowania momentu, mamy teraz dyfuzję strumienia.
Zgodnie z prawem Fick-a jest ona proporcjonalna do gradientu zachowanej własności na jednostkę masy. Ze względu na to, że dyfuzja strumienia F⃗ D jest określona dla płynu w spoczynku, bierzemy pod uwagę tylko energię wewnętrzną i otrzymujemy:
F⃗ D = −γρκ∇e (7.22)
gdzie: γ = 𝐶𝑝
𝐶𝑣 – stosunek ciepła właściwego przy stałym ciśnieniu i stałej objętości, κ – współczynnik wyrównania temperatur.
Dyfuzja strumienia jest jedną formą przepływu ciepła do objętości kontrolnej, związanej z dyfuzją ciepła przez przewodnictwo termiczne cząsteczek – przepływ ciepła związany z gradientem temperatury. Z tego powodu równanie (7.22) jest zapisywane zgodnie z równaniem Fouriera o przewodnictwie cieplnym:
F⃗ D = −k∇T (7.23)
gdzie: k – współczynnik przewodzenia ciepła, T – temperatura absolutna statyczna.
Kolejną składową przepływu ciepła jest podgrzewanie objętości wskutek absorpcji, promieniowania lub chemiczny reakcji. Źródło ciepła oznaczane jest jako qℎ̇ (w wielkości transferu ciepła na jednostkę masy w czasie). Dodając pracę wykonaną przez siły masowe f e uzyskuje się kompletne źródła objętości Qv:
QV= ρf e∙ v⃗ + qḣ (7.24)
Ostatnią składową są źródła powierzchniowe Q⃗⃗ S. Odpowiadają one pracy wykonanej w czasie przez ciśnienie oraz naprężenia ścinające i normalne:
Q⃗⃗ S = −pv⃗ + τ̿ ∙ v⃗ (7.25) Dodając powyższe składniki otrzymujemy równanie zachowania energii:
∂
∂t∫ ρEdΩΩ + ∮ ρE(v⃗ ∙ n∂Ω ⃗ )dS= ∮ k(∇T ∙ n∂Ω ⃗ )dS+ ∫ (ρf Ω e∙ v⃗ + qḣ )dΩ−
∮ p(v∂Ω ⃗ ∙ n⃗ )dS+ ∮ (τ̿ ∙ v⃗ ) ∙ n∂Ω ⃗ dS
(7.26)
Równanie (7.26) zapisuje się w trochę innej formie. Używając zależności pomiędzy całkowitą entalpią, energią i ciśnieniem:
H = h +|v⃗ |2
2 = E +p
ρ (7.27)
gdzie: H – entalpia całkowita.
Podstawiając człon konwekcji (ρEv⃗ ) i cieśnienia (pv⃗ ) do równania (7.26) i stosując zależność (7.27), można zapisać ostatecznie równanie zachowania energii:
∂
∂t∫ ρEdΩΩ + ∮ ρH(v⃗ ∙ n∂Ω ⃗ )dS= ∮ k(∇T ∙ n∂Ω ⃗ )dS+ ∫ (ρf Ω e∙ v⃗ + qḣ )dΩ+
∮ (τ̿ ∙ v∂Ω ⃗ ) ∙ n⃗ dS
(7.28)
Naprężenia lepkościowe wynikają z tarcia pomiędzy płynem a powierzchnią objętości kontrolnej i opisane są tensorem naprężeń 𝜏̿. W układzie kartezjańskim jego zapis wygląda następująco: działające na sześcienny element płynu. Naprężenia normalne (a) „przesuwają” ścianki sześcianu wzdłuż osi współrzędnych, natomiast naprężenia ścinające (b) „ścinają”
element.
Naprężenia lepkościowe zależą od własności dynamicznych płynu. Dla powietrza czy wody Newton określił, że naprężenia ścinające są proporcjonalne do gradientu prędkości. Taki płyn nazywa się newtonowskim. Inne płyny, jak np. tworzywo sztuczne w stanie płynnym czy krew, zachowują się w innych sposób – są to płyny nienewtonowskie. Dla większości przypadków praktycznych, badany płyn jest newtonowski i stąd składowe tensora naprężeń lepkościowych można obliczyć z wzorów:
gdzie: λ – objętościowy współczynnik lepkości, μ – współczynnik lepkości dynamicznej [Pa·s].
Współczynnik lepkości kinematycznej v wyrażany jest wzorem:
ν = μ
ρ (7.31)
Rys. 7.3. Naprężenia działające na objętość płynu: (a) normalne, (b) ścinające [9, 61]
Równania (7.30) zostały zdefiniowane przez George Stokes’a. Wyrażenia typu 𝜇 (∂u
∂x) w nieprężeniach normalnych reprezentują rozszerzalność liniową – zmiany w kształcie. Z kolei wyrażenia typu λdivv⃗ reprezentują rozszerzalność objętościową – wielkość zmiany objętości, co jest tak naprawdę zmianą gęstości. Stokes wprowadził także zależność:
λ +2
3μ = 0 (7.32)
Powyższe równanie (7.32) jest określane jako lepkość objętościowa. Odpowiada ona rozproszeniu energii w płynie przy stałej temperaturze i podczas określonej zmiany objętości. Za wyjątkiem bardzo wysokich temperatur i ciśnień nie ma żadnego dowodu doświadczalnego, że hipoteza Stockes’a nie stosuję się do równania (7.32) i dlatego jest wykorzystywana do usunięcia z tych zależności współczynnika λ. Dla naprężeń lepkościowych otrzymywane są zatem:
τxx = 2μ (∂u
∂x−1 3divv⃗ ) τyy= 2μ (∂v
∂y−1 3divv⃗ ) τzz = 2μ (∂w
∂z −1 3divv⃗ )
(7.33)
Rozważając przypadek taki, że do równania (7.28) w postaci wektorowej wprowadzi się dwie wielkości wektorowe strumienia – F⃗ C i F⃗ V. Pierwsza z nich, F⃗ C, odnosi się do przepływu konwekcyjnego w płynie i jest nazywana wektorem strumienia konwekcyjnego. Dla równań zachowania pędu i energii zawiera także wyrażenia związane z ciśnieniem – p(n⃗ ) i p(v⃗ ∙ n⃗ )). Drugi wektor strumienia pola (nazwany jako wektor lepkości strumienia F⃗ V) określa naprężenia lepkościowe i dyfuzję ciepła.
Dodatkowo wprowadzając wyrażenie Q⃗⃗ , które zawiera wszystkie źródła objętości związane z siłami masowymi i ogrzewaniem objętościowym. Biorąc powyższe
założenia i wprowadzając iloczyn skalarny z jednostkowym wektorem normalnym n⃗ można złączyć równania razem we wzór:
∂
Wektor wielkości zachowawczych W⃗⃗⃗ w przestrzeni trójwymiarowej zawiera pięć składowych:
Wektor strumienia konwekcyjnego F⃗⃗⃗⃗ zawiera składowe: C
FC
gdzie: V – kontrawariantna prędkość normalna do elementu powierzchniowego dS, definiowana jako iloczyn skalarny wektora prędkości i jednostkowego wektora normalnego:
V = v⃗ ∙ n⃗ = nxu + nyv + nzw (7.37) Entalpia całkowita H w wektorze (7.36) jest obliczana wg wzoru (7.27). Dla wektora lepkości strumienia razem z równaniem (7.29) otrzymuje się:
F⃗ V =
są wyrażeniami opisującymi pracę naprężeń lepkościowych i konwekcji ciepła w płynie. Wektor źródeł jest równy:
Q⃗⃗ = lepkościowych (7.30), zależności (7.31–7.37) są nazywane równaniami Naviera-Stokesa. Opisują one wymianę (strumień) masy, pędu i energii przez granicę dΩ objętości kontrolnej Ω, która jest umiejscowiona na stałe w przestrzeni (rys. 7.1).
Wyprowadzone powyżej równania Naviera-Stokesa przedstawione są w zapisie całkowym, zgodnie z prawami zachowania. Używa się także innych zapisów tych równań, np. w zapisie różniczkowym z zastosowaniem teorii Gauss’a [8, 12, 25, 29].
Równania Naviera-Stokesa w przestrzeni trójwymiarowej składają się z pięciu równań dla pięciu zmiennych zachowawczych ρ, ρu, ρv, ρw i ρE. Zawierają także siedem nieokreślonych zmiennych pola prądu: ρ, u, v, w, E, p i T. Dlatego konieczne jest przedstawienie dwóch dodatkowych równań, w których należy określić zależności termodynamiczne między dwoma parametrami stanu. Przykładowo, określenie ciśnienia jako funkcji gęstości i temperatury oraz wewnętrznej energii lub entalpii będącej funkcją ciśnienia i temperatury. Oprócz tego, należy przedstawić współczynnik lepkości dynamicznej μ i współczynnik przewodzenia ciepła k jako funkcje stanu płynu, w celu zamknięcia całego systemu równań. Powyższe zależności zależą głównie od typu płynu branego pod uwagę.
W przypadku symulacji numerycznych (opisanych w niniejszej pracy zajmujących się przepływem) powietrze traktuje się jako gaz doskonały. Dla takiego płynu równanie stanu ma postać [30, 34, 38]:
p = ρRT (7.41)
gdzie: R – indywidualna stała gazowa [J/kgK].
Entalpię h oblicza się zgodnie z wzorem:
h = cpT (7.42)
Ciśnienie oblicza się jak funkcję zmiennych zachowawczych. Do wyprowadzenia odpowiedniego wzoru należy połączyć równanie (7.27) odnoszące całkowitą entalpię do całkowitej energii z równaniem stanu (7.41). Podstawiając do równania (7.42) na entalpię definicje [1, 4, 26, 31]:
R = cp− cV, γ =cp
cV, (7.43)
otrzymuje się ostatecznie wzór na ciśnienie:
p = (γ − 1)ρ [E −u2 + v2+ w2
2 ] (7.44)
W tym przypadku temperatura jest obliczana zgodnie z wzorem 7.41.
Współczynnik lepkości dynamicznej μ dla gazu doskonałego zależy w dużym stopniu od temperatury a w małym stopniu od ciśnienia. Stosuje się w tym przypadku często wzór Sutherlanda, który dla powietrza wynosi:
μ =1,45T3/2
T + 110 ∙ 10−6 (7.45)
gdzie temperatura T podana jest w stopniach Kelwina (K). Stąd przy T = 288K otrzymuje się μ = 1,78 ∙ 10−5kg/ms. Zależność temperatury od współczynnika przewodzenia ciepła κ dla gazów wynosi:
k = cp μ
Pr (7.46)