FLIP Fluid: частки + сітка

1. Навіщо поєднувати частки та сітку?

У симуляції рідин домінують дві родини методів. Чисто ейлерові сіткові розв'язувачі (Navier–Stokes через скінченні різниці / об'єми) дискретизують швидкість і тиск на фіксованій ґратці комірок. Вони роблять розв'язання тиску для нестисливості тривіальним — це просто рівняння Пуассона на сітці, — але адвекція самого поля швидкостей вносить числову дифузію: рідина втрачає енергію та дрібні деталі на кожному кроці, а вільні поверхні (бризки, краплі) складно відстежувати.

Чисто лагранжеві розв'язувачі часток, як-от SPH, ідеально відстежують масу, імпульс і вільні поверхні (частки просто є там, де рідина), але обчислення тиску з неструктурованої хмари часток дороге й схильне до шумної, «стисливої» поведінки, якщо не налаштувати ретельно.

FLIP, запропонований Бреккбіллом і Раппелом (Brackbill & Ruppel, 1986) для фізики плазми та адаптований для графічних рідин Чжу і Бріджсоном (Zhu & Bridson, 2005), поєднує сильні сторони обох підходів:

• Частки несуть швидкість, позицію та форму вільної поверхні — без числової дифузії, без згладжування бризок
• Фонова сітка розв'язує проекцію тиску точно, за один лінійний розв'язок, гарантуючи нестисливе (бездивергентне) поле швидкостей
• Інформація переноситься частка → сітка → частка на кожному кроці, тож кожне представлення робить лише те, у чому воно найкраще

2. PIC проти FLIP проти суміші

Схема перенесення має дві крайності, і виробничі симулятори використовують суміш.

PIC (Particle-In-Cell)

Після розв'язання тиску на сітці PIC просто копіює нову сіткову швидкість назад на кожну частку, якої вона торкається:

Це безумовно стабільно, але повторна інтерполяція на сітку і назад щокадру діє як фільтр низьких частот — кінетична енергія витікає, і рідина за кілька десятків кадрів виглядає густою, загальмованою й «каламутною».

FLIP (Fluid-Implicit-Particle)

Замість того, щоб взяти нову швидкість напряму, FLIP переносить лише зміну, яку дав сітковий розв'язок, і додає цю дельту до власної старої швидкості частки:

Оскільки крізь дисипативну сіткову інтерполяцію проходить лише поправка тиску, майже вся кінетична енергія частки зберігається. Недолік: чистий FLIP шумніший і може накопичувати хибну високочастотну швидкість («шум часток») за довгі прогони.

Суміш PIC/FLIP

Майже кожна реальна реалізація лінійно змішує обидва методи з параметром α (типово 0.95–0.99, тобто переважно FLIP):

Невелика частка демпфування PIC не дає симуляції втратити стабільність, зберігаючи водночас живий, енергійний вигляд FLIP.

3. Перенесення частка → сітка (P2G)

Сіткові величини живуть на MAC-сітці (Marker-and-Cell): тиск зберігається в центрах комірок, а кожна компонента швидкості — на відповідній грані комірки (u на x-гранях, v на y-гранях, w на z-гранях). Таке зміщення запобігає «шаховим» артефактам тиску і робить оператори градієнта тиску / дивергенції простими центральними різницями.

Щоб перенести швидкості часток на сітку, кожна частка розсіює свій імпульс на сусідні вузли сітки за допомогою гладкого ядра інтерполяції — зазвичай трилінійної «наметової» функції, іноді квадратичного B-сплайна для додаткової гладкості:

Ділення на суму ваг (а не кількість часток) нормалізує нерівну кількість часток у кожній комірці — саме це не дає ділянкам з малою кількістю часток отримувати хибно занижені швидкості.

4. Проекція тиску

Після розсіювання сіткове поле швидкостей зазвичай не є бездивергентним — здається, що рідина всередині комірок створюється або зникає. Нестисливість вимагає ∇·u = 0 всюди. Ми досягаємо цього, віднімаючи градієнт поля тиску p, що точно компенсує дивергенцію (розклад Гельмгольца–Ходжа):

Взявши дивергенцію обох частин і вимагаючи, щоб результат дорівнював нулю, отримуємо рівняння Пуассона для тиску:

Дискретизоване на MAC-сітці, це стає одним лінійним рівнянням на кожну комірку FLUID, що пов'язує її тиск із (до 6) сусідами — великою розрідженою симетричною додатно визначеною системою:

A — дискретний лапласіан (7-точковий шаблон у 3D), b — (від'ємна, масштабована) дивергенція швидкості до проекції. Систему розв'язують ітеративно — метод спряжених градієнтів з передобумовленням (PCG) з передобумовлювачем Modified Incomplete Cholesky (MIC) — галузевий стандарт, що сходиться за кілька десятків ітерацій навіть для сіток із мільйонами комірок.

Коли p відоме, кожна грань швидкості, суміжна з комірку FLUID, коригується:

Отримане сіткове поле бездивергентне з точністю до похибки розв'язувача — рідина не стискається і не розширюється.

5. Перенесення сітка → частка (G2P)

Скориговані сіткові швидкості тепер інтерполюються назад на кожну частку, використовуючи ті самі трилінійні ваги, що й P2G (ця симетрія важлива — використання різних ядер для розсіювання й збирання порушує збереження імпульсу). Як описано в розділі 2, виробничий код обчислює як значення PIC, так і дельту FLIP, і змішує їх з α:

де Δv_grid = v_grid^new − v_grid^old, обидва трилінійно інтерпольовані в позиції самої частки.

6. Адвекція та межові умови

Маючи кінцеві швидкості часток, позиції інтегрують вперед — зазвичай схемою Рунге-Кутти 2-го чи 3-го порядку (RK2/RK3), а не простим Ейлером, оскільки частки поблизу бризок можуть переміщатися на кілька комірок за крок:

Тверді межі (стінки контейнера, колайдери) забезпечуються двічі: один раз на сітці, обнуляючи нормальну компоненту швидкості на гранях SOLID перед розв'язанням тиску (умова Неймана), і один раз на частках, повертаючи будь-яку частку, що «протунелювала» крізь стінку, назад усередину й гасячи її внутрішню компоненту швидкості.

Тонкість, унікальна для FLIP: оскільки частки ніколи не переприсаджуються в рівномірний патерн, комірки можуть з часом залишатися без часток (числове «протікання» крізь тонкі щілини) або переповнюватися (злипання). Більшість реалізацій відстежують легковаговий цільовий підрахунок часток на комірку й нагодою переприсаджують/видаляють частки, щоб залишатися, скажімо, у межах 4–12 часток на рідинну комірку.

7. Далі за FLIP: APIC та MLS-MPM

Звичайний FLIP втрачає кутовий момент під час перенесення — закрутіть частку на місці, і цикл P2G/G2P погасить обертання. APIC (Affine Particle-In-Cell), представлений Цзяном та ін. (Jiang et al., 2015), виправляє це, зберігаючи для кожної частки невелику афінну матрицю швидкості C_p (локальний градієнт швидкості) поряд із самою швидкістю, і використовуючи її для відновлення локально-лінійного поля швидкостей під час обох перенесень:

APIC точно зберігає імпульс і кутовий момент у самому кроці перенесення, тому він значною мірою замінив звичайне змішування FLIP у сучасних виробничих розв'язувачах. MLS-MPM (Hu et al., 2018) узагальнює цей підхід ще далі за допомогою формулювання методу рухомих найменших квадратів і поширює той самий механізм частка/сітка на пружні та зернисті матеріали (сніг, пісок, тканина) — це вже поза межами чистої рідини, але побудовано на точно такому ж каркасі P2G/розв'язок/G2P, описаному вище.

8. Псевдокод

Один крок симуляції FLIP, спрощено:

function stepFLIP(particles, grid, dt, alpha):

  // 1. Перенесення частка → сітка (P2G)
  grid.clear()
  for each particle p:
    splatVelocity(grid, p.pos, p.vel)
    grid.markFluidCell(p.pos)
  grid.normalizeByWeights()
  gridVelOld = grid.velocity.copy()

  // 2. Зовнішні сили (гравітація) на сітці
  grid.applyGravity(dt)
  grid.enforceSolidBoundaries()

  // 3. Проекція тиску: A·p = b, розв'язується PCG
  b = grid.computeDivergence()
  p = conjugateGradient(grid.laplacianA, b, iters=100)
  grid.subtractPressureGradient(p, dt)

  // 4. Перенесення сітка → частка (G2P), суміш PIC/FLIP
  for each particle p:
    vPIC  = interpolate(grid.velocity, p.pos)
    dv    = interpolate(grid.velocity - gridVelOld, p.pos)
    vFLIP = p.vel + dv
    p.vel = alpha * vFLIP + (1 - alpha) * vPIC

  // 5. Адвекція часток (RK2) та обробка зіткнень
  for each particle p:
    mid = p.pos + 0.5 * dt * p.vel
    p.pos += dt * interpolate(grid.velocity, mid)
    resolveCollisions(p)

  // 6. Переприсадження/видалення для стабільної густини часток
  grid.rebalanceParticles(particles, targetPerCell=8)

Зверніть увагу на два повних проходи по сітці: P2G має завершитись для кожної частки перед запуском розв'язання тиску, а розв'язання тиску має завершитись перед будь-якою інтерполяцією G2P — сітка є точкою синхронізації, яка робить усю схему зручною для паралелізації на GPU, оскільки частки ніколи не повинні знати одна про одну напряму.