The Role of Mass Distribution, Higher – Order Acceleration Energies and Generalized Forces in the Advanced Dynamics

1. Introduction

In recent years, the development of dynamic models for mechanical systems has evolved toward more accurate and universally adaptable formulations. This shift reflects the increasing complexity of engineering applications, particularly in fields like robotics, aerospace, and flexible multibody systems. In such systems, traditional simplifications – such as modeling discrete masses and neglecting higher – order acceleration – fail to capture critical effects inherent to complex behaviors. Within this context, three central aspects have emerged as essential for accurate modeling: mass distribution, higher-order acceleration energies, and generalized forces.

Mass distribution significantly affects system dynamics, especially flexible structures or spatial linkages. Unlike models based solely on discrete parameters, mass distribution models allow for spatial variation of inertia and stiffness. This capability leads to improved accuracy in analyzing the local dynamics, internal wave propagation, and modal coupling effects.

Mizhidon [1] introduced a unified mathematical framework that incorporates both discrete and distributed parameters yielding a more realistic representation of mechanical systems. His model enhances predictions of deformation and energy transfer across structural elements with variable density and elasticity.

Similarly, Yang et al. [2] proposed a multiscale approach for flexible multibody systems, integrating spatially distributed mass and stiffness parameters. Their findings demonstrate that neglecting distributed effects introduces systematic errors in dynamic simulations, particularly in high-frequency or transient regimes.

In classical dynamics, kinetic energy is typically defined based on velocities, i.e., first-order derivatives. However, in systems subjected to sudden accelerations, shocks, or control input variations, higher-order derivatives such as jerk or snap become relevant. Negrean, Crișan, and Vlase [3,4] have introduced theoretical formulations that incorporate these acceleration energies within the analytical mechanics framework. Specifically, they extended the classical Lagrangian and Gibbs-Appell equations to include second- and third-order time derivatives of the acceleration energies.

In a subsequent study [5], the same authors extend the formulation on kinetic and acceleration energy expressions and propose new formulations for generalized variables based on higher-order derivatives. The concept of generalized forces is proposed as an extension of d’Alembert Principle, adapted for systems characterized by higher order accelerations. Similarly, de Jong et al. [6] derived dynamic equations for rigid body spatial mechanisms incorporating higher-order derivatives, demonstrating that dynamic balancing – often addressed heuristically – can be treated rigorously handled using analytical models based on jerk and snap. Generalized forces thus became central in advanced dynamic. These are not limited to some functions of displacement and velocity but can involve accelerations and their derivatives. Jafari et al. [7] presented a flexible-link robot model that integrates generalized forces derived from both structural compliance and control feedback mechanisms. Their model improves stability and accuracy, particularly under high-speed motion or external influences.

Complementarily, Condurache [8] conducted a theoretical study on the effects of higher-order accelerations in the dynamic behavior of mechanical linkages. He demonstrated that forces transmitted through joints and links must be defined with respect to time derivatives of acceleration, especially in high-speed motion regimes.

This perspective creates opportunities for developing different control strategies and optimization algorithms based on higher-order acceleration energies.

In the present paper, the simultaneous evaluation of mass distribution, higher-order acceleration energies, and generalized forces lead to a more comprehensive and reliable modeling approach. The work aims to synthesize recent theoretical advancements and develop a unified dynamic model applicable to mechanical systems with distributed mass and inertia, while preserving rigid body assumptions.

2. Mass Distribution in the Dynamics of Mechanical Systems

In multibody systems (MBS), such as robotic mechanical structures, it is common to adopt simplifying assumptions concerning mass properties. A frequently applied hypothesis is that the mass is continuously distributed throughout the entire mechanical structure—from the fixed base to the final kinetic assembly. Based on this assumption, the present study adopts the term “mass distribution” as a more accurate descriptor than the traditional expression “mass geometry,” which typically refers to ideal rigid bodies.

According to classical Newtonian dynamics, inertia in the translational motion is iis characterized by the system’s mass and the position of its center of mass while rotational inertia is defined by the mechanical moments of inertia—originally formalized by Leonhard Euler in 1758 — and further extended through the inertia tensor and the pseudo-inertia tensor [9]-[14].

These matrix representations are fundamental in defining dynamic quantities such as angular momentum, kinetic energy, higher-order acceleration energies, and their respective time derivatives, as required by the formulation of higher-order differential equations characteristic of analytical dynamics [15,16,17,18,19].

In this study, the authors further develop the concept of mass distribution (MD – type) based on prior work presented in [4,5]. To facilitate the analysis, several simplifying assumptions are introduced. Each kinetic subassembly, within the multibody system (MBS), is modeled as a rigid body, (see Figure 1).

At the same time, MD-type properties are using the mass and geometric integrals which are applicable exclusively to homogeneous bodies with simple or regular geometric shapes. A homogeneous body is, by definition, characterized by a constant density over an infinite number of infinitesimal mass elements$\left(dm\right)$, continuously spread throughout its geometry. In most cases, each kinetic subassembly is a compound body with an irregular geometric shape, which makes the direct application of mass integrals impossible.

As an illustrative example, the kinetic ensemble $i=1\to n$, presented in Figure 1 as part of the multibody system (MBS) is considered.

Its geometrical decomposition is detailed in Figure 2.

Based on the considerations above, each kinetic ensemble $i=1\to n$ is decomposed into a finite number of homogeneous bodies with simple or regular geometrical shapes, denoted as $\left(j\right)\in \left(i\right)$, with $j=1\to p_i$, where $p_i\in N (natural numbers)$.

The kinetic subassembly presented in Figure 2 is decomposed into simpler components $\left(j\right)\in \left(i\right)$, each assumed to have regular geometry and uniform density. This decomposition facilitates the computation of mass distribution properties such as center of mass, inertia tensor, and pseudo-inertial tensor for the full subassembly. To each component body a local reference frame centered at its geometric centroid is attached. Additionally, for each $\left(i\right)$ kinetic ensemble, a dedicated coordinate frame$\left\{iG\right\}$is defined at the geometric center of its corresponding joint, enabling consistent formulation of dynamic equations within the multibody system.

Based on these assumptions and the authors’ previous formulations, this section develops the mass distribution properties associated with such systems. These properties include: the total mass and the position of the center of masses, mechanical moments of inertia (axial, planar, polar, and products of inertia); the inertia tensor together with its generalized variation law; and the pseudo-inertia tensor, as detailed in [3,4,5].

The total mass and the position of the center of mass are first calculated for each component body$\left(j\right)\in \left(i\right)$as illustrated by the example in Figure 3:

A consistent formulation of mass distribution properties requires a clear definition of the differential mass element $\left(dm\right)$ which integrate density and geometry into a unified framework applicable to different body types, including linear, surface, and volumetric structures.

To this end, the concept of $\left(dm\right)$ is generalized as follows:

$$dm=\left\{{\rho }_V\cdot dV ; {\rho }_A\cdot dA ; {\rho }_L\cdot dL\right\}$$

(1)

where ${\rho }_L$ is the linear density characteristic to one – dimensional (1D) components such as rods or beams; ${\rho }_A$ denotes the surface density applicable to two – dimensional (2D) elements such as plates or shells and ${\rho }_V$represents the volumetric density, and is used for three – dimensional (3D) bodies.

This unified definition allows the computation of mass – related properties to be expressed in an integral form applicable to a wide variety of bodies with varying geometry and dimensions (1D, 2D or 3D).

According to [3,4,5], the position vector of the center of masses is defined by:

$${}^j\overline{r}_{C_j}={\displaystyle \frac{{\displaystyle \int {}^j\overline{r}_j\cdot dm}}{m_j}}$$

(2)

where $\left\{j\right\}$is the system attached to the mass center of the simple homogenous body $\left(j\right)$part of kinetic element $\left(i\right)$ and $m_j$ is the mass of the body $\left(j\right)$, given by:

$$m_j={\displaystyle \int dm}$$

(3)

and $\left(dm\right)$ is substituted by (1) as a function of body’s geometry.

The position of the center of mass $C_j$ must be expressed relative to the frames $\left\{i\right\}$ or $\left\{0\right\}$ which are associated with the kinematic structure of the multibody system (MBS). Consequently, the expressions are rewritten as follows:

$$\left[\begin{array}{c}{}^{\left(0\right)i}\overline{r}_{C_j}\\ 1\end{array}\right]={}_j{}^{\left(0\right)i}\left[T\right]\cdot \left[\begin{array}{c}{}^j\overline{r}_{C_j}\\ 1\end{array}\right];$$

(4)

The global center of masses for the kinetic ensemble $\left(i\right)$, denoted ${}^{\left(0\right)i}\overline{r}_{C_i}$, is computed as a mass weighted average of its component bodies’ centers of mass.

The defining expression considers the additive and subtractive contributions of each component $\left(j\right)$through the sign of operator ${\sigma }_j$:

${\sigma }_j=1$ if the component $\left(j\right)$ belongs to subassembly $\left(i\right)$ and

${\sigma }_j=-1$ when the component $\left(j\right)$ is subtracted from subassembly $\left(i\right)$ (e.g., a cavity or excluded region). The position vector ${}^{\left(0\right)i}\overline{r}_{C_i}$ defines the positioning of subassembly $\left(i\right)$relative to fixed frame$\left\{0\right\}$ or $\left\{i\right\}$, as presented below:

$${}^{\left(0\right)i}\overline{r}_{C_i}={\displaystyle \frac{{\displaystyle \sum _{j=1}^{p_i}{\sigma }_j\cdot {}^{\left(0\right)i}\overline{r}_{C_j}\cdot m_j}}{M_i}}.$$

(5)

where sometimes ${\left\{j\right\}}_{OR}\equiv {\left\{i\right\}}_{OR}$ is recommended.

The symbol ${}_j{}^{\left(0\right)i}\left[T\right]$ expresses the rotation and respectively locating matrix between frames $\left\{j\right\}$and $\left\{i\right\}$ or $\left\{0\right\}$(see also Figure 4).

Considering that the subassembly $\left(i\right)$ has been decomposed into $p_i$ simple bodies $\left(j\right)$, the application of expressions (2)-(4) to each component allows the computation of the total mass and position of the mass center of the kinetic ensemble $i=1\to n$.

The parameter $p_i$ represents the number of simple bodies resulting from the decomposition of subassembly $\left(i\right)$ and the final expressions are given as follows:

$$M_i={\displaystyle \sum _{j=1}^{p_i}{\sigma }_j\cdot m_j};$$

(6)

$$\mathrm{where}{\sigma }_j=\left\{\left(1, j\in i\right);\left(0; j\notin i\right)\right\};$$

(7)

Equation (6) defines the total mass of the kinetic ensemble, while Equation (5) provides the position vector of mass center.

2.1. The Inertial Tensor in Case of a Rigid Body

When a kinetic ensemble$i=1\to n$, from a multibody system is characterized by resultant rotational motion, its inertia properties are primarily expressed through mechanical moments of inertia and their extended representations, namely the inertial tensor and pseudo-inertia tensor [3,4]. Following the same approach used in the previous section, these inertia – related properties will first be determined for each homogeneous component body$\left(j\right)\in \left(i\right)$. An illustrative example is shown in Figure 3.

Each homogeneous body$\left(j\right)\in \left(i\right)$ has already been analyzed in terms of its mass (see Equation (2)) and position of the mass center (Equations (3)–(5)). According to Figure 3, the mass center of the body (denoted $C_j$) serves as the origin of three concurrent frames, which are described below:

$$C_j\in \left\{\left\{j\right\}; {\left\{i^{\ast }\right\}}_{OR}\equiv {\left\{i\right\}}_{OR}; {\left\{0^{\ast }\right\}}_{OR}\equiv {\left\{0\right\}}_{OR}\right\}$$

(8)

Considering the formulations of the authors, a few symbols and notations are used:

$$u_j=\left\{x_j; y_j; z_j\right\}; v_j=\left\{y_j; z_j; x_j\right\}; w_j=\left\{z_j; x_j; y_j\right\}$$

(9)

$${}^{\left(0\right)i}\overline{u}_j=\left\{{}^{\left(0\right)i}\overline{x}_j; {}^{\left(0\right)i}\overline{y}_j; {}^{\left(0\right)i}\overline{z}_j\right\},{}^j\overline{u}_{i\left(0\right)}=\left\{{}^j\overline{x}_{i\left(0\right)}; {}^j\overline{y}_{i\left(0\right)}; {}^j\overline{z}_{i\left(0\right)}\right\}$$

(10)

$${}^{\left(0\right)i}\overline{v}_j=\left\{{}^{\left(0\right)i}\overline{y}_j; {}^{\left(0\right)i}\overline{z}_j; {}^{\left(0\right)i}\overline{x}_j\right\},{}^j\overline{v}_{i\left(0\right)}=\left\{{}^j\overline{y}_{i\left(0\right)}; {}^j\overline{z}_{i\left(0\right)}; {}^j\overline{x}_{i\left(0\right)}\right\}$$

(11)

$${}^{\left(0\right)i}\overline{w}_j=\left\{{}^{\left(0\right)i}\overline{z}_j; {}^{\left(0\right)i}\overline{x}_j; {}^{\left(0\right)i}\overline{y}_j\right\}, {}^j\overline{w}_{i\left(0\right)}=\left\{{}^j\overline{z}_{i\left(0\right)}; {}^j\overline{x}_{i\left(0\right)}; {}^j\overline{y}_{i\left(0\right)}\right\}$$

(12)

where in Equation(9), the symbols refer to the Cartesian coordinates or axes of the local reference frames. Expressions (10)-(12) define the unit vectors that describe the orientation between the frames $\left\{j\right\}$, $\left\{i\right\}$ and $\left\{0\right\}$ or between$\left\{i\right\}$,$\left\{0\right\}$, and$\left\{j\right\}$respectively.

Based on these definitions, the rotation matrix can be expressed as follows:

$${}_j{}^{\left(0\right)i}\left[R\right]=\left[\begin{array}{ccc}{}^{\left(0\right)i}\overline{x}_j& {}^{\left(0\right)i}\overline{y}_j& {}^{\left(0\right)i}\overline{z}_j\end{array}\right]=\left[\begin{array}{l}{}^j\overline{x}_{i\left(0\right)}^T\\ {}^j\overline{y}_{i\left(0\right)}^T\\ {}^j\overline{z}_{i\left(0\right)}^T\end{array}\right],$$

(13)

To express the inertia properties for each homogenous body $\left(j\right)$ relative to its own mass center, the local coordinate system is positioned such that its origin coincides with the center of masses, thus leading to the following condition:

$${\displaystyle \int {}^{\left(0\right)\left(i\right)j}\overline{r}_j^{\ast }\cdot dm=0},\mathrm{where}{}^j\overline{r}_j^{\ast }={\left[\begin{array}{ccc}{}^{j}x_j^{\ast }& {}^{j}y_j^{\ast }& {}^{j}z_j^{\ast }\end{array}\right]}^T$$

(14)

where ${}^j\overline{r}_j^{\ast }$ represents the position vector of a mass element relative to the mass center of body $\left(j\right)$, relative to a local frame.

Additionally, the coordinate transformation between the frame attached to body $\left(j\right)$ and the reference frame $\left\{i\right\}$ corresponding to of ensemble $\left(i\right)$ or the frame $\left\{0\right\}$attached to the fixed basis of the robot, is defined by the following relations:

$${}^{\left(0\right)i}\overline{r}_j^{\ast }={}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^j\overline{r}_j^{\ast }, {}^{\left(0\right)i}\overline{r}_{C_j}={\displaystyle \frac{{\displaystyle \int {}^{\left(0\right)i}\overline{r}_j\cdot dm}}{m_j}}, {}^j\overline{r}_j^{\ast }={}_j{}^{\left(0\right)i}\left[R\right]_T\cdot {}^{\left(0\right)i}\overline{r}_j^{\ast },$$

(15)

where Equation (15) defines the position of body$\left(j\right)$relative to frames: $\left\{\left\{j\right\}; \left\{i^{\ast }\right\}; \left\{0^{\ast }\right\}\right\}\in C_j$.

Equation (15) is rewritten below as follows:

$$\left[\begin{array}{l}{}^{\left(0\right)i}x_j^{\ast }\\ {}^{\left(0\right)i}y_j^{\ast }\\ {}^{\left(0\right)i}z_j^{\ast }\end{array}\right]=\left[\begin{array}{l}{}^j\overline{x}_{i\left(0\right)}^T\\ {}^j\overline{y}_{i\left(0\right)}^T\\ {}^j\overline{z}_{i\left(0\right)}^T\end{array}\right]\cdot {}^j\overline{r}_j^{\ast }.$$

(16)

Another notation, essential in the present formulation, is the skew symmetric matrix associated with the position vector defined in Equation (14), as well as its transpose:

$$\left[{}^j\overline{r}_j^{\ast }\times \right]=\left[\begin{array}{ccc}0& -{}^{j}z_j^{\ast }& {}^{j}y_j^{\ast }\\ {}^{j}z_j^{\ast }& 0& -{}^{j}x_j^{\ast }\\ -{}^{j}y_j^{\ast }& {}^{j}x_j^{\ast }& 0\end{array}\right],\mathrm{and}\left[{}^j\overline{r}_j\times \right]=\left[\begin{array}{ccc}0& -{}^{j}z_j& {}^{j}y_j\\ {}^{j}z_j& 0& -{}^{j}x_j\\ -{}^{j}y_j& {}^{j}x_j& 0\end{array}\right]$$

(17)

To determine the inertia properties in either frame$\left\{i\right\}$or$\left\{0\right\}$, the position of a mass element is given by the following expressions:

$${}^{\left(0\right)i}\overline{r}_j={}^{\left(0\right)i}\overline{r}_{C_j}+{}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^j\overline{r}_j^{\ast };$$

(18)

$$\left[\begin{array}{c}{}^{\left(0\right)i}\overline{r}_j\\ 1\end{array}\right]={}_j{}^{\left(0\right)i}\left[T\right]\cdot \left[\begin{array}{c}{}^j\overline{r}_j^{\ast }\\ 1\end{array}\right]=\left[\begin{array}{cc}{}_j{}^{\left(0\right)i}\left[R\right]& {}^{\left(0\right)i}\overline{r}_{C_j}\\ \begin{array}{ccc}0& 0& 0\end{array}& 1\end{array}\right]\cdot \left[\begin{array}{c}{}^j\overline{r}_j^{\ast }\\ 1\end{array}\right]$$

(19)

where Equation(19) shows the position of elementary mass in homogeneous coordinates.

Since each homogeneous body$\left(j\right)\in \left(i\right)$ is assumed to have simple geometrical shape, its moment of inertia can be expressed as a mass integral, as established in references [3,4,5,6,7,8,9]. This integral involves the product of the differential mass element$\left(dm\right)$ and the square of its distance from a given axis or plane, and is defined as follows:

$${}^{j}I_u^{\ast }={\displaystyle \int \left({}^{j}v_j^{\ast 2}+{}^{j}w_j^{\ast 2}\right)\cdot dm}; {}^{j}I_{uv}^{\ast }={\displaystyle \int {}^{j}u_j^{\ast }\cdot }{}^{j}v_j^{\ast }\cdot dm; {}^{j}I_{uu}^{\ast }={\displaystyle \int {}^{j}u_j^{\ast 2}}\cdot dm.$$

(20)

$$\mathrm{and}{}^{\left(0\right)i}I_{ju}^{\ast }={\displaystyle \int {\delta }_{ju}^{\ast 2}}\cdot dm={\displaystyle \int \left({}^{\left(0\right)i}v_j^{\ast 2}+{}^{\left(0\right)i}w_j^{\ast 2}\right)\cdot dm}$$

(21)

,

$${}^{\left(0\right)i}I_{juv}^{\ast }={\displaystyle \int {}^{\left(0\right)i}u_j^{\ast }\cdot }{}^{\left(0\right)i}v_j^{\ast }\cdot dm; {}^{\left(0\right)i}I_{juu}^{\ast }={\displaystyle \int {}^{\left(0\right)i}u_j^{\ast 2}}\cdot dm.$$

Accordingly, the axial, centrifugal and planar mechanical moments of inertia relative to frame $\left\{j\right\}\in C_j$ are known, and are included in the following set:

$$\left\{{}^{j}I_u^{\ast }; {}^{j}I_{uv}^{\ast }; {}^{j}I_{uu}^{\ast }\right\}, \mathrm{where} j=1\to p_i \mathrm{and} \left(j\right)\in \left(i\right)$$

(22)

In the next step, the axial, centrifugal and planar mechanical moments of inertia with respect to $\left\{ \left\{i^{\ast }\right\}, \left\{0^{\ast }\right\} \right\}\in C_i$, must be computed by considering the following:

$$\left\{{}^{\left(0\right)i}I_{ju}^{\ast }, {}^{\left(0\right)i}I_{juv}^{\ast }, {}^{\left(0\right)i}I_{juu}^{\ast }\right\}, \mathrm{where} j=1\to p_i \mathrm{and} \left(j\right)\in \left(i\right)$$

(23)

The squared distance from a mass element $\left(dm\right)$ to an axis defined by the unit vector ${}^j\overline{u}_{i\left(0\right)}$ (see Figure 3) is determined according to the following expression:

$${\delta }_{ju}^{\ast 2}={\left({}^j\overline{u}_{i\left(0\right)}\times {}^j\overline{r}_j^{\ast }\right)}^T\cdot \left({}^j\overline{u}_{i\left(0\right)}\times {}^j\overline{r}_j^{\ast }\right)={}^j\overline{u}_{i\left(0\right)}^T\cdot \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot {}^j\overline{u}_{i\left(0\right)}$$

(24)

As a result, based on Equations (21)-(24), the axial mechanical moment of inertia is given by:

$${}^{\left(0\right)i}I_{ju}^{\ast }={\displaystyle \int {\delta }_{ju}^{\ast 2}}\cdot dm={}^j\overline{u}_{i\left(0\right)}^T\cdot \left\{{\displaystyle \int \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot dm}\right\}\cdot {}^j\overline{u}_{i\left(0\right)}$$

(25)

However, the mechanical moments of inertia, defined in Equation (21) can also be computed using the coordinates included in the transfer matrix formulation (16), as follows:

$$\begin{array}{c}{}^{\left(0\right)i}u_j^{\ast 2}={}^j\overline{u}_{i\left(0\right)}^T\cdot {}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{u}_{i\left(0\right)}, {}^{\left(0\right)i}v_j^{\ast 2}={}^j\overline{v}_{i\left(0\right)}^T\cdot {}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{v}_{i\left(0\right)},\\ {}^{\left(0\right)i}w_j^{\ast 2}={}^j\overline{w}_{i\left(0\right)}^T\cdot {}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{w}_{i\left(0\right)}\end{array}$$

(26)

The squared distance between a differential mass element $\left(dm\right)$ and the mass center $C_j$ is mathematically defined as:

$${}^{\left(0\right)i}u_j^{\ast 2}+{}^{\left(0\right)i}v_j^{\ast 2}+{}^{\left(0\right)i}w_j^{\ast 2}={}^{\left(0\right)i}\overline{r}_j^{\ast T}\cdot {}^{\left(0\right)i}\overline{r}_j^{\ast }\equiv {}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{r}_j^{\ast }$$

(27)

In Equation (26) the sum of the squared components of the position vector is reformulated as a scalar product allowing compact matrix representation:

$${}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{\ast T}={}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{r}_j^{\ast }\cdot I_3-\left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T$$

(28)

Substituting Equation(28) in Equation(26), expressions become:

$${}^{\left(0\right)i}u_j^{\ast 2}={}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{r}_j^{\ast }-{}^j\overline{u}_{i\left(0\right)}^T\cdot \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot {}^j\overline{u}_{i\left(0\right)}$$

(29)

$${}^{\left(0\right)i}v_j^{\ast 2}={}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{r}_j^{\ast }-{}^j\overline{v}_{i\left(0\right)}^T\cdot \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot {}^j\overline{v}_{i\left(0\right)}$$

(30)

$${}^{\left(0\right)i}w_{\ast 2}={}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{r}_j^{\ast }-{}^j\overline{w}_{i\left(0\right)}^T\cdot \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot {}^j\overline{w}_{i\left(0\right)}$$

(31)

$${}^{\left(0\right)i}v_{\ast 2}+{}^{\left(0\right)i}w_{\ast 2}={\delta }_{ju}^{\ast 2}={}^{\left(0\right)i}\overline{r}_j^{\ast T}\cdot {}^{\left(0\right)i}\overline{r}_j^{\ast }-{}^{\left(0\right)i}u_{\ast 2}={}^j\overline{u}_{i\left(0\right)}^T\cdot \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot {}^j\overline{u}_{i\left(0\right)}$$

(32)

Consequently, the mechanical moment of inertia in Equation (25) can be reformulated using the relations provided in Equation (24) or (32), as follows:

$${}^{\left(0\right)i}I_{ju}^{\ast }={}^j\overline{u}_{i\left(0\right)}^T\cdot \left\{{\displaystyle \int \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot }dm\right\}\cdot {}^j\overline{u}_{i\left(0\right)}$$

(33)

The corresponding mass integral from (25) or (33) is expressed symbolically as follows:

$$\begin{array}{c}{}^{j}I_j^{\ast }={\displaystyle \int \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot }dm\\ =\left[\begin{array}{ccc}{\displaystyle \int \left({}^{j}y_j^{\ast 2}+{}^{j}z_j^{\ast 2}\right)\cdot dm}& -{\displaystyle \int {}^{j}x_j^{\ast }\cdot {}^{j}y_j^{\ast }\cdot dm}& -{\displaystyle \int {}^{j}z_j^{\ast }\cdot {}^{j}x_j^{\ast }\cdot dm}\\ -{\displaystyle \int {}^{j}x_j^{\ast }\cdot {}^{j}y_j^{\ast }\cdot dm}& {\displaystyle \int \left({}^{j}z_j^{\ast 2}+{}^{j}x_j^{\ast 2}\right)\cdot dm}& -{\displaystyle \int {}^{j}y_j^{\ast }\cdot {}^{j}z_j^{\ast }\cdot dm}\\ -{\displaystyle \int {}^{j}z_j^{\ast }\cdot {}^{j}x_j^{\ast }\cdot dm}& -{\displaystyle \int {}^{j}y_j^{\ast }\cdot {}^{j}z_j^{\ast }\cdot dm}& {\displaystyle \int \left({}^{j}x_j^{\ast 2}+{}^{j}y_j^{\ast 2}\right)\cdot dm}\end{array}\right]\end{array}$$

(34)

In the matrix (34), the main diagonal contains the axial moments of inertia, while off – diagonal terms – symmetric and with negative sign – correspond to centrifugal moments of inertia (products of inertia). According to [2,3,4], the matrix in Equation (34) is known as the inertia tensor (axial and centrifugal) of body $\left(j\right)\in \left(i\right)$ with respect to frame $\left\{j\right\}$ evaluated at its center of masses $C_j$. Based on this inertia tensor, symbolized by (34), the axial mechanical moment of inertia relative to frames$\left\{\left\{i^{\ast }\right\}; \left\{0^{\ast }\right\}\right\}\in C_j$ is determined as:

$${}^{\left(0\right)i}I_{ju}^{\ast }={}^j\overline{u}_{i\left(0\right)}^T\cdot {}^{j}I_j^{\ast }\cdot {}^j\overline{u}_{i\left(0\right)}$$

(35)

$$\underset{\left(3\times 3\right)}{Diag}\left\{\left[\begin{array}{c}{}^j\overline{x}_{i\left(0\right)}^T\\ {}^j\overline{y}_{i\left(0\right)}^T\\ {}^j\overline{z}_{i\left(0\right)}^T\end{array}\right]\cdot {}^{j}I_j^{\ast }\cdot \left[\begin{array}{ccc}{}^j\overline{x}_{i\left(0\right)}& {}^j\overline{y}_{i\left(0\right)}& {}^j\overline{z}_{i\left(0\right)}\end{array}\right]\right\}=\underset{\left(3\times 3\right)}{Diag}\left\{{}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_j^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T\right\}$$

(36)

Based on the definitions provided in Equations (23) and (21),the centrifugal moment of inertia is computed as follows:

$${}^{\left(0\right)i}I_{juv}^{\ast }={\displaystyle \int {}^{\left(0\right)i}u_j^{\ast }{}_{\cdot }}{}^{\left(0\right)i}v_j^{\ast }\cdot dm={}^j\overline{u}_{i\left(0\right)}^T\cdot {\displaystyle \int {}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{\ast T}\cdot } dm\cdot {}^j\overline{v}_{i\left(0\right)}$$

(37)

The $\left(3\times 3\right)$ matrix, defined by Equation (28) is rewritten as presented below:

$${}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{\ast T}={}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{r}_j^{\ast }\cdot I_3-\left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T$$

(38)

When substituted in Equation (37), this expression leads to the following transformation:

$${}^{\left(0\right)i}I_{juv}^{\ast }={\displaystyle \int \left[{}^j\overline{r}_j^{\ast T}\cdot {}^j\overline{r}_j^{\ast }\cdot \right.}{}^j\overline{u}_{i\left(0\right)}^T\cdot {}^j\overline{v}_{i\left(0\right)}-\left.{}^j\overline{u}_{i\left(0\right)}^T\cdot \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot {}^j\overline{v}_{i\left(0\right)}\right]\cdot dm$$

(39)

$$-{}^j\overline{u}_{i\left(0\right)}^T\cdot \left\{{\displaystyle \int \left[{}^j\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^j\overline{r}_j^{\ast }\times \right]}^T\cdot dm}\right\}\cdot {\overline{v}}_{i\left(0\right)}=-{}^j\overline{u}_{i\left(0\right)}^T\cdot {}^{j}I_j^{\ast }\cdot {}^j\overline{v}_{i\left(0\right)}={}^{\left(0\right)i}I_{juv}^{\ast }$$

(40)

where (34) has been substituted.

Using Equation (13) and its transpose, the centrifugal moments of inertia become:

$${}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_j^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T-\underset{\left(3\times 3\right)}{Diag}\left\{{}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_j^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T\right\}=\left[\begin{array}{ccc}& -{}^{\left(0\right)i}I_{jxy}^{\ast }& -{}^{\left(0\right)i}I_{jxz}^{\ast }\\ -{}^{\left(0\right)i}I_{jyx}^{\ast }& & -{}^{\left(0\right)i}I_{jyz}^{\ast }\\ -{}^{\left(0\right)i}I_{jzx}^{\ast }& -{}^{\left(0\right)i}I_{jzy}^{\ast }& \end{array}\right]$$

(41)

Equation (36) and (41) are included in the matrix of axial and centrifugal moments of inertia:

$${}^{\left(0\right)i}I_j^{\ast }={\displaystyle \int \left[{}^{\left(0\right)i}\overline{r}_j^{\ast }\times \right]\cdot {\left[{}^{\left(0\right)i}\overline{r}_j^{\ast }\times \right]}^T\cdot }dm={}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_j^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T$$

(42)

$${}^{\left(0\right)i}I_j^{\ast }={}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_j^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T=\left[\begin{array}{ccc}{}^{\left(0\right)i}I_{jx}^{\ast }& -{}^{\left(0\right)i}I_{jxy}^{\ast }& -{}^{\left(0\right)i}I_{jxz}^{\ast }\\ -{}^{\left(0\right)i}I_{jyx}^{\ast }& {}^{\left(0\right)i}I_{jy}^{\ast }& -{}^{\left(0\right)i}I_{jyz}^{\ast }\\ -{}^{\left(0\right)i}I_{jzx}^{\ast }& -{}^{\left(0\right)i}I_{jzy}^{\ast }& {}^{\left(0\right)i}I_{jz}^{\ast }\end{array}\right]$$

(43)

The matrices defined in Equations (42) and (43), presented in their expanded form, constitute the inertia tensor axial and centrifugal of the body $\left(j\right)\in \left(i\right)$ relative to frame $\left\{\left\{i^{\ast }\right\}; \left\{0^{\ast }\right\}\right\}\in C_j$ and evaluated at the mass center$C_j$. Simultaneously, these expressions represent the variation law of the inertia tensor with respect to concurrent reference frames located at the mass center $C_j$. Based on Equation (22) and (23), the following section addresses the computation of the planar moment of inertia, as defined by:

$${}^{\left(0\right)i}I_{juu}^{\ast }={\displaystyle \int {}^{\left(0\right)i}u_j^{\ast 2}\cdot } dm={}^j\overline{u}_{i\left(0\right)}^T\cdot {\displaystyle \int {}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{\ast T}\cdot } dm\cdot {}^j\overline{u}_{i\left(0\right)}$$

(44)

By substituting Equation (38) into Equation (44) and evaluating the mass integrals, the following result is obtained:

$${}^{j}I_{pj}^{\ast }={\displaystyle \int {}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{T\ast }\cdot } dm=\left[\begin{array}{ccc}{}^{j}I_{xx}^{\ast }& {}^{j}I_{xy}^{\ast }& {}^{j}I_{xz}^{\ast }\\ {}^{j}I_{yx}^{\ast }& {}^{j}I_{yy}^{\ast }& {}^{j}I_{yz}^{\ast }\\ {}^{j}I_{zx}^{\ast }& {}^{j}I_{zy}^{\ast }& {}^{j}I_{zz}^{\ast }\end{array}\right]$$

(45)

The matrix of inertia moments defined in Equation (45) is referred to as the planar and centrifugal inertia tensor of body $\left(j\right)\in \left(i\right)$relative to reference frame $\left\{j\right\}$ and applied at the mass center$C_j$. Based on this tensor representation, the planar mechanical moments of inertia of body $\left(j\right)\in \left(i\right)$relative to frame$\left\{\left\{i^{\ast }\right\}; \left\{0^{\ast }\right\}\right\}\in C_j$ are determined as:

$${}^{\left(0\right)i}I_{juu}^{\ast }={}^j\overline{u}_{i\left(0\right)}^T\cdot {}^{j}I_{pj}^{\ast }\cdot {}^j\overline{u}_{i\left(0\right)}$$

(46)

$$\underset{\left(3\times 3\right)}{Diag}\left\{{}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_{pj}^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T\right\}=\underset{\left(3\times 3\right)}{Diag}\left\{\left[\begin{array}{c}{}^j\overline{x}_{i\left(0\right)}^T\\ {}^j\overline{y}_{i\left(0\right)}^T\\ {}^j\overline{z}_{i\left(0\right)}^T\end{array}\right]\cdot {}^{j}I_{pj}^{\ast }\cdot \left[\begin{array}{ccc}{}^j\overline{x}_{i\left(0\right)}& {}^j\overline{y}_{i\left(0\right)}& {}^j\overline{z}_{i\left(0\right)}\end{array}\right]\right\}$$

(47)

Based on Equation (41), in which the inertia tensor from Equation (45) is substituted, the centrifugal moments of inertia are determined using the following expressions:

$$\left[\begin{array}{ccc}& {}^{\left(0\right)i}I_{jxy}^{\ast }& {}^{\left(0\right)i}I_{jxz}^{\ast }\\ {}^{\left(0\right)i}I_{jyx}^{\ast }& & {}^{\left(0\right)i}I_{jyz}^{\ast }\\ {}^{\left(0\right)i}I_{jzx}^{\ast }& {}^{\left(0\right)i}I_{jzy}^{\ast }& \end{array}\right]={}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_{pj}^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T-\underset{\left(3\times 3\right)}{Diag}\left\{{}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_{pj}^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T\right\}$$

(48)

Therefore, Equations (47) and (48) are included into a matrix that represents the planar and centrifugal moments of inertia:

$${}^{\left(0\right)i}I_{pj}^{\ast }={\displaystyle \int {}^{\left(0\right)i}\overline{r}_j^{\ast }\cdot {}^{\left(0\right)i}\overline{r}_j^{T\ast }\cdot } dm={}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_{pj}^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T$$

(49)

$${}^{\left(0\right)i}I_{pj}^{\ast }={}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_{pj}^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T=\left[\begin{array}{ccc}{}^{\left(0\right)i}I_{jxx}^{\ast }& {}^{\left(0\right)i}I_{jxy}^{\ast }& {}^{\left(0\right)i}I_{jxz}^{\ast }\\ {}^{\left(0\right)i}I_{jyx}^{\ast }& {}^{\left(0\right)i}I_{jyy}^{\ast }& {}^{\left(0\right)i}I_{jyz}^{\ast }\\ {}^{\left(0\right)i}I_{jzx}^{\ast }& {}^{\left(0\right)i}I_{jzy}^{\ast }& {}^{\left(0\right)i}I_{jzz}^{\ast }\end{array}\right]$$

(50)

The matrices defined in Equations (49) and (50), in their expanded form, represent the planar – centrifugal inertia tensor of body $\left(j\right)\in \left(i\right)$ relative to reference frame$\left\{\left\{i^{\ast }\right\}; \left\{0^{\ast }\right\}\right\}\in C_j$ applied at its center of masses $C_j$. At the same time, Equations (42) and (43) describe the variation law of the inertia tensor with respect to concurrent reference frames located at the mass center $C_j$.

In many cases, the fundamental concepts and theorems of analytical dynamics are expressed in a matrix form. Therefore, the position vectors from Equation (19) are rewritten using homogeneous coordinates, which leads to the following matrix representation:

$${}^{j}I_{psj}^{\ast }={\displaystyle \int \left(\begin{array}{c}{}^j\overline{r}_j^{\ast }\\ 1\end{array}\right)\cdot \left(\begin{array}{cc}{}^j\overline{r}_j^{T\ast }& 1\end{array}\right)\cdot } dm=\left[\begin{array}{cc}{\displaystyle \int {}^j\overline{r}_j^{\ast }\cdot {}^j\overline{r}_j^{\ast T}\cdot dm}& {\displaystyle \int {}^j\overline{r}_j^{\ast }\cdot dm}\\ {\displaystyle \int {}^j\overline{r}_j^{\ast T}\cdot dm}& m_j\end{array}\right]=\left[\begin{array}{cc}{}^{j}I_{pj}^{\ast }& {\left[0\right]}_{\left(3\times 1\right)}\\ {\left[0\right]}_{\left(1\times 3\right)}& m_j\end{array}\right]$$

(51)

In this formulation, Equation (14) is substituted, incorporating the planar – centrifugal inertia tensor (from Equation (45)), the static moments and the mass of the body. The resulting matrix defines the pseudo-inertia tensor of the body relative to reference frame $\left\{j\right\}$, applied at the mass center$C_j$. Accordingly, the expression of the pseudo-inertia tensor relative to frame $\left\{\left\{i^{\ast }\right\}; \left\{0^{\ast }\right\}\right\}\in C_j$ is given by:

$$\begin{array}{c}{}^{\left(0\right)i}I_{psj}^{\ast }={\displaystyle \int \left(\begin{array}{c}{}^{\left(0\right)i}\overline{r}_j^{\ast }\\ 1\end{array}\right)\cdot \left(\begin{array}{cc}{}^{\left(0\right)i}\overline{r}_j^{T\ast }& 1\end{array}\right)\cdot } dm\\ =\left[\begin{array}{cc}{\displaystyle \int {}^{\left(0\right)i}\overline{r}_j^{\ast }\cdot {}^{\left(0\right)i}\overline{r}_j^{\ast T}\cdot dm}& {\displaystyle \int {}^{\left(0\right)i}\overline{r}_j^{\ast }\cdot dm}\\ {\displaystyle \int {}^{\left(0\right)i}\overline{r}_j^{\ast T}\cdot dm}& m_j\end{array}\right]=\left[\begin{array}{cc}{}^{\left(0\right)i}I_{pj}^{\ast }& {\left[0\right]}_{\left(3\times 1\right)}\\ {\left[0\right]}_{\left(1\times 3\right)}& m_j\end{array}\right]\end{array}$$

(52)

where the conditions (14) are substituted again. But, (52) can be written in another form:

$${}^{\left(0\right)i}I_{psj}^{\ast }={}_j{}^{\left(0\right)i}\left[T\right]^{\ast }\cdot {}^{j}I_{pj}^{\ast }\cdot {}^j{}_{\left(0\right)i}\left[T\right]_{\ast }^T$$

(53)

$$\mathrm{where}{}_j{}^{\left(0\right)i}\left[T\right]^{\ast }=\left[\begin{array}{cc}{}_j{}^{\left(0\right)i}\left[R\right]& {\left[0\right]}_{\left(3\times 1\right)}\\ {\left[0\right]}_{\left(1\times 3\right)}& 1\end{array}\right]$$

(54)

The square matrix defined in Equation (53), symmetric and positive defined, represents the variation law of the pseudo-inertial tensor with respect to concurrent reference frames located at the mass center $C_j$. The matrix defined with (54) is the homogenous transformation matrix which is used to determine a pure rotation from frame $\left\{j\right\}$to frame $\left\{i\right\}$. The components of this matrix are represented by ${}_j{}^{\left(0\right)i}\left[R\right]$ which is the rotation matrix that characterizes the rotation of frame $\left\{j\right\}$ relative to frames $\left\{i\right\}$ or $\left\{0\right\}$. The rightmost column of the homogenous transformation and the bottom row contain zeros and a scalar 1. The matrix does not include translation, making it a pure rotational transformation into homogenous space.

2.2. The Generalized Variation of the Inertial Tensor

The mass properties for any homogeneous body $\left(j\right)\in \left(i\right)$ assumed to have a simple geometric configuration [2], can be determined from the following predefined parameters:

$$\left\{m_j; {}^j\overline{r}_{C_j}{}^{j}I_j^{\ast }; {}^{j}I_j^{\ast }; {}^{j}I_{pj}^{\ast }; {}^{j}I_{psj}^{\ast }; j=1\to p_i; \left(j\right)\in \left(i\right)\right\}$$

(55)

Based on the expressions provided in Section 2 of this paper, the mass properties are computed relative to reference frame$\left\{\left\{i^{\ast }\right\}; \left\{0^{\ast }\right\}\right\}\in C_j$, as follows:

$\left\{m_j; {}^{\left(0\right)i}\overline{r}_{C_j}^{\ast }{}^{\left(0\right)i}I_j^{\ast }; {}^{\left(0\right)i}I_j^{\ast }; {}^{\left(0\right)i}I_{pj}^{\ast }; {}^{\left(0\right)i}I_{psj}^{\ast }; j=1\to p_i\right\}$.

In the following steps, the axial, centrifugal and planar mechanical moments of inertia for each homogeneous body $\left(j\right)\in \left(i\right)$ are to be computed with respect to frames $\left\{i\right\}$ and $\left\{0\right\}$, in accordance with the kinematic structure of MBS (as presented in Figure 1):

$$\left\{{}^{\left(0\right)i}I_u; {}^{\left(0\right)i}I_{uv}; {}^{\left(0\right)i}I_{uu}; where j=1\to p_i; \left(j\right)\in \left(i\right)\right\}$$

(56)

These quantities are included in the generalized variation law of the inertial and pseudo-inertial tensors, as presented below:

$$\left\{{}^{\left(0\right)i}I_j; {}^{\left(0\right)i}I_{pj}; {}^{\left(0\right)i}I_{psj}; where j=1\to p_i; \left(j\right)\in \left(i\right)\right\}$$

(57)

To begin with, the generalized variation law of the axial and centrifugal inertia tensor for body $\left(j\right)\in \left(i\right)$ relative to reference frame$\left\{\left\{i\right\}; \left\{0\right\}\right\}$is determined as follows:

$${}^{\left(0\right)i}I_j={\displaystyle \int \left[{}^{\left(0\right)i}\overline{r}_j\times \right]\cdot {\left[{}^{\left(0\right)i}\overline{r}_j\times \right]}^T\cdot } dm={}^{\left(0\right)i}I_{C_j}+{}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_j^{\ast }\cdot {}^j{}_{\left(0\right)i}\left[R\right]_T={}^{\left(0\right)i}I_{C_j}+{}^{\left(0\right)i}I_j^{\ast }$$

(58)

$${}^{\left(0\right)i}I_{C_j}=m_j\cdot \left[{}^{\left(0\right)i}\overline{r}_{C_j}\times \right]\cdot {\left[{}^{\left(0\right)i}\overline{r}_{C_j}\times \right]}^T=\left[\begin{array}{ccc}{}^{\left(0\right)i}y_{C_j}^2+{}^{\left(0\right)i}z_{C_j}^2& -{}^{\left(0\right)i}x_{C_j}\cdot {}^{\left(0\right)i}y_{C_j}& -{}^{\left(0\right)i}z_{C_j}\cdot {}^{\left(0\right)i}x_{C_j}\\ -{}^{\left(0\right)i}x_{C_j}\cdot {}^{\left(0\right)i}y_{C_j}& {}^{\left(0\right)i}z_{C_j}^2+{}^{\left(0\right)i}x_{C_j}^2& -{}^{\left(0\right)i}y_{C_j}\cdot {}^{\left(0\right)i}z_{C_j}\\ -{}^{\left(0\right)i}z_{C_j}\cdot {}^{\left(0\right)i}x_{C_j}& -{}^{\left(0\right)i}y_{C_j}\cdot {}^{\left(0\right)i}z_{C_j}& {}^{\left(0\right)i}x_{C_j}^2+{}^{\left(0\right)i}y_{C_j}^2\end{array}\right]$$

(59)

The expression given in Equation (59) is referred to as the axial and centrifugal inertia matrix of the mass center $C_j$ relative to reference frames$\left\{\left\{i\right\}; \left\{0\right\}\right\}$.

Based on Equation (57), the generalized variation law for the planar and centrifugal inertia tensor of the body $\left(j\right)\in \left(i\right)$ relative to frame $\left\{\left\{i\right\}; \left\{0\right\}\right\}$is formulated as follows:

$${}^{\left(0\right)i}I_{pj}={\displaystyle \int {}^{\left(0\right)i}\overline{r}_j\cdot {}^{\left(0\right)i}\overline{r}_j^T\cdot } dm={}^{\left(0\right)i}I_{p{C}_j}+{}_j{}^{\left(0\right)i}\left[R\right]\cdot {}^{j}I_{pj}^{\ast }\cdot {}_j{}^{\left(0\right)i}\left[R\right]^T={}^{\left(0\right)i}I_{p{C}_j}+{}^{\left(0\right)i}I_{pj}^{\ast }$$

(60)

$$\mathrm{where},{}^{\left(0\right)i}I_{p{C}_j}=m_j\cdot {}^{\left(0\right)i}\overline{r}_{C_j}\cdot {}^{\left(0\right)i}\overline{r}_{C_j}^T=\left[\begin{array}{ccc}{}^{\left(0\right)i}x_{C_j}^2& {}^{\left(0\right)i}x_{C_j}\cdot {}^{\left(0\right)i}y_{C_j}& {}^{\left(0\right)i}z_{C_j}\cdot {}^{\left(0\right)i}x_{C_j}\\ {}^{\left(0\right)i}x_{C_j}\cdot {}^{\left(0\right)i}y_{C_j}& {}^{\left(0\right)i}y_{C_j}^2& {}^{\left(0\right)i}y_{C_j}\cdot {}^{\left(0\right)i}z_{C_j}\\ {}^{\left(0\right)i}z_{C_j}\cdot {}^{\left(0\right)i}x_{C_j}& {}^{\left(0\right)i}y_{C_j}\cdot {}^{\left(0\right)i}z_{C_j}& {}^{\left(0\right)i}z_{C_j}^2\end{array}\right]$$

(61)

The expression (84) defines the planar and centrifugal inertia matrix evaluated at the mass center $C_j$, with respect to the frames$\left\{\left\{i\right\}; \left\{0\right\}\right\}$.

The pseudo – inertia tensor variation law of the body $\left(j\right)\in \left(i\right)$ with respect to frames $\left\{\left\{i\right\}; \left\{0\right\}\right\}$, obtained from Equations (51)-(53), is given below:

$$\begin{array}{c}{}^{\left(0\right)i}I_{psj}={\displaystyle \int \left[\begin{array}{c}{}^{\left(0\right)i}\overline{r}_j\\ 1\end{array}\right]\cdot \left[\begin{array}{cc}{}^{\left(0\right)i}\overline{r}_j^T& 1\end{array}\right]\cdot } dm\\ ={}_{\left(0^{\ast }\right)i^{\ast }}{}^{\left(0\right)i}\left[T\right]\cdot {}^{\left(0\right)i}I_{psj}^{\ast }\cdot {}_{\left(0^{\ast }\right)i^{\ast }}{}^{\left(0\right)i}\left[T\right]_T=\left[\begin{array}{cc}{}^{\left(0\right)i}I_{pj}& {}^{\left(0\right)i}\overline{r}_{C_j}\cdot m_j\\ {}^{\left(0\right)i}\overline{r}_{C_j}^T\cdot m_j& m_j\end{array}\right]\end{array}$$

(62)

$$\mathrm{where}{}_{\left(0^{\ast }\right)i^{\ast }}{}^{\left(0\right)i}\left[T\right]=\left[\begin{array}{cc}{I}_{3\times 3}& {}^{\left(0\right)i}\overline{r}_{C_j}\\ {\left[0\right]}_{\left(1\times 3\right)}& 1\end{array}\right]$$

(63)

is the locating matrix between frames $\left\{0^{*}\right\}/\left\{i^{*}\right\}$ versus $\left\{0\right\}/\left\{i\right\}$.

2.3. The Inertial Tensor for Multibody Systems (MBS)

The generalized variation laws of the inertial and pseudo-inertial tensors, as defined in Equation (57), have been defined for every homogeneous body $\left(j\right)\in \left(i\right)$ relative to frames$\left\{\left\{i\right\}; \left\{0\right\}\right\}$. In this section, based on the previously established expressions, the inertia properties for each kinetic ensemble$\left(i=1\to n\right)$will be determined. These include the axial, centrifugal and planar mechanical moments of inertia, in accordance with the models presented in references [4] and [5]:

$${}^{\left(0\right)i}I_u; {}^{\left(0\right)i}I_{uv}; {}^{\left(0\right)i}I_{uu},\mathrm{where} i=1\to n$$

(64)

These values, already defined, are now included in the expressions for the inertial and pseudo-inertial tensors:

$${}^{\left(0\right)i}I_i; {}^{\left(0\right)i}I_{pi}; {}^{\left(0\right)i}I_{psi},\mathrm{where} i=1\to n$$

(65)

Therefore, the generalized variation law for the axial and centrifugal inertia tensor of the kinetic ensemble $\left(i\right)$ relative to frame$\left\{\left\{i\right\}; \left\{0\right\}\right\}$is determined below:

$${}^{\left(0\right)i}I_i={\displaystyle \sum _{j=1}^{p_i}{\sigma }_j}\cdot {}^{\left(0\right)i}I_j={\displaystyle \sum _{j=1}^{p_i}{\sigma }_j}\cdot {}^{\left(0\right)i}I_{C_j}+{\displaystyle \sum _{j=1}^{p_i}{\sigma }_j}\cdot {}^{\left(0\right)i}I_j^{\ast }={}^{\left(0\right)i}I_{C_i}+{}^{\left(0\right)i}I_i^{\ast }$$

(66)

$${}^{\left(0\right)i}I_i^{\ast }={\displaystyle \sum _{j=1}^{p_i}{\sigma }_j}\cdot {}^{\left(0\right)i}I_j^{\ast }=\left[\begin{array}{ccc}{}^{\left(0\right)i}I_x^{\ast }& -{}^{\left(0\right)i}I_{xy}^{\ast }& -{}^{\left(0\right)i}I_{xz}^{\ast }\\ -{}^{\left(0\right)i}I_{yx}^{\ast }& {}^{\left(0\right)i}I_y^{\ast }& -{}^{\left(0\right)i}I_{yz}^{\ast }\\ -{}^{\left(0\right)i}I_{zx}^{\ast }& -{}^{\left(0\right)i}I_{zy}^{\ast }& {}^{\left(0\right)i}I_z^{\ast }\end{array}\right]$$

(67)

The generalized variation law of the planar and centrifugal inertia tensor of the kinetic ensemble$\left(i\right)$, relative to reference frame$\left\{\left\{i\right\}; \left\{0\right\}\right\}$, is established with expressions:

$${}^{\left(0\right)i}I_{pi}={\displaystyle \sum _{j=1}^{p_i}{\sigma }_j}\cdot {}^{\left(0\right)i}I_{pj}={}^{\left(0\right)i}I_{p{C}_i}+{}^{\left(0\right)i}I_{pi}^{\ast }={\displaystyle \sum _{j=1}^{p_i}{\sigma }_j}\cdot {}^{\left(0\right)i}I_{p{C}_j}+{\displaystyle \sum _{j=1}^{p_i}{\sigma }_j}\cdot {}^{\left(0\right)i}I_{pj}^{\ast }$$

(68)

$${}^{\left(0\right)i}I_{pi}={\displaystyle \sum _{j=1}^{p_i}{\sigma }_j}\cdot {}^{\left(0\right)i}I_{pj}=\left[\begin{array}{ccc}{}^{\left(0\right)i}I_{xx}& {}^{\left(0\right)i}I_{xy}& {}^{\left(0\right)i}I_{xz}\\ {}^{\left(0\right)i}I_{yx}& {}^{\left(0\right)i}I_{yy}& {}^{\left(0\right)i}I_{yz}\\ {}^{\left(0\right)i}I_{zx}& {}^{\left(0\right)i}I_{zy}& {}^{\left(0\right)i}I_{zz}\end{array}\right]$$

(69)

Following the formulation presented in Equation (65), the pseudo-inertial tensor corresponding to the kinetic ensemble $\left(i\right)$ relative to frames$\left\{\left\{i\right\}; \left\{0\right\}\right\}$is determined as:

$$\begin{array}{c}{}^{\left(0\right)i}I_{psi}=\left[\begin{array}{cc}{}^{\left(0\right)i}I_{pi}& M_i\cdot {}^{\left(0\right)i}\overline{r}_{C_i}\\ M_i\cdot {}^{\left(0\right)i}\overline{r}_{C_i}^T& M_i\end{array}\right]\\ ={\displaystyle \sum _{j=1}^{p_i}{\sigma }_j\cdot {}^{\left(0\right)i}I_{psj}}=\left[\begin{array}{cc}{\displaystyle \sum _{j=1}^{p_i}{\sigma }_j\cdot {}^{\left(0\right)i}I_{pj}}& {\displaystyle \sum _{j=1}^{p_i}{\sigma }_j\cdot {}^{\left(0\right)i}\overline{r}_{C_j}\cdot m_j}\\ {\displaystyle \sum _{j=1}^{p_i}{\sigma }_j\cdot {}^{\left(0\right)i}\overline{r}_{C_j}^T\cdot m_j}& {\displaystyle \sum _{j=1}^{p_i}{\sigma }_j\cdot m_j}\end{array}\right]\end{array}$$

(70)

The expression defined above contains the planar and centrifugal inertia tensor (as presented in Equation (69)), along with the static moments (Equation (5)) and the total mass (from Equation (6)) of the kinetic ensemble$\left(i\right)$.

As demonstrated in [3,4,5], both the inertial and pseudo-inertial tensors, are essential components in the matrix formulations of fundamental concepts in analytical dynamics, including angular momentum, kinetic energy and acceleration energy.

3. The Acceleration Energies of Higher Order

Unlike other research papers addressing these advanced concepts, the main objective of this section is to present, in a matrix form, the fundamental expressions defining the higher – order acceleration energies. By applying the differential principle in generalized form (a generalization of the Lagrange–D’Alembert principle), the dynamic equations of fast-moving mechanical systems naturally incorporate these higher-order acceleration energies. As documented in the scientific literature [20,21,22], the foundations of these equations were laid in 1879 by J. W. Gibbs who formulated the differential equations of motion. Subsequently, in 1899, Paul Appell extended this work through a more elaborate analysis [23,24]. The outcome of this development is a set of equations now known as Gibbs–Appell equations. These equations are applied to holonomic and non-holonomic systems.

The study presented in this section was carried out by considering the holonomic mechanical systems, for which, the Gibbs–Appell equations along with their higher-order derivatives, are refined to suit this type of systems. Furthermore, this section emphasizes the importance of higher acceleration energies as central functions, in the analysis of the dynamics of high – speed mechanical systems.

In this framework, the kinetic energy is substituted by the acceleration energy, also referred to as Appell’s function or “kinetic energy of acceleration” [3,5].

Unlike the studies mentioned above, the authors have developed explicit expressions for the first, second, and third – order acceleration energy, specific to mechanical systems characterized by high – speed motions [18].

To determine the general expressions for the acceleration energies of pth order, the corresponding differential matrices are introduced in the following subsection.

3.1. The Differential Matrices in the Advanced Kinematics

Based on the formulations of homogenous transformations and matrix exponentials, the differential matrices of these homogeneous transformations are developed. These matrices – also known as dynamics matrices – are essential components in the matrix representation of the acceleration energies, particularly for mechanical systems characterized by fast motions. As stated in [18], the dynamics matrices comprise first-, second- and higher – order differential matrices, which can be determined either by directly applying the partial derivatives to the homogenous transformations or by using the properties of matrix exponential functions.

The components of the differential matrices are represented by the submatrices corresponding to the rotation $\left(R\right)$ and position $\left(\overline{p}\right)$, respectively. The first-order differential matrices$A_{ij\left(k,m,p\right)}$, representing the time derivatives of the homogenous transformations are obtained as follows:

$$A_{ij\left(k,m,p\right)}=\left[\begin{array}{cccc}& A_{ij\left(k,m,p\right)}\left(R\right)& & A_{ij\left(k,m,p\right)}\left(\overline{p}\right)\\ 0& 0& 0& 0\end{array}\right]={\displaystyle \frac{\partial }{\partial q_j}}\left\{{}_i{}^0\left[T\right]\right\}={}_j{}^0\left[T\right]\cdot U_j\cdot {}_i{}^j\left[T\right].$$

(71)

By employing the matrix exponentials formalism, the two components represented by the rotation transformation $\left(R\right)$, and position $\left(\overline{p}\right)$, from Equation (71), are further defined:

$$\begin{array}{c}{A}_{ij}\left(R\right)={\displaystyle \frac{\partial }{\partial q_j}}\left\{{}_i{}^0\left[R\right]\right\}\\ =\left\{\exp \left\{{\displaystyle \sum _{k=0}^{j-1}\left[{\overline{k}}_k^{\left(0\right)}\times \right]\cdot q_k\cdot {\Delta }_k}\right\}\right\}\cdot \left[{\overline{k}}_j^{\left(0\right)}\times \right]\cdot {\Delta }_j\cdot \exp \left\{{\displaystyle \sum _{l=j}^i\left[{\overline{k}}_l^{\left(0\right)}\times \right]\cdot q_l\cdot {\Delta }_l}\right\}\cdot R_{i0}^{\left(0\right)},\end{array}$$

(72)

$$\begin{array}{c} A_{ij}\left(\overline{p}\right)={\displaystyle \frac{\partial {\overline{p}}_i}{\partial q_j}} \\ =\left\{\exp \left[{\displaystyle \sum _{k=0}^{j-1}\left({\overline{k}}_k^{(0)}\times \right)\cdot q_k\cdot {\Delta }_k}\right]\right\}\cdot \left[\left({\overline{p}}_j^{\left(0\right)} \times \right) \cdot {\overline{k}}_j^{\left(0\right)}\cdot {\Delta }_j+\left(1-{\Delta }_j\right)\cdot {\overline{k}}_j^{\left(0\right)}\right]+\\ +\exp \left[{\displaystyle \sum _{l=j}^i\left({\overline{k}}_l^{(0)}\times \right)\cdot q_l\cdot {\Delta }_l}\right]\cdot {\overline{p}}_i^{\left(0\right)}+{\Delta }_j\cdot \exp \left[{\displaystyle \sum _{k=0}^{j-1}\left({\overline{k}}_k^{(0)}\times \right)\cdot q_k\cdot {\Delta }_k}\right]\cdot A_{ij}^{\ast }\left(\overline{p}\right)\end{array}$$

(73)

$$\begin{gathered} \mathrm{where} A_{ij}^{\ast }\left(\overline{p}\right) ={\displaystyle \sum _{l=j}^i\left\{\exp \left[{\displaystyle \sum _{m=j-1}^{l-1}\left({\overline{k}}_m^{(0)}\times \right)\cdot q_m\cdot {\Delta }_m\cdot {\delta }_m}\right]\right\}}\cdot {\overline{b}}_l \\ {\delta }_m = \left\{\left(0 , m=j-1\right) , \left(1 , m\ge i\right)\right\}. \end{gathered}$$

(74)

$$\mathrm{and}{\Delta }_u=\left\{\left(1, q_u-\mathrm{rotation}\right); \left(0, q_u-\mathrm{translation}\right)\right\},u=\left\{ i; j; k; m; p \right\}$$

(75)

The differential matrix of second order is defined with matrix exponential functions:

$$A_{ijk}=\left[\begin{array}{cccc}& A_{ijk}\left(R\right)& & A_{ijk}\left(\overline{p}\right)\\ 0& 0& 0& 0\end{array}\right]={\displaystyle \frac{{\partial }^2}{\partial q_j\cdot \partial q_k}}\left\{{}_i{}^0\left[T\right]\right\}={}_k{}^0\left[T\right]\cdot U_k\cdot {}_j{}^k\left[T\right]\cdot U_j\cdot {}_i{}^j\left[T\right]$$

(76)

$$A_{ijk}\left(R\right)={\displaystyle \frac{{\partial }^2}{\partial q_j\cdot \partial q_k}}\left\{{}_i{}^0\left[R\right]\right\}=\left\{\exp \left\{{\displaystyle \sum _{l=0}^{k-1}\left({\overline{k}}_l^{\left(0\right)}\times \right)\cdot q_l\cdot {\Delta }_l}\right\}\right\}\cdot \left({\overline{k}}_k^{\left(0\right)}\times \right)\cdot {\Delta }_k\cdot A_{ijk}^{\ast }\left(R\right)$$

(77)

$$\mathrm{where}{A}_{ijk}^{\ast }\left(R\right)=\left\{\exp \left\{{\displaystyle \sum _{m=k}^{j-1}\left({\overline{k}}_m^{\left(0\right)}\times \right)\cdot q_m\cdot {\Delta }_m}\right\}\right\}\cdot \left({\overline{k}}_m^{\left(0\right)}\times \right)\cdot {\Delta }_m\cdot \left\{\exp \left\{{\displaystyle \sum _{p=m}^i\left({\overline{k}}_p^{\left(0\right)}\times \right)\cdot q_p\cdot {\Delta }_p}\right\}\right\}\cdot R_{i0}^{\left(0\right)},$$

(78)

$$A_{ijk}\left(\overline{p}\right)={\displaystyle \frac{\partial }{\partial q_k}}\left.\left[A_{ij}\left(\overline{p}\right)\right.\right]={\displaystyle \frac{{\partial }^2}{\partial q_k\cdot \partial q_m}}\left.\left[A_{ij}\left(\overline{p}\right)\right.\right].$$

(79)

The sub-matrices included in Equation (72) and (79) are determined according to [3,4,5], by using matrix exponential functions. The following symbols are explained:

$U_{j(k,m)}$ represents the derivative matrix operator (Uicker operator).

${}^0\overline{k}_u$ - it represents the unit vector of the driving axis in the initial configuration (zero)

$q_u$ - it represents the generalized variable of the driving axis.

${\overline{p}}_i^{\left(0\right)}$ and $R_{i0}^{\left(0\right)}$ defines the position vector in the initial configuration and the orientation matrix of the system $\left\{i\right\}$ relative to $\left\{0\right\}$, respectively;

${\overline{b}}_l$ - it represents the 4th column of the exponential of a homogeneous transformation and is defined according to [25], as follows:

$${\overline{b}}_l=\left\{I_3\cdot s{q}_l+\left\{{\overline{k}}_l^{\left(0\right)}\times \right\}\left[1-c\left(q_l\cdot {\Delta }_l\right)\right]+{\overline{k}}_l^{\left(0\right)}\cdot {\overline{k}}_l^{\left(0\right)T}\cdot \left[q_l-s\left(q_l\cdot {\Delta }_l\right)\right]\right\}\cdot {\overline{v}}_l^{\left(0\right)}$$

(80)

$$\mathrm{where}{\overline{v}}_l^{\left(0\right)}=\left\{{\overline{p}}_l\times \right\}{\overline{k}}_l^{\left(0\right)}\cdot {\Delta }_l+\left(1-{\Delta }_l\right)\cdot {\overline{k}}_l^{\left(0\right)}$$

(81)

${\overline{v}}_l^{\left(0\right)}$ - is a screw parameter or one of the homogeneous coordinates of the driving axis.

The first- and second- order differential matrices, previously introduced as reference expressions, play an essential in establishing dynamics matrices.

These dynamics matrices are key components in the matrix formulations that define higher – order acceleration energies.

3.2. The Matrix Expressions for the Acceleration Energies

The following relations are the starting point for the formulation of acceleration energy:

$$\begin{array}{c}{E}_A^{\left(p\right)}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\dots ;\stackrel{\left(p+1\right)}{\overline{\theta }\left(t\right)}\right]\\ ={\displaystyle \frac{1}{2}}\cdot {\displaystyle \sum _{i=1}^{n}Trace \left\{\stackrel{\left(p+1\right)}{{}_i{}^0\left[R\right]}\cdot \left[{}^{i}I_{pi}^{\ast }+M_i\cdot {}^i\overline{r}_{C_i}\cdot {}^i\overline{r}_{C_i}^T\right]\cdot \stackrel{\left(p+1\right)}{{}^i{}_0\left[R\right]_T}\right\}} +{\displaystyle \frac{1}{2}}\cdot {\displaystyle \sum _{i=1}^{n}Trace \left[\stackrel{\left(p+1\right)}{{\overline{p}}_i}\cdot \stackrel{\left(p+1\right)}{{\overline{p}}_i^T}\right]\cdot M_i}\\ ={\displaystyle \frac{1}{2}}\cdot {\displaystyle \sum _{i=1}^{n}Trace \left\{\stackrel{\left(p+1\right)}{{}_i{}^0\left[R\right]}\cdot \left[{\displaystyle \int {}^i\overline{r}_i^{\ast }\cdot {}^i\overline{r}_i^{\ast T}\cdot dm}\right.\right.}\left.+\left.{}^i\overline{r}_{C_i},\cdot ,{}^i\overline{r}_{C_i}^T,\cdot ,{\displaystyle \int dm}\right],\cdot ,\stackrel{\left(p+1\right)}{{}^i{}_0\left[R\right]_T}\right\}+{\displaystyle \frac{1}{2}}\cdot {\displaystyle \sum _{i=1}^{n}Trace \left[\stackrel{\left(p+1\right)}{{\overline{p}}_i}\cdot \stackrel{\left(p+1\right)}{{\overline{p}}_i^T}\right]\cdot {\displaystyle \int dm}}\end{array}$$

(82)

Equation (82) defines the acceleration energy of order $"p=1,2,3,\mathrm{...}"$ for the entire mechanical system. The notations used in this equation have the following meaning:

$(p)$ and $(p+1)$ represent the order of the absolute time derivatives,${}^i\overline{r}_i^{\ast }$is the position vector of the elementary mass $dm$, relative to a reference frame $\left\{i^{\ast }\right\}$ located at the mass center; ${}^{\left(i\right)}\overline{r}_{C_i}$is the position vector of the mass center, projected onto the fixed$\left\{0\right\}$ or moving $\left\{i\right\}$ frame; ${}_i{}^0\left[R\right]$ is the orientation matrix between the two reference systems.

The mechanical system is characterized by$\text{"}\mathrm{n}\text{"}$degrees of freedom (d.o.f.) or generalized coordinates, which are included in the column matrix $\overline{\theta }(t)={(q_i(t),fori=1\to n)}^T$.

It should be noted that the explicit expressions for the higher – order acceleration energies were previously determined in [4], using mass integrals.

In this section, only the matrix form of the higher orders acceleration energies is presented.

An essential aspect is that the dynamics matrices must be included in these expressions. These matrices comprise the differential matrices from the advanced kinematics, expressed through matrix exponentials [26] and developed in [5,18].

3.3. TheAcceleration Energies of the First Order

Starting from Equation (82) and applying a series of successive matrix transformations, the defining expression for the first – order acceleration energy of is obtained. This expression can be represented in the following matrix form:

$$\begin{array}{c}{E}_A^{(1)}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\ddot{\overline{\theta }}\left(t\right)\right]={\displaystyle \frac{1}{2}}\cdot \left\{{\ddot{\overline{\theta }}}^T\left(t\right)\cdot M\left[\overline{\theta }\left(t\right)\right]\cdot \ddot{\overline{\theta }}\left(t\right)\right.\\ \left.+{\ddot{\overline{\theta }}}^T\left(t\right)\cdot V\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]+\left[{\dot{\overline{\theta }}}^T\left(t\right)\cdot D\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]\cdot \dot{\overline{\theta }}\left(t\right)\right]\right\}. \end{array}$$

(83)

Equation (83) includes the column vector of the generalized velocities and accelerations, whose components correspond to the first- and second – order time derivatives of the column matrix of generalized variables $\left\{\overline{\theta }\left(t\right), \dot{\overline{\theta }}\left(t\right), \ddot{\overline{\theta }}\left(t\right)\right\}$.

Additionally, within Equation (83) a set of dynamics matrices can be identified.

These matrices are defined as follows:

$$\underset{\left(n\times n\right)}{M\left[\overline{\theta }\left(t\right)\right]}=\left[M_{ij}={\displaystyle \sum _{k=\max \left(i;j\right)}^{n}Trace\left[A_{ki}\cdot {}^{k}I_{psk}\cdot A_{kj}^T\right]} , i=1\to n, j=1\to n\right],$$

(84)

$$\underset{\left(n\times 1\right)}{V\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]}=\left[V_i\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right],i=1\to n\right].$$

(85)

Equation (84) defines a $\left(n\times n\right)$ matrix, commonly referred to as the mass matrix or the inertia matrix associated with the first – order acceleration energy. Its components denoted by $M_{ij}$ are explicitly defined in the same context. The column matrix of centrifugal and Coriolis terms of the first-order is defined by Equation(86), and is represented as an$\left(n\times 1\right)$ matrix, with the following components:

$$V_i\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]={\dot{\overline{\theta }}}^T\cdot \left[V_{ijm}={\displaystyle \sum _{k=\max \left(j;m\right)}^{n}Trace\left[A_{ki}\cdot {}^{k}I_{psk}\cdot A_{kjm}^T\right]} \right]\cdot \dot{\overline{\theta }}.\mathrm{where}j=1\to n, m=1\to n$$

(86)

The pseudo inertial matrix corresponding to the acceleration energies of first order is defined as follows:

$$D\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]={\dot{\overline{\theta }}}^T\cdot \left[D_{ijlm}={\displaystyle \sum _{k=\max \left(i;j;l;m\right)}^{n}Tr\left[A_{kij}\cdot {}^{k}I_{psk}\cdot A_{klm}^T\right]}\right]\cdot \dot{\overline{\theta }}\mathrm{where}\left\{i=1\to n, j=1\to n, l=1\to n, m=1\to n\right\}$$

(87)

The formulations in Equations (83) and (87) include the differential matrices of first and second-order denoted with $A_{ki}$ and$A_{kjm}$, respectively. These are computed via the Equations (71)–(74) in the case of the differential matrices of first-order ($A_{ki}$), while for the differential matrices of second order, the Equations (76)–(79) are applied.

3.3. The Acceleration Energy of Second-Order

As demonstrated in the authors’ previous research [3,4,5], sudden motions, transient motion phases, and the mechanical systems subjected to the action of a system of external forces, with a time variation law, are characterized by higher order linear and angular accelerations.

More specifically, fast – moving mechanical systems, which are subjected to the action of external forces with time variation law, are characterized by higher order linear and angular accelerations.

Unlike the formulations presented in [4], where the second – order acceleration energy of the was established using mass integrals, the following expression introduces a matrix – based formulation of the second – order acceleration energy, as:

$$\begin{array}{c}{E}_A^{\left(2\right)}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\ddot{\overline{\theta }}\left(t\right);\stackrel{⃛}{\overline{\theta }}\left(t\right)\right]\\ ={\displaystyle \frac{1}{2}}\cdot {\stackrel{⃛}{\overline{\theta }}}^T\left(t\right)\cdot M\left[\overline{\theta }\left(t\right)\right]\cdot \stackrel{⃛}{\overline{\theta }}\left(t\right)+3\cdot {\stackrel{⃛}{\overline{\theta }}}^T\left(t\right)\cdot V\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\ddot{\overline{\theta }}\left(t\right)\right]\\ +\mathrm{...}+\left\{{\dot{\overline{\theta }}}^T\left(t\right)\cdot \left[{\dot{\overline{\theta }}}^T\left(t\right)\cdot D\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]\cdot \dot{\overline{\theta }}\left(t\right)\right]\cdot \dot{\overline{\theta }}\left(t\right)\right\}\end{array}$$

(88)

The components $E_A^{\left(2\right)}\left[\overline{\theta }\left(t\right); \dot{\overline{\theta }}\left(t\right); {\ddot{\overline{\theta }}}^2\left(t\right)\right]$and $E_A^{\left(2\right)}\left[\overline{\theta }\left(t\right); {\dot{\overline{\theta }}}^6\left(t\right); \ddot{\overline{\theta }}\left(t\right)\right]$are not developed in this paper. Their components in explicit form can be found in the research presented in [4]. Equation (88) includes the dynamics matrices of the second order. These are defined as:

$$V\left[\overline{\theta }(t);\dot{\overline{\theta }}(t);\ddot{\overline{\theta }}(t)\right]=\left[V_i^{*}\left[\overline{\theta }(t);\dot{\overline{\theta }}(t);\ddot{\overline{\theta }}(t)\right], \mathrm{where} i=1\to n\right],\mathrm{and} \underset{\left(n\times 1\right)}{V_i^{*}}\left[\overline{\theta }(t);\dot{\overline{\theta }}(t);\ddot{\overline{\theta }}(t)\right]={\ddot{\overline{\theta }}}^T\cdot \left[\left\{V_{ijm}=V_{imj}\right\},\mathrm{where} j=1\to n \mathrm{and} m=1\to n\right]\cdot \dot{\overline{\theta }},$$

(89)

Therefore, both the inertia matrix and the matrix of Coriolis and centrifugal terms, are components in the expression of the second – order acceleration energy.

3.4. The Acceleration Energy of Third-Order

The dynamic analysis is further extended to include the third – order acceleration energy. As previously stated, the explicit integral formulation of the third – order acceleration energy was presented in [4]. According to the approach introduced in [3], the matrix form of the third – order, defined as a function which depends exclusively on $\stackrel{⃜}{\overline{\theta }}={({\stackrel{⃜}{q}}_{ i}, i=1\to n)}^T$, is given as follows:

$$\begin{array}{c}{E}_A^{\left(3\right)}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\ddot{\overline{\theta }}\left(t\right);\stackrel{⃛}{\overline{\theta }}\left(t\right);\stackrel{⃜}{\overline{\theta }}\left(t\right)\right]\\ ={\displaystyle \frac{1}{2}}\cdot {\stackrel{⃜}{\overline{\theta }}}^T\left(t\right)\cdot M\left[\overline{\theta }\left(t\right)\right]\cdot \stackrel{⃜}{\overline{\theta }}\left(t\right)+4\cdot {\stackrel{⃜}{\overline{\theta }}}^T\left(t\right)\cdot V\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\stackrel{⃛}{\overline{\theta }}\left(t\right)\right]+\\ +3\cdot {\stackrel{⃜}{\overline{\theta }}}^T\left(t\right)\cdot V^{\ast }\left[\overline{\theta }\left(t\right);{\ddot{\overline{\theta }}}^2\left(t\right)\right]\\ +\mathrm{...}+\left\{{\dot{\overline{\theta }}}^T\left(t\right)\cdot \left\{{\dot{\overline{\theta }}}^T\left(t\right)\cdot \left[{\dot{\overline{\theta }}}^T\left(t\right)\cdot D\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]\cdot \dot{\overline{\theta }}\left(t\right)\right]\cdot \dot{\overline{\theta }}\left(t\right)\right\}\cdot \dot{\overline{\theta }}\left(t\right)\right\}.\end{array}$$

(90)

The dynamics matrices of the third order, included in the acceleration energy of third order, are defined by the following expressions:

$$\begin{array}{c}V\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\stackrel{⃛}{\overline{\theta }}\left(t\right)\right]\\ =\underset{\left(n\times 1\right)}{Matrix}{\left\{ {\stackrel{⃛}{\overline{\theta }}}^T\cdot \left[V_{ijm}, \mathrm{where} i=1\to n, j=1\to n, m=1\to n\right]\cdot \dot{\overline{\theta }} \right\}}^T\end{array}$$

(91)

$V^{\ast }\left[\overline{\theta }\left(t\right);{\ddot{\overline{\theta }}}^2\left(t\right)\right]=\underset{\left(n\times 1\right)}{Matrix}{\left\{ {\ddot{\overline{\theta }}}^T\cdot \left[V_{ijm}, \mathrm{where} i=1\to n, j=1\to n, m=1\to n\right]\cdot \ddot{\overline{\theta }} \right\}}^T$It is important to note that, in case of third – order acceleration energy, the mass properties are characterized by the pseudo – inertial tensor and higher-order differential matrices.

3.5. The Acceleration Energy of p-th-Order

For multibody mechanical systems, the fourth-order acceleration energy becomes a central function $\stackrel{\left(5\right)}{\overline{\theta }}\left(t\right)$, in the dynamic formulation. The corresponding expression is given as follows:

$$\begin{array}{l}{E}_A^{\left(4\right)}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\ddot{\overline{\theta }}\left(t\right);\stackrel{⃛}{\overline{\theta }}\left(t\right);\stackrel{⃜}{\overline{\theta }}\left(t\right);\stackrel{\left(5\right)}{\overline{\theta }}\left(t\right)\right]={\displaystyle \frac{1}{2}}\cdot \stackrel{\left(5\right)}{{\overline{\theta }}^T}\left(t\right)\cdot \left(t\right)\cdot M\left[\overline{\theta }\left(t\right)\right]\cdot \stackrel{\left(5\right)}{\overline{\theta }}\left(t\right)+\mathrm{...}\\ +\mathrm{...}+\left\{{\dot{\overline{\theta }}}^T\left(t\right)\cdot \left\{{\dot{\overline{\theta }}}^T\left(t\right)\cdot \left\{{\dot{\overline{\theta }}}^T\left(t\right)\cdot \left[{\dot{\overline{\theta }}}^T\left(t\right)\cdot D\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]\cdot \dot{\overline{\theta }}\left(t\right)\right]\cdot \dot{\overline{\theta }}\left(t\right)\right\}\cdot \dot{\overline{\theta }}\left(t\right)\right\}\cdot \dot{\overline{\theta }}\left(t\right)\right\}.\end{array}$$

(92)

The acceleration energy of pth order can be established According to the models outlined in [5,18], its corresponding expression is given below:

$$\begin{array}{c}{E}_A^{\left(p\right)}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\ddot{\overline{\theta }}\left(t\right);\stackrel{⃛}{\overline{\theta }}\left(t\right);\stackrel{⃜}{\overline{\theta }}\left(t\right);\mathrm{...};\stackrel{\left(p\right)}{\overline{\theta }}\left(t\right)\right]={\displaystyle \frac{1}{2}}\cdot \stackrel{\left(p\right)}{{\overline{\theta }}^T}\left(t\right)\cdot \left(t\right)\cdot M\left[\overline{\theta }\left(t\right)\right]\cdot \stackrel{\left(p\right)}{\overline{\theta }}\left(t\right)+\mathrm{...}\\ \mathrm{...}+\left\{{\displaystyle \prod _{j=1}^p{{\dot{\overline{\theta }}}_j}^T\left(t\right)\cdot \left[{\dot{\overline{\theta }}}^T\left(t\right)\cdot D\left[\overline{\theta }\left(t\right);{\dot{\overline{\theta }}}^2\left(t\right)\right]\cdot \dot{\overline{\theta }}\left(t\right)\right]\cdot {\dot{\overline{\theta }}}_j\left(t\right) }\right\}.\end{array}$$

(93)

In the expression of higher-order acceleration energy, the contribution of mass properties is reflected through the presence pseudo – inertial tensor. The tensor is a key component of the dynamics matrices, which were previously introduced and defined in Equations (84), (87) and (89).

4. The Generalized Forces in the Dynamics of Systems

In accordance with differential principles typical to analytical dynamics of systems, the study of a system’s dynamical behavior fundamentally relies on the concept of generalized forces. These forces are mechanically developed in direct relation to the generalized variables, also referred to as independent parameters or degrees of freedom (d.o.f.) which uniquely describe the absolute motion of holonomic mechanical systems. From a mechanical perspective, the generalized forces are due to various sources including actuation mechanisms, gravitational fields, manipulating loads, as well as complex frictional effects associated with physical connections between kinetic subassemblies in a MBS – as illustrated in Figure 4 and Figure 5.

According to [3,4,18], each kinetic ensemble $\left(i=1\to n\right)$, that is a part of the mechanical structure of a robot, and implicitly integrated within the multibody system (MBS), is subject to a system of external and active forces, manipulating loads, and complex frictional effects, as seen in Figure 4.

Depending on the behavior (whether static or dynamic) of each physical link – particularly within driving joints of fifth order – a corresponding generalized force is generated. This may take the form of a generalized static force or a generalized actuation force, depending on the context of motion and interaction.

4.1. The Generalized Static Forces

To begin with and by following the classical algorithm described in [4,5], the generalized forces corresponding to static equilibrium are determined. Based on the fundamental theorems of statics applied to mechanical systems, and after performing a few transformations, the following expressions are obtained for the constraint force ${}^i\overline{f}_i$ and for the corresponding moment of the constraint forces ${}^i\overline{n}_i$:

$${}^i\overline{f}_i={\Delta }_{ m}^{ 2}\cdot {\displaystyle \sum _{j=i}^n{M}_j}\cdot {}^i{}_0\left[R\right]^T\cdot \overline{g}+{\left(-1\right)}^{{\Delta }_{ M}}\cdot {\displaystyle \frac{1-{\Delta }_m}{1+{\Delta }_m}}\cdot {}_{n+1}{}^i\left[R\right] \cdot {}^{n+1}\overline{f}_{n+1}$$

(94)

$$\begin{array}{c}{}^i\overline{n}_i={\Delta }_{ m}^{ 2}\cdot {\displaystyle \sum _{j=i}^n{M}_j}\cdot {}^i{}_0\left[R\right]^T\cdot \left[\left({}^0\overline{r}_{C_j}-{\overline{p}}_i\right) \times \overline{g}\right]\\ +{\left(-1\right)}^{{\Delta }_{ m}}\cdot {\displaystyle \frac{1-{\Delta }_{ m}}{1+{\Delta }_{ m}}}\cdot \left\{{}^i{}_0\left[R\right]^T\cdot \left(\overline{p}-{\overline{p}}_i\right)\times {}_{n+1}{}^i\left[R\right] \cdot {}^{n+1}\overline{f}_{n+1}\right.+\left.{}_{n+1}{}^i\left[R\right] \cdot {}^{n+1}\overline{n}_{n+1}\right\}\end{array}$$

(95)

$$\mathrm{where}{\Delta }_{ m}=\left\{\left\{-1;\left(SU; M_i\right)\right\};\left\{0;SU\right\};\left\{1; M_i\right\} \right\}$$

(96)

$$\overline{g}=\tau \cdot g\cdot {\overline{k}}_0, \mathrm{and} {\overline{k}}_0=\left\{{\overline{x}}_0;{\overline{y}}_0;{\overline{z}}_0\right\}\in \left\{0\right\}$$

(97)

$$\mathrm{and}\tau =-{\overline{k}}_g^T\cdot {\overline{k}}_0=\left\{\left(-1, \mathrm{if} {\overline{k}}_g^T\cdot {\overline{k}}_0=1\right), \left(1, \mathrm{if} {\overline{k}}_g^T\cdot {\overline{k}}_0=-1\right)\right\}, {\overline{k}}_g={\displaystyle \frac{{}_0\overline{g}}{\left|{}_0\overline{g}\right|}},$$

(98)

where ${\overline{k}}_g$ is the unit vector of ${}_0\overline{g}$ (the gravity acceleration vector).

The operator defined in Equation (96) serves to highlight two types of external loads: gravitational loads denoted by$\left(M_i\right)$ and manipulating loads represented by symbol$\left(SU\right)$. The equations (94) and (95) are established through outward iterations $\left(i=1\to n\right)$, following the model described in [4,5]. These expressions are then substituted in the formulation of the generalized static force, yielding to the following result:

$$Q_S^i=\left\{{}^i\overline{f}_i^T\cdot \left(1-{\Delta }_{ i}\right)+{}^i\overline{n}_i^T\cdot {\Delta }_{ i}\right\} \cdot {}^i\overline{k}_i=Q_g^i+Q_{SU}^i$$

(99)

where the quantities$Q_g^i$ and $Q_{SU}^i$ represent the generalized gravitational and manipulating forces, respectively, as defined in [25]. They are commonly referred to as generalized active forces. Unlike the classical formulation based on the principle of virtual work, the approach adopted here – as detailed in [5,18,25] – determines these forces using transfer matrices for linear and angular velocities. Based on Equations (94), (95) and (99), the starting equation for the generalized gravitational force is given by:

$$\begin{array}{c}{Q}_g^i={\displaystyle \sum _{j=i}^n{M}_j}\cdot \left\{{\overline{g}}^T\cdot {}_i{}^0\left[R\right]\cdot \left(1-{\Delta }_{ i}\right)+\right.\left.{\overline{g}}^T\cdot {\left[ \left({}^0\overline{r}_{C_j}-{\overline{p}}_i\right) \times \right]}^T\cdot {}_i{}^0\left[R\right]\cdot {\Delta }_{ i}\right\} \cdot {}^i\overline{k}_i\\ ={\displaystyle \sum _{j=i}^n{M}_j}\cdot {\overline{g}}^T\cdot \left\{{}^0\overline{k}_i\cdot \left(1-{\Delta }_{ i}\right)+\right.{}^0\overline{k}_i\times \left({}^0\overline{r}_{C_j}-{\overline{p}}_i\right)\left.\cdot {\Delta }_{ i}\right\}, \mathrm{where} {}_i{}^0\left[R\right]\cdot {}^i\overline{k}_i={}^0\overline{k}_i\end{array}$$

(100)

Considering the transfer matrices for linear and angular velocities, as detailed in [25], a set of transformations is carried out to determine the generalized gravitational force. These transformations are presented as follows:

$${}^0\overline{k}_i\cdot \left(1-{\Delta }_{ i}\right)+{}^0\overline{k}_i\times \left({\overline{p}}_n-{\overline{p}}_i\right)\cdot {\Delta }_{ i}={\displaystyle \frac{\partial {\overline{p}}_n}{\partial q_i}}={\overline{V}}_i\subset {}^{0}J_i\subset {}_{0}J\left[\overline{\theta }\left(t\right)\right]$$

(101)

$${\overline{g}}^T\cdot \left[{}^0\overline{k}_i\times \left({}^0\overline{r}_{C_j}-{\overline{p}}_n\right)\right]\cdot {\Delta }_{ i}={}^0\overline{k}_i^T\cdot {\Delta }_{ i}\cdot \left[\left({}^0\overline{r}_{C_j}-{\overline{p}}_n\right)\times \overline{g}\right]$$

(102)

$${}^0\overline{k}_i^T\cdot {\Delta }_{ i}= {\overline{\Omega }}_i\subset {}^{0}J_i^T\subset {}^{0}J^T\left[\overline{\theta }\left(t\right)\right]$$

(103)

The column vectors (101) – (103) are components of the Jacobian matrix, according with:

$${}_{0}J\left[\overline{\theta }\left(t\right)\right] = \left[{}^{0}J_i\left[{\overline{\theta }}_{i0}\left(t\right)\right] =\left[ \begin{array}{c}{V}_i\\ {\Omega }_i\end{array}\right]=\left[ \begin{array}{c}{}^0\overline{d}_i\left[{\overline{\theta }}_{i0}\left(t\right)\right]\\ {}^0\overline{\delta }_i\left[{\overline{\theta }}_{i0}\left(t\right)\right]\end{array}\right] , i=1\to n\right]$$

(104)

where

$${}^0\overline{d}_i = {\overline{V}}_i = {}^0\overline{k}_i\times \left({\overline{p}}_n-{\overline{p}}_i\right)\cdot {\Delta }_i + {}^0\overline{k}_i\cdot \left(1-{\Delta }_i\right) = {\displaystyle \frac{\partial {\overline{p}}_n}{\partial q_i}} = A_{ni}\left(\overline{p}\right)$$

(105)

$${}^0\overline{\delta }_i = {\overline{\Omega }}_i = {}^0\overline{k}_i\cdot {\Delta }_i = {}_i{}^0\left[R\right]\cdot {}^i\overline{k}_i\cdot {\Delta }_i$$

(106)

The same components of the Jacobian matrix can alternatively be determined using matrix exponentials, as detailed in [25,26]. The corresponding expressions are provided below:

$$\begin{array}{c}{\overline{V}}_i={\displaystyle \frac{\partial }{\partial q_i}}\left({\overline{p}}_n\right)=A_{ni}\left(\overline{p}\right)={\displaystyle \prod _{j=0}^{i-1}\exp \left\{\left\{ {\overline{k}}_j^{\left(0\right)}\times \right\} q_j \cdot {\Delta }_j\right\}}\cdot {\overline{v}}_i^{\left(0\right)} +\\ +{\displaystyle \prod _{k=i}^n\exp \left\{\left\{ {\overline{k}}_k^{\left(0\right)}\times \right\} q_k \cdot {\Delta }_k\right\}}\cdot {\overline{p}}_n^{\left(0\right)}+\\ +{\Delta }_i \cdot {\displaystyle \prod _{j=0}^{i-1}\exp \left\{\left\{ {\overline{k}}_j^{\left(0\right)}\times \right\} q_j \cdot {\Delta }_j\right\}}\cdot \left\{{\overline{k}}_i^{\left(0\right)}\times \right\}\cdot \left\{{\displaystyle \sum _{k=i}^n\left\{{\displaystyle \prod _{m=i-1}^{k-1}\exp \left\{\left\{ {\overline{k}}_m^{\left(0\right)}\times \right\} q_m \cdot {\delta }_m\right\}}\right\}}\right. \cdot {\overline{b}}_k\end{array}$$

(107)

$$\underset{\left(3\times 1\right)}{{\overline{\Omega }}_i} ={\displaystyle \prod _{j=0}^{i-1} \left\{ \exp \left\{\left\{ {\overline{k}}_j^{\left(0\right)}\times \right\} q_j\cdot {\Delta }_j\right\} \right\}}\cdot {\overline{k}}_i^{\left(0\right)}$$

(108)

Parameters from Equations (107) and (108) maintain the same definitions provided in § 3.1.

Based on these, the defining formulation for the generalized gravitational force results:

$$Q_g^i={\displaystyle \sum _{j=i}^n{M}_j}\cdot {\overline{g}}^T\cdot {\displaystyle \frac{\partial {}^0\overline{r}_{C_j}}{\partial q_i}}={}^{\left(n\right)0}J_i^T\cdot {}^{\left(n\right)0}\overline{\mathcal{F}}_{X_i}$$

(109)

$$\mathrm{where}\underset{\left(6\times 1\right)}{{}^{\left(\mathcal{n}\right)0}\overline{\mathcal{F}}_{{\mathcal{X}}_{\mathcal{i}}}}={\left[\begin{array}{cc}{}^{\left(n\right)0}\overline{F}_{X_i}^T& {}^{\left(n\right)0}\overline{N}_{X_i}^T\end{array}\right]}^T=\left[\begin{array}{c}\left(1-{\Delta }_{ i}\right)\cdot {\displaystyle \sum _{j=i}^n{M}_j}\cdot {}^n{}_{\left(0\right)n}\left[R\right]^T\cdot \overline{g}\\ \mathrm{................................................................}\\ {\Delta }_{ i}\cdot {\displaystyle \sum _{j=i}^n{M}_j}\cdot {}^n{}_{\left(0\right)n}\left[R\right]^T\cdot \left[\left({}^0\overline{r}_{C_j}-{\overline{p}}_n\right) \times \overline{g}\right]\end{array}\right]$$

(110)

$$\mathrm{and}{Q}_g\left(\overline{\theta }\right)={\left[Q_g^i={}^{\left(n\right)0}J_i^T\cdot {}^{\left(n\right)0}\overline{\mathcal{F}}_{X_i}; i=1\to n\right]}^{ T}$$

(111)

The column vector defined in Equation (110), expressed in Cartesian space, is mechanically equivalent to the wrench (i.e., the combined resultant force and moment) generated by the gravitational forces acting on the interval $\left[i , n\right]$ relative to the $\left\{n\right\}$ moving reference frame, attached to the geometric center of the final driving joint from MBS, as illustrated in Figure 4. Following the expressions (94), (95) and (99), the starting equation for the generalized manipulating force is given as:

$$\begin{array}{c}{Q}_{SU}^i={\left\{{}_{n+1}{}^i\left[R\right] \cdot {}^{n+1}\overline{f}_{n+1}\right\}}^T\cdot \left(1-{\Delta }_{ i}\right) \cdot {}^i\overline{k}_i\\ +\left\{{}^i{}_0\left[R\right]^T\cdot \left(\overline{p}-{\overline{p}}_i\right)\times {}_{n+1}{}^i\left[R\right] \cdot {}^{n+1}\overline{f}_{n+1}\right.{\left.+{}_{n+1}{}^i\left[R\right] \cdot {}^{n+1}\overline{n}_{n+1}\right\}}^T \cdot {}^i\overline{k}_i\cdot {\Delta }_{ i}\end{array}$$

(112)

$$\begin{array}{c}{Q}_{SU}^i={\left\{{}_{n+1}{}^0\left[R\right]\cdot {}^{n+1}\overline{f}_{n+1}\right\}}^T\cdot \left\{{}^0\overline{k}_i\cdot \left(1-{\Delta }_{ i}\right)+\right.\left.{}^0\overline{k}_i\times \left({\overline{p}}_n-{\overline{p}}_i\right)\cdot {\Delta }_{ i}\right\}+\\ +{\left\{{}_{n+1}{}^0\left[R\right]\cdot {}^{n+1}\overline{f}_{n+1}\right\}}^T\cdot \left[{}^0\overline{k}_i\times \left(\overline{p}-{\overline{p}}_n\right)\right]+{\left\{{}_{n+1}{}^0\left[R\right]\cdot {}^{n+1}\overline{n}_{n+1}\right\}}^T\cdot {}^0\overline{k}_i\cdot {\Delta }_{ i}\end{array}$$

(113)

$$Q_{SU}^i={}^{\left(n\right)0}J_i^T\cdot {}^{\left(n\right)0}\overline{\mathcal{F}}_X$$

(114)

$${}_0\underset{\left(6\times 1\right)}{{\overline{\mathcal{F}}}_X}={\left[\begin{array}{cc}{}^0\overline{F}_X^T& {}^0\overline{N}_X^T\end{array}\right]}^T=\left[\begin{array}{cc}{}_{n+1}{}^{\left(n\right)0}\left[R\right]& \left[0\right]\\ \left\{{}^n{}_{\left(0\right)n}\left[R\right]^T\cdot {\overline{p}}_{n+1n}\times \right\}{}_{n+1}{}^n\left[R\right]& {}_{n+1}{}^{\left(0\right)n}\left[R\right]\end{array}\right]\cdot \left[\begin{array}{c}{}^{n+1}\overline{f}_{n+1}\\ {}^{n+1}\overline{n}_{n+1}\end{array}\right]$$

(115)

$$Q_{SU}\left(\overline{\theta }\right)={}_{\left(n\right)0}J\left(\overline{\theta }\right) \cdot {}^{\left(n\right)0}\overline{\mathcal{F}}_X={\left[Q_{SU}^i={}^{\left(n\right)0}J_i^T\cdot {}^{\left(n\right)0}\overline{\mathcal{F}}_X; i=1\to n\right]}^T$$

(116)

The Cartesian column vector defined in Equation (115) is mechanically equivalent, as shown in [5,25] to the wrench associated with the manipulating load, with respect to the$\left\{n\right\}$ reference frame.

4.2. The Generalized Dynamics Forces

Considering the aspects presented in [25], in the context of the dynamical behavior of multibody systems (MBS), each driving joint is subjected not only to the active forces, discussed in the previous subsection, but also to the effects of generalized inertia and driving forces. The determination of these forces constitutes the main objective of the current subsection. As a starting point, the iterative algorithm associated to the classical formulation of the dynamic equations, is applied. The traditional approach is based on the D’Alembert principle of analytical dynamics which accounts for inertia by transforming a dynamic system into a quasi-static one through the inclusion of inertial forces.

According to [5,6,7,8] and similar to (94) and (95), the dynamical actions are rewritten as:

$$\begin{array}{c}{}^i\overline{f}_i={\Delta }_m^2\cdot \left\{{\displaystyle \sum _{j=i}^n{}_j{}^i\left[R\right]\cdot {}^j\overline{F}_j^{\ast }}+{\displaystyle \sum _{j=i}^n{M}_j\cdot {}^j{}_0\left[R\right]^T}\cdot \overline{g}\right\}\\ +{\left(-1\right)}^{{\Delta }_m}\cdot {\displaystyle \frac{1-{\Delta }_{ m}}{1+3\cdot {\Delta }_{ m}}}\cdot {}_{n+1}{}^i\left[R\right]\cdot {}^{n+1}\overline{f}_{n+1}\end{array}$$

(117)

$$\begin{array}{c}{}^i\overline{n}_i={\Delta }_m^2\cdot \left\{{\Delta }_{\theta }\cdot \left\{{\displaystyle \sum _{j=i}^n{}^i{}_0\left[R\right]^T\cdot \left({}^0\overline{r}_{C_j}-{\overline{p}}_i\right)\times {}_j{}^i\left[R\right]\cdot {}^j\overline{F}_j^{\ast }}\right.\right.\left.+{}_j{}^i\left[R\right]\cdot {}^j\overline{N}_j^{\ast }\right\}\left.+{\displaystyle \sum _{j=i}^n{M}_j}\cdot {}^i{}_0\left[R\right]^T\cdot \left[\left({}^0\overline{r}_{C_j}-{\overline{p}}_i\right) \times \overline{g}\right]\right\}+\\ {\left(-1\right)}^{{\Delta }_{ m}}\cdot {\displaystyle \frac{1-{\Delta }_m}{1+{\Delta }_m}}\cdot \left\{{}^i{}_0\left[R\right]^T\cdot \left(\overline{p}-{\overline{p}}_i\right)\times {}_{n+1}{}^i\left[R\right] \cdot {}^{n+1}\overline{f}_{n+1}\right.\left.+{}_{n+1}{}^i\left[R\right] \cdot {}^{n+1}\overline{n}_{n+1}\right\}\end{array}$$

(118)

In Equations (117) and (118) one can observe the application of two fundamental theorems in dynamics of mechanical systems: the theorem of the motion of the mass center, expressed in Equation (119) and theorem of the angular momentum formulated in Equation (120):

$${}^i\overline{F}_i^{\ast }={{\displaystyle M}}_i\cdot {}_i{\dot{\overline{v}}}_{C_i}$$

(119)

$${}^i\overline{N}_i^{\ast }={}^{i}I_i^{\ast }\cdot {}^i\dot{\overline{\omega }}_i+{}^i\overline{\omega }_i\times {}^{i}I_i^{\ast }\cdot {}^i\overline{\omega }_i$$

(120)

$${}_i{\dot{\overline{v}}}_{C_i}={}_i{\dot{\overline{v}}}_i+{}^i\dot{\overline{\omega }}_i\times {}^i\overline{r}_{C_i}+{}^i\overline{\omega }_i\times {}^i\overline{\omega }_i\times {}^i\overline{r}_{C_i}$$

(121)

The above expressions apply to each kinetic ensemble in the multibody system.

The involved symbols have the following physical interpretations: ${}^i\dot{\overline{v}}_{C_i}$ denotes the acceleration of mass center; ${}^i\overline{F}_i^{\ast }$ and ${}^i\overline{N}_i^{\ast }$ represent the resultant force and resultant moment, respectively, and ${}^{i}I_i^{\ast }$ is the axial and centrifugal inertia tensor defined relative to the frame$\left\{i^{\ast }\right\}$ which is applied at the mass center of the kinetic ensemble.

Starting from Equations (117) and (118), the following expression for the generalized driving force is obtained:

$$Q_m^i\left(t\right)={\Delta }_m^2\cdot \left[{\Delta }_{\theta }\cdot Q_{i\mathcal{F}}^i\left(t\right)+Q_g^i\left(t\right)\right]+{\left(-1\right)}^{{\Delta }_m}\cdot {\displaystyle \frac{1-{\Delta }_m}{1+3\cdot {\Delta }_m}}\cdot Q_{SU}^i\left(t\right)$$

(122)

where $Q_g^i\left(t\right)$ and $Q_{SU}^i\left(t\right)$ denote the generalized active forces, $Q_{i\mathcal{F}}^i\left(t\right)$ represents the generalized inertia force, and ${\Delta }_{\theta }$ is a scalar operator used to distinguish between static and dynamic operating conditions of the system, with ${\Delta }_{\theta }=0$ corresponding to static conditions and ${\Delta }_{\theta }=1$ corresponding to dynamic behavior.

$${\Delta }_{\theta }=\left\{\left[1; \left\{\dot{\overline{\theta }};\ddot{\overline{\theta }}\right\}\ne 0\right];\left[0; \left\{\dot{\overline{\theta }};\ddot{\overline{\theta }}\right\}=0\right]\right\}$$

(123)

The starting equation for determining the generalized inertia force $Q_{iö}^i\left(t\right)$is:

$$Q_{i\mathcal{F}}^i={\displaystyle \sum _{j=i}^n\left\{{\left\{{}_j{}^i\left[R\right]\cdot {}^j\overline{F}_j^{\ast }\right\}}^T\cdot \left(1-{\Delta }_{ i}\right)+\right.}{\Delta }_{ i}\cdot {\left\{{}^i{}_0\left[R\right]^T\cdot \left({}^0\overline{r}_{C_j}-{\overline{p}}_i\right) \times {}_j{}^i\left[R\right]\cdot {}^j\overline{F}_j^{\ast }\right\}}^T\left.+{\Delta }_{ i}\cdot {\left[{}_j{}^i\left[R\right]\cdot {}^j\overline{N}_j^{\ast }\right]}^T\right\}\cdot {}^i\overline{k}_i$$

(124)

Following the same approach as in Equation (113), the final expression of the generalized inertia force is written as follows:

$$Q_{i\mathcal{F}}^i ={}^{\left(n\right)0}J_i^T\cdot {}^{\left(n\right)0}\overline{\mathcal{F}}_{X_i}^{\ast }$$

(125)

$${}^0\overline{\mathcal{F}}_{X_i}^{\ast }={\left[\begin{array}{cc}{}^0\overline{F}_{X_i}^{*T}& {}^0\overline{N}_{X_i}^{*T}\end{array}\right]}^T=\left[\begin{array}{c}\left(1-{\Delta }_{ i}\right)\cdot {\displaystyle \sum _{j=i}^n{}_j{}^0\left[R\right]\cdot {}^j\overline{F}_j^{\ast }}\\ {\Delta }_{ i}\cdot {\displaystyle \sum _{j=i}^n\left\{\left({}^0\overline{r}_{C_j}-{\overline{p}}_n\right) \times {}_j{}^0\left[R\right]\cdot {}^j\overline{F}_j^{\ast }+{}_j{}^0\left[R\right]\cdot {}^j\overline{N}_j^{\ast }\right\}}\end{array}\right]$$

(126)

$$\underset{\left(n\times 1\right)}{Q_{i\mathcal{F}}}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\ddot{\overline{\theta }}\left(t\right)\right]={\left[Q_{i\mathcal{F}}^i\left(t\right)={}^{0}J_i^T\left(t\right)\cdot {}^0\overline{\mathcal{F}}_{X_i}^{\ast }\left(t\right) ; i=1\to n\right]}^{ T}$$

(127)

The column vector defined in Equation (110), expressed in Cartesian space, is mechanically equivalent to the wrench of the inertia forces acting on the interval $\left[i ;n\right]$, with respect to the $\left\{n\right\}$ moving frame attached to the geometric center of the last driving joint in the MBS structure (see Figure 4 and Figure 5).

Consequently, the generalized driving force corresponding to each driving axis in the multibody system (MBS) is finally expressed as follows:

$$\begin{array}{c}{Q}_m^i\left(t\right)\\ ={}^{\left(n\right)0}J_i^T\left(t\right)\cdot \left\{{\Delta }_{ m}^2\cdot \left.\left[{\Delta }_{ \theta }\cdot {}^{\left(n\right)0}\overline{\mathcal{F}}_{X_i}^{\ast }\left(t\right)+{}^{\left(n\right)0}\overline{\mathcal{F}}_{X_i}\left(t\right)\right]+{\left(-1\right)}^{{\Delta }_{ m}}\cdot {\displaystyle \frac{1-{\Delta }_{ m}}{1+3\cdot {\Delta }_{ m}}}\cdot {}^{\left(n\right)0}\overline{\mathcal{F}}_X\left(t\right)\right\}\right\}\end{array}$$

(128)

The generalized active and inertia forces share the same mathematical structure, which represents a significant advantage in formulating the dynamic’s equations (128) corresponding to each kinetic ensemble in the multibody system (MBS) (see Figure 5).

Considering the presence of complex friction forces acting within each driving joint of the robotic system, the corresponding generalized friction force, is defined as follows:

$$Q_{fd}^i=b_i\cdot {\dot{q}}_i+{\mu }_{iT}\cdot \left(1-{\Delta }_{ i}\right)\cdot \left|{}^i\overline{k}_i\times {}^i\overline{f}_i\right|\cdot \mathrm{sgn}\left({\dot{q}}_i \cdot {\Delta }_{\theta }\right)+{\mu }_{ i R}\cdot {\displaystyle \frac{d_i}{2}}\cdot {\Delta }_{ i}\cdot \left|{}^i\overline{k}_i\times {}^i\overline{f}_i\right|\cdot \mathrm{sgn}\left({\dot{q}}_i \cdot {\Delta }_{\theta }\right)$$

(129)

where the symbols: $b_i$ and ${\mu }_{i \left(T\right)R}$ represents the viscous friction coefficient and the dry friction coefficient, respectively both defined as a function of the joint type. Within Equation (129) the following operator is defined:

$${\Delta }_{ f}=\left\{\begin{array}{ccc}\left(-1, {\Delta }_{ m}=-1\right);& \left[0, {\Delta }_{ m}=\left(0, 1\right)\right];& \left(1, Q_{fd}^i\right)\end{array}\right\}$$

(130)

where ${\Delta }_{ f}$ highlights the contribution of the external loads by ${\Delta }_{ m}$ as well as the influence of the complex frictions.

In accordance with the mathematical aspects, previously discussed, the resultant force vector included in the generalized friction force (129) is determined as follows:

$${}^i\overline{f}_i={}^i{}_0\left[R\right]^T\cdot \left\{{}^0\overline{F}_{X_i}+{}^0\overline{F}_{X_i}^{\ast }+{}^0\overline{F}_X\right\}\subset Q_{fd}^i$$

(131)

Consequently, Equation (122) is reformulated as follows:

$$Q_{mf}^i\left(t\right)={\left(-1\right)}^{{\Delta }_{ f}}\cdot {\displaystyle \frac{1-{\Delta }_f}{1+3\cdot {\Delta }_f}}\cdot Q_m^i\left(t\right)+{\Delta }_{ f}^2\cdot Q_{fd}^i \left(t\right)$$

(132)

As a result, Equation (132) for $\left(i=1\to n\right)$ defines the system of $\left(n\right)$ generalized driving forces which are identical to the dynamic equations of the multibody system (MBS). These equations include the generalized active and inertia forces, as well as the complex friction effects acting within each joint. It can also evident that, in both cases – generalized gravitational forces $Q_g\left(\overline{\theta }\right)$,$Q_{SU}\left(\overline{\theta }\right)$ and generalized inertial forces $Q_{i\mathcal{F}}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\ddot{\overline{\theta }}\left(t\right)\right]$ – the influence of mass properties play an important role in shaping the dynamic behavior of the robot’s mechanical structure.

Furthermore, Equations (109), (114) and (128) emphasize one of the fundamental properties in analytical dynamics, namely the one related to the Jacobian matrix.

This matrix ensures, from both a mathematical and physical perspectives, the transfer of the gravitational loads, manipulation forces, and inertial wrenches – originally applied at the end-effector — to each driving axis of the robot.

5. The Higher Order Dynamic’s Equations

When mechanical systems (MBS) are subject to sudden or transitory motions, both theoretical and experimental research [25] confirms the existence and relevance of higher – order accelerations energy. These energies are incorporated into the higher – order dynamic equations, where the time variation of generalized forces becomes explicit. In the specific case of the robotic mechanical structure (see Figure 5 ) which is characterized by the sudden motions, both higher – order generalized accelerations, and the corresponding generalized forces are developed. Based on the studies in references [25,26,27,28,29,30], the following categories of higher – order derivatives are identified:

$$\stackrel{\left(k\right)}{Q_{i\mathcal{F}}^i\left(t\right)}={\displaystyle \sum _{m=1}^k\frac{\left(k-1\right)!}{\left(m-1\right)!\left(k-m\right)!}}\cdot \stackrel{\left(m-1\right)}{{}^{0}J_i\left[\overline{\theta }\left(t\right)\right]}\cdot \stackrel{\left[k-\left(m-1\right)\right]}{{}^0\overline{\mathcal{F}}_{X_i}^{\ast }}$$

(133)

$$\stackrel{\left(k\right)}{Q_g^i\left(t\right)}={\displaystyle \sum _{m=1}^k\frac{\left(k-1\right)!}{\left(m-1\right)!\left(k-m\right)!}}\cdot \stackrel{\left(m-1\right)}{{}^{0}J_i\left[\overline{\theta }\left(t\right)\right]}\cdot \stackrel{\left[k-\left(m-1\right)\right]}{{}^0\overline{\mathcal{F}}_{X_i}}$$

(134)

$$\stackrel{\left(k\right)}{Q_{SU}\left[\overline{\theta }\left(t\right)\right]}={\displaystyle \sum _{m=1}^k\frac{\left(k-1\right)!}{\left(m-1\right)!\left(k-m\right)!}}\cdot \stackrel{\left(m-1\right)}{{}_{0}J\left[\overline{\theta }\left(t\right)\right]}\cdot \stackrel{\left[k-\left(m-1\right)\right]}{{}^0\overline{\mathcal{F}}_X}$$

(135)

The meaning of the terms used in Equations (133) – (135) has been clearly defined in the previous sections, where $\left(k\ge 1\right)$ is the order of time derivation. However, in the context of dynamic equations, the notation$\left(k-1\right)$ is used instead of $\left(k\right)$.

When$\left(k=1\right)$, the equations (133) – (135) reduce to the classical forms: Equation (125), Equation (109), and finally Equation (114), respectively. By applying Lagrange, Gibbs-Appell and author’s own formulations [18,25], the generalized differential equations of higher order are obtained for mechanical systems that are dynamically characterized by sudden and transient motions

$$\begin{array}{c}\stackrel{\left(k-1\right)}{Q_i^j}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\dots ;\stackrel{\left(m\right)}{\overline{\theta }\left(t\right)}\right]\\ ={\displaystyle \frac{1}{m}}\cdot \left[{\displaystyle \frac{\partial \stackrel{\left(m+k-1\right)}{KE}}{\partial \stackrel{\left(m\right)}{q_j}}}-\left(m+1\right)\cdot {\displaystyle \frac{\partial \stackrel{\left(k-1\right)}{KE}}{\partial q_j}}\right]={\displaystyle \frac{\partial \stackrel{\left(m+k-3\right)}{E_A^{\left(1\right)}}}{\partial \stackrel{\left(m\right)}{q_j}}}={\displaystyle \frac{1}{m+1}}\cdot {\displaystyle \frac{\partial }{\partial \stackrel{\left(m\right)}{q_j}}}\left[2\cdot \stackrel{\left(m+k-5\right)}{E_A^{\left(2\right)}}+\stackrel{\left(m+k-3\right)}{E_A^{\left(1\right)}}\right]=\\ ={\displaystyle \frac{2}{\left(m+1\right)\cdot \left(m+2\right)}}\cdot {\displaystyle \frac{\partial }{\partial \stackrel{\left(m\right)}{q_j}}}\left[5\cdot \stackrel{\left(m+k-5\right)}{E_A^{\left(3\right)}}+2\cdot \stackrel{\left(m+k-3\right)}{E_A^{\left(2\right)}}+\right.\mathrm{...}\left.+\stackrel{\left(m+k-1\right)}{E_A^{\left(1\right)}}\right]\end{array}$$

(136)

$$\begin{array}{c}\mathrm{where} \stackrel{\left(0\right)}{E_A^{\left(1\right)}}=E_A^{\left(1\right)}\text{} j=1\to n, k=1 \\ m\ge \left[\left(k+1\right)=2\right] , \mathrm{and} \left(k\right) \mathrm{are} \mathrm{time} \mathrm{deriving} \mathrm{orders}\end{array}$$

(137)

$$\begin{gathered} \stackrel{\left(0\right)}{E_A^{\left(2\right)}}=E_A^{\left(2\right)}, j=1\to n, k=2, m\ge \left[\left(k+1\right)=3\right] \\ \mathrm{where} \stackrel{\left(0\right)}{E_A^{\left(3\right)}}=E_A^{\left(3\right)}, j=1\to n , k=3, m\ge \left[\left(k+1\right)=4\right], m=4,5,6,\dots \end{gathered}$$

(138)

The acceleration energies of first, second and third order have been defined both in the present study and in previous works [4,5,18,25].

Authors have proposed in [18], a formulation for the generalized higher – order differential equations, applicable to the mechanical systems (MBS), dynamically characterized by sudden and transient motions.

In these equations, the central function becomes the pth order acceleration energy, to which successive time derivatives of higher – order is applied, as follows:

$$\left\{{\displaystyle \frac{\left(k-1\right)!\cdot m!}{\left(m+k-1\right)!}}\cdot {\displaystyle \frac{\partial }{\partial \stackrel{\left(m\right)}{q_j}}}\left\{\left({\displaystyle \sum _{p=1}^k {\Delta }_p}\right)\cdot \stackrel{\left(m+k\right)-\left(2\cdot p+1\right)}{E_A^{\left(p\right)}}\right\}=\stackrel{\left(k-1\right)}{Q_{i\mathcal{F}}^j}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\dots ;\stackrel{\left(m\right)}{\overline{\theta }\left(t\right)}\right]\right\}$$

(139)

$$\mathrm{where} E_A^{\left(p\right)}=E_A^{\left(p\right)}\left[\overline{\theta }\left(t\right);\dot{\overline{\theta }}\left(t\right);\dots ;\stackrel{\left(p+1\right)}{\overline{\theta }\left(t\right)}\right] \mathrm{and} {\displaystyle \sum _{p=1}^k {\Delta }_p}={\displaystyle \sum _{p=1}^k\left[\frac{p\cdot \left(p+1\right)}{2}-{\delta }_p\right]}$$

(140)

The necessary conditions in (139) and (140) are the following:

$$\begin{array}{c}p=1\to k ; {\delta }_p=\left\{\left\{0 ; p=1\right\} ; \left\{1 ; p>1\right\}\right\}\\ \mathrm{and} k\ge 1; k=\left\{1;2; 3; 4; 5; \mathrm{...}\right\}; m\ge \left(k+1\right); m=\left\{2; 3; 4; 5; \mathrm{...}\right\}\end{array}$$

(141)

The expression (139) which is identical in form to Equation (133), represents the $k\ge 1$ order time derivative of the generalized inertia force. In the case of rapid motions, the Equation (122) is modified as follows:

$$\stackrel{\left(k-1\right)}{Q_m^i}\left(t\right)={\Delta }_m^2\cdot \left[{\Delta }_{\theta }\cdot Q_{i\mathcal{F}}^i\left(t\right)+Q_g^i\left(t\right)\right]+{\left(-1\right)}^{{\Delta }_m}\cdot {\displaystyle \frac{1-{\Delta }_m}{1+3\cdot {\Delta }_m}}\cdot \stackrel{\left(k-1\right)}{Q_{SU}^i}\left(t\right)$$

(142)

Thus, a system of n higher – order differential equations is obtained, describing the dynamic behavior of robotic mechanical structures subjected to rapid and transient motions.

The research presented in this paper introduces a novel and unified formulation of the higher – order dynamic equations modeling the dynamic behavior of multibody mechanical systems, particularly of those subjected to sudden accelerations and transient motions. Unlike classical approaches limited to first – order formulations, this study systematically integrates the acceleration energies of higher – order into a matrix – based model, enabling a more accurate and physically consistent modeling of robotic structures under dynamic stress. The introduction of generalized equations of p – order, together with the use of pseudo inertial tensors and advanced dynamic matrices, represents a significant original contribution to analytical dynamics and robotic modeling.

These formulations allow the direct connection between mass properties, higher – order kinematics and generalized forces, offering a clear perspective on the dynamic performance of robotic structures.

6. Conclusions

The paper is structured into four main sections, each addressing a fundamental aspect of advanced dynamics of mechanical systems: the distribution of mass properties, the formulation of higher-order acceleration energies, the generalized forces in multibody system dynamics, and the development of higher-order dynamic equations.

The first section introduces a novel and consistent framework for analytical modeling of multibody systems (MBS), particularly applicable to robotic structures subject to rapid dynamic transitions (transient and fast motions).

A key contribution is the formulation of the mass distribution (MD) concept, which generalizes the classical notion of mass geometry and emphasizes the dual role of mass and energy in defining matter and expressing its gravitational and inertial behavior.

This leads to a rigorous definition of the essential mass properties—total mass, center of masses, inertia tensor, its generalized variation law, and the pseudo-inertial tensor—initially for homogeneous bodies with simple geometry and later extended to complex subassemblies. These properties are essential in defining key dynamic quantities such as angular momentum, kinetic energy, and higher-order acceleration energies, which play a central role in the proposed higher-order differential equations.

In the second section, higher-order acceleration energies are introduced as a main component of the proposed dynamic model. The novelty lies in the formulation of these energies in matrix form, using homogeneous transformation matrices and matrix exponentials—an approach that is distinct from classical methods. Unlike traditional approaches, this study addresses holonomic systems undergoing general motion and integrates higher-order accelerations to improve the accuracy of dynamic models.

The third section focuses on the generalized forces such as gravitational, manipulation, inertial and friction forces, formulated based on d’Alembert’s principle and Jacobian based on transfer matrices of linear and angular velocities. These forces are explicitly expressed in a generalized form, for each driving joint and adapted to include friction effects at joint level. The final section extends the authors’ previous research by introducing generalized higher – order forces and dynamic equations obtained applying the analytical dynamics principles. These results are intended to model the behavior of complex mechanical structures subject to nonstationary regimes and transient inputs.

In conclusion, the mass distribution concept, the matrix – based formulation of higher – order acceleration energies, and the extension of Gibbs – Appell equations together define a unified and original theoretical framework for advanced dynamic modeling of multibody systems (MBS).

Compared to classical Newton – Euler or Lagrangian formulations, this model offers improved capability for capturing the influence of higher – order dynamics in fast and complex motion scenarios. Moreover, the higher – order differential equations developed here are capable of time – dependent behavior, showing strong potential for future applications in real – time control and high – speed robotics.

It is important to emphasize that this is a theoretical investigation focused exclusively on mathematical formulations and symbolic generalizations without experimental validation or numerical simulations. However, the presented formulations provide a solid foundation for future validation through computational simulations and experimental testing in robotic systems.