International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics

journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / n l m

On the resonance response of an asymmetric Duffing oscillator

Ivana Kovacic a,∗ , Michael J. Brennan b , Benjamin Lineton b
a
Department of Mechanics, Faculty of Technical Sciences, University of Novi Sad, 21000 Novi Sad, Serbia
b
Institute of Sound and Vibration Research, University of Southampton, Southampton SO17 1BJ, UK

A R T I C L E I N F O A B S T R A C T

Article history: The primary resonance response of an asymmetric Duffing oscillator with no linear stiffness term and
Received 12 November 2007 with hardening characteristic is investigated in this paper. An approximate solution corresponding to the
Received in revised form 26 March 2008 steady-state response is sought by applying the harmonic balance method. Its stability is also studied. It is
Accepted 17 May 2008
found that different shapes of frequency--response curves can exist. Multiple-valued solutions, indicating
the occurrence of jump phenomena, are observed analytically and confirmed numerically. The influence
Keywords: of the system parameters on the primary resonance response is also examined.
Primary resonance © 2008 Elsevier Ltd. All rights reserved.
Frequency--response curve (FRC)
Multiple jumps
Hysteretic behaviour

1. Introduction He also observed that when the damping ratio is increased, the num-
ber of the steady-state solutions decreases, so that for high damping
Non-linear systems can exhibit a wide range of phenomena that the system has just one steady-state solution. Rahman and Burton
are not found in linear systems, such as different kinds of resonance [7] showed numerically what is thought to be the first example of a
(main, subharmonic, superharmonic and ultrasubharmonic) [1,2]; "double jump” phenomenon in a single degree of freedom, harmon-
hysteresis [1,3]; co-existence of attractors, symmetry breaking, cas- ically excited Duffing oscillator with softening non-linearity, which
cade of period-doubling bifurcations and chaos [4,5]. corresponds to Eq. (1) with  < 0 and f0 = 0.
The aim of this paper is to investigate some of these phenomena The Helmholtz--Duffing equation
for an asymmetric Duffing oscillator, such as the main (primary)
resonance, hysteresis and jump phenomena. Its motion is governed z̈ + 2ż + 20 z + 2 z2 + 3 z3 = f1 cos t, (2)
by the non-dimensional equation
can also be recast in the form of Eq. (1) by using appropriate sub-
stitutions (see [8] or [9], for example). After such a transformation
ÿ + 2ẏ + y3 = f0 + f1 cos t, (1)
Ravindra and Mallik [8] applied the harmonic-balance method and
concluded that for a system with a hardening characteristic, the
where y is the displacement,  is the damping ratio,  is a parameter
jump in the transmissibility curve extends towards high frequen-
related to the non-linearity, f0 is a constant force, and f1 and  are
cies. The same authors [10] considered the Duffing equation of
the magnitude and frequency of harmonic excitation, respectively.
the hardening type with combined Coulomb and viscous damping,
Overdots denote derivatives with respect to time t. All variables are
and observed an anomalous jump in the system subject to base
non-dimensional.
excitation. The Helmholtz--Duffing equation was also analysed as
The model described by Eq. (1) belongs to the general class of
the equation of motion for a diaphragm air spring [11]. The jump
the Duffing oscillator, the origin of which was in electronics [6].
phenomena and hysteresis were studied by the bifurcation set of
Being associated with an iron core, a resonance jump phenomenon
a cusp catastrophe, which gives the relationship between the fre-
in electronics is called ferroresonance. Hayashi [6] considered
quency and amplitude of excitation yielding a discontinuous change
Eq. (1) for  =  = 1. He found that for a fixed value of f0 there can
in the amplitude of the motion. Lee et al. [12] considered microcan-
be five steady-state solutions for certain values of  and f1 , which
tilevers in tapping mode atomic microscopy when the tip-sample
implies the occurrence of multiple ferroresonance in the system.
distance was relatively large and showed numerically and experi-
mentally that four types of saddle-node bifurcations can occur in
* Corresponding author. Tel.: +381 214852241; fax: +381 21458133.
this system. A pair of stable periodic orbits and an unstable orbit
E-mail addresses: ivanakov@uns.ns.ac.yu (I. Kovacic), mjb@isvr.soton.ac.uk merged or disappeared, which led to multiple jumps and hysteretic
(M.J. Brennan), bl@isvr.soton.ac.uk (B. Lineton). behaviour. Yagasaki also studied the dynamics of two systems

0020-7462/$ - see front matter © 2008 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ijnonlinmec.2008.05.008
 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867 859

governed by the Helmholtz--Duffing equation: a microcantilever in

tapping mode atomic force microscopy when the tip-sample distance
was relatively small [13] and a pendulum subjected to linear feed-
back control [14]. He obtained the values of the magnitude of excita-
tion at which subcritical and supercritical saddle-node bifurcations
occur. Considering the case of the primary resonance, he analysed
bifurcation behaviour by using a second-order averaging method,
developed by the author in Ref. [15]. The results of this paper are
also included in Ref. [16], where this method was further developed
for analysing higher-order ultrasubharmonic, subharmonic and su-
perharmonic behaviour. In Ref. [15], Yagasaki also confirmed a fact
derived earlier by the multiple scale method [1]---for a specific or-
dering of the system parameters due to quadratic non-linearity, the
system modelled by Eq. (2) can behave as hardening or softening
depending on the sign of the expression 3 − 1022 /(920 ). Degen-
erate resonance behaviour in the Helmholtz--Duffing oscillator was
studied in Ref. [17] by using an extended version of the subharmonic
Melnikov method [18], leading to the conclusion that degenerate
resonance generally yields cusp bifurcations.
This study extends Hayashi's findings [6], by taking into account Fig. 1. The maximum number of the steady-state amplitudes: one, three or five, as
the influence of frequency on the response of the system. More- a function of the non-dimensional constant force f0 and magnitude of the harmonic
over, weakly and strongly non-linear oscillators are considered; two force f1 for the damping ratio  =0.025 and the coefficient of non-linearity  =0.0783.
values of the parameter of non-linearity are used:  = 0.0783 and
 = 3.7033, which are related to two configurations of a non-linear
Rule of Signs [20], the number of positive roots of the real algebraic
isolator [19]. In Section 2, typical frequency--response curves (FRCs)
Eq., (5), is either equal to the number of sign changes in the sequence
obtained analytically are distinguished. Their stability is analysed in
of the coefficients of the polynomial, where vanishing terms are dis-
Section 3. Numerical confirmation for the response of the system is
regarded, or it is less than that number by a positive even integer.
provided in Section 4. Section 5 is a study into the way in which
This implies that the system can have a maximum number of one,
the stiffness of the system changes during a cycle of vibration. In a
three, or five steady-state values. The way in which the maximum
linear system this is constant, but in the system investigated in this
number of the steady-state values of A0 and A1 depends on f0 and
paper it changes dramatically depending on the frequency of exci-
f1 is shown in Fig. 1 for the case when  = 0.025 and  = 0.0783.
tation. In Section 6 the influence of different system parameters on
It can be seen that for small values of the magnitude of harmonic
its response is examined.
excitation f1 , there are no multivalued amplitudes, and hence the
corresponding FRC is similar to that for a linear harmonically ex-
2. Frequency--response equations and curves cited system. To illustrate this case, the FRCs for the bias term and
the fundamental harmonic are plotted as red solid curves in Fig. 2
In the first step of the analysis of the dynamic behaviour of the for  = 0.025,  = 0.0783, f0 = 0.25 and f1 = 0.01. Also in this figure
system whose equation of motion is given by Eq. (1), an approximate are the results of the analysis from Section 2.1 and the results of
solution corresponding to the steady-state response in the region of numerical simulations described in Section 4. However, for the ma-
the primary resonance is sought. The harmonic balance method is jority of combinations shown in Fig. 1, there are three steady-state
applied and the approximate solution is assumed to be values. This multivaluedness, as well as that relating to the existence
of five steady-state amplitudes, implies the occurrence of a multiple
y(t) = A0 + A1 cos(t + ). (3) jump phenomenon. The FRCs with the maximum number of three
steady-state values corresponding to  = 0.025,  = 0.0783, f0 = 0.2
Substituting Eq. (3) into Eq. (1) and equating constant terms and the
and f1 = 0.1 are plotted in Fig. 3. Those for which five steady-state
coefficients of the terms containing cos t and sin t separately to
amplitudes occur, corresponding to  = 0.025,  = 0.0783, f0 = 0.4
zero, the system of coupled non-linear algebraic equations, in terms
and f1 = 0.1 are shown in Fig. 4. Again the results discussed in
of a bias term A0 , the amplitude of the harmonic term A1 and phase
Sections 2.1 and 4 are also plotted.
 is found to be

A30 + 32 A0 A21 = f0 , −A1 2 + 3A20 A1 + 34 A31 = f1 cos , 2.1. Peak amplitudes

− 2A1  = f1 sin . (4a-c) In order to find the values of the peak amplitudes of the FRCs
and the values of the frequency at which they occur, the case when
Combining Eqs. (4a-c) gives the implicit equation for the amplitude damping is negligible, i.e.  = 0, is first considered (although a trough
of the bias term A0 : occurs in the FRC of the bias term when a peak occurs in the FRC of
the harmonic term, the word peak is used here in both cases for the
253 A90 − 202 2 A70 − 152 f0 A60 sake of brevity). According to Eqs. (4b-c), the phase is either I = 0
or II = . The corresponding branches of the FRC for the bias term
+ 42 (2 + 42 )A50 + 16f0 2 A40
are given by
+ 3(2f12 −3f02 )A30 −4f0 2 (2 +42 )A20 +4f02 2 A0 −f03 = 0. (5) 
5 2 f 3A0
2I = A0 + 0 − f1 ,
For given values of , , f0 and f1 , this can be solved to give A0 . 2 2A0 2(f0 − A3 ) 0
Consequently, the amplitude of the harmonic term can be found from 
Eq. (4a). The dependence of these amplitudes on the frequency , 5 2 f 3 A 0
2II = A + 0 + f1 . (6)
i.e. the so-called FRC can then be plotted. According to Descartes's 2 0 2A0 2(f0 − A30 )
860 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867

1.5 2

1.5
1.48

A0
1
A0

0.5
1.46
0
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85
1.44
0.6 0.65 0.7 0.75 0.8 0.85 0.9

3
2.5
0.4 2

A1
1.5
0.3 1
0.5
A1

0.2 0
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85
0.1

0 Fig. 4. Frequency--response curves of the bias term A0 and harmonic term A1

for  = 0.025,  = 0.0783, f0 = 0.4 and f1 = 0.1: `--' solution of Eqs. (4a-c); `-- --'
0.6 0.65 0.7 0.75 0.8 0.85 0.9
unstable solution of Eqs. (4a-c); `*' numerical solution of Eq. (1) when increasing
the frequency; `o' numerical solution of Eq. (1) when decreasing the frequency; `--.'
Fig. 2. Frequency--response curves of the bias term A0 and harmonic term A1 for backbone curves and `' peak amplitudes.
 = 0.025,  = 0.0783, f0 = 0.25 and f1 = 0.01: `--' stable solution of Eqs. (4a-c);
`*' numerical solution of Eq. (1) when increasing the frequency; `--.' backbone curves
and `' peak amplitudes.
By using Eqs. (4a-c) the complete implicit equation for the FRC of
the bias term can be written as

 2
82 (f0 − A30 ) 5 2 f
1.5 2 + A21 −2 + A0 + 0 = f12 . (8)
3 A 0 2 2A0
1
The peak amplitude A0p lies on the curve (7) and, consequently,
A0

satisfies the previous equation for 2 = 2b0 . Thus, the following

0.5
implicit equation for the peak amplitude of the bias term is derived:

0
4f0 3 3f12 f2
0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 A60p − A0p + A20p − 0 = 0. (9)
5 202 5 2

In order to find the backbone curve and the peak amplitude of the
3 harmonic term, Eq. (4a) can be treated as a cubic equation in A0 .
2.5 Due to the form of the coefficients, it has one real root. Substituting
2 for A0 from Eq. (4a) into Eq. (4b) and combining this with Eq. (4c),
gives the implicit equation for the FRC of the harmonic term:
A1

1.5
1 ⎛ ⎛ 
0.5 2
⎜ ⎜ f 6
4 A1  + A1 ⎜
2 2 2 2 − 2
+ 3  ⎜
3 f0

+
0 + A1
0 ⎝ ⎝ 2
0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 4 2 8

 ⎞2 ⎞2

Fig. 3. Frequency--response curves of the bias term A0 and harmonic term A1 2
f 3 2 ⎟
6
A ⎟
3 f ⎟
for  = 0.025,  = 0.0783, f0 = 0.2 and f1 = 0.1: `--' solution of Eqs. (4a-c); `-- --' +
0 −
0 + 1⎟ + A ⎟ = f12 . (10)
unstable solution of Eqs. (4a-c); `*' numerical solution of Eq. (1) when increasing 2 4 2 8 ⎠ 4 1⎠
the frequency; `o' numerical solution of Eq. (1) when decreasing the frequency; `--.'
backbone curves and `' peak amplitudes.

The backbone curve of the peak amplitudes of the harmonic term

A1p , is
The backbone curve, which is the locus of the peak amplitudes of
the bias term A0p , is ⎛  
 ⎞2
2 2
⎜
3 f0
f0
6
A1
3 f0
f0
6
A 1 ⎟ 3 2
2b1 = 3⎜

⎟
5 2 f ⎝ 2 + 4 2 + 8 + 2 − 4 2 + 8 ⎠ + 4 A1 . (11)
2b0 = A + 0 . (7)
2 0 2A0
 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867 861

The peak amplitude of the harmonic term satisfies the following one derives
equation ⎡ ⎤
2 + 0 1 cos  − 1 sin 
⎢ ⎥
⎛  ⎞2 ⎢ −2 sin  −2 − 2 sin 2 2 −2 + 0 − 2 cos 2 ⎥




 ⎣ 1 ⎦
2 ⎜
3 f0 f02 A6 3 f0

f02 A61 ⎟ f12
A41p + 4A1p ⎝
2
+
4 2
+ 1 +
8 2
−
4 2
+
8
⎠ −
3
2
= 0. (12) 2 1 cos  2 −2 + 0 + 2 cos 2 2 − 2 sin 2
⎧ ⎫
⎪ B ⎪
⎨ ⎬
The peak amplitude of the harmonic term can also be calculated × sin  = 0. (18)
⎪
⎩ ⎪
⎭
from Eq. (4a) if the value of A0p is known. A blue dashed-dotted line cos 
in Figs. 2--4 shows the backbone curves for the bias and harmonic
term defined by Eq. (7) and (11). Squares plotted in these figures Non-trivial solutions exist only when the determinant of the coeffi-
denote the values of their peak amplitudes calculated on the basis cient matrix
 
of Eq. (9) and (4a).  2 + 0
 1 cos  − 1 sin  

 
() ≡  −2 1 sin  −2− 2 sin 2 2 −2 + 0 − 2 cos 2  (19)
 
 2 cos   − + 0 + 2 cos 2
2 2
2− 2 sin 2 
1
3. Stability of the approximate harmonic balance solution
vanishes. The stability conditions lead to () = 0 on the boundary
In the case when several stationary values exist, it is to be ex- between the stable and unstable regions and () > 0 in a stable
pected that not all of them will correspond to stable motion. Thus, region. Thus, it follows that
stability analysis of the approximate harmonic balance solution given
by Eq. (3) is necessary. To carry out this analysis, a small perturba- 0 22 − 2 21 2 + 2 21 ( 0 − 2 ) − 0 ( 0 − 2 )2 > 0. (20)
tion u(t) is introduced into Eq. (1) as follows:
The results from this analysis are shown in Figs. 3 and 4, where the
dashed parts of the FRCs represent the unstable regions. When three
y(t) = A0 + A1 cos(t + ) + u(t). (13) steady states occur in the system for a single frequency, two of them
are stable and one unstable (Fig. 3). If there are five steady states,
The corresponding linearized variational equation is then three of them are stable and two unstable (Fig. 4). The condition
defined by Eq. (20) is obviously equivalent to the criterion that a
stability limit is given by a vertical tangent to the resonance curves.
ü + 2u̇ + 3(A0 + A1 cos(t + ))2 u = 0. (14)
4. Numerical simulations and confirmations
Using the substitution u(t) = e−t v(t), Eq. (14) transforms to Hill's
equation [6] In this section, the results previously obtained are compared with
the results of direct numerical integration of the equation of motion,
⎛ ⎞
2 Eq. (1). It should be emphasized that the system parameters are

v̈ + ⎝ 0 + 2 n cos n(t + )⎠ v = 0, (15) chosen in such a way that the response is predominantly at the
n=1 same frequency as the harmonic excitation, and that all other Fourier
components can be neglected in comparison to the bias term and
where the first harmonic. The numerical results were obtained using the
MATLAB ode45 function and calculating the first harmonic from the
3 2 Fourier series coefficients of the steady-state response. These are
0 = 3A20 + A − 2 , also shown in Figs. 2--4 using the symbols `*' and `o'. The former
2 1
corresponds to the results obtained as frequency is increased and the
1 = 3 A 0 A 1 , 2 = 34 A21 . (16) latter as frequency is decreased. It can be seen that in the frequency
range considered, both the bias term and the first harmonic are
According to Floquet theory, a system driven parametrically, such predicted reasonably well using the harmonic balance method.
as the system modelled by Eq. (15), can exhibit resonance when- As already mentioned and confirmed by the numerical results,
√ √
ever the driving frequency is equal to 2 0 /n, where 0 is the multivaluedness of the response is responsible for the jump phe-
the normalized natural frequency of the system and n is an in- nomenon. This can occur at several frequencies as shown in Figs. 3
teger [1,6]. For the stability analyses of oscillations having the and 4. To investigate these cases more thoroughly and to emphasize
same frequency as the approximate harmonic balance solution, the occurrence of non-linear hysteresis, the corresponding FRCs for
the second unstable region is of interest, i.e. n = 2. By virtue of A0 and A1 are studied further.
Floquet theory, the solution to Eq. (15) can be assumed to be of The excitation force and the damping ratio are held fixed, while
the form the excitation frequency is slowly varied up and down, i.e. it is varied
quasi-statically. The starting point corresponds to point 1 on the FRCs
in Fig. 5. As the frequency is slowly increased, the amplitudes for A0
v(t) = e
(B + sin(t + )), (17)
and A1 follow the path labelled 1--2--3--4--5--6. Point 2 is a jump-
down point for the bias term and a jump-up point for the amplitude
where  is the characteristic Floquet exponent. On the bound- of the harmonic term; on the other hand, point 4 is a jump-up point
ary between the stable and unstable region, the real part of for the bias term and a jump-down point for the amplitude of the
the term (− ± ) is equal to zero. The stability of the approx- harmonic term. If the frequency is slowly decreased starting from
imate solution of Eq. (1) is determined by the condition that point 6, the path 6--7--8--9--10 is followed. For the bias terms (and,
the real part of the term (− ± ) should be negative. Since respectively, the amplitude of the harmonic term) point 7 is a jump-
the characteristic Floquet exponent can be either real or imagi- down (jump-up) point, while point 9 is a jump-up (jump-down)
nary, this condition is equivalent to  > 0 and 2 > 2 . Substituting point. In all case, when the bias term experiences a jump-down, the
Eq. (17) into Eq. (15) and applying the harmonic balance method, oscillatory term experiences a jump-up and vice versa.
862 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867

1.5 4
1 6
5
10 2
1 3
A0

9
7
0.5 3
2
8 4
0
0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85
1

y
0
3
4
2.5
2 3 8 -1
A1

1.5 7
9
1 -2
2
5
0.5 1 6
10
0
-3
0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85
0 900 1800 2700 3600
t
Fig. 5. Jump phenomena for the system defined in Fig. 3.

The overall time-history response of the system corresponding

2
to this hysteretic behaviour, obtained by carrying out numerical in-
tegration of the equation of motion given by Eq. (1), is illustrated
in Figs. 6a and b for slowly increasing and decreasing frequency, re-
1
spectively.
If the system shown in Fig. 4 is subject to slowly varying fre-
quency, similar jump phenomena can be observed. An illustration of
0
y

this is presented in Fig. 7. As  is increased, the responses follow the

path 1--2--3--4--5--6. If the process is reversed starting from point 6,
the route 6--7--8--9--10 is taken. In both directions two jumps oc-
-1
cur. The peak amplitude is reached when approached from a lower
frequency. It is interesting to note the difference between this re-
sponse and that one shown in Fig. 5. In the case depicted in Fig. 7,
-2
a part of the FRC corresponding to a stable steady state is unattain-
able when frequency is increased or decreased (unstable parts are
always unattainable). This is the part between points 8 and 3. These
-3
steady states can be achieved if the frequency is adjusted to some 0 900 1800 2700 3600
value in the region between 8 and 3 and the initial conditions are t
such that the attractor corresponding to point P exists. The domains
of attractors for  = 0.66 for three coexisting periodic attractors P, Fig. 6. The overall time--history diagram as the illustration of a double jump phe-
Q and R shown in Fig. 7 are presented in Fig. 8, together with the nomenon for the system defined in Fig. 3 for the frequency rate d/dt=1.3389×10−4 :
(a) obtained by increasing the frequency from (t=0)=0.35 to (t=3600)=0.85 and
corresponding phase projections and Poincaré points.
(b) obtained by decreasing the frequency from (t = 0) = 0.85 to (t = 3600) = 0.35.
The FRCs shown in Figs. 3b and 4b exhibit a softening type char-
acteristic at low amplitudes and a hardening characteristic at higher
amplitudes as frequency is increased. Carnegie and Reift [21] re-
ported similar behaviour in the undamped Helmholtz--Duffing os- jump-up frequencies and the FRC between them, and confirmed it
cillator excited by a centrifugal exciting force (the magnitude of this numerically. The response curves obtained in all the aforementioned
force is not constant but it increases with the square of frequency), studies are characterized by the bend, first to the left and then to
which was derived from the original asymmetric Duffing equation; the right. In addition, the branches intersect at a point that is not
unlike Eq. (1), however, this contained a linear term. In this system a bifurcation point (see, for example, Figure 3 in [21], Figure 10 in
the softening then hardening behaviour was found to be due to the [12] and Figure 14b in [13]). However, in the study presented here,
second, rather than the first harmonic, which occurred as a predom- in which the system motion is governed by Eq. (1) with no linear
inant component of motion for a relatively high amplitude of exci- term and its response is defined by Eqs. (4a-c), such a distinguishing
tation. Lee et al. [12] observed the initial softening and subsequent feature is not observed.
hardening in the overall non-linear response of microcantilevers in
tapping mode atomic microscopy, which represent non-smooth sys- 5. Stiffness variation
tems with simultaneous parametric and base excitation. Yagasaki
analysed that system for the case when a parametric excitation term It can be seen from the equation of motion for the system given
can be neglected [13]. Utilizing some of his previous results regard- in Eq. (1) that the only cause of non-linearity is the stiffness. Unlike
ing Melnikov's method [17,18], he determined the jump-down and a linear system, this stiffness changes during each cycle of vibration.
 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867 863

This is illustrated in Figs. 9a and b for the system response given in Each plot in Fig. 9 shows the variation of the non-dimensional stiff-
Figs. 5a and b, for increasing and decreasing frequency, respectively. ness K = 3y2 as a function of the non-dimensional displacement y.
The plots marked points 1--10 correspond to the frequencies indi-
cated in Figs. 5a and b. Also marked on each plot is the static dis-
2 placement given by `o', the bias term of the displacement, A0 , given
1 by `∇' and the peak harmonic displacement either side of this dis-
10 2 6
1.5 R
placement A0 ± A1 , given by `*'. It is worth mentioning that A0 and
5
Q A1 in Fig. 9 were obtained by the harmonic balance approximation.
A0

1
9 There are several notable features about these graphs. One is that
7
0.5 the vibration at points 1 and 6 (at low and high frequencies, respec-
8 P 4
Ω8
3 tively) are very similar, with the bias position (A0 ) being almost the
Ω3
0 same as the static equilibrium position, and the dynamic responses
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 about this position being roughly the same. As the frequency is in-
creased to point 2, the bias position decreases a little and the vibra-
tion about this position increases. However, although the stiffness
changes during a cycle, it is always greater than zero. When the
3
8 P 3 4 system experiences a jump-up in A1 (from point 2 to 3) there is a
2.5
dramatic reduction in the bias position and hence the associated in-
2 7
9 Q stantaneous stiffness, and a dramatic increase in the vibration about
A1

1.5 5
this position. It can be seen that the stiffness changes dramatically
1 6 during a cycle from a large value to zero and back to a large value.
2
0.5 1 At this frequency and at the frequency corresponding to point 4, the
Ω8 R Ω3
0 10 vibration is almost symmetric about the origin (y = 0), which indi-
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 cates that the constant force has very little effect on the dynamic
response at this frequency.
Fig. 7. Jump phenomena for the system defined in Fig. 4. When a jump-down in the A1 FRC occurs (from point 4 to point 5),
there is a dramatic increase in the bias position and a dramatic

2 P

Q
y(0)

0
R
-1

-2

-3
-3 -2 -1 0 1 2 3
y(0)

y y y
1.5 0.1
1
1
0.2
0.5 0.5
y y y
-2 -1 1.4 1.6 1.8
-0.5 1 2 3 -1 1 2 3 -0.1
-0.5
-1 -0.2
-1
-1.5

Fig. 8. (a) Domains of periodic attractors P, Q and R from Fig. 7 when  = 0.66; phase projections and Poincaré points corresponding to the attractor; (b) P for y(0) = 1,
ẏ(0) = 0; (c) Q for y(0) = 1, ẏ(0) = −3 and (d) R for y(0) = 1, ẏ(0) = −1.
864 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867

Fig. 9. Non-dimensional stiffness K = 3y2 , as a function of non-dimensional displacement y, for the asymmetric Duffing oscillator;  = 0.025,  = 0.0783, f0 = 0.2 and f1 = 0.1.
Points 1--10 correspond to the points in Fig. 5: `o' static displacement; `∇' displacement A0 due to the constant force f0 and `*' peak harmonic displacement either side of A0 .

decrease in the vibration about this position. Thus the vibration reduction in the stiffness at the bias position but a much greater
is asymmetric again, with the influence of the constant force be- change in stiffness during a cycle. Conversely, a jump-down in
ing felt by the system. In summary, it can be seen that the stiff- the A1 FRC generally results in an increase in the stiffness at the
ness variation during a cycle of vibration is heavily dependent bias position and a much smaller change in the stiffness during a
upon frequency. A jump-up in the A1 FRC generally results in a cycle.
 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867 865

6. Effects of the system parameters on the resonance response

0.8
In this section, the influence of the system parameters given by
Eqs. (4a-c) on the FRCs is examined. 0.6

A0
0.4
6.1. Effect of the magnitude of the forces
0.2
A series of graphs is presented in Figs. 10--12 that illustrate the 0
effects of the magnitude of the constant force and the harmonic force 0.3 0.4 0.5 0.6 0.7 0.8 0.9
on the system resonance response.
If the magnitude of the harmonic excitation is small and fixed,
leading to a linear-like FRC (Fig. 10), an increase in the magnitude
3
of the constant force yields a decrease in peak amplitude of A1 and
a shift in the peaks towards higher frequencies. In this case, the
coefficient in front of the quadratic term in Eq. (9) can be neglected 2

A1
so that
 the peak of the FRC of the bias term flattens to a value of
A0 = 3 f0 / . According to Eq. (7), the peak occurs at higher frequencies 1
as f0 is increased.
In the case when three possible steady-state values exist, two 0
different shapes of the FRCs are observed for one fixed value of the 0.3 0.4 0.5 0.6 0.7 0.8 0.9
magnitude of the harmonic excitation. When the constant force is
small (Fig. 11), the FRC for the harmonic term have the same shape as
that of the hardening cubic Duffing oscillator [1]. One jump appears Fig. 11. Frequency--response curves of the bias term A0 and harmonic term A1 for
 =0.025,  =0.0783 and f1 =0.1 for different values of the magnitude of the constant
when increasing frequency, and one when decreasing frequency. As
force: f0 = 0.01 (red solid line) and f0 = 0.03 (blue dashed-dotted line).
the magnitude of the constant force becomes larger, the frequency
corresponding to the jump point for decreasing frequency slightly
increases, while the frequency corresponding to the peak amplitudes
remains the same.
For higher values of the constant force, the FRCs bend several 2
times (Fig. 12), so that multiple jump-up and jump-down points oc- 1.5
cur. It is noticeable that even in this case, the frequency at which the
peak amplitude occurs does not change with the magnitude of the
A0

1
constant force. However, the values of the frequency corresponding
to other jump points do change. 0.5

0
0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9

3
2.5
2
A0

1.5
A1

1.5
1
0.5
1 0
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9

Fig. 12. Frequency--response curves of the bias term A0 and harmonic term A1 for
0.4  =0.025,  =0.0783 and f1 =0.1 for different values of the magnitude of the constant
force: f0 = 0.25 (red solid line) and f0 = 0.3 (blue dashed-dotted line).
0.3
A1

0.2 The effect of increasing the magnitude of the constant force on the
shape of the FRCs for the case when five steady states exist is shown
0.1 in Fig. 13. As it increases, the frequency at which jump-up and jump-
down points occur also becomes higher. In addition, for some values
0 of the magnitude of the harmonic force, the frequency corresponding
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 to the jump-up point when increasing the frequency is higher than
the frequency corresponding to the peak (see green dashed line in
Fig. 13). Consequently, the peak amplitude is unattainable with a
Fig. 10. Frequency--response curves of the bias term A0 and harmonic term A1 for
 = 0.025,  = 0.0783 and f1 = 0.01 for different values of the magnitude of the slowly increasing frequency, and there is only one jump during that
constant force: f0 =0.1 (red solid line); f0 =0.25 (blue dashed-dotted line) and f0 =0.5 process. The critical value of f1 yielding such behaviour, as well as
(green dashed line). the critical value of f0 that changes the FRC from that typical of
866 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867

1.5
A0

0.5

0
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85

3
2.5
2
A1

1.5
1
0.5
0
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85
Fig. 15. The maximum number of the steady-state amplitudes: one, three or five, as
a function of the non-dimensional constant force f0 and magnitude of the harmonic
Fig. 13. Frequency--response curves of the bias term A0 and harmonic term A1 for
force f1 for the damping ratio  =0.025 and the coefficient of non-linearity  =3.7033.
 =0.025,  =0.0783 and f1 =0.1 for different values of the magnitude of the constant
force: f0 = 0.35 (red solid line); f0 = 0.4 (blue dashed-dotted line) and f0 = 0.45 (green
dashed line).
a hardening system (Fig. 11) to one with multiple bends shown in
Fig. 13, will be addressed in future studies.

0.1 6.2. Effect of damping

ζcr3-1
0.08
Numerical calculations were carried out to examine the influence
0.06 of damping on the different shapes of the FRCs. Two curves corre-
ζcr

0.04 ζcr5-3 sponding to the critical values of the parameter cr are shown in
Fig. 14a as a function of f0 with f1 = 0.1. One of curves cr5−3 corre-
0.02
sponds to the damping ratios that change the response of the sys-
0.05 tem from that with five steady-state amplitudes to one with three.
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 The other curve cr3−1 gives the values of damping that changes the
f0 response with three steady states to a single-valued response. Illus-
trations of the way in which the damping ratio affects the response
of the system when the other parameters are held fixed are given
2
in Figs. 14b and c. They confirm that damping can be used to avoid
1.5 the appearance of jumps in the system.
A0

1
6.3. Effect of the coefficient of non-linearity
0.5
Fig. 15 shows the maximum number of the steady-state am-
0
plitudes for the case when the damping ratio is  = 0.025, while
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85
 = 3.7033. Comparing this graph with the one in Fig. 1, which cor-
responds to the weakly non-linear system  = 0.0783, it can be seen
that this strongly non-linear system has fewer possibilities for yield-
3 ing single-valued FRCs. It can also be seen that FRCs with five steady-
2.5 state amplitudes occur for smaller magnitudes of the constant force
2 and harmonic force. This makes the configuration of the strongly
non-linear system less favourable than the weakly non-linear one.
A1

1.5
1
0.5 7. Conclusions
0
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 The primary resonance response of a non-linear oscillatory sys-
tem that is under simultaneous constant and harmonic excitation
has been investigated. The corresponding equation of motion is the
Fig. 14. (a) The critical values of the damping ration cr5−3 and cr3−1 as the functions asymmetric Duffing equation with no linear term and with hard-
of the parameter f0 for  =0.0783 and f1 =0.1. (b) and (c) Frequency--response curves
ening non-linearity. The harmonic balance method was applied to
of the bias term A0 and harmonic term A1 for f0 = 0.4 and for different values of
the damping ratio:  = 0.025 (red solid line);  = 0.05 (blue dashed-dotted line) and determine the frequency--response equations. The stability of the
 = 0.075 (green dashed line). system was studied by applying Floquet theory. It was found that,
 I. Kovacic et al. / International Journal of Non-Linear Mechanics 43 (2008) 858 -- 867 867

as with the effect of a different ordering of the asymmetry in the [4] J.M.T. Thompson, H.B. Stewart, Nonlinear Dynamics and Chaos, Wiley,
applied external force with respect to the values of the system pa- Chichester, 1986.
[5] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and
rameters, the FRCs can have markedly different shapes. First, they Bifurcation of Vector Fields, Springer, New York, 1983.
can be single-valued, having the shape similar to the response of a [6] C. Hayashi, Nonlinear Oscillations in Physical Systems, McGraw-Hill, New York,
linear system. Second, they can be multivalued with the maximum 1964.
[7] Z. Rahman, T.D. Burton, Application of catastrophe theory to forced vibration
of three steady-state values. In this group, different shapes can be of a diaphragm air spring, J. Sound Vib. 110 (1986) 363--380.
distinguished: one in which the peak is bent in a single direction, as [8] B. Ravindra, A.K. Mallik, Performance of nonlinear vibration isolators under
it is for the response of a symmetric Duffing oscillator; and the other, harmonic excitation, J. Sound Vib. 170 (3) (1994) 325--337.
[9] W. Szemplinska-Stupnicka, J. Bajkowski, The 1/2 subharmonic resonance and
in which the peak has a double bend---first towards lower frequen-
its transition to chaotic motion in a non-linear oscillator, Int. J. Non-Linear
cies, and then towards higher frequencies. Finally, the FRCs can have Mech. 21 (1986) 401--419.
double bends with five steady-state values. The multivaluedness of [10] B. Ravindra, A.K. Mallik, Hard Duffing-type vibration isolator with combined
Coulomb and viscous damping, Int. J. Non-Linear Mech. 28 (4) (1994) 427--440.
these curves causes the occurrence of the multiple jumps in the sys-
[11] A. Murata, Y. Kume, F. Hashimoto, Application of catastrophe theory to forced
tem. It is possible for the system to exhibit one or two jumps when vibration of a diaphragm air spring, J. Sound Vib. 112 (1987) 31--44.
increasing frequency and one or two when decreasing frequency. The [12] S.I. Lee, S.W. Howell, A. Raman, R. Reifenberger, Nonlinear dynamics of
effects of different parameters on the system behaviour have been microcantilevers in tapping mode atomic force microscopy: a comparison
between theory and experiment, Phys. Rev. B 66 (2002) 115409.
analysed. The analysis has shown that a suitable choice of the damp- [13] K. Yagasaki, Bifurcations and chaos in vibrating microcantilevers of tapping
ing ratio can entirely eliminate the possibility of the appearance of mode atomic force microscopy, Int. J. Non-Linear Mech. 42 (2007) 658--672.
jumps. [14] K. Yagasaki, A simple feedback control system: bifurcations of periodic orbits
and chaos, Nonlinear Dyn. 9 (1996) 391--417.
[15] K. Yagasaki, Second-order averaging and Melnikov analyses for forced non-
Acknowledgement linear oscillators, J. Sound Vib. 190 (1996) 587--609.
[16] K. Yagasaki, Higher-order averaging and ultra-subharmonics in forced
oscillators, J. Sound Vib. 210 (1998) 529--553.
This study has been financially supported by the Ministry of Sci- [17] K. Yagasaki, Degenerate resonances in forced oscillators, Discrete Continuous
ence and Environmental Protection, Republic of Serbia (Project no. Dyn. Syst. B 3 (2003) 423--438.
144008). [18] K. Yagasaki, The Melnikov theory for subharmonics and their bifurcations in
forced oscillations, SIAM J. Appl. Math. 56 (1996) 1720--1765.
[19] I. Kovacic, M.J. Brennan, T.P. Waters, A study of a non-linear vibration isolator
References with a quasi-zero stiffness characteristic, J. Sound Vib. 315 (2008) 700--711.
[20] G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers,
[1] A. Nayfeh, D.T. Mook, Nonlinear Oscillations, Willey, New York, 1979. McGraw-Hill, New York, 1961.
[2] J. Awrejcewicz, Deterministic Vibrations of Discrete System, WNT, Warszawa, [21] W. Carnegie, Z.F. Reift, Ultraharmonic resonance of a system with an
1996 (in Polish). asymmetrical restoring force characteristic, J. Mech. Eng. Sci. 11 (1969)
[3] K. Magnus, Vibrations, Blackie and Son, London, 1965. 592--597.