Bifurcation and chaos analysis of the floating raft vibration isolation system with quasi-zero-stiffness isolators

1. Introduction

Power equipment usually causes vibration and noise. Passive vibration isolation is widely used to solve the problem of vibration transmission in power equipment due to its economical and reliable characteristics [1]. At present, passive vibration isolation devices have good vibration isolation effects in the medium and high frequency ranges, but the vibration isolation effects in the low frequency and even ultra-low frequency ranges are not ideal [2]. The shortcomings of the vibration isolation system can be compensated by introducing nonlinear vibration isolation elements with high static and low dynamic stiffness [3]. Among them, the quasi-zero stiffness vibration isolator introduces a cubic stiffness term, which can significantly improve the vibration isolation performance of the vibration isolation system in the low-frequency region.

Related research on quasi-zero-stiffness vibration isolators focuses on their geometric configuration and vibration isolation performance. References [4-6] mainly improved the vibration isolation performance in the low-frequency region by improving the structure of the single-degree-of-freedom quasi-zero-stiffness vibration isolator. Jazer et al. [7] studied a single-degree-of-freedom vibration isolation system with dual nonlinearities in stiffness and damping and found that the effect of nonlinear stiffness on the frequency response of the system is opposite to that of nonlinear damping. Xu et al. [8] proposed a multi-directional quasi-zero-stiffness isolator with time-delay active control, which can suppress low-frequency vibrations in multiple directions. Kawana et al. [9] designed a three-degree-of-freedom vibration isolator that can suppress the vertical, horizontal, and rotational vibrations of a rigid body, and found that the transient nonlinear vibration of the rigid body has a 2:1:1 internal resonance. Lu et al. [10] and Tang et al. [11] respectively designed a quasi-zero-stiffness vibration isolation platform with a 6-degree-of-freedom rigid body. The results showed that the introduction of quasi-zero-stiffness vibration isolators can improve the low-frequency vibration isolation characteristics. The quasi-zero-stiffness vibration isolators improve the vibration isolation performance remarkably in the low-frequency region. However, the chaotic phenomena make the vibration state more complex.

The quasi-zero-stiffness isolators makes the system more susceptible to chaotic phenomena that are unique to nonlinear systems, especially under the influence of large low-frequency disturbances [12]. Lou et al. [13] evaluated the vibration isolation performance of nonlinear vibration isolators in a chaotic state and found that the line spectrum amplitude was significantly reduced. Li et al. [14] and Zhang et al. [15] used time-delay control methods to deal with chaos in a two-dimensional vibration isolation raft system and studied the stability of different time-delay feedback control methods. Based on a simplified two-dimensional raft system model, Zhang et al. [16] fully presented the line spectrum chaos method of nonlinear time-delay feedback control and explored the influence of control parameters on the chaotic effect of the system under multi-source excitation conditions. Chai et al. [17] used double-delay feedback control to make a high-static and low-dynamic vibration isolation system chaotic, while weakening the line spectrum characteristics and intensity. Zuo et al. [18] used the open-loop plus nonlinear closed-loop coupling method to realize the generalized chaotic synchronization, which can achieve the dual functions of concealing the line spectrum information and maintaining the vibration isolation performance. Obviously, when quasi-zero-stiffness isolators are applied to vibration isolation systems, it is inevitable to pay attention to their nonlinear vibration phenomena.

Existing research on quasi-zero-stiffness vibration isolation systems focuses on single-degree-of-freedom vibration isolation or double-layer vibration isolation [19-23]. Lu et al. [24] placed quasi-zero-stiffness isolators on the upper and lower layers of a double-layer vibration isolation system and found that the stiffness nonlinearity of the lower layer had a greater impact on the vibration isolation performance of the system than that of the upper layer. Ji et al. [25] explored the vibration isolation performance of two-degree-of-freedom quasi-zero stiffness isolators under different parameters. The study found that the damping ratio and the vertical stiffness ratio can improve the low-frequency vibration isolation performance. Zhao et al. [26] established a vertical vehicle-quasi-zero stiffness floating plate track coupling dynamic model and found that excessively high nonlinear stiffness levels would weaken the vibration isolation effect in the resonance and high-frequency regions. Zhang et al. [27] established the raft dynamic equation considering the foundation flexibility and the nonlinear stiffness of the isolator, and analyzed the vertical vibration of the system under single-frequency and multi-frequency excitation. Zhao et al. [28] established the dynamic equations of the airbag raft-rotor system considering the vertical and horizontal vibrations of the raft body, and analyzed the nonlinear characteristics of the system such as bifurcation and chaos under a single vibration excitation source. Li et al. [29] introduced the quasi-zero-stiffness isolator into the floating raft vibration isolation system and established a high-dimensional coupled dynamic equation with 12 degrees of freedom. They found that its vibration isolation performance was significantly improved, but did not analyze the nonlinear characteristics of the system. Li et al. [30] established a 6-degree-of-freedom vibration isolation model composed of eight isolators to solve the complex problem of rolling and pitching motions of ocean-going vessels. Wei et al. [31] designed a 2-degree-of-freedom bottom-springs grounded nonlinear vibration isolator model, and carried out numerical simulation under typical harmonic excitation. The results showed that the model has a good vibration isolation effect. Sun et al. [32] adopted the improved Incremental Harmonic Balance Method to discuss the coupling relationship between different degrees of freedom of the quasi-zero stiffness Gough-Stewart vibration isolation platform, considering the coupling relationship between each degree of freedom.

However, it is worth mentioning that researchers’ works mainly concentrated on either single-degree-of-freedom or double-layer vibration isolation while the floating raft vibration isolation with quasi-zero-stiffness isolators is not referred. In addition, the excitation resource from their research is always a single vibration excitation source, and the definition of excitation was not accurate enough to describe the real movement of the floating raft.

In this paper, an investigation into the nonlinear dynamic behaviors of the floating raft isolation system coupled with quasi-zero-stiffness isolators under multiple disturbance sources is presented. First, the coupling effects between the excitation source and isolation system are considered, and movements in the vertical, horizontal, and rotational directions are also taken into account. Based on this, the floating raft isolation model subjected to multiple disturbance sources and its motion equation is developed, and then its dynamic responses are mainly investigated by using the techniques of time history diagram, spectrum diagram, phase diagram, Poincaré diagram, and the bifurcation diagram. Last, the bifurcations of the floating raft isolation system with different parameters are analyzed through numerical methods, especially the effect of center distance and mass ratio. In particular, this work aims to present the theoretic supporting for the dynamics, vibration control and its parametric optimization of the floating raft isolation system coupled with quasi-zero-stiffness isolators.

2. Modeling of the floating raft vibration isolation system

2.1. Dynamics modeling

A plane motion model is established for the floating raft vibration isolation system considering multi-disturbance sources, as shown in Fig. 1. The plane motion model includes five parts: machinery A, machinery B, upper vibration isolator I, lower vibration isolator II, and raft C.

Fig. 1Plane motion model of the floating raft isolation system

The floating raft vibration isolation system is modeled by the multi-rigid body dynamic method. The machinery and raft are regarded as rigid bodies. A couple of quasi-zero-stiffness vibration isolators are vertically installed under each machinery, and three quasi-zero-stiffness vibration isolators are vertically installed under the raft. The foundation is rigidly fixed. As shown in Fig. 1, the mass centers of machinery A, machinery B and raft C are set as the origins of their local coordinate systems $o_A$, $o_B$, $o_C$, respectively. $o_A{x}_A$, $o_B{x}_B$, $o_C{x}_C$ are the ordinate axes of the local coordinate systems, and their directions are perpendicular to the geodetic plane where the system is located. $o_A{y}_A$, $o_B{y}_B$, $o_C{y}_C$ are the horizontal axes of the local coordinate systems, and their directions are parallel to the geodetic plane. $\varphi A$, $\varphi B$, $\varphi C$ represents the rotation angles of machinery A, machinery B and raft C in the plane, respectively. Set the origin of the global coordinate system at point $\mathrm{o}$. $o_x$, $o_y$ are the ordinate axis and horizontal axis of the overall coordinate system.

The dynamic equations of machinery A, machinery B and raft C are established respectively based on Newton's second law, which is formulated in Eq. (1):

1

$\left\{\begin{array}{l}{\mathbf{M}}_{\mathrm{A}}{\ddot{\mathbf{r}}}_{\mathrm{A}}={\mathbf{R}}_{\mathrm{A}}{\mathbf{F}}_{\mathrm{A}}-\sum _{i=1}^2{\mathbf{F}}_{\mathrm{{\rm I}}i}^t,\\ {\mathbf{M}}_{\mathrm{B}}{\ddot{\mathbf{r}}}_{\mathrm{B}}={\mathbf{R}}_{\mathrm{B}}{\mathbf{F}}_{\mathrm{B}}-\sum _{i=3}^4{\mathbf{F}}_{\mathrm{{\rm I}}i}^t,\\ {\mathbf{M}}_{\mathrm{C}}{\ddot{\mathbf{r}}}_{\mathrm{C}}=-\sum _{j=1}^3{\mathbf{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t-\sum _{i=1}^4{\mathbf{F}}_{\mathrm{{\rm I}}i}^b.\end{array}\mathrm{}\right.$

In the equilibrium equation of machinery A, ${\mathbf{M}}_{\mathrm{A}}=\mathrm{d}\mathrm{i}\mathrm{g}[M_{\mathrm{A}},M_{\mathrm{A}},J_{\mathrm{A}}]$ represents the mass matrix of machinery A. $M_{\mathrm{A}}$ and $J_{\mathrm{A}}$ represent the mass and planar moment of inertia, respectively. ${\ddot{\mathbf{r}}}_{\mathrm{A}}=[{\ddot{x}}_{\mathrm{A}},{\ddot{y}}_{\mathrm{A}},{\ddot{\varphi }}_{\mathrm{A}}{]}^T$ represents the acceleration vector of the center of machinery A. ${\mathbf{R}}_{\mathrm{A}}$ represents the transformation matrix of the local coordinate system to the overall coordinate system for machinery A, which can be expressed as Eq. (2). ${\mathbf{F}}_{\mathrm{A}}=[F_{\mathrm{A}}\mathrm{c}\mathrm{o}\mathrm{s}{\omega }_{\mathrm{A}}t,\mathrm{0,0}{]}^T$ represents the simple harmonic excitation vector of machinery A. Assuming the machinery is only subjected to vertical excitation, the excitation amplitude is $F_{\mathrm{A}}$ and the excitation frequency is ${\omega }_{\mathrm{A}}$. ${\mathbf{F}}_{\mathrm{{\rm I}}i}^t$ represents the restoring forces at the top ends of the upper vibration isolator I. The subscript $i$ indicates the $i$th vibration isolator in the upper level I. Superscript $t$ represents the top end of the isolator:

2

${\mathbf{R}}_{\mathrm{A}}=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_{\mathrm{A}}& -\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_{\mathrm{A}}& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_{\mathrm{A}}& \mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_{\mathrm{A}}& 0\\ 0& 0& 1\end{array}\right].$

Similarly, in the equilibrium equation of machinery B, ${\mathbf{M}}_{\mathrm{B}}=\mathrm{d}\mathrm{i}\mathrm{g}[M_{\mathrm{B}},M_{\mathrm{B}},J_{\mathrm{B}}]\text{,}$${\ddot{\mathbf{r}}}_{\mathrm{B}}=[{\ddot{x}}_{\mathrm{B}},{\ddot{y}}_{\mathrm{B}},{\ddot{\varphi }}_{\mathrm{B}}{]}^T$, ${\mathbf{F}}_{\mathrm{B}}=\left[F_{\mathrm{B}}\mathrm{c}\mathrm{o}\mathrm{s}\right({\omega }_{\mathrm{B}}t+\phi ),\mathrm{0,0}{]}^T$ represent the mass matrix, acceleration array and simple harmonic excitation array of machinery B, respectively. $\phi$ represents the phase difference between motivation ${\mathbf{F}}_{\mathrm{B}}$ and ${\mathbf{F}}_{\mathrm{A}}$. ${\mathbf{R}}_{\mathrm{B}}$ represents the transformation matrix of the local coordinate system to the overall coordinate system for machinery B, and it is formulated in the same form as ${\mathbf{R}}_{\mathrm{A}}$.

Similarly, in the equilibrium equation of raft C, ${\mathbf{M}}_{\mathrm{C}}=\mathrm{d}\mathrm{i}\mathrm{g}[M_{\mathrm{C}},M_{\mathrm{C}},J_{\mathrm{C}}]$ and ${\ddot{\mathbf{r}}}_{\mathrm{C}}=[{\ddot{x}}_{\mathrm{C}},{\ddot{y}}_{\mathrm{C}},{\ddot{\varphi }}_{\mathrm{C}}{]}^T$ represent the mass matrix and acceleration vector of raft C, respectively. ${\mathbf{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t$ represents the restoring forces acting on the top end of the lower vibration isolator II. ${\mathbf{F}}_{\mathrm{{\rm I}}i}^b$ represents the restoring forces at the bottom end of the upper vibration isolator I. The subscript $j$ indicates the $j$th vibration isolator in the lower level II. Superscript $b$ represents the bottom end of the isolator.

Both the upper vibration isolator I and the lower vibration isolator II are selected as quasi-zero-stiffness vibration isolators. Their restoring forces ${\mathbf{F}}_{\mathrm{{\rm I}}i}^t$, ${\mathbf{F}}_{\mathrm{{\rm I}}i}^b$ and ${\mathbf{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t$ can be expressed as $k\mathrm{\Delta }x(1+\epsilon \mathrm{\Delta }{x}^2)$, where $\mathrm{\Delta }x$ is the deformation, $k$ is the stiffness coefficient, $\epsilon$ is the cubic nonlinear stiffness coefficient.

The stiffness coefficient matrix and damping coefficient matrix of the upper vibration isolator I are ${\mathbf{K}}_{\mathrm{{\rm I}}}=\mathrm{d}\mathrm{i}\mathrm{g}[k_{\mathrm{{\rm I}}x},k_{\mathrm{{\rm I}}y},k_{\mathrm{{\rm I}}\varphi }]$, ${\mathbf{C}}_{\mathrm{{\rm I}}}=\mathrm{d}\mathrm{i}\mathrm{g}[c_{\mathrm{{\rm I}}x},c_{\mathrm{{\rm I}}y},c_{\mathrm{{\rm I}}\varphi }]$. $k_{\mathrm{{\rm I}}x}$, $k_{\mathrm{{\rm I}}y}$ and $k_{\mathrm{{\rm I}}\varphi }$ are the stiffness coefficients of the upper vibration isolator I in the vertical, horizontal and rotational directions. Accordingly, $c_{\mathrm{{\rm I}}x}$, $c_{\mathrm{{\rm I}}y}$ and $c_{\mathrm{{\rm I}}\varphi }$ are the damping coefficients of the upper vibration isolator I in the vertical, horizontal and rotational directions.

The stiffness coefficient matrix and damping coefficient matrix of the lower vibration isolator II are ${\mathbf{K}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}=\mathrm{d}\mathrm{i}\mathrm{g}[k_{\mathrm{{\rm I}}\mathrm{{\rm I}}x},k_{\mathrm{{\rm I}}\mathrm{{\rm I}}y},k_{\mathrm{{\rm I}}\mathrm{{\rm I}}\varphi }]$, ${\mathbf{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}=\mathrm{d}\mathrm{i}\mathrm{g}[c_{\mathrm{{\rm I}}\mathrm{{\rm I}}x},c_{\mathrm{{\rm I}}\mathrm{{\rm I}}y},c_{\mathrm{{\rm I}}\mathrm{{\rm I}}\varphi }]$. $k_{\mathrm{{\rm I}}\mathrm{{\rm I}}x}$, $k_{\mathrm{{\rm I}}\mathrm{{\rm I}}y}$ and $k_{\mathrm{{\rm I}}\mathrm{{\rm I}}\varphi }$ are the stiffness coefficients of the lower vibration isolator II in the vertical, horizontal and rotational directions. Accordingly, $c_{\mathrm{{\rm I}}\mathrm{{\rm I}}x}$, $c_{\mathrm{{\rm I}}\mathrm{{\rm I}}y}$ and $c_{\mathrm{{\rm I}}\mathrm{{\rm I}}\varphi }$ are the damping coefficients of the lower vibration isolator II in the vertical, horizontal and rotational directions.

At the same time, the frequency parameters corresponding to the stiffness matrix can be obtained, as shown in Table 1.

Table 1Frequency parameters

Parameter	Expression	Parameter	Expression
${\omega }_{1x}$	$\sqrt{2k_{\mathrm{{\rm I}}x}/M_{\mathrm{A}}}$	${\omega }_{1y}$	$\sqrt{2k_{\mathrm{{\rm I}}y}/M_{\mathrm{A}}}$
${\omega }_{2x}$	$\sqrt{2k_{\mathrm{{\rm I}}x}/M_{\mathrm{B}}}$	${\omega }_{2y}$	$\sqrt{2k_{\mathrm{{\rm I}}y}/M_{\mathrm{B}}}$
${\omega }_{3x}$	$\sqrt{3{\mathrm{k}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}x}/M_{\mathrm{C}}}$	${\mathrm{\omega }}_{3\mathrm{y}}$	$\sqrt{3{\mathrm{k}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}y}/M_{\mathrm{C}}}$

The equilibrium equations of the upper vibration isolator I and the lower vibration isolator II are formulated. And their restoring forces are expressed as:

3

$\left\{\begin{array}{l}{\mathbf{F}}_{\mathrm{{\rm I}}i}^t={\mathbf{T}}_{\mathrm{{\rm I}}i}^t{\mathbf{K}}_{\mathrm{{\rm I}}}[\left(\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^t-\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^b\right)+\epsilon \left(\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^t-\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^b{)}^3\right]+{\mathbf{T}}_{\mathrm{{\rm I}}i}^t{\mathbf{C}}_{\mathrm{{\rm I}}}\left(\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i}^t-\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i}^b\right),\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}i=\mathrm{1,2},\mathrm{3,4},\\ {\mathbf{F}}_{\mathrm{{\rm I}}i}^b=-{\mathbf{T}}_{\mathrm{C}i}^b{\mathbf{K}}_{\mathrm{{\rm I}}}[\left(\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^t-\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^b\right)+\epsilon \left(\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^t-\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^b{)}^3\right]-{\mathbf{T}}_{\mathrm{C}i}^b{\mathbf{C}}_{\mathrm{{\rm I}}}\left(\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i}^t-\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i}^b\right),\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}i=\mathrm{1,2},\mathrm{3,4},\\ {\mathbf{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t={\mathbf{T}}_{\mathrm{C}j}^t{\mathbf{K}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}[\left(\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t-\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^b\right)+\epsilon \left(\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t-\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^b{)}^3\right]+{\mathbf{T}}_{\mathrm{C}j}^t{\mathbf{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\left(\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t-\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^b\right),\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}j=\mathrm{1,2},3.\end{array}\right.$

Among them:

4

${\mathbf{T}}_{\mathrm{{\rm I}}i}^t=\left\{\begin{array}{l}{\mathbf{T}}_{\mathrm{A}i}^t,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}i=\mathrm{1,2},\\ {\mathbf{T}}_{\mathrm{B}i}^t,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}i=\mathrm{3,4},\end{array}\right.$

where, ${\mathbf{T}}_{\mathrm{A}i}^t$ and ${\mathbf{T}}_{\mathrm{B}i}^t$ represents the translation transformation matrix of the top end of the upper vibration isolator I in the global coordinate system. The specific form is expressed as follows:

5

${\mathbf{T}}_{\mathrm{A}1}^t=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ -y_{\mathrm{{\rm I}}i}^{\mathrm{t}}& x_{\mathrm{{\rm I}}i}^{\mathrm{t}}& 1\end{array}\right]\mathrm{},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}i=\mathrm{1,2},$

where $x_{\mathrm{{\rm I}}i}^t$ and $y_{\mathrm{{\rm I}}i}^t$ represents the coordinates of the upper end of the upper isolator I in the local coordinate system. The meaning and form of ${\mathbf{T}}_{\mathrm{B}i}^t$, ${\mathbf{T}}_{\mathrm{C}i}^b$, and ${\mathbf{T}}_{\mathrm{C}j}^t$ are similar to ${\mathbf{T}}_{\mathrm{A}i}^t$, which needs no further elaboration:

6

$\left\{\begin{array}{l}\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^t={\mathbf{r}}_{\mathrm{{\rm I}}i}^t-{\mathbf{r}}_{\mathrm{{\rm I}}i0}^t,\\ \mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i}^t={\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i}^t-{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i0}^t,\end{array}\right.\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}i=\mathrm{1,2},\mathrm{3,4},$

where ${\mathbf{r}}_{\mathrm{{\rm I}}i0}^t$ and ${\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i0}^t$ represents the deformation and velocity when the coordinates of the top end of the upper isolator I are in the local coordinate system. $\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^b$, $\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t$, $\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^b$, $\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i}^b$, $\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^t$, and $\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j}^b$ have similar meaning and form as $\mathrm{\Delta }{\mathbf{r}}_{\mathrm{{\rm I}}i}^t$ and $\mathrm{\Delta }{\dot{\mathbf{r}}}_{\mathrm{{\rm I}}i}^t$, which needs no further elaboration.

3. Bifurcation characteristics of the floating raft vibration isolation system

Parameter analysis is carried on bifurcation characteristics of the floating raft vibration isolation system with quasi-zero-stiffness. Since machineries A and B are arranged symmetrically, only the movement of machinery A and raft C is considered. The fourth-order Runge-Kutta method is used to solve the nonlinear dynamic equations, with the time step is specified as $T/125$.

3.1. The influence of the distance between the centroids

Centroid distance $s$ is a dimensionless parameter reflecting the installation distance between machinery A and machinery B. Fig. 2 shows the bifurcation diagram of the steady-state displacement response of machinery A as the centroid distance varies.

Fig. 2Bifurcation diagrams of machinery A with s

a) Vertical response bifurcation diagram of machinery A

b) Horizontal response bifurcation diagram of machinery A

c) Rotation response bifurcation diagram of machinery A

As shown in Fig. 2, the motions of machinery A in vertical and horizontal direction have close trend with the centroid distance $s$. Since the excitation direction is vertical, the vibration amplitude of vertical motion is larger. When the centroid distance $s$> 19.68, the motion of machinery A shows an alternating state of periodic motion and chaotic motion. When the centroid distance$s$is between 19.94 and 20.02, the motion of machinery A in three directions transforms into quasi-periodic motion. As the centroid distance$s$continues to increase, the motion goes through a period of alternating periodic motion and chaotic motion, and then transforms into periodic motion. It is worth noting that machinery A has no single-period motion. That is because the vertical and horizontal directions of machinery A provide coupling stiffness for the rotation direction. The rotation response of machinery A is affected by the movement in the other two directions, and the movement state is more complicated.

In general, under the motivation of vertical equal-amplitude double excitation with a phase difference of $\pi$/2, the motion of machinery A conforms to the law of period-chaos-period according to the centroid distance$s$.

Fig. 3Bifurcation diagrams of raft C with s

a) Vertical response bifurcation diagram of floating raft C

b) Horizontal response bifurcation diagram of floating raft C

c) Rotation response bifurcation diagram of floating raft C

Fig. 3 is the bifurcation diagram of the steady-state displacement response of raft C as the distance from the mass center varies. Similar to the movement of machinery A, the movement of raft C in all directions varies according to the centroid distance$s$. The variation rules are basically similar With the centroid distance$s$, the motion of the raft C in all directions first shows a stable periodic motion, followed by an alternating state of periodic motion, quasi-periodic motion, multi-periodic motion, and chaotic motion. Compared with machinery A, the vibration amplitude of raft C in all directions is smaller. However, in the unstable quasi-periodic motion state, the bifurcation of the raft C is more complex and the nonlinear vibration characteristics are more obvious.

In summary, when the value of the centroid distance$s$ approaches 20, the movement of the floating raft vibration isolation system will be greatly affected. The system movement will produce more obvious quasi-periodic motion and chaotic motion, and the movement in the three directions will become more unstable. Especially for horizontal movement and rotation, the vibration amplitude will be greatly increased.

Therefore, in order to avoid the nonlinear motion state in the process of system motion, this value range $s=$19.6-20.2 should be excluded in the design of floating raft isolation system. It can be seen from the above that $s=L/x_0$, $x_0=M_{\mathrm{A}}g/2k_{\mathrm{{\rm I}}x}$, so $s$ can be optimized by designing the vertical stiffness coefficient $k_{\mathrm{{\rm I}}x}$ of the upper vibration isolator.

3.2. The influence of the mass ratio

Mass ratio $\eta$ is a dimensionless parameter reflecting the mass ratio between raft C and machinery A. Fig. 4 is the bifurcation diagram of the steady-state displacement response of machinery A as the mass ratio varies.

Fig. 4Bifurcation diagrams of machinery A with η

a) Vertical response bifurcation diagram of machinery A

b) Horizontal response bifurcation diagram of machinery A

c) Rotation response bifurcation diagram of machinery A

It can be seen from the Fig. 4 that the displacement response of machinery A in each direction varies with the mass ratio $\eta$. There are multiple discontinuous quasi-periodic motion intervals and chaotic motion intervals in the variation. When the mass ratio $\eta <$0.890, the motion state of machinery A is stable. The motion of machinery A comes into a quasi-periodic state when the mass ratio $\eta$ is between 0.716 and 0.750 and at 0.766. It is in periodic motion in the rest of the value range. When the mass ratio $\eta$ is in the range of 0.892 and 1.398, the motion of machinery A in all directions entered an obvious nonlinear motion state. In particular, its horizontal response entered chaos through period-doubling bifurcation, and the vibration amplitudes in all directions oscillated over a large range. Quasi-periodic motion states and chaotic motion states appears alternately, and the vibration amplitude in the horizontal direction reaches the same order of magnitude as that in the vertical direction. After the mass ratio $\eta$ exceeding 1.400, quasi-periodic motion and chaotic motion states only appear between 1.530-1.606 and 1.944-2.000, and the rest of the value range is periodic motion.

In general, the motion of machinery A shows a pattern of alternating periodic-quasi-periodic-chaotic appearance according to the increment of mass ratio $\eta$.

Fig. 5 is the bifurcation diagram of the steady-state displacement response of raft C as the mass ratio $\eta$ changes. Similar to the movement of machinery A, the movement of raft C in all directions varies with the mass ratio. As the mass ratio $\eta$ increases, the motion of the raft C in all directions first shows a stable periodic motion, followed by periodic motion, quasi-periodic motion, multi-periodic motion, and chaotic motion. Compared with machinery A, the vibration amplitude of raft C in all directions is smaller. At the same time, the rotation response of raft C also has the phenomenon of entering a chaotic state by period-doubling bifurcation.

Fig. 4(c) and Fig. 5(c) reveal that the rotational response of machinery A and raft C varies with the mass ratio $\eta$. In the bifurcation characteristics, there is also a period-doubling bifurcation phenomena, which is also caused by the coupling stiffness in the rotation direction.

In summary, in order to avoid the nonlinear motion state during the system motion, the mass ratio $\eta$ should avoid the value ranges of 0.892-1.398, 1.530-1.606 and 1.944-2.000. The mass ratio $\eta$ is optimized by designing the raft mass $M_{\mathrm{C}}$ and the machinery A mass $M_{\mathrm{A}}$.

4. Chaotic characteristics of the floating raft vibration isolation system

In nonlinear systems, the Lyapunov exponent is a quantitative analysis method to identify chaotic motion. When the maximum Lyapunov exponent is positive, the system is in a chaotic state; otherwise, the system is in a periodic or quasi-periodic motion state. From the analysis in Section 3.2, we can see that the mass ratio, $\eta$ has more obvious and complex impact on the vibration characteristics of the system. Therefore, it is necessary to use the Lyapunov exponent to further analyze the nonlinear characteristics of the system.

Fig. 5Bifurcation diagrams of raft C with η

a) Vertical response bifurcation diagram of floating raft C

b) Horizontal response bifurcation diagram of floating raft C

c) Rotation response bifurcation diagram of floating raft C

This section adopts the Benettin method [33]to calculate the maximum Lyapunov index of the floating raft vibration isolation system based on the orbit perturbation principle. The initial condition of the phase trajectory is ${\mathbf{z}}_0=0_{9\times 1}$, ${\mathbf{z}\mathit{\text{'}}}_0=0_{9\times 1}$, ${\mathbf{z}}_1=[10^{-7},0_{1\times 8}{]}^T$. Fig. 6 plots the curve of maximum Lyapunov exponent versus mass ratio $\eta$ . By comparing Fig. 6 with Fig. 4, it can be found that the maximum Lyapunov exponent of the system is basically consistent with the period-doubling bifurcation diagram. As the mass ratio increases, when the system motion is in a periodic motion state, the maximum Lyapunov exponent is negative, when the system motion is in a chaotic motion state, the maximum Lyapunov exponent is positive.

In general, it also shows the evolution law of nonlinear vibration characteristics with alternating periodic motion, quasi-periodic motion, and chaotic motion. The local stability of the system can be further analyzed based on the Lyapunov exponent.

Fig. 6LLE diagram of floating raft system with η

Fig. 7 shows the time history diagram, spectrum diagram, phase diagram and Poincaré map of the vertical vibration response of machinery A when mass ratio $\eta =$ 1.002. The time history diagram in Fig. 7(a) depicts the steady-state vibration response, and the vibration amplitude fluctuates periodically within a certain range. The spectrum diagram in 7(b) shows that there are obvious double frequency components in the system. The phase diagram in Fig. 7(c) plots the steady-state vibration response and the velocity numerical solution, and the phase trajectory forms a ring with a certain width. Fig. 7(d) Poincaré map with 2000 intervals $T$. The points of the curve are drawn, and their Poincaré mapping points form a closed curve. Based on the results in Fig. 7, machinery A is in a quasi-periodic motion state at this time.

Fig. 7Machinery A is in a quasi-periodic motion at η= 1.002

a) Time history

b) Frequency spectrum

c) Phase portrait

d) Poincaré map

Fig. 8 shows time history diagram, spectrum diagram, phase diagram and Poincaré map of the vertical vibration response of machinery A when the mass ratio $\eta =$ 1.050. Fig. 8(a) depicts the steady-state vibration response, and the vibration amplitude fluctuates periodically within a certain range. The spectrum diagram in 8(b) shows that there are complex frequency components in the system and continuous frequency. Fig. 8(c) plots the steady-state vibration response and the velocity numerical solution, and the phase trajectory is messy which indicates unstable state. Fig. 8(d) plots Poincaré map with 2000 intervals $T$. The points scattered in a certain area. Based on the results in Fig. 8, machinery A is currently exhibiting chaotic motion state.

Fig. 8Machinery A is in chaotic motion at η= 1.050

a) Time history

b) Frequency spectrum

c) Phase portrait

d) Poincaré map

5. Conclusions

The nonlinear dynamic behavior of a floating raft isolation system coupled with quasi-zero-stiffness isolators under multiple disturbance sources is carried out. In this study, based on the multi-rigid body dynamic theory, the mathematical model of the 9-degree-of-freedom system is developed. And its nonlinear dynamic characteristics are analyzed by applying the Runge-Kutta method. The dynamic responses are mainly investigated by using time history diagram, spectrum diagram, phase diagram, Poincaré diagram, and the bifurcation diagram. The results of this research significantly extend and enrich the understanding of the behavior of the system. Based on the results of the numerical simulations, the following conclusions can be summarized. First, the floating raft exhibits alternate phenomenon of periodic motion, quasi-periodic motion, and chaotic motion when the center distance and mass ratio vary. The motion state of the floating raft vibration isolation system is more sensitive to the mass ratio than to the center distance. Second, the horizontal and rotational responses of the system become very intense in the chaotic state, with the response amplitudes in the horizontal and vertical directions reaching the same order of magnitude. Above all, the dynamic characteristics can provide theoretic supporting for the dynamics, vibration control, and parametric optimization of the floating raft isolation system coupled with quasi-zero-stiffness isolators. This study is based on the basic assumption that the external excitation phase difference is $\pi$/2, which cannot fully reflect the real system working scenario. In the future, it is planned to carry out research on different phase differences.