21. 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 acting on the mass point \(i\).

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 multiplied by elongation, in the direction of the vector \(x_j-x_i\). The sum runs over all mass points \(j\) connected to mass point \(i\) with a spring.

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\) as \(\dot {x}\).

21.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{\tau}{2} (v_n + v_{n+1} ) \\ M (v_{n+1} - v_n) & = & \frac{\tau}{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{\tau}{2} (a_n + a_{n+1}) \\ x_{n+1} & = & x_n + \frac{\tau}{2} (v_n + v_{n+1}) = x_n + \tau v_n + \frac{\tau^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 + \tau v_n + \tfrac{\tau^2}{4} (a_n + a_{n+1} ) \big) \]

The dimension of the nonlinear system is the same as the number of state variables (and not twice the size, as it looks after introducing the first order system).

21.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 Generalized \(\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 + \tau v_n + \tau^2 ( ( \frac{1}{2} - \beta) a_n + \beta a_{n+1} ) \\ v_{n+1} & = & v_n + \tau ( (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.

21.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 |^2 - l^2, \]

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 directly 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 \\ 0 \end{array} \right) = \nabla L (x, \lambda) \end{split}\]

21.4. Exercises

21.4.1. Mass-spring system

Merge branch Newmark from ASC-ODE into your main branch.

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.cpp 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.

• 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