Solver
The constrained equations of motion are defined by
$$\begin{split}
\mathbf{M}(\mathbf{q})\,\ddot{\mathbf{q}}+\mathrm{\mathbf{C}}_{\mathbf{q}}^{\mathsf{T}}(\mathbf{q})\lambda + \mathrm{\mathbf{G}}_{\mathbf{q}}^{\mathsf{T}}(\mathbf{q})\mu &= \mathbf{Q}(\mathbf{q},\dot{\mathbf{q}},\mathbf{u},t)\\
\mathbf{C}(\mathbf{q}) & = 0 \\
\mathbf{G}(\mathbf{q}, t) & = 0
\end{split}$$
in which \(\mathbf{q}\in\mathbb{R}^n\) denotes the vector of redundant generalized coordinates, \(\mathbf{M}\in\mathbb{R}^{(n\times n)}\) the symmetric mass matrix and \(\mathbf{Q}\in\mathbb{R}^n\) the vector of generalized and gyroscopic forces. Due to the \(m\) algebraic constraints \(\mathbf{C}(\mathbf{q}) = 0\), the equations of motion are extended by constraint forces of the form \(\mathrm{\mathbf{C}}_{\mathbf{q}}^{\mathsf{T}}\lambda\), where \(\mathrm{\mathbf{C}}_{\mathbf{q}}\in\mathbb{R}^{(m\times n)}\) represents the constraint Jacobian and the vector \(\lambda\in\mathbb{R}^m\) includes the Lagrange multipliers. Another \(r\) time dependent algebraic constraint equation \(\mathbf{G}(\mathbf{q}, t)=0\) can be used to impose a user defined motion, which again leads to a constraint force \(\mathrm{\mathbf{G}}_{\mathbf{q}}^{\mathsf{T}}\mu\), where \(\mathrm{\mathbf{G}}_{\mathbf{q}}\in\mathbb{R}^{(r\times n)}\) is the Jacobian of the motion constraint and \(\mu\in\mathbb{R}^r\) a Lagrange multiplier. Moreover, \(\mathbf{u}\) is a vector of time dependent control signals actuating the system.
In multibody dynamics, Euler parameters are often used for the numerical simulation of rigid body rotations because they lead to a relatively simple form of the rotation matrix which avoids the evaluation of trigonometric functions and can thus save computational time. The Hilber-Hughes-Taylor (HHT) method is widely employed for solving the equations of motion of mechanical systems. However, in the classical versions, the use of these integration schemes have a very unfavorable impact on the Euler parameter description of rotational motions. This effect, which does not appear for Euler angles, can be even observed if the numerical damping parameter \alpha in the HHT-method is set to zero. To circumvent this problem without loosing the advantage of Euler parameters, we use a modified HHT-method which reduces the damping effect on the angular velocity significantly and eliminates it completely for \(\alpha = 0\).