Jak działa symulacja płynów SPH: matematyka stojąca za cząstkami płynu

Problem z symulowaniem płynów

Zachowanie płynu jest rządzone równaniami Naviera-Stokesa — zestawem cząstkowych równań różniczkowych wiążących prędkość, ciśnienie, gęstość i lepkość w ośrodku ciągłym. Problem tkwi w słowie ciągłym. Prawdziwa woda zawiera około 3 × 10²⁵ cząsteczek na litr. Symulowanie każdej z nich jest niemożliwe.

Istnieją dwie szerokie strategie radzenia sobie z tym. Metody siatkowe (podejścia eulerowskie) dzielą przestrzeń na stałe komórki i śledzą, jak płyn przepływa między nimi. Metody cząsteczkowe (podejścia lagrangeowskie) poruszają się razem z płynem, śledząc porcje masy, gdy przemieszczają się przez przestrzeń.

Hydrodynamika wygładzonych cząstek (SPH) to najszerzej stosowana metoda cząsteczkowa. Została pierwotnie wynaleziona w 1977 roku przez Lucy'ego i Gingolda & Monaghana do symulacji astrofizycznych, a później zaadaptowana do dynamiki płynów. Jej kluczową zaletą jest naturalna obsługa powierzchni swobodnych — granic między płynem a pustą przestrzenią — które są notorycznie trudne dla metod siatkowych.

Jądro SPH: rozprowadzanie mas punktowych na objętości

Główną matematyczną sztuczką w SPH jest funkcja jądra W(r, h), zwana też jądrem wygładzającym. Każda cząstka ma długość wygładzania h — promień wpływu. Jądro determinuje, jak właściwość niesiona w punkcie jest rozprowadzana po otaczającej objętości.

Jądro musi spełniać kilka właściwości:

• Normalizacja: całkuje się do 1 po całej przestrzeni.
• Nośnik zwarty: jest zerowe poza promieniem 2h, więc odległe cząstki nie wnoszą nic (kluczowe dla wydajności).
• Symetria: W(r, h) zależy tylko od odległości |r|, nie od kierunku.
• Zbieżność deltowa: gdy h → 0, W dąży do funkcji delta Diraca.

Jądro sześciennego splajnu, wprowadzone przez Monaghana, jest najczęstszym wyborem. Jest gładkie, tanie obliczeniowo i dobrze zachowuje się numerycznie. Jądra Poly6 i Spiky, wprowadzone przez Müllera i in. w ich przełomowej pracy z 2003 roku o interaktywnym SPH, są używane specjalnie odpowiednio do obliczeń gęstości i ciśnienia — kluczowy wgląd, który dramatycznie poprawia stabilność symulacji.

Szacowanie wielkości z cząstek

Mając jądro, możesz oszacować dowolną wielkość skalarną lub wektorową A w punkcie r, sumując wkłady od wszystkich pobliskich cząstek j:

A(r) = Σⱼ mⱼ (Aⱼ / ρⱼ) W(r − rⱼ, h)

Gdzie mⱼ to masa cząstki j, Aⱼ to wartość A niesiona przez tę cząstkę, ρⱼ to jej gęstość, a W to jądro. Sama gęstość jest szacowana w ten sam sposób:

ρ(r) = Σⱼ mⱼ W(r − rⱼ, h)

Pochodne są obliczane analitycznie przez różniczkowanie jądra — bez potrzeby różnic skończonych. Gradient A w cząstce i wynosi:

∇A(rᵢ) = Σⱼ mⱼ (Aⱼ / ρⱼ) ∇W(rᵢ − rⱼ, h)

To jest matematyczna maszyneria, która przekłada ciągłe równania Naviera-Stokesa na dyskretny system cząsteczkowy.

Siła ciśnienia

Płyn opiera się kompresji. Gdy cząstki zbliżają się do siebie zbyt mocno, gęstość rośnie powyżej wartości spoczynkowej ρ₀, a siła przywracająca ciśnienie rozpycha je z powrotem. Jest to obliczane poprzez równanie stanu — najprostsze to prawo gazu doskonałego:

p = k(ρ − ρ₀)

Gdzie k to stała sztywności. Wyższe k czyni płyn bardziej nieściśliwym, ale może powodować niestabilność numeryczną (cząstki gwałtownie odbijają się). Prawdziwie nieściśliwe SPH używa bardziej wyrafinowanych solverów ciśnienia (PCISPH lub DFSPH), które iteracyjnie wymuszają zerową dywergencję, kosztem wyższej złożoności obliczeniowej.

Siła ciśnienia działająca na cząstkę i od jej sąsiadów j wynosi:

fᵢ_pressure = −mᵢ Σⱼ mⱼ ((pᵢ + pⱼ) / (2ρⱼ)) ∇W(rᵢ − rⱼ, h)

Symetryzowany człon ciśnienia (pᵢ + pⱼ)/2 jest ważny — zapewnia zachowanie trzeciej zasady dynamiki Newtona oraz zachowanie całkowitego pędu układu.

Siła lepkości

Prawdziwe płyny opierają się ścinaniu — to jest lepkość. Bez niej cząstki płynu przechodzące blisko siebie nie przekazują energii, a symulacja wygląda jak bąbelki, a nie woda. Człon lepkości SPH dodaje siłę tłumiącą proporcjonalną do różnicy prędkości między sąsiednimi cząstkami:

fᵢ_viscosity = μ Σⱼ mⱼ ((vⱼ − vᵢ) / ρⱼ) ∇²W(rᵢ − rⱼ, h)

Gdzie μ to współczynnik lepkości dynamicznej, a ∇²W to Laplasjan jądra. Müller i in. zaproponowali specjalne jądro dla tego członu, którego Laplasjan jest zawsze dodatni — unikając niestabilności związanych ze zmianą znaku, które dręczyły wcześniejsze implementacje.

Napięcie powierzchniowe

Woda tworzy krople, ponieważ molekuły wewnętrzne są przyciągane równomiernie we wszystkich kierunkach, podczas gdy molekuły powierzchniowe odczuwają netto skierowane do wewnątrz ciągnięcie. W SPH jest to modelowane za pomocą pola koloru — skalara, który wynosi 1 wewnątrz płynu i 0 na zewnątrz. Gradient tego pola wskazuje w kierunku powierzchni płynu; jego krzywizna determinuje siłę napięcia powierzchniowego.

Napięcie powierzchniowe jest tym, co czyni symulacje płynów SPH piękne: woda tworzy krople, krople się łączą, a cienkie warstwy naturalnie się rozrywają. Bez niego cząstki płynu rozpraszają się w bezkształtną chmurę.

Dlaczego SPH jest dobre do przepływów z powierzchnią swobodną

Metody siatkowe śledzą, które komórki siatki zawierają płyn, a które są puste. Gdy powierzchnia się porusza, komórki muszą być klasyfikowane, wypełniane i opróżniane — wyzwanie księgowe, które wymaga specjalnych metod, takich jak Volume of Fluid (VOF) czy podejścia Level Set.

SPH nie ma siatki. Cząstki są płynem. Gdziekolwiek cząstki idą, powierzchnia płynu naturalnie podąża za nimi. Rozpryskiwanie, łączenie się, rozrywanie, formowanie kropel — wszystko to wyłania się automatycznie z dynamiki cząstek bez jakiegokolwiek jawnego śledzenia powierzchni.

To czyni SPH metodą z wyboru dla:

• Symulacji pęknięcia tamy — klasycznego benchmarku, w którym kolumna wody zapada się i zalewa basen.
• Piany oceanicznej i łamiących się fal — gdzie metody siatkowe zmagają się z drobnoskalową dynamiką kropel.
• Interakcji płyn-struktura — cząstki naturalnie oddziałują z ciałami sztywnymi i odkształcalnymi.
• Produkcji przemysłowej — formowanie wtryskowe, odlewanie, przepływy powlekające.

Wyzwanie wydajnościowe

Naiwny algorytm SPH sprawdza każdą cząstkę względem każdej innej cząstki, by znaleźć sąsiadów w promieniu 2h. To praca O(n²) — podwojenie liczby cząstek poczwórnie zwiększa obliczenia.

W praktyce wszystkie implementacje SPH w czasie rzeczywistym używają haszowania przestrzennego: dziel przestrzeń na komórki o rozmiarze 2h i przypisz każdą cząstkę do komórki. Zapytania o sąsiadów muszą sprawdzać tylko 27 otaczających komórek (w 3D). To redukuje złożoność do O(n) przeciętnie dla jednorodnie rozłożonych cząstek.

Mimo to tysiące cząstek na klatkę są bliskie limitowi CPU dla symulacji w czasie rzeczywistym. Nowoczesne wysokowierne SPH używa shaderów obliczeniowych GPU, by uruchamiać przeszukiwanie sąsiadów i obliczanie sił równolegle na tysiącach wątków GPU jednocześnie, umożliwiając miliony cząstek w tempie interaktywnym.