Mass spring algorithm:
This example shows a mass spring system that consists of particles linked by damped springs:
- compute spring forces
- 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.