# Nonlinear Dynamics Analysis of the Semiactive Suspension System with Magneto-Rheological Damper

Authors: Hailong Zhang; Enrong Wang; Fuhong Min; Ning Zhang; Chunyi Su; Subhash Rakheja
Journal: Shock and Vibration
 (2015)

---

## Abstract

This paper examines dynamical behavior of a nonlinear oscillator which models a quarter-car forced by the road profile. The magneto-rheological (MR) suspension system has been established, by employing the modified Bouc-Wen force-velocity (F-v) model of magneto-rheological damper (MRD). The possibility of chaotic motions in MR suspension is discovered by employing the method of nonlinear stability analysis. With the bifurcation diagrams and corresponding Lyapunov exponent (LE) spectrum diagrams detected through numerical calculation, we can observe the complex dynamical behaviors and oscillating mechanism of alternating periodic oscillations, quasiperiodic oscillations, and chaotic oscillations with different profiles of road excitation, as well as the dynamical evolutions to chaos through period-doubling bifurcations, saddle-node bifurcations, and reverse period-doubling bifurcations.


---

## Body

## 1. Introduction

Magneto-rheological fluid (MRF) is a suspension of micron-sized, magnetic particles in a carrier fluid. When exposed to a magnetic field, the rheology of MRF reversibly and consecutive changes from a free-flowing Newton liquid to a semisolid Bingham with controllable yield strength, which is known as magneto-rheological effect [1–3]. Magneto-rheological damper (MRD) based on MRF has significant promise for effective vibration damping in many applications, such as MRD-based semiactive suspension system, which has attracted much attention to improve the ride comfort and handling safety of the vehicle [4–6]. Nevertheless, there are still numerous challenges in controlling a MRD to get the superior performances, because it has highly nonlinear characteristics and chaotic motion due to typical hysteretic characteristics.In the past decades, many studies of MR suspension focused on the semiactive control algorithms. Many new ideas were proposed, such as “Skyhook” control or artificial intelligence based control, and had been implemented in a practical vehicle. However, chaotic behavior may exist due to the nonlinearities in MR suspension, which was not studied well yet. Li et al. [7] investigated the chaotic motion in nonlinear suspension system with hysteretic characteristics and verified the path from quasiperiodic to chaos by deriving Melnikov method. Siewe Siewe [8] applied the method of multiple scales to analyze local bifurcation in the quarter-car system with periodically excited road profile, and a variety of nonlinear behaviors were found, such as resonance and antiresonance phenomena and saddle-node bifurcation. Litak et al. [9] used the analytical Melnikov theory and predicted the lowest critical amplitude that a single degree of freedom (DoF) vehicle model may transit to a chaotic motion, under a road surface profile consisting of harmonic and noisy components. Luo and Rajendran [10] carried out the periodic motion and stability of a single DoF semiactive suspension model by developing a mapping structure, and the model of MRD was formulated with piecewise linear equations. However, all of above models are simple as single DoF model which is very different from the practical situation owing to the neglect of wheel movement, and this may explain that the vibrations appear for impossible amplitude of 0.4 m in [9]. Consequently, Borowiec and Litak [11] studied a 2-DoF quarter-car and find out the transition to chaos, by applying Melnikov theory and recurrence approach. In earlier studies, the suspension system was mostly considered as single DoF system for its simplicity. What is more, calculation models of MRD were mostly made employing polynomials or piecewise linear model in [9], while the Bouc-Wen model is considered appropriate to describe well the dynamic performance of MRD [12–14]. So far, there are no systematic nonlinear dynamics of the Bouc-Wen model based MR suspension system, owing to the complex structure of model.In this paper, from the experimental results, such a commercial MRD is modeled using a modified Bouc-Wen model proposed by the author. A 2-DoF model is established to express the MR suspension system, by employing the identification results. The stability of the system is analyzed according to the stability criterion, and all possible motions of the MR suspension system are determined by plotting the frequency response, bifurcation diagrams, and phase portraits under different road profiles. The Lyapunov exponent (LE) is calculated to detect chaotic motion. Time series with combination of power spectrum density is used as assisted means for the special system parameters.The paper is organized as follows. In Section2, the dynamic model of 2-DoF MR suspension system is formulated, employing the modified Bouc-Wen calculation model of MRD. In Section 3, the stability of the system is analyzed by calculating the eigenvalue of the Jacobian matrix of fixed point. Next, the numerical calculation is conducted to determine the dynamic behavior and to confirm the stability analysis as well. Comprehensive numerical results include frequency response, bifurcation diagrams, phase plane portraits, Poincare map, and time series; thus, the process of the transition to chaotic motion is revealed. Finally, the conclusions of the research are presented in Section 4.


## 2. Dynamic Model of Quarter-Car Semiactive Suspension System

### 2.1. Mechanical Model and Formulation

In Figure1, a classic dynamic model of a quarter MR suspension system [15] is presented, which involves main vehicle components such as the car body, suspension spring, MRD, and wheel. Amongst, only the MRD has strong force hysteresis and saturation nonlinearities. The other components have been addressed in linearization. m
s and m
u are defined as the sprung and unsprung masses. k
s, i
d, and F
d represent the suspension stiffness, direct drive current, and the yielded damping force of the MRD. k
t and c
t denote equivalent stiffness and damping coefficient of the tire, and x
i, x
s, and x
u are defined as the road excitation and vertical displacements of the sprung and unsprung masses, respectively. Following the principle of Newton’s second law of motion, the dynamic equation is formulated as(1)
m
s
x
¨
s
+
k
s
x
s
-
x
u
+
F
d
=
0
,
m
u
x
¨
u
+
k
t
x
u
-
x
i
+
c
t
x
˙
u
-
x
˙
i
-
k
s
x
s
-
x
u
-
F
d
=
0
.Figure 1
Schematic diagram of 2-DoF MR suspension system.


### 2.2. Modified Bouc-Wen Calculation Model of MRD

The accurate and practical MRD dynamic model is crucial for application. Nevertheless, there is no recognized MRD dynamic model yet. Amongst, relative effective models include Bouc-Wen model by Spencer and the phenomenon model based on Bouc-Wen [14]. Figure 2 shows the structure of typical Bouc-Wen phenomenon model, and it accurately describes the MRD inherent hysteresis nonlinear properties. However, it cannot describe nonlinear response and saturation characteristics of the magnetic field, because the linear output item in the model is just to represent the relationship between the damping force and the control current. Therefore, the author has imported the sigmoid function to improve the conventional Bouc-Wen phenomenon model [16]. The issue is effectively solved by decoupling the hysteretic characteristics with the current modulation separation:(2)
F
d
=
f
i
d
=
c
i
d
F
h
x
r
,
v
r
,
0
≤
i
d
≤
I
m
,
(3)
c
i
d
=
1
+
k
2
1
+
exp
⁡
-
a
0
i
d
+
I
0
-
k
2
1
+
exp
⁡
-
a
0
I
0
,
(4)
F
h
x
r
,
v
r
=
c
1
y
˙
+
k
1
x
r
-
x
0
,
(5)
x
r
=
x
s
-
x
u
,
y
˙
=
1
c
0
+
c
1
α
z
+
c
0
x
˙
r
+
k
0
x
r
-
y
,
z
˙
=
-
γ
x
˙
r
-
y
˙
z
z
n
-
1
-
β
x
˙
r
-
y
˙
z
n
+
A
x
˙
r
-
y
˙
,where i
d and I
m express the direct current and its maximum value for driving the MRD, respectively, and 0
≤
i
d
≤
I
m. c
(
i
d
) denotes the saturated nonlinear direct current control function proposed by the authors and c
(
i
d
)
≥
1, c
(
i
d
)
=
1, for i
d
=
0. F
h
(
v
r
) denotes yielded passive damping force with hysteresis depending on the piston relative displacement velocity (
v
r
) of the MRD for i
d
=
0. x
r is the piston travel of the MRD, y and z are inner variables without units, and k
0, k
1, k
2, a
0, I
0, α, β, γ, c
0, c
1, n, A, and x
0 are constants, respectively.Figure 2
Bouc-Wen mechanical model.As is shown in Figure3(a), a CARRERA MagneShock MRD of the vehicle suspension is further employed from [16], which permits maximum control direct current 0.5 A at 12 V [16]. On the basis of the actual measured characteristic data of MRD, the model parameters are identified as k
0
=
184.1, k
1
=
1528.1, k
2
=
10.092, a
0
=
7.526, I
0
=
0.069, α
=
20373.7, β
=
233849.1, γ
=
8816.9, c
0
=
1368.7, c
1
=
6222.7, n
=
2, A
=
20.6, and x
0
=
-
0.004. Figure 3(b) shows the comparison of the modified Bouc-Wen hysteresis F
-
v model and the tested data of employed candidate MRD, under the harmonic excitation with amplitude 12.5 mm at 1.5 Hz and different direct drive current (0–0.4 A), which exhibits ideal coordination.Figure 3
(a) A pictorial view of the candidate MRD. (b) Comparison of tested and calculated results.
(a)
(b)MRD has clear strong nonlinear property from the above results. Due to the damping force dependency on driven current and excitation frequency, theF
-
v curve varies evidently under different driven and excitation profiles [17]. As a result, nonlinear behavior caused by hysteresis cannot be ignored in engineering application.


## 2.1. Mechanical Model and Formulation

In Figure1, a classic dynamic model of a quarter MR suspension system [15] is presented, which involves main vehicle components such as the car body, suspension spring, MRD, and wheel. Amongst, only the MRD has strong force hysteresis and saturation nonlinearities. The other components have been addressed in linearization. m
s and m
u are defined as the sprung and unsprung masses. k
s, i
d, and F
d represent the suspension stiffness, direct drive current, and the yielded damping force of the MRD. k
t and c
t denote equivalent stiffness and damping coefficient of the tire, and x
i, x
s, and x
u are defined as the road excitation and vertical displacements of the sprung and unsprung masses, respectively. Following the principle of Newton’s second law of motion, the dynamic equation is formulated as(1)
m
s
x
¨
s
+
k
s
x
s
-
x
u
+
F
d
=
0
,
m
u
x
¨
u
+
k
t
x
u
-
x
i
+
c
t
x
˙
u
-
x
˙
i
-
k
s
x
s
-
x
u
-
F
d
=
0
.Figure 1
Schematic diagram of 2-DoF MR suspension system.


## 2.2. Modified Bouc-Wen Calculation Model of MRD

The accurate and practical MRD dynamic model is crucial for application. Nevertheless, there is no recognized MRD dynamic model yet. Amongst, relative effective models include Bouc-Wen model by Spencer and the phenomenon model based on Bouc-Wen [14]. Figure 2 shows the structure of typical Bouc-Wen phenomenon model, and it accurately describes the MRD inherent hysteresis nonlinear properties. However, it cannot describe nonlinear response and saturation characteristics of the magnetic field, because the linear output item in the model is just to represent the relationship between the damping force and the control current. Therefore, the author has imported the sigmoid function to improve the conventional Bouc-Wen phenomenon model [16]. The issue is effectively solved by decoupling the hysteretic characteristics with the current modulation separation:(2)
F
d
=
f
i
d
=
c
i
d
F
h
x
r
,
v
r
,
0
≤
i
d
≤
I
m
,
(3)
c
i
d
=
1
+
k
2
1
+
exp
⁡
-
a
0
i
d
+
I
0
-
k
2
1
+
exp
⁡
-
a
0
I
0
,
(4)
F
h
x
r
,
v
r
=
c
1
y
˙
+
k
1
x
r
-
x
0
,
(5)
x
r
=
x
s
-
x
u
,
y
˙
=
1
c
0
+
c
1
α
z
+
c
0
x
˙
r
+
k
0
x
r
-
y
,
z
˙
=
-
γ
x
˙
r
-
y
˙
z
z
n
-
1
-
β
x
˙
r
-
y
˙
z
n
+
A
x
˙
r
-
y
˙
,where i
d and I
m express the direct current and its maximum value for driving the MRD, respectively, and 0
≤
i
d
≤
I
m. c
(
i
d
) denotes the saturated nonlinear direct current control function proposed by the authors and c
(
i
d
)
≥
1, c
(
i
d
)
=
1, for i
d
=
0. F
h
(
v
r
) denotes yielded passive damping force with hysteresis depending on the piston relative displacement velocity (
v
r
) of the MRD for i
d
=
0. x
r is the piston travel of the MRD, y and z are inner variables without units, and k
0, k
1, k
2, a
0, I
0, α, β, γ, c
0, c
1, n, A, and x
0 are constants, respectively.Figure 2
Bouc-Wen mechanical model.As is shown in Figure3(a), a CARRERA MagneShock MRD of the vehicle suspension is further employed from [16], which permits maximum control direct current 0.5 A at 12 V [16]. On the basis of the actual measured characteristic data of MRD, the model parameters are identified as k
0
=
184.1, k
1
=
1528.1, k
2
=
10.092, a
0
=
7.526, I
0
=
0.069, α
=
20373.7, β
=
233849.1, γ
=
8816.9, c
0
=
1368.7, c
1
=
6222.7, n
=
2, A
=
20.6, and x
0
=
-
0.004. Figure 3(b) shows the comparison of the modified Bouc-Wen hysteresis F
-
v model and the tested data of employed candidate MRD, under the harmonic excitation with amplitude 12.5 mm at 1.5 Hz and different direct drive current (0–0.4 A), which exhibits ideal coordination.Figure 3
(a) A pictorial view of the candidate MRD. (b) Comparison of tested and calculated results.
(a)
(b)MRD has clear strong nonlinear property from the above results. Due to the damping force dependency on driven current and excitation frequency, theF
-
v curve varies evidently under different driven and excitation profiles [17]. As a result, nonlinear behavior caused by hysteresis cannot be ignored in engineering application.


## 3. Nonlinear Dynamics Characteristics of MR Suspension

Stability is undoubtedly the key requirements in the field of vehicle suspension control. In previous studies, comprehensive performances of the MR suspension were more highly valued, by comparing the transmissibility of vibration and vibration amplitude with the traditional passive one. Here, the nonlinear dynamic response and the mechanism of the instability process are the problem to be solved.In this section, the stability analysis of the MR suspension is conducted in view of modern control theory and nonlinear dynamics theory. Stability of the MR suspension system is examined at the equilibrium of the dimensionless nonlinear model. Further, the bifurcation analysis, combining LE assessment, is carried out to describe nonlinear dynamical evolution process and to find out the route to chaos. In addition, the phase portrait and Poincare map diagrams are plotted to give the more intuitive response of the system under different excitation conditions.


### 3.1. Stability Analysis

For definition of time coefficientτ
=
ω
·
t, in which ω
2
=
k
s
/
m
s, the corresponding dimensionless equation of the motion is written as [7, 10](6)
x
¨
s
+
k
s
m
s
ω
2
x
s
-
x
u
+
1
m
s
ω
F
d
=
0
,
x
¨
u
-
k
s
m
u
ω
2
x
s
-
x
u
-
1
m
u
ω
F
d
+
k
t
m
u
ω
2
x
u
-
X
i
+
c
t
m
u
ω
x
˙
u
-
X
i
'
=
0
,where X
i denotes dimensionless form of road profile. The state space is defined as x
=
[
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
], where x
1
=
x
s, x
2
=
x
˙
s, x
3
=
x
u, x
4
=
x
˙
u, x
5
=
y, and x
6
=
z. The state space is divided into four sections from (5), which are(7)
D
0
=
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
∣
α
x
6
+
c
0
x
2
-
x
4
+
k
0
x
1
-
x
3
-
x
5
<
0
∩
x
6
<
0
,
D
1
=
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
∣
α
x
6
+
c
0
x
2
-
x
4
+
k
0
x
1
-
x
3
-
x
5
>
0
∩
x
6
>
0
,
D
2
=
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
∣
α
x
6
+
c
0
x
2
-
x
4
+
k
0
x
1
-
x
3
-
x
5
>
0
∩
x
6
<
0
,
D
3
=
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
∣
α
x
6
+
c
0
x
2
-
x
4
+
k
0
x
1
-
x
3
-
x
5
<
0
∩
x
6
>
0
,respectively.Considering no excitation input, we get the state equation of the system further as follows from (1) to (7):(8)
x
˙
1
=
x
2
,
x
˙
2
=
a
21
x
1
+
a
22
x
2
+
a
23
x
3
+
a
24
x
4
+
a
25
x
5
+
a
26
x
6
,
x
˙
3
=
x
4
,
x
˙
4
=
a
41
x
1
+
a
42
x
2
+
a
43
x
3
+
a
44
x
4
+
a
45
x
5
+
a
46
x
6
,
x
˙
5
=
a
51
x
1
+
a
52
x
2
+
a
53
x
3
+
a
54
x
4
+
a
55
x
5
+
a
56
x
6
,
x
˙
6
=
A
+
B
x
6
2
·
a
61
x
1
+
a
62
x
2
+
a
63
x
3
+
a
64
x
4
+
a
65
x
5
+
a
66
x
6
,where(9)
a
21
=
-
k
s
m
s
ω
2
-
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
-
c
i
d
k
1
m
s
ω
,
a
22
=
-
c
i
d
m
s
ω
c
0
c
1
c
0
+
c
1
,
a
23
=
k
s
m
s
ω
2
+
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
+
c
i
d
k
1
m
s
ω
,
a
24
=
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
,
a
25
=
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
,
a
26
=
-
c
i
d
m
s
ω
α
c
1
c
0
+
c
1
,
a
41
=
k
s
m
u
ω
2
+
c
i
d
m
u
ω
k
0
c
1
c
0
+
c
1
+
c
i
d
k
1
m
u
ω
,
a
42
=
c
i
d
m
u
ω
c
0
c
1
c
0
+
c
1
,
a
43
=
-
k
s
m
u
ω
2
-
c
i
d
m
u
ω
k
0
c
1
c
0
+
c
1
-
c
i
d
k
1
m
u
ω
-
k
t
m
u
ω
2
,
a
44
=
-
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
+
c
t
m
u
ω
,
a
45
=
-
c
i
d
m
u
ω
k
0
c
1
c
0
+
c
1
,
a
46
=
c
i
d
m
u
ω
α
c
1
c
0
+
c
1
,
a
51
=
k
0
c
0
+
c
1
,
a
52
=
c
0
c
0
+
c
1
,
a
53
=
-
k
0
c
0
+
c
1
,
a
54
=
-
c
0
c
0
+
c
1
,
a
55
=
-
1
c
0
+
c
1
,
a
56
=
α
c
0
+
c
1
,
a
61
=
-
k
0
c
0
+
c
1
,
a
62
=
c
1
c
0
+
c
1
,
a
63
=
k
0
c
0
+
c
1
,
a
64
=
-
c
1
c
0
+
c
1
,
a
65
=
k
0
c
0
+
c
1
,
a
66
=
-
α
c
0
+
c
1
,
B
=
-
β
-
γ
for
x
∈
D
0
∪
D
1
,
B
=
-
β
+
γ
for
x
∈
D
2
∪
D
3
.The parameters of the system are selected based on actual vehicle [15]: m
s
=
562.5 kg, m
u
=
90 kg, k
s
=
57000 N/m, k
t
=
285000 N/m, and c
t
=
100 N/m·s
-
1. Regardless of semiactive control of MRD, the driven current is constant i
d
=
0.5 A, and other model parameters are as mentioned above. By setting the left side of (8) to zero, the equilibrium is obtained. Obviously, the system has a fixed point X
0
(
0,0
,
0,0
,
0,0
), the stability of which is determined by the characteristic equation of the Jacobian matrix (10) at X
0:(10)
J
=
0
-
λ
1
0
0
0
0
a
21
a
22
-
λ
a
23
a
24
a
25
a
26
0
0
0
-
λ
1
0
0
a
41
a
42
a
43
a
44
-
λ
a
45
a
46
a
51
a
52
a
53
a
54
a
55
-
λ
a
56
A
a
61
A
a
62
A
a
63
A
a
64
A
a
65
A
a
66
-
λ
.The eigenvalues of matrix (10) are calculated as λ
1
=
-
1.5176
+
2.3913
i, λ
2
=
-
1.5176
-
2.3913
i, λ
3
=
-
0.0099
+
0.2057
i, λ
4
=
-
0.0099
-
0.2057
i, λ
5
=
0.0052, and λ
6
=
0. Note that λ
1
,
2 and λ
3
,
4 are two couples of conjugate complex roots, while λ
5 is positive and λ
6 is zero, and therefore the system is instable at X
0.Furthermore, due to quadratic and absolute terms in the above six dimensions system, analytical solution is scarcely to obtain. Therefore, the numerical methods are usually employed to analyze such a system. Considering the application background, the range of frequencies and amplitudes of the road profile during vehicle running is certain, which are normally below 15 Hz and 10 cm, respectively [4–6]. We use the common harmonic excitation as the road surface. It is expressed as x
i
=
A
m
sin
⁡
(
Ω
·
t
), in which A
m represents roughness of the road surface and Ω is angular frequency Ω
=
2
π
f. The dimensionless form of road excitation is X
i
=
Ω
2
A
sin
⁡
(
Ω
/
ω
·
τ
). The following analysis methods of nonlinear dynamics are applied to such a special complex system:(1)
By obtaining the frequency band response, we find out the area that the system is sensitive to the corresponding road profiles.
(2)
By drawing the bifurcation diagrams, we study the nonlinear dynamics of the system under different road parameters. Combing the calculated LE spectrum, we can further determine the chaotic motion of the system under the corresponding road profiles.
(3)
Based on the above analysis, the dynamics evolution process is vividly portrayed using phase plane portraits, time series, and power spectrum of the system response at critical parameters.


### 3.2. Numerical Results

It is known that the dynamics of a vibration system may be analyzed through the frequency response diagram [18]. Therefore, for the studied system, the frequency response is obtained by plotting the vibration amplitude of sprung mass m
s expressed by v
s. In addition, the LE spectrum [19] is used to reveal the detail of the system frequency response. Figure 4(a) shows the frequency response of the model, which covers pass-band of the road, and Figure 4(b) presents the LE spectrum as well. The road frequency f is slowly increased by increment of 0.001 Hz. As is illustrated in Figure 4(a) there exists a critical jump of the system response for f
=
1.752 Hz [20] near the resonance point. The phenomenon of the jump causes the motion to change, but the system remains stable, as is shown in gray shadow. However, the diagram exhibits a more complicated and different behavior while f is increased to 2.71 Hz, with that restoration of stability in a short time. With f being increased, it is shown that the system falls into instability area as 2.75 Hz–4 Hz, which indicates that the chaotic motion may appear when f is within or near this area, as is shown in red shadow. This is confirmed by the LE spectrum showed in Figure 4(b). The diagram illustrates that there exist positive LE among the LE spectra in instable area shaded in red, which indicates the existence of the chaos. Next, the frequency response goes back to normal and all the LEs are less than zero, which indicates that the system remains stable up to 15 Hz.Figure 4
Frequency response diagram and LE spectrum diagram.
(a)
(b)The bifurcation and max LE diagrams under parameter variations are efficient methods for analyzing the nonlinear dynamic behavior. Figure5 shows the bifurcation diagram and the corresponding LE diagram, with f varying in the above-mentioned instable area. The bifurcation diagram is obtained by plotting the stroboscopic point of the displacement x
s, because x
s is important for the safety of the vehicle. As is shown in Figure 5(a), when f
=
2.08 Hz, the system loses period-1 stability, bifurcating from period-1 into period-2 motion. For f
∈
(2.08–2.74) Hz, the period-2 motion of the system develops into period-8 through the period-doubling bifurcation. Moreover, it is interesting that complex motion occurs at f
∈
(2.439–2.488) from the bifurcation diagram, and the max LEs are positive. At the same time, note that the points in the bifurcation diagram are limited in enveloping curve, not wide distribution. Consequently, the exact dynamic behavior in this area will be identified later. Next, after short stay in the strange area, the system gets back to periodic motion until f
=
2.719 Hz, through series of inverse period-doubling bifurcation from period-8 to period-2. Nevertheless, the system directly enters the chaotic state after period-2 motion, as is shown in Figure 5(b), and the max LE turns positive. When f
∈
(2.78–4) Hz, periodic motion, period-doubling bifurcation, and chaotic motions appear alternately and the system suffers the saddle-node bifurcation. Beyond the threshold for the onset of chaotic motion, there are some “periodic windows,” which could be the feature of the transient chaos. The max LEs before and after the saddle-node bifurcation point are opposite in sign. Then, for larger frequency over 3.73 Hz, there are no longer “periodic windows” presented up to the cutoff value of f in simulation. In Figure 5(c), we can observe that the system escapes from chaos to periodic motions through reverse period-doubling bifurcation.Figure 5
Bifurcation diagram and max LEs diagram ofx
s by varying f: (a) f
∈
(2.0–2.7) Hz; (b) f
∈
(2.7–4.0) Hz; (c) f
∈
(4.0–5.6) Hz.
(a)
(b)
(c)In order to give a rather clear presentation of the dynamical behavior, we depict phase portraits (with lines) and Poincare maps (with points) [21]. Figures 6(a)–6(h) show results for the following road frequency f: (a) f
=
1.5 Hz; (b) f
=
2.2 Hz; (c) f
=
2.39 Hz; (d) f
=
2.245 Hz; (e) f
=
2.471 Hz; (f) f
=
2.8 Hz; (g) f
=
3.071 Hz; (h) f
=
3.2 Hz; (i) f
=
3.75 Hz, and amplitude A
m is fixed at 0.08 m. The displacement x
s and velocity v
s are adopted. A periodic dynamic response is exhibited, including period-1, period-2, period-4, and period-8. In Figure 6(a), the period-1 motion is presented, according to the one loop circle as well as a single point of the Poincare map. The period-1 motion is replaced by period-2, period-4, and period-8 in turn, which confirms the analysis of process of bifurcating. Figure 6(e) shows the phase plane at f
=
2.471 Hz in the above-mentioned strange area. Note that the phase plane is keeping regular even consist of amount of closed curves, and the Poincare map contains limited points. In order to verify the dynamic behavior, in the above strange area, we plot the timing diagram and calculate the power spectrum density (PSD) diagram of x
s for f
=
2.4617 Hz. It is indicated that periodic and chaotic coexisting state exists when the system falls in the strange area [22], as is shown in Figure 7. Figure 6(f) reveals the chaotic attractor, because both of the phase plane and Poincare map distribute irregularly throughout the phase space. In the chaotic area, the saddle-node bifurcation is confirmed by phase plane portraits and Poincare maps, as is shown in Figures 6(g), 6(h), and 6(i). Note, as shown in Figure 6(i), a new chaotic attractor exists as f is increased, which covers a larger area compared with Figure 6(f). It means that the vibration amplitude rapidly increases with a small increase of road excitation frequency, which is harmful to driving safety.Figure 6
Phase portraits and corresponding Poincare maps (velocityv
s versus displacement x
s) for different f: (a) f
=
1.5 Hz; (b) f
=
2.2 Hz; (c) f
=
2.39 Hz; (d) f
=
2.245 Hz; (e) f
=
2.471 Hz; (f) f
=
2.8 Hz; (g) f
=
3.071 Hz; (h) f
=
3.2 Hz; (i) f
=
3.75 Hz.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 7
Timing history and PSD ofx
s for f
=
2.4617 Hz, A
m
=
0.08 m.
(a)
(b)The road surface amplitude is also essential in vehicle dynamics analysis. Figure8 illustrates the influence of road amplitude on the system dynamics, of the global bifurcation graph and LE spectrum under different amplitude of road excitation. The amplitude varies from 0.005 m to the max 0.1 m, containing the normal road condition, and frequency f is fixed at 5 Hz. It is observed that the displacement x
s and velocity v
s exhibit a periodic response, including period-1, period-2, and period-4 motion, until the road amplitude is increased to A
m
=
0.0315 m. Figure 9 expresses the phase portraits and Poincare map corresponding to A
m
=
0.01 m, 0.02316 m, 0.02720 m, and 0.06 m, respectively. The phase plane and Poincare map is plotted in Figures 9(a), 9(b), and 9(c), which confirms the process of period-doubling bifurcation shown in Figure 8(a). Note, for A
m
=
0.02720 m, a limit circle occurs, which corresponds to the phenomenon of jump in the process of period-2 shown in Figure 8(a). When the amplitude of road excitation is increased larger than A
m
=
0.0315 m, the system enters the chaotic region from Figure 8 because the max LE turns positive. In addition, just as is shown in Figure 9(d), the phase plane and Poincare map of A
m
=
0.06 m become irregular. The system maintains chaotic motion and vibration amplitude of suspension is getting far larger than the sustainable limit, which seriously threatens the handling safety of the vehicle.Figure 8
Bifurcation diagram and LE spectrum diagram ofx
s by varying A
m
∈
(0.005–0.1) m.
(a)
(b)Figure 9
Phase plane portraits and corresponding Poincare maps (velocityv
s versus displacement x
s) for different A
m: (a) A
m
=
0.01 m; (b) A
m
=
0.02316 m; (c) A
m
=
0.02720 m; (d) A
m
=
0.06 m.
(a)
(b)
(c)
(d)The results of analysis above are complete descriptions of the set of parameters where chaos occurs under road profiles. The proposed analytical method is useful for estimating the suspension parameters so that chaos does not occur as desired, and thus avoiding the unexpected dangers in the running process of the vehicle.


## 3.1. Stability Analysis

For definition of time coefficientτ
=
ω
·
t, in which ω
2
=
k
s
/
m
s, the corresponding dimensionless equation of the motion is written as [7, 10](6)
x
¨
s
+
k
s
m
s
ω
2
x
s
-
x
u
+
1
m
s
ω
F
d
=
0
,
x
¨
u
-
k
s
m
u
ω
2
x
s
-
x
u
-
1
m
u
ω
F
d
+
k
t
m
u
ω
2
x
u
-
X
i
+
c
t
m
u
ω
x
˙
u
-
X
i
'
=
0
,where X
i denotes dimensionless form of road profile. The state space is defined as x
=
[
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
], where x
1
=
x
s, x
2
=
x
˙
s, x
3
=
x
u, x
4
=
x
˙
u, x
5
=
y, and x
6
=
z. The state space is divided into four sections from (5), which are(7)
D
0
=
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
∣
α
x
6
+
c
0
x
2
-
x
4
+
k
0
x
1
-
x
3
-
x
5
<
0
∩
x
6
<
0
,
D
1
=
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
∣
α
x
6
+
c
0
x
2
-
x
4
+
k
0
x
1
-
x
3
-
x
5
>
0
∩
x
6
>
0
,
D
2
=
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
∣
α
x
6
+
c
0
x
2
-
x
4
+
k
0
x
1
-
x
3
-
x
5
>
0
∩
x
6
<
0
,
D
3
=
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
∣
α
x
6
+
c
0
x
2
-
x
4
+
k
0
x
1
-
x
3
-
x
5
<
0
∩
x
6
>
0
,respectively.Considering no excitation input, we get the state equation of the system further as follows from (1) to (7):(8)
x
˙
1
=
x
2
,
x
˙
2
=
a
21
x
1
+
a
22
x
2
+
a
23
x
3
+
a
24
x
4
+
a
25
x
5
+
a
26
x
6
,
x
˙
3
=
x
4
,
x
˙
4
=
a
41
x
1
+
a
42
x
2
+
a
43
x
3
+
a
44
x
4
+
a
45
x
5
+
a
46
x
6
,
x
˙
5
=
a
51
x
1
+
a
52
x
2
+
a
53
x
3
+
a
54
x
4
+
a
55
x
5
+
a
56
x
6
,
x
˙
6
=
A
+
B
x
6
2
·
a
61
x
1
+
a
62
x
2
+
a
63
x
3
+
a
64
x
4
+
a
65
x
5
+
a
66
x
6
,where(9)
a
21
=
-
k
s
m
s
ω
2
-
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
-
c
i
d
k
1
m
s
ω
,
a
22
=
-
c
i
d
m
s
ω
c
0
c
1
c
0
+
c
1
,
a
23
=
k
s
m
s
ω
2
+
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
+
c
i
d
k
1
m
s
ω
,
a
24
=
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
,
a
25
=
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
,
a
26
=
-
c
i
d
m
s
ω
α
c
1
c
0
+
c
1
,
a
41
=
k
s
m
u
ω
2
+
c
i
d
m
u
ω
k
0
c
1
c
0
+
c
1
+
c
i
d
k
1
m
u
ω
,
a
42
=
c
i
d
m
u
ω
c
0
c
1
c
0
+
c
1
,
a
43
=
-
k
s
m
u
ω
2
-
c
i
d
m
u
ω
k
0
c
1
c
0
+
c
1
-
c
i
d
k
1
m
u
ω
-
k
t
m
u
ω
2
,
a
44
=
-
c
i
d
m
s
ω
k
0
c
1
c
0
+
c
1
+
c
t
m
u
ω
,
a
45
=
-
c
i
d
m
u
ω
k
0
c
1
c
0
+
c
1
,
a
46
=
c
i
d
m
u
ω
α
c
1
c
0
+
c
1
,
a
51
=
k
0
c
0
+
c
1
,
a
52
=
c
0
c
0
+
c
1
,
a
53
=
-
k
0
c
0
+
c
1
,
a
54
=
-
c
0
c
0
+
c
1
,
a
55
=
-
1
c
0
+
c
1
,
a
56
=
α
c
0
+
c
1
,
a
61
=
-
k
0
c
0
+
c
1
,
a
62
=
c
1
c
0
+
c
1
,
a
63
=
k
0
c
0
+
c
1
,
a
64
=
-
c
1
c
0
+
c
1
,
a
65
=
k
0
c
0
+
c
1
,
a
66
=
-
α
c
0
+
c
1
,
B
=
-
β
-
γ
for
x
∈
D
0
∪
D
1
,
B
=
-
β
+
γ
for
x
∈
D
2
∪
D
3
.The parameters of the system are selected based on actual vehicle [15]: m
s
=
562.5 kg, m
u
=
90 kg, k
s
=
57000 N/m, k
t
=
285000 N/m, and c
t
=
100 N/m·s
-
1. Regardless of semiactive control of MRD, the driven current is constant i
d
=
0.5 A, and other model parameters are as mentioned above. By setting the left side of (8) to zero, the equilibrium is obtained. Obviously, the system has a fixed point X
0
(
0,0
,
0,0
,
0,0
), the stability of which is determined by the characteristic equation of the Jacobian matrix (10) at X
0:(10)
J
=
0
-
λ
1
0
0
0
0
a
21
a
22
-
λ
a
23
a
24
a
25
a
26
0
0
0
-
λ
1
0
0
a
41
a
42
a
43
a
44
-
λ
a
45
a
46
a
51
a
52
a
53
a
54
a
55
-
λ
a
56
A
a
61
A
a
62
A
a
63
A
a
64
A
a
65
A
a
66
-
λ
.The eigenvalues of matrix (10) are calculated as λ
1
=
-
1.5176
+
2.3913
i, λ
2
=
-
1.5176
-
2.3913
i, λ
3
=
-
0.0099
+
0.2057
i, λ
4
=
-
0.0099
-
0.2057
i, λ
5
=
0.0052, and λ
6
=
0. Note that λ
1
,
2 and λ
3
,
4 are two couples of conjugate complex roots, while λ
5 is positive and λ
6 is zero, and therefore the system is instable at X
0.Furthermore, due to quadratic and absolute terms in the above six dimensions system, analytical solution is scarcely to obtain. Therefore, the numerical methods are usually employed to analyze such a system. Considering the application background, the range of frequencies and amplitudes of the road profile during vehicle running is certain, which are normally below 15 Hz and 10 cm, respectively [4–6]. We use the common harmonic excitation as the road surface. It is expressed as x
i
=
A
m
sin
⁡
(
Ω
·
t
), in which A
m represents roughness of the road surface and Ω is angular frequency Ω
=
2
π
f. The dimensionless form of road excitation is X
i
=
Ω
2
A
sin
⁡
(
Ω
/
ω
·
τ
). The following analysis methods of nonlinear dynamics are applied to such a special complex system:(1)
By obtaining the frequency band response, we find out the area that the system is sensitive to the corresponding road profiles.
(2)
By drawing the bifurcation diagrams, we study the nonlinear dynamics of the system under different road parameters. Combing the calculated LE spectrum, we can further determine the chaotic motion of the system under the corresponding road profiles.
(3)
Based on the above analysis, the dynamics evolution process is vividly portrayed using phase plane portraits, time series, and power spectrum of the system response at critical parameters.


## 3.2. Numerical Results

It is known that the dynamics of a vibration system may be analyzed through the frequency response diagram [18]. Therefore, for the studied system, the frequency response is obtained by plotting the vibration amplitude of sprung mass m
s expressed by v
s. In addition, the LE spectrum [19] is used to reveal the detail of the system frequency response. Figure 4(a) shows the frequency response of the model, which covers pass-band of the road, and Figure 4(b) presents the LE spectrum as well. The road frequency f is slowly increased by increment of 0.001 Hz. As is illustrated in Figure 4(a) there exists a critical jump of the system response for f
=
1.752 Hz [20] near the resonance point. The phenomenon of the jump causes the motion to change, but the system remains stable, as is shown in gray shadow. However, the diagram exhibits a more complicated and different behavior while f is increased to 2.71 Hz, with that restoration of stability in a short time. With f being increased, it is shown that the system falls into instability area as 2.75 Hz–4 Hz, which indicates that the chaotic motion may appear when f is within or near this area, as is shown in red shadow. This is confirmed by the LE spectrum showed in Figure 4(b). The diagram illustrates that there exist positive LE among the LE spectra in instable area shaded in red, which indicates the existence of the chaos. Next, the frequency response goes back to normal and all the LEs are less than zero, which indicates that the system remains stable up to 15 Hz.Figure 4
Frequency response diagram and LE spectrum diagram.
(a)
(b)The bifurcation and max LE diagrams under parameter variations are efficient methods for analyzing the nonlinear dynamic behavior. Figure5 shows the bifurcation diagram and the corresponding LE diagram, with f varying in the above-mentioned instable area. The bifurcation diagram is obtained by plotting the stroboscopic point of the displacement x
s, because x
s is important for the safety of the vehicle. As is shown in Figure 5(a), when f
=
2.08 Hz, the system loses period-1 stability, bifurcating from period-1 into period-2 motion. For f
∈
(2.08–2.74) Hz, the period-2 motion of the system develops into period-8 through the period-doubling bifurcation. Moreover, it is interesting that complex motion occurs at f
∈
(2.439–2.488) from the bifurcation diagram, and the max LEs are positive. At the same time, note that the points in the bifurcation diagram are limited in enveloping curve, not wide distribution. Consequently, the exact dynamic behavior in this area will be identified later. Next, after short stay in the strange area, the system gets back to periodic motion until f
=
2.719 Hz, through series of inverse period-doubling bifurcation from period-8 to period-2. Nevertheless, the system directly enters the chaotic state after period-2 motion, as is shown in Figure 5(b), and the max LE turns positive. When f
∈
(2.78–4) Hz, periodic motion, period-doubling bifurcation, and chaotic motions appear alternately and the system suffers the saddle-node bifurcation. Beyond the threshold for the onset of chaotic motion, there are some “periodic windows,” which could be the feature of the transient chaos. The max LEs before and after the saddle-node bifurcation point are opposite in sign. Then, for larger frequency over 3.73 Hz, there are no longer “periodic windows” presented up to the cutoff value of f in simulation. In Figure 5(c), we can observe that the system escapes from chaos to periodic motions through reverse period-doubling bifurcation.Figure 5
Bifurcation diagram and max LEs diagram ofx
s by varying f: (a) f
∈
(2.0–2.7) Hz; (b) f
∈
(2.7–4.0) Hz; (c) f
∈
(4.0–5.6) Hz.
(a)
(b)
(c)In order to give a rather clear presentation of the dynamical behavior, we depict phase portraits (with lines) and Poincare maps (with points) [21]. Figures 6(a)–6(h) show results for the following road frequency f: (a) f
=
1.5 Hz; (b) f
=
2.2 Hz; (c) f
=
2.39 Hz; (d) f
=
2.245 Hz; (e) f
=
2.471 Hz; (f) f
=
2.8 Hz; (g) f
=
3.071 Hz; (h) f
=
3.2 Hz; (i) f
=
3.75 Hz, and amplitude A
m is fixed at 0.08 m. The displacement x
s and velocity v
s are adopted. A periodic dynamic response is exhibited, including period-1, period-2, period-4, and period-8. In Figure 6(a), the period-1 motion is presented, according to the one loop circle as well as a single point of the Poincare map. The period-1 motion is replaced by period-2, period-4, and period-8 in turn, which confirms the analysis of process of bifurcating. Figure 6(e) shows the phase plane at f
=
2.471 Hz in the above-mentioned strange area. Note that the phase plane is keeping regular even consist of amount of closed curves, and the Poincare map contains limited points. In order to verify the dynamic behavior, in the above strange area, we plot the timing diagram and calculate the power spectrum density (PSD) diagram of x
s for f
=
2.4617 Hz. It is indicated that periodic and chaotic coexisting state exists when the system falls in the strange area [22], as is shown in Figure 7. Figure 6(f) reveals the chaotic attractor, because both of the phase plane and Poincare map distribute irregularly throughout the phase space. In the chaotic area, the saddle-node bifurcation is confirmed by phase plane portraits and Poincare maps, as is shown in Figures 6(g), 6(h), and 6(i). Note, as shown in Figure 6(i), a new chaotic attractor exists as f is increased, which covers a larger area compared with Figure 6(f). It means that the vibration amplitude rapidly increases with a small increase of road excitation frequency, which is harmful to driving safety.Figure 6
Phase portraits and corresponding Poincare maps (velocityv
s versus displacement x
s) for different f: (a) f
=
1.5 Hz; (b) f
=
2.2 Hz; (c) f
=
2.39 Hz; (d) f
=
2.245 Hz; (e) f
=
2.471 Hz; (f) f
=
2.8 Hz; (g) f
=
3.071 Hz; (h) f
=
3.2 Hz; (i) f
=
3.75 Hz.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 7
Timing history and PSD ofx
s for f
=
2.4617 Hz, A
m
=
0.08 m.
(a)
(b)The road surface amplitude is also essential in vehicle dynamics analysis. Figure8 illustrates the influence of road amplitude on the system dynamics, of the global bifurcation graph and LE spectrum under different amplitude of road excitation. The amplitude varies from 0.005 m to the max 0.1 m, containing the normal road condition, and frequency f is fixed at 5 Hz. It is observed that the displacement x
s and velocity v
s exhibit a periodic response, including period-1, period-2, and period-4 motion, until the road amplitude is increased to A
m
=
0.0315 m. Figure 9 expresses the phase portraits and Poincare map corresponding to A
m
=
0.01 m, 0.02316 m, 0.02720 m, and 0.06 m, respectively. The phase plane and Poincare map is plotted in Figures 9(a), 9(b), and 9(c), which confirms the process of period-doubling bifurcation shown in Figure 8(a). Note, for A
m
=
0.02720 m, a limit circle occurs, which corresponds to the phenomenon of jump in the process of period-2 shown in Figure 8(a). When the amplitude of road excitation is increased larger than A
m
=
0.0315 m, the system enters the chaotic region from Figure 8 because the max LE turns positive. In addition, just as is shown in Figure 9(d), the phase plane and Poincare map of A
m
=
0.06 m become irregular. The system maintains chaotic motion and vibration amplitude of suspension is getting far larger than the sustainable limit, which seriously threatens the handling safety of the vehicle.Figure 8
Bifurcation diagram and LE spectrum diagram ofx
s by varying A
m
∈
(0.005–0.1) m.
(a)
(b)Figure 9
Phase plane portraits and corresponding Poincare maps (velocityv
s versus displacement x
s) for different A
m: (a) A
m
=
0.01 m; (b) A
m
=
0.02316 m; (c) A
m
=
0.02720 m; (d) A
m
=
0.06 m.
(a)
(b)
(c)
(d)The results of analysis above are complete descriptions of the set of parameters where chaos occurs under road profiles. The proposed analytical method is useful for estimating the suspension parameters so that chaos does not occur as desired, and thus avoiding the unexpected dangers in the running process of the vehicle.


## 4. Conclusions

This paper studies the nonlinear dynamics problems in the practical application of MRD. The 2-DoF MR vehicle suspension was set up based on the identified data of a commercial MRD. By calculating the eigenvalues of the Jacobian matrix at fixed point, the possibility of chaotic movement in the system was discovered. The theoretical analysis was then confirmed by numerical simulations. Under the single frequency harmonic excitation, the nonlinear dynamic evolution process was analyzed with frequency varying. The main conclusion would be that the loss of stability of the system appears near the midfrequency band and high amplitude of road surface, through the dynamical evolution to chaos by period-doubling bifurcations, saddle-node bifurcations, and reverse period-doubling bifurcations. These show the importance of parametric excitation in the control of vibrations in a MR suspension system. The research provides the reference for nonlinear dynamic analysis and control method in engineering application.



---

*Source: 971731-2015-11-08.xml*