Симулятор частинок із пружинами — Position Based Dynamics (PBD)

1. Чому PBD, а не силові пружини

Класична пружинно-масова система обчислює силу Гука для кожної пружини, F = -k·(|d| - L₀)·d̂, підсумовує сили для кожної частинки і інтегрує їх явним або напівнеявним методом Ейлера. Проблема — жорсткість: щоб пружина здавалась твердою, потрібно збільшувати k, але велике k у поєднанні з фіксованим кроком часу dt робить систему числово жорсткою — інтегратор "перестрибує" рівновагу, енергія зростає з кожним кроком, і вся мережа перетворюється на коливальний, вибуховий безлад. Можна зменшити dt або перейти на неявну інтеграцію (розв'язання лінійної системи щокадру), але обидва варіанти додають реальну складність для браузерної демонстрації на 60 fps.

Position Based Dynamics (Мюллер та ін., 2007) переформульовує задачу. Замість питання "яка сила утримає цю пружину в спокої?", вона запитує "де мають бути ці дві точки, щоб обмеження виконувалося прямо зараз?" Кожна пружина стає обмеженням відстані, яке напряму зміщує позиції обох частинок до цільової довжини, з вагою, пропорційною оберненій масі. Оскільки метод оперує позиціями, а не силами, тут немає вибуху через жорсткість — натомість жорсткість керується тим, скільки разів ви повторюєте проєкцію за кадр.

2. Частинки у стилі Верле

Частинки PBD не зберігають швидкість як окрему змінну під час фази розв'язання обмежень — вони використовують інтегрування Штермера-Верле: швидкість неявно виражена різницею між поточною та попередньою позицією. Це робить проєкцію обмежень тривіальною (достатньо перемістити pos), при цьому рух залишається коректним після відновлення швидкості з позицій наприкінці кроку.

class Particle {
  constructor(x, y, mass = 1) {
    this.pos     = { x, y };       // поточна позиція (спочатку прогнозована, потім скоригована)
    this.prevPos = { x, y };       // позиція наприкінці минулого кадру
    this.vel     = { x: 0, y: 0 };
    this.invMass = mass > 0 ? 1 / mass : 0; // 0 = закріплена/статична
  }
}

function predictPosition(p, dt, gravity) {
  // зовнішні сили (гравітація) діють на швидкість, потім прогнозуємо нову позицію
  p.vel.x += 0 * dt;         // горизонтальна сила відсутня
  p.vel.y += gravity * dt;
  p.prevPos.x = p.pos.x;
  p.prevPos.y = p.pos.y;
  p.pos.x += p.vel.x * dt;
  p.pos.y += p.vel.y * dt;
}

function reconcileVelocity(p, dt) {
  // після того, як обмеження змістили p.pos, виводимо швидкість, що пояснює цей рух
  p.vel.x = (p.pos.x - p.prevPos.x) / dt;
  p.vel.y = (p.pos.y - p.prevPos.y) / dt;
}

3. Обмеження відстані

Пружина між частинками A і B з довжиною спокою L₀ стає функцією обмеження:

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

function solveDistanceConstraint(a, b, restLength, stiffness = 1) {
  const dx = b.pos.x - a.pos.x;
  const dy = b.pos.y - a.pos.y;
  const dist = Math.hypot(dx, dy) || 1e-6; // уникаємо ділення на нуль
  const C = dist - restLength;

  const wSum = a.invMass + b.invMass;
  if (wSum === 0) return; // обидві закріплені, нічого робити

  const nx = dx / dist, ny = dy / dist;
  const corr = (C / wSum) * stiffness; // жорсткість у [0,1], див. §5

  a.pos.x += a.invMass * corr * nx;
  a.pos.y += a.invMass * corr * ny;
  b.pos.x -= b.invMass * corr * nx;
  b.pos.y -= b.invMass * corr * ny;
}

4. Цикл розв'язувача PBD

Повний алгоритм на кадр має три фази: прогноз, проєкція (повторювана), узгодження.

1. Прогноз: застосувати гравітацію та зовнішні сили до швидкостей, потім спрогнозувати нові позиції (обмеження ще не застосовуються).
2. Проєкція: пройтися по всіх обмеженнях N разів, підштовхуючи позиції до задоволення кожного з них. Оскільки переміщення A для виправлення обмеження 1 може порушити обмеження 2, це релаксація у стилі Гаусса-Зейделя — вона збігається до розв'язку, а не розв'язує точно.
3. Узгодження: вивести швидкість зі зміни позиції та зберегти її для прогнозу наступного кадру й для інтерполяції під час рендерингу.

function simulate(particles, springs, dt) {
  const GRAVITY = 9.8;
  const SOLVER_ITERATIONS = 8;

  // 1. Прогноз
  for (const p of particles) {
    if (p.invMass === 0) continue; // закріплені частинки не рухаються
    predictPosition(p, dt, GRAVITY);
  }

  // 2. Проєкція обмежень, кілька проходів для збіжності
  for (let iter = 0; iter < SOLVER_ITERATIONS; iter++) {
    for (const s of springs) {
      solveDistanceConstraint(s.a, s.b, s.restLength, s.stiffness);
    }
  }

  // 3. Узгодити швидкість зі зміни позиції
  for (const p of particles) {
    if (p.invMass === 0) continue;
    reconcileVelocity(p, dt);
  }
}

5. Жорсткість через кількість ітерацій та SOR

Сприйнята "пружність" у PBD залежить від двох параметрів, а не від силової константи:

• Кількість ітерацій: більше проходів Гаусса-Зейделя за кадр наближають обмеження до повного задоволення (C → 0), що сприймається як жорсткіша, жорсткіша пружина. Менше ітерацій залишають залишкове розтягнення — м'якшу, пружнішу пружину.
• Коефіцієнт жорсткості для кожного обмеження k ∈ [0, 1]: масштабує корекцію на кожному проході коефіцієнтом k замість повного застосування. Один прохід із k=1 і N ітераціями — це не те саме, що один прохід із k=1/N — класична стаття про PBD показує, що жорсткість на одну ітерацію треба скоригувати до 1 − (1 − k)^(1/N), щоб сприйнята жорсткість залишалась незалежною від кількості ітерацій.

Корисна ментальна модель: ітерації міняють CPU на жорсткість. Вільна тканина може обійтись 4 ітераціями; натягнуте трамплінне полотно чи майже жорсткий стрижень, апроксимований обмеженнями, вимагає 15–30.

6. Демпфування швидкості та обмеження розтягу

Два вдосконалення суттєво покращують вигляд і поведінку мережі:

Глобальне демпфування швидкості

Без демпфування енергія, внесена перетягуванням частинки, ніколи не згасає, і мережа тремтить нескінченно. Масштабуйте швидкість коефіцієнтом демпфування щокадру:

for (const p of particles) {
  p.vel.x *= 0.995;
  p.vel.y *= 0.995;
}

Обмеження розтягу

Тканина, симульована лише м'якими пружинами, може розтягнутись у кілька разів більше довжини спокою під навантаженням. Поширене рішення — жорстке обмеження максимального розтягу, застосоване після м'яких пружин: якщо пружина перевищує, скажімо, 10% над довжиною спокою, вона примусово повертається назад (жорсткість = 1) незалежно від звичайної м'якості матеріалу.

function limitStretch(a, b, restLength, maxStretch = 1.1) {
  const dx = b.pos.x - a.pos.x;
  const dy = b.pos.y - a.pos.y;
  const dist = Math.hypot(dx, dy) || 1e-6;
  const maxLen = restLength * maxStretch;
  if (dist <= maxLen) return;
  solveDistanceConstraint(a, b, maxLen, 1); // жорстке обмеження, повна жорсткість
}

7. Побудова пружинної мережі (сітка тканини)

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

function buildGrid(cols, rows, spacing) {
  const particles = [];
  const springs = [];

  for (let y = 0; y < rows; y++)
    for (let x = 0; x < cols; x++)
      particles.push(new Particle(x * spacing, y * spacing));

  const idx = (x, y) => y * cols + x;
  const link = (i, j, stiffness) => {
    const a = particles[i], b = particles[j];
    const rest = Math.hypot(a.pos.x - b.pos.x, a.pos.y - b.pos.y);
    springs.push({ a, b, restLength: rest, stiffness });
  };

  for (let y = 0; y < rows; y++)
    for (let x = 0; x < cols; x++) {
      // структурні
      if (x < cols - 1) link(idx(x,y), idx(x+1,y), 0.9);
      if (y < rows - 1) link(idx(x,y), idx(x,y+1), 0.9);
      // зсувні (діагоналі)
      if (x < cols - 1 && y < rows - 1) {
        link(idx(x,y),   idx(x+1,y+1), 0.6);
        link(idx(x+1,y), idx(x,y+1),   0.6);
      }
      // згинальні (через 2 клітинки, запобігають різким складкам)
      if (x < cols - 2) link(idx(x,y), idx(x+2,y), 0.3);
      if (y < rows - 2) link(idx(x,y), idx(x,y+2), 0.3);
    }

  return { particles, springs };
}

8. Закріплення, перетягування та колізії

Мас-зважена корекція PBD робить взаємодію майже безкоштовною:

• Закріплення: встановіть invMass = 0 для кутових частинок, які мають бути нерухомими у просторі (наприклад, верхній ряд підвішеної тканини). Уся вага корекції автоматично переходить на незакріплених сусідів.
• Перетягування: поки кнопка миші натиснута, щокадру напряму перезаписуйте pos перетягуваної частинки на позицію курсора (еквівалентно invMass = 0 для цього кадру). Відпускання миші відновлює звичайний invMass, і reconcileVelocity підхоплює швидкість кидка з руху останнього кадру.
• Колізія з підлогою / іншими тілами: після проходів обмежень пружин виконайте дешевий прохід обмеження позиції — if (p.pos.y > floorY) p.pos.y = floorY; — розглядаючи колізію просто як ще одне обмеження, розв'язане у тому самому циклі.

function applyGroundCollision(particles, floorY, restitution = 0.3) {
  for (const p of particles) {
    if (p.invMass === 0 || p.pos.y <= floorY) continue;
    p.pos.y = floorY;
    // відображаємо збережену prevPos, щоб узгоджена швидкість відскакувала, а не зупинялась миттєво
    p.prevPos.y = floorY + (floorY - p.prevPos.y) * restitution;
  }
}

9. Типові пастки

• Пружина нульової довжини / ділення на нуль: дві частинки, що опинилися в одній точці, дають dist = 0. Завжди захищайтесь малим епсилоном (|| 1e-6) перед діленням.
• Забута мас-зважена корекція: застосування повної корекції до обох частинок порівну (ігноруючи invMass) призводить до дрейфу закріплених частинок і поведінки важких об'єктів так, ніби вони невагомі.
• Замало ітерацій для натягнутих мереж: мотузка, симульована лише з 1–2 ітераціями, виглядає як гума, а не мотузка — довгі ланцюги обмежень потребують пропорційно більше проходів, щоб корекція поширилась з одного кінця на інший.
• Фіксований крок часу проти змінного часу кадру: передавання змінного dt від requestAnimationFrame напряму в прогнозатор робить жорсткість непослідовною при провалах частоти кадрів. Накопичуйте реальний час і крокуйте симуляцію фіксованими частками dt (наприклад, 1/60 с) для відтворюваної поведінки.
• Збіжність, залежна від порядку: проєкція Гаусса-Зейделя чутлива до порядку обмежень, що може внести тонке спрямоване зміщення у великих сітках. Перемішування порядку обмежень щокадру або перехід на схему Якобі з усередненою корекцією прибирає це зміщення ціною трохи повільнішої збіжності.