Physics-Based Simulation

Stiffness and Damping in Constraint Solving

Stiffness

Material stiffness is controlled by a parameter $k\in [0,1]$, which controls how strongly a constraint is enforced. The most straightforward way to incorporate stiffness is to scale the calculated correction by $k$: $\Delta {\mathbf{p}}_i\leftarrow k\cdot \Delta {\mathbf{p}}_i$. However, this makes the effective stiffness dependent on the number of solver iterations, $n_s$. If a constraint is projected $n_s$ times, the remaining error after the projections would be proportional to $(1-k{)}^{n_s}$. A more robust formulation that decouples stiffness from the iteration count is: $$k^{\prime }=1-(1-k{)}^{1/n_s}\text{\hspace{1em}(32.5.1)}$$ By scaling the correction by $k^{\prime }$ instead of $k$, the remaining error after $n_s$ iterations becomes $(1-k^{\prime }{)}^{n_s}=1-k$, which is independent of $n_s$. This allows artists to tune the material stiffness $k$ without worrying about how many solver iterations are being used.

Damping

Although the PBD method is generally stable, simulations can exhibit excessive oscillations or kinetic energy gain, particularly with stiff constraints. The quality of dynamic simulations can generally be improved by the incorporation of an appropriate damping scheme. It can improve the stability by reducing temporal jittering of the point positions of an object. It also allows for larger time steps which increases the performance of a dynamic simulation.

A naive damping of all velocities can undesirably affect the global motion of an object, slowing down its overall translation and rotation. A more sophisticated approach is to damp only the relative motions of particles while preserving the total linear and angular momentum of the system.

Method 32.5.1 (Momentum-Conserving Damping).

1. Compute Center of Mass and Velocity: ${\mathbf{p}}_{cm}=\frac{{\sum }_i{m}_i{\mathbf{p}}_i}{{\sum }_i{m}_i}$, ${\mathbf{v}}_{cm}=\frac{{\sum }_i{m}_i{\mathbf{v}}_i}{{\sum }_i{m}_i}$

2. Compute Angular Momentum and Inertia Tensor: Let ${\mathbf{r}}_i={\mathbf{x}}_i-{\mathbf{p}}_{cm}$.

   $\mathbf{L}={\sum }_i{\mathbf{r}}_i\times (m_i{\mathbf{v}}_i)$

   $\mathbf{I}={\sum }_i{m}_i(({\mathbf{r}}_i\cdot {\mathbf{r}}_i)\mathbf{Id}-{\mathbf{r}}_i\otimes {\mathbf{r}}_i)$

3. Compute Angular Velocity: $\mathbf{\omega }={\mathbf{I}}^{-1}\mathbf{L}$

4. Apply Damping: For each particle $i$:

   a. Calculate the velocity deviation from rigid body motion: $\Delta {\mathbf{v}}_i=({\mathbf{v}}_{cm}+\mathbf{\omega }\times {\mathbf{r}}_i)-{\mathbf{v}}_i$

   b. Apply damping to the deviation: ${\mathbf{v}}_i\leftarrow {\mathbf{v}}_i+k_{damping}\Delta {\mathbf{v}}_i$

This method effectively isolates the non-rigid components of motion and damps only them. In the extreme case where $k_{damping}=1$, all relative motion is eliminated, and the object behaves as a perfect rigid body.

Extended Position Based Dynamics (XPBD)

While the stiffness formulation in Equation (32.5.1) decouples the effective stiffness from the solver iteration count, a fundamental limitation of standard PBD remains: the resulting material stiffness is still dependent on the simulation time step $\Delta t$. An extension known as Extended Position Based Dynamics (XPBD) addresses this issue, enabling true parameter-independent stiffness.

XPBD is derived from a compliant constraint formulation and introduces the concept of compliance as the inverse of stiffness.

Definition 32.5.1 (Compliance). Compliance, $\alpha \ge 0$, is the inverse of a material's stiffness $k$, with $\alpha =1/k$. It describes a material's propensity to deform under load. In XPBD, a compliance of $\alpha =0$ corresponds to an infinitely stiff, or hard, constraint.

The core idea of XPBD is to treat the Lagrange multiplier $\lambda$ not as a temporary value recalculated in each iteration, but as a physical quantity representing the accumulated impulse that is incrementally updated. In each solver iteration, we calculate an impulse increment, $\Delta \lambda$, and add it to the total impulse $\lambda$ for that constraint. The formula for this increment modifies Equation (32.4.3) as follows: $$\Delta \lambda =\frac{-C(\mathbf{p})-\hat{\alpha }\lambda }{{\sum }_j{w}_j|{\nabla }_{{\mathbf{p}}_j}C(\mathbf{p}){|}^2+\hat{\alpha }}\text{\hspace{1em}(32.5.2)}$$ Here, $\lambda$ is the total accumulated Lagrange multiplier for the constraint from previous iterations within the current time step. The term $\hat{\alpha }$ is the time-step scaled compliance, defined as: $$\hat{\alpha }=\frac{\alpha }{(\Delta t{)}^2}$$ This scaling ensures that the compliance parameter $\alpha$ has physically consistent units within the dynamical system. After computing $\Delta \lambda$, the solver updates both the particle positions and the accumulated Lagrange multiplier for that constraint: $$\Delta {\mathbf{p}}_i=-\Delta \lambda w_i{\nabla }_{{\mathbf{p}}_i}C(\mathbf{p})\text{and}\lambda \leftarrow \lambda +\Delta \lambda$$

Remark 32.5.1 (Interpretation of XPBD). The $+\hat{\alpha }$ term in the denominator of Equation (32.5.2) acts to limit the magnitude of the corrective impulse $\Delta \lambda$. As compliance $\alpha$ increases (i.e., the material becomes softer), $\hat{\alpha }$ grows, thus reducing the impulse applied per iteration. In the case of zero compliance ($\alpha =0$), the $\hat{\alpha }$ terms vanish, and the formula for $\Delta \lambda$ reduces to the standard PBD formulation in Equation (32.4.3).

The primary benefit of XPBD is not faster convergence, but convergence to a physically consistent state that correctly reflects the user-defined compliance $\alpha$, independent of the time step or iteration count. If the solver is terminated early, the system still exhibits the desired softness to some degree rather than an artificial, uncontrolled compliance. This makes material behavior more predictable and robust.