Проблема симуляції рідин
Поведінка рідини описується рівняннями Нав'є-Стокса — набором диференціальних рівнянь у частинних похідних, що пов'язують швидкість, тиск, густину та в'язкість у неперервному середовищі. Проблема — саме слово неперервному. Справжня вода містить приблизно 3 × 1025 молекул на літр. Симулювати кожну з них неможливо.
Для розв'язання цієї проблеми існують дві широкі стратегії. Сітчасті методи (ейлерів підхід) ділять простір на фіксовані комірки й відстежують, як рідина тече між ними. Частинкові методи (лагранжів підхід) рухаються разом із рідиною, відстежуючи пакунки маси в міру їхнього переміщення в просторі.
Метод згладжених частинок гідродинаміки (SPH) — найпоширеніший частинковий метод. Він був винайдений у 1977 році Люсі та Гінголдом і Монаганом для астрофізичних симуляцій, а згодом адаптований для гідродинаміки. Його ключова перевага — природна обробка вільних поверхонь — меж між рідиною і порожнім простором, які, як відомо, надзвичайно складні для сітчастих методів.
Ядро SPH: розподіл точкових мас в об'єми
Основний математичний прийом у SPH — це функція ядра W(r, h), також звана ядром згладжування. Кожна частинка має довжину згладжування h — радіус впливу. Ядро визначає, як властивість, що несеться в точці, розподіляється по навколишньому об'єму.
Ядро повинне задовольняти кілька властивостей:
- Нормалізація: інтеграл по всьому простору дорівнює 1.
- Компактний носій: воно дорівнює нулю поза радіусом 2h, тому віддалені частинки не роблять жодного внеску (критично для продуктивності).
- Симетричність: W(r, h) залежить лише від відстані |r|, а не від напрямку.
- Збіжність до дельта-функції: при h → 0 W наближається до дельта-функції Дірака.
Кубічне сплайн-ядро, запроваджене Монаганом, — найпоширеніший вибір. Воно гладке, обчислювально дешеве і добре поводиться чисельно. Ядра Poly6 і Spiky, запроваджені Мюллером та ін. у їхній знаковій статті 2003 року про інтерактивний SPH, використовуються спеціально для обчислення густини й тиску відповідно — ключова ідея, що суттєво покращує стабільність симуляції.
Оцінка величин за частинками
Маючи ядро, можна оцінити будь-яку скалярну чи векторну величину A в точці r, підсумовуючи внески від усіх сусідніх частинок j:
A(r) = Σⱼ mⱼ (Aⱼ / ρⱼ) W(r − rⱼ, h)
Де mⱼ — маса частинки j, Aⱼ — значення A, яке несе ця частинка, ρⱼ — її густина, а W — ядро. Сама густина оцінюється так само:
ρ(r) = Σⱼ mⱼ W(r − rⱼ, h)
Похідні обчислюються аналітично шляхом диференціювання ядра — скінченні різниці не потрібні. Градієнт A у частинці i:
∇A(rᵢ) = Σⱼ mⱼ (Aⱼ / ρⱼ) ∇W(rᵢ − rⱼ, h)
Це і є математичний механізм, який перетворює неперервні рівняння Нав'є-Стокса на дискретну систему частинок.
Сила тиску
Рідина чинить опір стисканню. Коли частинки надто зближуються, густина зростає вище рівноважного значення ρ₀, і відновлювальна сила тиску розштовхує їх. Це обчислюється через рівняння стану — найпростіше з них це закон ідеального газу:
p = k(ρ − ρ₀)
Де k — стала жорсткості. Більше k робить рідину більш нестисливою, але може спричинити чисельну нестабільність (частинки різко «підстрибують»). Справжній нестисливий SPH використовує складніші розв'язувачі тиску (PCISPH або DFSPH), які ітеративно забезпечують нульову дивергенцію — ціною вищої обчислювальної вартості.
Сила тиску на частинку i від сусідів j дорівнює:
fᵢ_pressure = −mᵢ Σⱼ mⱼ ((pᵢ + pⱼ) / (2ρⱼ)) ∇W(rᵢ − rⱼ, h)
Симетризований член тиску (pᵢ + pⱼ)/2 важливий — він забезпечує виконання третього закону Ньютона і збереження повного імпульсу системи.
Сила в'язкості
Справжні рідини чинять опір зсуву — це і є в'язкість. Без неї частинки рідини, що проходять поруч одна з одною, не передають жодної енергії, і симуляція виглядає радше як бульбашки, ніж як вода. В'язкісний член SPH додає демпфувальну силу, пропорційну різниці швидкостей між сусідніми частинками:
fᵢ_viscosity = μ Σⱼ mⱼ ((vⱼ − vᵢ) / ρⱼ) ∇²W(rᵢ − rⱼ, h)
Де μ — коефіцієнт динамічної в'язкості, а ∇²W — лапласіан ядра. Мюллер та ін. запропонували спеціальне ядро для цього члена, чий лапласіан завжди додатний — це уникає нестабільностей зі зміною знаку, які турбували попередні реалізації.
Поверхневий натяг
Вода утворює краплі, тому що внутрішні молекули притягуються однаково в усіх напрямках, тоді як поверхневі молекули відчувають сумарну силу, спрямовану всередину. У SPH це моделюється за допомогою кольорового поля — скаляра, що дорівнює 1 усередині рідини й 0 поза нею. Градієнт цього поля спрямований до поверхні рідини; його кривизна визначає силу поверхневого натягу.
Поверхневий натяг — те, що робить симуляції рідини SPH красивими: вода збирається в краплі, краплі зливаються, а тонкі плівки природно розриваються. Без нього частинки рідини розсіюються в безформну хмару.
Чому SPH добре підходить для течій із вільною поверхнею
Сітчасті методи відстежують, які комірки сітки містять рідину, а які порожні. У міру руху поверхні комірки потрібно класифікувати, заповнювати й спорожняти — завдання обліку, що вимагає спеціальних методів, таких як Volume of Fluid (VOF) або підходи Level Set.
У SPH немає сітки. Частинки і є рідиною. Куди рухаються частинки, туди природно й слідує поверхня рідини. Розбризкування, злиття, розпад, утворення крапель — усе це виникає автоматично з динаміки частинок без будь-якого явного відстеження поверхні.
Це робить SPH методом вибору для:
- Симуляцій прориву дамби — класичного тесту, де стовп води обвалюється і затоплює резервуар.
- Океанських бризок і хвиль, що ламаються — де сітчасті методи мають труднощі з дрібномасштабною динамікою крапель.
- Взаємодії рідини й конструкції — частинки природно взаємодіють із твердими тілами й деформованими твердими об'єктами.
- Промислового виробництва — литво під тиском, лиття, потоки покриття.
Подивіться на SPH у дії прямо у браузері — інтерактивна симуляція на /fluid/ запускає тисячі частинок у реальному часі. Спробуйте наливати рідину під різними кутами або додавати перешкоди, щоб побачити, як потік розділяється і возз'єднується.
Виклик продуктивності
Наївний алгоритм SPH перевіряє кожну частинку відносно кожної іншої, щоб знайти сусідів у межах радіуса 2h. Це O(n²) роботи — подвоєння кількості частинок вчетверо збільшує обчислення.
На практиці всі реалізації SPH реального часу використовують просторове хешування: простір ділиться на комірки розміром 2h, і кожна частинка призначається комірці. Запити сусідів потребують перевірки лише 27 навколишніх комірок (у 3D). Це знижує складність до O(n) у середньому для рівномірно розподілених частинок.
Навіть так, тисячі частинок на кадр близькі до межі можливостей CPU для симуляції реального часу. Сучасний високоточний SPH використовує обчислювальні шейдери GPU, щоб виконувати пошук сусідів і обчислення сил паралельно на тисячах потоків GPU одночасно, дозволяючи мільйонам частинок працювати в інтерактивному темпі.