Wewnątrz SPH: algorytm, który napędza naszą symulację płynu

Gdy przeciągasz kursorem po naszej symulacji płynu i obserwujesz, jak ciecz faluje, chlupie i pryska, widzisz działającą w twojej przeglądarce hydrodynamikę wygładzonych cząstek (Smoothed Particle Hydrodynamics, SPH). SPH została wynaleziona w 1977 roku niezależnie przez Lucy'ego oraz przez Gingolda i Monaghana do symulowania przepływów astrofizycznych. Dziś stanowi podstawę efektów płynów w czasie rzeczywistym w grach, potokach efektów wizualnych i symulacjach inżynierskich. Oto matematyka, która za tym stoi.

Podstawowa idea: cząstki jako próbki płynu

Klasyczna obliczeniowa dynamika płynów rozwiązuje równania Naviera-Stokesa na stałej siatce. Każda komórka siatki ma gęstość, prędkość i ciśnienie, a równania opisują, jak płyn przepływa między sąsiednimi komórkami. Sprawdza się to dobrze dla przepływów pozostających w stałej domenie — na przykład opływu powietrza wokół skrzydła — ale staje się złożone, gdy sama powierzchnia płynu się porusza: fale, plamy, krople.

SPH przyjmuje inne podejście. Zamiast stałej siatki, płyn jest reprezentowany przez dużą liczbę cząstek, z których każda niesie masę, pozycję, prędkość i inne wielkości. Ciągłe pole płynu jest rekonstruowane z cząstek w dowolnym punkcie przestrzeni przez ważone uśrednianie po pobliskich cząstkach. Cząstki poruszają się wraz z płynem, automatycznie obsługując swobodne powierzchnie — gdzie kończy się płyn, a zaczyna powietrze — bez żadnego specjalnego traktowania.

To lagranżowskie (podążające za cząstkami) podejście jest kluczową przewagą SPH nad metodami eulerowskimi (zakotwiczonymi na siatce) dla przepływów ze swobodną powierzchnią: powierzchnia jest śledzona niejawnie przez to, gdzie znajdują się cząstki.

Funkcje jądra: matematyka wygładzania

„Wygładzenie" w SPH pochodzi od funkcji jądra W(r, h) — dzwonowatej funkcji, która przypisuje wagę każdej sąsiedniej cząstce na podstawie odległości r i długości wygładzania h. Jądro musi być znormalizowane (całkuje się do 1), symetryczne (zależy tylko od |r|) i zanikać poza pewną odległością odcięcia (zwykle 2h).

Najczęściej używanym jądrem w SPH czasu rzeczywistego jest jądro Poly6 do estymacji gęstości i jądro Spiky do sił ciśnienia. Jądro Spiky, zaproponowane przez Müllera i współautorów w ich artykule SIGGRAPH z 2003 roku, który wprowadził SPH do gier i grafiki komputerowej, ma niezerowy gradient w samej cząstce — kluczowy do obliczania odpychających sił ciśnienia, które zapobiegają zapadaniu się cząstek razem.

Długość wygładzania h jest zazwyczaj ustawiana tak, by obejmować 20–50 cząstek w 2D. Większe h oznacza płyn bardziej gładki, ale mniej szczegółowy; mniejsze h zachowuje drobne cechy, ale wymaga większej liczby cząstek dla stabilności.

Gęstość i ciśnienie: krok po kroku

Każdy krok symulacji obejmuje sekwencję obliczeń dla każdej cząstki:

• 1 Znajdź sąsiadów. Dla każdej cząstki i, zbierz wszystkie cząstki j w odległości 2h. Siatka haszująca przestrzennie sprawia, że jest to O(N), a nie O(N²).
• 2 Oblicz gęstość. ρ_i = Σ_j m_j W(|r_i − r_j|, h) — ważona suma mas sąsiednich cząstek.
• 3 Oblicz ciśnienie. Używając równania stanu: P_i = k(ρ_i − ρ₀), gdzie k to stała sztywności, a ρ₀ to gęstość spoczynkowa.
• 4 Oblicz siłę ciśnienia. F_pressure,i = −Σ_j m_j (P_i + P_j) / (2ρ_j) ∇W(|r_i − r_j|, h)
• 5 Oblicz siłę lepkości. F_visc,i = μ Σ_j m_j (v_j − v_i) / ρ_j ∇²W(|r_i − r_j|, h)
• 6 Całkuj. Dodaj grawitację, zastosuj całkowitą siłę, zaktualizuj prędkość i pozycję za pomocą integratora leapfrog lub symplektycznego Eulera.

Siły ciśnienia i równanie stanu

Prawdziwa woda jest niemal nieściśliwa — jej gęstość ledwo się zmienia pod wpływem ciśnienia. SPH z prostym równaniem stanu P = k(ρ − ρ₀) jest słabo ściśliwe: pozwala na niewielkie fluktuacje gęstości w zamian za uniknięcie kosztownego rozwiązywania, którego wymaga ścisła nieściśliwość. Parametr sztywności k musi być na tyle wysoki, by utrzymać zmienność gęstości poniżej kilku procent, ale nie na tyle wysoki, by krok czasowy musiał być niepraktycznie mały dla zachowania stabilności symulacji.

Dla dokładniejszego, nieściśliwego SPH — stosowanego w zastosowaniach inżynierskich — algorytmy takie jak IISPH (Implicit Incompressible SPH) czy DFSPH (Divergence-Free SPH) rozwiązują ciśnienie niejawnie w każdym kroku czasowym. Zbiegają się do niemal zerowego odchylenia gęstości, ale wymagają iteracyjnych solverów, które są zbyt wolne dla renderowania w czasie rzeczywistym.

Napięcie powierzchniowe i lepkość

Napięcie powierzchniowe w SPH jest modelowane albo metodą pola koloru (obliczanie normalnej interfejsu na podstawie gradientów gęstości i zastosowanie siły opartej na krzywiźnie), albo siłami kohezji między cząstkami. Nasza symulacja używa uproszczonego modelu kohezji, który tworzy charakterystyczną tendencję kropel wody do ściągania się z powrotem w kule.

Lepkość tłumi względną prędkość między sąsiednimi cząstkami — wyraz siły lepkości sprawia, że szybko poruszające się cząstki spowalniają sąsiadów i odwrotnie. Wysoka lepkość daje przepływ przypominający miód; zerowa lepkość (przepływ nielepki) daje ostre, nierealistyczne nieciągłości prędkości. Jądro laplasjanowe używane dla lepkości — z jego dodatnią drugą pochodną — zapewnia, że siła lepkości jest zawsze dyssypatywna, zapobiegając wprowadzaniu energii.

Dlaczego SPH bije metody siatkowe dla swobodnych powierzchni

Decydującą przewagą SPH dla przepływów ze swobodną powierzchnią jest brak konieczności śledzenia interfejsu. W metodach siatkowych swobodna powierzchnia musi być jawnie śledzona metodami takimi jak VOF (Volume of Fluid) czy Level Set. Dodają one znaczącą złożoność i mogą powodować numeryczną dyfuzję, która z czasem rozmywa ostre interfejsy.

Powierzchnia SPH jest wszędzie tam, gdzie są cząstki. Krople naturalnie się oddzielają, łączą i odbijają, ponieważ cząstki po prostu podążają za siłami. To sprawia, że SPH jest dominującą metodą w symulacji płynów czasu rzeczywistego dla gier i aplikacji interaktywnych, a także silnym konkurentem w offline'owych efektach wizualnych — symulacja oceanu w wielu współczesnych filmach to SPH lub hybrydowa metoda SPH-siatka.

Kompromisem jest to, że SPH ma trudności z cienkimi warstwami płynu i turbulentnymi przepływami o niskiej lepkości, gdzie metody siatkowe błyszczą. Prawdziwe produkcyjne symulatory płynów często wykorzystują podejścia hybrydowe — SPH blisko powierzchni, gdzie liczy się śledzenie swobodnej powierzchni, metody siatkowe w głębi, gdzie nieściśliwość musi być egzekwowana precyzyjnie.