18. Mechanical Systems

We consider a system of point masses \(m_i\) at time-dependent positions \(x_i(t)\).

There is potential energy \(U(x)\) stored in that system. Examples are gravitational energy

\[ \sum m_i g x_{i,z}, \]

or energy stored in springs \(k_{ij}\) connecting masses \(m_i\) and \(m_j\):

\[ \sum_{k_{ij}} \frac{1}{2} k_{ij} ( |x_i - x_j | - l_{ij} )^2, \]

where \(l_{ij}\) is the length of the spring without load, and \(|x_i - x_j|\) is the Euclidean distance.

Taking the negative derivative

\[ F_i := -\frac{\partial U}{\partial x_i}(x) \]

gives the force onto the mass.

In case of gravity it is \(-m_i g e_z\), in case of springs it is

\[ -\sum_j k_{ij} ( |x_i - x_j | - l_{ij} ) \frac{x_i - x_j}{|x_i - x_j|} \]

i.e. spring constant multiplid by elongation, in the direction of the vector \(x_j-x_i\).

By Newton’s principle, the force leads to acceleration:

\[ m_i \ddot x_i = -\frac{\partial U}{\partial x_i} \]

Now we denote the derivative w.r.t. time \(t\) via \(\dot {x}\).

18.1. The Newmark method

We start with reducing the second order equation to a first order system. Therefore, we introduce velocity \(v_i = \dot x_i\):

\[\begin{split} \begin{array}{rcl} \dot x_i & = & v_i \\ m_i \dot v_i & = & -F_i(x) \end{array} \end{split}\]

Applying the Crank-Nicolson time-stepping scheme we obtain

\[\begin{split} \begin{array}{rcl} x_{n+1} - x_n & = & \frac{h}{2} (v_n + v_{n+1} ) \\ M (v_{n+1} - v_n) & = & \frac{h}{2} (F(x_n) + F(x_{n+1})) \end{array} \end{split}\]

By renaming \(a_n := M^{-1} F(x_n)\), we can substitute

\[\begin{split} \begin{array}{rcl} v_{n+1} & = & v_n + \frac{h}{2} (a_n + a_{n+1}) \\ x_{n+1} & = & x_n + \frac{h}{2} (v_n + v_{n+1}) = x_n + h v_n + \frac{h^2}{4} (a_n + a_{n+1} ) \end{array} \end{split}\]

and obtain the equation for the new \(a_{n+1}\):

\[ M a_{n+1} = F \big( x_n + h v_n + \tfrac{h^2}{4} (a_n + a_{n+1} ) \big) \]

18.2. Generalized-\(\alpha\) method

Although the Newmark method preserves energy exactly for linear ODEs, it leads to instabilities for non-linear equations. For this reason, the following generalization is usually used in mechanical simulation codes:

By adding a bunch of parameters, one ends up with the Generalied \(\alpha\) method:

\[ M a_{n+1-\alpha_m} = F(x_{n+1-\alpha_f}) \]

with

\[\begin{split} \begin{array}{rcl} x_{n+1} & = & x_n + h v_n + h^2 ( ( \frac{1}{2} - \beta) a_n + \beta a_{n+1} ) \\ v_{n+1} & = & v_n + h ( (1-\gamma) a_n + \gamma a_{n+1} ) \\ x_{n+1-\alpha_f} & = & (1-\alpha_f) x_{n+1} + \alpha_f x_n \\ a_{n+1-\alpha_m} & = & (1-\alpha_m) a_{n+1} + \alpha_m a_n \end{array} \end{split}\]

Depending of one free parameter \(\rho^\infty\), parameters tuned for optimal accuracy and stabilty are obtained by

\[\begin{split} \begin{array}{rcl} \alpha_m & = & \frac{2 \rho^\infty - 1}{\rho^\infty + 1} \\ \alpha_f & = & \frac{\rho^\infty}{\rho^\infty + 1} \\ \beta & = & \frac{1}{4} (1 - \alpha_m + \alpha_f)^2 \\ \gamma & = & \frac{1}{2} - \alpha_m + \alpha_f \end{array} \end{split}\]

The parameter \(\rho^\infty\) specifies the damping for high frequency functions. The reference is J. Chung and G.M. Hulbert A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized-α Method, J. Appl. Mech., 1993, see An easy to read derivation is here.

18.3. Systems with constraints

Think of a pendulum, where a mass \(m\) is attached to a (mass-less) stick of length \(l\). One can describe the system by the angle \(\alpha\) between the stick and the vertical axis. By some high-school mathematics, one splits the total gravitational force into the direction of the stick, and the orthogonal complement, and ends up with the equation for the angle \(\alpha\):

\[ \ddot \alpha = -\frac{g}{l} \sin \alpha \]

However, such a direct modeling with angles (so called generalized coordinates) becomes difficult for more complex systems. A simpler possibility is to pose the length constraint

\[ 0 = g(x) := | x - x_0 | - l, \]

where \(x_0\) is the anchor point of the pendulum, and define the Lagrange function

\[ L(x, \lambda) = -U(x) + \left< \lambda , g(x) \right>, \]

with \(U(x) = m g x_z\), and pose the equation

\[\begin{split} m_i \ddot x_i & = & \frac{\partial}{\partial x_i} L(x, \lambda) \\ 0 & = & \nabla_\lambda L(x, \lambda) \end{split}\]

The physical meaning of the Lagrange parameter \(\lambda\) is the longitudinal force in the pendulum (precisely: \(\lambda \nabla_x g\) is the force vector). This is an ordinary differential equation, with algebraic constraints, a so called differential-algebraic equation (DAE).

One can diretly use the generalized-\(\alpha\) method for this DAE. The right hand side is the gradient of the Lagrange function, the mass matrix on the left hand side is extended by \(0\) for the Lagrange parameters:

\[\begin{split} \left( \begin{array}{cc} M & 0 \\ 0 & 0 \end{array} \right) \left( \begin{array}{cc} \ddot x \\ \ddot \lambda \end{array} \right) = \nabla L (x, \lambda) \end{split}\]

18.4. Exercises

18.4.1. Test Newmark and Generalized-\(\alpha\) solvers

• compile and run the cods test_newmark and test_alpha. Plot the solution functions

• test the equation for the pendulum using the angle \(\alpha\)

• implement a double pendulum as DAE

18.4.2. Mass-spring system

For small system it is reasonable to implement directly the right-hand-side function. For larger mechanical system it is of advantage to build some data structures describing the system. See files here. The class MassSpringSystem has containers for lists of masses, fixations and spring connections. By Python-bindings one can setup such a system by adding these components. The MSS_Functions implements the right-hand-side function for the ODE based on the MassSpringSystem model. There is

• file mass_spring.cc for testing a system from C++

• file test_mass_spring.py for testing the system from Python

• Jupyter-notebook file mass_spring.ipynb for testing the system including 3D visualization

Exercise:

• build and run the examples in the mass_spring directory

• add distance constraints to the MassSpring system. Enforce them using Lagrange multipliers, similar as in test_alpha.cc.

• the class MSS_Function implements the derivative using numerical differentiation. Implement the exact derivative.

• Experiment with different mechanical structures:

  • Extend the double-pendulum to a chain

  • Build more complex structures using beams (for example a crane), and simulate vibration

  • Replace the springs of the double-pendulum by distance-constraints

  • Build a spinning top (German: Kreisel). We know angular momentum keeps it spinning. Connect three masses by distance constraints modeling the spinning top. Start with velocities consistent to a rotating motion