**3. Analysis of Bifurcation Characteristic**

As below, the bifurcation characteristic of load vertical vibration is analyzed in detail. Moreover, the multi-scale method is utilized to solve the nonlinear equation of load vertical vibration.

Firstly, the small parameter factor ε and a series of slow-change time scales *Tn* are employed [49,50].

$$T\_n = \varepsilon^n t \quad n = 0, 1, 2 \cdots \tag{7}$$

Then, *Tn* are regarded as independent variables; there are:

$$\frac{\mathbf{d}}{\delta t} = \frac{\partial}{\partial T\_0} + \varepsilon \frac{\partial}{\partial T\_1} + \dots + \varepsilon^n \frac{\partial}{\partial T\_n} = D\_0 + \varepsilon D\_1 + \dots + \varepsilon^n D\_n \tag{8}$$

$$\begin{array}{rcl} \frac{\mathbf{d}^2}{\mathbf{d}t^2} &=& \frac{\mathbf{d}}{\mathbf{d}t} \Big( \frac{\partial}{\partial T\_0} + \varepsilon \frac{\partial}{\partial T\_1} + \dots + \varepsilon^n \frac{\partial}{\partial T\_n} \Big) \\ &=& \left( D\_0 + \varepsilon D\_1 + \dots + \varepsilon^n D\_n \right)^2 \\ &= D\_0^2 + 2\varepsilon D\_0 D\_1 + \varepsilon^2 (D\_1^2 + 2D\_0 D\_2) + \dotsb \end{array} \tag{9}$$

where *Dn* is a symbol of partial differential operator:

$$D\_n = \frac{\partial}{\partial T\_n} \quad n = 0, 1, 2 \cdots \tag{10}$$

Equation (1) can be further sorted as:

$$\ddot{y} + \alpha\_0^2 y = \frac{1}{m\_1} \left[ F \cos(\omega t) - \alpha F\_k(y) - \beta F\_c(\dot{y}) - c\_1 \dot{y} \right] \tag{11}$$

where <sup>ω</sup><sup>0</sup> <sup>=</sup> <sup>√</sup> *k*1/*m*1.

Then, the right side of Equation (11) is multiplied by ε; there is:

$$
\ddot{y} + \omega\_0^2 y = \frac{\varepsilon}{m\_1} \left[ F \cos(\omega t) - \alpha F\_k(y) - \beta F\_c(\dot{y}) - c\_1 \dot{y} \right] \tag{12}
$$

Suppose that ω is close to ω0, that is:

$$
\omega = a\wp + \varepsilon\sigma \tag{13}
$$

where ω is perturbation frequency, ω<sup>0</sup> is natural frequency, σ is the frequency tuning factor.

Equation (13) is substituted into Equation (12); there is:

$$
\ddot{y} + \omega\_0^2 y = \varepsilon \left[ \frac{F}{m\_1} \cos(\alpha y + \varepsilon \sigma) t + f(y, \dot{y}) \right] \tag{14}
$$

Among them:

$$f(y, \dot{y}) = -\frac{1}{m\_1} [aF\_k(y) + \beta F\_c(\dot{y}) + c\_1 \dot{y}] \tag{15}$$

Assume that Equation (14) has the solution with the following form when it is under the external excitation:

$$y = y\_0(T\_{0\prime}T\_1) + y\_1(T\_{0\prime}T\_1) \tag{16}$$

Substituting Equation (16) into Equation (14), Equations (8) and (9) are introduced. Then, the partial differential equations can be obtained:

$$\begin{cases} D\_0^2 y\_0 + \omega\_0^2 y\_0 = 0\\ D\_0^2 y\_1 + \omega\_0^2 y\_1 = -2D\_0 D\_1 y\_0 + f(y\_0, D\_0 y\_0) + \frac{F}{m\_1} \cos(\omega\_0 T\_0 + \sigma T\_1) \end{cases} \tag{17}$$

The solution for the first equation of Equation (17) is set as:

$$y\_0(T\_{0\prime}T\_1) = a(T\_1)\cos[\omega\_0 T\_0 + \psi(T\_1)] = A(T\_1)e^{i\omega\_0 T\_0} + \text{cc} \tag{18}$$

$$A(T\_1) = \frac{a(T\_1)}{2} \mathfrak{e}^{\mathrm{i}\psi(T\_1)}\tag{19}$$

where *a*(*T*1) and ψ(*T*1) are the slow-change functions of vibration amplitude and phase angle, respectively. cc indicates the conjugated plural.

Substituting Equation (18) into the second equation of Equation (17), there is:

$$\begin{aligned} D\_0^2 y\_1 + \omega\_0^2 y\_1 &= \mathrm{i}\omega\_0 \Big( -2D\_1 - \frac{\varepsilon\_1}{m\_1} \Big) \mathrm{Ae}^{i\omega\_0 T\_0} + \frac{F}{2m\_1} \mathrm{e}^{i(\omega\_0 T\_0 + \sigma T\_1)} \\ -\frac{a\varrho\_e}{m\_1} \Big[ \Big( \frac{A\_F}{L\_1} + \frac{A\_b}{L - L\_1} \Big) \mathrm{Ae}^{i\omega\_0 T\_0} + \Big( \frac{A\_F}{L\_1^3} + \frac{A\_b}{(L - L\_1)^3} \Big) \big( A^3 \mathrm{e}^{i3\omega\_0 T\_0} + 3A^2 \overline{A} \mathrm{e}^{i\omega\_0 T\_0} \Big) \Big] \\ -\frac{\beta F\_N}{m\_1} \Big[ \mu\_8 \mathrm{sgn} \left( \dot{y} \right) - \mathrm{i}\omega\_0 \mu\_1 A \mathrm{e}^{i\omega\_0 T\_0} - \mathrm{i}\omega\_0^3 \mu\_2 \big( A^3 \mathrm{e}^{i3\omega\_0 T\_0} - 3A^2 \overline{A} \mathrm{e}^{i\omega\_0 T\_0} \Big) \Big] + \infty \end{aligned} \tag{20}$$

*Processes* **2019**, *7*, 718

The secular term is further eliminated. In Equation (20), the coefficient of ei<sup>ω</sup>0*T*<sup>0</sup> is set to zero; there is:

$$\begin{split} \|\omega\_{0}(2D\_{1}+\frac{c\_{1}}{m\_{1}})A-\frac{F}{2m\_{1}}\mathbf{e}^{\mathbf{i}\alpha T\_{1}}+\frac{\alpha\theta\_{c}}{m\_{1}}\Big[\Big(\frac{A\_{p}}{L\_{1}}+\frac{A\_{b}}{L-L\_{1}}\Big)A+\Big(\frac{A\_{p}}{L\_{1}^{3}}+\frac{A\_{b}}{(L-L\_{1})^{3}}\Big)\big(3A^{2}\overline{A}\Big)\Big] \\ +\frac{\beta F\chi}{m\_{1}}\Big[-\mathrm{i}\alpha\_{0}\mu\_{1}A+\mathrm{i}\mathrm{i}3\alpha\_{0}^{3}\mu\_{2}A^{2}\overline{A}\Big]=0 \end{split} \tag{21}$$

Then, Equation (19) is substituted into Equation (21). Moreover, the real and imaginary parts are separated; there are:

$$\begin{cases} D\_1 a = -\frac{c\_1}{2m\_1} a + \frac{\beta F\_N}{m\_1} \left( \frac{\mu\_1}{L} a - \frac{\beta g^3}{8} \omega\_0^2 \mu\_2 \right) + \frac{F}{2m\_1 a \wp} \sin(\sigma T\_1 - \psi) \\\ D\_1 \psi = \frac{a \beta\_r}{2m\_1 a \wp} \left[ \left( \frac{A\_p}{L\_1} + \frac{A\_b}{L - L\_1} \right) + \frac{3 \varrho^2}{4} \left( \frac{A\_p}{L\_1} + \frac{A\_b}{(L - L\_1)^3} \right) \right] - \frac{F}{2m\_1 a \wp a} \cos(\sigma T\_1 - \psi) \end{cases} \tag{22}$$

Suppose that ϕ = σ*T*<sup>1</sup> − ψ, then Equation (22) can be changed into:

$$\begin{cases} D\_1 a = -\frac{\varepsilon\_1}{2m\_1} a + \frac{\beta F\_N}{m\_1} \left( \frac{\mu\_1}{2} a - \frac{9\mu^3}{8} \alpha\_0^2 \mu\_2 \right) + \frac{F}{2m\_1 \omega\_0} \sin \varphi\\ D\_1 \varphi = \sigma - \frac{a\mu\_r}{2m\_1 \omega\_0} \left[ \left( \frac{A\_p}{L\_1} + \frac{A\_b}{L - L\_1} \right) + \frac{3\mu^2}{4} \left( \frac{A\_p}{L\_1} + \frac{A\_b}{(L - L\_1)^3} \right) \right] + \frac{F}{2m\_1 \omega\_0 a} \cos \varphi \end{cases} \tag{23}$$

In Equation (23), when *D*1*a* = 0 and *D*1ϕ = 0, the system has a stable vibration amplitude and frequency; there are:

$$\begin{cases} \frac{c\_1}{m\_1}a - \frac{\beta \Gamma \wedge a}{m\_1} \left(\mu\_1 - \frac{g\_0^2}{4} \alpha\_0^2 \mu\_2\right) = \frac{F}{m\_1 \omega \alpha\_0} \sin \varphi\\ 2\sigma a - \frac{\alpha \beta \varepsilon a}{m\_1 \omega \alpha\_0} \left[ \left(\frac{A\_p}{L\_1} + \frac{A\_b}{L - L\_1}\right) + \frac{3a^2}{4} \left(\frac{A\_p}{L\_1^3} + \frac{A\_b}{\left(L - L\_1\right)^3}\right) \right] = -\frac{F}{m\_1 \omega \alpha\_0} \cos \varphi \end{cases} \tag{24}$$

In Equation (24), the square sum is performed to eliminate the parameter ϕ; there are:

$$\begin{aligned} \left\{ \frac{c\_1}{m\_1} a - \frac{\beta \mathcal{F}\_N a}{m\_1} \left( \mu\_1 - \frac{g\_0 \mathcal{F}}{4} \omega\_0^2 \mu\_2 \right) \right\}^2 + \\ \left\{ 2 \sigma a - \frac{a \beta\_0 a}{m\_1 \omega\_0} \left[ \left( \frac{A\_p}{L\_1} + \frac{A\_b}{L - L\_1} \right) + \frac{3 \sigma^2}{4} \left( \frac{A\_p}{L\_1^3} + \frac{A\_b}{\left( L - L\_1 \right)^3} \right) \right] \right\}^2 = \left( \frac{F}{m\_1 \omega\_0} \right)^2 \end{aligned} \tag{25}$$

Equation (25) is unfolded and organized; there is:

$$a^6 + \lambda\_1 a^4 + \lambda\_2 a^2 + \mu = 0\tag{26}$$

Among them:

$$\begin{split} \lambda\_{1} = \quad - \begin{bmatrix} \frac{\alpha \beta\_{r}}{3\omega \eta} \left[ \frac{A\_{p}}{L\_{1}^{3}} + \frac{A\_{b}}{\left(L - L\_{1}\right)^{3}} \right] \left[ m\_{1}\sigma - \frac{\alpha \beta\_{r}}{2\omega \eta} \left( \frac{A\_{p}}{L\_{1}} + \frac{A\_{b}}{L - L\_{1}} \right) \right] - \frac{\beta F\_{N}\mu\_{2}\omega\_{0}^{2}}{2} \left( c\_{1} - \beta F\_{N}\mu\_{1} \right) \end{bmatrix} \\ \quad + \begin{Bmatrix} \frac{\alpha^{2} \beta\_{r}^{2}}{16\omega\_{0}^{2}} \left[ \frac{A\_{p}}{L\_{1}^{3}} + \frac{A\_{b}}{\left(L - L\_{1}\right)^{3}} \right]^{2} + \frac{9\beta^{2} F\_{N}^{2}\mu\_{2}^{2}\omega\_{0}^{4}}{16} \end{Bmatrix} \tag{27}$$

$$\begin{split} \lambda\_{2} &= \quad \left\{ \left[ 2\sigma - \frac{\alpha \beta\_{r}}{m\_{1}\omega\wp} \left( \frac{A\_{p}}{L\_{1}} + \frac{A\_{b}}{L - L\_{1}} \right) \right]^{2} + \frac{\xi \mathbbm{1}}{m1} - \frac{\beta F\_{N}\mu\_{1}}{m1} \right)^{2} \right\} \\ &\overset{\perp}{\leftarrow} \left\{ \frac{9\alpha^{2} \beta\_{r}^{2}}{16m\_{1}^{2}\omega\_{0}^{2}} \left[ \frac{A\_{p}}{L\_{1}^{3}} + \frac{A\_{b}}{(L - L\_{1})^{3}} \right]^{2} + \frac{81\beta^{2} F\_{N}^{2} \mu\_{2}^{2} \omega\_{0}^{4}}{16m\_{1}^{2}} \right\} \end{split} \tag{28}$$

$$\mu = -\left(\frac{F}{m\_1}\right)^2 + \left\{\frac{9\alpha^2 \beta\_\varepsilon^2}{16m\_1^2} \left[\frac{A\_p}{L\_1^3} + \frac{A\_b}{\left(L - L\_1\right)^3}\right]^2 + \frac{81\beta^2 F\_N^2 \mu\_2^2 \omega\_0^6}{16m\_1^2}\right\} \tag{29}$$

Equation (26) can be further written as:

$$a^7 + \lambda\_1 a^5 + \lambda\_2 a^3 + \mu a = 0\tag{30}$$

According to the singularity theory [51], Equation (30) is the universal unfolding of the paradigm *a*<sup>7</sup> + μ*a* = 0, λ<sup>1</sup> and λ<sup>2</sup> are the open fold parameters, and μ is the external disturbance quantity. If the values of the open fold parameters are different, the system will have different bifurcation patterns.

Assume:

$$G(a,\mu,\lambda) = a^{\mathcal{T}} + \lambda\_1 a^5 + \lambda\_2 a^3 + \mu a \tag{31}$$

Then, the derivations about *a* and μ in Equation (31) are respectively performed, and the following can be obtained: .

$$\dot{G}\_a(a,\mu,\lambda) = 7a^6 + 5\lambda\_1 a^4 + 3\lambda\_2 a^2 + \mu \tag{32}$$

$$
\dot{G}\_{\mu}(a,\mu,\lambda) = a \tag{33}
$$

$$\ddot{\mathcal{G}}\_a(a,\mu,\lambda) = 42a^5 + 20\lambda\_1 a^3 + 6\lambda\_2 a \tag{34}$$

In the singularity theory, the transition set is a very important concept. This is the set of folding parameters corresponding to the non-persistent bifurcation diagram of the universal unfolding *G*(*a*, μ, λ). The sets of bifurcation point, lag point, and double limit point correspond to the three types of the non-persistence bifurcation diagram [52,53].

(1) Bifurcation point set

*<sup>B</sup>* <sup>=</sup> {<sup>λ</sup> <sup>∈</sup> **<sup>R</sup>***n*<sup>|</sup> exist (*a*, <sup>μ</sup>) <sup>∈</sup> **<sup>R</sup>** <sup>×</sup> **<sup>R</sup>**, make *<sup>G</sup>* <sup>=</sup> . *Ga* <sup>=</sup> . *G*<sup>μ</sup> = 0 at (*a*, μ, λ)}; (2) Lag point set *<sup>H</sup>* <sup>=</sup> {<sup>λ</sup> <sup>∈</sup> **<sup>R</sup>***n*<sup>|</sup> exist (*a*, <sup>μ</sup>) <sup>∈</sup> **<sup>R</sup>** <sup>×</sup> **<sup>R</sup>**, make *<sup>G</sup>* <sup>=</sup> . *Ga* <sup>=</sup> .. *Ga* = 0 at (*a*, μ, λ)}; (3) Double limit point set *<sup>D</sup>* <sup>=</sup> {<sup>λ</sup> <sup>∈</sup> **<sup>R</sup>***n*<sup>|</sup> exist (*ai*, <sup>μ</sup>)(*<sup>i</sup>* = 1, 2) <sup>∈</sup> **<sup>R</sup>** <sup>×</sup> **<sup>R</sup>**, *<sup>a</sup>*<sup>1</sup> *<sup>a</sup>*2, make *<sup>G</sup>* <sup>=</sup> . *Ga* = 0 at (*ai*, μ, λ)}; (4) Transition set: " = *B* ∪ *H* ∪ *D*

The transition set can divide the neighborhood of the origin of the real space **R***<sup>n</sup>* into several subregions, and the bifurcation diagram of the universal unfolding *G*(*a*, μ, λ) in each subregion is persistent. In the same subregion, the bifurcation diagrams corresponding to different folding parameters λ are equivalent. Then, all the persistent bifurcation maps of *G*(*a*, μ, λ) can be enumerated according to these subregions.

Next, singularity theory is applied to solve the transition set of bifurcation Equation (30), and the bifurcation characteristics of the system are explored. The subscripts *R* and *I* represent the new transition set generated by the conventional transition set and nonlinear action, respectively [54].

Bifurcation point set:

$$B\_R = \bigotimes \quad B\_I = \bigotimes \tag{35}$$

Lag point set:

$$H\_R = \left\{\frac{3\lambda\_2^2}{\lambda\_1^2} - \lambda\_2 = 0\right\} \quad H\_I = \mathcal{Q} \tag{36}$$

Double limit point set:

$$D\_R = \oslash \qquad D\_I = \left\{ 6a^4 + 4\lambda\_1 a^2 + 2\lambda\_2 = 0 \right\} \tag{37}$$

Transition set:

$$
\sum = B\_R \cup H\_R \cup D\_R \cup B\_I \cup H\_I \cup D\_I \tag{38}
$$

At this point, the transition set of the system under the open fold parameters λ<sup>1</sup> and λ<sup>2</sup> is shown in Figure 3. The topological structures of the bifurcation curves in different subregions divided by the transition set are displayed in Figure 4.

**Figure 4.** Topological structure of bifurcation curves in different subregions.

According to Figures 3 and 4, the system plane is divided into four subregions (I, II, III, IV) by the transition set under the open fold parameters λ<sup>1</sup> and λ2. The bifurcation curves in each subregion have their own topological structure. Furthermore, with the change of λ<sup>1</sup> and λ2, the topological structure changes at the transition set. It indicates that the system has different vibration behaviors in diverse subregions and will exhibit different bifurcation behaviors under various parameter combinations. Therefore, by analyzing the bifurcation characteristics of the system, the parameter region that causes the system to be unstable can be determined. Meanwhile, from the obtained topological structures of the bifurcation curve in different subregions, the bifurcation state of the system can be changed by the disturbance parameter μ. In a certain bifurcation form, the change of μ will lead to the change of vibration amplitude, which in turn changes the vibration behavior of the system.
