Devlog #84 – Fala 63: Przepływ Stokesa, upakowanie kul, krzywa Hilberta, spirala Ulama, przepływ maksymalny i halo lodowe

Nowe symulacje

🌊 Przepływ Stokesa — fizyka ruchu pełzającego

Równania Stokesa (Re → 0)

Przy bardzo niskiej liczbie Reynoldsa Re = ρUL/μ ≪ 1 siły bezwładności są pomijalne w porównaniu z siłami lepkości, a równania Naviera-Stokesa redukują się do liniowych równań Stokesa:

∇p = μ ∇²u       (pęd: gradient ciśnienia równoważy dyfuzję lepką)
∇ · u = 0        (ciągłość: przepływ nieściśliwy)

W postaci funkcji prądu 2D (u = ∂ψ/∂y, v = −∂ψ/∂x):
  ∇⁴ψ = 0         (równanie biharmoniczne)

Ponieważ równania są liniowe i nie zawierają pochodnych czasowych (czas pojawia się jedynie poprzez warunki brzegowe), przepływ Stokesa jest quasi-statyczny: pole prędkości dostosowuje się natychmiastowo do warunków brzegowych. Ma to głęboką konsekwencję — odwracalność.

Funkcja prądu dla przepływu wokół kuli

Dla jednorodnego przepływu U wokół kuli o promieniu R, dokładna funkcja prądu Stokesa we współrzędnych sferycznych (r, θ) wynosi:

ψ(r, θ) = U · sin²θ · [r²/2  -  3R·r/4  +  R³/(4r)]

Składowe prędkości:
  u_r = (1/r² sinθ) · ∂ψ/∂θ = U cosθ · [1 - 3R/(2r) + R³/(2r³)]
  u_θ = -(1/r sinθ) · ∂ψ/∂r = -U sinθ · [1 - 3R/(4r) - R³/(4r³)]

Siła oporu (opór Stokesa):
  F_drag = 6π μ R U

Linie prądu są symetryczne z przodu i z tyłu kuli — właściwość unikalna dla przepływu Stokesa. Dla porównania, przy wysokim Re ślad przepływu łamie symetrię, a opór drastycznie rośnie.

Demonstracja odwracalności

Najbardziej uderzającą cechą symulacji jest demonstracja odwracalności. Barwne kropelki barwnika są wstrzykiwane i unoszone przez pole prędkości Stokesa. Gdy przepływ jest odwrócony (U → −U), cząstki barwnika dokładnie odtwarzają swoje ścieżki i ponownie zbiegają się do pierwotnych pozycji — system nie ma efektywnego „mieszania”, ponieważ liniowość przepływu Stokesa oznacza brak chaotycznego rozciągania i składania. To słynny eksperyment Taylora z syropem kukurydzianym z 1966 roku. Ten sam pokaz przy przepływie Naviera-Stokesa przy Re = 100 pokazałby kropelki nieodwracalnie rozmazane przez mieszanie turbulentne.

Twierdzenie grzebienia

Pojedyncza pozioma linia barwnika jest dzielona na N odcinków, z których każdy jest zmuszony przejść przez osobny równoległy cylinder. Po przejściu przez N cylindrów, N osobnych nitek barwnika łączy się z powrotem w jedną poziomą linię — identyczną z oryginałem. Jest to konsekwencja liniowości i odwracalności przepływu Stokesa: ścieżka każdej cząstki płynu przez przeszkodę w przepływie Stokesa jest jednoznacznie determinowana pozycją wejścia i odwracalna, więc pozycje wyjścia gwarantowanie odtwarzają się w porządku.

🔵 Upakowanie kul — hipoteza Keplera i współczynniki upakowania

Współczynniki upakowania 2D

W 2D znane są dwa regularne (sieciowe) upakowania. Upakowanie kwadratowe umieszcza środki okręgów na siatce kwadratu jednostkowego; upakowanie heksagonalne (trójkątne) przesuwa naprzemienne rzędy o połowę średnicy:

Upakowanie kwadratowe:
  środki w (i, j) dla i,j ∈ ℤ; promień r = 0,5
  pole okręgu / pole komórki jednostkowej = π r² / 1² = π/4 ≈ 0,7854

Upakowanie heksagonalne:
  środki w (i + j/2, j·√3/2); promień r = 0,5
  pole komórki jednostkowej = √3/2 (romb)
  współczynnik upakowania = π r² / (√3/2) = π/(2√3) ≈ 0,9069

Losowe sekwencyjne dodawanie (RSA):
  Umieszczaj dyski po jednym w losowych pozycjach, jeśli nie ma nakładania.
  Współczynnik upakowania zakleszczonego ≈ 0,547 (stała parkingowa Renyiego)

Upakowanie heksagonalne jest dowodnie optymalne w 2D — żaden układ jednakowych okręgów nie może osiągnąć wyższej gęstości. Udowodnił to Thue w 1910 roku. Symulacja 2D pokazuje bieżący współczynnik upakowania i charakterystyczny heksagonalny wzór pustek między okręgami.

Sieci krystaliczne 3D i hipoteza Keplera

Problem 3D jest znacznie trudniejszy. Cztery klasyczne sieciowe upakowania zaimplementowane w symulacji to:

Sieć prosta regularna (SC):
  a = 2r;  współczynnik upakowania = π/6 ≈ 0,5236;  liczba kontaktów Z = 6

Sieć centrowana przestrzennie (BCC):
  a = 4r/√3;  współczynnik upakowania = π√3/8 ≈ 0,6802;  Z = 8

Sieć ściennie centrowana (FCC):
  a = r√8;  współczynnik upakowania = π/(3√2) ≈ 0,7405;  Z = 12

Upakowanie heksagonalne zwarte (HCP):
  Ten sam współczynnik upakowania co FCC = π/(3√2) ≈ 0,7405;  Z = 12
  Różni się od FCC jedynie sekwencją ułożenia warstw (ABAB vs ABCABC)

Kepler przypuszczał w 1611 roku, że π/(3√2) ≈ 74,05% to maksymalny możliwy współczynnik upakowania dla jednakowych kul. Zostało to formalnie udowodnione przez Thomasa Halesa w 1998 roku za pomocą dowodu wspomaganego komputerowo liczącego 250 stron plus 3 GB kodu — jednego z najdłuższych dowodów w matematyce. Zarówno FCC, jak i HCP osiągają to maksimum.

Projekcja 3D i wizualizacja kontaktów

Prawy panel pokazuje animowaną projekcję 3D wybranej sieci. Każda kula renderowana jest półprzezroczyście, aby widoczna była wewnętrzna struktura. Wiązania kontaktowe są rysowane między stykającymi się parami kul (odległość = 2r w granicach tolerancji zmiennoprzecinkowej). Liczba kontaktów Z — liczba najbliższych sąsiadów stykających się z daną kulą — wizualnie potwierdza typ sieci i jest równa liczbie pocałunkowej w danej geometrii.

🌀 Krzywa Hilberta — krzywe wypełniające przestrzeń i lokalność

Konstrukcja przez rekurencyjne podstawienie

Krzywa Hilberta jest generowana przez gramatykę podstawieniową zastosowaną do wzorca początkowego. W rzędzie 1 jest to kształt ∪ łączący 4 punkty. W rzędzie n każdy z 4 pod-kwadratów z rzędu n−1 jest zastępowany obróconą/odbitą kopią krzywej rzędu n−1, połączoną segmentami mostkującymi. Krzywa rzędu n odwiedza 4ⁿ punktów i ma całkowitą długość (4ⁿ − 1) × (rozmiar komórki jednostkowej), która rozbiega się, gdy n → ∞, podczas gdy krzywa pozostaje ciągła i zawarta w kwadracie jednostkowym.

// mapowanie d-na-(x,y) (manipulacja bitami, wg Skillinga 2004):
function d2xy(n, d):
    rx, ry, s = 0, 0, 1
    x, y = 0, 0
    while s < n:
        rx = (d >> 1) & 1
        ry = (d ^ rx) & 1
        rotate(s, x, y, rx, ry)
        x += s * rx
        y += s * ry
        d >>= 2
        s <<= 1
    return x, y

// Odwrotność (x,y)-na-d:
function xy2d(n, x, y):
    rx, ry, s = 0, 0, n >> 1
    d = 0
    while s > 0:
        rx = (x & s) > 0 ? 1 : 0
        ry = (y & s) > 0 ? 1 : 0
        d += s * s * ((3 * rx) ^ ry)
        rotate(s, x, y, rx, ry)
        s >>= 1
    return d

Zachowanie lokalności

Kluczową właściwością, która czyni krzywą Hilberta praktycznie użyteczną (w zarządzaniu pamięcią podręczną CPU, drzewach R w bazach danych, indeksowaniu przestrzennym GIS), jest to, że punkty blisko siebie na krzywej 1D są również blisko siebie w przestrzeni 2D, i odwrotnie. To zachowanie lokalności można skwantyfikować: dla dwóch punktów 2D o odległości euklidesowej d, ich indeksy Hilberta różnią się co najwyżej o O(d). Dla porównania, uporządkowanie wierszowe („skanowanie rastrowe”) ma w najgorszym przypadku pary na przeciwnych krańcach przestrzeni indeksów, które w 2D dzieli tylko jeden rząd, a Z-order (Morton) jest pośredni.

Krzywe Peana i Mortona

Krzywa Peana (1890, pierwsza odkryta krzywa wypełniająca przestrzeń) dzieli każdy kwadrat na 9 pod-kwadratów w wzorze wężowym, odwiedzając 9ⁿ punktów w rzędzie n. Krzywa Z-order Mortona przeplata bity współrzędnych x i y:

// kodowanie Mortona (Z-order):
morton(x, y) = spread_bits(x) | (spread_bits(y) << 1)

spread_bits(v):  // wstaw bit 0 między każdy bit v
    v = (v | (v << 8)) & 0x00FF00FF
    v = (v | (v << 4)) & 0x0F0F0F0F
    v = (v | (v << 2)) & 0x33333333
    v = (v | (v << 1)) & 0x55555555
    return v

Dymek po najechaniu pokazuje zarówno mapowanie w przód (odległość na krzywej d → pozycja siatki (x,y)), jak i odwrotność, pozwalając na jeden rzut oka zobaczyć, jak lokalność jest zachowana lub naruszona dla dowolnie wybranego punktu.

🔢 Spirala Ulama — wzory liczb pierwszych i grzbiety wielomianowe

Konstrukcja spirali Ulama

W 1963 roku Stanislaw Ulam, znudzony podczas wykładu konferencyjnego, zaczął zapisywać liczby całkowite w kwadratowej spirali i zakreślać liczby pierwsze. Zauważył, że liczby pierwsze mają tendencję do skupiania się na przekątnych — znacznie częściej niż wynikałoby to z przypadku. Spirala przypisuje liczbę całkowitą n do pozycji siatki, wijąc się na zewnątrz od centrum:

n=1 → (0,0)
n=2 → (1,0)  n=3 → (1,1)  n=4 → (0,1)  n=5 → (-1,1)
n=6 → (-1,0) n=7 → (-1,-1) n=8 → (0,-1) n=9 → (1,-1) ...

Ogólne wicie (ruch w prawo 1, w górę 1, w lewo 2, w dół 2, w prawo 3, ...):
  długości segmentów: 1,1,2,2,3,3,4,4,...
  kierunki: E,N,W,S,E,N,W,S,...

Symulacja oblicza pozycje spiralne dla wszystkich N do ~90 000, korzystając z wcześniej obliczonej tabeli wicia — wystarczająco szybko, aby renderować pełną siatkę przy 60 fps z interaktywnym powiększaniem i przesuwaniem.

Dlaczego przekątne? Wielomiany kwadratowe

Przekątne odpowiadają wielomianom kwadratowym w n. Wzdłuż jednej głównej przekątnej: wartości to f(k) = 4k² + 2k + 1 dla k = 0, 1, 2, … Wzdłuż drugiej: f(k) = 4k² − 2k + 1. Każdy wielomian kwadratowy dający wiele liczb pierwszych rozświetli się jako jasna przekątna na spirali Ulama. Najsłynniejszym przykładem jest wielomian Eulera f(n) = n² + n + 41, który jest pierwszy dla n = 0, 1, 2, …, 39 (40 kolejnych wartości) — rekord, który przetrwał wieki. Na spirali pojawia się to jako niewątpliwie jasny grzbiet diagonalny, przełączalny w symulacji.

Spirala Sacksa

Spirala Sacksa (Robert Sacks, 1994) mapuje liczbę całkowitą n na współrzędne biegunowe, umieszczając ją pod kątem θ = 2π√n i promieniem r = √n — tak że kwadraty doskonałe leżą na dodatniej osi x, a spirala wykonuje dokładnie jeden obrót na każdy kwadrat doskonały:

θ(n) = 2π √n
r(n) = √n
x(n) = r cosθ,  y(n) = r sinθ

Liczby pierwsze na spirali Sacksa tworzą zakrzywione łuki zamiast prostych przekątnych, ale skupienie jest jeszcze bardziej uderzające wizualnie. Każdy łuk odpowiada innej klasie reszt kwadratowych modulo małe liczby pierwsze. Symulacja renderuje obie spirale obok siebie z zsynchronizowanym powiększeniem.

📦 Przepływ maksymalny — Edmonds-Karp i twierdzenie o przekroju minimalnym

Sieci przepływowe i problem przepływu maksymalnego

Sieć przepływowa to graf skierowany G = (V, E), gdzie każda krawędź (u, v) ma nieujemną pojemność c(u, v). Przepływ f(u, v) na każdej krawędzi musi spełniać: (1) ograniczenie pojemności f(u, v) ≤ c(u, v), oraz (2) zachowanie przepływu: przepływ wchodzący do dowolnego węzła wewnętrznego równa się przepływowi wychodzącemu. Problem przepływu maksymalnego pyta o maksymalny całkowity przepływ ze źródła s do ujścia t.

Twierdzenie przepływ maksymalny = przekrój minimalny (Ford-Fulkerson, 1956):
  Maksymalny przepływ s-t równa się minimalnej pojemności przekroju s-t.

Przekrój s-t: podział (S, T) zbioru V z s ∈ S, t ∈ T
Pojemność przekroju = Σ c(u,v) dla krawędzi (u,v) z u ∈ S, v ∈ T

Algorytm Edmondsa-Karpa

Algorytm Edmondsa-Karpa (1972) to Ford-Fulkerson z konkretnym wyborem ścieżki powiększającej: zawsze używaj najkrótszej ścieżki (najmniej krawędzi) z s do t w grafie rezydualnym, znajdowanej metodą BFS. Ten wybór gwarantuje wielomianową złożoność czasową O(VE²), podczas gdy ogólny Ford-Fulkerson może wymagać O(E × max_flow) kroków przy niefortunnym wyborze ścieżek.

Edmonds-Karp:
1. Zbuduj graf rezydualny G_f: dla każdej krawędzi (u,v) o pojemności c i przepływie f,
       dodaj krawędź w przód o pojemności rezydualnej c - f
       dodaj krawędź wsteczną o pojemności rezydualnej f
2. Znajdź najkrótszą ścieżkę powiększającą s→t przez BFS w G_f
3. Wąskie gardło Δ = minimalna pojemność rezydualna wzdłuż ścieżki
4. Powiększ: zwiększ przepływ na każdej krawędzi w przód o Δ, zmniejsz wsteczną o Δ
5. Powtarzaj, aż nie istnieje ścieżka powiększająca (BFS nie znajduje ścieżki s→t)

Złożoność: O(VE²) — co najwyżej VE/2 powiększeń, każde BFS w O(E)

Identyfikacja przekroju minimalnego i animacja

Po zakończeniu algorytmu przekrój minimalny jest identyfikowany przez BFS w grafie rezydualnym startujący z s: wszystkie wierzchołki osiągalne z s tworzą zbiór S; wszystkie pozostałe tworzą T. Każda krawędź przechodząca z S do T przy pełnej pojemności jest krawędzią przekroju — rysowaną na czerwono. Symulacja pokazuje każdy krok powiększający BFS jeden po drugim: wybrana ścieżka podświetlona jest na żółto, wartości przepływu na każdej krawędzi aktualizują się liczbowo, a pojemności rezydualne animują się. Po zbieżności animacja przekroju minimalnego rysuje granicę podziału i etykietuje ją pojemnością (równą przepływowi maksymalnemu).

• Preset 1 — klasyczny przykład z podręcznika; przepływ maksymalny = przekrój minimalny = 23.
• Preset 2 — sieć diamentowa; demonstruje przekrój o rozmiarze 3 oddzielający źródło od ujścia.
• Preset 3 — rurociąg z krawędzią wąskiego gardła; pojedyncza krawędź przekroju jest trywialnie przekrojem minimalnym.
• Preset 4 — sformułowanie dopasowania dwudzielnego; przepływ maksymalny daje maksymalny rozmiar dopasowania.
• Preset 5 — duża losowa sieć; demonstruje wydajność przy 20 węzłach i 40 krawędziach.

🌈 Halo lodowe — optyka Monte Carlo atmosferycznych kryształów lodu

Fizyka halo lodowych

Atmosferyczne halo lodowe powstają wskutek załamania i odbicia światła słonecznego na ściankach heksagonalnych pryzmatów lodowych (struktura krystaliczna Ih). Pryzmat heksagonalny ma dwie heksagonalne ścianki podstawowe i sześć prostokątnych ścianek bocznych. Kluczowy kąt załamania zależy od kąta pryzmatu α między dwiema ściankami i współczynnika załamania n lodu dla światła widzialnego (n ≈ 1,31 przy 550 nm, dyspersyjnie od ~1,307 dla czerwieni do ~1,317 dla fioletu).

Prawo Snella i minimalna dewiacja

Dla promienia wchodzącego przez ściankę 1 i wychodzącego przez ściankę 2 pryzmatu o kącie klina α, kąt dewiacji δ wynosi:

Przy wejściu:   n_air · sin(θ_i) = n_ice · sin(θ_r)
Przy wyjściu:    n_ice · sin(θ_r') = n_air · sin(θ_e)

Geometria:   θ_r + θ_r' = α   (kąt między ściankami pryzmatu)
Dewiacja całkowita: δ = θ_i + θ_e - α

Dewiacja minimalna (δ_min) występuje, gdy θ_i = θ_e (promień symetryczny):
  sin((δ_min + α)/2) = n · sin(α/2)

Dla halo 22°: α = 60° (naprzemienne ścianki pryzmatu), n = 1,31
  sin((δ_min + 60°)/2) = 1,31 · sin(30°) = 0,655
  (δ_min + 60°)/2 = arcsin(0,655) ≈ 40,9°
  δ_min ≈ 21,8° ≈ 22°

Kąt dewiacji minimalnej definiuje wewnętrzną krawędź pierścienia halo: żaden promień nie może wyjść z dewiacją mniejszą niż δ_min dla danej geometrii pryzmatu. Promienie kumulują się w pobliżu δ_min (faza stacjonarna), tworząc jasny pierścień przy tym promieniu kątowym.

Raytracing metodą Monte Carlo

Zamiast całkować analitycznie rozkład dewiacji, symulacja wystrzeliwuje N = 100 000 promieni na klatkę, każdy przez losowo zorientowany (lub ograniczony) kryształ lodu:

dla każdego promienia:
    1. Próbkuj orientację kryształu:
         losowa → wszystkie promienie halo 22°
         płytka (przechylenie σ ≈ 1°) → słońce poboczne / parhelium
         kolumna (przechylenie σ ≈ 1°) → łuk okołozenitalny
    2. Oblicz trafienie ścianki wejściowej: przecięcie promienia z geometrią pryzmatu
    3. Zastosuj prawo Snella przy wejściu → kierunek załamany wewnątrz kryształu
    4. Propaguj do ścianki wyjściowej; zastosuj prawo Snella ponownie
    5. Oblicz (azymut, elewację) promienia wyjściowego
    6. Akumuluj w histogramie kątowym 2D (projekcja nieba)
    7. Zastosuj dyspersję długości fali: powtórz dla λ = 400, 550, 700 nm

Wynikowy histogram jest normalizowany i renderowany jako obraz nieba w kątowym układzie współrzędnych. Halo 22° pojawia się jako rozproszony jasny pierścień o promieniu kątowym 22° od symulowanego słońca, jaśniejszy na wewnętrznej krawędzi (dewiacja minimalna). Czerwień znajduje się wewnątrz pierścienia (mniejsze n → większe δ_min), a błękit na zewnątrz — odwrotnie niż w tęczy.

Słońca poboczne i łuk okołozenitalny

Słońca poboczne (parhelia) powstają wskutek kryształów płytkowych, które osadzają się z osią c pionową na skutek oporu aerodynamicznego. Promienie wchodzą przez pionową ściankę pryzmatu i wychodzą przez inną pionową ściankę pod kątem 60°, ale teraz tylko kryształy o właściwym kącie azymutalnym przyczyniają się, koncentrując światło przy ±22° od słońca wzdłuż horyzontu:

Warunek parhelium: światło wchodzi/wychodzi przez pionowe ścianki pryzmatu
  Przechylenie kryształu od poziomu: σ ≈ 0°–2° (osadzanie płytek)
  Pozycja kątowa: przy wysokości słońca h, słońce poboczne jest przesunięte o
    Δaz = arctan(tan(22°)/cos(h))
  Przy h = 0°: Δaz = 22°;  przy h = 30°: Δaz ≈ 26°

Łuk okołozenitalny (CZA) powstaje, gdy światło wchodzi przez górną ściankę podstawową i wychodzi przez pionową ściankę pryzmatu kryształów płytkowych. Geometria daje zakrzywiony łuk 46° powyżej słońca, z czerwienią na zewnątrz (na dole) i fioletem wewnątrz — tworząc najbardziej intensywnie zabarwiony łuk w optyce atmosferycznej, często bardziej nasycony niż tęcza.

Co dalej

Fala 64 jest w trakcie tworzenia, z kandydatami obejmującymi aerodynamikę współczynnika oporu, geometrię kopuły geodezyjnej, stabilne dopasowanie (Gale-Shapley), wzory interferencyjne hologramów, algorytm wykrywania społeczności Louvain oraz strukturę ziaren / krystalografię Woronoja. Biblioteka jest na dobrej drodze, aby przekroczyć 600 symulacji do lata.