Nauka #26 – Dynamika płynów i turbulencja: Navier-Stokes, rozkład Reynoldsa, warstwy przyścienne i kaskada Kołmogorowa

Równania Naviera-Stokesa zapisano w latach 1820–1840 i opisują całą klasyczną mechanikę płynów. Milion dolarów z Nagrody Millenijnej Instytutu Matematycznego Claya czeka na kogoś, kto udowodni (lub obali) tezę, że gładkie rozwiązania zawsze istnieją dla trójwymiarowego przepływu nieściśliwego. Mimo to turbulencja — która wymaga takich rozwiązań — jest już rozumiana na tyle dobrze, by projektować samoloty, okręty podwodne i turbiny gazowe. To produktywne napięcie między matematyczną niekompletnością a inżynierskim sukcesem sprawia, że dynamika płynów jest tak fascynująca.

1. Równania Naviera-Stokesa

Równania Naviera-Stokesa wyrażają drugą zasadę dynamiki Newtona dla elementu płynu: tempo zmiany pędu równa się sumie gradientów ciśnienia, naprężeń lepkich i sił masowych. Dla nieściśliwego płynu newtonowskiego (stała gęstość ρ, lepkość dynamiczna μ) przyjmują one zwartą postać wektorową:

Ciągłość (nieściśliwość):
  ∇ · u = 0    (ρ = const)

Pęd:
  ρ(∂u/∂t + (u·∇)u) = −∇p + μ∇²u + ρg
  ┌┐               ┌┐        ┌┐      ┌┐         ┌┐
  bezwładność    advekcja   ciśnienie   lepkość     siła masowa

Postać bezwymiarowa (skalowanie przez U, L):
  ∂u*/∂t* + (u*·∇*)u* = −∇*p* + (1/Re)∇*²u*

Liczba Reynoldsa Re = ρUL/μ = UL/ν    (ν = μ/ρ = lepkość kinematyczna)
  Re << 1: dominują siły lepkie → przepływ laminarny, pełzający (przepływ Stokesa)
  Re ~ 100: laminarny ze stałą recyrkulacją (ślad za walcem)
  Re ~ 10^3: okresowe odrywanie wirów (ulica Kármána, St = fD/U ≈ 0.2)
  Re > 10^4: przejściowy;
  Re > 10^5: w pełni turbulentny

Równanie Bernoulliego (nielepkie, wzdłuż linii prądu):
  p + ½ρu² + ρgz = const   (energia na jednostkę objętości)
  Pochodzenie: ∇(½u²) = (u·∇)u + u×(∇×u)

Nieliniowy człon advekcji (u·∇)u jest źródłem niemal całej złożoności mechaniki płynów. Przekazuje on energię kinetyczną między skalami w sposób, którego liniowy człon lepki nie potrafi cofnąć, dając początek wirom, śladom aerodynamicznym i ostatecznie turbulencji. Przy niskich liczbach Reynoldsa człon ten jest pomijalny w porównaniu z lepkością, a przepływy są gładkie i przewidywalne. Przy wysokich liczbach Reynoldsa lepkość może rozpraszać jedynie najdrobniejsze skale; ruchy wielkoskalowe są zasadniczo nielepkie, a ich energia musi przejść przez coraz mniejsze wiry, zanim zostanie zamieniona w ciepło.

2. Teoria warstwy przyściennej

Gdy płyn opływa powierzchnię ciała stałego, warunek braku poślizgu wymaga, aby prędkość płynu przy powierzchni zrównała się z prędkością ściany. Szybkie przejście od zera (przy ścianie) do prędkości swobodnego strumienia U zachodzi w cienkiej warstwie przyściennej, po raz pierwszy przeanalizowanej przez Ludwiga Prandtla w 1904 roku. Ta obserwacja otworzyła drogę do praktycznej aerodynamiki: poza cienką warstwą przyścienną przepływ jest w przybliżeniu nielepki; wewnątrz niej dominują efekty lepkości.

Laminarna warstwa przyścienna Blasiusa (płaska płyta, zerowy gradient ciśnienia):
  δ(x) / x = 5.0 / √Re_x       (grubość przy 99%)
  δ*(x) / x = 1.72 / √Re_x      (grubość wyparcia: brakujący strumień masy)
  C_f(x) = 0.664 / √Re_x            (lokalny współczynnik tarcia powierzchniowego)

Przejście do turbulencji (płaska płyta):
  Re_x,krit ≈ 5×10^5 (poziom turbulencji swobodnego strumienia ~0,1%)
  Wyższa turbulencja swobodnego strumienia lub niekorzystny gradient ciśnienia → wcześniejsze przejście

Turbulentna warstwa przyścienna (profil potęgowy):
  u/U ≈ (y/δ)^{1/7}    (prawo 1/7, Re > 10^5)
  C_f = 0.0592 Re_x^{-1/5}

Prawo ściany (warstwa wewnętrzna):
  u^+ = y^+                    (podwarstwa lepka, y^+ < 5)
  u^+ = (1/κ) ln(y^+) + B   (obszar prawa logarytmicznego, 30 < y^+ < 200)
  y^+ = y u_τ / ν,   u_τ = √(τ_w/ρ) (prędkość tarcia)
  κ ≈ 0.41 (stała von Kármána),  B ≈ 5.0

Oderwanie warstwy przyściennej:
  Występuje, gdy gradient ciśnienia dP/dx > 0 (niekorzystny gradient ciśnienia)
  Płyn przy ścianie zwalnia do zera; odwrócenie → bąbel odrywania → ślad
  Kontrola oderwania: ssanie, generatory wirów, turbulatory (dołki na piłce golfowej)

3. Niestabilności hydrodynamiczne

Przepływ laminarny jest niestabilny, jeśli małe zaburzenia rosną w czasie zamiast zanikać. Liniowa analiza stabilności określa, czy dany przepływ bazowy jest stabilny, wprowadzając małe zaburzenie i sprawdzając, czy zlinearyzowane równania dopuszczają rosnące rozwiązania. Dwie niestabilności są szczególnie ważne dla przejścia do turbulencji i dynamiki warstwy mieszania: Kelvina-Helmholtza i Rayleigha-Taylora.

Niestabilność Kelvina-Helmholtza (ścinanie na granicy ośrodków):
  Dwa nielepkie płyny o prędkościach U_1, U_2 i gęstościach ρ_1, ρ_2 stykające się przy y=0
  Tempo wzrostu dla liczby falowej k:
    σ^2 = −gk(ρ_1−ρ_2)/(ρ_1+ρ_2) + k²ρ_1ρ_2(U_1−U_2)²/(ρ_1+ρ_2)²
  Czyste ścinanie (ta sama gęstość): zawsze niestabilne dla dowolnego k > 0
  Stabilizowane przez grawitację i stratyfikację gęstości (liczba Richardsona Ri = N²/S²)
  Ri < 1/4 (kryterium Milesa-Howarda): ścinanie przezwycięża stratyfikację → wiry KH

Niestabilność Rayleigha-Taylora (ciężki płyn nad lekkim):
  Ciężki płyn (ρ_2) spoczywający na lżejszym płynie (ρ_1) w polu grawitacyjnym g
  Tempo wzrostu: σ = √[Atkg]    (A = liczba Atwooda = (ρ_2−ρ_1)/(ρ_2+ρ_1))
  Wszystkie długości fal są niestabilne, ale krótkie fale rosną szybciej; napięcie powierzchniowe stabilizuje k > k_c
  Nieliniowo: grzybkowate pióropusze ciężkiego płynu penetrujące w dół (kolce) i lekkiego płynu unoszącego się (bąble)
  Zastosowania: wyrzuty materii supernowych, implozja kapsuły w fuzji inercyjnej, mieszanie oceaniczne

Konwekcja Rayleigha-Bénarda (napędzana wyporem):
  Warstwa płynu ogrzewana od dołu, chłodzona od góry
  Liczba Rayleigha: Ra = gβΔT L³ / (να)  (β = rozszerzalność cieplna, α = dyfuzyjność)
  Początek (przewodnictwo → wałki konwekcyjne): Ra_c = 1708  (Chandrasekhar 1961)
  Ra > 10^4: oscylacyjny; Ra > 10^6: turbulentne pióropusze
  Liczba Nusselta Nu = h L / λ ~ Ra^{1/3} (reżim turbulentny)

4. Rozkład Reynoldsa i równania RANS

W pełni turbulentne przepływy zawierają wiry w ogromnym zakresie skal. Bezpośrednie symulowanie każdej skali długości i czasu (symulacja bezpośrednia, DNS) wymaga liczby punktów siatki proporcjonalnej do Re^{9/4} — dla turbiny silnika odrzutowego przy Re ~ 10^7 to ~10^{15} punktów, znacznie przekraczające obecną moc obliczeniową. Modele inżynierskie zamiast tego rozkładają prędkość na część średnią i fluktuującą: rozkład Reynoldsa.

Rozkład Reynoldsa:
  u_i(x,t) = U_i(x) + u_i'(x,t)
  <u_i> = U_i  (średnia czasowa),  <u_i'> = 0

Uśrednione Reynoldsem równanie pędu:
  ρ U_j ∂U_i/∂x_j = −∂P/∂x_i + ∂/∂x_j [μ∂U_i/∂x_j − ρ<u_i'u_j'>]
                                                └──────────────────────└
                                                tensor naprężeń Reynoldsa τ^R_ij

Problem domknięcia:
  τ^R_ij ma 6 niezależnych niewiadomych na punkt → więcej niewiadomych niż równań
  Hipoteza Boussinesqa: τ^R_ij = −2ν_t S_ij  (ν_t = turbulentna lepkość wirowa)
  Analogia z lepkością molekularną; zawodzi w silnie anizotropowych przepływach

Model k-ε (Launder i Spalding, 1974 — najczęściej stosowany model RANS):
  Równanie transportu turbulentnej energii kinetycznej: k = ½<u_i'u_i'>
  Równanie transportu tempa dyssypacji: ε = ν<∂u_i'/∂x_j ∂u_i'/∂x_j>
  Lepkość wirowa: ν_t = C_μ k² / ε  (C_μ = 0.09)
  Skala długości turbulencji swobodnego strumienia: L_t = k^{3/2}/ε

Model k-ω SST (Menter 1994 — standard w aerodynamice):
  Łączy k-ω blisko ściany (dokładny w podwarstwie) z k-ε w swobodnym strumieniu
  Lepsze przewidywanie oderwania przy niekorzystnym gradiencie ciśnienia, przeciągnięcia profilu
  Stosowany w ANSYS Fluent, OpenFOAM, SU2

5. Teoria kaskady energii Kołmogorowa

W 1941 roku Andrieij Kołmogorow sformułował statystyczną teorię w pełni rozwiniętej turbulencji, która mimo swojej prostoty pozostaje fundamentem tej dziedziny. Centralną ideą jest kaskada energii Richardsona: energia jest wprowadzana w dużych skalach (skala całkowa L), przekazywana do coraz mniejszych wirów bez istotnej dyssypacji i ostatecznie rozpraszana jako ciepło w lepkiej skali Kołmogorowa η. W pośrednim zakresie inercyjnym statystyki fluktuacji prędkości są uniwersalne.

Mikroskale Kołmogorowa (zależą tylko od ν i ε):
  Długość: η = (ν³/ε)^{1/4}
  Czas:   τ_η = (ν/ε)^{1/2}
  Prędkość: u_η = (νε)^{1/4}

Rozdzielenie skal:
  L / η ~ Re^{3/4}       (od skali całkowej do skali Kołmogorowa)
  Dla atmosferycznej warstwy przyściennej (Re ~ 10^8): L/η ~ 10^6

Pierwsza hipoteza podobieństwa (lokalna izotropia dla małych skal):
  W skalach r << L statystyki przyrostów prędkości są uniwersalne
  (zależą wyłącznie od ν i ε)

Druga hipoteza podobieństwa (zakres inercyjny, η << r << L):
  Statystyki zależą wyłącznie od ε (lepkość nieistotna)
  Funkcja strukturalna: S_2(r) = <|u(x+r) − u(x)|²> = C_2 (εr)^{2/3}
  C_2 ≈ 2.0 (stała Kołmogorowa)

Widmo energii (prawo −5/3 Kołmogorowa):
  E(k) = C_K ε^{2/3} k^{−5/3}   (zakres inercyjny)
  C_K ≈ 1.5 (spektralna stała Kołmogorowa)
  Pierwsza weryfikacja eksperymentalna: Grant i in. (1962) w kanale pływowym

Poprawki intermitencji (K62, Obuchow 1962):
  Rzeczywista turbulencja jest przestrzennie przerywana (intermitentna); dyssypacja jest fraktalna
  Wykładnik skalowania: ζ_p = p/3 − μp(p−3)/18   (μ ≈ 0.25)
  Formalizm multifraktalny (Parisi i Frisch 1985)

Możliwości DNS (2024):
  Maksymalne DNS Re_λ ~ 2000 (liczba Reynoldsa w skali Taylora)
  Siatka: 12 288³ ≈ 2×10^12 punktów  (działa na 10 najlepszych superkomputerach)
  Cel: Re_λ ≈ 10^4 (porównywalne z warstwą przyścienną atmosfery) w ciągu ~15 lat

6. Symulacja dużych wirów i zastosowania

Pomiędzy siłową dokładnością DNS a szybkością RANS mieści się symulacja dużych wirów (LES): bezpośrednio oblicza się wiry większe niż siatka, a modeluje się jedynie efekty podsiatkowe. LES rejestruje niestacjonarne mieszanie, struktury spójne i aeroakustykę, które RANS całkowicie pomija, przy ułamku kosztu DNS. Nowoczesne przemysłowe CFD coraz częściej wykorzystuje hybrydowe metody RANS-LES, takie jak Detached Eddy Simulation (DES) i Scale-Adaptive Simulation (SAS).

Filtr siatkowy (kreska oznacza wielkość przefiltrowaną):
  ū(x) = ∫ G(x − x') u(x') dx'    (G = filtr prostokątny lub Gaussa, szerokość Δ)
  ū zawiera skale > Δ;  naprężenie podsiatkowe τ_ij^sgs = ​u_i u_j – ū_i ū_j

Model Smagorinsky'ego (1963):
  τ_ij^sgs = −2(C_s Δ)² |S̄| S̄_ij
  C_s ≈ 0.1–0.18 (uniwersalne dla turbulencji izotropowej; wymaga redukcji blisko ścian)

Dynamiczny model Smagorinsky'ego (Germano 1991):
  C_s obliczane lokalnie ze stosunku naprężeń filtru testowego do filtru siatkowego
  Eliminuje potrzebę ręcznego strojenia; automatycznie obsługuje przejście

LES modelowana przy ścianie (WMLES):
  Modelowanie obszaru blisko ściany za pomocą RANS lub algebraicznego prawa ściany → grubsza siatka blisko powierzchni
  N_siatka ~ Re^{13/7} (WMLES) vs Re^{37/14} (LES rozwiązująca ścianę) vs Re^{9/4} (DNS)
  WMLES całego samolotu przy Re ~ 10^8 wykonalna z ~10^9 punktami

Zastosowania:
  Spalanie:          LES rozwiązuje turbulentne szczotkowanie płomienia, palniki wirowe przedmieszania
  Energetyka wiatrowa: LES atmosferycznej warstwy przyściennej + śladów za turbinami (SOWFA, AMR-Wind)
  Aeroakustyka:       fluktuacje ciśnienia z LES + równanie FW-H → hałas w polu dalekim
  Klimat:             mezoskalowe wiry oceaniczne (~30 km) rozwiązywane obecnie w modelach CMIP7 o wysokiej rozdzielczości

Powiązane symulacje