SPH: fluid on the GPU

1. Why particles?

Традиційні методи симуляції рідин — наприклад, Euler-grid CFD — розбивають простір на фіксовану сітку комірок і відстежують значення тиску та швидкості в кожній комірці. Це добре для труб і каналів, але погано для вільних поверхонь: хвилі, бризки, краплі.

SPH (Smoothed Particle Hydrodynamics) — безсіткове лагранжеве поле. Замість комірок у нас є рухомі частки, кожна з масою, позицією та швидкістю. Рідина описується тим, де знаходяться частки, а не де вони «мали бути». Це дозволяє природно симулювати:

• Вільні поверхні та бризки
• Злиття та розділення крапель
• Взаємодію з довільними твердими тілами
• Різні рідини одночасно (вода + масло)

2. Kernel function (Smoothing Kernel)

Основна ідея SPH: будь-яку фізичну величину A у точці r можна відновити як зважену суму значень сусідніх часток. Ваги задаються smoothing kernel W(r, h), де h — радіус впливу (smoothing length).

Ядро W має бути нормованим (∫W dr = 1), гладким і рівним нулю за межами радіуса h. Популярні варіанти:

Poly6 (Müller et al., 2003)

Широко використовується для обчислення густини. Дешево обчислювати, бо залежить лише від r².

Spiky (Müller et al., 2003)

Використовується для градієнта тиску. Poly6 має нульовий градієнт при r→0, що призводить до «злиплих» часток. Spiky не має цього недоліку.

Wendland C2 (сучасний стандарт)

Компактне ядро вищого порядку гладкості. Не має від'ємних значень, стабільніше за Spiky при великій кількості часток.

3. Pressure and equation of state

Перш ніж рахувати сили, для кожної частки i обчислюємо локальну густину ρᵢ як суму вкладів усіх сусідів j у радіусі h:

З відомою густиною отримуємо тиск через рівняння стану. Для рідини, близької до нестисливої, зазвичай використовують формулу Тайта (Tait's equation):

Pressure force на частку i від усіх сусідів j визначається через симетризований градієнт:

Симетризована форма (pᵢ + pⱼ) / 2ρⱼ гарантує, що дія дорівнює протидії (третій закон Ньютона), запобігаючи нефізичному дрейфу.

4. Viscosity

Viscosity «розмазує» швидкості між сусідніми частками — вода «тягне» за собою повільніших сусідів. Стандартна формула SPH для в'язкісної сили:

де μ — динамічна в'язкість (0.01 для води, 10–1000 для меду), а W_visc — ядро, лапласіан якого гарантовано невід'ємний:

5. Motion integration

Маючи всі сили, інтегруємо рівняння руху Ньютона. Загальне прискорення частки:

де g — зовнішнє прискорення (гравітація, вітер). Для інтеграції в часі зазвичай використовують Leapfrog або Verlet — вони зберігають кінетичну енергію краще, ніж простий метод Ейлера:

Часовий крок Δt обмежений умовою стабільності CFL: частка не повинна переміщатися більше ніж h за один крок. Типово Δt = 0.001 – 0.005 с.

6. Spatial hashing

Наївна реалізація: для кожної частки перебираємо всіх N сусідів → O(N²). При N = 10 000 це 100 мільйонів перевірок на кадр — занадто повільно.

Uniform grid hashing ділить простір на комірки розміром h × h. Кожна частка хешується в комірку за формулою:

Тепер для пошуку сусідів достатньо перевірити лише 9 комірок (2D) або 27 комірок (3D) довкола поточної — O(N · 27) ≈ O(N).

7. Pseudocode

Один кроку симуляції SPH (спрощено):

function stepSPH(particles, dt):

  // 1. Оновити просторовий хеш
  updateHashGrid(particles, h)

  for each particle i:
    // 2. Знайти сусідів у радіусі h
    neighbors = getNeighbors(i, h)

    // 3. Обчислити густину
    ρ[i] = Σ m[j] * W_poly6(dist(i,j), h)
              for j in neighbors

    // 4. Обчислити тиск (Tait)
    p[i] = k * ((ρ[i] / ρ₀)^7 - 1)

  for each particle i:
    neighbors = getNeighbors(i, h)

    // 5. Pressure force
    f_press = -m[i] * Σ m[j] * (p[i]+p[j]) / (2*ρ[j])
                     * gradW_spiky(r[i]-r[j], h)

    // 6. В'язкісна сила
    f_visc = μ * Σ m[j] * (v[j]-v[i]) / ρ[j]
                * lapW_visc(dist(i,j), h)

    // 7. Прискорення
    a[i] = (f_press + f_visc) / ρ[i] + gravity

  for each particle i:
    // 8. Leapfrog інтеграція
    v[i] += a[i] * dt
    r[i] += v[i] * dt

    // 9. Граничні умови (відбиття)
    handleBoundary(i)

Кроки 2–4 (густина/тиск) і 5–7 (сили) обов'язково розділяють у два цикли: при першому проході по всіх частках ρ і p мають бути вже відомі для всіх сусідів.