⚙️ Симуляція фізики · Механіка суцільного середовища
📅 Квітень 2026⏱ 15 хв🔴 Просунутий · Останнє оновлення: 3 липня 2026 р.

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

Сніжинки, що хрумтять під ногами, піщані дюни, що реалістично обвалюються, і скло, що розбивається з фізично точними візерунками тріщин — усі три стали можливими в кіно та іграх завдяки методу матеріальної точки. MPM — це гібридна техніка симуляції суцільного середовища, що представляє матеріал як лагранжеві частинки, які несуть фізичний стан, водночас розв'язуючи рівняння руху на ейлеровій фоновій сітці. Результат — метод, який справляється з великими деформаціями, фрагментацією та фазовими переходами, що розірвали б FEM-сітку й затопили б частинкову симуляцію SPH.

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

1. Навіщо MPM? Обмеження FEM та SPH

Два традиційні підходи до симуляції суцільного середовища мають фундаментальні обмеження, які MPM покликаний подолати:

Метод скінченних елементів (FEM)

FEM дискретизує матеріал у сітку елементів і розв'язує диференціальне рівняння в частинних похідних (PDE) на цій сітці. Фізична точність відмінна. Проблема: сітка прив'язана до топології матеріалу. Коли матеріал руйнується, рветься або зазнає екстремальної деформації (розтяг у 10× чи більше), елементи сітки інвертуються (від'ємний якобіан), і симуляція катастрофічно зривається. Перебудова сітки за кожної події руйнування дорога й вимагає складної топологічної хірургії.

Згладжена гідродинаміка частинок (SPH)

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

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

MPM спочатку розробили в 1994 році Сульскі, Чен і Шраєр як узагальнення методу Fluid-Implicit-Particle (FLIP) на механіку твердих тіл. Він лишався дослідницьким методом, доки Disney Research не розробила пружно-пластичну модель снігу для фільму «Крижане серце» (2013) й детально опублікувала алгоритм — спричинивши величезний інтерес у спільноті комп'ютерної графіки.

2. Алгоритм MPM: частинка-сітка-частинка

Кожен часовий крок MPM складається з трьох фаз — перенесення стану з частинок на сітку, розв'язання фізики на сітці та перенесення оновленого стану назад на частинки:

// Часовий крок MPM: Частинка → Сітка → Частинка // === ФАЗА 1: Частинка на сітку (P2G) === // Скинути сітку for all grid nodes i: m_i = 0; v_i = 0 // Розкидати масу й імпульс частинки на сітку через ваги інтерполяції for each particle p: for each grid node i in p's support (зазвичай 3×3×3 вузли): w = weight(x_p, x_i) // вага інтерполяції m_i += w * m_p v_i += w * m_p * v_p // зважений імпульс f_i -= w * V_p * σ_p ∇w // дивергенція напружень → сила на сітці // Обчислити швидкість на сітці з накопиченого імпульсу for all grid nodes i with m_i > 0: v_i = v_i / m_i // імпульс → швидкість // === ФАЗА 2: Розв'язок на сітці === for all grid nodes i: v_i_new = v_i + Δt * (f_i / m_i + g) // розв'язати Ньютона (g = гравітація) застосувати граничні умови до v_i_new // === ФАЗА 3: Сітка на частинку (G2P) === for each particle p: v_p_new = 0; ∇v_p = 0 for each grid node i in p's support: w = weight(x_p, x_i) v_p_new += w * v_i_new // зібрати швидкість ∇v_p += v_i_new ⊗ ∇w // градієнт швидкості (для деформації) // Оновити градієнт деформації F F_p = (I + Δt * ∇v_p) · F_p // F відстежує сумарну деформацію // Застосувати оновлення пластичності / конститутивне до F_p // Перенести позицію частинки x_p += Δt * v_p_new

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

3. Інтерполяція: зв'язування частинок із сіткою

Ваги інтерполяції w(x_p, x_i) — ядро алгоритму MPM: вони визначають, наскільки кожна частинка робить внесок у кожен вузол сітки й навпаки. Стандартний MPM використовує ядра B-сплайнів, визначені на кроці комірки фонової сітки h:

// 1D вага квадратичного B-сплайна (найпоширеніша для MPM) // Для відстані |d| = |x_p - x_i| / h (в одиницях сітки): w(d) = 0.75 − d² для |d| < 0.5 w(d) = 0.5 * (1.5 − |d|)² для 0.5 ≤ |d| < 1.5 w(d) = 0 для |d| ≥ 1.5 // 3D: w(x_p, x_i) = w(dₓ) * w(d_y) * w(d_z) (роздільна) // Радіус підтримки: ±1.5 комірки → 3×3 вузли у 2D, 3×3×3 у 3D (27 вузлів) // Градієнт ваги (потрібен для обчислення сили): ∇w = (dw/dₓ * w_y * w_z, w_x * dw/d_y * w_z, w_x * w_y * dw/d_z) / h

Кубічні B-сплайни вищого порядку дають гладкіші результати (використовуються в деяких кіно- виробництвах), але розширюють підтримку до 4×4×4 вузлів (64 вузли на частинку у 3D), приблизно потроюючи вартість P2G/G2P. Квадратичний B-сплайн — стандартний вибір для більшості симуляцій.

Компроміс роздільності сітки: дрібніші комірки сітки дають вищу просторову точність, але вимагають пропорційно більше частинок (ціль: 4–8 частинок на комірку сітки), щоб уникнути дірок у представленні. Сітка 128³ з 8 млн частинок — типова MPM-симуляція середньої роздільності. Сніг кіноякості у «Крижаному серці» використовував кількість частинок у сотнях мільйонів.

4. Конститутивні моделі: сніг, пісок, пружне тіло

Конститутивна модель визначає, який матеріал представляє симуляція, задаючи, як напруження σ пов'язане з деформацією F. Кожна частинка несе власний градієнт деформації F (матрицю 3×3, що представляє локальний розтяг і обертання відносно вихідної конфігурації), а конститутивна модель відображає F у тензор напружень Коші σ.

Нео-гуківська пружність (гума, м'які тканини)

// Нео-гуківська гіперпружна модель // F = градієнт деформації (оновлюється щокроку як F ← (I + ΔtL)·F) J = det(F) // коефіцієнт зміни об'єму // Напруження Піоли-Кірхгофа (перше PK): P = μ(F − F⁻ᵀ) + λ log(J) F⁻ᵀ Напруження Коші: σ = (1/J) P Fᵀ Параметри матеріалу: μ = модуль зсуву (жорсткість до зміни форми) λ = модуль об'ємного стиснення (жорсткість до зміни об'єму) Обидва пов'язані з модулем Юнга E та коефіцієнтом Пуассона ν: μ = E / (2(1+ν)), λ = Eν / ((1+ν)(1−2ν))

Пружно-пластична модель снігу (Stomakhin та ін. 2013, Disney)

// Сніг: пружна деформація + пластичний потік через сингулярний розклад // Полярний розклад: F = R·S (обертання × розтяг) // SVD: F = U·Σ·Vᵀ де Σ = diag(σ₁, σ₂, σ₃) — головні розтяги // Пробний пружний крок: вважати всю деформацію пружною F_E_trial = (I + ΔtL)·F_E_prev // Пластичність: обмежити сингулярні значення до [1-θ_c, 1+θ_s] // θ_c = критичне стиснення (сніг ущільнюється: ~2.5%) // θ_s = критичний розтяг (сніг розділяється: ~7.5%) SVD(F_E_trial) → U, Σ, Vᵀ Σ_clamped = clamp(Σ, 1−θ_c, 1+θ_s) F_E = U·Σ_clamped·Vᵀ // пружна частина (обмежена) F_P = Vᵀ·Σ⁻¹·U^ · F_E_trial // пластична частина поглинає залишок J_E = det(F_E) // Модифіковане нео-гуківське напруження з пластичним зміцненням: μ_eff = μ · exp(ξ(1 − J_P)) // ξ = коефіцієнт зміцнення P = 2μ_eff(F_E − F_E⁻ᵀ) + λ(J_E−1)·J_E·F_E⁻ᵀ

Модель піску Друкера-Прагера

// Пісок: гранульований, без зчеплення, текучість залежить від тиску // Умова текучості (Друкер-Прагер): f(σ) = ‖dev(σ)‖ + sin(φ) · tr(σ) ≤ 0 dev(σ) = σ − (1/3)tr(σ)I (девіаторне напруження) φ = кут тертя піску (~30° для сухого піску) // Якщо f > 0 (текучість настала): спроєктувати напруження назад на поверхню текучості // — див. Klar та ін. 2016 щодо алгоритму зворотного відображення // При стисненні: дозволити ущільнення (частинки можуть перекриватися на сітці) // При розтязі: частинки розділяються вільно (без розтягового зчеплення)

5. Руйнування та зміна топології

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

// Руйнування через механіку пошкоджень (найпростіший підхід) // Кожна частинка несе змінну пошкодження d ∈ [0, 1] // Зростання пошкодження, коли розтягова деформація перевищує поріг: if (ε_max > ε_crit): d += Δd // накопичити пошкодження if d > 1.0: d = 1.0 // повністю зруйновано // Застосувати пошкодження до напруження: σ_effective = (1 − d) · σ_elastic // Повністю пошкоджені частинки (d→1): більше не передають розтягове напруження // Вони все ще роблять внесок у масу й можуть стикатися як фрагменти // Шлях тріщини виникає автоматично з топології поля пошкоджень

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

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

6. APIC та MLS-MPM: сучасні варіанти

APIC: Affinely Particle-In-Cell (Jiang та ін. 2015)

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

// APIC: кожна частинка зберігає (v_p, C_p) — швидкість + афінна матриця // Перенесення P2G (з афінним коригувальним членом): momentum_i += w_ip * m_p * (v_p + C_p · (x_i − x_p)) // ↑ афінна корекція: точніший сплеск імпульсу // Перенесення G2P: v_p = Σᵢ w_ip * v_i_grid C_p = (4/h²) Σᵢ w_ip * v_i_grid ⊗ (x_i − x_p) // C_p несе градієнт швидкості → не потрібне окреме обчислення ∇v // Переваги: точно зберігає момент імпульсу (виправляє витік обертання у FLIP) // нульова числова дисипація без шуму в стилі FLIP

MLS-MPM: MPM на основі рухомих найменших квадратів (Hu та ін. 2018)

MLS-MPM ще більше спрощує APIC, використовуючи апарат рухомих найменших квадратів для уніфікації перенесень P2G та G2P. Результат — формулювання, де весь крок MPM можна записати у ≈60 рядках коду без окремого обчислення дивергенції напружень — напруження Коші входить природно через оновлення афінної матриці. MLS-MPM — основа популярного фреймворку симуляції taichi і того, що реалізує більшість сучасних навчальних матеріалів з MPM.

// MLS-MPM P2G (повний): напруження входить як частина афінного імпульсу // D_p = (1/4) h² I (для квадратичних B-сплайнів) momentum_i += w_ip * (m_p * v_p + (m_p * C_p − Δt * V_p * σ_p D_p⁻¹)(x_i−x_p)) // ↑ член напружень вкладено в афінний імпульс // Жодного окремого циклу сил! Дивергенція напружень вбудована в афінний сплеск. // Увесь 2D MPM-рушій: ~50 рядків коду

7. Начерк реалізації на JavaScript

Мінімальне 2D пружне тіло MLS-MPM можна реалізувати менш ніж у 100 рядках. Ключові структури даних — плоскі типізовані масиви для продуктивності:

// 2D MLS-MPM основні структури даних const N_GRID = 64; // роздільність сітки (64×64) const DX = 1.0 / N_GRID; const INV_DX = N_GRID; const N_P = 8000; // кількість частинок // Стан частинки (розкладка SoA для ефективності кешу) const px = new Float32Array(N_P); // позиція x const py = new Float32Array(N_P); // позиція y const vx = new Float32Array(N_P); // швидкість x const vy = new Float32Array(N_P); // швидкість y const Cxx = new Float32Array(N_P); // афінна матриця C (4 компоненти) const Cxy = new Float32Array(N_P); const Cyx = new Float32Array(N_P); const Cyy = new Float32Array(N_P); const Jdet = new Float32Array(N_P).fill(1); // коефіцієнт об'єму (J = det F) // Сітка (скидається щокроку) const gm = new Float32Array((N_GRID+1)**2); // маса сітки const gvx = new Float32Array((N_GRID+1)**2); // імпульс/швидкість сітки x const gvy = new Float32Array((N_GRID+1)**2); // імпульс/швидкість сітки y function step(dt) { gm.fill(0); gvx.fill(0); gvy.fill(0); // скидання перед P2G for (let p = 0; p < N_P; p++) { // P2G const bx = Math.floor(px[p] * INV_DX - 0.5); const by = Math.floor(py[p] * INV_DX - 0.5); // ... ваги квадратичного B-сплайна ... // ... розкидати масу + імпульс + напруження ... } for (let i = 0; i < (N_GRID+1)**2; i++) { // оновлення сітки if (gm[i] > 0) { gvx[i] = gvx[i]/gm[i]; gvy[i] = gvy[i]/gm[i] + dt * -9.8; // гравітація } // граничні умови: нульова швидкість на краях сітки } for (let p = 0; p < N_P; p++) { // G2P // зібрати швидкість і матрицю C із сітки // оновити J: J *= 1 + dt*(dv_x/dx + dv_y/dy) // перенести: px[p] += dt*vx[p]; py[p] += dt*vy[p] } }
Продуктивність: 2D MLS-MPM симуляція з 8 000 частинок працює на ~120 кадрах/с у JavaScript. 3D симуляція зі 100 000 частинок вимагає реалізації обчислювального шейдера WebGL чи WebGPU, щоб працювати в реальному часі — цикли P2G/G2P ганебно паралельні за частинками, а оновлення сітки паралельне за вузлами сітки.

8. Застосування: Disney, ігри, інженерія

Кіно-VFX

Ігри (MPM у реальному часі)

Нещодавні досягнення в GPU-прискореному MPM дозволяють застосування в реальному часі в іграх. Ключова оптимізація — запуск сітки на грубішій роздільності (наприклад, 32³ комірки) з меншою кількістю частинок (10 000–50 000), із використанням обчислювальних шейдерів WebGPU чи Vulkan для частинково-сіткових перенесень.

Наукові та інженерні застосування

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