# Decomposition of the Equations of Motion in the Analysis of Dynamics of a 3-DOF Nonideal System

Authors: Jan Awrejcewicz; Roman Starosta; Grażyna Sypniewska-Kamińska
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The dynamic response of a nonlinear system with three degrees of freedom, which is excited by nonideal excitation, is investigated. In the considered system the role of a nonideal source is played by a direct current motor, where the central axis of the rotor is not coincident with the axis of rotation. This translation generates a torque whose magnitude depends on the angular velocity. During the system operation a general coordinate assigned to the nonideal source grows rapidly as a result of rotation. We propose the decomposition of the equations of motion in such a way to extract the solution which is directly related to the rotation of an unbalanced rotor. The remaining part of the solution describes pure oscillation depending on the dynamical behaviour of the whole system. The decomposed equations are solved numerically. The influence of selected system parameters on the rotor vibration is examined. The presented approach can be applied to separate vibration and rotation of motions in many other engineering systems.


---

## Body

## 1. Introduction

The behaviour of mechanical systems subjected to the external excitation periodically changing in time belongs to classical problems of multibody dynamics. When the excitation does not depend on the system response, it is said to be an ideal loading. These problems are widely discussed in the literature. In real problems the motion of a system less or more affects the source of energy, especially in the neighbourhood of resonance. It is caused by a limited power supply to external loading. Such source of energy is called nonideal. The problem of nonideal vibrations of multibody systems leads to a rather sophisticated mathematical description, especially when nonlinearities are taken into account. In that case an additional equation, which describes how the energy source supplies energy to the system, should be added. It causes that an additional degree-of-freedom appears. Therefore, the nonideal system has one degree of freedom more than the ideal counterpart.In our investigations we used results published in a few papers concerning such systems [1–4] and in the references included therein. It should be emphasized that first ideas and modelling concepts of nonideal vibrations were published by Kononenko [5] and the first book entirely devoted to that problem by Sommerfeld [6]. An overview of investigations in this field is described by Balthazar et al. in [7], and hence it is omitted here for the paper conciseness. In that paper the authors reviewed main properties of nonideal vibrating systems. The so-called Sommerfeld effect [1, 7] connected with jump phenomena near the resonance is discussed among other topics. This phenomenon suggests that the vibrational response provides an energy sink. The parametric resonance in the nonideal system with DC motor was analyzed in [1, 2]. The authors of paper [1] used an asymptotic approach to investigate the behaviour of the system near internal resonance.The system investigated in the paper is shown in Figure1. The electric DC motor with an eccentrically mounted rotor is assumed to be a nonideal source of vibrations. In this case the additional degree of freedom is connected to its rotations. The generalized coordinate describing the rotor motion grows in time, so the whole process cannot be considered as vibrational. A new idea proposed in the paper consists in decomposition of the equation of motion related to the rotor and separation of the rotations and vibrations. After this operation a new set of equations of motion is derived where all generalized coordinates describe the vibrations.Figure 1
The investigated system.


## 2. Formulation of the Problem

Let us consider planar motion of the system composed of three elements: the support of massm
1 which can move only vertically, the mathematical pendulum of mass m
2 and length L, and the DC motor. The examined system is shown in Figure 1. It is assumed that the stator of the motor whose mass is equal to m
3 does not move with respect to the base. The rotor of mass m
0 is eccentrically mounted on the axis of rotation and its mass centre lies at the crossing point of the axes of symmetry (as shown in Figure 1). The eccentricity of the rotor is denoted by r
e. The moment of inertia of the rotor about its central axis is equal I
0. The pendulum is suspended from a support by a joint. The support, in turn, is connected to the basis via a viscous damped elastic suspension. Both the spring and the damper are linear. The elastic constant of the spring is denoted by k, whereas the damping coefficient is C
1. Length of the nonstretched spring is L
0. Moreover, at the joint connecting the pendulum with the support, there is viscous damping denoted by C
2. Force F acting in the direction of X-axis and torque M are changing harmonically. Namely, the magnitude of F is given by F
=
F
0
cos
⁡
⁡
(
Ω
1
t
) and M
=
M
0
cos
⁡
⁡
(
Ω
2
t
). Owing to the commonly used characteristic of a DC motor, the produced torque depends linearly on the angular velocity, so we can write that M
p
=
U
1
-
U
1
Λ
˙
(
t
), where U
1 and U
2 depend on the electromagnetic field of the DC motor. Both quantities are constant for each model of the motor considered [1]. Angular velocity Λ
˙
(
t
) oscillates as a result of the rotor unbalancing as well as the interaction with support. Variation of the angular speed influences the torque, so even during stationary operation of the electric motor the torque is not constant.The system has three degrees of freedom (DOFs). Total spring elongationZ
(
t
) and angles Φ
(
t
) and Λ
(
t
) are chosen as generalized coordinates of the system. The kinetic energy relative to the motionless frame and written in the Cartesian coordinate systemOXY is as follows:
(1)
T
=
1
2
m
1
(
X
˙
1
2
+
Y
˙
1
2
)
+
1
2
m
2
(
X
˙
2
2
+
Y
˙
2
2
)
+
1
2
m
3
(
X
˙
3
2
+
Y
˙
3
2
)
+
1
2
m
0
(
X
˙
0
2
+
Y
˙
0
2
)
+
1
2
I
0
Λ
˙
2
.All conservative forces which are acting in the system are represented by the potential energy(2)
V
=
-
m
0
g
X
0
-
m
2
g
X
1
-
m
1
g
X
2
-
m
3
g
X
3
+
1
2
k
Z
2
,
where X
1
(
t
)
=
(
L
0
+
Z
(
t
)
+
h
/
2
), Y
1
=
0 are the coordinates of the movable support, X
3
=
X
2
-
H, Y
3
=
0 are the coordinates of the stator mass centre, X
2
=
X
1
+
L
cos
⁡
⁡
(
Φ
(
t
)
), Y
2
=
L
sin
⁡
(
Φ
(
t
)
) are the coordinates of the suspended pendulum, and X
0
=
X
1
-
H
+
r
e
cos
⁡
⁡
(
Λ
(
t
)
), Y
0
=
r
e
sin
⁡
(
Λ
(
t
)
) are the coordinates of the mass centre of the eccentrically mounted rotor and g is the gravity of Earth.The position of stable equilibrium of the system is determined by the following values of the generalized coordinates:(3)
Z
r
=
(
m
0
+
m
1
+
m
2
+
m
3
)
g
k
,
Φ
r
=
0
,
Λ
r
=
0
.The equations of motion around the equilibrium position are obtained using Lagrange equations of the second type, wherein the nonconservative forces acting on the system are introduced as the generalized forces. They are as follows:(4)
Z
¨
1
+
ω
1
2
Z
1
+
c
1
m
c
Z
˙
1
-
m
0
r
e
m
c
Λ
˙
2
cos
⁡
⁡
Λ
-
L
m
2
m
c
Φ
˙
2
cos
⁡
⁡
Φ
-
m
0
r
e
m
c
Λ
¨
sin
⁡
Λ
-
L
m
2
m
c
Φ
¨
sin
⁡
Φ
-
F
0
m
c
cos
⁡
⁡
(
Ω
1
t
)
=
0
,
Φ
¨
+
ω
2
2
sin
⁡
Φ
+
c
2
L
2
m
2
Φ
˙
-
Z
¨
1
L
sin
⁡
Φ
-
M
0
L
2
m
2
cos
⁡
⁡
(
Ω
2
t
)
=
0
,
(
1
+
m
0
r
e
2
I
0
)
Λ
¨
+
ω
3
2
sin
⁡
Λ
-
U
1
I
0
+
U
2
I
0
Λ
˙
-
m
0
r
e
I
0
Z
¨
1
sin
⁡
Λ
=
0
,
where Z
1
=
Z
-
Z
r is the dynamic elongation of the spring (measured from the equilibrium position), m
c
=
m
0
+
m
1
+
m
2
+
m
3 is the total mass of the system, and ω
1
=
k
/
m
c, ω
2
=
g
/
L, ω
3
=
g
m
0
r
e
/
I
0 are the characteristic frequencies.Let us introduce the dimensionless coordinatez
~
=
Z
1
/
L, time τ
=
t
ω
1, and the following dimensionless parameters:(i)
eccentricityβ
e
=
r
e
/
L,
(ii)
mass fractionsμ
2
=
m
2
/
m
c and μ
0
=
m
0
/
m
c,
(iii)
frequenciesw
2
=
ω
2
/
ω
1, w
3
=
ω
3
/
ω
1, p
1
=
Ω
1
/
ω
1, p
2
=
Ω
2
/
ω
1,
(iv)
amplitudes of external excitationsf
1
=
F
0
/
L
m
c
ω
1
2, f
2
=
M
0
/
L
2
m
2
ω
1
2,
(v)
parameters of the DC motor characteristicsu
1
=
U
1
/
I
0
ω
1
2, u
2
=
U
2
/
I
0
ω
1,
(vi)
damping coefficientsc
1
=
C
1
/
m
c
ω
1, c
2
=
C
2
/
L
2
m
2
ω
1.Then the equations of motion take the nondimensional form:(5)
z
~
¨
+
z
~
+
c
1
z
~
˙
-
β
e
μ
0
λ
~
˙
2
cos
⁡
⁡
λ
~
-
μ
2
φ
~
˙
2
cos
⁡
⁡
φ
~
-
β
e
μ
0
λ
~
¨
sin
⁡
λ
~
-
μ
2
φ
~
¨
sin
⁡
φ
~
-
f
1
cos
⁡
⁡
(
p
1
τ
)
=
0
,
φ
~
¨
+
w
2
2
sin
⁡
φ
~
+
c
2
φ
~
-
z
~
¨
sin
⁡
φ
~
-
f
2
cos
⁡
⁡
(
p
2
τ
)
=
0
,
(
1
+
w
3
2
β
e
w
2
2
)
λ
~
¨
+
w
3
2
sin
⁡
λ
~
-
u
1
+
u
2
λ
~
˙
-
w
3
2
w
2
2
z
~
¨
sin
⁡
λ
~
=
0
.Now functions denoted byz
~, φ
~, λ
~ (equivalent to generalized coordinates Z, Φ, Λ) depend on τ. Equations (5) are supplemented by the initial conditions: z
~
(
0
)
=
z
0, z
~
˙
(
0
)
=
v
0, φ
~
(
0
)
=
φ
0, φ
~
˙
(
0
)
=
ω
0
φ, λ
~
(
0
)
=
λ
0, and λ
~
˙
(
0
)
=
ω
0
λ where z
0, v
0, φ
0, ω
0
φ, λ
0, and ω
0
λ are the known quantities.


## 3. Decomposition of the Governing Equations

The coordinateλ
~
(
τ
) is equal to the angle measured from the initial position. After each rotation of the rotor, the angle increases by 2
π. This increase is dominant in time history of λ
~
(
τ
) which is shown in Figure 2. This graph and several subsequent graphs are made for chosen data included in set SET1 =
{
f
1
=
0, f
2
=
0.01, p
2
=
0.77, c
1
=
0.0001, c
2
=
0.0001, β
e
=
0.001, μ
2
=
0.3, μ
0
=
0.15, u
1
=
0.04, u
2
=
0.02, w
2
=
0.06, w
3
=
0.004, z
0
=
0, v
0
=
0, φ
0
=
0.1, ω
0
φ
=
0, λ
0
=
0.1, ω
0
λ
=
0
}. However, in the time history of λ
~
(
τ
) there are also oscillations yielded by the rotor unbalance. They are invisible in the scale of the graph due to their small amplitude relative to the values of λ
~
(
τ
). Oscillations of the rotor play a significant role in dynamics of the whole system. In order to study the interactions between particular parts of the system, for instance, near resonance, oscillations of the rotor should be separated from its rotation.Figure 2
Time history ofλ
~ for SET1.Thus, it is desirable to split functionλ
~
(
τ
) into a component describing unlimited increase and the second one relative to pure oscillations. We propose decomposition of this function in the following manner:
(6)
λ
~
(
τ
)
=
α
0
(
τ
)
+
α
1
(
τ
)
,
where function α
0
(
τ
) satisfies the initial value problem of the linear differential equation and inhomogeneous initial conditions:
(7)
(
1
+
w
3
2
β
e
w
2
2
)
α
¨
0
-
u
1
+
u
2
α
˙
0
=
0
,
α
0
(
0
)
=
λ
0
,
α
˙
0
(
0
)
=
ω
0
λ
.Decomposition (6) causes that the problem is described now by the system of four differential equations. The Cauchy problem governed by (7) describes the dynamics of the rotation of the rotor under the action of linearly changing torque and can be solved analytically.Substituting (6)-(7) into (5), we obtain a new form of the governing equations with unknown functions z, φ, and α
1:
(8)
z
¨
+
z
+
c
1
z
˙
-
β
e
μ
0
(
α
˙
0
+
α
˙
1
)
2
cos
⁡
⁡
(
α
0
+
α
1
)
-
μ
2
φ
˙
2
cos
⁡
⁡
φ
-
β
e
μ
0
(
α
¨
0
+
α
¨
1
)
sin
(
α
0
+
α
1
)
-
μ
2
φ
¨
sin
⁡
φ
-
f
1
cos
⁡
⁡
(
p
1
τ
)
=
0
,
(9)
φ
¨
+
w
2
2
sin
⁡
φ
+
c
2
φ
˙
-
z
¨
sin
⁡
φ
-
f
2
cos
⁡
⁡
(
p
2
τ
)
=
0
,
(10)
(
1
+
w
3
2
β
e
w
2
2
)
α
¨
1
+
w
3
2
sin
⁡
(
α
0
+
α
1
)
+
u
2
α
˙
1
-
w
3
2
w
2
2
z
¨
sin
⁡
(
α
0
+
α
1
)
=
0
.The initial conditions for this set of differential equations are(11)
z
(
0
)
=
z
0
,
z
˙
(
0
)
=
v
0
,
φ
(
0
)
=
φ
0
,
φ
˙
(
0
)
=
ω
0
φ
,
α
1
(
0
)
=
0
,
α
˙
1
(
0
)
=
0
.Taking into account the introduced decomposition (6), the original initial value problem (5) is transformed into two problems, that is, the problem of pure rotation of the rotor and the problem of pure oscillation in the whole system. In the first case the term “pure” is used to emphasize that rotation α
0
(
τ
) does not depend on the rotor unbalancing or interactions between system parts. In the second case, it is used to note that oscillations have been completely separated from rotations. Let us note that in equations (8)–(10) function α
0
(
τ
) is a solution of the problem (7), and it has the following form:
(12)
α
0
=
-
u
1
u
2
2
+
ω
0
λ
u
2
+
(
(
(
u
2
ω
0
λ
-
u
1
)
(
w
3
2
β
e
-
(
w
2
2
+
w
3
2
β
e
)
exp
⁡
(
-
u
2
w
2
2
τ
w
2
2
+
w
3
2
β
e
)
)
×
(
w
3
2
β
e
-
(
w
2
2
+
w
3
2
β
e
)
exp
⁡
(
-
u
2
w
2
2
τ
w
2
2
+
w
3
2
β
e
)
)
)
×
(
u
2
2
w
2
2
)
-
1
×
(
w
3
2
β
e
-
(
w
2
2
+
w
3
2
β
e
)
exp
⁡
(
-
u
2
w
2
2
τ
w
2
2
+
w
3
2
β
e
)
)
)
)
+
λ
0
+
u
1
τ
u
2
.It should be noticed that functionα
1
(
τ
) governs oscillations of the rotor around its rotation described by function α
0
(
τ
). In (8) and (10), trigonometric functions appear yielded by sum α
0
(
τ
)
+
α
1
(
τ
). We expand them using the classical trigonometric identities, which clearly show that the coefficients in differential equations (8) and (10) are the functions of time.Equations (8)–(10), yielded by the proposed decomposition of the original problem, describe pure oscillation of the system parts. The associated time histories of z, φ, and α
1 are presented in Figures 3, 4, and 5.Figure 3
Time history ofz for data from SET1.Figure 4
Time history ofφ for data from SET1.Figure 5
Time history ofα
1 for data from SET1.Coincidence of the solutions of the initial problems (5) and (8)–(10) together with the solution (12) is confirmed based on the results presented in Figure 6 (Δ
z
=
z
~
-
z, Δ
φ
=
φ
~
-
φ, Δ
α
=
λ
~
-
(
α
0
+
α
1
)), which validates our decomposition proposal.Identity of solutions:z
~ versus z, φ
~ versus φ, and λ
~ versus (
α
0
+
α
1
) for data from SET1.
(a)
(b)
(c)Very high conformity of solutions of the initial problems (5) and (8)–(10) has been achieved numerically for various parameters, despite the fact that the postulated form of solution (6) is not a general solution of the problem (5).Functionα
0
(
τ
) for data set SET1 is presented graphically in Figure 7. Let us observe that when time τ tends to infinity, function α
0
(
τ
) approaches its asymptote given by
(13)
lim
⁡
τ
→
∞
(
α
0
(
τ
)
)
=
-
u
1
u
2
2
+
ω
0
λ
u
2
+
(
u
2
ω
0
λ
-
u
1
)
w
3
2
β
e
u
2
2
w
2
2
+
λ
0
+
u
1
τ
u
2
=
α
^
0
(
τ
)
.Figure 7
Asymptotic convergence ofα
0 and α
^
0 for data from SET1.Therefore, the exponential transient componentexp
⁡
⁡
(
-
(
(
u
2
w
2
2
τ
)
/
(
w
2
2
+
w
3
2
β
e
)
)
), which disappears in time, can be omitted, apart from the initial stage of motion. Owing to this observation we can approximate solution (8) in the following manner:
(14)
α
0
(
τ
)
≈
α
^
0
(
τ
)
.


## 4. Properties of Functionα
1

Functionα
1
(
τ
) describes, in principle, the oscillation of the nonideal source (beyond the transitional period), being an oscillation around the rotational motion. This effect is much more spectacular for generalized velocities. Time histories of derivative α
˙
0
(
τ
) and the sum of derivatives α
˙
0
(
τ
)
+
α
˙
1
(
τ
) are shown in Figure 8.Figure 8
Time histories ofα
˙
0 (black line) and (
α
˙
0
+
α
˙
1
) (red line) for data from SET1.The character of vibrations described by functionα
1
(
τ
) strongly depends on the values of mechanical parameters. Depending on the parameters used, the system oscillation can be either quasiperiodic or chaotic. For the values of parameters collected in SET1 they are quasiperiodic. In order to detect the frequencies of α
1
(
τ
), the Fourier analysis was performed.The used sampling covered time interval (0, 2500) with the step equal to 0.005 of dimensionless time units. Figure9 shows a few peaks corresponding to the dimensionless frequencies equal to 1, 2, 3, and 4. The first peak represents the frequency of support oscillation. The largest peak corresponds to the value of slope u
1
/
u
2 of asymptote α
^
0
(
τ
) given by (13). The remaining ones correspond to harmonics of the higher order.Figure 9
Discrete Fourier transform ofα
1 for SET1.Figure10 displays how the interval between two neighbouring extrema of function α
1 changes in time. Four regular series of the “period,” which is determined in this way, can be observed. The series are drawn in different colours. Each such series consists of many points. The ordinate of each of these points is equal to the distance between two instants at which function α
1 achieves its adjacent extrema. The abscissa represents the instant at which the first extremum of the given pair occurs. Time on the horizontal axis (in Figure 10) is measured from the start of the simulation process. Apart from the beginning of the motion, we can observe that the series oscillate around the value Δ
T
≈
π, which is equal to 2
π
/
(
u
1
/
u
2
). Results were performed for the following set of the values of parameters SET2 =
{
f
1
=
0, f
2
=
0, c
1
=
0.0001, c
2
=
0.0001, β
e
=
0.001, μ
2
=
0.3, μ
0
=
0.15, u
1
=
0.04, u
2
=
0.02, w
2
=
0.06, w
3
=
0.004, z
0
=
0, v
0
=
0, φ
0
=
0.1, and ω
0
φ
=
0, λ
0
=
0.1, ω
0
λ
=
0
}. A fragment of the time history of function α
1, corresponding to Figure 10, is shown in Figure 11.Figure 10
Time interval between two adjacent maxima of functionα
1 for SET2.Figure 11
Time history ofα
1 for SET1.The character of change of “period”Δ
T depends, among others, on the mass of the pendulum, and the dependence is illustrated in Figure 12. Calculations were carried out assuming the following values of the mass fraction: μ
2
=
0.05, μ
2
=
0.1, μ
2
=
0.2, and μ
2
=
0.4. All other parameters for the simulations are the same as in the SET2. The values of the interval between the adjacent extrema of function α
1 change very slowly when the mass fraction of the pendulum is small. The rate of these changes clearly accelerates with the increasing mass of the pendulum.Time interval between two adjacent maxima of functionα
1 for various values of μ
2.
(a)
(b)
(c)
(d)The changing “period”Δ
T behaves less regularly when the system is exposed to other external loadings. Figure 13 shows the results of numerical simulations in which the action of force F and torque M is taken into account. The calculations were obtained for the following set of the values of parameters SET3 =
{
f
1
=
0, f
2
=
0.01, p
2
=
0.77, c
1
=
0.0001, c
2
=
0.0001, β
e
=
0.002, μ
2
=
0.3, μ
0
=
0.15, u
1
=
0.04, u
2
=
0.02, w
2
=
0.06, w
3
=
0.004, z
0
=
0, v
0
=
0, φ
0
=
0.1, ω
0
φ
=
0, λ
0
=
0.1, and ω
0
λ
=
0
}.Figure 13
Time interval between two adjacent maxima of functionα
1 for SET3.Although four series are still observable, however, they are not so regular as in the case when only the DC motor excited oscillations. Point clouds forming each of these series are more fuzzy. Moreover, the effect of external loading causes that the regions of irregular changes of “period”Δ
T appear. In the case presented in Figure 13, such a irregular behaviour is observed for a few time intervals. The first of them can be regarded as a transient state. In the middle of the second one a gap in the point cloud is observable. It means that the values of Δ
T are bigger than the values of ordinates on the vertical axis. The frequency of vibrations of the unbalanced rotor is then smaller than the frequency determined by slope u
1
/
u
2. The next interval is marked in Figure 13 as (
τ
1
,
τ
2
). The degree of irregularity in this range is smaller than the one studied previously. Fragments of the time history of function α
1, corresponding to Figure 13, are shown in Figure 14.Time interval between two adjacent maxima of functionα
1 for SET3.
(a)
(b)The graph on the left-hand side shows the time history of coordinateα
1
(
τ
) when changes of Δ
T are irregular (for τ
∈
(
τ
1
,
τ
2
)). A more regular motion is presented in the graph on the right-hand side (for τ
>
τ
2).


## 5. Conclusions

The nonlinear system with three degrees of freedom and excited by the DC motor with the unbalanced rotor was studied both analytically and numerically. Many novel nonlinear effects were detected, illustrated, and described. First, the general coordinate assigned to the nonideal source grew rapidly as a result of the rotation of the rotor. The oscillation caused by interaction with the excited mechanical system was imperceptible in that scale. Second, the decomposition of the system of equations of motion was proposed in order to separate the infinitely growing component. Though the decomposition carried out was an identity in the mathematical sense, it was not completely effective in each studied case. However, it allowed the rotor oscillation to be examined separately which essentially simplified the analysis.Third, the equations of motion, decomposed in this way, were solved numerically. Properties of the solution governing pure oscillation of unbalanced rotor were analyzed. The character of this phenomenon strongly depended on the values of parameters of the system. In some cases the quasiperiodic oscillations were exhibited and their main frequency oscillated around the value equal to the slope of the asymptote of functionα
0 representing a pure rotation.



---

*Source: 816840-2014-08-05.xml*