Фізична симуляція · Чисельні методи · Комп’ютерна графіка
📅 Березень 2026 ⏱ ≈ 14 хв читання 🎯 Просунутий · Останнє оновлення: 23 червня 2026 р.

Метод матеріальної точки (MPM) — як симулювати сніг, пісок і м’які тіла

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

Автор: Команда MySimulator · Редакційна перевірка: Редакція MySimulator

1. Чому MPM? Лагранж проти Ейлера

Симуляції механіки суцільного середовища поділяються на дві родини:

MPM перебуває посередині: частинки несуть історію матеріалу (градієнт деформації, напруження, швидкість, густину), тоді як свіжа фонова сітка перестворюється на кожному часовому кроці для опрацювання контакту, проєкції тиску та граничних умов. Сітка відкидається після кожного кроку — заплутування неможливе.

Витоки: MPM був запроваджений Сульскі та ін. (1994) як розширення методу частинок у комірках (PIC) з фізики плазми. Стаття Disney Research для SIGGRAPH 2013 року («A Material Point Method for Snow Simulation») принесла його у кіновиробництво, а відкрита мова Taichi (2019) зробила MPM на GPU доступним для дослідників у всьому світі.

2. Огляд алгоритму: чотири кроки на кадр

Кожен крок симуляції виконує чотири ключові операції:

  1. P2G (частинка → сітка): розмазування маси та імпульсу частинки на сусідні вузли сітки за допомогою інтерполяційних ваг
  2. Сили на сітці: обчислення та застосування сил (пружне напруження, гравітація, імпульси зіткнень) до швидкостей вузлів сітки
  3. Граничні умови: накладання стін, липких/слизьких поверхонь, об’єктів на сітці
  4. G2P (сітка → частинка): оновлення швидкостей та положень частинок із сітки; оновлення градієнта деформації
Для частинки p, для вузла сітки i: w_{ip} = вага вузла i, що впливає на частинку p (квадратичний B-сплайн або кубічний, на основі відстані в сітці) P2G: m_i += w_{ip} · m_p mv_i += w_{ip} · m_p · v_p (+ кутовий доданок APIC) f_i += −w_{ip} · vol_p · σ_p · ∇w_{ip} (пружна сила) Сітка: v_i^new = (mv_i + dt · f_i) / m_i G2P: v_p += Σ_i w_{ip} · v_i^new (реконструкція PIC) x_p += dt · v_p F_p = (I + dt · Σ_i v_i ⊗ ∇w_{ip}) · F_p (градієнт деформації)

3. P2G — перенесення частинка-сітка

Кожна частинка впливає на найближчі 3×3 (2D) або 3×3×3 (3D) вузли сітки. Інтерполяційне ядро — це квадратичний B-сплайн — обраний за компактний носій (2 комірки), C¹-неперервність та властивість розкладу одиниці:

1-D вага квадратичного B-сплайна для відстані d = |x_p − x_i| / h: N(d) = 3/4 − d² якщо d ≤ 1/2 (3/2 − d)² / 2 якщо 1/2 < d ≤ 3/2 0 інакше 2-D: w_{ip} = N(|(x_p − x_i)/h|) · N(|(y_p − y_i)/h|) Розклад одиниці: Σ_i w_{ip} = 1 (сумарна вага дорівнює 1) Повна інтерполяція: x_p = Σ_i w_{ip} · x_i (відновлює положення частинки)

Внесок пружного напруження у силу на сітці походить від першого напруження Піоли–Кірхгофа P = ∂ψ/∂F, де ψ — це густина енергії деформації, а F — градієнт деформації:

Сила на сітці від частинки p: f_i += −vol_p^0 · P_p · F_p^T · ∇w_{ip} де vol_p^0 = початковий об’єм частинки (зберігається як маса / ρ₀) P_p = тензор першого напруження Піоли–Кірхгофа F_p = градієнт деформації (матриця 3×3, спочатку I)

4. Оновлення сітки — сили та граничні умови

Щойно маса та імпульс растеризовані на сітку, сили інтегруються явно:

v_i^{n+1} = v_i^n + dt · (f_i / m_i + g) g = гравітаційне прискорення [0, −9.8] м/с² Граничні умови (липка стіна при x = 0): if (x_i ≤ 0 && v_i.x < 0): v_i.x = 0 // без проникнення За потреби також: v_i.y = 0 // тертя (липке) Або просто: (залишити v_i.y без змін) // без тертя

Накладання граничних умов на сітку тривіальне порівняно з лагранжевим МСЕ, де контакт потребує складних формулювань функцій зазору. Це одна з ключових практичних переваг MPM.

5. G2P — перенесення сітка-частинка

Після оновлення сітки інформація повертається до частинок. Критична операція — оновлення градієнта деформації F, який кодує те, як матеріал розтягнувся та зазнав зсуву від початку симуляції:

Градієнт швидкості в частинці p: (∇v)_p = Σ_i v_i^new ⊗ ∇w_{ip} (зовнішній добуток) Оновлення градієнта деформації: F_p^{n+1} = (I + dt · (∇v)_p) · F_p^n Це мультиплікативне оновлення градієнта деформації. Для гіперпружних матеріалів J = det(F) відстежує зміну об’єму: J = 1: нестислива течія J < 1: стиск J > 1: розширення (розтяг)

Положення частинок оновлюються останніми, після того як пластичне зворотне відображення обмежило F до поверхні текучості (для моделей снігу/піску):

x_p^{n+1} = x_p^n + dt · v_p^{n+1}

6. Конститутивні моделі

Градієнт деформації F керує обчисленням напружень. Різні матеріали використовують різні густини енергії деформації:

Неогуківське пружне тверде тіло

ψ(F) = μ/2 · (tr(FᵀF) − 3) − μ·ln(J) + λ/2·ln(J)² Перше напруження PK: P = μ(F − F⁻ᵀ) + λ·ln(J)·F⁻ᵀ Параметри Ламе з модуля Юнга E та коефіцієнта Пуассона ν: μ = E / (2(1+ν)) λ = Eν / ((1+ν)(1−2ν))

Сніг (Disney 2013)

Сніг руйнується при розтягу та ущільнюється при стиску. Це моделюється через мультиплікативну пластичність:

F = F_E · F_P (пружно-пластичний розклад) SVD F_E: F_E = U · Σ · Vᵀ Обмеження сингулярних значень (межі розтягу): Σ_i = clamp(Σ_i, 1 − θ_c, 1 + θ_s) θ_c = критичний стиск (напр., 0.025) θ_s = критичний розтяг (напр., 0.0075) Зміцнення: μ, λ масштабуються на exp(ξ(1 − J_p)) ξ = коефіцієнт зміцнення ≈ 10 (свіжий сніг) до 0 (мокра сльота)

Пісок Друкера–Прагера

Пісок має критерій текучості за кутом тертя — він витримує стиск, але не розтяг, а міцність на зсув залежить від тиску:

Поверхня текучості (Друкер–Прагер): f(σ) = ||dev(σ)|| + sin(φ)·tr(σ) ≤ 0 φ = кут тертя (пісок: ~30°) Дилатансія: зміна об’єму при пластичному зсуві Зворотне відображення: проєкція напруження на найближчу точку поверхні текучості Якщо в розтязі (tr(σ) > 0): руйнування → F = I, нульове напруження

Слабостислива рідина (тиск у стилі SPH)

Тиск за рівнянням стану: p = k₀ · ((J⁻γ) − 1) k₀ = модуль об’ємного стиску, γ = 7 (вода) Напруження Коші: σ = −pI (лише ізотропний тиск) Без напруження зсуву → наближення нев’язкої рідини

7. Варіанти APIC та FLIP

Базова схема PIC (реконструкція швидкості суто із сітки) є чисельно дифузійною — кутовий момент втрачається. Існують два покращення:

APIC P2G (частинка → сітка): mv_i += w_{ip} · m_p · (v_p + C_p · (x_i − x_p)) APIC G2P (сітка → частинка): v_p = Σ_i w_{ip} · v_i C_p = (4/h²) · Σ_i w_{ip} · v_i ⊗ (x_i − x_p) (градієнт швидкості)
MLS-MPM (Moving Least Squares, Hu та ін. 2018): Додатково об’єднує APIC з обчисленням сил на основі квадратури. Доданок сили на сітці природно випливає з перенесення імпульсу APIC, удвічі зменшуючи складність коду. Реалізація MLS-MPM на Taichi у 88 рядках може симулювати пружне тіло, сніг, рідину та пісок лише заміною конститутивної моделі.

8. Основний цикл MLS-MPM на JavaScript

Нижче показано основний 2-D цикл MLS-MPM — спрощений для ясності. Повна робоча реалізація з рендерингом на canvas потребує ~200 рядків: 

// Ядро MLS-MPM: сітка 64×64, ваги квадратичного B-сплайна
const n = 64; // роздільність сітки
const dt = 1e-4;
const dx = 1 / n;
const inv_dx = n;

// Стан частинок
let pos, vel, C, J; // положення, швидкість, афінна матриця, відношення об’ємів

function mpmStep() {
  // 1. Скидання сітки
  const gm  = new Float32Array(n * n);      // маса
  const gmv = new Float32Array(n * n * 2); // імпульс xy

  // 2. P2G
  for (let p = 0; p < pos.length; p++) {
    const [px, py] = pos[p];
    const base = [Math.floor(px * inv_dx - .5), Math.floor(py * inv_dx - .5)];

    // Ваги квадратичного B-сплайна  (3 вузли на вісь)
    const fx = [px * inv_dx - base[0], px * inv_dx - base[0] - 1, px * inv_dx - base[0] - 2];
    const fy = [py * inv_dx - base[1], py * inv_dx - base[1] - 1, py * inv_dx - base[1] - 2];
    const wx = [.5 * (1.5-fx[0])**2, .75-fx[1]**2, .5*(.5+fx[2])**2];
    const wy = [.5 * (1.5-fy[0])**2, .75-fy[1]**2, .5*(.5+fy[2])**2];

    // Доданок напруження MLS-MPM (неогуківська рідина: P = k₀(1 - 1/J) I)
    const stress = -dt * 4 * inv_dx * inv_dx * 400 * (1 - 1 / J[p]);

    for (let di = 0; di < 3; di++)
      for (let dj = 0; dj < 3; dj++) {
        const w = wx[di] * wy[dj];
        const ix = base[0] + di, iy = base[1] + dj;
        if (ix < 0 || ix >= n || iy < 0 || iy >= n) continue;
        const idx = ix * n + iy;
        const dpos = [(di - fx[di]) * dx, (dj - fy[dj]) * dx];
        gm[idx]        += w;
        gmv[idx*2]   += w * (vel[p][0] + C[p][0]*dpos[0] + stress * dpos[0]);
        gmv[idx*2+1] += w * (vel[p][1] + C[p][3]*dpos[1] + stress * dpos[1]);
      }
  }
  // 3. Оновлення сітки: нормалізація за масою, додавання гравітації, накладання межі
  // 4. G2P: реконструкція швидкості та оновлення F, x  (цикл по частинках)
}
Нотатка щодо продуктивності: наївна реалізація на JS опрацьовує ~5 000 частинок при 60 fps. Обчислювальні шейдери WebGL (Transform Feedback або обчислення WebGPU) можуть підняти це до 100 000–1 000 000 частинок, оскільки всі цикли P2G та G2P легко паралелізуються для кожної частинки.
🏖️ Спробувати симуляцію піску →