Some Properties of Motion Equations Describing the Nonlinear Dynamical Response of a Multibody System with Flexible Elements

Abstract

The industrial applications use instruments and machines operating at high speeds, developing high forces, low temperatures, corrosive environments, extreme pressures, and so forth. Under these conditions, the elasticity of elements such a machine is built of cannot be ignored anymore, and models are needed to more accurately “grasp” the mechanical phenomena accompanying the operation. The vibrations and the loss of stability are the main effects occurring under these conditions. For the study on this kind of systems with rigid motion and elastic elements, numerous models have been elaborated, the main idea being the discretization of the elements and the use of finite element method. Finally, second-order differential equations with variable coefficients are obtained; these equations are strong nonlinear ones due to the time-dependent values of angular speed and acceleration, and they can be linearized considering a very short period of time, in which the motion is considered to be “frozen.” The aim of this paper is to present some characteristic properties of these systems.

1. Introduction

A mechanical system, a machine or instrument, is made up of elastic elements, the elasticity manifesting itself more or less. The rigid elements assumption generally made when studying such technical systems represents a first approximation leading to rapid results closer to reality. Depending on the given instrument operation conditions, this assumption may lead to correct results or to results considerably deviating from the real situation. If the instrument or the machine works with low operation speeds or if it is subjected to lower loads, then the model built based on the rigid elements assumption may lead to excellent results. The elasticity becomes a significant element if the loads occurred are high and/or the operation speeds are high (Figure 1). In this case, the deformations of the machine element will influence, usually in a negative way, the correct operation of the system. The resonance and the loss of stability represent the main forms of manifestation of elasticity. They will occur in the case of an inadequate design leading to a fast machine damage. The main method of approach of such a problem is the method of finite elements, a method used in a lot of works for elaborating models describing the behavior of machines containing elastic elements.

Figure 1 

Mechanical system with elastic elements [1].

As a result of modeling a set of differential equations having a complex nonlinear form is obtained. It contains a series of additional terms due to the relative movements of nodal coordinates and the general motion, and the operation motion is called “rigid motion” of the system. As the geometrical configuration of a mechanical system with rigid motion changes from one moment to the other, the motion equations are valid for that very moment, in which the system motion is to be considered to be “frozen”. In this case, Δ𝑡 periods of time, when the coefficients can be approximate as constant, are being considered. It is difficult to estimate the length of this interval Δ𝑡 that depends on the topology and geometry of the multibody system. We consider that it is short enough in order to regard the motion as “frozen.” The analysis made is an incremental analysis. The motion equations written in the local coordinate system may be expressed under the form [1–9] 𝐦_𝐞𝜹̈_𝐞+2𝐜_𝐞(𝜔)𝜹̇_𝐞+𝐤_𝐞+𝐤_𝐞(𝜀)+𝐤_𝐞𝜔²+𝐤^𝐆_𝐞𝜹_𝐞=𝐪_𝐞+𝐪^∗_𝐞−𝐪^𝐢_𝐞(𝜀)−𝐪^𝐢_𝐞𝜔²−𝐦^𝐢_𝐄𝐞𝐈𝜺_𝐋−𝐦^𝐢_𝐨𝐞𝐑^𝐓𝐫̈_𝐨𝐆, (1.1) where 𝐦_𝐞: inertia matrix, symmetric; 𝐜_𝐞(𝜔): matrix of Coriolis terms, skew symmetric; 𝐤_𝐞: classical rigidity matrix, symmetric; 𝐤_𝐞(𝜀): modification of rigidity due to the angular acceleration of the mobile reference system; 𝐤_𝐞(𝜔²): modification of rigidity due to the angular speed of the mobile reference system; 𝐤^𝐆_𝐞: rigidity due to the second-order effects; 𝐪_𝐞+𝐪^∗_𝐞: the external forces concentrated and distributed; 𝐪^𝐢_𝐞(𝜀): inertial forces due to the angular acceleration of the mobile reference system; 𝐪^𝐢_𝐞(𝜔²): inertial forces due to the angular speed of the mobile reference system; 𝐦^𝐢_𝐄𝐞𝐈𝜺_𝐋: inertial forces due to rotation only at bars; 𝐦^𝐢_𝐨𝐞𝐑^𝐓𝐫̈_𝐨𝐆: inertial forces due to the shifting of the mobile reference system.

After passing to a common global coordinate system (neglecting the second-order effects), finally, the second-order differential equations system with variable coefficients is obtained [10–12] 𝐌Δ̈+𝐂Δ̇+𝐊+𝐊(𝜀)+𝐊𝜔²Δ=𝝔, (1.2) where 𝐌: inertia matrix, symmetric; 𝐂: matrix of Coriolis terms, skew symmetric; 𝐊: rigidity matrix, symmetric; 𝐊(𝜀): modification of rigidity due to the angular acceleration of the mechanical system components, skew symmetric; 𝐊(𝜔²): modification of rigidity due to the angular speed of the mechanical system components, symmetric; 𝐐: matrix of the total nodal loads.

To summarize, matrices 𝐌,𝐊,𝐊(𝜔²) of the whole structure are symmetric and 𝐂,𝐊(𝜀) are skew symmetric. New researches increase the complexity of the problem of the study of the multibody systems with flexible elements [13–17], but there are not many useful results concerning solving the equations. Some properties of such a system will be presented later.

2. Properties of Motion Equations of the Mechanical System

The motion equations of such system have properties allowing an easier solution of the equation system obtained but also a qualitative interpretation of the dynamic response of the system. We are presenting these properties as follows.

2.1. P1. In the Rayleigh Quotient, the Eigenvalues Do Not Depend Directly on the Damping Matrix

We consider a very short period of time in which the variation of the matrix coefficients of motion equations is of no importance (very low). The system solution of the following form is selected: Δ=𝐀cos𝜔𝑡+𝜑. (2.1) By derivation, the following is being successively obtained: Δ̇=−𝜔𝐀sin𝜔𝐭+𝜑,Δ̈=−𝜔²𝐀cos𝜔𝐭+𝜑. (2.2)

Introducing into the homogeneous system associated to (1.2), we get −𝜔²𝐌𝐀cos𝜔𝐭+𝜑−𝜔𝐂𝐀sin𝜔𝐭+𝜑+𝐊+𝐊(𝜀)+𝐊𝜔²𝐀cos𝜔𝐭+𝜑=0. (2.3)

By premultiplying (2.3) by 𝐀, we get −𝜔²𝐀^𝐓𝐌𝐀cos𝜔𝐭+𝜑−𝜔𝐀^𝐓𝐂𝐀sin𝜔𝐭+𝜑+𝐀^𝐓𝐊+𝐊(𝜀)+𝐊𝜔²𝐀cos𝜔𝐭+𝜑=0. (2.4)

Since 𝐂 and 𝐊(𝜀) are skew-symmetric, we have 𝐀^𝐓𝐂𝐀=0,𝐀^𝐓𝐊(𝜀)𝐀=0, (2.5) and (2.4) becomes −𝜔²𝐀^𝐓𝐌𝐀+𝐀^𝐓𝐊+𝐊𝜔²𝐀cos𝜔𝐭+𝜑=0, (2.6) where (see Rayleigh [18]) 𝜔²=𝐀^𝐓𝐊+𝐊𝜔²𝐀𝐀^𝐓.𝐌𝐀 (2.7) Moreover, these values are real, since 𝐊+𝐊(𝜔²) and 𝐌 are symmetric; thus there is no damping in the system. It results that the matrix 𝐂, having in equations the significance of a viscous damping, does not introduce a damping term in the obtained solutions. Its skew symmetry denotes the system property not dissipating the energy to the exterior.

Let us now write the relation (2.1) as follows: Δ=𝐀₁cos𝜔𝐭+𝐀₂sin𝜔𝐭. 2.1^′ We can enounce the following property.

2.2. P2. Matrix A Defined by Relation (2.11) Has Two Null Eigenvalues

The matrices 𝐊(𝜀) and 𝐊(𝜔²) are neglected; they have generally a small contribution to the rigidity matrix. We propose ourselves to solve the homogeneous system associated to (1.2) which becomes under this assumption 𝐌Δ̈+𝐂Δ̇+𝐊Δ=0, (2.8) with 𝐌 and 𝐊 symmetric and 𝐂 skew-symmetric. With the classical substitution Δ=𝐱₁,Δ̇=𝐱₂, (2.9) the first-order linear system is obtained 𝐱̇₁𝐱̇₂=0𝐄−𝐌⁻¹𝐊−𝐌⁻¹𝐂𝐱₁𝐱₂ (2.10) or 𝐱̇=𝐀𝐱, 2.10^′ where 𝐱𝐱=¹𝐱₂;𝐀=0𝐄−𝐌⁻¹𝐊−𝐌⁻¹𝐂, (2.11) and the unit matrix was named 𝐄.

If 𝑛 is the system dimension the dimension of matrix 𝐀 will be 2𝑛×2𝑛.

Let 𝜆₁,𝜆₂,…,𝜆_2𝑛 be the eigenvalues for matrix 𝐀 and 𝐕₁,𝐕₂,…,𝐕_2𝐧 the corresponding eigenvectors. We write 𝐐=𝐕₁𝐕₂…𝐕_2𝐧. (2.12) The matrix of eigenvectors is clear. By transforming 𝐱=𝐐𝐪, the system is reduced to 𝐪̇=𝜆𝐪, (2.13) where 𝜆𝜆=¹𝜆₂0𝜆0⋱_2𝐧,𝜆_𝑗∈𝐶. (2.14)

The equations 𝑞̇_𝑖=𝜆_𝑖𝑞_𝑖 allow the solutions 𝑞_𝑖=𝐶_𝑖exp(𝜆_𝑖𝑡).

If in the condition det(𝐀−𝜆𝐄)=0 we plug in 𝜆=0, we get the condition det 𝐀=0.

As a proof, Developing on rows, we get finally det𝐀=(−1)^𝐧det𝐌⁻¹𝐊=(−1)^𝐧det𝐌⁻¹det𝐊=0, (2.15) because 𝐊 is singular.

We will show in what follows that 𝐀 allows another eigenvalue 0. If an eigenvalue 𝜆₁ is known, matrix 𝐀 can be reduced to a matrix of dimension (2𝑛−1)×(2𝑛−1) which has 2𝑛−1 eigenvalues equal to the eigenvalues of 𝐀 where 𝜆₁ has been eliminated. It is known [19] that matrix 𝐁=𝐓⁻¹𝐀𝐓 has the same eigenvalues as 𝐀, but the eigenvectors are 𝐕^′=𝐓⁻¹𝐕. (2.16) Then, if a matrix 𝐁 has the form 𝜆𝐁=¹⋮𝑏₁₂⋯𝑏_1𝑛⋯⋮⋯⋯⋯⋯⋮𝑏₂₂⋯𝑏_2𝑛0⋮⋮𝑏_𝑛2⋯𝑏_𝑛𝑛==𝜆₁⋮𝐁₁₂⋯⋯⋯0⋮𝐁₂₂ (2.17) with 𝐁₂₂=𝑏₂₂…𝑏_2𝑛𝑏_𝑛2…𝑏_𝑛𝑛;𝐁₁₂=𝑏₁₂…𝑏_1𝑛, (2.18) then the problem of eigenvalues det𝐁₂₂−𝜆𝐄=0 (2.19) is to be written (𝜆−𝜆₁)det(𝐁₂₂−𝜆𝐄)=0 resuming to determine the eigenvalues of 𝐁₂₂.

It is checked by direct computation that the transformation 𝐓 is as follows: ⋮𝐓=⋮0⋮⋯⋯⋯𝐕⋮𝐄_2𝐧−1⋮⋮=𝐞1⋮0⋯⋯⋯⋯⋯₂⋮𝐞₃⋮𝐄_2𝐧−1𝐞⋮⋮_2𝐧⋮. (2.20)

By simple computation, we get 𝐓⁻¹=1⋮0⋯⋮⋯⋯⋯⋯−𝐞₂⋮−𝐞₃⋮𝐄_2𝐧−1⋮−𝐞_2𝐧⋮. (2.21)

It results in the possibility of avoiding the computation of the eigenvalues for matrix 𝐀 with det 𝐀 = 0, considering the fact that we know an eigenvalue 𝜆₁=0, and, by the transformation presented, we will compute the other eigenvalues as eigenvalues of 𝐁₂₂.

For 𝜆₁=0, the eigenvector 𝐕₁ corresponding to matrix 𝐀 will be computed. We will have 𝐀𝐕₁=0 (2.22) or ⋮0_𝐧⋮𝐄_𝐧⋮⋮⋮⋯⋯⋯⋯⋮⋯⋯⋯−𝐌⁻¹𝐊⋮−𝐌⁻¹𝐂⋮𝐞₁𝐞₂⋮𝐞_𝑛⋯𝐞_𝑛1⋮𝐞_2𝐧=00⋮⋮⋮⋮⋮0. 2.22^′ The 𝑛 first rows give us 𝑒_𝑛+1=0;𝑒_𝑛+2=0,…,𝑒_2𝑛=0.

Using these results, the other 𝑛 rows give us ⋮−𝐌⁻¹𝐊⋮−𝐌⁻¹𝐂⋮𝐞₁𝐞₂⋮𝐞_𝑛⋯0⋮0=0. (2.23) or, carrying out the multiplications of the matrix blocks, −𝐌⁻¹𝐊𝐞₁𝐞₂⋮𝐞_𝐧=0. 2.23^′ The result will be that the first 𝑛 components of the eigenvector are components of the vector corresponding to value 𝜆₁=0 for matrix −𝐌⁻¹𝐊 which characterizes the system 𝐌Δ̈+𝐊Δ=0. (2.24) With these components, we build 𝐞𝐓=1⋮0⋯⋯⋯⋯⋯₂⋮𝐞⋮⋮_𝐧⋮⋯⋮𝐄_2𝐧−1⋮⋮0⋮,𝐓⁻¹=1⋮0⋯⋯⋯⋯⋯−𝐞₂⋮⋮⋮−𝐞_𝐧⋮⋯⋮𝐄_2𝐧−1⋮⋮,0⋮ (2.25) where 𝑒₁=1. We compute 𝐁=𝐓⁻¹𝐀𝐓.

Consider 𝐁=1⋮0_{1𝐱(𝐧−1)}⋯⋯⋯⋯⋯−𝐞₂⋮⋮⋮−𝐞_𝐧⋮⋯⋮𝐄_2𝐧−1⋮⋮⋅⋮⋮00⋮^𝐧×1⋮𝐄_𝐧⋮⋮⋮⋯⋮⋯−𝐌⁻¹𝐊⋮−𝐌⁻¹𝐂⋮⋅𝐞1⋮0⋯⋯⋯⋯⋯₂⋮𝐞⋮⋮_𝐧⋮⋯⋮𝐄_2𝐧−1⋮⋮=0⋮1⋮0⋯⋯⋯⋯⋯−𝐞₂⋮⋮⋮−𝐞_𝐧⋮⋯⋮𝐄_2𝐧−1⋮⋮.⋮0⋮0⋮𝐄_𝐧⋮⋮⋮𝛼⋯⋯⋯⋯⋯0⋮_∗⋮−𝐌⁻¹𝐂=⋮⋮0⋮⋮1⋮−𝐞₂10_𝐧×𝐧⋮−𝐞₃1⋮⋮⋯⋮−𝐞_𝐧10⋯⋯⋯⋮⋯⋯⋯⋯⋯_𝐧×1⋮𝛼^∗⋮−𝐌⁻¹𝐂,⋮⋮ (2.26) where [𝛼^∗] is the matrix having the dimensions 2𝑛×(2𝑛−1). It is obtained from −𝐌⁻¹𝐊 having the dimension 2𝑛×2𝑛 by eliminating the first column.

We have 𝐁₂₂=⋮−𝐞₂1⋮−𝐞₃010⋮⋮⋱⋮−𝐞_𝐧⋮𝛼⋯1⋯⋯⋯⋮⋯⋯⋯⋯⋯_∗⋮−𝐌⁻¹𝐂⋮=⋮−𝐞₂1⋮−𝐞₃0⋮⋮⋱⋮−𝐞_𝐧1⋮𝛼⋯⋯⋯⋮⋯⋯⋯⋯_∗⋮𝐂₁𝐂₂𝐂_𝐧⋮, (2.27) where the matrix columns −𝐌⁻¹𝐂 have been named 𝐂₁,𝐂₂,…,𝐂_𝐧.

By developing the matrix determinant 𝐁₂₂ on rows, we obtain det𝐁₂₂=(−1)^𝐧𝐞₂⋮𝛼⋮01⋮00⋮⋮⋱⋮01⋯⋯⋯⋮⋯⋯⋯⋯_∗⋮𝐂₁𝐂₂𝐂_𝐧⋮+(−1)^𝐧×⋮−𝐞₃1⋮−𝐞₄0⋮⋮⋱⋮−𝐞_𝐧1⋮𝛼⋯⋯⋯⋮⋯⋯⋯⋯_∗⋮𝐂₁𝐂_𝟐𝐂_𝐧⋮=𝑒₂det𝛼^∗⋮𝐶₂+𝑒₃det𝛼^∗⋮𝐶₃+⋯+𝑒_𝑛det𝛼^∗⋮𝐶_𝑛+𝑒₁det𝛼^∗⋮𝐶₁=det𝛼^∗⋮𝐞₁𝐂₁+𝐞₂𝐂₂+⋯+𝐞_𝐧𝐂_𝐧 (2.28) (we took into account that 𝑒₁=1).

The following relation has been considered: det𝐀𝐁=det𝐀det𝐁. (2.29) We have 1⋱11=det1001=det𝑒₁01𝑒₂⋱01𝑒_𝑛⋅𝑒₁0𝑒₂⋱0𝑒_𝑛. (2.30) Taking into account the previous relation, we may write det𝐁₂₂1=det𝑒₁01𝑒₂⋱01𝑒_𝑛⋅𝑒₁0𝑒₂⋱0𝑒_𝑛⋮det𝛼^∗⋮𝐞₁𝐂₁+𝐞₂𝐂₂+⋯+𝐞_𝐧𝐂_𝐧⋮1=det𝑒₁01𝑒₂⋱01𝑒_𝑛𝛼⋅det^∗𝐋1𝐞₁⋮𝐞²₁𝐂₁₁+𝐞₁𝐞₂𝐂₁₂+⋯+𝐞₁𝐞_𝐧𝐂_1𝐧𝛼^∗_𝐋2𝐞₂⋮𝐞₁𝐞₂𝐂₂₁+𝐞²₂𝐂₂₂+⋯+𝐞₂𝐞_𝐧𝐂_2𝐧⋮𝛼^∗_𝐋𝐧𝐞_𝐧⋮𝐞₁𝐞_𝐧𝐂_𝐧1+𝐞_𝐧𝐞₂𝐂_𝐧2+⋯+𝐞²_𝐧𝐂_𝐧2=0, (2.31) since 𝛼^∗_𝐋1𝐞₁+𝛼^∗_𝐋2𝐞₂+⋯+𝛼^∗_𝐋𝐧𝐞_𝐧=𝐕₁^𝐓𝐕₂𝐕₃⋯𝐕_𝐧𝑒=0,₂₁𝑐₁₁+𝑒₂𝑒₁𝑐₁₂+⋯+𝑒_𝑛𝑒₁𝑐_𝑛1+𝑒₂𝑒₁𝑐₂₁+𝑒²₂𝑐₂₂+⋯+𝑒₂𝑒_𝑛𝑐_𝑛2+⋯+𝑒_𝑛𝑒₁𝑐_𝑛1+𝑒_𝑛𝑒₂𝑐_𝑛2+⋯+𝑒²_𝑛𝑐_𝑛𝑛=_𝑖,𝑗𝑒_𝑖𝑒_𝑗𝑐_𝑖𝑗+𝑐_𝑗𝑖=0,𝑐_𝑖𝑗=−𝑐_𝑗𝑖. (2.32) with 𝛼^∗_𝐿1,𝛼^∗_𝐿2,…,𝛼^∗_𝐿𝑛, the lines 1,2,…,𝑛 of the matrix 𝛼^∗ have been named. It results in det 𝐁₂₂=0, thus, the matrix 𝐁₂₂ has a null eigenvalue. Matrix 𝐁 will have two null eigenvalues.

To these two eigenvalues corresponds the nonharmonic solution 𝑞=𝐶₁𝑡+𝐶₂ which will represent the rigid motion of the multibody system in a first approximation (𝐶₁ and 𝐶₂ are integration constants). The other values are different from null for a multibody system with only one degree of mobility generally being different from each other.

3. Conclusions

The paper presents a few properties of motion equations in the case of mechanical systems having elastic elements. The properties are due to the existence of the skew-symmetric matrix 𝐂, by which the relative motion of nodal coordinates is manifested in case of applying the finite element method by the Coriolis effects. These properties allow a qualitative analysis of the obtained motion equations. Thus, the Coriolis effects due to the relative motions will determine a modification (generally small) of the systems eigenvalues. The Coriolis damping is not dissipative, meaning that the systems energy is not being influenced by the terms in which the skew-symmetric matrix 𝐂 occurs. The last conclusion is that in case of modeling with the finite element method, the rigid motion, considered in a first approximation an a uniform motion of the system, can be eliminated from the motion equations written.

This fact also suggests the incremental solution of the problem on small periods of time in which the motion may be considered as “frozen” or uniform.

Acknowledgments

This paper is supported by the Sectoral Operational Programme Human Resources Development (SOP HRD), financed from the European Social Fund, and by the Romanian Government under contract no. POSDRU/89/1.5/S/59323.