Sage Journals: Discover world-class research

Abstract

The Euler–Poinaré principle is a reduced Hamilton’s principle under Lie group framework. In this article, it is applied to derive a hybrid set of dynamical equations of rigid multibody systems, which include four parts: the classical Euler–Lagrange equations of rigid bodies in their translational coordinates of mass center; Euler–Poinaré equations via orientation matrices and their related angular velocities; the constraint equations due to different joints in Cartesian coordinates and Lie groups; and the reconstruction equations between special orthogonal groups and their Lie algebras. The generalized mass matrices of dynamical equations are constant, which is computationally efficient. All the equations can be constructed systematically and can be solved easily. The construction equations can be used to design Lie group integrators of multibody system dynamics. The procedure presented in this article can be extended easily to flexible multibody systems, systems with non-holonomic constraints, and so on.

Keywords

Multibody system dynamics Euler–Lagrange equations Lie group Euler–Poinaré equations reconstruction equations constraint equations

Introduction

Multibody systems are mechanical systems that include interconnected bodies through different joints driven by some types of forces; their dynamical equations are usually denoted by nonlinear highly coupled ordinary differential equations or differential/algebraic equations, which must be solved numerically. The computational efficiency and stability of these dynamical equations depend strongly on the configuration description of multibody systems^1–9 and especially on the expressions of the orientations of bodies.

The classical Euler angles are usually used to describe the rotations of bodies in multibody systems,¹ but the resulting dynamical equations of motion are nonlinear inertia coupled with singular problems. If the Euler angles are replaced by the Euler parameters,^2–6 the singular problems can be overcome. However, the nonlinear inertia coupling still exists. In addition, we must introduce new constraints in order to satisfy the definition of Euler parameters. If the three directors from direction cosine matrix of every body are used to express the orientation,¹⁰ the generalized mass matrices of the dynamical equations are constant, which contributes to high computational efficiency, but six new constraint equations must be introduced for each body which enlarges the computational scales. If we use the orientation matrices as special orthogonal groups to describe the rotations,^11–14 the motions of bodies in multibody systems can be directly denoted by the Lie groups and their corresponding Lie algebras. Meanwhile, we do not need any additional constraints for the rotation description. Noticeably, the resulting generalized mass matrices are constant. This Lie group modeling method has become one of foundations for structure-preserving algorithms.

Technically, different descriptions of system configuration lead to different types of dynamical equations even though we use same modeling principle. In vectorized parameterization space, the Euler–Lagrange equations or Hamilton’s equations of multibody system dynamics can be derived under the Hamilton’s principle. While under the framework of Lie group and Lie algebra, these equations are the reduced Euler–Poinaré and Lie–Poisson equations.^15–17 By making use of the integrators of Lie group,^18–21 the Euler–Poinaré and Lie–Poisson equations can be easily applied to design the structure-preserving integrators of Lie group. These integrators have been implemented for long-term simulation of free rigid body as well as full rigid body.²² Additionally, by combining Lie group methods and the variational principle of discrete mechanics,^23,24 Lee et al.²⁵ and Bou-Rabee and Marsden²⁶ proposed the more stable algorithms of symplecticity–momentum–energy-Lie group preserving. Mizzi,²⁷ Liu,²⁸ and Boyer and Primault²⁹ used Lie group methods to construct equations for the motion of multibody systems. Unfortunately, the motions of rotation are parameterized finally. Park and Chung,³⁰ Brüls et al.,³¹ and Sonneville and Brüls^32,33 focused on Lie group integrators of multibody system dynamics, but they did not present the systematic derivations of dynamic equations under Lie group framework. In this article, based on Euler–Poinaré theory, we focus on how to derive the dynamical equations of rigid multibody systems by using Lie group and Lie algebra expressions of rotations.

The rest of this article is as follows. The classic Euler–Poinaré equations of free rigid body are introduced in section “The classical Euler–Poinaré equations of a rigid body.” In section “The Euler–Poinaré equations of rigid multibody system dynamics,” we derive the Euler–Poinaré equations of rigid multibody systems with joint constraints. In section “Generalized external forces in Lie groups and Lie algebras,” we analyze some generalized external forces under Lie group and Lie algebra. Section “Examples” presents two examples using the proposed method. In the last section, the contributions of this article are summarized and some future research topics are also presented.

The classical Euler–Poinaré equations of a rigid body

Lie group is a smooth differentiable manifold on which composition operations such as matrix product can be defined. The Euler–Poinaré principle is based on matrix Lie group. If Lie group $g$ is used to describe the configuration space of a dynamic system, its corresponding Lie algebras $ξ$ and $η$ can be denoted as^15,16

ξ = g^{- 1} \overset{\cdot}{g}

(1)

η = g^{- 1} δ g

(2)

where $ξ$ and $η$ are skew-symmetric matrices in vector space and $δ g$ is the variation in $g$ . For the problem of fixed time terminals, $δ g (0) = δ g (T) = 0$ , $η (0) = η (T) = 0$ , and $t \in [0, T]$ . $ξ$ and $η$ satisfy the following relationship of constrained variation

δ ξ = \overset{\cdot}{η} + a d_{ξ} η

(3)

where $a d_{ξ} η = [ξ, η] = ξ η - η ξ$ is the commutation operator between two Lie algebras $ξ$ and $η$ . Using Lie group and Lie algebra, the reduced Hamilton action is given as

\begin{matrix} S & = \int_{0}^{T} L (g, \overset{\cdot}{g}) dt = \int_{0}^{T} L (g^{- 1} g, g^{- 1} \overset{\cdot}{g}) dt \\ = \int_{0}^{T} l (e, ξ) dt = \int_{0}^{T} l (ξ) dt \end{matrix}

(4)

where $e$ is the unit element of Lie group, and $L$ and $l$ are the traditional Lagrangian and reduced Lagrangian, respectively. Hamilton’s principle based on Lie group and Lie algebra is

\begin{matrix} δ S = δ \int_{0}^{T} l (ξ) dt = \int_{0}^{T} 〈 \frac{\partial l}{\partial ξ}, δ ξ 〉 dt = \int_{0}^{T} 〈 \frac{\partial l}{\partial ξ}, \overset{\cdot}{η} + a d_{ξ} η 〉 dt \\ = \int_{0}^{T} 〈 - \frac{d}{dt} (\frac{\partial l}{\partial ξ}) + {ad}_{ξ}^{*} (\frac{\partial l}{\partial ξ}), η 〉 dt = 0 \end{matrix}

(5)

which is the reduced Euler–Poinaré principle, and its classic Euler–Poinaré equation reads

\frac{d}{dt} (\frac{\partial l}{\partial ξ}) - {ad}_{ξ}^{*} (\frac{\partial l}{\partial ξ}) = 0

(6)

where ${ad}_{ξ}^{*} η = [η, ξ] = η ξ - ξ η$ is the conjugate of $a d_{ξ} η$ . After $ξ$ is solved from equation (6), the corresponding Lie group $g$ can be recovered by the following reconstruction equation which is the classic Poisson equation

\overset{\cdot}{g} = g ξ

(7)

For a free rigid body, its orientation matrix can be selected as the special orthogonal group, that is, $g = R \in S (3)$ , $R^{- 1} = R^{T}$ , and $det (R) = 1$ . Its corresponding Lie algebras $ξ, η \in s (o)$ are the skew-symmetric matrices of the angular velocity and the variation in quasi rotation angles, respectively, that is, $ξ = \tilde{Ω}$ , $η = δ \tilde{Θ}$ . Thus, equations (1)–(3) can be rewritten as

\tilde{Ω} = R^{T} \overset{\cdot}{R}

(8)

δ \tilde{Θ} = R^{T} δ R

(9)

δ \tilde{Ω} = δ \overset{\cdot}{\tilde{}} Θ + a d_{\tilde{Ω}} δ \tilde{Θ}

(10)

So, we can rewrite equations (6) and (7) as follows

\frac{d}{dt} (\frac{\partial l}{\partial \tilde{Ω}}) - {ad}_{\tilde{Ω}}^{*} (\frac{\partial l}{\partial \tilde{Ω}}) = 0

(11)

\overset{\cdot}{R} = R \tilde{Ω}

(12)

The equivalent vector form of equation (10) is

δ Ω = δ \overset{\cdot}{Θ} + Ω \times δ Θ

(13)

Based on equation (13), we rewrite the Euler–Poinaré principle (5) using vector expression as follows

\begin{matrix} δ S & = δ \int_{0}^{T} l (\tilde{Ω}) dt = δ \int_{0}^{T} l (Ω) dt = \int_{0}^{T} 〈 \frac{\partial l}{\partial Ω}, δ Ω 〉 dt \\ = \int_{0}^{T} 〈 \frac{\partial l}{\partial Ω}, δ \overset{\cdot}{Θ} + Ω \times δ Θ 〉 dt \\ = \int_{0}^{T} 〈 - \frac{d}{dt} (\frac{\partial l}{\partial Ω}) + \frac{\partial l}{\partial Ω} \times Ω, δ Θ 〉 dt \end{matrix}

(14)

The corresponding Euler–Poinaré equation of a free rigid body is in the following vector form

\frac{d}{dt} (\frac{\partial l}{\partial Ω}) + Ω \times \frac{\partial l}{\partial Ω} = 0

(15)

The Lagrangian of a free rigid body is

l = \frac{1}{2} Ω \cdot J Ω

(16)

Thus, equation (15) can be transformed into the following classical Euler equation of rotation

J \overset{\cdot}{Ω} + Ω \times J Ω = 0

(17)

The Euler–Poinaré equations of rigid multibody system dynamics

A typical multibody system is depicted in Figure 1, which contains $n$ bodies and $m$ constraints. For every body $B_{i} (i = 1, 2, \dots, n)$ in the system, we set up the body-fixed frame at its mass center $C_{i}$ and its axes coincide with its three inertia principle axes. The position of a typical point $p$ on $B_{i}$ in the inertia space is

r_{i}^{p} = r_{i} + R_{i} s'_{i}

(18)

where $r_{i}$ is the position coordinate of the mass center of body $B_{i}$ in inertia space, and $R_{i}$ is the direction cosine matrix of $B_{i}$ . $s'_{i}$ is the relative position coordinate of point $p$ with respect to the mass center of $B_{i}$ expressed in the body-fixed frame. The derivative of equation (18) with respect to time $t$ is

{\overset{\cdot}{r}}_{i}^{p} = {\overset{\cdot}{r}}_{i} + {\overset{\cdot}{R}}_{i} s' = {\overset{\cdot}{r}}_{i} + R_{i} R_{i}^{T} {\overset{\cdot}{R}}_{i} s'_{i} = {\overset{\cdot}{r}}_{i} + R_{i} {\tilde{Ω}}_{i} s'_{i}

(19)

where

{\tilde{Ω}}_{i} = R_{i}^{T} {\overset{\cdot}{R}}_{i}

(20)

Figure 1.

A typical multibody system.

The kinetic energy of $B_{i}$ is

T_{i} = \frac{1}{2} \int_{V_{i}} ρ {\overset{\cdot}{r}}_{i}^{pT} {\overset{\cdot}{r}}_{i}^{p} dV = \frac{1}{2} m_{i} {\overset{\cdot}{r}}_{i}^{T} \overset{\cdot}{r} + \frac{1}{2} tr ({\overset{\cdot}{R}}_{i} {\hat{J}}_{i} {\overset{\cdot}{R}}_{i}^{T})

(21)

where ${\hat{J}}_{i}$ is the non-standard moment of inertia of $B_{i}$ , which is

{\hat{J}}_{i} = \int_{V_{i}} ρ {s'}_{i} {s'}_{i}^{T} dV = E \int_{V_{i}} {‖ {s'}_{i} ‖}^{2} dV - J_{i}

(22.1)

where $J_{i}$ is the standard moment of inertia of $B_{i}$

J_{i} = \int_{V_{i}} ρ {\tilde{s}'}_{i}^{T} {\tilde{s}'}_{i} dV

(22.2)

Equation (21) can be rewritten using angular velocity as

\begin{matrix} T_{i} = & \frac{1}{2} m_{i} ‖ {\overset{\cdot}{r}}_{i} ‖^{2} + \frac{1}{2} \int_{V_{i}} ρ {‖ {\tilde{Ω}}_{i} {s'}_{i} ‖}^{2} dV = \frac{1}{2} m_{i} ‖ {\overset{\cdot}{r}}_{i} ‖^{2} \\ + \frac{1}{2} tr ({\tilde{Ω}}_{i} {\hat{J}}_{i} {\tilde{Ω}}_{i}^{T}) \end{matrix}

(23.1)

T_{i} = \frac{1}{2} m_{i} ‖ {\overset{\cdot}{r}}_{i} ‖^{2} + \frac{1}{2} Ω_{i} \cdot J_{i} Ω_{i}

(23.2)

Here, we assume that the potential energy of the system is the following form

U = U (r, R) = U (r_{1}, r_{2}, \dots, r_{n}, R_{1}, R_{2}, \dots, R_{n})

(24)

and express the holonomic constraints of the system in the following compact form

ϕ = ϕ (r, R) = (\begin{matrix} ϕ_{1} (r_{1}, r_{2}, \dots, r_{n}, R_{1}, R_{2}, \dots, R_{n}) \\ ϕ_{2} (r_{1}, r_{2}, \dots, r_{n}, R_{1}, R_{2},, R_{n}) \\ \dots \\ ϕ_{m} (r_{1}, r_{2}, \dots, r_{n}, R_{1}, R_{2}, \dots, R_{n}) \end{matrix}) = 0

(25)

where $r = {[\begin{matrix} r_{1}^{T} & r_{2}^{T} & \dots & r_{n}^{T} \end{matrix}]}^{T}$ and $R = (R_{1}, R_{2}, \dots, R_{n})$ .

Next, we separately denote the translation velocities, angular velocities, and their skew-symmetric matrices as $v = \overset{\cdot}{r}$ , $v = {[\begin{matrix} v_{1}^{T} & v_{2}^{T} & \dots & v_{n}^{T} \end{matrix}]}^{T}$ , $Ω = {[\begin{matrix} Ω_{1}^{T} & Ω_{2}^{T} & \dots & Ω_{n}^{T} \end{matrix}]}^{T}$ , and $\tilde{Ω} = dig [\begin{matrix} {\tilde{Ω}}_{1} & {\tilde{Ω}}_{2} & \dots & {\tilde{Ω}}_{n} \end{matrix}]$ . Thus, the reduced Lagrangian of the system is

\begin{matrix} L & = T - U = L (r, \overset{\cdot}{r}, R, \overset{\cdot}{R}) = l (r, v, R, \tilde{Ω}) \\ = T (r, v, \tilde{Ω}) - U (r, R) \end{matrix}

(26)

where $T = \sum_{i = 1}^{n} T_{i}$ is the kinetic energy. Based on these notations, the Hamilton principle can be stated as the following reduced form

\begin{matrix} δ S & = δ \int_{0}^{T} (L (r, v, R, \overset{\cdot}{R}) - λ^{T} ϕ (r, R)) dt \\ = δ \int_{0}^{T} (l (r, v, R, \tilde{Ω}) - λ^{T} ϕ (r, R)) dt = 0 \end{matrix}

(27)

According to equation (10), equation (27) can be manipulated further as

\begin{matrix} δ S = \int_{0}^{T} (〈 \frac{\partial l}{\partial r}, δ r 〉 + 〈 \frac{\partial l}{\partial v}, δ \overset{\cdot}{r} 〉 + 〈 λ^{T} \frac{\partial ϕ}{\partial r}, δ r 〉) dt \\ + \sum_{i = 1}^{n} \int_{0}^{T} (〈 \frac{\partial l}{\partial R_{i}}, δ R_{i} 〉 + 〈 - \frac{d}{dt} (\frac{\partial l}{\partial {\tilde{Ω}}_{i}}) + {ad}_{{\tilde{Ω}}_{i}}^{*} (\frac{\partial l}{\partial {\tilde{Ω}}_{i}}), δ {\tilde{Θ}}_{i} 〉) dt \\ - \int_{0}^{T} (〈 Φ, δ λ 〉 + 〈 \frac{\partial (λ^{T} ϕ (r, R))}{\partial R_{i}}, δ R_{i} 〉) dt = 0 \end{matrix}

(28)

For any two matrices $A_{n \times n}$ and $B_{n \times n}$ , the following properties are applicable

〈 A, B 〉 : = \sum_{p, q = 1}^{n} A_{p, q} B_{p, q} = tr (A B^{T})

(29.1)

A = \frac{1}{2} (A + A^{T}) + \frac{1}{2} (A - A^{T}) = A_{Sym} + A_{ssy}

(29.2)

tr (AB) = tr (BA) = tr (B^{T} A^{T}) = tr (A^{T} B^{T})

(29.3)

where $A_{Sym} = 1 / 2 (A + A^{T})$ and $A_{ssy} = 1 / 2 (A - A^{T})$ are, respectively, symmetric matrix and skew-symmetric matrix. So, we have

\begin{matrix} 〈 \frac{\partial l}{\partial R_{i}}, δ R_{i} 〉 & = 〈 \frac{\partial l}{\partial R_{i}}, R_{i} R_{i}^{T} δ R_{i} 〉 = 〈 \frac{\partial l}{\partial R_{i}}, R_{i} δ {\tilde{Θ}}_{i} 〉 \\ = tr (\frac{\partial l}{\partial R_{i}} {(R_{i} δ {\tilde{Θ}}_{i})}^{T}) = tr (\frac{\partial l}{\partial R_{i}} δ {\tilde{Θ}}_{i}^{T} R_{i}^{T}) \\ = tr (R_{i}^{T} \frac{\partial l}{\partial R_{i}} δ {\tilde{Θ}}_{i}^{T}) = 〈 R_{i}^{T} \frac{\partial l}{\partial R_{i}}, δ {\tilde{Θ}}_{i} 〉 \\ = 〈 \frac{1}{2} (R_{i}^{T} \frac{\partial l}{\partial R_{i}} + {(\frac{\partial l}{\partial R_{i}})}^{T} R_{i}) \\ + \frac{1}{2} (R_{i}^{T} \frac{\partial l}{\partial R_{i}} - {(\frac{\partial l}{\partial R_{i}})}^{T} R_{i}), δ {\tilde{Θ}}_{i} 〉 \end{matrix}

(30)

As $1 / 2 (R_{i}^{T} (\partial l / \partial R_{i}) + {(\partial l / \partial R_{i})}^{T} R_{i})$ is a symmetric matrix and $δ {\tilde{Θ}}_{i}$ is a skew-symmetric matrix, we have

〈 \frac{1}{2} (R_{i}^{T} \frac{\partial l}{\partial R_{i}} + {(\frac{\partial l}{\partial R_{i}})}^{T} R_{i}), δ {\tilde{Θ}}_{i} 〉 = 0

(31.1)

Because $1 / 2 (R_{i}^{T} (\partial l / \partial R_{i}) + {(\partial l / \partial R_{i})}^{T} R_{i})$ is a skew-symmetric matrix, we can define

{\overset{⌣}{Ξ}}_{i} = \frac{1}{2} (R_{i}^{T} \frac{\partial l}{\partial R_{i}} - {(\frac{\partial l}{\partial R_{i}})}^{T} R_{i}) = - \frac{1}{2} {\tilde{Ξ}}_{i}

(31.2)

For any vectors $x, y$ , $- 1 / 2 〈 \tilde{x}, \tilde{y} 〉 = 〈 x, y 〉 = x^{T} y$ , equation (30) can be written as

〈 \frac{\partial l}{\partial R_{i}}, δ R_{i} 〉 = 〈 {\overset{⌣}{Ξ}}_{i}, δ {\tilde{Θ}}_{i} 〉 = - \frac{1}{2} 〈 {\tilde{Ξ}}_{i}, δ {\tilde{Θ}}_{i} 〉 = 〈 Ξ_{i}, δ Θ_{i} 〉

(31.3)

where

\frac{\partial l}{\partial R_{i}} = \frac{\partial U}{\partial R_{i}}

(31.4)

Based on the identity $\tilde{(x \times y)} = \tilde{x} \tilde{y} - \tilde{y} \tilde{x} = y x^{T} - x y^{T}$ , we can give the physical explanation of ${\tilde{Ξ}}_{i}$ and $Ξ_{i}$ as follows

\begin{matrix} R_{i}^{T} \frac{\partial U}{\partial R_{i}} & - {(\frac{\partial U}{\partial R_{i}})}^{T} R_{i} = [\begin{matrix} R_{i 1}^{T} & R_{i 2}^{T} & R_{i 3}^{T} \end{matrix}] [\begin{matrix} U_{R_{i 1}} \\ U_{R_{i 2}} \\ U_{R_{i 3}} \end{matrix}] \\ - [\begin{matrix} U_{R_{i 1}}^{T} & U_{R_{i 2}}^{T} & U_{R_{i 3}}^{T} \end{matrix}] [\begin{matrix} R_{i 1} \\ R_{i 2} \\ R_{i 3} \end{matrix}] \\ = (R_{i 1}^{T} U_{R_{i 1}} - U_{R_{i 1}}^{T} R_{i 1}) + (R_{i 2}^{T} U_{R_{i 2}} - U_{R_{i 2}}^{T} R_{i 2}) \\ + (R_{i 3}^{T} U_{R_{i 3}} - U_{R_{i 3}}^{T} R_{i 3}) = \tilde{- (R_{i 1}^{T} \times U_{R_{i 1}}^{T})} \\ - \tilde{(R_{i 2}^{T} \times U_{R_{i 2}}^{T})} - \tilde{(R_{i 3}^{T} \times U_{R_{i 3}}^{T})} \end{matrix}

(31.5)

So, we have

- {\tilde{Ξ}}_{i} = - \tilde{(R_{i 1}^{T} \times U_{R_{i 1}}^{T})} - \tilde{(R_{i 2}^{T} \times U_{R_{i 2}}^{T})} - \tilde{(R_{i 3}^{T} \times U_{R_{i 3}}^{T})}

(31.6)

- Ξ_{i} = - (R_{i 1}^{T} \times U_{R_{i 1}}^{T}) - (R_{i 2}^{T} \times U_{R_{i 2}}^{T}) - (R_{i 3}^{T} \times U_{R_{i 3}}^{T})

(31.7)

Compared with equation (31.3), we can obtain

Ξ_{i} = \frac{\partial U}{\partial Θ_{i}}

(31.8)

$Ξ_{i}$ is the partial derivative of potential energy with respect to quasi-angles, while ${\tilde{Ξ}}_{i}$ is its skew-symmetric matrix. Using similar method, we can derive the variation in constraint terms in equation (28) with respect to orientation matrix as follows

\begin{matrix} 〈 \frac{\partial (λ^{T} ϕ (r, R))}{\partial R_{i}}, δ R_{i} 〉 = \frac{1}{2} 〈 (R_{i}^{T} \frac{\partial (λ^{T} ϕ (r, R))}{\partial R_{i}} - {(\frac{\partial (λ^{T} ϕ (r, R))}{\partial R_{i}})}^{T} R_{i}), δ {\tilde{Θ}}_{i} 〉 \\ = 〈 \frac{1}{2} (R_{i}^{T} \frac{\partial (\sum_{j = 1}^{m} λ_{j} ϕ_{j} (r, R))}{\partial R_{i}} - {(\frac{\partial (\sum_{j = 1}^{m} λ_{j} ϕ_{j} (r, R))}{\partial R_{i}})}^{T} R_{i}), δ {\tilde{Θ}}_{i} 〉 \\ = \sum_{j = 1}^{m} 〈 λ_{j}, 〈 \frac{1}{2} (R_{i}^{T} \frac{\partial ϕ_{j} (r, R)}{\partial R_{i}} - {(\frac{\partial ϕ_{j} (r, R)}{\partial R_{i}})}^{T} R_{i}), δ {\tilde{Θ}}_{i} 〉 〉 \\ = \sum_{j = 1}^{m} 〈 λ_{j}, 〈 {\overset{⌣}{Φ}}_{j, i}, δ {\tilde{Θ}}_{i} 〉 〉 = - \sum_{j = 1}^{m} 〈 λ_{j}, 〈 \frac{1}{2} {\tilde{Φ}}_{j, i}, δ {\tilde{Θ}}_{i} 〉 〉 = \sum_{j = 1}^{m} 〈 λ_{j}, 〈 Φ_{j, i}, δ Θ_{i} 〉 〉 \\ = \sum_{j = 1}^{m} 〈 λ_{j}, 〈 ϕ_{j, Θ_{i}}, δ Θ_{i} 〉 〉 = λ^{T} ϕ_{Θ} δ Θ \end{matrix}

(31.9)

where

Φ_{j, i} = ϕ_{j, Θ_{i}} = \frac{\partial ϕ_{j}}{\partial Θ_{i}}

(31.9)

and its corresponding skew-symmetric matrix is

- {\tilde{Φ}}_{j, i} = R_{i}^{T} \frac{\partial ϕ_{j} (r, R)}{\partial R_{i}} - {(\frac{\partial ϕ_{j} (r, R)}{\partial R_{i}})}^{T} R_{i}

(31.9)

Substituting equations (31.3) and (31.9) into equation (28), we obtain

\begin{array}{l} δ S = \int_{0}^{T} (〈 \frac{\partial l}{\partial r}, δ r 〉 + 〈 - \frac{d}{d t} (\frac{\partial l}{\partial v}), δ r 〉 - 〈 λ^{T} \frac{\partial ϕ}{\partial r}, δ r 〉) d t - \int_{0}^{T} 〈 Φ, δ λ 〉 d t \\ + \sum_{i = 1}^{n} \int_{0}^{T} (〈 {\overset{⌣}{Ξ}}_{i}, δ {\tilde{Θ}}_{i} 〉 + 〈 - \frac{d}{d t} (\frac{\partial l}{\partial {\tilde{Ω}}_{i}}) + a d_{{\tilde{Ω}}_{i}}^{*} (\frac{\partial l}{\partial {\tilde{Ω}}_{i}}), δ {\tilde{Θ}}_{i} 〉 - 〈 λ^{T} {\overset{⌣}{Φ}}_{i}, δ {\tilde{Θ}}_{i} 〉) d t = 0 \end{array}

(32)

where the dynamical equations of multibody systems can be derived as

{\begin{matrix} \frac{d}{d t} (\frac{\partial l}{\partial v}) - \frac{\partial l}{\partial r} + ϕ_{r}^{T} λ = 0 \\ \frac{d}{d t} (\frac{\partial l}{\partial {\tilde{Ω}}_{i}}) - a d_{{\tilde{Ω}}_{i}}^{*} (\frac{\partial l}{\partial {\tilde{Ω}}_{i}}) - {\overset{⌣}{Ξ}}_{i} + λ^{T} {\overset{⌣}{Φ}}_{i} = 0 (i = 1, 2, \dots, n) \\ ϕ = 0 \\ \dot{r} = v \\ {\dot{R}}_{i} = R_{i} {\tilde{Ω}}_{i} (i = 1, 2, \dots, n) \end{matrix}

(33)

where the first equation is the traditional Euler–Lagrange equation on translational motion of mass center. The second one is the Euler–Poinaré equation on rotational motion of every body. The third one is constraint equation. Note that the last two equations are reconstruction equations.

By using equation (13), we can rewrite the first two equations of equation (33) as the following vector form

m \overset{\cdot}{v} - \frac{\partial U}{\partial r} + ϕ_{r}^{T} λ = 0

(34.1)

\frac{d}{dt} (\frac{\partial l}{\partial Ω}) + Ω \times (\frac{\partial l}{\partial Ω}) - U_{Θ} + ϕ_{Θ}^{T} λ = 0

(34.2)

For rigid bodies, we have

\frac{\partial l}{\partial Ω} = J Ω

(35)

With equation (35), equation (34.2) can be rewritten as the following simple form

J \overset{\cdot}{Ω} + Ω \times J Ω - U_{Θ} + ϕ_{Θ}^{T} λ = 0

(36)

The constraints for the interconnected bodies are important features of multimode systems. In this section, we employ the method proposed in Haug^5,6 to construct constraint library for the multibody systems according to the definition of basic constraints. The constraint equations can be expressed as the following general form

ϕ (r_{i}, R_{i}, r_{j}, R_{j}) = 0

(37.1)

Its variations, first-order derivatives, and second order-derivatives are, respectively

δ ϕ = ϕ_{r_{i}} δ r_{i} + ϕ_{Θ_{i}} δ Θ_{i} + ϕ_{Θ_{j}} δ Θ_{j} + ϕ_{r_{j}} δ r_{j}

(37.2)

\overset{\cdot}{ϕ} = ϕ_{r_{i}} {\overset{\cdot}{r}}_{i} + ϕ_{Θ_{i}} Ω_{i} + ϕ_{Θ_{j}} Ω_{j} + ϕ_{r_{j}} {\overset{\cdot}{r}}_{j} = 0

(37.3)

\begin{matrix} \overset{\cdot\cdot}{ϕ} = & ϕ_{r_{i}} {\overset{\cdot\cdot}{r}}_{i} + ϕ_{Θ_{i}} {\overset{\cdot}{Ω}}_{i} + ϕ_{Θ_{j}} {\overset{\cdot}{Ω}}_{j} + ϕ_{r_{j}} {\overset{\cdot\cdot}{r}}_{j} + {\overset{\cdot}{ϕ}}_{r_{i}} {\overset{\cdot}{r}}_{i} \\ + {\overset{\cdot}{ϕ}}_{Θ_{i}} Ω_{i} + {\overset{\cdot}{ϕ}}_{Θ_{j}} Ω_{j} + {\overset{\cdot}{ϕ}}_{r_{j}} {\overset{\cdot}{r}}_{j} = 0 \end{matrix}

(37.4)

For the velocity level constraint equation (37.3), we present its derivation in detail as follows

\begin{matrix} \overset{\cdot}{ϕ} & = \overset{\cdot}{ϕ} (r_{i}, R_{i}, r_{j}, R_{j}) \\ = ϕ_{r_{i}} {\overset{\cdot}{r}}_{i} + ϕ_{r_{j}} {\overset{\cdot}{r}}_{j} + 〈 \frac{\partial ϕ}{\partial R_{i}}, {\overset{\cdot}{R}}_{i} 〉 + 〈 \frac{\partial ϕ}{\partial R_{j}}, {\overset{\cdot}{R}}_{j} 〉 \\ = ϕ_{r_{i}} {\overset{\cdot}{r}}_{i} + ϕ_{r_{j}} {\overset{\cdot}{r}}_{j} + 〈 \frac{\partial ϕ}{\partial R_{i}}, R_{i} {\tilde{Ω}}_{i} 〉 + 〈 \frac{\partial ϕ}{\partial R_{j}}, R_{j} {\tilde{Ω}}_{j} 〉 \\ = ϕ_{r_{i}} {\overset{\cdot}{r}}_{i} + ϕ_{r_{j}} {\overset{\cdot}{r}}_{j} + 〈 \frac{\partial ϕ}{\partial R_{i}}, R_{i} {\tilde{Ω}}_{i} 〉 + 〈 \frac{\partial ϕ}{\partial R_{j}}, R_{j} {\tilde{Ω}}_{j} 〉 \\ = ϕ_{r_{i}} {\overset{\cdot}{r}}_{i} + ϕ_{r_{j}} {\overset{\cdot}{r}}_{j} - tr ({\tilde{Ω}}_{i} R_{i}^{T} \frac{\partial ϕ}{\partial R_{i}}) - tr ({\tilde{Ω}}_{j} R_{j}^{T} \frac{\partial ϕ}{\partial R_{j}}) \\ = ϕ_{r_{i}} {\overset{\cdot}{r}}_{i} + ϕ_{r_{j}} {\overset{\cdot}{r}}_{j} - tr (R_{i}^{T} \frac{\partial ϕ}{\partial R_{i}} {\tilde{Ω}}_{i}) - tr (R_{j}^{T} \frac{\partial ϕ}{\partial R_{j}} {\tilde{Ω}}_{j}) \\ = ϕ_{r_{i}} {\overset{\cdot}{r}}_{i} + ϕ_{r_{j}} {\overset{\cdot}{r}}_{j} + 〈 R_{i}^{T} \frac{\partial ϕ}{\partial R_{i}}, {\tilde{Ω}}_{i} 〉 + 〈 R_{j}^{T} \frac{\partial ϕ}{\partial R_{j}}, {\tilde{Ω}}_{j} 〉 \\ = ϕ_{r_{i}} {\overset{\cdot}{r}}_{i} + ϕ_{r_{j}} {\overset{\cdot}{r}}_{j} + ϕ_{Θ_{i}} Ω_{i} + ϕ_{Θ_{j}} Ω_{j} \end{matrix}

(37.5)

Generalized external forces in Lie groups and Lie algebras

In this section, we focus on some conservative external forces, such as gravity, elastic force due to linear translational and twisted springs. The corresponding generalized forces are minus derivatives of energy with respect to quasi-coordinates.

For the generalized external forces due to gravity, we let $\bar{g} = {[\begin{matrix} 0 & 0 & g \end{matrix}]}^{T}$ denote gravity acceleration matrix. Based on the assumption that the origin of body-fixed frame is located at mass center and its axes are along the three inertia principle axes, we have the following three equations

U_{i} = r_{i}^{T} m_{i} \bar{g}

(38.1)

δ U_{i} = δ r_{i}^{T} m_{i} \bar{g}

(38.2)

\frac{\partial U_{i}}{\partial r_{i}} = m_{i} \bar{g}, \frac{\partial U_{i}}{\partial Θ_{i}} = 0

(38.3)

Equations (38.1), (38.2), and (38.3) are potential energy of every body, its variation, and partial derivative with quasi-coordinates, respectively.

By taking account into the linear translational spring between point $p$ on $B_{i}$ and point $q$ on $B_{j}$ , we denote $d_{ij}$ and $k$ as the relative position vector between points $p$ and $q$ and the elastic modulus, respectively. The elastic potential energy and its variation are given as follows

U_{ij} = \frac{1}{2} k {| | d_{ij} | - | d_{ij 0} | |}^{2}

(39.1)

δ U_{ij} = F_{ij} \frac{d_{ij}^{T}}{| d_{ij} |} δ d_{ij}

(39.2)

where the variation in $d_{ij}$ is

\begin{matrix} δ d_{ij} & = δ r_{j} + δ R_{j} s'_{j}^{q} - δ r_{i} - δ R_{i} s'_{i}^{p} \\ = δ r_{j} + R_{j} δ {\tilde{Θ}}_{j} s'_{j}^{q} - δ r_{i} - R_{i} δ {\tilde{Θ}}_{i} s'_{i}^{p} \\ = δ r_{j} - R_{j} \tilde{s}'_{j}^{q} δ Θ_{j} - δ r_{i} + R_{i} \tilde{s}'_{i}^{p} δ Θ_{i} \end{matrix}

(39.3)

According to equation (39.3), we have $F_{ij} = k | | d_{ij} | - | d_{ij 0} | |$ in equation (39.2).

From equation (39.2), we can derive the derivatives of energy with respect to quasi-coordinates as follows

\begin{matrix} \frac{\partial U_{ij}}{\partial r_{j}} = F_{ij} \frac{d_{ij}^{T}}{| d_{ij} |}, \frac{\partial U_{ij}}{\partial r_{i}} = - F_{ij} \frac{d_{ij}^{T}}{| d_{ij} |}, \\ \frac{\partial U_{ij}}{\partial Θ_{j}} = - F_{ij} \frac{d_{ij}^{T}}{| d_{ij} |} R_{j} \tilde{s}'_{j}^{q}, \frac{\partial U_{ij}}{\partial Θ_{i}} = F_{ij} \frac{d_{ij}^{T}}{| d_{ij} |} R_{i} \tilde{s}'_{i}^{p} \end{matrix}

(39.4)

For a twisted spring between $B_{i}$ and $B_{j}$ , the elastic energy and its variation are

U_{ij} = \frac{1}{2} k {| | Θ_{ij} | - | Θ_{ij 0} | |}^{2}

(40.1)

δ U_{ij} = M_{ij} \frac{Θ_{ij}^{T}}{| Θ_{ij} |} δ Θ_{ij}

(40.2)

where $M_{ij} = k | | Θ_{ij} | - | Θ_{ij 0} | |$ . According to the following equations

R_{j} = R_{i} R_{ij}, R_{ij} = R_{i}^{T} R_{j}

(40.3)

δ R_{ij} = δ R_{i}^{T} R_{j} + R_{i}^{T} δ R_{j}

(40.4)

R_{ij} R_{ij}^{T} δ R_{ij} = - R_{i}^{T} δ R_{i} R_{i}^{T} R_{j} + R_{i}^{T} R_{j} R_{j}^{T} δ R_{j}

(40.5)

R_{ij} δ {\tilde{Θ}}_{ij} = - δ {\tilde{Θ}}_{i} R_{ij} + R_{ij} δ {\tilde{Θ}}_{j}

(40.6)

One can obtain

δ {\tilde{Θ}}_{ij} = δ {\tilde{Θ}}_{j} - R_{ji} δ {\tilde{Θ}}_{i} R_{ij} = δ {\tilde{Θ}}_{j} - A d_{R_{ji}} δ {\tilde{Θ}}_{i}

(40.7)

δ Θ_{ij} = δ Θ_{j} - R_{ji} δ Θ_{i}

(40.8)

Finally, the partial derivatives of energy with respect to quasi-coordinates are written as follows

\begin{matrix} \frac{\partial U_{ij}}{\partial r_{i}} = 0, \frac{\partial U_{ij}}{\partial r_{j}} = 0, \frac{\partial U_{ij}}{\partial Θ_{j}} = M_{ij} \frac{Θ_{ij}^{T}}{| Θ_{ij} |}, \\ \frac{\partial U_{ij}}{\partial Θ_{j}} = - M_{ij} \frac{Θ_{ij}^{T}}{| Θ_{ij} |} R_{ji} \end{matrix}

(40.9)

Examples

A heavy top

For a heavy top presented in Figure 2, we set up its body-fixed frame according to section “The Euler–Poinaré equations of rigid multibody system dynamics.” Its kinetic energy and potential energy are, respectively

T = \frac{1}{2} m v^{T} v + \frac{1}{2} Ω^{T} J Ω

U = r^{T} m \bar{g}

Figure 2.

Heavy top.

The constraint equation at joint O and its variation are, respectively

ϕ (r, R) = r + Rs' = 0

\begin{matrix} δ ϕ (r, R) & = δ r + δ Rs' = δ r + R R^{T} δ Rs' \\ = δ r + R δ \tilde{Θ} s' = δ r - R \tilde{s}' δ Θ \end{matrix}

where $s'$ is the coordinate of $O$ in the body-fixed frame. The dynamical equations of the heavy top can be derived as follows according to the results in section “The Euler–Poinaré equations of rigid multibody system dynamics”

{\begin{matrix} m \overset{\cdot}{v} - m \bar{g} + λ = 0 \\ J \overset{\cdot}{Ω} + Ω \times J Ω + Rs' \times λ = 0 \\ r + Rs' = 0 \\ \overset{\cdot}{r} = v \\ \overset{\cdot}{R} = R \tilde{Ω} \end{matrix}

A spatial pendulum with three rigid bodies

A spatial pendulum shown in Figure 3 consists of three bodies interconnected through three spherical joints. The kinetic energy, potential energy, constraint equations, and their variations are, respectively

T = \frac{1}{2} \sum_{i = 1}^{3} m_{i} v_{i}^{T} v_{i} + \frac{1}{2} \sum_{i = 1}^{3} Ω_{i}^{T} J_{i} Ω_{i}

U = \sum_{i = 1}^{3} r_{i}^{T} m_{i} \bar{g}

ϕ (r, R) = [\begin{matrix} ϕ_{1} (r, R) \\ ϕ_{2} (r, R) \\ ϕ_{3} (r, R) \end{matrix}] = [\begin{matrix} r_{1} + R_{1} s'_{1}^{1} \\ r_{1} + R_{1} s'_{1}^{2} - r_{2} - R_{2} s'_{2}^{2} \\ r_{2} + R_{2} s'_{2}^{3} - r_{3} - R_{3} s'_{3}^{3} \end{matrix}] = 0

\begin{matrix} δ ϕ = & [\begin{matrix} 1 & 0 & 0 \\ 1 & - 1 & 0 \\ 0 & 1 & - 1 \end{matrix}] [\begin{matrix} δ r_{1} \\ δ r_{2} \\ δ r_{3} \end{matrix}] \\ + [\begin{matrix} - R_{1} \tilde{s}'_{1}^{1} & 0 & 0 \\ - R_{1} \tilde{s}'_{1}^{2} & R_{2} \tilde{s}'_{2}^{2} & 0 \\ 0 & - R_{2} \tilde{s}'_{2}^{3} & R_{3} \tilde{s}'_{3}^{3} \end{matrix}] [\begin{matrix} δ Θ_{1} \\ δ Θ_{2} \\ δ Θ_{3} \end{matrix}] \\ = ϕ_{r} δ r + ϕ_{Θ} δ Θ = 0 \end{matrix}

Figure 3.

A spatial pendulum.

The dynamic equations of the system are derived as follows according to the results in section “The Euler–Poinaré equations of rigid multibody system dynamics”

{\begin{matrix} m \overset{\cdot}{v} - m \bar{g} + ϕ_{r}^{T} λ = 0 \\ J \overset{\cdot}{Ω} + Ω \times J Ω + ϕ_{Θ}^{T} λ = 0 \\ ϕ (r, R) = 0 \\ {\overset{\cdot}{r}}_{i} = v_{i} \\ {\overset{\cdot}{R}}_{i} = R_{i} {\tilde{Ω}}_{i} \\ (i = 1, 2, 3) \end{matrix}

Obviously, numerical methods of differential/algebraic equations in Lie group appearing in this article will be a new research subject in multibody system dynamics, which are not discussed here due to scope of this research. Please refer to Uhlar and Betsch,¹⁰ Brüls et al.,¹¹ Brüls and Cardona,¹² Celledoni et al.,¹³ Celledoni and Owren,¹⁴ and Müller and Terze³⁴ for some new developments on this topic.

Conclusion

A new type of dynamical equations of rigid multibody systems is derived using Euler–Poinaré principle. The results include four typical equations: the Euler–Lagrange equations with constraints on translations of mass centers; the Euler–Poinaré equations with constraints on rotations of every body; the reconstruction equations used to recover the orientation matrices from Lie algebras; and the constraint equations describing inter-connections of bodies. All these equations have simple forms which are easy for numerical implementation. The inertia matrices are constant, leading to high computational efficiency. In addition, the Lie group structure-preserving algorithms can be developed from the Lie group integrators, which is very important for long-term simulations of multibody systems.

This work can be extended to other new researches naturally, including the following: (1) the discrete Euler–Poinaré equations and Lie–Poisson equations of multibody systems based on the variational principles for discrete mechanics; (2) continuous and discrete Euler–Poinaré equations and Lie–Poisson equations of flexible multibody systems; and (3) Lie group variational integrators of multibody systems and their applications in the field of control and optimization and so forth.

Footnotes

Academic Editor: Chuanzeng Zhang

Declaration of conflicting interests

The authors declare that there is no conflict of interest.

Funding

This research was supported by the Natural Science Foundation of China (Grant Nos 11472144, 11002075, 11272166, 11472143).

References

Liu

Pan

. Dynamics of multibody systems. Beijing, China: Higher Education Press, 2014.

Hong

. Computational multibody system dynamics. Beijing, China: Higher Education Press, 1999.

Shabana

. Dynamics of multibody systems. New York: Wiley, 1989.

Nikravesh

. Computer-aided analysis of mechanical systems. Englewood Cliffs, NJ: Prentice Hall, 1988.

Haug

. Computer aided kinematics and dynamics of mechanical systems, volume 1: basic methods. Needham Heights, MA: Allyn & Bacon, 1989.

Haug

. Intermediate dynamics. Englewood Cliffs, NJ: Prentice Hall, 1992.

Géradin

Cardona

. Flexible multibody dynamics: a finite element approach. Chichester: John Wiley & Sons, 2001.

Bauchau

. Flexible multibody dynamics. Berlin: Springer, 2011.

Simeon

. Computational flexible multibody dynamics—a differential algebraic approach. Berlin: Springer, 2013.

10.

Uhlar

Betsch

. A rotationless formulation of multibody dynamics: modeling of screw joints and incorporation of control constraints. Multibody Syst Dyn 2009; 22(1): 69–95.

11.

Brüls

Cardona

Arnold

. Lie group generalized-α time integration of constrained flexible multibody systems. Mech Mach Theor 2012; 48: 121–137.

12.

Brüls

Cardona

. On the use of Lie group time integrators in multibody dynamics. J Comput Nonlin Dyn 2010; 5: 1–13.

13.

Celledoni

Marthinsen

Owren

. An introduction to Lie group integrators-basics, new developments and applications. J Comput Phys 2014; 257B: 1040–1061.

14.

Celledoni

Owren

. Lie group methods for rigid body dynamics and time integration on manifolds. Comput Meth Appl Mech Eng 2003; 192: 421–438.

15.

Holm

Marsden

Ratiu

. The Euler-Poinaré equations and semidirect products with applications to continuum theories. Adv Math 1998; 137: 1–81.

16.

Holm

Schmah

Stoica

. Geometric mechanics and symmetry: from finite to infinite dimensions. Oxford: Oxford University Press, 2009.

17.

Rowley

. Lie-Poisson integrators: a Hamiltonian, variational approach. Int J Numer Meth Eng 2010; 82: 1609–1644.

18.

Crouch

Grossman

. Numerical integration of ordinary differential equations on manifolds. J Nonlinear Sci 1993; 3: 1–33.

19.

Lewis

Simo

. Conserving algorithms for the dynamics of Hamiltonian systems on Lie groups. J Nonlinear Sci 1994; 4: 253–299.

20.

Munthe-Kaas

. High order Runge-Kutta methods on manifolds. Appl Numer Math 1999; 29: 115–127.

21.

Iserles

Munthe-Kaas

Nørsett

. Lie-group methods. Acta Numer 2000; 9: 215–365.

22.

Hairer

Lubich

Wanner

. Geometric numerical integration: structure-preserving algorithms for ordinary differential equations. 2nd ed. Berlin: Springer-Verlag, 2006.

23.

Marsden

Pekarsky

Shkoller

. Discrete Euler-Poinaré and Lie-Poisson equations. Nonlinearity 1999; 12: 1647–1662.

24.

Marsden

West

. Discrete mechanics and variational integrators. Acta Numer 2001; 10: 357–514.

25.

Lee

Leok

McClamroch

. Lie group variational integrators for the full body problem. Comput Meth Appl Mech Eng 2007; 196: 2907–2924.

26.

Bou-Rabee

Marsden

. Hamilton-Pontryagin integrators on Lie groups. Part I: introduction and structure preserving properties. Found Comput Math 2009; 9: 197–219.

27.

Mizzi

. Equations of motion of a multibody system and Lie groups. Meccanica 1995; 30: 201–219.

28.

Liu

. A Lie group formulation of Kane’s equations for multibody systems. Multibody Syst Dyn 2008; 20: 29–49.

29.

Boyer

Primault

. The Poincare-Chetayev equations and flexible multibody systems. J Appl Math Mech 2005; 69: 925–942.

30.

Park

Chung

W-K

. Geometric integration on Euclidean group with application to articulated multibody systems. IEEE Trans Robot 2005; 21: 850–863.

31.

Brüls

Arnold

Cardona

. Two Lie group formulations for dynamic multibody systems with large rotations. In: Proceedings of the ASME 2011, international design engineering technical conferences and computers and information in engineering conference, IDETC/CIE 2011, Washington, DC, 28–31 August 2011, DETC2011/MSNDC-48132.

32.

Sonneville

Brüls

. Formulation of kinematic joints and rigidity constraints in multibody dynamics using a Lie group approach. In: The 2nd joint international conference on multibody system dynamics, Stuttgart, 29 May–1 June 2012.

33.

Sonneville

Brüls

. A formulation on the special Euclidean group for dynamic analysis of multibody systems. J Comput Nonlin Dyn 2014; 9: 041002 (8 pp.).

34.

Müller

Terze

. The significance of the configuration space Lie group for the satisfaction in numerical time integration of multibody systems. Mech Mach Theory 2014; 82: 173–202.

Dynamical equations of multibody systems on Lie groups

Abstract

Keywords

Introduction

The classical Euler–Poinaré equations of a rigid body

The Euler–Poinaré equations of rigid multibody system dynamics

Generalized external forces in Lie groups and Lie algebras

Examples

A heavy top

A spatial pendulum with three rigid bodies

Conclusion

Footnotes

Declaration of conflicting interests

Funding

References