Mass spring algorithm:

This example shows a mass spring system that consists of particles linked by damped springs:

1. compute spring forces
2. time integration to get new particle positions and velocities

Note that the simulation is only conditionally stable since a conditionally stable explicit time integration method is used.

1. Compute spring forces

In general a spring force can be obtained for any holonomic constraint. In this example we use distance constraints $$C_i(\mathbf{x}_{i_1}, \mathbf{x}_{i_2}) = \| \mathbf x_{i_1} -\mathbf x_{i_2} \|-d,$$ where $d$ is the rest length between particles $\mathbf{x}_{i_1}$ and $\mathbf{x}_{i_2}$.

Potential energy

For a scalar constraint we can define a potential energy as: $$E(\mathbf x) = \frac k 2 C(\mathbf x)^2,$$ where $k$ is the stiffness of the spring.

General spring force

The spring force for a particle $j$ is then determined by the negative gradient of the potential energy function: $$\mathbf f_j = - \frac{\partial E(\mathbf x)}{\partial \mathbf x_j} = -k \frac{\partial C(\mathbf x)}{\partial \mathbf x_j} C(\mathbf x).$$

General damping force

The corresponding damping force is obtained by using the time derivative of the constraint function: $$\mathbf f^D_j = -\mu \frac{\partial C(\mathbf x)}{\partial \mathbf x_j} \dot{C}(\mathbf x).$$

Constraint gradients:

To compute the spring and damping forces, the constraint gradients are required which are computed as:

$$\begin{align*} \frac{\partial C_i}{\partial \mathbf x_{i_1}} &= \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|} \\ \frac{\partial C_i}{\partial \mathbf x_{i_2}} &= - \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|} \end{align*}$$

Spring force

So finally we get the following spring forces for a distance constraint $C_i(\mathbf{x}_{i_1}, \mathbf{x}_{i_2})$: $$\begin{align*} \mathbf f_{i_1} = -k (\| \mathbf x_{i_1} -\mathbf x_{i_2} \|-d) \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|} \\ \mathbf f_{i_2} = +k (\| \mathbf x_{i_1} -\mathbf x_{i_2} \|-d) \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|}. \end{align*}$$

Damping force

The damping forces for a distance constraint $C_i(\mathbf{x}_{i_1}, \mathbf{x}_{i_2})$ are determined as: $$\begin{align*} \mathbf f^D_{i_1} = -k \left ((\mathbf v_{i_1}- \mathbf v_{i_2}) \cdot \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|} \right ) \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|} \\ \mathbf f^D_{i_2} = +k \left ((\mathbf v_{i_1}- \mathbf v_{i_2}) \cdot \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|} \right ) \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|}. \end{align*}$$

2. Time integration

Finally, the particles are advected by numerical time integration. In our case we use a symplectic Euler method: $$\begin{align*} \mathbf v(t + \Delta t) &= \mathbf v(t) + \frac{\Delta t}{m} \left (\mathbf f(t) + \mathbf f^D + \mathbf f^{\text{ext}} \right ) \\ \mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t), \end{align*}$$ where $\mathbf f^{\text{ext}}$ are the external forces.