SPH: Budowanie płynu z cząstek

Płyn z cząstek

Istnieją dwa sposoby symulacji płynu. Metoda Eulera umieszcza siatkę w przestrzeni i pyta, co przechodzi przez każdy komórkowy obszar; metoda Lagrangea umieszcza cząstki w płynie i pozwalają je na transport płynu. Hydrodynamika Przysmakiwana (SPH) to metoda Lagrange’a. Została zaproponowana w 1977 roku przez Gingold, Monaghan i (samodzielnie) Lucy dla astrofizyki, gdzie nie ma siatki do umieszczenia, a następnie przeniosła się do grafiki komputerowej i inżynierii, ponieważ obsługuje splaszczania, wolne powierzchnie i oddzielające się kropli bez żadnego kodu dla przypadków specjalnych.

Podstawowym pomysłem jest interpolaacja z użyciem jądra przysmakiwania. Jakiekolwiek pole A może być odtworzone w dowolnej punkcie x sumując wkład cząstek bliskich do tego punktu, każdy ważony przez jądro W, które maleje ze wzrostem odległości i znikne po przekroczeniu promienia przysmakiwania h:

A(x) = Σ_j m_j · (A_j / ρ_j) · W(x − x_j, h) Elegancja polega na tym, że pochodne pola stają się pochodne jądra — funkcji analitycznej, którą różniczkuje się raz na papierze. Brakuje tu siatki, ani terminu transportowego, ani numerycznego rozpraszania z powodu interpolaacji między komórkami; masa jest dokładnie zachowana, ponieważ przemieszcza się wraz z cząsteczkami.

A(x) = Σ_j  m_j · (A_j / ρ_j) · W(x − x_j, h)

Jądra

Jądro musi być normalizowane (integruje się do 1), mieć zbiór skończony (jest zerowe poza h, więc sumy są lokalne) i być wystarczająco gładkie, aby jego gradient był dobrze zachowany. W praktyce używa się wielu jąder dla różnych wyrażeń, o ile chodzi o formułę Müllera, Charypara i Grossa z 2003 roku, na której opierają się większość implementacji SPH w czasie rzeczywistym:

poly6 → gęstość gładkie, ekonomiczne, ale jego gradient zniknie w centrum spiky → siła ciśnienia gradient jest niezerowy dla r → 0, co uniemożliwia aglomerację cząstek viscosity → siła viskoznosc Laplacian jest pozytywny wszędzie, więc nie można wprowadzać energii Używanie poly6 do gradientu ciśnienia to klasyczny błąd początkujących: ponieważ jego gradient zniknie, gdy cząstki zbliżą się do siebie, siła odpychająca zniknie dokładnie wtedy, gdy jest ona potrzebna, i cząstki zaczną aglomerować. Jądro spikowe istnieje exactly to poprawić to. Jądro szpilkowatościa sześciennego jest drugim powszechnym wyborem i standardem w SPH astronomicznym.

poly6      → density        smooth, cheap, but its gradient
                             vanishes at the centre
spiky      → pressure force  non-zero gradient at r → 0, which is
                             what stops particles clumping
viscosity  → viscous force   Laplacian is positive everywhere, so it
                             cannot inject energy

Gęstość, ciśnienie i równanie Taita

Każde kroku zaczyna się od gęstości. Każda cząstka sumuje masę swoich sąsiadów przemnożoną przez jądro:

ρ_i = Σ_j m_j · W(x_i − x_j, h) // w tym przypadku uwzględnia się także samą cząstkę Następnie ciśnienie. Przegubowe przepływy wymagająbyłoby rozwiązywania globalnej równania Poissona na każdym kroku, co byłoby drogim procesem; słabo przegubowe SPH zamiast tego używa równania stanu, które nagradza kompresję tak silnie, że płyn pozostaje prawie nieprzegubowy. Standardowym równaniem jest równanie Taita:

p_i = B · ((ρ_i / ρ₀)^γ − 1) // γ = 7 dla wody B = ρ₀ · c² / γ // c = sztuczna prędkość dźwięku Wykładnik γ = 7 sprawia, że ciśnienie wzrasta brutalnie stromie, gdy gęstość przekracza gęstość spoczynkową ρ₀, co oznacza, że mała błąd w gęstości powoduje duży siły odzyskawcze. Stała sprężystości B ustawiana jest na podstawie sztucznej prędkości dźwięku c, wybranego nie jako rzeczywista 1500 m/s wody, ale około dziesięciokrotnie większej od najświeższego oczekiwanej przepływu — to pozwala utrzymać oscylacje gęstości pod około 1%, jednocześnie pozwalając na krok czasowy tysiące razy większy niż rzeczywista prędkość dźwięku. Prostej liniowej równania stanu, p = k(ρ − ρ₀), często wystarcza dla demonstracji w czasie rzeczywistym w celach wizualnych.

Siła ciśnienia na cząstce używa formy symetrycznej, która zachowuje dokładnie chwilową pęd, ponieważ siła z i na j jest równa i przeciwna:

ρ_i = Σ_j  m_j · W(x_i − x_j, h)      // include the particle itself

Gęstość i naprężenie powierzchniowe

Gęstość uwydatnia gładką pole prędkości: każda cząstka jest przyciągana ku średniej prędkości swoich sąsiadów, ważona przez Laplasjan jądra gęstości. Określa to, jak płynny wydaje się płyn, a jednocześnie cicho stabilizuje symulację, tłumacząc szum wysokofrekansowy cząsteczkowy, który tworzy termin napięcia ciśnienia. Astronomiczna SPH używa innego narzędzia – gęstości artystycznej (termin Monagha α–β), którego zadaniem jest pozwolenie na powstanie szoków bez przepływu cząstek przez siebie.

Naprężenie powierzchniowe wymaga małej maszyny, ponieważ cząstka wewnątrz nie wie, że znajduje się wewnątrz. Standardowa technika to pole kolorów: każdemu cząstce przypisuje wartość 1 i gładzi ją za pomocą jądra. Gradient jest prawie zerowy głęboko wewnątrz płynu i pokazuje na zewnątrz na wolnej powierzchni, gdzie suma sąsiadów jest jednostronna. Wartość gradientu działa jako detector powierzchniowy, jego kierunek to normalna do powierzchni, a odchylenie normalizowanej normalnej to krzywizna κ. Siła jest wtedy f = −σ · κ · n̂, która przyciąga powierzchnię do płaskości i jest to, co sprawia, że pętlę wolno upadająca staje się kuli, a małe kapcie pozostają w kształcie kulki.

Szukanie sąsiadów: część, która decyduje o szybkości klatek

Każda z tych sum biegnie po „sąsiadach wewnątrz h”. Wykonane prosto to jest skan O(n²) na każdym kroku i dominuje nad wszystkim innym. Rozwiązanie polega na użyciu uniformnej siatki hashowskiej z rozmiarem komórki dokładnie równym h: wtedy sąsiedztwo cząsteczki może znajdować się tylko w jej własnej komórce oraz 8 (2D) lub 26 (3D) sasiednich komórkach.

const c = 1 / h;                      // cell size = smoothing radius
const key = (x, y) => (Math.floor(x * c) * 92837111 ^
                       Math.floor(y * c) * 689287499) >>> 0;

// rebuild each step: counting sort into a flat bucket array
buckets.fill(0);
for (const p of particles) buckets[key(p.x, p.y) % NB]++;
// prefix-sum, scatter, then iterate only the 3×3 cell block per particle

Krok czasowy i warunek CFL

SPH jest explikite, więc krok jest ograniczony. Cząstka nie może przesunąć się dalej niż ułamek promienia rozsuwania w jednym kroku, a krok musi również rozwiązać prędkość sztuczna dźwięku i skali czasu viskozacji oraz sił. Standardowe ograniczenia CFL to:

Δt ≤ 0.25 · h / (c + v_max) // dźwięk / adwecja Δt ≤ 0.25 · sqrt(h / |a_max|) // przyspieszenie (np. grawitacja) Δt ≤ 0.125 · h² / ν // viskozacja rozsuwania Wybierz najmniejszą z tych trzech wartości, z uwzględnieniem czynnika bezpieczeństwa. Przewyższenie tych ograniczeń spowoduje wybuch wyrazu napięcia — który jest bardzo sztywny, dokładnie dlatego że γ = 7: cząstki się nadmiernie pokrywają, pojawiają się wzbudzenia gęstości, siła odzyskująca przekracza granice i symulacja eksploduje na strumień NaNów w kilku ramkach. To jest powodem, dla którego prędkość sztuczna dźwięku jest utrzymywana jak najszybsza dozwolona przez cel kompresyjności: c występuje w mianowniku pierwszego ograniczenia, więc każdy dodany czynnik do c kosztuje bezpośrednio rozmiar kroku. Integrator jest zazwyczaj Leapfrog dla tej samej przyczyny co w przypadku N-powłok — jest symplektyczny, kosztuje jedno obliczenie siły i tu siły są znacznie za drogie do oceny cztery razy na krok dla RK4.

Δt ≤ 0.25 · h / (c + v_max)           // sound / advection
Δt ≤ 0.25 · sqrt(h / |a_max|)         // acceleration (e.g. gravity)
Δt ≤ 0.125 · h² / ν                   // viscous diffusion

Często zadawane pytania

Wypróbuj na żywo

Wszystko powyżej działa bezpośrednio w Twojej przeglądarce — otwórz SPH Fluid i zmieniaj parametry podczas działania. Nic nie jest instalowane ani przesyłane na serwer, cały model działa w jednej karcie.