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

Nonlinear Dynamics, Jan 2017

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs11071-017-3355-y.pdf

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

Inertia forces and shape integrals in the floating frame of reference formulation Grzegorz Orzechowski Marko K. Matikainen Aki M. Mikkola 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. Flexible multibody dynamics; Floating frame of reference formulation; Inertia forces; Quadratic velocity vector; Velocity-dependent terms; Inertia shape integrals 1 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 strain– displacement 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. 2 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. Fig. 1 Coordinates used for kinematic description in the floating frame of reference formulation. X Y Z denotes the global reference frame, while x yz is a local, floating frame of reference. P is an arbitrary point in the current configuration that corresponds to the point P0 in the undeformed configuration 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˙ ⎤ q˙ = ⎣ ω¯ ⎦ . p ˙ 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 p˙ 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: dq = 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. Basic kinematic equations 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: rP = R+A u¯ = R + A u¯0 + u¯f u¯ ( p) is the position vector that points to P in the local coordinate frame, u¯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˙ = Aω˜¯ is used, and the property of the “tilde” operator is exploited, which for two arbitrary vectors a and b gives a˜b = − b˜a. The time independence of u¯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 velocities q˙ and matrix L (θ , p): L = I − A u¯˜ AS = I B AS . I denotes an identity matrix of an appropriate size (here 3 × 3), and B = − A u˜¯. 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: M = LT Ldm = 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 = ⎣ ⎡ mtt mtr mt f ⎤ mrr mr f ⎦ . sym m f f 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. 3 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: mtt = where I1 is the mass of the deformable body: I1 = Translational-rotational component The translationalrotational component mtr (θ , p) is expressed as follows: mtr = I j = I2 = ρ ( u¯0 + S p) dV = I2 + I3 p, where Eq. (8) is partially used, and the definitions for the shape integrals I2 and I3 are: I3 = Translational-flexible component The translationalflexible component mt f (θ ) is expressed as: ρu¯i2dV = Iii + 2Ii8i p + pT Ii9i p. 7 For j = 1, 2, 3: +u¯0,i S j p + u¯0, j Si p + pT SiT S j p. ⎡ 0 u¯3 −u¯2 ⎤ ⎡ u˜¯T u˜¯ = ⎣ −u¯3 0 u¯1 ⎦ ⎣ u¯2 −u¯1 0 ⎡ u¯22 + u¯32 −u¯1u¯2 −u¯1u¯3 ⎤ = ⎣ −−uu¯¯23uu¯¯11 u¯−12u¯+3u¯u¯232 u¯−21u¯+2u¯u¯322 ⎦ , where u¯i = u¯0,i + Si p for i = 1, 2, 3 is an i th component of the vector u¯, u¯0,i is an i th component of the vector u¯0, and Si is an i th row of the matrix S. The following expression is then evident: u¯i2 = u¯0,i + Si p 2 = u¯02,i + 2u¯0,i Si p + pT SiT Si p. (22) mt f = Rotational component The purely rotational component mrr ( p) is written as: mrr = Ii7j = Ii8j = the nine shape integrals Ii8j are: where the six shape integrals Ii7j are: and the six shape integrals Ii9j (the following three are symmetric) are: Ii9j = Next, the matrix Irr can be written using the shape integrals from Eqs. (26), (27) and (28): where p¯ = diag ( p, p, p), and which is a symmetric matrix. Moreover, I 8 ⎡ 2 I282 + I383 − I812 + I821 − I813 + I8381 ⎤ 8 8 8 8 8 ¯ = ⎣ − I12 + I21 2 I11 + I33 − I23 + I832 ⎦ 8 8 8 8 8 − I13 + I31 − I23 + I32 2 I11 + I22 ⎡ p 0 0 ⎤ ⎣ 0 p 0 ⎦ = I8 p¯, 0 0 p where I8 is a constant matrix of shape integrals: ⎡ I922 + I393 −I912 −I193 ⎤ 9 T I191 + I933 9 −− II11923 T − I293 T I911−+I23I292 ⎦⎥ p¯ = p¯T I9 p¯, where I9 is a constant and symmetric matrix: I9 = ⎢ ⎣ 9 T I11 + I33 9 T I11 + I22 Rotational-flexible component The rotational-flexible sub-matrix mr f ( p) of the mass matrix is: mr f = where the orthogonality of the rotation matrix and the skew-symmetric matrix property u˜¯T = − u¯˜ is used. Next, the following expression can be written as: ⎡ −u¯3 S2 + u¯2 S3 ⎤ u˜¯S = ⎣ u¯3 S1 − u¯1 S3 ⎦ . −u¯2 S1 + u¯1 S2 Assuming, as previously, that i = 1, 2, 3 and j = 1, 2, 3: Ir f = In addition, using the property I9ji = Ii9j T that directly results from definition (28) results in: ⎡ pT I293 − I293 T ⎤ I¯ 5 = ⎢⎢⎢⎣ ppTT II919123T−−I91I291T3 ⎥⎥⎥⎦ = ⎣ ⎡ pT 0 0 ⎤ ⎡ I923 − I923 T ⎤ 0 pT 0T ⎦ ⎣⎢ I913T −9I19T3 ⎥⎦ = p¯T I5. 0 0 p I192 − I12 Flexible component The purely flexible sub-matrix of the mass matrix, which is constant, is expressed as: m f f = I6 = 4 Quadratic velocity vector: Lagrangian approach The quadratic velocity vector can be evaluated based on differentiation of kinetic energy T = T (q, q˙ ) [7]: Qi − Qv = dt where Qi = Qi (q, q¨ ) is a vector of inertial forces, and Qv = Qv (q, q˙ ) is a vector of quadratic velocity forces. The kinetic energy of a flexible body may be expressed as: Qi = M q¨ , 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: 1 ∂ Qv = − M˙ q˙ + 2 ∂ q It should be noted that mass sub-matrices mtt and m f f are constant, while mtr (θ , p), mt f (θ ), mrr ( p) and mr f ( p) depend on system coordinates. The quadratic velocity vector based on Eq. (48) can be written as: m˙ tt R˙ + m˙tr ω¯ + m˙t f p˙ ⎤ m˙rt R˙ + m˙rr ω¯ + m˙r f p˙ ⎦ m˙ f t R˙ + m˙ f r ω¯ + m˙ f f p˙ ⎡ m˙tr ω¯ + m˙t f p˙ ⎤ = ⎣ m˙rt R˙ + m˙rr ω¯ + m˙r f p˙ ⎦ , m˙ f t R˙ + m˙ f r ω¯ qT M q˙ = R˙T mtt R˙ + R˙T mtr ω¯ + R˙T mt f p˙ ˙ Because the mass matrix is symmetric, the equation can be simplified further: qT M q˙ = R˙T mtt R˙ + 2 R˙T mtr ω¯ + 2 R˙T mt f p˙ ˙ In addition, the vector of the quadratic velocity may be divided into translational, rotational and elastic components: Qv = ⎣ In the following paragraphs, each of these components will be considered separately. 4.1 Translational component Considering Eqs. (48), (49) and (51), Qv,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: − A I˜ j ω¯ = − A˙ I˜ j ω¯ − A I˜˙ j ω¯ . Remembering that A˙ = ω˜ A = Aω˜¯ , and because I˙ j = I3 p˙ from Eq. (16): The second term in Eq. (53) is: Therefore, the value of the quadratic velocity vector for that component is: m˙rt R˙ = 4.2 Rotational component As previously, with Eqs. (48), (49) and (51): Even though mrr ( p) and mr f ( p) do not depend on rotational coordinates, the terms ω¯ T mrr ω¯ and ω¯ T mr f p˙ must be included, because the angular velocity vector is dependent on rotational coordinates. The following equation can also be shown to be valid: − m˙rt R˙ + = 0. (59) Because matrix mrt (θ , p) is a function of rotational and elastic coordinates, the expression can also be written as: Starting from the left-hand side, using Eq. (15) and remembering that mrt = mtTr = I˜ j AT : I˜ j AT R˙ = I˜ j = I˜ j AT R˜˙ω = − I˜ j AT ω˜ A AT R˙ where the equalities ∂∂π¯ AT v = AT v˜ A, ω = Aω¯ and ω˜¯ = AT ω˜ 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 aT = R˙T mtr (so a = mrt R˙) and b = ω¯ . According to Eq. (4), ∂ω¯ /∂π¯ = ω¯˜ . Therefore, T I˜ j AT R˙ = ω¯˜ T I˜ j AT R˙ + I˜ j AT R˜˙ A = −ω¯˜ I˜ j AT R˙ − I˜ j ω¯ AT R˙ where properties ( a˜b) = a˜ b˜− b˜ a˜ and AT a˜ A = AT 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, p˙ = The left-hand side of Eq. (66) is equal to: I˜ j AT R˙ = − = − R˙¯˜I3 p˙ = I˜˙ j AT R˙, p ˙ p ˙ p˙ = where the symbol R¯˙ = AT R˙ is introduced, and the equalities ∂ I j /∂ p = I3 and I3 p˙ = I˙ j are assumed. Next, the right-hand side is equal to: = R˙ = R˙T T ∂ A I˙ j − A I˜˙ j = I˜˙ j AT R˙, where the equality ∂ ( Av) /∂π¯ = − Av˜ (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 m˙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 Irr can be expressed as: I˙ rr = I8 ˙¯p + p¯T Iˆ9 ˙¯p. Thus, matrix I˙ rr is described with the shape integrals Ii8j and Ii9j . The second term in Eq. (69) is: m˙r f p˙ = I˙ r f p˙, where Ir f is given by Eq. (40), and I˙ r f p˙ = ¯˙pT I5 p˙ = I˙¯ 5 p˙ = 0 as: ⎡ p˙T I923 − I293 T p˙ ⎤ where once again the symmetry of the matrix Irr is assumed. The last term of the Qv,r is: It can be shown that: m˙ f t R˙ = ∂∂p R˙T mtr ω¯ T . Starting from the left-hand side: d m˙ f t R˙ = dt AI3 T R˙ = I3 T A˙T R˙ = I3 T ω˜¯ T AT R˙, The first of the remaining terms is: d m˙ f r ω¯ = dt The first term in parenthesis is: T IrTf ω¯ = I˙ rTf ω¯ = I5 T ˙¯pω¯ . = 21 ω¯ T I8ωˆ¯ + p¯T Iˆ9ωˆ¯ which results from the equality AT A˙ = ω¯˜ . On the right-hand side, = − 4.3 Flexible (elastic) component Using Eqs. (48), (49) and (51), the flexible (elastic) component of the quadratic velocity vector can be written: where I8 is defined in Eq. (32), Iˆ9 is given by Eq. (72) and: = − ω¯ T ˙¯pT I5 T I4 p˙ + p¯T I5 p˙ = − I5 T ˙¯pω¯ . (86) Finally, an expression for the quadratic velocity vector for the elastic component is: Qv, f = −2 I5 T ˙¯pω¯ + 21 ωˆ¯ T 4.4 Quadratic velocity force vector Qv = ⎣ 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. 5 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: where virtual change of the system parameters can be written as: r¨ P = L˙q˙ + Lq¨ . = − Therefore, the virtual work can be expressed as: where Qi (q, q¨ ) is a vector of inertial forces given by: Qi = The expression for the vector of inertial forces results in the definition used previously in Eq. (47). Qv (q, q˙ ) is a vector of quadratic velocity forces, which can be further defined as: Qv = − The integrant of the quadratic velocity force, using Eq. (10), is: = −LT = −LT − A˙ u˜¯ + A u˜˙¯ ω¯ + A˙S p˙ , Qiv = −LT = −LT Finally, the vector of the quadratic velocity forces is: Therefore, the first term in Eq. (102) becomes: = − treated separately. Aω˜ ω˜ u 2Aω˜ S p ¯ ¯ ¯ + ¯ ˙ T B Aω˜ ω˜ u 2Aω˜ S p ¯ ¯ ¯ + ¯ ˙ ⎥ dV, S A Aω˜ ω˜ u 2Aω˜ S p ¯ ¯ ¯ + ¯ ˙ where partitioning of the force vector into translational, rotational and elastic components, as in Eq. (52), is In the following section, each component will be 5.1 Translational component The translational component can be written: v,t = ρ Aω˜ ω˜ u 2Aω˜ S p dV − ¯ ¯ ¯ − ¯ ˙ Aω˜ ω˜ ρudV 2Aω˜ ρSdV p. = − ¯ ¯ ¯ − ¯ ˙ By substituting Eqs. (16) and (18): which is the same equation as in Eq. (57). 5.2 Rotational component The rotational component for the vector of quadratic velocity forces can be written as: v,r = − T ρ B Aω˜ ω˜ u 2Aω˜ S p dV ¯ ¯ ¯ + ¯ ˙ v,r = − ρ u˜ ω˜ ω˜ u 2ω˜ S p dV. ¯ ¯ ¯ ¯ + ¯ ˙ Using the equalities aa 0 and ab˜ b˜a (ab) or ˜ = ˜ − ˜ = ˜ ab˜ b˜a (ab), this expression can be written as: ˜ = ˜ + ˜ u˜ω˜ ω˜ u ω˜ u˜ u˜ω u˜ω ¯ ¯ ¯ ¯ = − ¯ ¯ + ¯ ¯ ¯ ¯ T T ω˜ u˜ u˜ω u˜ω u˜ω ω˜ u˜ u˜ω. (103) = ¯ ¯ ¯ ¯ − ¯ ¯ ¯ ¯ = ¯ ¯ ¯ ¯ ρ u˜ω˜ ω˜ udV ¯ ¯ ¯ ¯ = − T ω˜ ρ u˜ u˜dV ω = − ¯ ¯ ¯ ¯ and, using Eq. (20): T − ¯ ¯ ¯ ¯ = −ω˜¯ Irrω¯ . ω˜ ρ u˜ u˜dV ω Then, the second term can be written as: 2 ρ u˜ω˜ S pdV ¯ ¯ ˙ = − ⎡ u¯3u˙¯3 + u¯2u¯˙2 Assuming i 1, 2, 3 and j 1, 2, 3: = = Thus, by using Eqs. (27) and (28): 8 T 9 ρu¯iu˙¯ jdV I p p I p = ij ˙ + ij ˙ = ⎢ ⎡ pT I9 I9 p pT I9 T p 22 + 33 ˙ − 12 ˙ T 9 T 9 9 T 9 T p I p p I I p p I p ⎥ − 12 ˙ 11 + 33 ˙ − 23 ˙ ⎦ T 9 T 9 9 p I p p I I p − 23 ˙ 11 + 22 ˙ 8 8 8 8 I I I I p = ⎢ − 12 11 + 33 − 32 ⎥ ¯˙ I I I − 23 11 + 22 9 9 9 9 T I I I I ⎥ ¯˙p. − 12 11 + 33 − 23 ⎦ I I I − 23 11 + 22 Next, direct calculations reveal that: T ρ u˜ u˜˙dV I ¯ ¯ = ˙rr − be expanded: Using Eqs. (40), (41) and (42), the expression I p can r f ˙ ⎢ T 9 T 9 ⎥ p I I p. 13 − 13 ⎥ ˙ 8 8 8 8 ⎤ I I p I I p 21 − 12 ˙ 31 − 13 ˙ I I ⎣ 8 8 8 8 I I p I I p 0 13 − 31 ˙ 23 − 32 ˙ ⎢ ⎣ T 9 9 T T 9 9 T p I I p p I I p 13 − 13 ˙ 23 − 23 ˙ T 9 T 9 T 9 T 9 ⎤ p I I p p I I p 12 − 12 ˙ 13 − 13 ˙ Now inserting Eqs. (71) and (113) into Eq. (111) using Eqs. (33) and (72): ⎡ 8 8 8 8 8 8 ⎤ 2 I I I I I I 22 + 33 − 12 + 21 − 13 + 31 1 ⎢ 8 8 8 8 8 8 ⎥ I I 2 I I I I = ⎢ − 12 + 21 11 + 33 − 23 + 32 ⎥ ¯˙p 2 ⎢ I I I I 2 I I 1 8 8 8 8 − ⎢ I12 −I21 0 I32 −I23 ⎥ ¯˙p 2 ⎣ ⎦ 0 I8 I8 I8 I8 ⎤ 21 − 12 31 − 13 I I I I 0 2 I I I I I I I I I I 2 I I 9 T 9 9 T 9 ⎤ 0 I I I I 12 − 12 13 − 13 I I I I ⎡ I8 I8 I8 I8 ⎤ 22 + 33 − 21 − 31 8 8 8 8 = ⎢ −I12 I11 +I33 −I32 ⎥ ¯˙p ⎣ 8 8 8 8 ⎦ I I I I − 13 − 23 11 + 22 ⎡ I9 I9 I9 T I9 T 22 + 33 12 − 13 T 9 9 9 9 T +p¯ ⎢ I I I I ⎥ ¯˙p, ⎣ − 12 11 + 33 − 23 ⎦ 9 9 9 9 I I I I − 13 − 23 11 + 22 which shows that Eq. (111) is valid. Therefore, Eq. (102) can be written as: Q I p ω v,r = −ω¯˜ Irrω¯ − I˙rrω¯ + rf ˙ ¯ Equation (115) has the same form as Eq. (77), which is obtained by the Lagrangian approach. 5.3 Flexible (elastic) component The elastic component of Q can be determined as: v In addition, T T T T S ω˜ω˜u S S S ¯ ¯ ¯ = 1 2 3 T T ρS A Aω˜ω˜ u 2Aω˜ Sp dV (116) ¯ ¯ ¯ + ¯ ˙ T T ρ S ω˜ω˜ u 2S ω˜ Sp dV. (117) ¯ ¯ ¯ + ¯ ˙ 2 2 ⎡ −ω¯2 − ω¯3 ω¯1ω¯2 ω¯1ω¯3 ⎤ ⎡ u¯1 ⎤ 2 2 ω¯1ω¯2 −ω¯1 − ω¯3 ω¯2ω¯3 2 2 ω¯1ω¯3 ω¯2ω¯3 −ω¯1 − ω¯2 u¯3 2 2 ⎡ ω¯1ω¯2u¯2 + ω¯1ω¯3u¯3 − u¯1(ω¯2 + ¯3 ω )⎤ 2 2 ω¯1ω¯2u¯1 + ω¯2ω¯3u¯3 − u¯2(ω¯1 + ¯3 ⎦ ω ) 2 2 ω¯1ω¯3u¯1 + ω¯2ω¯3u¯2 − u¯3(ω¯1 + ¯2 ω ) T T S u S u = ω¯1ω¯2 1 ¯2 + 2 ¯1 T T S u S u + ω¯1ω¯3 1 ¯3 + 3 ¯1 T T S u S u + ω¯2ω¯3 2 ¯3 + 2 ¯1 T 2 2 T 2 2 S u ω S u ω − 1 ¯1 ω¯2 + ¯3 − 2 ¯2 ω¯1 + ¯3 T 2 2 S u ω − 3 ¯3 ω¯1 + ¯2 − I911 T p ω¯ 22 + ω¯ 32 − I292 T p ω¯ 12 + ω¯ 32 9 T p ω¯ 1 + ω¯ 22 . − I33 2 On the other hand, 2 I922 + I393 ⎢⎢⎢ − I912 T + I192 ⎣ − I13 + I193 − I923 T + I293 9 T 2 I911 + I33 +2 I181 T ω¯ 2 + ω¯ 3 2 2 − 21 2 I812 + I281 T ω¯ 1ω¯ 2 + 2 I813 + I381 T ω¯ 1ω¯ 3 + 21 2 I922 T p ω¯ 12 + ω¯ 32 + 2 I933 T p ω¯ 12 + ω¯ 2 2 +2 I191 T p ω¯ 22 + ω¯ 3 2 2 I911 + I292 p ¯ A direct comparison of the two above equations shows that: The next term in the elastic component of the Qv vector is: − 2 = −2 where the equality that u˙¯ f = u¯˙ = S p˙ is used. This can be expressed as: ⎡ 0 u˙¯3 −u¯˙2 ⎤ ST u˜˙¯T = S1T S2T S3T ⎣ −u¯˙3 0 u˙¯1 ⎦ u¯˙2 −u¯˙1 0 S3T S2 − S2T S3 p˙ S1T S3 − S3T S1 p˙ S2T S1 − S1T S2 p˙ . (124) − 2 ρ ST ω˜¯ S p˙dV = −2 I˙ r f T ω¯ = −2 I5 T ˙¯pω¯ . Finally, the vector of the elastic component of quadratic velocity forces derived from the principle of virtual work of the inertial forces is: p¯ ω¯ − 2 I5 T ˙¯pω¯ , (127) which is equivalent to Eq. (87) obtained using a direct differentiation of the kinetic energy. 5.4 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). 6 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 G¯ is velocity transformation matrix that depends only on rotational coordinates. Analogously, The coordinates can be transformed from q to qˆ : Therefore, the transformed inertial and quadratic velocity forces are: Qˇv,r = − G¯ Next, the inertial forces must be transformed. The equations that correspond with rotations should be multiplied by G¯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 vector q¨ . Therefore, those accelerations should be expressed as a function of the new rotational variables. Taking the time derivative of Eq. (128) gives: Finally, the vector Qˇi can be expressed as: Moreover, the value of the inertial forces can be written as follows: where Mˆ = T T M T , so 6.1 Quadratic velocity vector with Euler parameters Assuming that the vector θ is a set of the Euler parameters, the well-known equalities G˙¯θ˙ = 0, ω˜¯ = 21 G¯ G¯˙T and G¯T G¯ = 4 I − θ θ T can be used. The term G¯T ω˜¯ can be simplified as: G¯T ω˜¯ = 21 G¯T G¯ G¯˙T = 2 I − θ θ T With the equality θ T δθ = 0, the virtual work of the second part of the above equation can be shown to be 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: ⎡ 0 0 0 ⎤ ⎡ mtr G¯˙θ˙ ⎤ T T M ⎣ 0 G˙¯ 0 ⎦ q˙ˆ = ⎢ G¯T mrr G¯˙θ˙ ⎥ 0 0 0 ⎣ m f r G¯˙θ˙ ⎦ Therefore, the final expression of the quadratic velocity vector, which should be used for most rotational parameters, is: equal to zero, i.e., 2 G˙¯T θ θ T δθ = 0. Therefore, Eq. (137) can be written in the simplified form in terms of Euler parameters as: ⎢ = ⎢⎢ ⎣ −2 I5 T ˙¯pω¯ + 2 ¯ˆ 1 ωT 6.2 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 G˙¯θ˙ = 0 and G¯T ω˜¯ = 2 G˙¯T or equivalently −ω¯˜ G¯ = 2 G˙¯ , then the simplifications made for the Euler param eters 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⎤ G˙¯ = ⎣ −ψ˙ sin ψ sin θ + θ˙ cos ψ cos θ −ψ˙ cos ψ 0⎦ . −θ˙ sin θ 0 0 (145) Looking at the first term and Using the notations sψ and cψ to denote, respectively, sinus and co-sinus of the angle ψ : ⎡ φ˙ψ˙ cψ sθ + φ˙θ˙ sψ cθ − θ˙ψ˙ sψ ⎤ = ⎣ −φ˙ψ˙ sψ sθ + φ˙θ˙ cψ cθ − θ˙ψ˙ cψ ⎦ . −φ˙θ˙ sθ And because the parameters θ and the time derivatives θ˙ are generally independent, the above equation cannot be reduced to the zero vector. In addition: ⎡ 0 −φ˙cθ − ψ˙ φ˙cψsθ − θ˙sψ ⎤ ω˜¯ G¯ = ⎣ φ˙cθ + ψ˙ 0 −φ˙sψsθ − θ˙cψ ⎦ −φ˙cψsθ + θ˙sψ φ˙sψsθ + θ˙cψ 0 = ⎣ ⎡ −ψ˙ cψsθ − θ˙sψcθ φ˙cθ sψ + ψ˙ sψ φ˙cψsθ − θ˙sψ ⎤ ψ˙ sψsθ − θ˙cψcθ φ˙cθ cψ + ψ˙ cψ −φ˙sψsθ − θ˙cψ ⎦ θ˙sθ −φ˙sθ 0 ⎡ 0 φ˙cθ sψ φ˙cψsθ − θ˙sψ ⎤ = ⎣ 00 φ˙−cφθ˙csθψ −φ˙sψsθ − θ˙cψ ⎦ − G¯˙. 0 Therefore, G˙¯ is a part of ω˜¯ G¯ , but the first term on the right-hand side of the above equation is not equal to − G¯˙, so it cannot be generally stated that G¯T ω˜¯ is equal to 2 G˙¯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. 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 7 Conclusions 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. Acknowledgements This work was supported, in part, by the Polish National Science Centre under Grant No. DEC2012/07/B/ST8/03993. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 1. Baharudin , M.E. , Korkealaakso , P. , Rouvinen , A. , Mikkola , A. : Crane operators training based on the real-time multibody simulation . In: H. Gattringer , J. Gerstmayr (eds.) Multibody System Dynamics, Robotics and Control , pp. 213 - 229 . Springer, Vienna ( 2013 ) 2. Bauchau , O.A. : Flexible Multibody Dynamics. Springer, New York ( 2011 ) 3. Haug , E.J. : Computer Aided Kinematics and Dynamics of Mechanical Systems: Basic Methods. Prentice Hall , Boston ( 1989 ) 4. Lugrís , U. , Naya , M.A. , Luaces , A. , Cuadrado , J. : Efficient calculation of the inertia terms in floating frame of reference formulations for flexible multibody dynamics . Proc. Inst. Mech. Eng. Part K J. Multibody Dyn . 223 ( 2 ), 147 - 157 ( 2009 ) 5. Nikravesh , P.E. : Computer-Aided Analysis of Mechanical Systems . Prentice Hall, New Jersey ( 1988 ) 6. Shabana , A.A.: Flexible multibody dynamics: review of past and recent developments . Multibody Syst. Dyn . 1 ( 2 ), 189 - 222 ( 1997 ) 7. Shabana , A.A.: Dynamics of Multibody Systems, 4th edn. Cambridge University Press, New York ( 2013 ) 8. Sherif , K. , Nachbagauer , K. : A detailed derivation of the velocity-dependent inertia forces in the floating frame of reference formulation . J. Comput. Nonlin. Dyn . 9 ( 4 ), 044 , 501 ( 2014 ) 9. Wasfy , T.M. , Noor , A.K. : Computational strategies for flexible multibody systems . Appl. Mech. Rev . 56 ( 6 ), 553 - 613 ( 2003 ) 10. Yoo , W.S. , Haug , E.J. : Dynamics of articulated structures. Part I. Theory. J. Struct. Mech . 14 ( 1 ), 105 - 126 ( 1986 )


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs11071-017-3355-y.pdf

Grzegorz Orzechowski, Marko K. Matikainen, Aki M. Mikkola. Inertia forces and shape integrals in the floating frame of reference formulation, Nonlinear Dynamics, 2017, 1953-1968, DOI: 10.1007/s11071-017-3355-y