Heavy top constrained-rigid-body example
We provide here a simple application of the Kynema formulation for the
heavy-top problem, which is a rotating body fixed to the ground by a spherical
joint. It is a common benchmark problem for constrained-rigid-body dynamics and
for testing Lie-group time integrators like that used in Kynema. Our approach
follows much of what is described in in [@Bruls-etal:2012], but with the key
differences that we formulate the problem in inertial coordinates rather than
material coordinates and we include three translational DOFs for the spherical joint, \(\underline{u}_\mathrm{SJ}\). The generalized coordinates for this problem are
\[\begin{split}\underline{q} = \begin{bmatrix}
\underline{u}_\mathrm{HT} \\
\underline{\underline{R}}_\mathrm{HT} \\
\underline{u}_\mathrm{SJ}
\end{bmatrix}\end{split}\]
where \(\underline{q} = \mathbb{R}^{9}\).
We assume the heavy top is a thin disk with mass
\(m=15\) kg. The \(6\times6\) mass matrix in material coordinates is
\[\begin{split}\underline{\underline{M}}^*_\mathrm{HT} = \begin{bmatrix}
15 \mathrm{~kg}& 0 & 0 & 0 & 0 & 0\\
0 & 15 \mathrm{~kg} & 0 & 0 & 0 & 0\\
0 & 0 & 15 \mathrm{~kg} & 0 & 0 & 0\\
0 & 0 & 0 & 0.234375 \mathrm{~kg~m}^2 & 0 & 0\\
0 & 0 & 0 & 0 & 0.46875 \mathrm{~kg~m}^2 & 0\\
0 & 0 & 0 & 0 & 0 & 0.234375 \mathrm{~kg~m}^2 \\
\end{bmatrix}\end{split}\]
and the \(9\times9\) system mass matrix in material coordinates is
\[\begin{split}\underline{\underline{M}}^* = \begin{bmatrix}
\underline{\underline{M}}^*_\mathrm{HT} & \underline{\underline{0}}\\
\underline{\underline{0}} & \underline{\underline{0}}
\end{bmatrix}\end{split}\]
The heavy-top center of mass reference position and orientation (see
Eq. (5)) are given by
\[\begin{split}\underline{x}^\mathrm{r}_\mathrm{HT} = ( 0, 1 , 0 )^T\, \mathrm{m}, \quad
\underline{\underline{R}}^\mathrm{r}_\mathrm{HT} = \underline{\underline{I}} \,, \\\end{split}\]
respectively. The only component of external force (see
Eq. (6)) is gravity:
\[\underline{f} = [0,0,-g,0,0,0,0,0,0]^T\]
where \(g=9.81\) m/s\(^2\).
The problem is constrained such that the spherical joint is located at the global origin, the heavy-top center of mass is located 1 m from the
origin, and the heavy-top body-attached-coordinate-system y-axis is pointing away from the global origin. The Kynema implementation employs the following constraint:
\[\begin{split}\underline{\Phi} = \begin{bmatrix}
\underline{u}_\mathrm{HT} - \underline{u}_\mathrm{SJ}
- \left( \underline{\underline{R}}_\mathrm{HT} - \underline{\underline{I}}\right)
\underline{x}^r_\mathrm{HT} \\
\underline{u}_\mathrm{SJ}
\end{bmatrix}\end{split}\]
where \(\underline{\Phi} \in \mathbb{R}^6\), and the contraints \(\underline{\lambda} \in \mathbb{R}^6\) are denoted
\[\begin{split}\underline{\lambda} = \begin{bmatrix}
\underline{\lambda}_\mathrm{HT} \\
\underline{\lambda}_\mathrm{SJ}
\end{bmatrix}\end{split}\]
where \(\underline{\lambda}_\mathrm{HT},\, \underline{\lambda}_\mathrm{SJ}
\in \mathbb{R}^3\).
The constraint gradient matrix is
\[\begin{split}\underline{\underline{B}} = \begin{bmatrix}
\underline{\underline{I}} & \widetilde{ \underline{\underline{R}}_\mathrm{HT} \underline{x}^{\mathrm{r}}_\mathrm{HT} } & - \underline{\underline{I}} \\
\underline{\underline{0}} &
\underline{\underline{0}} &
\underline{\underline{I}}
\end{bmatrix}\end{split}\]
and the contribution to the iteration matrix in Eq. (4) is
\[\begin{split}\underline{\underline{K}}^\Phi = \begin{bmatrix}
\underline{\underline{0}} &
\underline{\underline{0}} &
\underline{\underline{0}} \\
\underline{\underline{0}} &
\widetilde{\lambda_\mathrm{HT}}
\widetilde{\underline{\underline{R}}_\mathrm{HT} \underline{x}^{\mathrm{r}}_\mathrm{HT} } & \underline{\underline{0}} \\
\underline{\underline{0}} &
\underline{\underline{0}} &
\underline{\underline{0}}
\end{bmatrix}\end{split}\]
Note
The term \(\underline{\underline{K}}^\Phi\) is include here for
completeness, but is not in the current Kynema implementation, and
nonlinear-system convergence has been satisfactory.
The Kynema regression test suite includes the spinning, heavy top
problem with the following initial conditions:
\[\begin{split}\begin{aligned}
\underline{u}^\mathrm{init}_\mathrm{HT} &= \left[ 0, 0, 0 \right]^T \,
\mathrm{m}\\
\underline{\underline{R}}_\mathrm{HT}^\mathrm{init} &= \underline{\underline{R}}^\mathrm{r}_\mathrm{HT} \\
\underline{u}^\mathrm{init}_\mathrm{SJ} &= \left[ 0, 0, 0 \right]^T \,
\mathrm{m}
\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned}
\dot{\underline{u}}^\mathrm{init}_\mathrm{HT} &=
\widetilde{\omega^\mathrm{init}_\mathrm{HT}}
\left(\underline{x}^\mathrm{r}_\mathrm{HT}
+\underline{u}^\mathrm{init}_\mathrm{HT}\right) \\
\omega^\mathrm{init}_\mathrm{HT} &= (0, 150, -4.61538)^T \,
\mathrm{rad/s}\\
\dot{\underline{u}}^\mathrm{init}_\mathrm{SJ} &= \left[ 0, 0, 0 \right]^T\,
\mathrm{m/s}
\end{aligned}\end{split}\]
Brüls, O., A. Cardona, and M. Arnold. 2012. “Lie Group
Generalized-\(\alpha\) Time Integration For Constrained
Flexible Multibody Systems.” Mechanism and Machine Theory,
121–37.