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

Authors: Maria Luminiţa Scutaru; Sorin Vlase
Journal: Journal of Applied Mathematics
 (2012)

---

## 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.


---

## Body

## 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 (Figure1). 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, Δt periods of time, when the coefficients can be approximate as constant, are being considered. It is difficult to estimate the length of this interval Δt 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]
(1.1)meδ¨e+2ce(ω)δ˙e+[ke+ke(ε)+ke(ω2)+keG]δe=qe+qe*-qei(ε)-qei(ω2)-mEeiIεL-moeiRTr¨oG,
where me: inertia matrix, symmetric; ce(ω): matrix of Coriolis terms, skew symmetric; ke: classical rigidity matrix, symmetric; ke(ε): modification of rigidity due to the angular acceleration of the mobile reference system; ke(ω2): modification of rigidity due to the angular speed of the mobile reference system; keG: rigidity due to the second-order effects; qe+qe*: the external forces concentrated and distributed; qei(ε): inertial forces due to the angular acceleration of the mobile reference system; qei(ω2): inertial forces due to the angular speed of the mobile reference system; mEeiIεL: inertial forces due to rotation only at bars; moeiRTr¨oG: 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)MΔ¨+CΔ˙+[K+K(ε)+K(ω2)]Δ=ϱ,
where M: inertia matrix, symmetric; C: matrix of Coriolis terms, skew symmetric; K: rigidity matrix, symmetric; K(ε): modification of rigidity due to the angular acceleration of the mechanical system components, skew symmetric; K(ω2): modification of rigidity due to the angular speed of the mechanical system components, symmetric; Q: matrix of the total nodal loads.To summarize, matricesM,K,K(ω2) of the whole structure are symmetric and C,K(ε) 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:(2.1)Δ=Acos(ωt+φ).
By derivation, the following is being successively obtained:
(2.2)Δ˙=-ωAsin(ωt+φ),Δ¨=-ω2Acos(ωt+φ).Introducing into the homogeneous system associated to (1.2), we get
(2.3)-ω2MAcos(ωt+φ)-ωCAsin(ωt+φ)+(K+K(ε)+K(ω2))Acos(ωt+φ)=0.By premultiplying (2.3) by A, we get
(2.4)-ω2ATMAcos(ωt+φ)-ωATCAsin(ωt+φ)+AT(K+K(ε)+K(ω2))Acos(ωt+φ)=0.SinceC and K(ε) are skew-symmetric, we have
(2.5)ATCA=0,ATK(ε)A=0,
and (2.4) becomes
(2.6)[-ω2ATMA+AT(K+K(ω2))A]cos(ωt+φ)=0,
where (see Rayleigh [18])
(2.7)ω2=AT[K+K(ω2)]AATMA.
Moreover, these values are real, since K+K(ω2) and M are symmetric; thus there is no damping in the system. It results that the matrix C, 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:
(2.1′)Δ=A1cosωt+A2sinωt.
We can enounce the following property.


### 2.2. P2. MatrixA Defined by Relation (2.11) Has Two Null Eigenvalues

The matricesK(ε) and K(ω2) 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
(2.8)MΔ¨+CΔ˙+KΔ=0,
with M and K symmetric and C skew-symmetric. With the classical substitution
(2.9)Δ=x1,Δ˙=x2,
the first-order linear system is obtained
(2.10){x˙1x˙2}=[0E-M-1K-M-1C]{x1x2}
or
(2.10′)x˙=Ax,
where
(2.11)x={x1x2};A=[0E-M-1K-M-1C],
and the unit matrix was named E.Ifn is the system dimension the dimension of matrix A will be 2n×2n.Letλ1,λ2,…,λ2n be the eigenvalues for matrix A and V1,V2,…,V2n the corresponding eigenvectors. We write
(2.12)Q=[V1V2…V2n].
The matrix of eigenvectors is clear. By transforming x=Qq, the system is reduced to
(2.13)q˙=λq,
where
(2.14)λ=[λ1λ200⋱λ2n],λj∈C.The equationsq˙i=λiqi allow the solutions qi=Ciexp(λit).If in the conditiondet(A-λE)=0 we plug in λ=0, we get the condition det A=0.As a proof, Developing on rows, we get finally(2.15)detA=(-1)ndet(M-1K)=(-1)ndetM-1detK=0,
because K is singular.We will show in what follows thatA allows another eigenvalue 0. If an eigenvalue λ1 is known, matrix A can be reduced to a matrix of dimension (2n-1)×(2n-1) which has 2n-1 eigenvalues equal to the eigenvalues of A where λ1 has been eliminated. It is known [19] that matrix B=T-1AT has the same eigenvalues as A, but the eigenvectors are
(2.16)V′=T-1V.
Then, if a matrix B has the form
(2.17)B=[λ1⋮b12⋯b1n⋯⋮⋯⋯⋯⋯⋮b22⋯b2n0⋮⋮bn2⋯bnn=]=[λ1⋮B12⋯⋯⋯0⋮B22]
with
(2.18)B22=[b22…b2nbn2…bnn];B12=[b12…b1n],
then the problem of eigenvalues
(2.19)det(B22-λE)=0
is to be written (λ-λ1)det(B22-λE)=0 resuming to determine the eigenvalues of B22.It is checked by direct computation that the transformationT is as follows:
(2.20)T=[⋮0⋮⋯⋯⋯⋮V⋮E2n-1⋮⋮]=[1⋮0⋯⋯⋯⋯⋯e2⋮e3⋮E2n-1⋮⋮e2n⋮].By simple computation, we get(2.21)T-1=[1⋮0⋯⋮⋯⋯⋯⋯-e2⋮-e3⋮E2n-1⋮-e2n⋮].It results in the possibility of avoiding the computation of the eigenvalues for matrixA with det A = 0, considering the fact that we know an eigenvalue λ1=0, and, by the transformation presented, we will compute the other eigenvalues as eigenvalues of B22.Forλ1=0, the eigenvector V1 corresponding to matrix A will be computed. We will have
(2.22)AV1=0
or
(2.22′)[⋮0n⋮En⋮⋮⋯⋯⋯⋯⋮⋯⋯⋯⋮-M-1K⋮-M-1C⋮]{e1e2⋮en⋯en1⋮e2n}={00⋮⋮⋮⋮⋮0}.
The n first rows give us en+1=0;en+2=0,…,e2n=0.Using these results, the othern rows give us (2.23)[⋮-M-1K⋮-M-1C⋮]{e1e2⋮en⋯0⋮0}=0.
or, carrying out the multiplications of the matrix blocks,
(2.23′)[-M-1K]{e1e2⋮en}=0.
The result will be that the first n components of the eigenvector are components of the vector corresponding to value λ1=0 for matrix -M-1K which characterizes the system
(2.24)MΔ¨+KΔ=0.
With these components, we build
(2.25)T=[1⋮0⋯⋯⋯⋯⋯e2⋮⋮⋮en⋮⋯⋮E2n-1⋮0⋮⋮],T-1=[1⋮0⋯⋯⋯⋯⋯-e2⋮⋮⋮-en⋮⋯⋮E2n-1⋮0⋮⋮],
where e1=1. We compute B=T-1AT.Consider(2.26)B=[1⋮01x(n-1)⋯⋯⋯⋯⋯-e2⋮⋮⋮-en⋮⋯⋮E2n-1⋮0⋮⋮]·[⋮⋮0n×1⋮En⋮⋮⋯⋮⋯⋮-M-1K⋮-M-1C⋮]·[1⋮0⋯⋯⋯⋯⋯e2⋮⋮⋮en⋮⋯⋮E2n-1⋮0⋮⋮]=[1⋮0⋯⋯⋯⋯⋯-e2⋮⋮⋮-en⋮⋯⋮E2n-1⋮0⋮⋮].[⋮0⋮En⋮⋮⋯⋯⋯⋯⋯0⋮⋮[α*]⋮-M-1C⋮⋮0⋮]=[⋮1⋮-e210n×n⋮-e31⋮⋮⋯⋮-en1⋯⋯⋯⋮⋯⋯⋯⋯⋯0n×1⋮[α*]⋮-M-1C⋮⋮],
where [α*] is the matrix having the dimensions 2n×(2n-1). It is obtained from -M-1K having the dimension 2n×2n by eliminating the first column.We have(2.27)B22=[⋮-e21⋮-e3010⋮⋮⋱⋮-en⋯1⋮⋯⋯⋯⋮⋯⋯⋯⋯⋯[α*]⋮-M-1C⋮]=[⋮-e21⋮-e30⋮⋮⋱⋮-en1⋯⋯⋯⋮⋯⋯⋯⋯⋮[α*]⋮C1C2Cn⋮],
where the matrix columns -M-1C have been named C1,C2,…,Cn.By developing the matrix determinantB22 on rows, we obtain
(2.28)detB22=(-1)ne2[⋮01⋮00⋮⋮⋱⋮01⋯⋯⋯⋮⋯⋯⋯⋯⋮[α*]⋮C1C2Cn⋮]+(-1)n×[⋮-e31⋮-e40⋮⋮⋱⋮-en1⋯⋯⋯⋮⋯⋯⋯⋯⋮[α*]⋮C1C2Cn⋮]=e2det[α*⋮C2]+e3det[α*⋮C3]+⋯+endet[α*⋮Cn]+e1det[α*⋮C1]=det[[α*]⋮e1C1+e2C2+⋯+enCn]
(we took into account that e1=1).The following relation has been considered:(2.29)detAB=detAdetB.
We have
(2.30)1=det[101⋱01]=det([1e101e2⋱01en]·[e10e2⋱0en]).
Taking into account the previous relation, we may write
(2.31)detB22=det([1e101e2⋱01en]·[e10e2⋱0en])det[⋮α*⋮e1C1+e2C2+⋯+enCn⋮]=det[1e101e2⋱01en]·det[αL1*e1⋮e12C11+e1e2C12+⋯+e1enC1nαL2*e2⋮e1e2C21+e22C22+⋯+e2enC2n⋮αLn*en⋮e1enCn1+ene2Cn2+⋯+en2Cn2]=0,
since
(2.32)αL1*e1+αL2*e2+⋯+αLn*en={V1}T[V2V3⋯Vn]=0,(e12c11+e2e1c12+⋯+ene1cn1)+(e2e1c21+e22c22+⋯+e2encn2)+⋯+(ene1cn1+ene2cn2+⋯+en2cnn)=∑i,jeiej(cij+cji)=0,(cij=-cji).
with αL1*,αL2*,…,αLn*, the lines 1,2,…,n of the matrix α* have been named. It results in det B22=0, thus, the matrix B22 has a null eigenvalue. Matrix B will have two null eigenvalues.To these two eigenvalues corresponds the nonharmonic solutionq=C1t+C2 which will represent the rigid motion of the multibody system in a first approximation (C1 and C2 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.


## 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:(2.1)Δ=Acos(ωt+φ).
By derivation, the following is being successively obtained:
(2.2)Δ˙=-ωAsin(ωt+φ),Δ¨=-ω2Acos(ωt+φ).Introducing into the homogeneous system associated to (1.2), we get
(2.3)-ω2MAcos(ωt+φ)-ωCAsin(ωt+φ)+(K+K(ε)+K(ω2))Acos(ωt+φ)=0.By premultiplying (2.3) by A, we get
(2.4)-ω2ATMAcos(ωt+φ)-ωATCAsin(ωt+φ)+AT(K+K(ε)+K(ω2))Acos(ωt+φ)=0.SinceC and K(ε) are skew-symmetric, we have
(2.5)ATCA=0,ATK(ε)A=0,
and (2.4) becomes
(2.6)[-ω2ATMA+AT(K+K(ω2))A]cos(ωt+φ)=0,
where (see Rayleigh [18])
(2.7)ω2=AT[K+K(ω2)]AATMA.
Moreover, these values are real, since K+K(ω2) and M are symmetric; thus there is no damping in the system. It results that the matrix C, 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:
(2.1′)Δ=A1cosωt+A2sinωt.
We can enounce the following property.


## 2.2. P2. MatrixA Defined by Relation (2.11) Has Two Null Eigenvalues

The matricesK(ε) and K(ω2) 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
(2.8)MΔ¨+CΔ˙+KΔ=0,
with M and K symmetric and C skew-symmetric. With the classical substitution
(2.9)Δ=x1,Δ˙=x2,
the first-order linear system is obtained
(2.10){x˙1x˙2}=[0E-M-1K-M-1C]{x1x2}
or
(2.10′)x˙=Ax,
where
(2.11)x={x1x2};A=[0E-M-1K-M-1C],
and the unit matrix was named E.Ifn is the system dimension the dimension of matrix A will be 2n×2n.Letλ1,λ2,…,λ2n be the eigenvalues for matrix A and V1,V2,…,V2n the corresponding eigenvectors. We write
(2.12)Q=[V1V2…V2n].
The matrix of eigenvectors is clear. By transforming x=Qq, the system is reduced to
(2.13)q˙=λq,
where
(2.14)λ=[λ1λ200⋱λ2n],λj∈C.The equationsq˙i=λiqi allow the solutions qi=Ciexp(λit).If in the conditiondet(A-λE)=0 we plug in λ=0, we get the condition det A=0.As a proof, Developing on rows, we get finally(2.15)detA=(-1)ndet(M-1K)=(-1)ndetM-1detK=0,
because K is singular.We will show in what follows thatA allows another eigenvalue 0. If an eigenvalue λ1 is known, matrix A can be reduced to a matrix of dimension (2n-1)×(2n-1) which has 2n-1 eigenvalues equal to the eigenvalues of A where λ1 has been eliminated. It is known [19] that matrix B=T-1AT has the same eigenvalues as A, but the eigenvectors are
(2.16)V′=T-1V.
Then, if a matrix B has the form
(2.17)B=[λ1⋮b12⋯b1n⋯⋮⋯⋯⋯⋯⋮b22⋯b2n0⋮⋮bn2⋯bnn=]=[λ1⋮B12⋯⋯⋯0⋮B22]
with
(2.18)B22=[b22…b2nbn2…bnn];B12=[b12…b1n],
then the problem of eigenvalues
(2.19)det(B22-λE)=0
is to be written (λ-λ1)det(B22-λE)=0 resuming to determine the eigenvalues of B22.It is checked by direct computation that the transformationT is as follows:
(2.20)T=[⋮0⋮⋯⋯⋯⋮V⋮E2n-1⋮⋮]=[1⋮0⋯⋯⋯⋯⋯e2⋮e3⋮E2n-1⋮⋮e2n⋮].By simple computation, we get(2.21)T-1=[1⋮0⋯⋮⋯⋯⋯⋯-e2⋮-e3⋮E2n-1⋮-e2n⋮].It results in the possibility of avoiding the computation of the eigenvalues for matrixA with det A = 0, considering the fact that we know an eigenvalue λ1=0, and, by the transformation presented, we will compute the other eigenvalues as eigenvalues of B22.Forλ1=0, the eigenvector V1 corresponding to matrix A will be computed. We will have
(2.22)AV1=0
or
(2.22′)[⋮0n⋮En⋮⋮⋯⋯⋯⋯⋮⋯⋯⋯⋮-M-1K⋮-M-1C⋮]{e1e2⋮en⋯en1⋮e2n}={00⋮⋮⋮⋮⋮0}.
The n first rows give us en+1=0;en+2=0,…,e2n=0.Using these results, the othern rows give us (2.23)[⋮-M-1K⋮-M-1C⋮]{e1e2⋮en⋯0⋮0}=0.
or, carrying out the multiplications of the matrix blocks,
(2.23′)[-M-1K]{e1e2⋮en}=0.
The result will be that the first n components of the eigenvector are components of the vector corresponding to value λ1=0 for matrix -M-1K which characterizes the system
(2.24)MΔ¨+KΔ=0.
With these components, we build
(2.25)T=[1⋮0⋯⋯⋯⋯⋯e2⋮⋮⋮en⋮⋯⋮E2n-1⋮0⋮⋮],T-1=[1⋮0⋯⋯⋯⋯⋯-e2⋮⋮⋮-en⋮⋯⋮E2n-1⋮0⋮⋮],
where e1=1. We compute B=T-1AT.Consider(2.26)B=[1⋮01x(n-1)⋯⋯⋯⋯⋯-e2⋮⋮⋮-en⋮⋯⋮E2n-1⋮0⋮⋮]·[⋮⋮0n×1⋮En⋮⋮⋯⋮⋯⋮-M-1K⋮-M-1C⋮]·[1⋮0⋯⋯⋯⋯⋯e2⋮⋮⋮en⋮⋯⋮E2n-1⋮0⋮⋮]=[1⋮0⋯⋯⋯⋯⋯-e2⋮⋮⋮-en⋮⋯⋮E2n-1⋮0⋮⋮].[⋮0⋮En⋮⋮⋯⋯⋯⋯⋯0⋮⋮[α*]⋮-M-1C⋮⋮0⋮]=[⋮1⋮-e210n×n⋮-e31⋮⋮⋯⋮-en1⋯⋯⋯⋮⋯⋯⋯⋯⋯0n×1⋮[α*]⋮-M-1C⋮⋮],
where [α*] is the matrix having the dimensions 2n×(2n-1). It is obtained from -M-1K having the dimension 2n×2n by eliminating the first column.We have(2.27)B22=[⋮-e21⋮-e3010⋮⋮⋱⋮-en⋯1⋮⋯⋯⋯⋮⋯⋯⋯⋯⋯[α*]⋮-M-1C⋮]=[⋮-e21⋮-e30⋮⋮⋱⋮-en1⋯⋯⋯⋮⋯⋯⋯⋯⋮[α*]⋮C1C2Cn⋮],
where the matrix columns -M-1C have been named C1,C2,…,Cn.By developing the matrix determinantB22 on rows, we obtain
(2.28)detB22=(-1)ne2[⋮01⋮00⋮⋮⋱⋮01⋯⋯⋯⋮⋯⋯⋯⋯⋮[α*]⋮C1C2Cn⋮]+(-1)n×[⋮-e31⋮-e40⋮⋮⋱⋮-en1⋯⋯⋯⋮⋯⋯⋯⋯⋮[α*]⋮C1C2Cn⋮]=e2det[α*⋮C2]+e3det[α*⋮C3]+⋯+endet[α*⋮Cn]+e1det[α*⋮C1]=det[[α*]⋮e1C1+e2C2+⋯+enCn]
(we took into account that e1=1).The following relation has been considered:(2.29)detAB=detAdetB.
We have
(2.30)1=det[101⋱01]=det([1e101e2⋱01en]·[e10e2⋱0en]).
Taking into account the previous relation, we may write
(2.31)detB22=det([1e101e2⋱01en]·[e10e2⋱0en])det[⋮α*⋮e1C1+e2C2+⋯+enCn⋮]=det[1e101e2⋱01en]·det[αL1*e1⋮e12C11+e1e2C12+⋯+e1enC1nαL2*e2⋮e1e2C21+e22C22+⋯+e2enC2n⋮αLn*en⋮e1enCn1+ene2Cn2+⋯+en2Cn2]=0,
since
(2.32)αL1*e1+αL2*e2+⋯+αLn*en={V1}T[V2V3⋯Vn]=0,(e12c11+e2e1c12+⋯+ene1cn1)+(e2e1c21+e22c22+⋯+e2encn2)+⋯+(ene1cn1+ene2cn2+⋯+en2cnn)=∑i,jeiej(cij+cji)=0,(cij=-cji).
with αL1*,αL2*,…,αLn*, the lines 1,2,…,n of the matrix α* have been named. It results in det B22=0, thus, the matrix B22 has a null eigenvalue. Matrix B will have two null eigenvalues.To these two eigenvalues corresponds the nonharmonic solutionq=C1t+C2 which will represent the rigid motion of the multibody system in a first approximation (C1 and C2 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 matrixC, 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 C 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.



---

*Source: 628503-2012-12-10.xml*