Interdisciplinary Description of Complex Systems 6(2), 117-131, 2008

RESPONSE AND DYNAMICAL STABILITY

OF OSCILLATORS WITH DISCONTINUOUS
OR STEEP FIRST DERIVATIVE
OF RESTORING CHARACTERISTIC
Hinko Wolf1, * Dubravko Banić2 and Željko Božić1

1
Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb
Zagreb, Croatia
2
Faculty of Graphic Arts, University of Zagreb
Zagreb, Croatia

Regular article Received: 2. August 2008. Accepted: 8. January 2009.

ABSTRACT
Response and dynamical stability of oscillators with discontinuous or steep first derivative of
restoring characteristic is considered in this paper. For that purpose, a simple single-degree-of-
freedom system with piecewise-linear force-displacement relationship subjected to a harmonic force
excitation is analysed by the method of piecing the exact solutions (MPES) in the time domain and by
the incremental harmonic balance method (IHBM) in the frequency domain. The stability of the
periodic solutions obtained in the frequency domain by IHBM is estimated by the Floquet-Lyapunov
theorem. Obtained frequency response characteristic is very complex and includes multi-frequency
response for a single frequency excitation, jump phenomenon, multi-valued and non-periodic
solutions. Determining of frequency response characteristic in the time domain by MPES is
exceptionally time consuming, particularly inside the frequency ranges of co-existence of multiple
stable solutions. In the frequency domain, IHBM is very efficient and very well suited for obtaining
wide range frequency response characteristics, parametric studies and bifurcation analysis. On the
other hand, neglecting of very small harmonic terms (which in-significantly influence the r.m.s.
values of the response and are very small in comparison to other terms of the spectrum) can cause
very large error in evaluation of the eigenvalues of the monodromy matrix, and so they can lead to
incorrect prediction of the dynamical stability of the solution. Moreover, frequency ranges are
detected inside which the procedure of evaluation of eigenvalues of the monodromy matrix does not
converge with increasing the number of harmonics included in the supposed approximate solution.

KEY WORDS
dynamical stability, response characteristic, non-linear vibrations, piecewise-linear system

CLASSIFICATION
ACM: J.2 Engineering
JEL: Z0

*Corresponding author, η: hwolf@fsb.hr; +385 (1) 6168 168;

Faculty of Mechanical Engineering & Naval Architecture, I. Lučića 5, HR – 10 000 Zagreb, Croatia.
 H. Wolf, D. Banić and Ž. Božić

INTRODUCTION
Among the great number of various types of non-linear dynamic systems a very specific
group constitutes non-linear systems described by differential equations which contain non-
linear restoring characteristic with discontinuous or steep first derivative (for example
systems with clearances, rolling bearings, gears, clutches, impacting oscillators, etc.).
Frequency response characteristics of these systems are usually very complex and include
multi-frequency response for a single frequency excitation, jump phenomenon, multi-valued
solutions, and possibility of non-periodic solutions. Both periodic and non-periodic responses
can be determined in the time domain by using digital simulation. But procedures of that kind
can be exceptionally time consuming, particularly inside the frequency ranges of co-existence
of multiple stable solutions (where many combinations of initial conditions have to be
examined for obtaining all possible steady-state solutions), for lightly damped systems (since
a great number of excitation periods must be simulated to obtain a steady-state response), and
when the state of the system is near to bifurcation. These methods are not suitable for
obtaining wide range frequency response characteristics, unstable solutions and for
bifurcation analysis also. A very efficient method for solving strong non-linear differential
equations in the frequency domain is the harmonic balance method (HBM) [1-6]. When the
assumption of dominance of primary resonance in the response is satisfied, the HBM (single
harmonic) is very accurate and numerically very efficient method for obtaining periodic
response of non-linear systems with harmonic excitation. But it becomes very inaccurate if
the influence of higher harmonics in the response is significant. Multi frequency harmonic
balance methods (e.g., Incremental harmonic balance method or Newton-Raphson harmonic
balance method [7-14]) provide the study of effects of superharmonics and subharmonics to
response. These methods become exceptionally efficient in combination with path following
techniques [15-17], and can be successfully applied to a wide range of non-linear problems.
They are very well suited for parametric studies because a new solution can be sought by these
methods, with the previous solution used as a very good approximation. Since these methods
enable obtaining both dynamically stable and unstable solutions, determining of dynamical
stability of these solutions should be reliable and numerically efficient. Since the estimation
of dynamical stability of the steady state response by Floquet-Lyapunov theorem [18] is a
sensitive procedure [19-21], the factors which can lead to incorrect prediction of the
dynamical stability must be taken into consideration.
Responses determined in the time domain (MPES) and in the frequency domain (HBM and
IHBM) are considered in this paper as well as problems which can occur in estimation of
dynamical stability of the periodic solutions obtained in the frequency domain. For that
purpose, a simple single-degree-of-freedom system with piecewise-linear force-displacement
relationship subjected to a harmonic excitation is analysed.

MODEL OF A MECHANICAL SYSTEM WITH A CLEARANCE

Model of a simple mechanical system with clearance is shown in Figure 1. It consists of an
inertia element m, a linear viscous damping parameter c, and a non-linear elastic element
defined by a piecewise-linear function g(x) and a coefficient k. When the system is excited by
a periodic harmonic force F(t), the motion of the system can be described by the non-linear
differential equation:
d2 x dx
m 2 + c + kg ( x ) = F (t ) = Fm + Fp cos ( Ωt + ϕ F ) = f 0 + f C cos ( Ωt ) + f S sin ( Ωt ) , (1)
dt dt

118
Response and dynamical stability of oscillators with discontinuous or steep first derivative …

F(t) x
kg(x)
c&x m

Figure 1. Model of vibration system.

where f0 = Fm represents mean transmitted force, Fp = f C2 + f S2 is the amplitude of the

vibratory component at frequency Ω, while fC and fS are force component amplitudes of the
corresponding harmonic terms and ϕF is the excitation phase angle.
The piecewise linear function g(x) and its derivative are shown in Figure 2(a) and Figure 2(b),
respectivelly. Parameter b denotes one-half of the clearance space. Since the procedure of
prediction of the dynamical stability is based on derivative of a non-linear function,
expressions for non-linear function and its derivative are given:
g(x) = h*(x – b*), (2)
∂ g( x ) *
=h , (3)
∂x
where
⎧1 , b<x ⎫ ⎧ b, b<x ⎫
∗ ⎪ ⎪ ∗ ⎪ ⎪
h = ⎨0 , −b ≤ x ≤ b ⎬ , b = ⎨ 0, −b ≤ x ≤ b ⎬ . (4)
⎪1 , x < −b⎪⎭ ⎪−b, x < −b ⎪⎭
⎩ ⎩

∂g ( x)
∂x
g(x) 1

-b
O b x

-b O b x
(a) (b)

Figure 2. Graphs showing a) non-linear function g ( x) and b) its derivative ∂g(x)/∂x.

BRIEF DESCRIPTION OF THE APPLIED METHODS

THE INCREMENTAL HARMONIC BALANCE METHOD (IHBM)
By introducing a non-dimensional time θ as a new independent variable, the differential
equation (1) can be rewritten in the non-dimensional form:
η2 d 2 x 2ζ η d x M

ν 2 dθ 2
+
ν dθ
+ g( x ) = F (θ ) = ∑
n=0
( fn cos n(νθ ) + gn sin n(νθ )), (5)

119
 H. Wolf, D. Banić and Ž. Božić

with
x b k c f fS
x= , b = , ω0 = , ζ = , fC = C 2 , f S = ,
l l m 2mω 0 mlω 0 mlω 20
Ω
η= , τ = ω0 ⋅ t , Ωt = ητ = νθ .
ω0
In this way, the period of the response (with ν subharmonics taken in consideration) is always
2π, making it possible (by using the IHBM) to consider any number of superharmonics and
subharmonics included in the supposed approximate solution. Any characteristic dimension
of the system is denoted by l here.
Supposed approximate solution is given by:
N
x = ∑ ai cos iθ + bi sin iθ = T a,, (6)
i=0
where
T = [1, cosθ , cos 2θ ,...,cos Nθ , sin θ , sin 2θ ,...,sin Nθ ] ,

a = [ a0 , a1 ,..., a N , b1 , b2 ,..., bN ] .
T

The equation N = νM represents the number of all harmonics included in the supposed
solution, ν is the number of subharmonics and M is the number of superharmonics. By
applying this method, which consists of two basic steps: incrementation and Galerkin's
procedure, the non-linear differential equation (5) is transformed into the system of 2N + 1
linearized incremental algebraic equations:
j+1
K ∆a = r ,
j j
(7)
j+1 j+1
a = a + ∆a ,
j
(8)
with Fourier coefficients (a0, ai, bi, i = 1, …, N) as unknowns. In equations (7) and (8), the
superscript j denotes the number of iterations. In each incremental step, only linear (i.e.,
linearized) algebraic equations have to be formed and solved. A solution is obtained from the
iteration process when the corrective vector norm ||r|| is smaller than a certain (arbitrary)
convergence criterion. The comprehensive description of the method, its application to
piecewise-linear systems and the way of determining elements of Jacobian matrix K and the
corrector r in explicit form is given by Wong et al. in [10]. Generally, accuracy of the
approximate solution obtained by using IHBM depends on the number of harmonics included
in the solution, accuracy of procedures used for determining elements of K and r, and a value
of convergence criterion. Since the IHBM described by Wong et al. ([10]) is used in this
work, accuracy of the procedure of determining elements of K and r depends only on the
precision of numerical determination of times θi in which the system changes stage stiffness
region (Fig. 3). Regarding the parameter b, the three stages in the problem are defined in
accordance with (4).

THE METHOD OF PIECING THE EXACT SOLUTIONS (MPES)

By introducing the non-dimensional time τ as an independent variable (what is convenient
when one determines response in the time domain), the differential equation (1) can be
rewritten in the non-dimensional form:
d2 x dx
+ 2ζ + g( x ) = Fm + Fp cos(ητ + ϕ F ) , (9)
dτ 2
dτ

120
Response and dynamical stability of oscillators with discontinuous or steep first derivative …

x
θ0 θ L+1
θ1 θL
b
0 θ L−1
−b
θ2 θ3 θ L−2

0 2π
θ

Figure 3. Solutions to the equations x = b and x = −b .

where
Fm Fp
Fm = , Fp = .
mlω 0
2
mlω 20
Force-displacement relationship g( x ), shown in Fig. 2(a), is piecewise-linear. Local solutions
of differential equation (9) are known explicitly inside each of the stage stiffness, and can be
repeatedly matched at x = b and x = −b , to obtain a global solution of (9). Piecing together
of these local solutions is not directly possible, because the times of flight in each stage
stiffness region cannot be found in a closed form. But, the matching of local solutions can be
numerically done very easily. Only approximation done by the applying of this procedure is
in the numerical determination of times in which the system changes stage stiffness region
( x = b , x = −b ). Effective amplitudes x p of the steady state time domain responses are
calculated by using:
T /2 T /2
1 1
x (τ ) dτ , x p = 2 ( x (τ ) − x m ) dτ ,
2

T − T∫/ 2 ∫
xm = (10)
T −T /2

where T denotes the period of the response. The solutions of equation (9) inside each of the
stage stiffness region are:
1. for −b ≤ x ≤ b
C1 −2ζ (τ −τ 0 ) Fm Q
x=− e + (τ − τ 0 ) + 1a sin(ητ + ϕ F − ϕ R ) + C2 , (11a)
2ζ 2ζ η
dx F
= C1e−2ζ (τ − τ 0 ) + m + Q1a cos(ητ + ϕ F − ϕ R ) , (11b)
dτ 2ζ
where
Fp ⎛η⎞
Q1a = , ϕ R = tan −1 ⎜ ⎟,
η + 4ζ
2 2 ⎝ 2ζ ⎠
⎛d x⎞ F
C1 = ⎜ ⎟ − m − Q1a cos(ητ 0 + ϕ F − ϕ R ) ,
⎝ d τ ⎠ 0 2ζ

C1 Q1a
C2 = x 0 + − sin(ητ 0 + ϕ F − ϕ R ) .
2ζ η

121
 H. Wolf, D. Banić and Ž. Božić

in which x 0 and (d x /dτ)0 are non-dimensional displacement and velocity at the initial time
τ = τ0, i.e. at time in which the motion is determined by the given equation.
2. for x ≥ b
[
x = e−ζ (τ − τ 0 ) A2 cos 1 − ζ 2 (τ − τ 0 ) + B2 sin 1 − ζ 2 (τ − τ 0 ) + Fm + ]
+(1 − α )b + Q2 cos(ητ + ϕ F − ϕ R ) , (12a)

⎜ ⎟ =e
⎝ dτ ⎠
(
⎛ d x ⎞ −ζ (τ −τ 0 ) ⎡
⎢
⎣ ) ( )
B2 1 − ζ 2 − ζA2 cos 1 − ζ 2 (τ − τ 0 ) − A2 1 − ζ 2 + ζB2 sin 1 − ζ 2 (τ − τ 0 )⎤⎥ −
⎦
−η Q2 sin(ητ + ϕ F − ϕ R ) , (12b)
where
Fp ⎛ 2ζ η ⎞
Q2 = , ϕ R = tan −1 ⎜ ⎟,
⎝ 1 − η2 ⎠
(1 − η ) + (2ζη)
2 2 2

A2 = x 0 − Fm − (1 − α )b − Q2 cos(ητ 0 + ϕ F − ϕ R ) ,
1 ⎡⎛ d x ⎞ ⎤
B2 = ⎢⎜ ⎟ + ζ A2 + ηQ2 sin(ητ 0 + ϕ F − ϕ R )⎥ .
1 − ζ 2 ⎢⎣⎝ d τ ⎠ 0 ⎥⎦
3. for x ≤ −b
[
x = e−ζ (τ −τ 0 ) A3 cos 1 − ζ 2 (τ − τ 0 ) + B3 sin 1 − ζ 2 (τ − τ 0 ) + Fm − ]
−(1 − α )b + Q3 cos(ητ + ϕ F − ϕ R ) , (13a)

⎜ ⎟ =e
⎝ dτ ⎠
(
⎛ d x ⎞ −ζ (τ −τ 0 ) ⎡
⎢
⎣ ) ( )
B3 1 − ζ 2 − ζA3 cos 1 − ζ 2 (τ − τ 0 ) − A3 1 − ζ 2 + ζB3 sin 1 − ζ 2 (τ − τ 0 )⎤⎥ −
⎦
−η Q3 sin(ητ + ϕ F − ϕ R ) , (13b)
where
Fp ⎛ 2ζ η ⎞
Q3 j = Q2 j = , ϕ R = tan −1 ⎜ ⎟,
⎝ 1 − η2 ⎠
(1 − η ) + (2ζη)
2 2 2

A3 = x 0 − Fm + (1 − α )b − Q3 cos(ητ 0 + ϕ F − ϕ R ) ,
1 ⎡⎛ d x ⎞ ⎤
B3 = ⎢⎜ ⎟ + ζ A3 + ηQ3 sin(ητ 0 + ϕ F − ϕ R )⎥ .
1 − ζ 2 ⎢⎣⎝ d τ ⎠ 0 ⎥⎦

THE STABILITY OF THE STEADY STATE SOLUTION

When the periodic solution is obtained, the stability of the given solution can be determined
by examining the perturbed solution x *:
x* = x + ∆ x* , (14)
*
where ∆ x is a small perturbation of a periodic solution x . By substitution of equation (14)
into equation (5), and after expanding the non-linear function g( x ) in Taylor's series about
the periodic solution with neglecting non-linear incremental terms, one obtains a linear
homogeneous differential equation with time changing periodic coefficients ∂g( x )/∂ x :

122
Response and dynamical stability of oscillators with discontinuous or steep first derivative …

η2 d ∆ x 2ζ η d∆ x ∂ g( x ) *
2 * *

+ + ∆x = 0 . (15)
ν dθ
2 2
ν dθ ∂x
When the steady state solution x (θ) is determined, the values of ∂g( x )/∂ x are known inside
the period of the response. A very efficient and very often used method for determining the
stability of the periodic solution is based on the Floquet-Lyapunov theorem [18, 22]. For that
purpose equation (15) can be rewritten in the state variable form as:
d X∗
= A (θ ) X * , (16)
dθ
where
⎧ ⎫ ⎧ d∆x * ⎫ ⎡ 0 1 ⎤
⎪
⎪ ∆ * ⎪
⎪ d ∗ ⎪
⎪ ⎪
⎪ ⎢ ⎥
= ⎨ 2dθ * ⎬ , A (θ ) = ⎢ ν 2 ⎛ ∂ g ( x ) ⎞
x X
X* = ⎨ * ⎬
, 2ν ζ ⎥ . (17)
⎪ d∆x ⎪ dθ ⎪ d ∆x ⎪ − 2⎜ ⎟ −
⎢ η ⎝ ∂x ⎠ η ⎥⎦
⎩⎪ dθ ⎭⎪ ⎩⎪ dθ 2 ⎭⎪ ⎣
Since the matrix A (θ ) is a periodic function of θ with period 2π, the stability criteria are
related to the eigenvalues of the monodromy matrix, which is defined as the state transition
matrix at the end of one period. According to Floquet-Lyapunov theorem, the solution is
stable if all the moduli of the eigenvalues of the monodromy matrix are less than unity.
Otherwise the solution is unstable. Bifurcation occurs when one of the moduli of the
eigenvalues of the monodromy matrix reaches unity. Generally, it is not possible to derive an
analytic expression for the transition matrix. But, if the non-linear force-displacement
relationship is piecewise-linear, its derivative ∂g(x)/∂x = h* is, according to (4), constant
inside each of the intervals [θi, θi+1]. Figure 3 shows a period of the response where θ0 = 0
and θL + 1 = 2π. There are L times denoted as θ1, θ2, ... , θL, in which the system undergoes a
stiffness change. Consequently, A(θi, θi+1) is also a constant matrix inside that interval.
According to [23], for the constant A(θi, θi+1) (inside the interval [θi, θi+1]), transition matrix
Φ(θi+1, θi) can be expressed as:
Φ (θi +1 ,θ i ) = e A(θi ,θi+1 )⋅ (θi+1 −θi ) ; (18)

and for the whole interval [0, 2π] according to [10] one obtains:
L
⎡⎣Φ ( 2π, 0 ) ⎤⎦ = ∏e
i=0
A(θi ,θi +1 )⋅ (θi+1 −θ i )
. (19)

Beside the precision of numerical determination of times θi in which the system changes
stage stiffness region ( x = b , x = −b ), the only approximation occurring in this procedure is
the accuracy of computation of the matrix exponential exp[A(θi, θi+1)⋅(θi+1 – θi)] and the
L
product of matrix exponentials ∏e
i=0
A(θi ,θi +1 )⋅ (θi+1 −θi )
.

NUMERICAL EXAMPLES
Figures 4 and 5 show effective amplitude-frequency plots x p = x p (η ) obtained by MPES
(both periodic and non-periodic solutions) for the parameter values: b = 1, ζ = 0,03,

123
 H. Wolf, D. Banić and Ž. Božić

⎛ ⎞
Figure 4. Effective amplitude-frequency plot x p = x p (η ) obtained by MPES: x0 = 0 , ⎜ d x ⎟ = 0 .
⎝ dτ ⎠0

Figure 5. Effective amplitude-frequency plot x p = x p (η ) obtained by MPES for x0 = 0 ,

⎛dx ⎞ ⎛dx ⎞ ⎛dx ⎞ ⎛dx ⎞
⎜ ⎟ = 0 ; x0 = 1 , ⎜ ⎟ = 1 ; x0 = 1 , ⎜ ⎟ = −1 and x0 = −1 , ⎜ ⎟ = −1 .
⎝ dτ ⎠0 ⎝ dτ ⎠0 ⎝ d τ ⎠0 ⎝ d τ ⎠0

f0 = Fm = 0.25, fC = Fp = 0.25, f S = 0 (ϕ F = 0 ) . Figure 4 shows 1990 effective amplitudes

x p of the time domain responses obtained at 1990 non-dimensional frequencies η for initial

124
Response and dynamical stability of oscillators with discontinuous or steep first derivative …

conditions: x0 = 0 and (d x /dτ)0 = 0. Figure 5 shows 7960 effective amplitudes x p obtained

at the same 1990 non-dimensional frequencies η, for four different initial conditions: x0 = 0 ,
⎛dx ⎞ ⎛dx ⎞ ⎛dx ⎞ ⎛dx ⎞
⎜ ⎟ = 0 ; x0 = 1 , ⎜ ⎟ = 1 ; x0 = 1 , ⎜ ⎟ = −1 and x0 = −1 , ⎜ ⎟ = −1 . Figure 6 shows
⎝ dτ ⎠0 ⎝ d τ ⎠0 ⎝ dτ ⎠0 ⎝ dτ ⎠0
comparison of results obtained by MPES and those obtained by IHBM in the case when
supposed approximate solution includes only a constant term and the first harmonic (single
harmonic balance method). As one can see, a good agreement of the results obtained by these
two methods is achieved, but only when the assumption of dominance of primary resonance
in the response is satisfied.

Figure 6. Comparison of the results obtained by MPES (dots) and by single harmonic balance
method (line).

In Figure 7 the numerical results obtained by IHBM are compared with those obtained by
MPES. Figure 7 shows excellent agreement between the results obtained by these methods.
Non-periodic responses obtained by MPES are not found by the incremental harmonic
balance method, because this method is limited only to consideration of periodic vibrations.
Also, frequency response characteristics obtained by MPES are incomplete, because the
results of MPES depend on given initial conditions, making it difficult to find all possible
solutions.
As it is shown in [21], the accuracy of determining the eigenvalues of the monodromy matrix
depend significantly on the number of harmonics included in the supposed approximate
solution. Consequently, neglecting of very small harmonic terms of the actual time domain
response can cause a very large error in evaluation of the eigenvalues of the monodromy
matrix and can lead to incorrect prediction of the dynamical stability of the solution (Fig. 8).
Figure 8(a) shows the relative differences of the effective amplitudes x p (the a branch of the
amplitude-frequency plot from Figure 7) obtained with N = 7, 12, 15, 18 and 30 harmonics
with respect to the effective amplitudes x p obtained with N = 100 harmonics

125
 H. Wolf, D. Banić and Ž. Božić

(x ) − (x )
(x ) =
p N p 100
, N = 7, 12, 15, 18, 30 . (20)
p N diff
(x )p 100

Figure 7. Comparison of the results obtained by MPES (dots) and by IHBM (line).

Figure 8(b) shows a corresponding plot of maximum modulus of the eigenvalues of the
monodromy matrix λmax = λmax (η , N ) . A very specific situation occurs at η = 0,176. In this
case the value of λmax is more precisely determined for N = 7 than for N = 12, N = 15 and
N = 18, what can lead to incorrect estimation of dynamical stability of the response and can
make bifurcation analysis difficult. The spectrum of the corresponding time domain response
is shown in Fig. 9. In this example amplitudes of the harmonics for N > 8 are exceptionally
small in comparison to other terms of the spectrum and in-significantly influence effective
amplitude x p (amplitudes of the higher harmonics (N > 8) are less than 0.7 % of the
amplitude of the dominant harmonic). Absence of convergence can be explained by essential
difference between IHBM (the method used for obtaining approximate steady state solutions)
and the procedure of evaluating the monodromy matrix (used for estimation of dynamical
stability of the steady state solution by Floquet-Lyapunov theorem). Since the IHBM is based
on the Galerkin’s procedure, Fourier’s coefficients of the supposed approximate solution
( a0 , ai , bi , i = 1,..., N ) are determined in the way that differential equation (5) is satisfied in
average, but not in every point of the response x = x (θ). On the other hand, stability
estimation based on evaluation of the eigenvalues of the monodromy matrix depends only on
the position of points θ1 , θ 2 , …, θ L in which the system undergoes a stiffness change (Fig. 3)
and estimation of stability is influenced only by the differences between the approximate and
exact positions of the points θ1 , θ 2 , …, θ L . Increasing of number of harmonics N decreases
average difference between the approximate and the exact solution and in this way increases
the probability of more accurate determination of points θ1 , θ 2 , …, θ L .

126
Response and dynamical stability of oscillators with discontinuous or steep first derivative …

0.005
a)
N=30 N=18

N=15
-0.005
Difference of the responses

-0.01

-0.015 N=12

N=7
-0.02

-0.025

-0.03
0.174 0.175 0.176 0.177 0.178 0.179 0.18
Non-dimensional frequency

1.15
b)
N=100
1.1
N=30

1.05
unstable
Maximum modulus of the eigenvalues

0.95 stable

0.9

N=12
0.85

N=15
0.8
N=18
0.75
N=7
0.7

0.65
0.174 0.175 0.176 0.177 0.178 0.179 0.18

Non-dimensional frequency

Figure 8. a) Differences of effective amplitudes ( x p ) (the a branch of the amplitude-frequency

N diff

plot from Figure 7) obtained with N = 7, 12, 15, 18, and 30, and b) corresponding plot of maximum
modulus of the eigenvalues of the monodromy matrix λ m ax = λ m ax (η , N ) .

127
 H. Wolf, D. Banić and Ž. Božić

Figure 9. The spectrum of the time domain response for η = 0,179 (branch a in Fig 7).

CONCLUSIONS
Response and dynamical stability of oscillators with discontinuous or steep first derivative of
restoring characteristic is considered in this paper. For that purpose, a simple single-degree-
of-freedom system with piecewise-linear force-displacement relationship subjected to a
harmonic force excitation is analysed by the method of piecing the exact solutions (MPES) in
the time domain and by the incremental harmonic balance method (IHBM) in the frequency

128
Response and dynamical stability of oscillators with discontinuous or steep first derivative …

domain. The stability of the periodic solutions obtained in the frequency domain by IHBM is
estimated by the Floquet-Lyapunov theorem.
The considerable advantage of using this piecewise-linear model is in the possibility of
expressing monodromy matrix exactly as a product of matrix exponentials, what is not
possible for a general non-linear function. In this way, the inaccuracy of evaluating
monodromy matrix can be caused only by insufficient precision of numerical determination
of the times in which the system changes stage stiffness region, and by numerical procedures
of evaluation matrix exponential and product of matrix exponentials. On the other hand, local
solutions of differential equation (1) are known explicitly inside each of the stage stiffness,
and can be repeatedly matched at x = b and x = −b , to obtain a global solution of (1) in the
time domain. Piecing together of these local solutions is not directly possible, because the
times of flight in each stage stiffness region cannot be found in a closed form. But, the
matching of local solutions can be numerically done very easily. Only approximation done by
applying this procedure is in the precision of numerical determination of times in which the
system changes stage stiffness region (x = b, x = –b).
Obtained frequency response characteristic is very complex and includes multi-frequency
response for a single frequency excitation, jump phenomenon, multi-valued and non-periodic
solutions. Determining of frequency response characteristic in the time domain by MPES is
exceptionally time consuming, particularly inside the frequency ranges of co-existence of
multiple stable solutions, where many combinations of initial conditions have to be examined
for obtaining complete frequency response characteristic. In the frequency domain, IHBM is
very efficient and very well suited for obtaining wide range frequency response
characteristics, parametric studies and bifurcation analysis. On the other hand, neglecting of
very small harmonic terms (which in-significantly influence the r.m.s. values of the response
and are very small in comparison to other terms of the spectrum) can cause very large error in
evaluation of the eigenvalues of the monodromy matrix, and so they can lead to incorrect
prediction of the dynamical stability of the solution. Moreover, frequency ranges inside
which the procedure of evaluation of eigenvalues of the monodromy matrix does not
converge with increasing the number of harmonics included in the supposed approximate
solution are detected.

ACKNOWLEDGMENTS
The data shown are result of a scientific project (Numerical and experimental research of
non-linear mechanical systems) that was conducted with the support of the Ministry of
science, education and sports, Republic of Croatia.

REFERENCES
[1] Comparin, R.J. and Singh, R.: Non-linear frequency response characteristics of an
impact pair.
Journal of Sound and Vibration 134, 59-90, 1989,
[2] Comparin, R.J. and Singh, R.: Frequency response characteristics of a multi-degree-of-
freedom system with clearances.
Journal of Sound and Vibration 142, 101-124, 1990,
[3] Kahraman, A. and Singh, R.: Non-linear dynamics of a spur gear pair.
Journal of Sound and Vibration 142, 49-75, 1990,
[4] Kahraman, A. and Singh, R.: Non-linear dynamics of geared rotor bearing system with
multiple clearances.
Journal of Sound and Vibration 144, 469-506, 1991,

129
 H. Wolf, D. Banić and Ž. Božić

[5] Chatterjee, S.; Mallik, A.K. and Ghosh, A.: On impact dampers for non-linear vibrating
systems.
Journal of Sound and Vibration 187, 403-420, 1995,
[6] Chatterjee, S.; Mallik, A.K. and Ghosh, A.: Periodic response of piecewise non-linear
oscillators under harmonic excitation.
Journal of Sound and Vibration 191, 129-144, 1996,
[7] Lau, S.L.; Cheung, Y.K. and Wu, S.W.: Incremental harmonic balance method with
multiple time scales for aperiodic vibration of nonlinear systems.
Journal of Applied Mechanics 50, 871-876, 1983,
[8] Pierre, C.; Ferri, A.A. and Dowell, E.H.: Multy-harmonic analysis of dry friction
damped systems using an incremental harmonic balance method.
Journal of Applied Mechanics 52, 958-964, 1985,
[9] Choi, Y.S. and Noah, S.T.: Forced Periodic Vibration of Unsymmetric Piecewise Linear
Systems.
Journal of Sound and Vibration 121, 117-126, 1988,
[10] Wong, C.W.; Zhang, W.S. and Lau, S.L.: Periodic forced vibration of unsymmetrical
piecewise-linear systems by incremental harmonic balance method.
Journal of Sound and Vibration 149, 91-105, 1991,
[11] Lau, S.L. and Yuen, S.W.: Solution diagram of non-linear dynamic systems by the IHB method.
Journal of Sound and Vibration 167, 303-316, 1993,
[12] Leung, A.Y. and Chui, S.K.: Non-linear vibration of coupled duffing oscillators by an
improved incremental harmonic balance method.
Journal of Sound and Vibration 181, 619-633, 1995,
[13] Kahraman, A. and Blankenship, G.W.: Interactions between commensurate parametric
and forcing excitations in a system with clearance.
Journal of Sound and Vibration 194, 317-336, 1996,
[14] Lerusse, A.; Cardona, A. and Gerardin, M.: Multy harmonic balance method applied to
nonlinear structures.
2nd European Nonlinear Oscillations Conference, September 9-13 1996. Prague Czech, 1996,
[15] Narayanan, S. and Sekar, P.: A frequency domain based numeric-analytical method for
non-linear dynamical systems.
Journal of Sound and Vibration 211, 409-424, 1998,
[16] Cardona, A.; Lerusse, A. and Gerardin, M.: Fast Fourier nonlinear vibration analysis.
Computational Mechanics 22, 128-142, 1998,
[17] Raghothama, A. and Narayanan, S.: Bifurcation and chaos in geared rotor bearing
system by incremental harmonic balance method.
Journal of Sound and Vibration 226, 469-492, 1999,
[18] Minorsky, N. Nonlinear oscillations.
D. Van Nostrand Company, Princeton, New Jersey, Toronto, London, New York, 1962,
[19] Wolf, H. and Stegić, M.: The influence of neglecting small harmonic terms on
estimation of dynamical stability of the response of non-linear oscillators.
Computational Mechanics 24, 230-237, 1999,
[20] Wolf, H.; Kodvanj, J. and Bjelovučić-Kopilović, S.: Effect of smoothing piecewise-linear
oscillators on their stability prediction.
Journal of Sound and Vibration 270, 917-932, 2004,
[21] Wolf, H.; Terze, Z. and Sušić, A.: Dynamical stability of the response of oscillators with
discontinuous or steep first derivative of restoring characteristic.
European Journal of Mechanics A/Solids 23, 1041-1050, 2004,

130
Response and dynamical stability of oscillators with discontinuous or steep first derivative …

[22] Nayfeh, A.H. and Balachandram, B.: Applied Nonlinear Dynamics: Analytical,
Computational and Experimental Methods.
Wiley, New York, 1995,
[23] D’Souza, A.F. and Garg, W.H.: Advanced Dynamics Modeling and Analysis.
Prentice-Hall, Englewood Cliffs, New Yersey, 1974.

ODZIV I DINAMIČKA STABILNOST VIBRACIJSKIH SUSTAVA

S PREKINUTOM ILI STRMOM DERIVACIJOM
POVRATNE KARAKTERISTIKE
H. Wolf1, D. Banić2 i Ž. Božić1
1
Fakultet strojarstva i brodogradnje Sveučilišta u Zagrebu
Zagreb, Hrvatska
2
Grafički fakultet Sveučilišta u Zagrebu
Zagreb, Hrvatska
SAŽETAK
U radu su razmatrani odziv i dinamička stabilnost vibracijskih sustava koji imaju prekinutu ili strmu prvu
derivaciju povratne karakteristike. U tu svrhu je analiziran jednostavni vibracijski sustav s jednim stupnjem
slobode gibanja i karakteristikom krutosti koja se sastoji od linearnih segmenata koji je uzbuđen s harmonijskom
silom. Odziv sustava je u vremenskoj domeni dobiven s metodom povezivanja egzaktnih rješenja po
segmentima (MPES), a u frekvencijskoj domeni s inkrementalnom metodom harmonijske ravnoteže (IHBM).
Procjena stabilnosti periodičnih rješenja dobivenih u frekvencijskoj domeni korištenjem IHBM izvršena je
primjenom Floquet-Lyapunovog teorema. Dobiveni graf funkcije povećanja je vrlo složen i sadrži
višefrekvencijske odzive uzrokovane jednofrekvencijskom uzbudom, tzv. skokove amplitude, te višestruka i
neperiodična rješenja. Određivanje grafa funkcije povećanja u vremenskoj domeni s MPES izuzetno je
dugotrajno a to je najizraženije u područjima frekvencija u kojima postoje višestruka stabilna rješenja. IHBM, s
kojom se odziv sustava određuje u frekvencijskoj domeni vrlo je efikasna i dobro prilagođena metoda za
određivanje cjelovitog grafa funkcije povećanja, kao i za parametarsku i bifurkacijsku analizu. S druge strane,
zanemarivanje vrlo malih harmonika (koji neznatno utječu na srednju vrijednost i efektivnu amplitudu odziva i
koji su vrlo mali u odnosu na ostale harmonike u spektru) može uzrokovati vrlo velike pogreške u određivanju
vlastitih vrijednosti prijenosne matrice i tako dovesti do pogrešne procjene dinamičke stabilnosti rješenja.
Štoviše, uočena su frekvencijska područja unutar kojih postupak određivanja vlastitih vrijednosti prijenosne
matrice ne konvergira s povećanjem broja harmonika uključenih u pretpostavljeno približno rješenje.
KLJUČNE RIJEČI
dinamička stabilnost, graf funkcije povećanja, nelinearne vibracije, sustav linearan po segmentima

131