Particle Simulator with Springs — Position Based Dynamics (PBD)

1. Why PBD Instead of Force Springs

A classic mass-spring system computes a Hookean force for every spring, F = -k·(|d| - L₀)·d̂, sums forces per particle, and integrates with explicit or semi-implicit Euler. The problem is stiffness: to make a spring feel rigid you raise k, but a large k combined with a fixed timestep dt makes the system numerically stiff — the integrator overshoots, energy grows every step, and the whole net turns into an oscillating, exploding mess. You can shrink dt or switch to implicit integration (solving a linear system every frame), but both add real complexity for a browser demo running at 60 fps.

Position Based Dynamics (Müller et al., 2007) reframes the problem. Instead of asking "what force keeps this spring at rest length?", it asks "where should these two points be so the constraint is satisfied right now?" Each spring becomes a distance constraint that directly displaces the two particle positions toward the target length, weighted by inverse mass. Because it operates on positions, not forces, there is no stiffness-induced blow-up — stiffness is instead controlled by how many times you repeat the projection per frame.

2. Verlet-Style Particles

PBD particles don't store velocity as a first-class variable during the constraint phase — they use Störmer-Verlet integration: velocity is implicit in the difference between the current and previous position. This makes constraint projection trivial (just move pos) while still producing correct motion once you reconcile positions back into a velocity at the end of the step.

class Particle {
  constructor(x, y, mass = 1) {
    this.pos     = { x, y };       // current position (predicted, then corrected)
    this.prevPos = { x, y };       // position at the end of last frame
    this.vel     = { x: 0, y: 0 };
    this.invMass = mass > 0 ? 1 / mass : 0; // 0 = pinned/static
  }
}

function predictPosition(p, dt, gravity) {
  // external forces (gravity) act on velocity, then we predict a new position
  p.vel.x += 0 * dt;         // no horizontal force
  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) {
  // after constraints moved p.pos, derive the velocity that explains the move
  p.vel.x = (p.pos.x - p.prevPos.x) / dt;
  p.vel.y = (p.pos.y - p.prevPos.y) / dt;
}

3. The Distance Constraint

A spring between particles A and B with rest length L₀ becomes a constraint function:

This is exactly the same shape as the impulse formula in a force-based solver, except it corrects position directly instead of accumulating a velocity change. Heavier particles (smaller invMass) move less; lighter particles move more — mass ratios are respected automatically.

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; // avoid divide-by-zero
  const C = dist - restLength;

  const wSum = a.invMass + b.invMass;
  if (wSum === 0) return; // both pinned, nothing to do

  const nx = dx / dist, ny = dy / dist;
  const corr = (C / wSum) * stiffness; // stiffness in [0,1], see §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. The PBD Solver Loop

The full per-frame algorithm has three phases: predict, project (repeatedly), reconcile.

1. Predict: apply gravity and external forces to velocities, then predict new positions (no constraints applied yet).
2. Project: loop over all constraints N times, nudging positions toward satisfying each one. Because moving A to fix constraint 1 can un-satisfy constraint 2, this is a Gauss-Seidel-style relaxation — it converges toward a solution rather than solving exactly.
3. Reconcile: derive velocity from the position delta and store it for the next frame's prediction and for rendering-time interpolation.

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

  // 1. Predict
  for (const p of particles) {
    if (p.invMass === 0) continue; // pinned particles don't move
    predictPosition(p, dt, GRAVITY);
  }

  // 2. Project constraints, several passes for convergence
  for (let iter = 0; iter < SOLVER_ITERATIONS; iter++) {
    for (const s of springs) {
      solveDistanceConstraint(s.a, s.b, s.restLength, s.stiffness);
    }
  }

  // 3. Reconcile velocity from the position change
  for (const p of particles) {
    if (p.invMass === 0) continue;
    reconcileVelocity(p, dt);
  }
}

5. Stiffness via Iteration Count and SOR

Perceived "springiness" in PBD comes from two knobs, not from a force constant:

• Iteration count: more Gauss-Seidel passes per frame push the constraint closer to fully satisfied (C → 0), which reads as a stiffer, more rigid spring. Fewer iterations leave residual stretch — a softer, bouncier spring.
• Per-constraint stiffness factor k ∈ [0, 1]: scale the correction each pass by k instead of applying it fully. A single pass with k=1 and N iterations is not the same as one pass with k=1/N — the classic PBD paper shows the per-iteration stiffness must be corrected to 1 − (1 − k)^(1/N) so that the perceived stiffness stays independent of the iteration count.

A useful mental model: iterations trade CPU for rigidity. A loose cloth might use 4 iterations; a taut trampoline sheet or a near-rigid rod approximated with constraints wants 15–30.

6. Velocity Damping and Stretch Limiting

Two refinements make the network look and feel much better:

Global velocity damping

Without damping, energy injected by dragging a particle never decays and the net jitters forever. Scale velocity by a damping factor each frame:

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

Stretch limiting

Cloth simulated with soft springs alone can stretch to several times its rest length under load. A common fix is a hard maximum-stretch constraint applied after the soft springs: if a spring exceeds, say, 10% over rest length, clamp it back exactly (stiffness = 1) regardless of the material's normal softness.

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); // hard clamp, full stiffness
}

7. Building a Spring Network (Cloth Grid)

A grid of particles connected with three kinds of springs produces a convincing cloth or net: structural springs (horizontal/vertical neighbours) hold the shape, shear springs (diagonal neighbours) resist skewing, and bend springs (two cells away) resist folding.

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++) {
      // structural
      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);
      // shear (diagonals)
      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);
      }
      // bend (2 cells away, stops sharp folding)
      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. Pinning, Dragging, and Collisions

PBD's mass-weighted correction makes interaction almost free:

• Pinning: set invMass = 0 on the corner particles you want fixed in space (e.g. the top row of a hanging cloth). All correction weight automatically shifts to the unpinned neighbours.
• Dragging: while the mouse is down, each frame directly overwrite the dragged particle's pos to the cursor position (equivalent to invMass = 0 for that frame). Releasing the mouse restores its normal invMass and lets reconcileVelocity pick up a throw velocity from the last frame's motion.
• Collision with the ground / other bodies: after the spring-constraint passes, run a cheap position-clamp pass — if (p.pos.y > floorY) p.pos.y = floorY; — treating collision as just another constraint solved in the same loop.

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;
    // reflect stored prevPos so reconciled velocity bounces instead of stopping dead
    p.prevPos.y = floorY + (floorY - p.prevPos.y) * restitution;
  }
}

9. Common Pitfalls

• Zero-length spring / divide-by-zero: two particles landing on the exact same point make dist = 0. Always guard with a small epsilon (|| 1e-6) before dividing.
• Forgetting mass weighting: applying the full correction to both particles equally (ignoring invMass) makes pinned particles drift and heavy objects behave like they're weightless.
• Too few iterations for taut networks: a rope simulated with only 1–2 iterations looks like rubber, not rope — long constraint chains need proportionally more passes to propagate correction from one end to the other.
• Fixed timestep vs. variable frame time: feeding a variable dt from requestAnimationFrame straight into the predictor makes stiffness feel inconsistent across frame-rate dips. Accumulate real time and step the simulation in fixed dt chunks (e.g. 1/60 s) for reproducible behaviour.
• Order-dependent convergence: Gauss-Seidel projection is sensitive to constraint iteration order, which can introduce a subtle directional bias in large grids. Shuffling constraint order per frame, or switching to a Jacobi scheme with averaged corrections, removes the bias at the cost of slightly slower convergence.