3. SPH
Интро
SPH (Smoothed Particle Hydrodynamics) – метод сглаженных частиц для гидродинамики. Супер древний метод, который придумали в 1977 году для астрофизики. Астрофизика нас не особо интересует, но метод оказался очень крутым и его стали использовать для симуляции жидкостей в компьютерной графике.
Основная идея заключается в том, что мы моделируем жидкость с помощью частиц. Каждая частица представляет собой небольшую порцию жидкости и взаимодействует с другими частицами в своем окружении. Для того чтобы придать им ощущение жидкости, нужно вычислить плотность и давление в каждой частице, а затем использовать эти значения для вычисления сил, действующих на частицы. Проблема в том, что эти величины не определены для отдельных частиц. Они должны быть непрерывными. Основная идея заключается в использовании сглаживающей функции (kernel function), которая определяет, как частицы влияют на плотность и давление в своей окрестности.
Каждая частица образует вокруг себя “облако” влияния, и чем ближе к частице, тем оно сильнее.
Функция распространения влияния называется kernel function. Она должна удовлетворять нескольким свойствам:
- Влияние должно быть ограничено радиусом сглаживания. Вообще, влияние частицы может быть и бесконечным, просто это будет выглядеть странно, если частица будет влиять на все.
- Быть убывающей с расстоянием: влияние частицы должно уменьшаться с увеличением расстояния от нее. Нам ничего не мешает сделать так, чтобы частица влияла одинаково на все расстояния, но тогда мы получим проблемы на границе радиуса сглаживания.
- Быть нормированной: интеграл по всему пространству должен быть равен 1. Это нужно чтобы увеличение радиуса сглаживания не влияло на вычисляемое значение.
Математическое обоснование
Если у нас есть некоторая величина $A(\mathbf{r})$, определенная в пространстве, то мы можем записать ее с помощью дельта-функции: $$ A(\mathbf{r}) = (A * \delta)(\mathbf{r}) = \int A(\mathbf{r’}) \delta(\mathbf{r} - \mathbf{r’}) d\mathbf{r'}
Где $\delta(\mathbf{r})$ – дельта-функция Дирака. Тк с дельта-функцией работать неудобно, мы заменяем ее на kernel function $W(\mathbf{r}, h)$, которая приближает дельта-функцию при $h \to 0$:
И тогда мы можем записать аппроксимацию величины $A(\mathbf{r})$ как свертку с kernel function:
$$ A(\mathbf{r}) \approx (A * W)(\mathbf{r}) = \int A(\mathbf{r’}) W(\mathbf{r} - \mathbf{r’}, h) d\mathbf{v’} $$ где $\mathbf{v’}$ – объемный элемент в точке $\mathbf{r’}$.
Отсюда вытекают требования к kernel function:
- $\lim_{h \to 0} W(\mathbf{r}, h) = \delta(\mathbf{r})$ (приближение дельта-функцию)
- $\int W(\mathbf{r}, h) d\mathbf{r} = 1$ (нормировка)
- $W(\mathbf{r}, h) \ge 0$ (неотрицательность)
- $W(\mathbf{r}, h) = W(-\mathbf{r}, h)$ (симметрия)
- $W(\mathbf{r}, h) = 0$ при $\lVert \mathbf{r} \rVert > kh$ для некоторого $k$ (ограниченность влияния)
Ну в общем она должна выглядеть как-то так:
Например типичная kernel function – кубический сплайн:
$$ W(\mathbf{r},h)=\sigma_{d} \begin{cases} 6\bigl(q^{3}-q^{2}\bigr)+1, & 0 \le q \le \tfrac{1}{2},\\ 2(1-q)^{3}, & \tfrac{1}{2} < q \le 1,\\ 0, & \text{otherwise}, \end{cases} $$
где $h$ – радиус сглаживания, $r$ – расстояние от частицы до точки, в которой мы вычисляем значение, $q=\dfrac{\lVert\mathbf{r}\rVert}{h}$, а $\sigma_{d}$ – нормировочный коэффициент, зависящий от размерности пространства:
$$ \sigma_{d}= \begin{cases} \frac{8}{\pi h^{3}}, & d=3,\\ \frac{40}{7 \pi h^{2}}, & d=2,\\ \frac{4}{3 h}, & d=1. \end{cases} $$
Вычисление значения
Математическое обоснование
Теперь мы можем вычислить значение любой величины $A$ в точке $\mathbf{r}$ с помощью следующей формулы:
$$ (A * W)(\mathbf{r}) = \int A(\mathbf{r’}) W(\mathbf{r} - \mathbf{r’}, h) d\mathbf{v’} $$
где $\mathbf{v’}$ – объемный элемент в точке $\mathbf{r’}$. И для точечной массы $m$ и плотности $\rho$ мы можем записать $d\mathbf{v’} = \frac{m}{\rho} d\mathbf{r’}$.
Если мы моделируем жидкость с помощью частиц, то мы можем аппроксимировать интеграл суммой по всем частицам:
$$ A(\mathbf{r}) = \sum_{j} A_{j} \frac{m_{j}}{\rho_{j}} W(\mathbf{r} - \mathbf{r_{j}}, h) $$
Получается формула для вычисления плотности в точке $\mathbf{r}$:
$$ \rho(\mathbf{r}) = \sum_{j} m_{j} W(\mathbf{r} - \mathbf{r_{j}}, h) $$
В итоге очень простая идея: суммируем вклад всех частиц, взвешенный по их массе и kernel function.
Пора визуализировать:
Здесь я для каждой ячейки сетки вычисляю плотность; чем выше плотность, тем ярче цвет. Kernel function – кубический сплайн. Зная плотность, мы можем вычислить давление. А зная давление, можно вычислить силы. В итоге можно симулировать движение жидкости.
Получаем псевдокод алгоритма SPH:
function SPH_Simulation(particles) {
for (let i = 0; i < particles.length; i++) {
// Вычисляем плотность в частице i
rho_i = calculateDensity(i, particles)
}
for (let i = 0; i < particles.length; i++) {
// Вычисляем давление в частице i
P_i = calculatePressure(rho_i)
// Вычисляем силы на частицу i
F_i = calculateForces(i, particles, P_i, rho_i)
}
for (let i = 0; i < particles.length; i++) {
// Обновляем положение и скорость частицы i
updateParticle(i, F_i)
}
}В основном дальше каждая статья по SPH отличается определением функций для вычисления давления и сил и их использование