💨 Dym VFX — Stable Fluids

Stable Fluids to metoda opracowana przez Josa Stama (1999) do bezwarunkowo stabilnej symulacji płynów w czasie rzeczywistym. Wykorzystuje niejawne całkowanie czasowe dla dyfuzji oraz krok adwekcji semi-lagranżowskiej, który nigdy nie staje się numerycznie niestabilny, niezależnie od wielkości kroku czasowego.

Nieściśliwe równania Naviera-Stokesa opisują ruch płynu: ∂u/∂t = -(u·∇)u + ν∇²u + f, przy ∇·u = 0 dla nieściśliwości. Poszczególne wyrazy odpowiadają odpowiednio adwekcji (transportowi przez przepływ), dyfuzji (lepkości) i siłom zewnętrznym.

Adwekcja semi-lagranżowska śledzi każdą komórkę siatki wstecz w czasie wzdłuż pola prędkości, aby znaleźć pochodzenie płynu, a następnie interpoluje gęstość/prędkość w tym punkcie. Jest bezwarunkowo stabilna, ale wprowadza dyfuzję numeryczną, przez co przepływ z czasem wygląda na rozmyty.

Krok projekcji wymusza nieściśliwość (∇·u = 0), rozwiązując równanie Poissona dla ciśnienia i odejmując gradient ciśnienia od prędkości. Zapobiega to gromadzeniu się lub zanikaniu płynu, zachowując jego całkowitą ilość.

Wzmocnienie wirowości przeciwdziała dyfuzji numerycznej adwekcji semi-lagranżowskiej. Wykrywa obszary o dużej wirowości (rotacji prędkości) i przykłada siły skupiające ruch wirowy, dając bardziej turbulentne i realistyczne efekty dymu i ognia.

Dyfuzja numeryczna to sztuczne wygładzanie wprowadzane przez dyskretyzację równań adwekcji. Metody semi-lagranżowskie z interpolacją dwuliniową są na nią szczególnie podatne, przez co dym z czasem wygląda na rozmyty bez technik korekcyjnych, takich jak wzmocnienie wirowości.

Gęstość dymu jest skalarem pasywnym transportowanym przez pole prędkości. Pełny krok czasowy: (1) dodanie sił, (2) dyfuzja prędkości, (3) projekcja do pola bezźródłowego, (4) adwekcja prędkości, (5) ponowna projekcja, (6) adwekcja gęstości. Gęstość jest mnożona przez współczynnik zanikania, by symulować zanikanie dymu.

Symulacje Stable Fluids w czasie rzeczywistym zwykle używają siatek od 64×64 do 256×256 dla 60 kl./s w JavaScript. Implementacje na GPU (WebGL) mogą działać z siatką 512×512 lub większą. Wyższa rozdzielczość oddaje więcej szczegółów turbulencji, ale wymaga więcej obliczeń na klatkę.

Iteracja Gaussa-Seidla rozwiązuje układy liniowe wynikające z niejawnej dyfuzji i projekcji ciśnienia. Każda komórka jest aktualizowana z użyciem najnowszych wartości sąsiadów w powtarzanych przebiegach aż do zbieżności — zwykle 4-20 iteracji na klatkę dla dobrego kompromisu jakości i szybkości.

Wirowość ω = ∇×u jest obliczana w każdej komórce. Gradient |ω| wskazuje kierunek ku centrom wirów. Siła wzmacniająca f = ε(η×ω), gdzie η = ∇|ω|/|∇|ω||, popycha płyn do rotacji wokół istniejących wirów, przywracając rotację utraconą wskutek dyfuzji numerycznej.