6.3 Rigid Body Dynamics
Rigid Body Rotation
Exponential Map
Conventionally, the rigid body rotation tensor \(\boldsymbol{\Lambda}\) corresponding to a rotation of angle \(\chi\) about the unit vector \(\mathbf{n}\) may be expressed in terms of the vector \(\boldsymbol{\chi}=\chi\mathbf{n}\) as
\[ \begin{equation} \boldsymbol{\Lambda}\left(\boldsymbol{\chi}\right)=\cos\chi\,\mathbf{I}-\sin\chi\,\boldsymbol{\mathcal{E}}\cdot\mathbf{n}+\left(1-\cos\chi\right)\mathbf{n}\otimes\mathbf{n}\label{eq:rbr-rotation-tensor} \end{equation} \]
where \(\boldsymbol{\mathcal{E}}\) is the third-order permutation pseudo-tensor with Cartesian components \(\varepsilon_{ijk}\). Making use of the trigonometric identity,
\[ \begin{equation} \cos\chi=1-2\sin^{2}\frac{1}{2}\chi\,,\label{eq:rbr-identity-1} \end{equation} \]
this expression may be rearranged as
\[ \begin{equation} \boldsymbol{\Lambda}\left(\boldsymbol{\chi}\right)=\mathbf{I}-\frac{\sin\chi}{\chi}\,\boldsymbol{\mathcal{E}}\cdot\boldsymbol{\chi}+\frac{2}{\chi^{2}}\sin^{2}\frac{1}{2}\chi\left(\boldsymbol{\mathcal{E}}\cdot\boldsymbol{\chi}\right)^{2}\,,\label{eq:rbr-rotation-tensor-alt-1} \end{equation} \]
where we have made use of the identity
\[ \begin{equation} \left(\boldsymbol{\mathcal{E}}\cdot\boldsymbol{\chi}\right)^{2}=\boldsymbol{\chi}\otimes\boldsymbol{\chi}-\chi^{2}\mathbf{I}\,.\label{eq:rbr-identity-2} \end{equation} \]
Letting
\[ \begin{equation} \hat{\boldsymbol{\chi}}=-\boldsymbol{\mathcal{E}}\cdot\boldsymbol{\chi}\label{eq:rbr-antisymmetric-tensor} \end{equation} \]
represent the antisymmetric tensor with axial vector \(\boldsymbol{\chi}\), \(\boldsymbol{\Lambda}\left(\boldsymbol{\chi}\right)\) may now be represented as
\[ \begin{equation} \boldsymbol{\Lambda}\left(\boldsymbol{\chi}\right)\equiv\exp\left[\hat{\boldsymbol{\chi}}\right]=\mathbf{I}+\frac{\sin\chi}{\chi}\hat{\boldsymbol{\chi}}+\frac{2}{\chi^{2}}\sin^{2}\left(\frac{1}{2}\chi\right)\,\hat{\boldsymbol{\chi}}^{2}\,,\label{eq:rbr-rotation-tensor-alt-2} \end{equation} \]
where \(\exp\left[\hat{\boldsymbol{\chi}}\right]\) is known as the exponential map. Thus, the exponential map provides the rotation tensor for a rotation \(\chi\) about the unit vector \(\mathbf{n}\). Note that \(\boldsymbol{\Lambda}\cdot\boldsymbol{\chi}=\boldsymbol{\chi}\), since \(\hat{\boldsymbol{\chi}}\cdot\boldsymbol{\chi}=\boldsymbol{\chi}\times\boldsymbol{\chi}=\mathbf{0}\).
Let \(\mathbf{Q}\) be any orthogonal transformation, then
\[ \begin{equation} \begin{aligned}\mathbf{Q}\cdot\exp\left[\hat{\boldsymbol{\chi}}\right]\cdot\mathbf{Q}^{T} & =\mathbf{I}+\frac{\sin\chi}{\chi}\mathbf{Q}\cdot\hat{\boldsymbol{\chi}}\cdot\mathbf{Q}^{T}+\frac{2}{\chi^{2}}\sin^{2}\frac{1}{2}\chi\left(\mathbf{Q}\cdot\hat{\boldsymbol{\chi}}\cdot\mathbf{Q}^{T}\right)^{2}\\ & =\exp\left[\mathbf{Q}\cdot\hat{\boldsymbol{\chi}}\cdot\mathbf{Q}^{T}\right]\equiv\exp\left[\hat{\boldsymbol{\theta}}\right] \end{aligned} \label{eq:rbr-orthogonal-transformation} \end{equation} \]
where \(\hat{\boldsymbol{\theta}}=\mathbf{Q}\cdot\hat{\boldsymbol{\chi}}\cdot\mathbf{Q}^{T}\) and its corresponding axial vector is \(\boldsymbol{\theta}=\mathbf{Q}\cdot\boldsymbol{\chi}\), implying that \(\theta=\chi\). This property of the exponential map is used in the next derivation.
Consider a vector \(\mathbf{Z}\) in the reference configuration of a rigid body. Upon rigid body rotation, this vector is currently at
\[ \begin{equation} \mathbf{z}\left(t\right)=\boldsymbol{\Lambda}\left(t\right)\cdot\mathbf{Z}\,.\label{eq:rbr-rotation-t} \end{equation} \]
The corresponding axial vector of \(\boldsymbol{\Lambda}\left(t\right)\) is \(\boldsymbol{\chi}\left(t\right)\). At a subsequent time \(t^{\prime}\), we would similarly have
\[ \begin{equation} \mathbf{z}\left(t^{\prime}\right)=\boldsymbol{\Lambda}\left(t^{\prime}\right)\cdot\mathbf{Z}=\exp\left[\hat{\boldsymbol{\theta}}\right]\cdot\boldsymbol{\Lambda}\left(t\right)\cdot\mathbf{Z}\,,\label{eq:rbr-orthogonal-rotation-t-prime} \end{equation} \]
where here, \(\boldsymbol{\theta}\) is the incremental (finite) rotation from \(t\) to \(t^{\prime}\). Alternatively, we may choose to write
\[ \begin{equation} \mathbf{z}\left(t^{\prime}\right)=\boldsymbol{\Lambda}\left(t^{\prime}\right)\cdot\mathbf{Z}=\boldsymbol{\Lambda}\left(t\right)\cdot\exp\left[\hat{\boldsymbol{\Theta}}\right]\cdot\mathbf{Z}\,,\label{eq:rbr-incremental-finite-rotation} \end{equation} \]
such that
\[ \begin{equation} \begin{aligned}\boldsymbol{\Lambda}\left(t^{\prime}\right) & =\exp\left[\hat{\boldsymbol{\theta}}\right]\cdot\boldsymbol{\Lambda}\left(t\right)\\ & =\boldsymbol{\Lambda}\left(t\right)\cdot\exp\left[\hat{\boldsymbol{\Theta}}\right] \end{aligned} \,,\label{eq:rbr-material-spatial-increments} \end{equation} \]
implying that
\[ \begin{equation} \begin{aligned}\exp\left[\hat{\boldsymbol{\theta}}\right] & =\exp\left[\boldsymbol{\Lambda}\left(t\right)\cdot\hat{\boldsymbol{\Theta}}\cdot\boldsymbol{\Lambda}^{T}\left(t\right)\right]\\ \boldsymbol{\theta} & =\boldsymbol{\Lambda}\left(t\right)\cdot\boldsymbol{\Theta} \end{aligned} \,.\label{eq:rbr-material-spatial-redux} \end{equation} \]
Note from these relations that \(\theta=\Theta\). Thus, \(\boldsymbol{\Theta}\) is the material representation of the incremental rotation from \(t\) to \(t^{\prime}\), while \(\boldsymbol{\theta}\) is the corresponding spatial representation.
An alternative to the exponential map is the Cayley transform,
\[ \begin{equation} \boldsymbol{\Lambda}\left(\boldsymbol{\chi}\right)=\cay\left[\hat{\boldsymbol{\chi}}\right]=\mathbf{I}+\frac{2}{1+\left(\frac{1}{2}\chi\right)^{2}}\left(\frac{1}{2}\hat{\boldsymbol{\chi}}+\frac{1}{4}\hat{\boldsymbol{\chi}}^{2}\right)\label{eq:rbr-Cayley-transform} \end{equation} \]
which is a second order approximation to the exponential map. This formula is a correction to that appearing in (which has \(\frac{1}{2}\chi^{2}\) in the denominator). According to Puso , the Cayley transform must be used to enforce conservation of momentum and energy in a midpoint rule discretization scheme, whenever the rigid body is connected to a deformable body, or whenever two rigid bodies are connected by a joint. Comparing the above expression to the exponential map \(\exp\left[\hat{\boldsymbol{\theta}}\right]\) using \eqref{eq:rbr-rotation-tensor-alt-2}, we find that \(\boldsymbol{\chi}\) and \(\boldsymbol{\theta}\) are related via
\[ \begin{equation} \boldsymbol{\chi}=2\tan\frac{\theta}{2}\,\mathbf{n}\,.\label{eq:exponential-Cayley-relation} \end{equation} \]
Linearization Along Rotational Increment
Let \(\boldsymbol{\theta}\) represent a spatial rotational increment, such that a rotation tensor compounded by an infinitesimal incremental rotation is given by
\[ \begin{equation} \boldsymbol{\Lambda}_{\varepsilon}=\cay\left[\varepsilon\hat{\boldsymbol{\theta}}\right]\cdot\boldsymbol{\Lambda}\,.\label{eq:rbr-spatial-rotation-increment} \end{equation} \]
Using the Cayley transform for illustration, the linearization of \(\boldsymbol{\Lambda}\) along the increment \(\boldsymbol{\theta}\) is obtained from
\[ \begin{equation} \begin{aligned}D\boldsymbol{\Lambda}\left[\boldsymbol{\theta}\right] & =\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\cay\left[\varepsilon\hat{\boldsymbol{\theta}}\right]\cdot\boldsymbol{\Lambda}\\ & =\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\left(\mathbf{I}+\frac{2}{1+\frac{1}{4}\varepsilon^{2}\theta^{2}}\left(\frac{1}{2}\varepsilon\hat{\boldsymbol{\theta}}+\frac{1}{4}\varepsilon^{2}\hat{\boldsymbol{\theta}}^{2}\right)\right)\cdot\boldsymbol{\Lambda}\\ & =\hat{\boldsymbol{\theta}}\cdot\boldsymbol{\Lambda} \end{aligned} \,.\label{eq:rbr-spatial-rotation-linearization} \end{equation} \]
The same result may be obtained with the exponential map. Similarly, using an infinitesimal material rotational increment such that \(\boldsymbol{\Lambda}_{\varepsilon}=\boldsymbol{\Lambda}\cdot\cay\left[\varepsilon\hat{\boldsymbol{\Theta}}\right]\), we may find
\[ \begin{equation} D\boldsymbol{\Lambda}\left[\boldsymbol{\Theta}\right]=\boldsymbol{\Lambda}\cdot\hat{\boldsymbol{\Theta}}\,.\label{eq:rbr-material-rotation-increment} \end{equation} \]
General Rigid Body Motion
If the point \(\mathbf{x}\) is connected to a rigid body, its motion is given by
\[ \begin{equation} \mathbf{x}=\mathbf{r}\left(t\right)+\boldsymbol{\Lambda}\left(t\right)\cdot\mathbf{Z}=\mathbf{r}\left(t\right)+\mathbf{z}\left(t\right)\label{eq:rbm-position} \end{equation} \]
where \(\mathbf{r}\left(t\right)\) is the position of the rigid body center of mass and \(\boldsymbol{\Lambda}\left(t\right)\) is the body's rotation tensor, which satisfies \(\boldsymbol{\Lambda}\left(t_{0}\right)=\mathbf{I}\) at the initial time \(t_{0}\); here, \(\mathbf{z}\left(t\right)=\boldsymbol{\Lambda}\left(t\right)\cdot\mathbf{Z}\) is the distance of the point from the body's center of mass, and \(\mathbf{X}=\mathbf{r}\left(t_{0}\right)+\mathbf{Z}\) is the initial position. The velocity of that point is
\[ \begin{equation} \begin{aligned}\dot{\mathbf{x}} & =\dot{\mathbf{r}}\left(t\right)+\dot{\boldsymbol{\Lambda}}\left(t\right)\cdot\mathbf{Z}\end{aligned} \,,\label{eq:rbm-velocity} \end{equation} \]
where
\[ \begin{equation} \dot{\boldsymbol{\Lambda}}\left(t\right)=\hat{\boldsymbol{\omega}}\left(t\right)\cdot\boldsymbol{\Lambda}\left(t\right)=\boldsymbol{\Lambda}\left(t\right)\cdot\hat{\mathbf{W}}\left(t\right)\,.\label{eq:rbm-lambda-dot} \end{equation} \]
Here, \(\hat{\boldsymbol{\omega}}\) is an antisymmetric tensor with axial vector \(\boldsymbol{\omega}\) which represents the spatial angular velocity vector; similarly, \(\hat{\mathbf{W}}\) is an antisymmetric tensor with axial vector \(\mathbf{W}\) (the material angular velocity), such that \(\boldsymbol{\omega}=\boldsymbol{\Lambda}\cdot\mathbf{W}\) and
\[ \begin{equation} \hat{\boldsymbol{\omega}}=\boldsymbol{\Lambda}\cdot\hat{\mathbf{W}}\cdot\boldsymbol{\Lambda}^{T}\,.\label{eq:rbm-material-spatial-angular-velocity} \end{equation} \]
We may now rewrite
\[ \begin{equation} \begin{aligned}\dot{\mathbf{x}} & =\dot{\mathbf{r}}\left(t\right)+\hat{\boldsymbol{\omega}}\left(t\right)\cdot\mathbf{z}\left(t\right)\\ & =\dot{\mathbf{r}}\left(t\right)+\boldsymbol{\Lambda}\left(t\right)\cdot\hat{\mathbf{W}}\left(t\right)\cdot\mathbf{Z} \end{aligned} \,,\label{eq:rbm-rigid-point-velocity} \end{equation} \]
so that the acceleration of the point is
\[ \begin{equation} \begin{aligned}\ddot{\mathbf{x}} & =\ddot{\mathbf{r}}\left(t\right)+\left(\hat{\boldsymbol{\alpha}}\left(t\right)+\hat{\boldsymbol{\omega}}^{2}\left(t\right)\right)\cdot\mathbf{z}\\ & =\ddot{\mathbf{r}}\left(t\right)+\boldsymbol{\Lambda}\left(t\right)\cdot\left(\hat{\mathbf{A}}\left(t\right)+\hat{\mathbf{W}}^{2}\left(t\right)\right)\cdot\mathbf{Z} \end{aligned} \,,\label{eq:rbm-rigid-point-acceleration} \end{equation} \]
where \(\boldsymbol{\alpha}=\dot{\boldsymbol{\omega}}=\boldsymbol{\Lambda}\cdot\mathbf{A}\) is the spatial angular acceleration vector, \(\mathbf{A}=\dot{\mathbf{W}}\) is the material angular acceleration vector. As shown below, the time discretization is performed in the material frame.
Rigid Body Momentum Balance
For a rigid body, the conservation of linear momentum is given by
\[ \begin{equation} \frac{d}{dt}\left(m\dot{\mathbf{r}}\right)=\dot{\mathbf{p}}=\mathbf{f}^{ext}\left(t\right)\label{eq:rbm-linear-momentum-balance} \end{equation} \]
where \(m\) is the mass of the rigid body, \(\dot{\mathbf{r}}\) is the velocity of the center of mass, \(\mathbf{p}=m\dot{\mathbf{r}}\) is the linear momentum, and \(\mathbf{f}^{ext}\left(t\right)\) represents the sum of external forces acting on the body. Here, \(m\) is constant for a rigid body. There are typically four contributions to \(\mathbf{f}^{ext}\left(t\right)\): Body forces \(\mathbf{f}_{b}^{ext}\left(t\right)=m\mathbf{b}\left(t\right)\) (where \(\mathbf{b}\) represents the body force per mass, such as gravitational acceleration), other user-prescribed forces \(\mathbf{f}_{p}^{ext}\left(t\right)\) (which act at the center of mass), forces \(\mathbf{f}_{c}^{ext}\left(t\right)\) produced by rigid body connectors (such as revolute and prismatic joints, or contact forces), and forces \(\mathbf{f}_{f}^{ext}\left(t\right)\) produced by rigid-flexible connections (where deformable materials interface with the rigid body), in which case \(\mathbf{f}_{f}^{ext}\) is evaluated from the traction \(\mathbf{t}=\boldsymbol{\sigma}\cdot\mathbf{n}\) over that interface, with \(\boldsymbol{\sigma}\) representing the stress in the deformable material.
The conservation of angular momentum is similarly given by
\[ \begin{equation} \frac{d}{dt}\left(\mathbf{J}\cdot\boldsymbol{\omega}\right)=\dot{\mathbf{h}}=\boldsymbol{\omega}\times\mathbf{h}+\mathbf{J}\cdot\boldsymbol{\alpha}=\mathbf{m}^{ext}\left(t\right)\label{eq:rbm-angular-momentum-balance} \end{equation} \]
where \(\mathbf{J}\) is the rigid body mass moment of inertia about its center of mass, \(\boldsymbol{\omega}\) is its angular velocity, \(\mathbf{h}=\mathbf{J}\cdot\boldsymbol{\omega}\) is its angular momentum, \(\boldsymbol{\alpha}=\dot{\boldsymbol{\omega}}\) is the rigid body angular acceleration, and \(\mathbf{m}^{ext}\left(t\right)\) is the sum of moments acting on the rigid body. External moments include contributions from user-prescribed moments/torques \(\mathbf{m}_{p}^{ext}\left(t\right)\), from rigid body connectors, \(\mathbf{m}_{c}^{ext}\left(t\right)=\mathbf{z}_{c}\left(t\right)\times\mathbf{f}_{c}^{ext}\left(t\right)\) where \(\mathbf{z}_{c}\left(t\right)\) is the connector insertion relative to the rigid body center of mass, and rigid-flexible interfaces, \(\mathbf{m}_{f}^{ext}\left(t\right)=\mathbf{z}_{f}\left(t\right)\times\mathbf{f}_{f}^{ext}\left(t\right)\) where \(\mathbf{z}_{f}\left(t\right)\) is the position of the interface point relative to the rigid body center of mass. Since body forces \(\mathbf{f}_{b}^{ext}\) and user-prescribed forces \(\mathbf{f}_{p}^{ext}\) act at the center of mass, they do not contribute to \(\mathbf{m}^{ext}\left(t\right)\). Note that
\[ \begin{equation} \mathbf{J}\left(t\right)=\boldsymbol{\Lambda}\left(t\right)\cdot\mathbf{J}_{r}\cdot\boldsymbol{\Lambda}^{T}\left(t\right)\,,\label{eq:rbm-moi-rotation} \end{equation} \]
where \(\mathbf{J}_{r}\) is the mass moment of inertia about the center of mass in the reference configuration and \(\boldsymbol{\Lambda}\left(t\right)\) is the rotation tensor representing the orientation of the rigid body at time \(t\), with \(\boldsymbol{\Lambda}=\mathbf{I}\) in the reference configuration.
The virtual work statement is given by
\[ \begin{equation} \delta W=\delta\mathbf{r}\cdot\left(\mathbf{f}^{ext}\left(t\right)-\dot{\mathbf{p}}\right)+\delta\boldsymbol{\theta}\cdot\left(\mathbf{m}^{ext}\left(t\right)-\dot{\mathbf{h}}\right)\,,\label{eq:rbm-virtual-work} \end{equation} \]
where \(\delta\mathbf{r}\) is the virtual velocity of the center of mass and \(\delta\boldsymbol{\theta}\) is the virtual angular velocity of the rigid body.
Time Discretization
Newmark Integration for Rigid Body Dynamics
Let \(t_{n}\) and \(t_{n+1}\) represent consecutive time points. According to the Newmark integration scheme, the rigid body center of mass velocity and acceleration at \(t_{n+1}\) may be expressed in terms of their values at \(t_{n}\) as
\[ \begin{equation} \begin{aligned}\dot{\mathbf{r}}_{n+1} & =\dot{\mathbf{r}}_{n}+\Delta t\left[\left(1-\gamma\right)\ddot{\mathbf{r}}_{n}+\gamma\ddot{\mathbf{r}}_{n+1}\right]\\ \ddot{\mathbf{r}}_{n+1} & =\frac{1}{\beta\Delta t}\left[\frac{1}{\Delta t}\left(\mathbf{r}_{n+1}-\mathbf{r}_{n}\right)-\dot{\mathbf{r}}_{n}\right]+\left(1-\frac{1}{2\beta}\right)\ddot{\mathbf{r}}_{n} \end{aligned} \,,\label{eq:td-Newmark-COM} \end{equation} \]
where \(\beta\) and \(\gamma\) are Newmark parameters that satisfy \(0\le2\beta\le1\) and \(0\le\gamma\le1\).
Let the rigid body rotation tensor \(\boldsymbol{\Lambda}\left(t\right)\) be expressed as \(\boldsymbol{\Lambda}\left(t\right)=\exp\left[\boldsymbol{\xi}\left(t\right)\right]\), and \(\boldsymbol{\xi}\left(t\right)\) is the material rotation of the rigid body from its reference configuration. Thus, \(\boldsymbol{\Lambda}_{n}=\exp\left[\boldsymbol{\xi}_{n}\right]\) and \(\boldsymbol{\Lambda}_{n+1}=\exp\left[\boldsymbol{\xi}_{n+1}\right]\) respectively represent the rigid body rotation tensors at \(t_{n}\) and \(t_{n+1}\). (In practice, \(\boldsymbol{\xi}\) is stored as a quaternion to facilitate the multiplication of rotation tensors.) These tensors are related by the incremental spatial rotation \(\boldsymbol{\theta}\) or material rotation \(\boldsymbol{\Theta}\) from \(t_{n}\) to \(t_{n+1}\) according to
\[ \begin{equation} \boldsymbol{\Lambda}_{n+1}=\cay\left[\boldsymbol{\theta}\right]\cdot\boldsymbol{\Lambda}_{n}=\boldsymbol{\Lambda}_{n}\cdot\cay\left[\boldsymbol{\Theta}\right]\,.\label{eq:td-incremental-rotation} \end{equation} \]
Here, it should be understood that the material frame for this incremental rotation is the configuration at time \(t_{n}\), while the spatial frame is the configuration at \(t_{n+1}\). For rotational motion, the Newmark scheme is applied in the material frame as
\[ \begin{equation} \begin{aligned}\mathbf{W}_{n+1} & =\frac{\gamma}{\beta\Delta t}\boldsymbol{\Theta}-\mathbf{W}_{n}+\left(2-\frac{\gamma}{\beta}\right)\left(\mathbf{W}_{n}+\frac{\Delta t}{2}\mathbf{A}_{n}\right)\\ \mathbf{A}_{n+1} & =\frac{1}{\gamma\Delta t}\left(\mathbf{W}_{n+1}-\mathbf{W}_{n}\right)+\left(1-\frac{1}{\gamma}\right)\mathbf{A}_{n}\\ & =\frac{1}{\beta\Delta t}\left(\frac{1}{\Delta t}\boldsymbol{\Theta}-\mathbf{W}_{n}\right)+\left(1-\frac{1}{2\beta}\right)\mathbf{A}_{n} \end{aligned} \,.\label{eq:td-Newmark-mrot} \end{equation} \]
Then, using the relations \(\boldsymbol{\omega}=\boldsymbol{\Lambda}\cdot\mathbf{W}\) and \(\boldsymbol{\alpha}=\boldsymbol{\Lambda}\cdot\mathbf{A}\) at \(t_{n}\) and \(t_{n+1}\), along with \eqref{eq:td-incremental-rotation}, we may express these relations in the spatial frame as
\[ \begin{equation} \begin{aligned}\boldsymbol{\omega}_{n+1} & =\cay\left[\boldsymbol{\theta}\right]\cdot\left(\frac{\gamma}{\beta\Delta t}\boldsymbol{\theta}-\boldsymbol{\omega}_{n}+\left(2-\frac{\gamma}{\beta}\right)\left(\boldsymbol{\omega}_{n}+\frac{\Delta t}{2}\boldsymbol{\alpha}_{n}\right)\right)\\ \boldsymbol{\alpha}_{n+1} & =\cay\left[\boldsymbol{\theta}\right]\cdot\left(\frac{1}{\beta\Delta t}\left(\frac{1}{\Delta t}\boldsymbol{\theta}-\boldsymbol{\omega}_{n}\right)+\left(1-\frac{1}{2\beta}\right)\boldsymbol{\alpha}_{n}\right) \end{aligned} \,.\label{eq:td-Newmark-srot} \end{equation} \]
In a nonlinear solution scheme we solve for \(\boldsymbol{\Theta}\) incrementally. According to \eqref{eq:rbr-material-rotation-increment}, the linearization of \(\cay\left[\boldsymbol{\Theta}\right]\) along an increment \(\Delta\boldsymbol{\Theta}\) is given by
\[ \begin{equation} D\left(\cay\left[\boldsymbol{\Theta}\right]\right)\left[\Delta\boldsymbol{\Theta}\right]=\cay\left[\boldsymbol{\Theta}\right]\cdot\widehat{\Delta\boldsymbol{\Theta}}\,,\label{eq:td-linearization-rot} \end{equation} \]
so that
\[ \begin{equation} D\boldsymbol{\Lambda}_{n+1}\left[\Delta\boldsymbol{\Theta}\right]=\boldsymbol{\Lambda}_{n+1}\cdot\widehat{\Delta\boldsymbol{\Theta}}\,.\label{eq:td-linearization-Lambda} \end{equation} \]
The linearizations of \(\mathbf{W}_{n+1}\) and \(\mathbf{A}_{n+1}\), as given in \eqref{eq:td-Newmark-mrot}, along an increment \(\Delta\boldsymbol{\Theta}\) requires us to first evaluate \(D\boldsymbol{\Theta}\left[\Delta\boldsymbol{\Theta}\right]\). According to Puso ,
\[ \begin{equation} D\boldsymbol{\Theta}\left[\Delta\boldsymbol{\Theta}\right]=\mathbf{T}\left(\boldsymbol{\Theta}\right)\cdot\Delta\boldsymbol{\Theta}\,,\label{td-Theta-linearization} \end{equation} \]
where
\[ \begin{equation} \mathbf{T}\left(\boldsymbol{\Theta}\right)=\mathbf{I}+\frac{1}{2}\hat{\boldsymbol{\Theta}}+\frac{1}{4}\boldsymbol{\Theta}\otimes\boldsymbol{\Theta}\,.\label{eq:td-T-Theta-relation} \end{equation} \]
Thus,
\[ \begin{equation} \begin{aligned}D\mathbf{W}_{n+1}\left[\Delta\boldsymbol{\Theta}\right] & =\frac{\gamma}{\beta\Delta t}\mathbf{T}\left(\boldsymbol{\Theta}\right)\cdot\Delta\boldsymbol{\Theta}\\ D\mathbf{A}_{n+1}\left[\Delta\boldsymbol{\Theta}\right] & =\frac{1}{\beta\Delta t^{2}}\mathbf{T}\left(\boldsymbol{\Theta}\right)\cdot\Delta\boldsymbol{\Theta} \end{aligned} \,.\label{eq:tf-W-A-linearizations} \end{equation} \]
Generalized-\(\alpha\) Method for Rigid Body Dynamics
In the generalized-\(\alpha\) method, we evaluate forces and moments at time \(t_{n+\alpha_{f}}=\left(1-\alpha_{f}\right)t_{n}+\alpha_{f}t_{n+1}\) and the time rate of change of linear and angular momenta at time \(t_{n+\alpha_{m}}=\left(1-\alpha_{m}\right)t_{n}+\alpha_{m}t_{n+1}\), where \(\alpha_{f}\) and \(\alpha_{m}\) may be evaluated from the spectral radius for an infinite time step, \(\rho_{\infty}\) (Section Generalized \(\alpha-\)Method). For second-order systems these parameters may be evaluated from
\[ \begin{equation} \alpha_{f}=\frac{1}{1+\rho_{\infty}}\,,\quad\alpha_{m}=\frac{2-\rho_{\infty}}{1+\rho_{\infty}}\,,\label{eq:ga-alphas-2nd} \end{equation} \]
Then, the Newmark parameters are given by
\[ \begin{equation} \begin{aligned}\beta & =\frac{1}{4}\left(1+\alpha_{m}-\alpha_{f}\right)^{2}\,,\\ \gamma & =\frac{1}{2}+\alpha_{m}-\alpha_{f}\,. \end{aligned} \label{eq:ga-Newmark-params} \end{equation} \]
Accordingly, to solve numerically for \(\delta W=0\) over the time domain, we express \eqref{eq:rbm-virtual-work} in the discretized time domain as
\[ \begin{equation} \delta\mathbf{r}\cdot\left(\mathbf{f}_{n+\alpha_{f}}^{ext}-\dot{\mathbf{p}}_{n+\alpha_{m}}\right)+\delta\boldsymbol{\theta}\cdot\left(\mathbf{m}_{n+\alpha_{f}}^{ext}-\dot{\mathbf{h}}_{n+\alpha_{m}}\right)=0\,,\label{eq:ga-virtual-work} \end{equation} \]
or equivalently,
\[ \begin{equation} \left[\begin{array}{cc} \delta\mathbf{r} & \delta\boldsymbol{\theta}\end{array}\right]\cdot\left[\begin{array}{c} \mathbf{f}_{n+\alpha_{f}}^{ext}-\dot{\mathbf{p}}_{n+\alpha_{m}}\\ \mathbf{m}_{n+\alpha_{f}}^{ext}-\dot{\mathbf{h}}_{n+\alpha_{m}} \end{array}\right]=0\,,\label{eq:ga-virtual-work-alt} \end{equation} \]
Thus, the residual vector is given by
\[ \begin{equation} \left[\mathbf{R}\right]=\left[\begin{array}{c} \mathbf{f}_{n+\alpha_{f}}^{ext}\\ \mathbf{m}_{n+\alpha_{f}}^{ext} \end{array}\right]-\left[\begin{array}{c} \dot{\mathbf{p}}_{n+\alpha_{m}}\\ \dot{\mathbf{h}}_{n+\alpha_{m}} \end{array}\right]\label{eq:ga-residual-vector} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\dot{\mathbf{p}}_{n+\alpha_{m}} & =\left(1-\alpha_{m}\right)\dot{\mathbf{p}}_{n}+\alpha_{m}\dot{\mathbf{p}}_{n+1}\\ \dot{\mathbf{h}}_{n+\alpha_{m}} & =\left(1-\alpha_{m}\right)\dot{\mathbf{h}}_{n}+\alpha_{m}\dot{\mathbf{h}}_{n+1} \end{aligned} \,.\label{eq:ga-momenta-interpolation} \end{equation} \]
According to the Newmark integration scheme,
\[ \begin{equation} \begin{aligned}\dot{\mathbf{p}}_{n+1} & =\frac{\mathbf{p}_{n+1}-\mathbf{p}_{n}}{\gamma\Delta t}+\left(1-\frac{1}{\gamma}\right)\dot{\mathbf{p}}_{n}\\ \dot{\mathbf{h}}_{n+1} & =\frac{\mathbf{h}_{n+1}-\mathbf{h}_{n}}{\gamma\Delta t}+\left(1-\frac{1}{\gamma}\right)\dot{\mathbf{h}}_{n} \end{aligned} \,,\label{eq:ga-momenta-Newmark} \end{equation} \]
where \(\mathbf{p}_{n+1}=m\dot{\mathbf{r}}_{n+1}\) and \(\mathbf{h}_{n+1}=\mathbf{J}_{n+1}\cdot\boldsymbol{\omega}_{n+1}\).
The nonlinear system \(\mathbf{R}=\mathbf{0}\) is solved using a Newton scheme that requires linearizing \(\mathbf{R}\) along increments \(\Delta\mathbf{r}\) and \(\Delta\boldsymbol{\theta}\). Thus,
\[ \begin{equation} \mathbf{R}+D\mathbf{R}\left[\Delta\mathbf{r}\right]+D\mathbf{R}\left[\Delta\boldsymbol{\theta}\right]\approx\mathbf{0}\,.\label{eq:ga-Newton-scheme} \end{equation} \]
The increments \(\Delta\mathbf{r}\) and \(\Delta\boldsymbol{\theta}\) are evaluated at \(t_{n+1}\) and the iterative Newton scheme requires updates of the form
\[ \begin{equation} \begin{aligned}\mathbf{r}_{n+1}^{j+1} & =\mathbf{r}_{n+1}^{j}+\Delta\mathbf{r}\\ \cay\left[\boldsymbol{\theta}^{j+1}\right] & =\cay\left[\Delta\boldsymbol{\theta}\right]\cdot\cay\left[\boldsymbol{\theta}^{j}\right] \end{aligned} \,,\label{eq:ga-Newton-iterations} \end{equation} \]
where \(j\) represents the Newton iteration. At each Newton iteration, the current value of \(\cay\left[\boldsymbol{\theta}^{j+1}\right]\) is used to perform the update
\[ \begin{equation} \boldsymbol{\Lambda}_{n+1}^{j+1}=\cay\left[\boldsymbol{\theta}^{j+1}\right]\cdot\boldsymbol{\Lambda}_{n}\,,\label{eq:ga-rotation-iterations} \end{equation} \]
until convergence is achieved.
In practice, it is convenient to store \(\boldsymbol{\theta}\) and \(\Delta\boldsymbol{\theta}\) in quaternions, recognizing that
\[ \cay\left[\theta\mathbf{n}\right]=\exp\left[\left(2\tan^{-1}\frac{\theta}{2}\right)\mathbf{n}\right]\,, \]
where \(\mathbf{n}\) is the unit vector along \(\boldsymbol{\theta}\) and \(\theta=\left\Vert \boldsymbol{\theta}\right\Vert\). Thus, it is \(\left(2\tan^{-1}\frac{\theta}{2}\right)\mathbf{n}\) which is stored in the quaternion, instead of \(\theta\mathbf{n}\).
In the linearization of \(\mathbf{R}\), the contributions from the rate of change of linear momentum \(\dot{\mathbf{p}}_{n+\alpha_{m}}\) reduce to
\[ \begin{aligned}D\dot{\mathbf{p}}_{n+\alpha_{m}}\left[\Delta\mathbf{r}\right] & =\frac{\alpha_{m}}{\beta\Delta t^{2}}m\Delta\mathbf{r}\\ D\dot{\mathbf{p}}_{n+\alpha_{m}}\left[\Delta\boldsymbol{\theta}\right] & =\mathbf{0} \end{aligned} \,. \]
To evaluate the contributions from the rate of change of angular momentum \(\dot{\mathbf{h}}_{n+\alpha_{m}}\), we start from \(D\dot{\mathbf{h}}_{n+\alpha_{m}}=\alpha_{m}D\dot{\mathbf{h}}_{n+1}\). Then, it becomes necessary to transform the variables to the material frame,
\[ \begin{aligned}\dot{\mathbf{h}}_{n+1} & =\boldsymbol{\omega}_{n+1}\times\mathbf{h}_{n+1}+\mathbf{J}_{n+1}\cdot\boldsymbol{\alpha}_{n+1}\\ & =\boldsymbol{\Lambda}_{n+1}\cdot\left(\hat{\mathbf{W}}_{n+1}\cdot\mathbf{J}_{r}\cdot\mathbf{W}_{n+1}+\mathbf{J}_{r}\cdot\mathbf{A}_{n+1}\right) \end{aligned} \,. \]
It follows that
\[ D\dot{\mathbf{h}}_{n+\alpha_{m}}\left[\Delta\mathbf{r}\right]=\mathbf{0}\,. \]
Then, using the relations in Section Newmark Integration for Rigid Body Dynamics, it can be shown that
\[ D\dot{\mathbf{h}}_{n+\alpha_{m}}\left[\Delta\boldsymbol{\Theta}\right]=\alpha_{m}\left[\frac{1}{\beta\Delta t}\left(\left(\gamma\hat{\boldsymbol{\omega}}_{n+1}+\frac{1}{\Delta t}\mathbf{I}\right)\cdot\mathbf{J}_{n+1}-\gamma\hat{\mathbf{h}}_{n+1}\right)\cdot\mathbf{T}\left(\boldsymbol{\theta}\right)-\hat{\dot{\mathbf{h}}}_{n+1}\right]\cdot\Delta\boldsymbol{\theta}\equiv\alpha_{m}\mathbf{K}\cdot\Delta\boldsymbol{\theta}\,. \]
Alternatively, we may use the discretization in \eqref{eq:ga-momenta-Newmark} to produce
\[ D\dot{\mathbf{h}}_{n+\alpha_{m}}\left[\Delta\boldsymbol{\Theta}\right]=\frac{\alpha_{m}}{\gamma\Delta t}D\mathbf{h}_{n+1}\left[\Delta\boldsymbol{\Theta}\right] \]
where
\[ D\mathbf{h}_{n+1}\left[\Delta\boldsymbol{\Theta}\right]=\left(\frac{\gamma}{\beta\Delta t}\mathbf{J}_{n+1}\cdot\mathbf{T}\left(\boldsymbol{\theta}\right)-\hat{\mathbf{h}}_{n+1}\right)\cdot\Delta\boldsymbol{\theta} \]
so that
\[ D\dot{\mathbf{h}}_{n+\alpha_{m}}\left[\Delta\boldsymbol{\Theta}\right]=\frac{\alpha_{m}}{\Delta t}\left(\frac{1}{\beta\Delta t}\mathbf{J}_{n+1}\cdot\mathbf{T}\left(\boldsymbol{\theta}\right)-\frac{1}{\gamma}\hat{\mathbf{h}}_{n+1}\right)\cdot\Delta\boldsymbol{\theta}\equiv\alpha_{m}\mathbf{K}\cdot\Delta\boldsymbol{\theta} \]
Therefore, the contribution to \(D\mathbf{R}\) from the linear and angular momenta produces a stiffness matrix called the mass matrix,
\[ D\left[\begin{array}{c} \dot{\mathbf{p}}_{n+\alpha_{m}}\\ \dot{\mathbf{h}}_{n+\alpha_{m}} \end{array}\right]=\alpha_{m}\left[\begin{array}{cc} \frac{m}{\beta\Delta t^{2}}\mathbf{I} & \mathbf{0}\\ \mathbf{0} & \mathbf{K} \end{array}\right]\left[\begin{array}{c} \Delta\mathbf{r}\\ \Delta\boldsymbol{\theta} \end{array}\right]\,. \]
Consider that the moments \(\mathbf{m}_{n+\alpha_{f}}^{ext}\) about the rigid body center of mass are produced by the forces \(\mathbf{f}_{n+\alpha_{f}}^{ext}\) according to
\[ \mathbf{m}_{n+\alpha_{f}}^{ext}=\mathbf{z}_{n+\alpha_{f}}\times\mathbf{f}_{n+\alpha_{f}}^{ext}=\hat{\mathbf{z}}_{n+\alpha_{f}}\cdot\mathbf{f}_{n+\alpha_{f}}^{ext}\,, \]
where \(\mathbf{z}_{n+\alpha_{f}}\) is the moment arm at \(t_{n+\alpha_{f}}\),
\[ \begin{aligned}\mathbf{z}_{n+\alpha_{f}} & =\left(1-\alpha_{f}\right)\mathbf{z}_{n}+\alpha_{f}\mathbf{z}_{n+1}\\ & =\left[\left(1-\alpha_{f}\right)\boldsymbol{\Lambda}_{n}+\alpha_{f}\boldsymbol{\Lambda}_{n+1}\right]\cdot\mathbf{Z}\\ & \equiv\boldsymbol{\Lambda}_{n+\alpha_{f}}\cdot\mathbf{Z} \end{aligned} \]
where \(\mathbf{Z}\) is the moment arm in the reference configuration. Note that
\[ \begin{aligned}D\mathbf{z}_{n+\alpha_{f}}\left[\Delta\mathbf{r}\right] & =\mathbf{0}\\ D\mathbf{z}_{n+\alpha_{f}}\left[\Delta\boldsymbol{\theta}\right] & =\alpha_{f}D\boldsymbol{\Lambda}_{n+1}\left[\Delta\boldsymbol{\theta}\right]\cdot\mathbf{Z}=-\alpha_{f}\hat{\mathbf{z}}_{n+1}\cdot\Delta\boldsymbol{\theta} \end{aligned} \]
Thus,
\[ \left[\begin{array}{c} D\mathbf{f}_{n+\alpha_{f}}^{ext}\left[\Delta\mathbf{r}\right]\\ D\mathbf{m}_{n+\alpha_{f}}^{ext}\left[\Delta\mathbf{r}\right] \end{array}\right]=\left[\begin{array}{c} D\mathbf{f}_{n+\alpha_{f}}^{ext}\left[\Delta\mathbf{r}\right]\\ \mathbf{z}_{n+\alpha_{f}}\times D\mathbf{f}_{n+\alpha_{f}}^{ext}\left[\Delta\mathbf{r}\right] \end{array}\right] \]
and
\[ \left[\begin{array}{c} D\mathbf{f}_{n+\alpha_{f}}^{ext}\left[\Delta\boldsymbol{\theta}\right]\\ D\mathbf{m}_{n+\alpha_{f}}^{ext}\left[\Delta\boldsymbol{\theta}\right] \end{array}\right]=\left[\begin{array}{c} D\mathbf{f}_{n+\alpha_{f}}^{ext}\left[\Delta\boldsymbol{\theta}\right]\\ \mathbf{z}_{n+\alpha_{f}}\times D\mathbf{f}_{n+\alpha_{f}}^{ext}\left[\Delta\boldsymbol{\theta}\right]-\left(\alpha_{f}\hat{\mathbf{z}}_{n+1}\cdot\Delta\boldsymbol{\theta}\right)\times\mathbf{f}_{n+\alpha_{f}}^{ext} \end{array}\right]\,. \]