Devlog #85 – Fala 64: Współczynnik oporu, kopuła geodezyjna, stabilne dopasowanie, hologram, Louvain i struktura ziaren

Nowe symulacje

🪂 Współczynnik oporu — aerodynamika w funkcji kształtu i liczby Reynoldsa

Równanie siły oporu

Siła oporu działająca na dowolne ciało poruszające się w płynie opisana jest pojedynczym bezwymiarowym współczynnikiem C_d, który pochłania całą geometryczną złożoność odrywania przepływu, formowania śladu i odzysku ciśnienia:

F_d = ½ · ρ · v² · C_d · A

gdzie:
  ρ  = gęstość płynu (kg/m³)
  v  = prędkość względna (m/s)
  A  = pole odniesienia (pole czołowe dla brył nieopływowych, pole rzutu dla profili lotniczych)
  C_d = współczynnik oporu (bezwymiarowy)

Dekompozycja:
  C_d = C_d,ciśnienia + C_d,tarcia
  Opór ciśnieniowy (formy) dominuje dla brył nieopływowych
  Opór tarciowy (lepki) dominuje dla brył opływowych i płaskich płyt

Symulacja oblicza F_d w czasie rzeczywistym podczas regulacji suwaków prędkości i gęstości płynu, z wzorem wyświetlanym liczbowo. Dostępnych jest sześć kształtów bryły, każdy z charakterystyczną krzywą C_d(Re) wyprowadzoną z danych eksperymentalnych.

Porównanie kształtów i wartości C_d

Przy pośrednich liczbach Reynoldsa (Re ≈ 10⁴–10⁵), sześć brył obejmuje niemal dwa rzędy wielkości współczynnika oporu:

Kształt bryły       C_d (Re ~ 10⁵)    Uwagi
─────────────────────────────────────────────────────────
Płaska płyta        ~1,17             Czysty opór ciśnieniowy, masywny ślad
Sześcian            ~1,05             Odrywanie na ostrych krawędziach
Cylinder            ~1,0              Szeroki ślad, odrywanie przy ~80°
Kula                ~0,47             Podkrytyczny (laminarna warstwa przyścienna)
Bryła opływowa      ~0,04–0,10        Przepływ przylegający, cienki ślad
Profil (NACA 0012)  ~0,006–0,012      Głównie opór tarciowy przy niskim kącie natarcia

Panel nakładki wielokształtowej wykreśla wszystkie sześć krzywych C_d(Re) jednocześnie na osi log-log, ułatwiając zobaczenie, gdzie krzywe się zbiegają, rozbiegają i gdzie występuje kryzys oporu.

Kryzys oporu

Najbardziej uderzającą cechą symulacji współczynnika oporu jest kryzys oporu dla kuli. Przy podkrytycznych liczbach Reynoldsa (Re < 3×10⁵), warstwa przyścienna jest laminarna i odrywa się wcześnie, dając szeroki ślad i C_d ≈ 0,47. Gdy Re przekracza wartość krytyczną, warstwa przyścienna przechodzi w turbulentną, zanim osiągnie punkt odrywania. Turbulentne warstwy przyścienne mają wyższy pęd i mogą pozostać przylegające dalej wokół kuli, drastycznie zwężając ślad:

Podkrytyczny (Re < 3×10⁵):
  Laminarna warstwa przyścienna → wczesne odrywanie (~80° od przedniego punktu stagnacji)
  Szeroki ślad, duży opór ciśnieniowy → C_d ≈ 0,47

Przejście krytyczne (Re ≈ 3×10⁵):
  Przejście warstwy przyściennej w turbulentną → opóźnione odrywanie (~120° od przodu)
  Wąski ślad, mały opór ciśnieniowy → C_d ≈ 0,10

Efekt w sporcie:
  Dołki w piłce golfowej: sztucznie wywołują turbulencję warstwy przyściennej przy niższym Re
  → kryzys oporu przy Re ≈ 4×10⁴ zamiast 3×10⁵
  → 2× dłuższy zasięg wobec gładkiej piłki przy typowych prędkościach golfowych

Wizualizacja odrywania przepływu rysuje punkt odrywania jako animowany znacznik na konturze bryły i szkicuje szerokość śladu jako zacieniony obszar, dając intuicyjne wyczucie tego, co liczba C_d oznacza fizycznie.

🔺 Kopuła geodezyjna — od dwudziestościanu do triangulacji sferycznej

Dwudziestościan jako punkt wyjścia

Regularny dwudziestościan ma 20 równobocznych ścian trójkątnych, 30 krawędzi i 12 wierzchołków, wszystkie na okręgu opisanym na sferze. Jest to bryła platońska o największej liczbie ścian, co czyni go naturalnym punktem wyjścia dla aproksymacji sfery. Współrzędne jego wierzchołków to:

Złoty podział φ = (1 + √5) / 2 ≈ 1,6180

12 wierzchołków (znormalizowanych do sfery jednostkowej):
  (0, ±1, ±φ)  →  znormalizowane: (0, ±1/√(1+φ²), ±φ/√(1+φ²))
  (±1, ±φ, 0)  →  znormalizowane podobnie
  (±φ, 0, ±1)  →  znormalizowane podobnie

Liczba ścian:  F = 20
Liczba krawędzi:  E = 30
Liczba wierzchołków: V = 12
Euler: V - E + F = 12 - 30 + 20 = 2  ✓

Podział o częstotliwości ν

Aby podzielić każdą ścianę trójkątną przy częstotliwości ν, każda krawędź oryginalnego dwudziestościanu jest dzielona na ν równych części, tworząc siatkę ν² mniejszych trójkątów na ścianę. Nowe wierzchołki wewnętrzne umieszczane są przez interpolację liniową, a następnie rzutowane na zewnątrz na sferę opisaną:

Dla każdego oryginalnego trójkąta (A, B, C) przy częstotliwości ν:
  Dla i = 0 do ν:
    Dla j = 0 do ν - i:
      k = ν - i - j
      Punkt P = (i·A + j·B + k·C) / ν    // interpolacja barycentryczna
      P_sphere = P / |P|                   // rzutowanie na sferę jednostkową

Liczba wyników:
  Nowych wierzchołków na ścianę:  (ν+1)(ν+2)/2
  Nowych trójkątów na ścianę: ν²
  Wierzchołków łącznie: V = 10ν² + 2
  Ścian łącznie:    F = 20ν²
  Krawędzi łącznie:   E = 30ν²
  Sprawdzenie Eulera:    V - E + F = 10ν² + 2 - 30ν² + 20ν² = 2  ✓ dla wszystkich ν

Wyświetlanie na żywo wzoru Eulera w symulacji aktualizuje się podczas przeciągania suwaka częstotliwości, potwierdzając V−E+F = 2 na każdym kroku — satysfakcjonująca kontrola spójności topologii.

Tryb kopuły i właściwości strukturalne

Tryb kopuły przycina podzieloną sferę na równiku, zachowując tylko dolną półkulę. Zmniejsza to w przybliżeniu o połowę liczbę wierzchołków, ale pierścień graniczny wymaga wzmocnienia zastrzałami. Symulacja oblicza maksymalną wariancję długości zastrzałów — kopuły geodezyjne zbudowane z pojedynczej częstotliwości mają wszystkie zastrzały niemal równej długości, co upraszcza budowę. Przy ν = 1 istnieje tylko jedna długość zastrzału; przy ν = 6 wciąż istnieje tylko 9 różnych długości dla całej struktury kopuły.

💍 Stabilne dopasowanie — Gale-Shapley i algorytm odroczonej akceptacji

Problem stabilnego dopasowania

Mając N proponentów i N akceptorów, każdy z pełnym ścisłym rankingiem preferencji względem drugiej strony, stabilne dopasowanie to dopasowanie doskonałe (każda osoba dopasowana dokładnie raz) bez pary blokującej: pary (p, a), gdzie proponent p preferuje akceptora a nad swoje obecne dopasowanie, a akceptor a preferuje p nad swoje obecne dopasowanie. Niestabilne dopasowanie w praktyce „rozpadłoby się” — para blokująca zdezerterowałaby.

Definicja formalna:
  Proponenci P = {p₁, …, pN}
  Akceptorzy A = {a₁, …, aN}
  Każdy pᵢ ma ranking: aσ(i,1) > aσ(i,2) > … > aσ(i,N)
  Każdy aⱼ ma ranking: pτ(j,1) > pτ(j,2) > … > pτ(j,N)

Para blokująca: (pᵢ, aⱼ) taka, że
  pᵢ preferuje aⱼ nad μ(pᵢ)   ORAZ
  aⱼ preferuje pᵢ nad μ⁻¹(aⱼ)
gdzie μ jest bieżącym dopasowaniem.

Odroczona akceptacja Gale-Shapleya

Algorytm (Gale i Shapley, 1962) działa w rundach. Każdy niedopasowany proponent proponuje swojemu najbardziej preferowanemu akceptorowi jeszcze nieodrzuconemu. Każdy akceptor tymczasowo przyjmuje najlepszą dotychczas otrzymaną propozycję i odrzuca wszystkie inne. Kluczowy wgląd polega na tym, że propozycje przesuwają się jedynie w dół listy preferencji proponenta, podczas gdy akceptorzy mogą jedynie „wymieniać się w górę” — więc algorytm zawsze się kończy i produkuje stabilne dopasowanie:

Gale-Shapley (optymalny dla proponentów):
  dopóki ∃ niedopasowany proponent pᵢ z niepustą listą:
    aⱼ = najlepszy pozostały akceptor pᵢ
    jeśli aⱼ jest niedopasowany:
      tymczasowo dopasuj (pᵢ, aⱼ)
    w przeciwnym razie jeśli aⱼ preferuje pᵢ nad obecnym posiadaczem pₖ:
      rozdopasuj pₖ (pₖ staje się ponownie wolny)
      tymczasowo dopasuj (pᵢ, aⱼ)
    w przeciwnym razie:
      aⱼ odrzuca pᵢ; pᵢ usuwa aⱼ z listy

Złożoność: O(N²) propozycji w najgorszym przypadku
Kończy się po co najwyżej N² rundach

Tryb krok po kroku symulacji podświetla bieżącego proponenta na żółto, docelowego akceptora na niebiesko, a odrzucone propozycje na czerwono. Licznik rund propozycji pokazuje, ile rund było potrzebnych, co waha się między N (najlepszy przypadek) a N² (najgorszy przypadek).

Optymalność i weryfikacja par blokujących

Gale-Shapley produkuje stabilne dopasowanie optymalne dla proponentów: każdy proponent otrzymuje najlepszego partnera, jakiego mógłby otrzymać w dowolnym stabilnym dopasowaniu. Akceptorzy otrzymują najgorszego stabilnego partnera. Detektor par blokujących symulacji wyczerpująco sprawdza wszystkie N² par po zbieżności, potwierdzając brak par blokujących. Test weryfikacyjny 2000 uruchomień losuje wszystkie listy preferencji i sprawdza stabilność każdej z nich — potwierdzając poprawność nawet dla adwersarialnych struktur preferencji.

🎭 Hologram — nagrywanie interferencyjne i rekonstrukcja czoła fali

Nagrywanie: interferencja na kliszy holograficznej

W przeciwieństwie do fotografii, która rejestruje jedynie natężenie, hologram rejestruje zarówno amplitudę, jak i fazę czoła fali obiektu poprzez interferowanie go ze spójną wiązką odniesienia. Wzór natężenia na kliszy holograficznej wynosi:

Geometria nagrywania:
  E_ref(x) = A_r · exp(i k x sin θ_r)     (płaska fala odniesienia)
  E_obj(x) = Σⱼ aⱼ / rⱼ · exp(i k rⱼ)    (fale sferyczne z punktów obiektu)

Zarejestrowane natężenie (transmitancja hologramu ∝ ekspozycja):
  I(x) = |E_ref + E_obj|²
       = |E_ref|² + |E_obj|² + E_ref* E_obj + E_ref E_obj*
         ───────────────────────────────────────────────────
         podkład DC   plamkowanie   człon sygnału   człon sprzężony

Człon sygnału E_ref* · E_obj koduje mapę fazy obiektu.

Symulacja renderuje I(x) jako wzór prążków w skali szarości na kliszy holograficznej w czasie rzeczywistym podczas przeciągania punktowych źródeł obiektu po scenie. Wiele punktów obiektu tworzy złożone nakładające się zestawy prążków, które razem kodują obiekt 3D.

Rekonstrukcja: rzędy dyfrakcji

Oświetlenie przetworzonego hologramu oryginalną wiązką odniesienia dyfrakcyjnie odczytuje zapisane czoło fali. Pole przepuszczone rozkłada się na trzy rzędy:

Transmitancja hologramu t(x) ∝ I(x)   (liniowy reżim nagrywania)

Zrekonstruowane pole = t(x) · E_ref(x)
  = (|E_ref|² + |E_obj|²) · E_ref    → rząd zerowy (DC, przejście proste)
  +  E_ref* · E_obj · E_ref           → rząd +1: E_obj zrekonstruowane  (obraz wirtualny)
  +  E_ref · E_obj* · E_ref           → rząd −1: fala sprzężona E_obj* (obraz rzeczywisty)

Rząd +1: fala rozbieżna wydaje się pochodzić z oryginalnych pozycji obiektu → obraz wirtualny
Rząd −1: fala zbieżna skupia się w obraz rzeczywisty (pseudoskopowy, odwrócona głębia)

Symulacja animuje rekonstrukcję, propagując fale z płaszczyzny hologramu. Obraz wirtualny rzędu +1 pojawia się za kliszą w prawidłowej pozycji 3D; obraz rzeczywisty rzędu −1 formuje się z przodu. Suwak reguluje długość fali rekonstrukcji, pokazując, jak niedopasowanie długości fali do długości fali nagrywania powoduje powiększenie i zniekształcenie chromatyczne.

🔵 Louvain — wykrywanie społeczności i modularność

Modularność jako cel

Modularność Q podziału grafu mierzy, o ile więcej krawędzi wewnątrzspołecznościowych istnieje w porównaniu z losowym modelem zerowym o tej samej sekwencji stopni:

Q = (1/2m) · Σᵢⱼ [Aᵢⱼ - kᵢkⱼ/(2m)] · δ(cᵢ, cⱼ)

gdzie:
  m     = całkowita liczba krawędzi
  Aᵢⱼ  = wpis macierzy sąsiedztwa
  kᵢ   = stopień węzła i
  cᵢ   = etykieta społeczności węzła i
  δ(a,b) = 1 jeśli a = b, 0 w przeciwnym razie

Zakres Q: [−0,5, 1)
  Q ~ 0: brak struktury społecznościowej (graf losowy)
  Q > 0,3: znacząca struktura społecznościowa
  Q > 0,7: bardzo silne społeczności

Dwufazowy algorytm Louvain

Algorytm Louvain (Blondel i in., 2008) naprzemiennie stosuje dwie fazy, aż nie jest możliwa dalsza poprawa:

Faza 1 — lokalna optymalizacja:
  Inicjalizacja: każdy węzeł we własnej społeczności
  Dla każdego węzła i (losowa kolejność):
    Oblicz ΔQ dla przeniesienia i do każdej sąsiedniej społeczności c:
      ΔQ = [Σ_in + 2k_i,in) / (2m) - ((Σ_tot + kᵢ) / (2m))²]
           − [Σ_in / (2m) − (Σ_tot / (2m))² − (kᵢ / (2m))²]
    gdzie Σ_in  = krawędzie wewnątrz c
          Σ_tot = całkowity stopień węzłów w c
          k_i,in = krawędzie od i do c
    Przenieś i do społeczności z maksymalnym ΔQ > 0
  Powtarzaj, aż żaden ruch nie zwiększa Q

Faza 2 — zgrubianie (agregacja super-węzłów):
  Zbuduj nowy graf: węzły = społeczności z fazy 1
  Wagi krawędzi = suma krawędzi między społecznościami
  Pętle własne = suma krawędzi wewnątrz każdej społeczności
  Powtórz fazę 1 na zgrubionym grafie

Złożoność: O(n log n) w praktyce na grafach rzadkich
Kończy się, gdy Q przestaje rosnąć między pełnymi przejściami

Generowanie grafu i weryfikacja

Symulacja wykorzystuje stochastyczny model blokowy do generowania grafów testowych o znanej rzeczywistej strukturze społecznościowej. SBM umieszcza n_{społeczności} grup o równym rozmiarze; prawdopodobieństwo krawędzi wewnątrz grupy p_in przewyższa prawdopodobieństwo między grupami p_out, tworząc zaplanowany podział. Po zbieżności Louvain symulacja oblicza znormalizowaną informację wzajemną (NMI) między wykrytymi społecznościami a rzeczywistością wzorcową i wyświetla ostateczną wartość Q. Dla dobrze rozdzielonych społeczności (p_in/p_out > 4) Louvain niezawodnie odzyskuje zaplanowaną strukturę z NMI > 0,95.

🔩 Struktura ziaren — krzepnięcie, zgrubianie i relacja Halla-Petcha

Nukleacja Woronoja

Krzepnięcie zaczyna się od nukleacji: zarodki krystaliczne powstają w losowych pozycjach w stopie. Każde ziarno rośnie na zewnątrz, aż spotka sąsiednie ziarno, tworząc granicę ziarna. Powstała tessellacja to diagram Woronoja — każda komórka ziarna zawiera wszystkie punkty bliższe swojemu jądru niż jakiemukolwiek innemu:

Komórka Woronoja dla jądra i:
  V(i) = { x ∈ ℝ² : |x − nᵢ| ≤ |x − nⱼ| dla wszystkich j ≠ i }

Granica ziarna między komórkami i oraz j:
  Symetralna odcinka nᵢnⱼ
  spotyka się w potrójnym złączu (trzy ziarna) pod kątem ~120°
  (minimalizuje energię granicy ziaren na jednostkę długości)

Symulacja zarodkuje ziarna jednostajnie losowo i renderuje tessellację Woronoja, przypisując każdemu ziarnu losową orientację krystalograficzną (kodowaną kolorystycznie odcieniem).

Wyżarzanie modelem Pottsa i wzrost ziaren

Wzrost ziaren napędzany jest krzywizną: granice ziaren migrują w stronę swojego środka krzywizny, redukując całkowitą długość granicy i energię. Model Pottsa dyskretyzuje to na siatce: każde miejsce sieciowe niesie spin (ID ziarna); krok Metropolisa Monte Carlo odwraca miejsca graniczne do orientacji sąsiedniego ziarna, jeśli redukuje to całkowitą energię granicy:

Energia modelu Pottsa:
  H = J · Σ_{⟨i,j⟩} (1 − δ(sᵢ, sⱼ))

  J = energia granicy ziaren na wiązanie
  Suma tylko po parach najbliższych sąsiadów
  δ(a,b) = 1 w tym samym ziarnie, 0 w poprzek granicy

Aktualizacja Metropolisa w miejscu i:
  Wybierz losową orientację sąsiedniego ziarna s'
  ΔH = H(s') − H(sᵢ)
  Zaakceptuj, jeśli ΔH ≤ 0, lub z prawdopodobieństwem exp(−ΔH / kT)

Prawo wzrostu ziaren (wzrost paraboliczny):
  ⟨D⟩² − ⟨D₀⟩² = K · t
  ⟨D⟩ ∝ tⁿ  z n ≈ 0,5 dla idealnego normalnego wzrostu ziaren

Relacja Halla-Petcha

Granica plastyczności metalu polikrystalicznego rośnie wraz ze zmniejszaniem się rozmiaru ziarna, ponieważ granice ziaren działają jako przeszkody dla ruchu dyslokacji. Relacja Halla-Petcha kwantyfikuje to:

Relacja Halla-Petcha:
  σ_y = σ₀ + k / √d

gdzie:
  σ_y = granica plastyczności (MPa)
  σ₀  = naprężenie tarcia sieciowego (granica plastyczności pojedynczego kryształu)
  k   = nachylenie Halla-Petcha (współczynnik umocnienia, MPa·m^½)
  d   = średnia średnica ziarna (m)

Przykładowe wartości (stal):
  σ₀ ≈ 70 MPa
  k  ≈ 0,74 MPa·m^½
  Przy d = 10 μm: σ_y ≈ 70 + 0,74/√(10×10⁻⁶) ≈ 70 + 234 = 304 MPa
  Przy d = 1  μm: σ_y ≈ 70 + 740 = 810 MPa

Symulacja śledzi średnią średnicę ziarna D(t) w miarę postępu wyżarzania i wykreśla σ_y vs D^−½ w czasie rzeczywistym obok mapy ziaren. Cztery predefiniowane materiały (żelazo czyste, stop aluminium, miedź, stal austenityczna) dostarczają realistyczne wartości σ₀ i k, dzięki czemu wykres Halla-Petcha ma fizycznie dokładne jednostki.

Co dalej

Fala 65 jest w trakcie tworzenia z kolejnymi sześcioma symulacjami obejmującymi działanie kapilarne (prawo Jurina i kształt menisku), przeszukiwanie ścieżek A* z trzema heurystykami, krzywe parametryczne od de Casteljau Béziera przez B-spliny aż po NURBS, wiązki wirowe Laguerre'a-Gaussa z orbitalnym momentem pędu, stabilność koloidalną przez teorię DLVO oraz protokoły trasowania internetowego, w tym Bellman-Forda i Dijkstry link-state. Biblioteka jest na dobrej drodze, aby przekroczyć 600 symulacji w czerwcu.