Inertia forces and shape integrals in the floating frame of reference formulation

Modeling and analysis of complex dynamical systems can be effectively performed using multibody system (MBS) simulation software. Many modern MBS packages are able to efficiently and reliably handle rigid and flexible bodies, often offering a wide choice of different formulations. Despite many advances in modeling of flexible systems, the most widely used formulation remains the well-established floating frame of reference formulation (FFRF). Although FFRF usually allows inclusion of only small elastic deformations, this assumption is adequate for many industrial applications. In addition, FFRF is computationally efficient if implemented with appropriate model order reduction techniques and effective handling of system inertia terms by utilization of so-called inertia shape integrals. However, derivation of the system of equations of motion for FFRF bodies is a complex and often error-prone task. The main goal of this paper is to provide a reliable, detailed, universal and clear set of inertia terms within the FFRF. The paper concentrates on detailed derivation of the inertia forces with focus on accurate determination and exploitation of the inertia shape integrals. Two standard methods are employed, namely the Lagrangian and Virtual Work approaches. Additionally, the introduced derivations are executed without selection of specific rotational parameters. Direct application of Euler parameters and Euler angles is presented. It is found that the derived expressions are well suited for direct implementation and simplify derivation of force components.


Introduction
Multibody system dynamics analysis is a numerical simulation tool used to construct and solve equations of motion for complex dynamic systems. In many practical engineering applications, a simplified multibody system dynamics approach can be used to analyze a system's motions and forces. Such an approach describes the complex dynamic system as a number of rigid interconnected bodies. However, in an increasing number of engineering applications, the current trend toward lightweight design is forcing researchers and designers to account for body deformation to ensure simulation accuracy and accurate input for strength-of-materials analysis. Deformation in multibody system dynamics can be defined with a number of different methods, of which the floating frame reference formulation (FFRF) is the approach most commonly used in industrial applications. In FFRF, body deformation is described with respect to the reference coordinate system, usually with conventional finite elements [6]. Consequently, FFRF can be applied to any arbitrary structural shape or type. Such an approach can, however, only be conveniently used for small elastic deformations, and large body reference motions must be described by translations and rotations of the reference coordinate system itself.
If the reference coordinate system is included in the deformation description and if a linear straindisplacement relationship is assumed, model order reduction techniques can be used to reduce the complexity of the calculations [9]. Model reduction enables computationally efficient description of structural flexibility such that real-time simulation becomes possible [1,4]. In the model reduction, the body flexibility is defined by assumed deformation shapes, which are commonly eigenvectors obtained from an eigenvalue analysis of the finite element model. As typical finite element models have many degrees of freedom, most practical applications utilize model reduction in FFRF.
In FFRF, body deformation and the motion of the reference coordinate system are coupled with inertial force descriptions. As a result, the inertial forces are nonlinear functions of body orientation and must be updated at each time step. The position, velocity and acceleration of every particle within the deformable body must thus be accounted for in the updating procedure. If body deformation with respect to the reference coordinate system is described using conventional finite elements, a flexible body is a discrete system comprising a finite number of nodal points. However, because the finite element models used in the flexibility description comprise a large number of nodal degrees of freedom, inertial force updating can become computationally expensive.
Identifying in advance inertial components that remain constant enables improvement in the computational efficiency of the updating procedure. Constant inertial components can be obtained from the deformation description, i.e., the finite element model, and they are referred to as the inertia shape integrals. The inertia shape integrals are added to formulate a discrete description of the deformation.
When discussing nonlinear inertial forces, many textbooks on multibody system dynamics are not providing a detailed description related to the inertia shape integrals. Moreover, many of these textbooks and many recent publications in the field describe rotational parameters as Euler angles or Euler parameters. However, a good understanding of the inertial shape integrals is essential for accurate and efficient simulation, and the inertia shape integrals must thus be carefully identified.
The main objective of this paper is to present a detailed derivation of the inertial forces for complex dynamic systems based on the floating frame of reference formulation. Appropriate formulas can be easily found in the literature, for example, in the papers of Yoo and Haug [10] and, more recently, Lugrís et al. [4] and Sherif and Nachbagauer [8] or, for a detailed derivation, in the well-known book by Shabana [7]. The present paper differs from its predecessors in the way it presents a detailed and clearly elaborated derivation of the velocity-dependent inertia forces, with a detailed and explicit use of the inertia shape integrals, and with the derivation accomplished using angular velocity and angular acceleration vectors. Using the presented transformation procedure for arbitrary rotational parameters, the provided expressions can be easily followed, verified and applied in computer code.
The use of angular velocity and acceleration results in a general expression of inertial forces because the well-known relationship between angular velocity vectors and the time derivation of rotational parameters can be used to express the inertial forces in any set of generalized rotational coordinates. The particular emphasis of the paper is the expression of inertia shape integrals in a form that simplifies implementation of FFRF. To verify the introduced expressions, inertial forces are derived using the concepts of Virtual Work and the Lagrangian approach.

System coordinates, basic kinematic equations and mass matrix
For simplicity and clarity, the body number index has been omitted from the variables in the following text.
Choice of the system coordinates In the current paper, q (t), where t is time, denotes a vector of generalized body coordinates that consists of the translational and rotational coordinates of the floating reference frame and the elastic coordinates. In addition, the following set of system velocities is used: R is the time derivative of the floating frame position with coordinates R (t), andω (t) is an angular velocity vector expressed in the local floating frame (a "bar" symbol indicates that the values are expressed in the local coordinate system of the body). The time derivativeṗ of a vector of elastic coordinates p (t) describes the deformation of the flexible body. The choice of the local angular velocity as a rotational velocity vector is common practice in the field of multibody dynamics. See Haug [3, see Eq. (11.3.11)] or Nikravesh [5, see formulation III, Eq. (11.24)]. However, there is no such concept as a rotation vector that can be differentiated to obtain the angular velocity. Instead, a differential rotation vector dπ (t) and a virtual rotation vector δπ (t) may be introduced to simplify the upcoming derivations. Using these vectors, the differential and virtual changes of generalized coordinates can be expressed as follows: Although the differential and virtual changes to the "true" coordinates (as position or elastic coordinates) are independent of the time derivatives, there is some relationship to the differential and virtual changes of the angular velocity vector for the differential rotational or virtual rotational vectors [2, see Sects. 4.12.4 and 7.3]: where A (θ ) is the orthogonal rotation matrix that can be expressed as a function of the rotational coordinates θ (t), and (˜) denotes the 3 × 3 skew-symmetric matrix associated with the vector used to express the vector (cross) product of two vectors in matrix notation. Similar expressions can be written in a global coordinate system with the global angular velocity vector ω (t): The choice of specific system coordinates q, in particular rotational coordinates, will be discussed in Sect. 6. Figure 1 presents the deformable body modeled using the floating frame of reference approach. In this approach, the position vector r P (q) of an arbitrary point P of the flexible body can be expressed as:

Basic kinematic equations
u ( p) is the position vector that points to P in the local coordinate frame,ū 0 is the position vector that points to P in the reference configuration (usually undeformed), u f ( p) is the deformation vector that accounts for (very often small) displacements of the point P with respect to its reference position, and S is a space-dependent matrix that relates the body deformations and elastic coordinates. For nodal coordinates and the finite element method, matrix S is a shape function matrix.
Using Eq. (8), the velocity vector of an arbitrary point P can be written: where the equalityȦ = Aω is used, and the property of the "tilde" operator is exploited, which for two arbitrary vectors a and b givesãb = −ba. The time independence ofū 0 and S is used to derive the third equation component. Finally, the equation can be rewritten in a form that employs the total vector of generalized velocitiesq and matrix L (θ, p): I denotes an identity matrix of an appropriate size (here 3 × 3), and B = −Aũ. It should be noted that quantities that do not explicitly depend on time, system coordinates or spatial variables, like the constant identity matrix I, are designated in the text with upright Roman style.
Flexible body mass matrix Using the L matrix, a mass matrix can be written as follows: where m, ρ and V are the mass, the density and the volume of the deformable body in its reference (initial) configuration. Note that dm = ρdV . The mass matrix can be expressed in the following compact form: m αβ for α, β = t, r, f denotes the blocks of the mass matrix split according to the system coordinates for t associated with translational coordinates, r with rotational coordinates and f with elastic (flexible) coordinates.

Mass components and inertia shape integrals
Aspects of the mass matrix components and inertia shape integrals that simplify the matrix updating process are discussed in the following paragraphs. The shape integrals of the mass matrix will be denoted as capital I with a superscript marker and with no additional designation. Matrices denoted as I with additional subscripts or designated with bars and hats are matrices that depend on the shape integrals. Additionally, they may depend also on the body coordinates.

Translational component
The simplest form has a purely translational component: where I 1 is the mass of the deformable body: Translational-rotational component The translationalrotational component m tr (θ , p) is expressed as follows: while where Eq. (8) is partially used, and the definitions for the shape integrals I 2 and I 3 are: and Rotational component The purely rotational component m rr ( p) is written as: where I rr = V ρũ Tũ dV is a symmetric matrix, and: whereū i =ū 0,i + S i p for i = 1, 2, 3 is an ith component of the vectorū,ū 0,i is an ith component of the vectorū 0 , and S i is an ith row of the matrix S.
The following expression is then evident: For j = 1, 2, 3: Therefore, where the six shape integrals I 7 i j are: the nine shape integrals I 8 i j are: and the six shape integrals I 9 i j (the following three are symmetric) are: Next, the matrix I rr can be written using the shape integrals from Eqs. (26), (27) and (28): which is a symmetric matrix. Moreover, where I 8 is a constant matrix of shape integrals: In addition, using the property I 9 ji = I 9 i j T that directly results from definition (28) results in: T p p T I 9 11 + I 9 33 p − p T I 9 23 p − p T I 9 where I 9 is a constant and symmetric matrix: where the orthogonality of the rotation matrix and the skew-symmetric matrix propertyũ T = −ũ is used. Next, the following expression can be written as: Assuming, as previously, that i = 1, 2, 3 and j = 1, 2, 3: thus, Flexible component The purely flexible sub-matrix of the mass matrix, which is constant, is expressed as: where

Quadratic velocity vector: Lagrangian approach
The quadratic velocity vector can be evaluated based on differentiation of kinetic energy T = T (q,q) [7]: where is a vector of quadratic velocity forces. The kinetic energy of a flexible body may be expressed as: where M = M (q) is a symmetric mass matrix of a system. Therefore, an expression for the inertial forces can be written as: while the quadratic velocity vector may be expressed as: It should be noted that mass sub-matrices m tt and m f f are constant, while m tr (θ , p), m t f (θ ), m rr ( p) and m r f ( p) depend on system coordinates. The quadratic velocity vector based on Eq. (48) can be written as: while: Because the mass matrix is symmetric, the equation can be simplified further: In addition, the vector of the quadratic velocity may be divided into translational, rotational and elastic components: In the following paragraphs, each of these components will be considered separately.

Translational component
Considering Eqs. (48), (49) and (51), Q v,t may be expressed as follows: None of the terms from Eq. (51) are included, because the mass matrix does not depend directly on R. It can be shown that: Remembering thatȦ =ω A = Aω, and becauseİ j = I 3ṗ from Eq. (16): The second term in Eq. (53) is: Therefore, the value of the quadratic velocity vector for that component is: Q v,r = − ṁ rtṘ +ṁ rrω +ṁ r fṗ Even though m rr ( p) and m r f ( p) do not depend on rotational coordinates, the termsω T m rrω andω T m r fṗ must be included, because the angular velocity vector is dependent on rotational coordinates.
The following equation can also be shown to be valid: Because matrix m rt (θ , p) is a function of rotational and elastic coordinates, the expression can also be written as: Firstly, Starting from the left-hand side, using Eq. (15) and where the equalities ∂ ∂π A T v = A Tṽ A, ω = Aω andω = A Tω A are used and orthogonality of the rotation matrix and some properties of skew-symmetric matrices. Note that: for arbitrary vectors a and b that are a function of the real variable vector q. Transposing both sides of Eq. (63) yields: Equation (64) can be applied to evaluate the right-hand side of Eq. (61) with a T =Ṙ T m tr (so a = m rtṘ ) and b =ω. According to Eq. (4), ∂ω/∂π =ω. Therefore, where properties (ãb) =ãb−bã and A Tã A = A T a (for arbitrary vectors a and b) are used. Notice that both sides of Eq. (61), that is Eqs. (62) and (65), are equal. Next, The left-hand side of Eq. (66) is equal to: where the symbolR = A TṘ is introduced, and the equalities ∂ I j /∂ p = I 3 and I 3ṗ =İ j are assumed. Next, the right-hand side is equal to: where the equality ∂ ( Av) /∂π = −Aṽ (for arbitrary vector v) is used. Both sides of the equation are equal. Eventually, because Eq. (59) is valid, Eq. (58) simplifies into: Starting with the termṁ rrω , and using the definition of the purely rotational sub-matrix of the mass matrix given in Eq. (20) leads to: The time derivative of the I rr can be expressed as: Therefore, Thus, matrixİ rr is described with the shape integrals I 8 i j and I 9 i j . The second term in Eq. (69) is: where I r f is given by Eq.
because for an arbitrary vector a and arbitrary square matrix A of proper sizes it can be shown that The next term is: where once again the symmetry of the matrix I rr is assumed.
The last term of the Q v,r is: Finally, It can be shown that: Starting from the left-hand side: which results from the equality A TȦ =ω. On the right-hand side, Therefore, both sides cancel to form a simplified quadratic velocity vector: The first of the remaining terms is: The first term in parenthesis is: where I 8 is defined in Eq. (32),Î 9 is given by Eq. (72) and: whereω T = ω 1ω2ω3 , I n×n is an n × n identity matrix, while n is the number of elastic coordinates. The last term of the quadratic velocity vector is: Finally, an expression for the quadratic velocity vector for the elastic component is: The reader should be able to easily calculate the value of the quadratic velocity vector based on the inertia shape integrals provided in Sect. 4, which will result in an efficient force evaluation without having to recalculate any mass integrals.

Quadratic velocity vector: virtual work approach
Next, the quadratic velocity vector will be developed using the virtual work of the inertial forces: Based on Eq. (9), the acceleration vector can be expressed as: Then, the virtual change of the position vector δr P is: where virtual change of the system parameters can be written as: Therefore, the virtual work can be expressed as: where Q i (q,q) is a vector of inertial forces given by: The expression for the vector of inertial forces results in the definition used previously in Eq. (47). Q v (q,q) is a vector of quadratic velocity forces, which can be further defined as: The integrant of the quadratic velocity force, using Eq. (10), is: based on equalitiesȦ = Aω andu =u f = Sṗ: Finally, the vector of the quadratic velocity forces is: where partitioning of the force vector into translational, rotational and elastic components, as in Eq. (52), is used.
In the following section, each component will be treated separately.

Translational component
The translational component can be written: By substituting Eqs. (16) and (18): which is the same equation as in Eq. (57).
(102) can be written as: Equation (115) has the same form as Eq. (77), which is obtained by the Lagrangian approach.

Flexible (elastic) component
The elastic component of Q v can be determined as: In addition, Therefore, A direct comparison of the two above equations shows that: The next term in the elastic component of the Q v vector is: where the equality thatu f =u = Sṗ is used. This can be expressed as: then, Finally, the vector of the elastic component of quadratic velocity forces derived from the principle of virtual work of the inertial forces is: which is equivalent to Eq. (87) obtained using a direct differentiation of the kinetic energy.

Quadratic velocity force vector
The translational, rotational and elastic components of the quadratic velocity force vector calculated based on the virtual work approach leads to expressions equivalent to those derived using the Lagrangian approach. Therefore, the final form of this force component is exactly the same as in Eq. (88).

Using orientation parameters
Utilization of the concepts of angular velocity and differential or virtual rotation vectors is convenient when developing a general form of the quadratic velocity vector such as in Eq. (88). However, specific coordinate types are used in many practical applications and implementations. Popular choices include Euler angles, Euler parameters, Rodriguez parameters, Cayley's parameters, and others [2,7]. To offer the option of selecting a specific coordinate type, the quadratic velocity vector must be expressed as a function of arbitrary rotation parameters θ (t). The following paragraphs describe how angular velocity and virtual rotations can be transformed to produce a vector of rotational parameters and its time derivatives. The transformation is carried out by applying the technique of virtual work, which was introduced in Sect. 5. According to Shabana [7], an angular velocity vector can be calculated as a function of the rotation parameters and derivative of the rotation parameters as: whereḠ is velocity transformation matrix that depends only on rotational coordinates. Analogously, The coordinates can be transformed from q toq: Using Eq. (93), the following equation can be written as: Therefore, the transformed inertial and quadratic velocity forces are: Because of the structure of the matrix T (θ), only the rotational component of the forces must be transformed. So, multiplying the left side of Eq. (77) bȳ G T yields the rotational component of the quadratic velocity forces after coordinate transformation: Next, the inertial forces must be transformed. The equations that correspond with rotations should be multiplied byḠ T , as in the case of the quadratic velocity vector. However, one should note that the inertial forces, given for example by (47), explicitly depend on the acceleration vectorq. Therefore, those accelerations should be expressed as a function of the new rotational variables. Taking the time derivative of Eq. (128) gives: Finally, the vectorQ i can be expressed as: The second part of the above expression is, in fact, a quadratic velocity term, which appears due to the differentiation of angular velocity, and it should be included in the definition of the quadratic velocity vector after the coordinate transformation. One can easily verify that: Therefore, the final expression of the quadratic velocity vector, which should be used for most rotational parameters, is: Moreover, the value of the inertial forces can be written as follows:

Quadratic velocity vector with Euler parameters
Assuming that the vector θ is a set of the Euler parameters, the well-known equalitiesĠθ = 0,ω = 1 2ḠĠ T andḠ TḠ = 4 I − θθ T can be used. The termḠ Tω can be simplified as: With the equality θ T δθ = 0, the virtual work of the second part of the above equation can be shown to be equal to zero, i.e., 2Ġ T θθ T δθ = 0. Therefore, Eq.
(137) can be written in the simplified form in terms of Euler parameters as:

Quadratic velocity vector with Euler angles
The literature does not offer many general derivations of Euler angles. Assuming that Euler angles describe consecutive rotations along the ZXZ axis, as is commonly done in the field of robotics, the current section vector of rotational parameters is: where φ, θ and ψ are proper angles. In local coordinates, the angular velocity vector can be described with Eq. (128). For the Euler angles, the velocity transformation matrix is [7]: IfĠθ = 0 andḠ Tω = 2Ġ T or equivalently −ωḠ = 2Ġ , then the simplifications made for the Euler parameters should also be applicable to the Euler angles. The direct calculation shows that: In addition, a direct differentiation with respect to time leads to: cos ψ sin θ +θ sin ψ cos θ −ψ sin ψ 0 −ψ sin ψ sin θ +θ cos ψ cos θ −ψ cos ψ 0 −θ sin θ 0 0 ⎤ ⎦ .
Looking at the first term and Using the notations sψ and cψ to denote, respectively, sinus and co-sinus of the angle ψ:Ḡθ Therefore,Ġ is a part ofωḠ, but the first term on the right-hand side of the above equation is not equal to −Ġ, so it cannot be generally stated thatḠ Tω is equal to 2Ġ T . When the Euler angles are used as a rotational parameters, Eq. (137) must be used directly without any further assumptions to obtain the correct value of the quadratic velocity vector. It is worth pointing out that for most choices of the rotational parameters Eq. (137) should be used for the velocity-dependent inertia forces and Eq. (141) is proven only for the Euler parameters.

Conclusions
FFRF is a well-known method in dynamic modeling of flexible bodies within a multibody framework. Yet, essential elements of its implementation, such as the representation of inertia or Coriolis, centrifugal and gyroscopic forces that result when kinetic energy is differentiated, are usually not described in sufficient detail. As a result, complex and error-prone derivations are often required, for example as in [7] and [8]. Moreover, most of the derivations presented in the literature are specific to a unique set of rotational coordinates and might not meet the requirements of a particular analysis. Equations derived for the Euler parameters can be especially misleading, because the Euler parameter identities often result in modified and simplified versions of force expressions that cannot be used for other rotational parameter choices. The current paper presents detailed derivations of the inertia and quadratic velocity force vectors within the FFRF. All force terms are presented as a function of inertia shape integrals to avoid costly integral updates at each time integration point. In addition, no specific assumptions are made for the choice of rotational parameters. Instead, the concepts of angular velocity and virtual or differential rotational vectors are applied. As a consequence, the presented expressions are well suited for direct implementation for any set of coordinates. In addition, the use of the angular velocity concept greatly simplifies the derivation of the force components.
In the description of inertial forces, the main challenge is often associated with the derivation of the quadratic velocity terms, especially the rotational component. Therefore, in this paper those forces are derived using two approaches-the Lagrangian approach and the virtual work principle. The final expressions perfectly agree for both attempts.
In an FFRF implementation, the choice of rotational parameters is often challenging. While this paper does not discuss how to make the choice, general transformation to an arbitrary set of rotational parameters is discussed in detail. In addition, two common choices are considered: using Euler parameters or using Euler angles. The simplifications made when using Euler parameters are described. When using Euler angles, the simplifications are not general, and in the most cases, generic force expressions must be applied.