SPH: Створення Рідини з Частинок

Гідродинаміка з використанням ядра згладжування

Існують два способи моделювання рідини. Еулерів підхід передбачає розміщення сітки в просторі та запит про те, що проходить через кожен елемент; лагранговий підхід передбачає розміщення частинок у рідині та дозволяє їм нести її. Smoothed Particle Hydrodynamics є лагранговим методом. Він був винайдений у 1977 році Гінгольдом, Магонаном і (незалежно) Люсі для астрофізики, де немає коробки, в якій можна було б розмістити сітку, і пізніше перейшов у комп'ютерну графіку та інженерію, оскільки він обробляє бризки, вільні поверхні та розділяються краплі без будь-якого спеціального коду.

Центральна ідея полягає в інтерполяції з використанням ядра згладжування. Будь-яке поле A може бути відновлено в точці x шляхом підсумовування внесків від сусідніх частинок, кожне з яких вагується ядром W, яке падає зі збільшенням відстані та зникає за межами радіусу згладжування h:

A(x) = Σ_j m_j · (A_j / ρ_j) · W(x − x_j, h) Елегантність полягає в тому, що похідні поля стають похідними ядра — аналітична функція, яку можна продиференціювати один раз, на папері. Немає сітки, немає терму адивекції, немає чисельного дифузії від інтерполяції між клітинами; маса точно зберігається, оскільки вона їде на частинках.

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

Ядра (Kernels)

Ядро должно быть нормализовано (интегрироваться до 1), иметь компактную поддержку (нулевое значение за пределами h, так что суммы локальны) и быть достаточно гладким, чтобы его градиент был хорошо сформирован. В реальности используется несколько ядер для различных терминов, следуя формулировке 2003 года Муллера, Чарипара и Гросса, которой опираются большинство реальных SPH-реализаций:

poly6 → плотность гладкая, дешевая, но ее градиент исчезает в центре

spiky → сила давления не нулевой градиент при r → 0, что предотвращает скопление частиц

viscosity → сила вязкости Лапласиан положительный везде, поэтому он не может вводить энергию

Использование poly6 для градиента давления является классической ошибкой для начинающих: поскольку его градиент стремится к нулю по мере приближения частиц друг к другу, отталкивающая сила исчезает точно тогда, когда она нужна, и частицы сжимаются в кластеры. Ядро spiky существует именно для того, чтобы это исправить. Кубическая сплайн-ядро является другим распространенным выбором и является стандартом в астрофизическом SPH.

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

Щільність, тиск та рівняння Tait

Кожен крок починається зі щільності. Кожна частинка сумує маси сусідів, зважені ядром:

ρ_i = Σ_j m_j · W(x_i − x_j, h) // включаючи саму частинку Потім тиск. Істинно незмінюваний потік вимагав би розв’язання глобальної рівняння Пуассона на кожному кроці, що є дорогим; замість цього слабко здатний SPH використовує рівняння стану, яке жорстко карає стиснення, щоб рідина залишалася майже незмінюваною. Рівняння Tait є стандартним:

p_i = B · ((ρ_i / ρ₀)^γ − 1) // γ = 7 для води B = ρ₀ · c² / γ // c = штучна швидкість звуку Піднесеність γ = 7 призводить до різкого підвищення тиску, коли щільність перевищує щільність спокою ρ₀, тому невелика похибка щільності створює велику силу відновлення. Постійна жорсткості B встановлюється відносно штучної швидкості звуку c, обраної не як справжні 1500 м/с води, але приблизно в десять разів швидше за очікувану максимальну швидкість потоку — це дозволяє щільності коливатися менше ніж на 1%, при цьому час кроку тисячу разів більший, ніж час, необхідний для звуку.

Сила тиску на частинку використовує симетричну форму, яка точно зберігає імпульс, оскільки сила від i до j рівна і протилежна за знаком:

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

В’язкість та поверхневий натяг

В’язкість згладжує поле швидкості: кожне частинка тягнеться до середньої швидкості її сусідів, зваженої лапласіаном ядра в’язкості. Це визначає, наскільки густим виглядає рідина, і водночас тихо стабілізує симуляцію, гамуючи високочастовий шум частинок, що створюється силою тиску. Астрофізичний SPH використовує інше пристрій – штучну в’язкість (термін α–β Монгана), яка дозволяє утворюватися ударним хвилям без перекриття частинками.

Поверхневий натяг потребує певних механізмів, оскільки частинка всередині не знає, що вона знаходиться всередині. Стандартний трюк – поле кольору: призначте кожній частинці значення 1 і згладьте його за допомогою ядра. Його градієнт майже нульовий глибоко всередині рідини та вказує назовні на вільному поверхні, де сума сусідів є односторонньою. Величина градієнту діє як детектор поверхні, а напрямок – нормалі поверхні, а дивергенція нормалізованої нормалі – кривина (κ). Сила тоді обчислюється як f = −σ · κ · n̂, яка вирівнює поверхню та змушує падаючий кульок набувати сферичної форми, а маленькі краплі залишаються округлими.

Пошук сусідів: частина, що визначає швидкість кадрів

Кожне з цих обчислень проходить через ‘сусідів у межах h’’. Якщо робити це наївним способом, то це O(n²) перебір, і це домінує в усьому. Виправлення – використання рівномірної просторової сітки комірок із розміром h; тоді сусіди частинки можуть знаходитися лише в її власній комірці та 8 (2D) або 26 (3D) сусідніх комірках.

const c = 1 / h; // розмір комірки = радіус згладжування const key = (x, y) => (Math.floor(x * c) * 92837111 ^ Math.floor(y * c) * 689287499) >>> 0; // перебудовувати кожен крок: сортування підрахунком у плоский масив ячеек buckets.fill(0); for (const p of particles) buckets[key(p.x, p.y) % NB]++; // префіксний суматор, розсіювання, потім ітерувати лише блок 3х3 навколо частинки Перебудова сітки займає O(n) з використанням сортування підрахунком, а запит на сусідів стає O(k), де k – середній кількість сусідів — зазвичай близько 20 у 2D та 30–40 у 3D шляхом налаштування h проти відстані між частинками. Занадто мало сусідів робить оцінку щільності шумною; занадто багато, і кожен крок стає квадратично дорожчим з h. Зберігання частинок відсортованих за коміркою також вирішує іншу половину проблеми продуктивності: локальність пам’яті, оскільки сусідні частинки знаходяться поруч один з одним в масиві.

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

Крок часу та умова CFL

SPH є явною, тому крок обмежений. Частинка не повинна рухатися далі ніж певна частка радіуса згладжування за один крок, і крок також повинен вирішувати штучну швидкість звуку та часові масштаби в'язкого та силового тиску. Стандартні межі CFL такі:

Δt ≤ 0.25 · h / (c + v_max) // звук / адвекція Δt ≤ 0.25 · sqrt(h / |a_max|) // прискорення (наприклад, гравітація) Δt ≤ 0.125 · h² / ν // в'язке дифузія Візьміть мінімум з трьох, з урахуванням коефіцієнта запасу. Порушення цих обмежень призводить до того, що термін тиску — який є дуже жорстким, особливо через γ = 7 — вибухає: частинки перекриваються, виникають піки щільності, сила відновлення перевищується, і симуляція перетворюється на розпилювання NaNs протягом кількох кадрів. Саме тому штучна швидкість звуку підтримується якомога нижчою, щоб відповідати цільовій стислимості: c з'являється в чисельнику першого обмеження, тому кожен доданий до нього фактор безпосередньо впливає на розмір кроку. Інтегратор зазвичай Leapfrog з тих же самих причин, що й у випадку N-тіла — він сінгулярний, оцінює один силовий вираз, і сили тут настільки дорогі, щоб їх оцінювати чотири рази на крок для 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

Часті запитання

Спробуйте наживо

Усе, що вище, працює прямо у вашому браузері — відкрийте SPH Fluid і змінюйте параметри під час роботи. Нічого не встановлюється, нічого не завантажується на сервер, уся модель живе в одній вкладці.