Periodic and torus motions of a two-degree-of-freedom dry friction vibration system

Vibration induced by dry friction is ubiquitous in various engineering fields. To explore the vibration characteristics for further studies and/or controls, it is of great theoretical and practical significances to investigate the non-linear dynamic behaviors of the friction systems. This study considers the slight vibration of a two-degree-of-freedom non-linear dry friction excitation system. The differential equations of system motion are established according to Newton’s law of motion. Moreover, the system’s non-linear dynamic is studied when the block velocity is always less than the friction surface velocity. The results indicate that the linearized matrix of the vibration system has a pair of purely imaginary eigenvalues for some critical values of the relevant parameters. The Poincaré-Birkhoff normal forms are utilized to simplify the motion equation under the non-resonant assumption to obtain a simplified equation with only the resonant terms. Furthermore, the truncated part of the simplified equation is analyzed in the case of only linear terms degeneration. Finally, numerical simulations reflect some qualitative conclusions about the system’s local dynamic properties, including equilibrium point, periodic motion, torus motion, and their stability.

1 Introduction

There are numerous dry friction phenomena in engineering practice, including wheel-rail contact, locomotive braking systems, machine tool guides, tool cutting, drilling, and friction damping. Indeed, dry friction can cause surface wear and fatigue failure of mechanical components, and the resulting dynamic behavior (e.g., generating noises) harms production and living environments. The phenomenon of machine tools, drill pipe chatter, brake whistle, and other phenomena caused by dry friction excitation will cause severe problems in engineering. However, dry friction can dissipate system energy to achieve vibration reduction or stimulate the required motion state. Therefore, studying the motion law of objects under the dry friction action is particularly important for finding effective engineering control methods and effectively utilizing the dry friction phenomenon.

Leine et al. [1] investigated the existence and stability of the periodic solutions of the vector field in the case of continuous, non-smooth, and unique solutions, corresponding to the existence and stability of the periodic solution of the dry friction system in the pure slip state. Luo [2] established the Poincaré mapping of the dry friction system in the pure slip state and analyzed the system’s periodic motion and the mapping bifurcation. A specific periodic solution was obtained in the presence of adhesion by adjusting the parameters to satisfy the functional relationship on the vector field interface. The dry friction system’s periodic motion and chaos were studied by numerical simulations [3–5]. Guo et al. [6] established the Poincaré mapping of the single-degree-of-freedom dry friction vibrator based on the series solution of the differential equation and investigated the existence and stability of the system under 1:4 strength resonance for periodic and torus motions. Guo [7] investigated the existence and stability of subharmonic periodic and torus motions of single-degree-of-freedom dry friction vibrator under 1:5 weak resonance. When the amplitude of the vibration induced by the dry friction reached to some motion constraints, collisions will occur, which may result in complex non-linear dynamic phenomena, such as chaos and chaos control [8]. Refs. [9–14] present the researches of the collision system or the friction collision coexistence system. By applying the stroboscopic controlled hybrid Poincaré map and OGY control method, Gritli and Belghith [9] investigated the non-linear dynamic characteristic and chaos control of a one-degree-of-freedom impact oscillator with a single rigid constraint, in which the border-collision bifurcation is explored for the use of OGY-based state-feedback control method. A two-degree-of-freedom impact oscillator with dry friction and external periodically forced excitation is considered by Li et al. [10] based on the flow switchability theory and G-functions. A one-degree-of-freedom flexible-impact oscillator is investigated by Stefani et al. [11], where some possible dynamic scenarios were obtained by experimental tests and numerical simulations. In another paper by Stefani et al. [13], a comprehensive numerical study is conducted, where, compared to Ref. [11], the range of selected parameters were extended. Based on the numerical results in Ref. [13], Stefani et al. [12] investigated further the effect of gap size on the response of a one-degree-of-freedom flexible-impact oscillator, in which the secondary resonances were observed for quite small gap and the number of resonances were analyzed by changing the gap size. Peng and Fan [14] studied a three-degree-of-freedom rigid-impact oscillator with dry friction by using the flow switchability theory, and divided the six-dimensional phase space of the system into different domains and boundaries/edges. It can be seen from the above mentioned researches that numerical simulation and experimental method are the main means to explore the non-linear dynamic phenomena of system with friction and/or collision. Furthermore, in the previous studies, the bifurcation analysis after impact generally began with some periodic solution. To the best of the author’s knowledge, when the amplitude of torus motion reached to motion constraints and then collision occurs, the non-linear vibration characteristic (i.e., the bifurcation analysis after impact begin with some torus motion) is an open problem, which motivated the current study. As an initial investigation, only the periodic and torus motions within the smooth case are considered in this paper, which provides a theoretical basis for further analysis collision with motion constraints.

Non-linear vibrations are widely encountered in engineering practice. As one of the important characteristics of non-linear vibrations, the periodic and/or torus motions have attracted much attention [15–17]. Hopf-Hopf bifurcation is an approach to investigate the periodic and torus motions of dynamical system. Wen et al. [18] proposed a new Hopf-Hopf bifurcation criterion based on the coefficients of the original equation of dynamics for a dry friction system and the corresponding feedback control method was investigated. Guo et al. [19] investigated Hopf-Hopf bifurcation for a simplified railway wheelset model in the case of non-resonance and near-resonance, in which the resonant coefficients and the truncated Poincaré-Birkhoff normal form were computed by use of MATCONT [20]. A numerical analysis of Hopf-Hopf bifurcation for a non-linear electric oscillator in the case of 1:2 resonance was performed by Revel et al. [21], where the resonance could arise more non-trivial mode interactions and lead to more complex dynamical phenomena. Hopf-Hopf bifurcation can also be found in other dynamical systems, for instance in vibro-impact system [22], high-dimensional maps [23] and infinite dynamical system [24]. It is noted that in the above literature of Refs. [18, 19, 21, 24] concerning differential equations, the non-linear resonant terms either were computed numerically in which the analytical expressions cannot be given or were simplified by use of symmetry which cannot apply to asymmetric systems. For this reason, the expressions of non-linear resonant terms for a system without symmetry will be given analytically in the current paper.

By means of the Hopf-Hopf bifurcation theory, the periodic and torus motions of a two-degree-of-freedom dynamical system with dry friction are considered both analytically and numerically in this paper. This paper is organized as follows: In Section 2, the differential equation of the system motion is established according to Newton’s laws of motion, and the problem is confined to the “smooth motion” case. This means the block velocity is considered smaller than the friction surface velocity. Section 3 employs the Poincaré-Birkhoff normal forms to process the dynamic equation in Section 2 to obtain a simplified equation equivalent to the original dynamic equation, which only contains the resonant terms. Moreover, the calculation formulas of the resonant terms coefficients and the eigenvalues’ rate of change on the perturbated parameters are given analytically. In Section 4, the amplitude equation is studied with specific examples, and the conclusions of the existence and stability of the equilibrium point and periodic and torus motions of the original vibration system are investigated. Some conclusions and discussions are made in Section 5.

2 Mechanical model and motion equations

Figure 1 shows the schematic of a two-degree-of-freedom dry friction system [25].

FIGURE 1

FIGURE 1. The schematic of a two-degree-of-freedom dry friction system.

The mass of each mass block is m, placed on two rotating circles. The speed of the circle’s edge is v. The remote ends of the two blocks are connected to the fixed end by springs with stiffness coefficients k₁ and k₂, respectively, and a spring-damping element is employed to connect the adjacent ends. The spring stiffness coefficient is k, and the damping coefficient is c. According to Newton’s law of motion, taking z₁ and z₂ as generalized coordinates, the following motion’s differential equation can be obtained:

$\{\begin{array}{l}m{\ddot{z}}_1=-k_1{z}_1-k\left(z_1+z_2\right)-c\left({\dot{z}}_1+{\dot{z}}_2\right)+F\left({\dot{z}}_1-v\right)\\ m{\ddot{z}}_2=-k_2{z}_2-k\left(z_1+z_2\right)-c\left({\dot{z}}_1+{\dot{z}}_2\right)+F\left({\dot{z}}_2-v\right)\end{array}(1)$

where $F\left({\dot{z}}_1-v\right)$ and $F\left({\dot{z}}_2-v\right)$ are dry frictions defined with the following general expression [1]:

$F\left(v_{\mathrm{r}}\right)=\left[-{\alpha }_0\mathrm{sign}\left(v_{\mathrm{r}}\right)+{\alpha }_1{v}_{\mathrm{r}}-{\alpha }_3{v}_{\mathrm{r}}^3\right]\cdot m(2)$

where v_r is the relative velocity of the mass block relative to the edge of the wheel, and ${\alpha }_0=1.5\mathrm{N}$, ${\alpha }_1=1.5\text{Ns}/\mathrm{m}$, ${\alpha }_3=0.45{\text{Ns}}^3{/\mathrm{m}}^3$. As shown by Figure 2.

FIGURE 2

FIGURE 2. The characteristic curve of the dry friction function described by Eq. 2.

This study considers that the relative velocity between the mass block and the edge of the wheel is always negative. Equivalently, a smooth vibration is considered. Thus, Eq. 2 can be written as:

$F\left(v_{\mathrm{r}}\right)=\left({\alpha }_0+{\alpha }_1{v}_{\mathrm{r}}-{\alpha }_3{v}_{\mathrm{r}}^3\right)\cdot m(3)$

Let ${\ddot{z}}_1={\ddot{z}}_2={\dot{z}}_1={\dot{z}}_2=0$, the system’s equilibrium position can be obtained as:

$z_1^0=\frac{F\left(-v\right)}{k_1+k\left(1+\frac{k_1}{k_2}\right)},z_2^0=\frac{F\left(-v\right)}{k_2+k\left(1+\frac{k_2}{k_1}\right)}$

To perform coordinate translation, let

$z_1=u_1+z_1^0,z_2=u_2+z_2^0(4)$

Hence, the origin is the system’s equilibrium position. Substituting Eq. 4 into Eq. 1 gives:

$\left\{\begin{array}{l}m{\ddot{u}}_1=-k_1{u}_1-k\left(u_1+u_2\right)-c\left({\dot{u}}_1+{\dot{u}}_2\right)+F\left({\dot{u}}_1-v\right)-F\left(-v\right)\\ m{\ddot{u}}_2=-k_2{u}_2-k\left(u_1+u_2\right)-c\left({\dot{u}}_1+{\dot{u}}_2\right)+F\left({\dot{u}}_2-v\right)-F\left(-v\right)\end{array}\right.(5)$

Let $u_1=y_1$, ${\dot{u}}_1=y_2$ , $u_2=y_3$ , ${\dot{u}}_2=y_4$, and expand the friction function according to Taylor’s formula. Equation 5 can be rewritten in the following matrix form

$\left[\begin{array}{c}{\dot{y}}_1\\ {\dot{y}}_2\\ {\dot{y}}_3\\ {\dot{y}}_4\end{array}\right]=\mathbf{A}\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]+\left[\begin{array}{c}0\\ \frac{1}{m}\left[\frac{1}{2}{F}^{″}\left(-v\right)y_2^2+\frac{1}{6}{F}^{‴}\left(-v\right)y_2^3\right]\\ 0\\ \frac{1}{m}\left[\frac{1}{2}{F}^{″}\left(-v\right)y_4^2+\frac{1}{6}{F}^{‴}\left(-v\right)y_4^3\right]\end{array}\right](6)$

where

$\mathbf{A}=\left[\begin{array}{cccc}0& 1& 0& 0\\ -\frac{k+k_1}{m}& \frac{F^{\prime }\left(-v\right)-c}{m}& -\frac{k}{m}& -\frac{c}{m}\\ 0& 0& 0& 1\\ -\frac{k}{m}& -\frac{c}{m}& -\frac{k+k_2}{m}& \frac{F^{\prime }\left(-v\right)-c}{m}\end{array}\right](7)$

let

$\mathbf{y}={\left[y_1,y_2,y_3,y_4\right]}^{\mathrm{T}}$

$\mathbf{F}\left(\mathbf{y}\right)=\left[\begin{array}{c}0\\ \frac{1}{m}\left[\frac{1}{2}{F}^{″}\left(-v\right)y_2^2+\frac{1}{6}{F}^{‴}\left(-v\right)y_2^3\right]\\ 0\\ \frac{1}{m}\left[\frac{1}{2}{F}^{″}\left(-v\right)y_4^2+\frac{1}{6}{F}^{‴}\left(-v\right)y_4^3\right]\end{array}\right]$

To be sure, the highest order terms of Eq. 3 and Eq. 5 are cubic. Therefore, Taylor’s formula (6) for Eq. 5 are actually cubic polynomial rather than infinite series.

3 Normal form theory

Generally, the following differential equations are considered [26]:

$\left\{\begin{array}{l}{\dot{z}}_1={\lambda }_1\left(\mathbf{\alpha }\right)z_1+g\left(z_1,{\overline{z}}_1,z_2,{\overline{z}}_2,\mathbf{\alpha }\right)\\ {\dot{z}}_2={\lambda }_2\left(\mathbf{\alpha }\right)z_2+h\left(z_1,{\overline{z}}_1,z_2,{\overline{z}}_2,\mathbf{\alpha }\right)\end{array}\right.(8)$

where $z_1,z_2,{\lambda }_1,{\lambda }_2$ are plurals, the parameter $\mathbf{\alpha }=\left({\alpha }_1,{\alpha }_2\right)$ is a two-dimensional real vector, and

${\lambda }_1={\mu }_1\left(\mathbf{\alpha }\right)+\mathrm{i}{\omega }_1\left(\mathbf{\alpha }\right),{\lambda }_2={\mu }_2\left(\mathbf{\alpha }\right)+\mathrm{i}{\omega }_2\left(\mathbf{\alpha }\right).$

where $\mathrm{i}$ is an imaginary unit. For all sufficiently small $‖\mathbf{\alpha }‖$, ${\mu }_1,{\mu }_2$, and ${\omega }_1,{\omega }_2$ are smooth functions of $\mathbf{\alpha }$ , meeting the following relations

${\mu }_1\left(\mathbf{0}\right)={\mu }_2\left(\mathbf{0}\right)=0,{\omega }_1\left(\mathbf{0}\right)>{\omega }_2\left(\mathbf{0}\right)>0$

Let $\mathbf{\mu }=\left({\mu }_1,{\mu }_2\right)=\left({\mu }_1\left(\mathbf{\alpha }\right),{\mu }_2\left(\mathbf{\alpha }\right)\right)$, ${\omega }_{10}={\omega }_1\left(\mathbf{0}\right),{\omega }_{20}={\omega }_2\left(\mathbf{0}\right)$.

The non-linear terms in Eq. 8 can be represented with the following Taylor expansions with respect to the first four arguments:

$g\left(z_1,{\overline{z}}_1,z_2,{\overline{z}}_2,\mathbf{\alpha }\right)={\displaystyle \sum }_{j+k+l+m\ge 2}{g}_{jklm}\left(\mathbf{\alpha }\right)z_1^j{\overline{z}}_1^k{z}_1^l{\overline{z}}_1^m$

$h\left(z_1,{\overline{z}}_1,z_2,{\overline{z}}_2,\mathbf{\alpha }\right)={\displaystyle \sum }_{j+k+l+m\ge 2}{h}_{jklm}\left(\mathbf{\alpha }\right)z_1^j{\overline{z}}_1^k{z}_1^l{\overline{z}}_1^m$

Lemma 1 [26]: For the above system, if.

$k{\omega }_{10}\ne l{\omega }_{20}\left(k,l>0,k+1\le 5\right),$

near the equilibrium point, Eq. 8 can be simplified to the following form through the reversible change of coordinates:

$\left\{\begin{array}{l}{\dot{w}}_1={\lambda }_1\left(\mathbf{\alpha }\right)w_1+G_{2100}\left(\mathbf{\alpha }\right)w_1{\left|w_1\right|}^2+G_{1011}\left(\mathbf{\alpha }\right)w_1{\left|w_2\right|}^2\\ +G_{3200}\left(\mathbf{\alpha }\right)w_1{\left|w_1\right|}^4+G_{2111}\left(\mathbf{\alpha }\right)w_1{\left|w_1\right|}^2{\left|w_2\right|}^2\\ +G_{1022}\left(\mathbf{\alpha }\right)w_1{\left|w_2\right|}^4+O\left({‖w_1,{\overline{w}}_1,w_2,{\overline{w}}_2‖}^6\right)\\ {\dot{w}}_2={\lambda }_2\left(\mathbf{\alpha }\right)w_2+H_{1110}\left(\mathbf{\alpha }\right)w_2{\left|w_1\right|}^2+H_{0021}\left(\mathbf{\alpha }\right)w_2{\left|w_2\right|}^2\\ +H_{2210}\left(\mathbf{\alpha }\right)w_2{\left|w_1\right|}^4+H_{1121}\left(\mathbf{\alpha }\right)w_2{\left|w_1\right|}^2{\left|w_2\right|}^2\\ +H_{0032}\left(\mathbf{\alpha }\right)w_2{\left|w_2\right|}^4+O\left({‖w_1,{\overline{w}}_1,w_2,{\overline{w}}_2‖}^6\right)\end{array}\right.(9)$

where $w_{\mathrm{1,2}}\in {\mathrm{C}}^1$ , ${‖\left(w_1,{\overline{w}}_1,w_2,{\overline{w}}_2\right)‖}^2=\left|w_1^2\right|+\left|w_2^2\right|$, and the complex-valued functions $G_{jklm}\left(\mathbf{\alpha }\right)$ and $H_{jklm}\left(\mathbf{\alpha }\right)$ are smooth and satisfy the following relation:

$\begin{array}{l}{G}_{2100}\left(\mathbf{0}\right)=\\ g_{2100}+\frac{\mathrm{i}}{{\omega }_{10}}{g}_{1100}{g}_{2000}+\frac{\mathrm{i}}{{\omega }_{20}}\left(g_{1010}{h}_{1100}-g_{1001}{\overline{h}}_{1100}\right)\\ -\frac{\mathrm{i}}{2{\omega }_{10}+{\omega }_{20}}{g}_{0101}{\overline{h}}_{0200}-\frac{\mathrm{i}}{2{\omega }_{10}-{\omega }_{20}}{g}_{0110}{h}_{2000}\\ -\frac{\mathrm{i}}{{\omega }_{10}}{\left|g_{1100}\right|}^2-\frac{2\mathrm{i}}{3{\omega }_{10}}{\left|g_{0200}\right|}^2\end{array}(10)$

$\begin{array}{l}{G}_{1011}\left(\mathbf{0}\right)=\\ g_{1011}+\frac{\mathrm{i}}{{\omega }_{20}}\left(g_{1010}{h}_{0011}-g_{1001}{\overline{h}}_{0011}\right)+\frac{\mathrm{i}}{{\omega }_{10}}(2g_{2000}{g}_{0011}-\\ g_{1100}{\overline{g}}_{0011}-g_{0011}{\overline{h}}_{0110}-g_{0011}{h}_{1010})-\frac{2\mathrm{i}}{{\omega }_{10}+2{\omega }_{20}}{g}_{0002}{\overline{h}}_{0101}-\\ \frac{2\mathrm{i}}{{\omega }_{10}-2{\omega }_{20}}{g}_{0020}{h}_{1001}-\frac{\mathrm{i}}{2{\omega }_{10}-{\omega }_{20}}{\left|g_{0110}\right|}^2-\frac{\mathrm{i}}{2{\omega }_{10}+{\omega }_{20}}{\left|g_{0101}\right|}^2\end{array}(11)$

$\begin{array}{l}{H}_{1110}\left(\mathbf{0}\right)=\\ h_{1110}+\frac{\mathrm{i}}{{\omega }_{10}}\left(g_{1100}{h}_{1010}-{\overline{g}}_{1100}{h}_{0110}\right)+\frac{\mathrm{i}}{{\omega }_{20}}(2h_{0020}{h}_{1100}-\\ h_{0011}{\overline{h}}_{1100}-g_{1010}{h}_{1100}-{\overline{g}}_{1001}{h}_{1100})+\frac{2\mathrm{i}}{2{\omega }_{10}-{\omega }_{20}}{g}_{0110}{h}_{2000}-\\ \frac{2\mathrm{i}}{2{\omega }_{10}+{\omega }_{20}}{\overline{g}}_{0101}{h}_{0200}-\frac{\mathrm{i}}{2{\omega }_{20}-{\omega }_{10}}{\left|h_{1001}\right|}^2-\frac{\mathrm{i}}{{\omega }_{10}+2{\omega }_{20}}{\left|h_{0101}\right|}^2\end{array}(12)$

$\begin{array}{l}{H}_{0021}\left(\mathbf{0}\right)=\\ h_{0021}+\frac{\mathrm{i}}{{\omega }_{10}}\left(g_{0011}{h}_{1010}-{\overline{g}}_{0011}{h}_{0110}\right)+\\ \frac{\mathrm{i}}{{\omega }_{20}}{h}_{0011}{h}_{0020}-\frac{\mathrm{i}}{2{\omega }_{20}-{\omega }_{10}}{g}_{0020}{h}_{1001}-\\ \frac{\mathrm{i}}{2{\omega }_{20}+{\omega }_{10}}{\overline{g}}_{0002}{h}_{0101}-\frac{\mathrm{i}}{{\omega }_{20}}{\left|h_{0011}\right|}^2-\frac{2\mathrm{i}}{3{\omega }_{20}}{\left|h_{0002}\right|}^2\end{array}(13)$

In relations (10–13), all $g_{jklm},h_{jklm}$ are taken value at $\mathbf{\alpha }=\mathbf{0}$.

Lemma 2 [26]: Based on Lemma 1, if the following relations are satisfied.

$\mathrm{R}\mathrm{e}{\mathit{G}}_{2100}\mathbf{(}\mathbf{0})\ne 0;\mathrm{R}\mathrm{e}{\mathit{G}}_{1011}(\mathbf{0})\ne 0;$

$\mathrm{R}\mathrm{e}{\mathit{H}}_{1110}(\mathbf{0})\ne 0;\mathrm{R}\mathrm{e}{\mathit{H}}_{0021}(\mathbf{0})\ne 0;$

$\mathrm{d}\mathrm{e}\mathrm{t}(\frac{\mathrm{\partial }\mathbf{\mu }}{\mathbf{\partial }\mathbf{\alpha }}){|}_{\mathbf{\alpha }=\mathbf{0}}\ne 0$

Then, Eq. 9 can be transformed into the following form through coordinate transformation and modular, complex angle substitution:

$\left\{\begin{array}{l}{\dot{r}}_1=r_1\left({\mu }_1+p_{11}\left(\mathbf{\mu }\right)r_1^2+p_{12}\left(\mathbf{\mu }\right)r_2^2+s_1\left(\mathbf{\mu }\right)r_2^4\right)+{\mathrm{\Phi }}_1\left(r_1,r_2,{\phi }_1,{\phi }_2,\mathbf{\mu }\right)\\ {\dot{r}}_2=r_2\left({\mu }_2+p_{21}\left(\mathbf{\mu }\right)r_1^2+p_{22}\left(\mathbf{\mu }\right)r_2^2+s_2\left(\mathbf{\mu }\right)r_1^4\right)+{\mathrm{\Phi }}_2\left(r_1,r_2,{\phi }_1,{\phi }_2,\mathbf{\mu }\right)\\ {\dot{\varphi }}_1={\omega }_1\left(\mathbf{\mu }\right)+{\mathrm{\Psi }}_1\left(r_1,r_2,{\phi }_1,{\phi }_2,\mathbf{\mu }\right)\\ {\dot{\varphi }}_2={\omega }_2\left(\mathbf{\mu }\right)+{\mathrm{\Psi }}_2\left(r_1,r_2,{\phi }_1,{\phi }_2,\mathbf{\mu }\right)\end{array}\right.(14)$

where $r_1,r_2$ are modulus variables, and ${\phi }_1,{\phi }_2$ are argument variables. The relevant coefficient satisfies the following relationship:$p_{11}\left(\mathbf{0}\right)=\mathrm{R}\mathrm{e}{G}_{2100}\left(\mathbf{0}\right),p_{12}\left(\mathbf{0}\right)=\mathrm{R}\mathrm{e}{G}_{1011}\left(\mathbf{0}\right)$$p_{21}\left(\mathbf{0}\right)=\mathrm{R}\mathrm{e}{H}_{1110}\left(\mathbf{0}\right)$, $p_{22}\left(\mathbf{0}\right)=\mathrm{R}\mathrm{e}{H}_{0021}\left(\mathbf{0}\right)$; $s_1,s_2$ are real numbers.

In the case of simple degeneration, considering the truncated part of the modular equation in Eq. 14, we have

$\left\{\begin{array}{l}{\dot{r}}_1=r_1\left({\mu }_1+p_{11}{r}_1^2+p_{12}{r}_2^2+s_1{r}_2^4\right)\\ {\dot{r}}_2=r_2\left({\mu }_2+p_{21}{r}_1^2+p_{22}{r}_2^2+s_2{r}_1^4\right)\end{array}\right.(15)$

Accordingly, the primary motion forms of the system near the equilibrium position and their mutual transformations with parameter changes can be obtained. The orbit structure on a torus in Eq. 14 is generically different from that in the truncated Equation 15 due to phase locking. Nevertheless, for the non-resonance case [i.e., kω₁₀≠lω₂₀, (k, l > 0, k + l ≤ 5)] in this paper, the qualitative dynamical characteristics for Eq. 14 are the same as Eq. 15 from the observability point of view ^{[26, p. 368]}.

In the proposed model, it can be seen that its linearized matrix $\mathbf{A}$ has two pairs of purely imaginary eigenvalues if $F^{\prime }\left(-v\right)$ and $c$ are zero simultaneously. Let $v_c$ be the critical velocity of $F^{\prime }\left(-v\right)=0$ and $\mathbf{\alpha }=\left({\alpha }_1,{\alpha }_2\right)=\left(v-v_c,c\right)$ as the perturbed parameter. Then, matrix $\mathbf{A}$ is a function of ${\alpha }_1,{\alpha }_2$, which can be written as:

$\mathbf{A}\left(\mathbf{\alpha }\right)=\left[\begin{array}{cccc}0& 1& 0& 0\\ -\frac{k+k_1}{m}& \frac{F^{\prime }\left(-v_c-{\alpha }_1\right)-{\alpha }_2}{m}& -\frac{k}{m}& -\frac{{\alpha }_2}{m}\\ 0& 0& 0& 1\\ -\frac{k}{m}& -\frac{{\alpha }_2}{m}& -\frac{k+k_2}{m}& \frac{F^{\prime }\left(-v_c-{\alpha }_1\right)-{\alpha }_2}{m}\end{array}\right](16)$

let $\frac{k_1}{m}={\omega }_1^2$, $\frac{k_2}{m}={\omega }_2^2$ , $\frac{k}{m}={\omega }^2$.

$\begin{array}{l}\mathbf{A}\left(\mathbf{0}\right)=\left[\begin{array}{cccc}0& 1& 0& 0\\ -\frac{k+k_1}{m}& 0& -\frac{k}{m}& 0\\ 0& 0& 0& 1\\ -\frac{k}{m}& 0& -\frac{k+k_2}{m}& 0\end{array}\right]\\ =\left[\begin{array}{cccc}0& 1& 0& 0\\ -\left({\omega }^2+{\omega }_1^2\right)& 0& -{\omega }^2& 0\\ 0& 0& 0& 1\\ -{\omega }^2& 0& -\left({\omega }^2+{\omega }_2^2\right)& 0\end{array}\right]\end{array}(17)$

$\mathbf{A}\left(\mathbf{0}\right)$ has two pairs of purely imaginary eigenvalues: $\pm \mathrm{i}{\omega }_{10}$, $\pm \mathrm{i}{\omega }_{20}$ (${\omega }_{10}>{\omega }_{20}>0$), assuming that it satisfies the non-resonance condition, namely:

$k{\omega }_{10}\ne l{\omega }_{20}\left(k,l>0,k+l\le 5\right).$

When ${\alpha }_1,{\alpha }_2$ are slightly perturbed, the eigenvalues and the eigenvectors of $\mathbf{A}\left(\mathbf{\alpha }\right)$ depend on ${\alpha }_1,{\alpha }_2$ smoothly. Let ${\lambda }_1\left(\mathbf{\alpha }\right)={\mu }_1\left(\mathbf{\alpha }\right)+\mathrm{i}{\omega }_1\left(\mathbf{\alpha }\right)$, ${\lambda }_2\left(\mathbf{\alpha }\right)={\mu }_2\left(\mathbf{\alpha }\right)+\mathrm{i}{\omega }_2\left(\mathbf{\alpha }\right)$ be the eigenvalues of $\mathbf{A}\left(\mathbf{\alpha }\right)$, satisfying ${\mu }_1\left(\mathbf{0}\right)={\mu }_2\left(\mathbf{0}\right)=\mathbf{0},{\omega }_1\left(\mathbf{0}\right)={\omega }_{10},{\omega }_2\left(\mathbf{0}\right)={\omega }_{20}$; And ${\mathbf{q}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_2\left(\mathbf{\alpha }\right)\left(\in C^4\right)$ are its corresponding eigenvectors, namely,

$\left\{\begin{array}{l}\mathbf{A}\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)={\lambda }_1\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)\\ \mathbf{A}\left(\mathbf{\alpha }\right){\mathbf{q}}_2\left(\mathbf{\alpha }\right)={\lambda }_2\left(\mathbf{\alpha }\right){\mathbf{q}}_2\left(\mathbf{\alpha }\right)\end{array}\right.(18)$

Consider the following transformation

$\mathbf{y}=\mathbf{B}\left(\mathbf{\alpha }\right)\mathbf{x}(19)$

where $\mathbf{B}\left(\mathbf{\alpha }\right)=\left[\mathrm{R}\mathrm{e}{\mathbf{q}}_1\left(\mathbf{\alpha }\right),-\mathrm{I}\mathrm{m}{\mathbf{q}}_1\left(\mathbf{\alpha }\right),\mathrm{R}\mathrm{e}{\mathbf{q}}_2\left(\mathbf{\alpha }\right),-\mathrm{I}\mathrm{m}{\mathbf{q}}_2\left(\mathbf{\alpha }\right)\right]$ (20)

By substituting Eq. 19 into Eq. 6, we have:

$\dot{\mathbf{x}}=\mathbf{T}\left(\mathbf{\alpha }\right)\mathbf{x}+\stackrel{\sim }{\mathbf{F}}\left(\mathbf{x},\mathbf{\alpha }\right)(21)$

$\mathbf{T}\left(\mathbf{\alpha }\right)={\mathbf{B}}^{-1}\left(\mathbf{\alpha }\right)\mathbf{A}\left(\mathbf{\alpha }\right)\mathbf{B}\left(\mathbf{\alpha }\right)=\left[\begin{array}{cccc}{\mu }_1\left(\mathbf{\alpha }\right)& -{\omega }_1\left(\mathbf{\alpha }\right)& & \\ {\omega }_1\left(\mathbf{\alpha }\right)& {\mu }_1\left(\mathbf{\alpha }\right)& & \\ & & {\mu }_2\left(\mathbf{\alpha }\right)& -{\omega }_2\left(\mathbf{\alpha }\right)\\ & & {\omega }_2\left(\mathbf{\alpha }\right)& {\mu }_2\left(\mathbf{\alpha }\right)\end{array}\right](22)$

$\stackrel{\sim }{\mathbf{F}}\left(\mathbf{x},\mathbf{\alpha }\right)={\mathbf{B}}^{-1}\left(\mathbf{\alpha }\right)\mathbf{F}\left(\mathbf{B}\left(\mathbf{\alpha }\right)\mathbf{x}\right)(23)$

Let $z_1=x_1+x_2\mathrm{i}$, $z_2=x_3+x_4\mathrm{i}$ , the following equation can be obtained:

$x_1=\frac{z_1+{\overline{z}}_1}{2},x_2=\frac{z_1-{\overline{z}}_1}{2\mathrm{i}},x_3=\frac{z_2+{\overline{z}}_2}{2},x_4=\frac{z_2-{\overline{z}}_2}{2\mathrm{i}}$

By substituting the above equations into Eq. 21, the following plural form similar to Eq. 8 can be obtained:

$\left\{\begin{array}{l}{\dot{z}}_1={\lambda }_1\left(\mathbf{\alpha }\right)z_1+F_1\left(z_1,{\overline{z}}_1,z_2,{\overline{z}}_2,\mathbf{\alpha }\right)\\ {\dot{z}}_2={\lambda }_2\left(\mathbf{\alpha }\right)z_2+F_2\left(z_1,{\overline{z}}_1,z_2,{\overline{z}}_2,\mathbf{\alpha }\right)\end{array}\right.(24)$

Since only the case of linear terms degeneration is considered, it is only necessary to consider the first-order approximate expression of the eigenvalues’ real part $\mathbf{\mu }$ of matrix $\mathbf{A}\left(\mathbf{\alpha }\right)$ near ${\alpha }_1={\alpha }_2=0$, and to calculate the third-order resonant terms’ coefficients of normal form at ${\alpha }_1={\alpha }_2=0$.

First, the coefficients of the cubic resonant terms are calculated. In order to obtain the resonant terms’ coefficients in normal form, take ${\alpha }_1={\alpha }_2=0$ in Eq. 24, and substitute its coefficients into Equations 10-13.

Let

${\mathbf{q}}_1\left(\mathbf{0}\right)=\left[\begin{array}{c}{q}_1^1\\ q_1^2\\ q_1^3\\ q_1^4\end{array}\right]=\left[\begin{array}{c}\mathrm{i}\\ -{\omega }_{10}\\ \frac{\mathrm{i}\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right)}{{\omega }^2}\\ \frac{-{\omega }_{10}\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right)}{{\omega }^2}\end{array}\right](25)$

${\mathbf{q}}_2\left(\mathbf{0}\right)=\left[\begin{array}{c}{q}_2^1\\ q_2^2\\ q_2^3\\ q_2^4\end{array}\right]=\left[\begin{array}{c}\mathrm{i}\\ -{\omega }_{20}\\ \frac{\mathrm{i}\left({\omega }_{20}^2-{\omega }^2-{\omega }_1^2\right)}{{\omega }^2}\\ \frac{-{\omega }_{20}\left({\omega }_{20}^2-{\omega }^2-{\omega }_1^2\right)}{{\omega }^2}\end{array}\right](26)$

Now, when ${\alpha }_1={\alpha }_2=0$, Eq. 19 can be written as:

$\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]=\left[\begin{array}{c}-\left(x_2+x_4\right)\\ -{\omega }_{10}{x}_1-{\omega }_{20}{x}_3\\ \frac{{\omega }^2+{\omega }_1^2-{\omega }_{10}^2}{{\omega }^2}{x}_2+\frac{{\omega }^2+{\omega }_1^2-{\omega }_{20}^2}{{\omega }^2}{x}_4\\ \frac{{\omega }^2+{\omega }_1^2-{\omega }_{10}^2}{{\omega }^2}{\omega }_{10}{x}_1+\frac{{\omega }^2+{\omega }_1^2-{\omega }_{20}^2}{{\omega }^2}{\omega }_{20}{x}_3\end{array}\right](27)$

Equation 27 is substituted into Eq. 21 and is transformed into the form of Eq. 24. According to Equations 10–13, one can obtain:

$\begin{array}{l}{p}_{11}=\frac{1}{16m}\frac{F^{‴}\left(-v_c\right)}{\left({\omega }_{10}^2-{\omega }_{20}^2\right){\omega }_{10}}[\left({\omega }^2+{\omega }_1^2-{\omega }_{20}^2\right){\omega }_{10}^3\\ \left.+\frac{{\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right)}^3{\omega }_{10}^3}{{\omega }^4}\right]\end{array}(28)$

$\begin{array}{l}{p}_{12}=\frac{1}{8m}\frac{F^{‴}\left(-v_c\right)}{\left({\omega }_{10}^2-{\omega }_{20}^2\right){\omega }_{10}}[\left({\omega }^2+{\omega }_1^2-{\omega }_{20}^2\right){\omega }_{10}{\omega }_{20}^2\\ \left.-\frac{\left({\omega }^2+{\omega }_1^2-{\omega }_{10}^2\right){\left({\omega }^2+{\omega }_1^2-{\omega }_{20}^2\right)}^2}{{\omega }^4}{\omega }_{10}{\omega }_{20}^2\right]\end{array}(29)$

$\begin{array}{l}{p}_{21}=\frac{1}{8m}\frac{F^{‴}\left(-v_c\right)}{\left({\omega }_{10}^2-{\omega }_{20}^2\right){\omega }_{20}}[\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right){\omega }_{10}^2{\omega }_{20}\\ \left.+\frac{{\left({\omega }^2+{\omega }_1^2-{\omega }_{10}^2\right)}^2\left({\omega }^2+{\omega }_1^2-{\omega }_{20}^2\right)}{{\omega }^4}{\omega }_{10}^2{\omega }_{20}\right]\end{array}(30)$

$\begin{array}{l}{p}_{22}=\frac{1}{16m}\frac{F^{‴}\left(-v_c\right)}{\left({\omega }_{10}^2-{\omega }_{20}^2\right){\omega }_{20}}[\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right){\omega }_{20}^3\\ +\left.\frac{{\left({\omega }^2+{\omega }_1^2-{\omega }_{20}^2\right)}^3{\omega }_{20}^3}{{\omega }^4}\right]\end{array}(31)$

In order to calculate the first-order approximate expression of the real part $\mathbf{\mu }$ of the eigenvalues of matrix $\mathbf{A}\left(\mathbf{\alpha }\right)$ near ${\alpha }_1={\alpha }_2=0$, the adjoint eigenvectors ${\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{p}}_2\left(\mathbf{\alpha }\right)\left(\in C^4\right)$ are introduced, which satisfy

$\left\{\begin{array}{l}{\mathbf{A}}^{\mathrm{T}}\left(\mathbf{\alpha }\right){\mathbf{p}}_1\left(\mathbf{\alpha }\right)={\overline{\lambda }}_1\left(\mathbf{\alpha }\right){\mathbf{p}}_1\left(\mathbf{\alpha }\right)\\ {\mathbf{A}}^{\mathrm{T}}\left(\mathbf{\alpha }\right){\mathbf{p}}_2\left(\mathbf{\alpha }\right)={\overline{\lambda }}_2\left(\mathbf{\alpha }\right){\mathbf{p}}_2\left(\mathbf{\alpha }\right)\end{array}\right.(32)$

The corresponding ${\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{p}}_2\left(\mathbf{\alpha }\right)$ are chosen such that the following relations are satisfied

$<{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_1\left(\mathbf{\alpha }\right)>=<{\mathbf{p}}_2\left(\mathbf{\alpha }\right),{\mathbf{q}}_2\left(\mathbf{\alpha }\right)>=1(33)$

where $<\cdot ,\cdot >$ is the inner product in $C^4$. Accordingly, if $\mathbf{\xi }={\left[{\xi }_1,{\xi }_2,{\xi }_3,{\xi }_4\right]}^{\mathrm{T}}$ and $\mathbf{\eta }={\left[{\eta }_1,{\eta }_2,{\eta }_3,{\eta }_4\right]}^{\mathrm{T}}\in C^4$, then $<\mathbf{\xi },\mathbf{\eta }>={\displaystyle \sum _{i=1}^4}{\overline{\xi }}_i{\eta }_i$. According to the inner product definition, the following relationship is satisfied between the eigenvector and the accompanying eigenvector:

$<{\mathbf{p}}_2\left(\mathbf{\alpha }\right),{\mathbf{q}}_1\left(\mathbf{\alpha }\right)>=<{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_2\left(\mathbf{\alpha }\right)>=0$

By differentiating Eq. 18 with respect to ${\alpha }_j\left(j=\mathrm{1,2}\right)$, we have [26]

$\left\{\begin{array}{l}{\mathbf{A}}_{{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)+\mathbf{A}\left(\mathbf{\alpha }\right){\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)={\lambda }_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)+{\lambda }_1\left(\mathbf{\alpha }\right){\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)\\ {\mathbf{A}}_{{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_2\left(\mathbf{\alpha }\right)+\mathbf{A}\left(\mathbf{\alpha }\right){\mathbf{q}}_{2{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)={\lambda }_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_2\left(\mathbf{\alpha }\right)+{\lambda }_1\left(\mathbf{\alpha }\right){\mathbf{q}}_{2{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)\end{array}\right.(34)$

Taking the inner product of ${\mathbf{p}}_1\left(\mathbf{\alpha }\right)$ with the first equation in Eq. 34, one can obtain

$〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{A}}_{{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉+〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),\mathbf{A}\left(\mathbf{\alpha }\right){\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉=〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\lambda }_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉+〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\lambda }_1\left(\mathbf{\alpha }\right){\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉(35)$

In consideration of the inner product definition, Eq. 35 can be written as

$〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{A}}_{{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉+〈{\mathbf{A}}^{\mathrm{{\rm T}}}\left(\mathbf{\alpha }\right){\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉=〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\lambda }_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉+〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\lambda }_1\left(\mathbf{\alpha }\right){\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉(36)$

It follows from Eq. 32 that

$〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{A}}_{{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉+〈{\overline{\lambda }}_1\left(\mathbf{\alpha }\right){\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉=〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\lambda }_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉+〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\lambda }_1\left(\mathbf{\alpha }\right){\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉(37)$

In view of the inner product definition, Eq. 37 can be written as

$〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{A}}_{{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉+{\lambda }_1\left(\mathbf{\alpha }\right)〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉={\lambda }_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉+{\lambda }_1\left(\mathbf{\alpha }\right)〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉(38)$

From Eq. 33 and eliminating the terms ${\lambda }_1\left(\mathbf{\alpha }\right)〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{q}}_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)〉$, we have:

$〈{\mathbf{p}}_1\left(\mathbf{\alpha }\right),{\mathbf{A}}_{{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_1\left(\mathbf{\alpha }\right)〉={\lambda }_{1{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)(39)$

Taking the inner product of ${\mathbf{p}}_2\left(\mathbf{\alpha }\right)$ with the second equation in Eq. 34 and repeating the processes of Eqs 35–38, then one obtains

$〈{\mathbf{p}}_2\left(\mathbf{\alpha }\right),{\mathbf{A}}_{{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right){\mathbf{q}}_2\left(\mathbf{\alpha }\right)〉={\lambda }_{2{\alpha }_j}^{\prime }\left(\mathbf{\alpha }\right)(40)$

Taken value at $\mathbf{\alpha }=\mathbf{0}$, Eqs 39, 40 can be written into the following unified form:

${\lambda }_{i{\alpha }_j}^{\prime }\left(0\right)=<{\mathbf{p}}_i\left(0\right),{\mathbf{A}}_{{\alpha }_j}^{\prime }\left(0\right){\mathbf{q}}_i\left(0\right)>\left(i=1,2;j=\mathrm{1,2}\right)(41)$

Taking its real part, the change rate of $\mathbf{\mu }$ at ${\alpha }_1={\alpha }_2=0$ can be obtained. Thus, its first-order approximate expression can be easily obtained as:

$\left\{\begin{array}{l}{\mu }_1\approx {\mu }_{1{\alpha }_1}^{\prime }\cdot {\alpha }_1+{\mu }_{1{\alpha }_2}^{\prime }\cdot {\alpha }_2=\mathrm{R}\mathrm{e}\left({\lambda }_{1{\alpha }_1}^{\prime }\left(0\right)\right)\cdot {\alpha }_1+\mathrm{R}\mathrm{e}\left({\lambda }_{1{\alpha }_2}^{\prime }\left(0\right)\right)\cdot {\alpha }_2\\ {\mu }_2\approx {\mu }_{2{\alpha }_1}^{\prime }\cdot {\alpha }_1+{\mu }_{2{\alpha }_2}^{\prime }\cdot {\alpha }_2=\mathrm{R}\mathrm{e}\left({\lambda }_{2{\alpha }_1}^{\prime }\left(0\right)\right)\cdot {\alpha }_1+\mathrm{R}\mathrm{e}\left({\lambda }_{2{\alpha }_2}^{\prime }\left(0\right)\right)\cdot {\alpha }_2\end{array}\right.(42)$

Specifically, let

${\mathbf{p}}_1\left(\mathbf{0}\right)=\left[\begin{array}{c}{p}_1^1\\ p_1^2\\ p_1^3\\ p_1^4\end{array}\right]=\frac{\mathrm{i}{\omega }^4}{2{\omega }_{10}\left[{\omega }^4+{\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right)}^2\right]}\left[\begin{array}{c}{\omega }_{10}\\ \mathrm{i}\\ \frac{{\omega }_{10}\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right)}{{\omega }^2}\\ \frac{\mathrm{i}\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right)}{{\omega }^2}\end{array}\right](43)$

${\mathbf{p}}_2\left(\mathbf{0}\right)=\left[\begin{array}{c}{p}_2^1\\ p_2^2\\ p_2^3\\ p_2^4\end{array}\right]=\frac{\mathrm{i}{\omega }^4}{2{\omega }_{20}\left[{\omega }^4+{\left({\omega }_{20}^2-{\omega }^2-{\omega }_1^2\right)}^2\right]}\left[\begin{array}{c}{\omega }_{20}\\ \mathrm{i}\\ \frac{{\omega }_{20}\left({\omega }_{20}^2-{\omega }^2-{\omega }_1^2\right)}{{\omega }^2}\\ \frac{\mathrm{i}\left({\omega }_{20}^2-{\omega }^2-{\omega }_1^2\right)}{{\omega }^2}\end{array}\right](44)$

From Eq. 16, we have:

${\mathbf{A}}_{{\alpha }_1}^{\prime }=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& {\left.\frac{1}{m}\frac{\partial F^{\prime }\left(-{\alpha }_1-v_c\right)}{\partial {\alpha }_1}\right|}_{{\alpha }_1=0}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& {\left.\frac{1}{m}\frac{\partial F^{\prime }\left(-{\alpha }_1-v_c\right)}{\partial {\alpha }_1}\right|}_{{\alpha }_1=0}\end{array}\right](45)$

${\mathbf{A}}_{{\alpha }_2}^{\prime }=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& -\frac{1}{m}& 0& -\frac{1}{m}\\ 0& 0& 0& 0\\ 0& -\frac{1}{m}& 0& -\frac{1}{m}\end{array}\right](46)$

Substituting equations (25-26) and (43–46) into Eq. 41, one can obtain:

$\begin{array}{l}{\mu }_{1{\alpha }_1}^{\prime }=\mathrm{R}\mathrm{e}{\lambda }_{1{\alpha }_1}^{\prime }=\mathrm{R}\mathrm{e}<{\mathbf{p}}_1\left(\mathbf{0}\right),{\mathbf{A}}_{{\alpha }_1}^{\prime }\left(\mathbf{0}\right){\mathbf{q}}_1\left(\mathbf{0}\right)>\\ =\frac{1}{2m}\frac{\partial F^{\prime }\left(-{\alpha }_1-v_c\right)}{\partial {\alpha }_1}\end{array}(47)$

$\begin{array}{l}{\mu }_{1{\alpha }_2}^{\prime }=\mathrm{R}\mathrm{e}{\lambda }_{1{\alpha }_2}^{\prime }=\mathrm{R}\mathrm{e}<{\mathbf{p}}_1\left(\mathbf{0}\right),{\mathbf{A}}_{{\alpha }_2}^{\prime }\left(\mathbf{0}\right){\mathbf{q}}_1\left(\mathbf{0}\right)>\\ =-\frac{1}{2m}\frac{{\left({\omega }_{10}^2-{\omega }_1^2\right)}^2}{\left[{\omega }^4+{\left({\omega }_{10}^2-{\omega }^2-{\omega }_1^2\right)}^2\right]}\end{array}(48)$

$\begin{array}{l}{\mu }_{2{\alpha }_1}^{\prime }=\mathrm{R}\mathrm{e}{\lambda }_{2{\alpha }_1}^{\prime }=\mathrm{R}\mathrm{e}<{\mathbf{p}}_2\left(\mathbf{0}\right),{\mathbf{A}}_{{\alpha }_1}^{\prime }\left(\mathbf{0}\right){\mathbf{q}}_2\left(\mathbf{0}\right)>\\ =\frac{1}{2m}\frac{\partial F^{\prime }\left(-{\alpha }_1-v_c\right)}{\partial {\alpha }_1}\end{array}(49)$

$\begin{array}{l}{\mu }_{2{\alpha }_2}^{\prime }=\mathrm{R}\mathrm{e}{\lambda }_{2{\alpha }_2}^{\prime }=\mathrm{R}\mathrm{e}<{\mathbf{p}}_2\left(\mathbf{0}\right),{\mathbf{A}}_{{\alpha }_2}^{\prime }\left(\mathbf{0}\right){\mathbf{q}}_2\left(\mathbf{0}\right)>\\ =-\frac{1}{2m}\frac{{\left({\omega }_{20}^2-{\omega }_1^2\right)}^2}{\left[{\omega }^4+{\left({\omega }_{20}^2-{\omega }^2-{\omega }_1^2\right)}^2\right]}\end{array}(50)$

4 Case study and Hopf-Hopf bifurcation

For the system shown in Figure 1, taking $m=1,k=2,k_1=9,k_2=19$, it can be calculated that:

${\omega }_{10}=4.6244,{\omega }_{20}=3.2580$

It can be easily verified that it satisfies the non-resonance condition. $v_c=1.0541$ can be calculated from the dry friction function expression. Substituting the relevant values into relations (28–31) and (47–50) gives:

$\left\{\begin{array}{l}{p}_{11}=-93.9517;p_{12}=-0.2562;\\ p_{21}=-13.9188;p_{22}=-1.7296.\end{array}\right.(51)$

$\left\{\begin{array}{l}{\mu }_{1{\alpha }_1}^{\prime }=-1.4230;{\mu }_{1{\alpha }_2}^{\prime }=-0.6857;\\ {\mu }_{2{\alpha }_1}^{\prime }=-1.4230;{\mu }_{2{\alpha }_2}^{\prime }=-0.3143.\end{array}\right.(52)$

Let

$\left\{\begin{array}{l}\theta =\frac{p_{12}}{p_{22}}=0.1481\\ \delta =\frac{p_{21}}{p_{11}}=0.1481\end{array}\right.(53)$

$\delta \cdot \theta =0.0219(54)$

From Eq. 52, we have

${\left.\det \left(\frac{\partial \mathbf{\mu }}{\partial \mathbf{\alpha }}\right)\right|}_{\mathbf{\mu }=0}\ne 0(55)$

Substituting Eq. 52 into Eq. 42, one can obtain:

$\left\{\begin{array}{l}{\mu }_1\approx {\mu }_{1{\alpha }_1}^{\prime }\cdot {\alpha }_1+{\mu }_{1{\alpha }_2}^{\prime }\cdot {\alpha }_2=-1.4230\cdot {\alpha }_1-0.6857\cdot {\alpha }_2\\ {\mu }_2\approx {\mu }_{2{\alpha }_1}^{\prime }\cdot {\alpha }_1+{\mu }_{2{\alpha }_2}^{\prime }\cdot {\alpha }_2=-1.4230\cdot {\alpha }_1-0.3143\cdot {\alpha }_2\end{array}\right.(56)$

Substituting Equations (51) and (56) into Eq. 15, the specific form of Eq. 15 can be obtained. It is noted from Eqs. (51), (54) and (55) that the system undergoes Hopf-Hopf bifurcation in the “Simple” case [26]. Furthermore, since $s_1,s_2$ in Eq. 15 do not affect its qualitative dynamical behavior, they are not calculated here. Let ${\rho }_k=r_k^2$. Substituting it into Eq. 15, one can obtain:

$\left\{\begin{array}{l}{\dot{\rho }}_1=2{\rho }_1\left({\mu }_1+p_{11}{\rho }_1+p_{12}{\rho }_2+s_1{\rho }_2^2\right)\\ {\dot{\rho }}_2=2{\rho }_2\left({\mu }_2+p_{21}{\rho }_1+p_{22}{\rho }_2+s_2{\rho }_1^2\right)\end{array}\right.(57)$

From Eq. 51, it is easy to see $p_{11}{p}_{22}>0$ , $p_{11}<0$ , and $p_{22}<0$. Let ${\xi }_1=-p_{11}{\rho }_1$, ${\xi }_2=-p_{22}{\rho }_2$ , and $\tau =2t$ , then Eq. 57 can be written as [26]

$\left\{\begin{array}{l}{\dot{\xi }}_1={\xi }_1\left({\mu }_1-{\xi }_1-\theta {\xi }_2+\mathrm{\Theta }{\xi }_2^2\right)\\ {\dot{\xi }}_2={\xi }_2\left({\mu }_2-\delta {\xi }_1-{\xi }_2+\Delta {\xi }_1^2\right)\end{array}\right.(58)$

where $\mathrm{\Theta }=\frac{s_1}{p_{22}^2}$, $\Delta =\frac{s_2}{p_{11}^2}$ (since $\mathrm{\Theta }$ and $\Delta$ do not affect the system’s qualitative dynamic behavior, they should not be calculated). In Eqs (25), (26), (43) and (44), the vectors q₁, q₂ and p₁, p₂ are eigenvectors, which means kq₁, lq₂ and mp₁, np₂ are also eigenvectors. Indeed, the special selection of these eigenvectors has effects on the coefficients p₁₁, p₁₂, p₂₁ and p₂₂. However, it has no effects upon the coefficients θ and δ, then which has no effects upon the final simplified Equation 58 after these transformations ξ₁ = -p₁₁ρ₁, ξ₂ = -p₂₂ρ₂ and τ = 2t. The process of proof can be found in Ref. [26, p. 383], which is too tedious to be suitable presenting it in this paper.

For all values of ${\mu }_1,{\mu }_2$, $E_0=\left(\mathrm{0,0}\right)$ is the equilibrium point of Eq. 58.

If ${\mu }_1,{\mu }_2$ are both negative, then Eq. 58 has only one stable equilibrium point $E_0=\left(\mathrm{0,0}\right)$, as described with region $①$ in Figure 3.

FIGURE 3

FIGURE 3. Parameter region division in Ref. [26] (corresponding to $p_{11}<0$; $p_{22}<0$; $\delta >0$; $\theta >0$; $\delta \cdot \theta <1$).

If at least one of ${\mu }_1,{\mu }_2$ is positive, then from the branch straight line

$H_1=\left\{\left({\mu }_1,{\mu }_2\right):{\mu }_1=0\right\}$

and

$H_2=\left\{\left({\mu }_1,{\mu }_2\right):{\mu }_2=0\right\}$

the trivial equilibrium points $E_1=\left({\mu }_1,0\right)$ (${\mu }_1>0$) and $E_2=\left(0,{\mu }_2\right)$ (${\mu }_2>0$) are branched at the origin, respectively, as described with region $\mathrm{②③⑪⑫}$ in Figure 3. Notably, the equilibrium point at this time is the equilibrium point of the amplitude equation: $E_1=\left({\mu }_1,0\right)$ and $E_2=\left(0,{\mu }_2\right)$ represent the periodic motion of amplitude ${\xi }_1\ne 0,{\xi }_2=0$ and ${\xi }_1=0,{\xi }_2\ne 0$, respectively.

$\theta \delta -1\ne 0$ is known from Eq. 54. Thus, for sufficiently small $‖\mu ‖$, Eq. 58 has a non-trivial equilibrium point in the neighborhood of origin of the phase space:

$E_3=\left(-\frac{{\mu }_1-\theta {\mu }_2}{\theta \delta -1}+O\left({‖\mu ‖}^2\right),\frac{\delta {\mu }_1-{\mu }_2}{\theta \delta -1}+O\left({‖\mu ‖}^2\right)\right)$

Specifically, $\theta \delta -1=-0.9781$, which is negative. Thus, the parameter region where a non-trivial equilibrium point $E_3$ exists is:

$-\left({\mu }_1-\theta {\mu }_2\right)<0\mathrm{a}\mathrm{n}\mathrm{d}\delta {\mu }_1-{\mu }_2<0(59)$

which is the area below $T_1$ and above $T_2$, described as the following (see the region ⑤ shown in Figure 3):

$T_1:{\mu }_1=\theta {\mu }_2,{\mu }_2>0;$

$T_2:{\mu }_2=\delta {\mu }_1,{\mu }_1>0.$

It should be noted that the “equilibrium point” at this time is the equilibrium point of the amplitude equation. Moreover, the non-trivial equilibrium point $E_3$ represents the coupling of two non-zero vibrations ${\xi }_1\ne 0,{\xi }_2\ne 0$, corresponding to the torus motion of the original Eq. 1 (generally quasi-periodic motion).

The phase trajectories corresponding to the parameter region $\mathrm{①②③⑤⑪⑫}$ are shown in Figure 4.

FIGURE 4

FIGURE 4. Phase trajectories for different regions in Figure 3 in Ref. [26].

As shown in Figure 4, there is only a stable equilibrium point in region $\mathrm{①}$. There are stable periodic motions in both regions $\mathrm{②⑪}$ and $\mathrm{③⑫}$; the difference is that $\mathrm{③⑫}$ has one more unstable periodic motion than $\mathrm{②⑪}$. The unstable periodic motion is saddle point type and cannot be observed by numerical simulations. There are stable torus motion and two unstable periodic motions (saddle point type) in region $\mathrm{⑤}$, and numerical simulations can observe stable torus motion.

According to equations (51) and (53), the specific form of Figure 3 can be drawn, as shown in Figure 5. Now, $T_1$ and $T_2$ are described as:

$T_1:{\mu }_1=0.1481{\mu }_2,{\mu }_2>0$

$T_2:{\mu }_2=0.1481{\mu }_1,{\mu }_1>0$

FIGURE 5

FIGURE 5. The specific form of Figure 3, determined by the specific parameters of this paper.

According to Eq. 56, ${\mu }_1,{\mu }_2$ can be determined from ${\alpha }_1$, ${\alpha }_2$. ${\alpha }_1$ is the increment of the wheel edge velocity, and ${\alpha }_2$ is the damping $c$, as shown in Eq. 16. For given values of ${\alpha }_1$, ${\alpha }_2$, ${\mu }_1,{\mu }_2$ can be calculated according to Eq. 56. The equilibrium position is stable if they belong to the region between the two blue lines in Figure 5. If they belong to the area sandwiched by the blue and red lines in Figure 5, the stable periodic motion can be observed through numerical simulations. Suppose they belong to the area sandwiched by the two red lines in Figure 5. In that case, the stable torus motion can be observed through numerical simulations.

Take the first set of parameters: ${\alpha }_1=0.04,{\alpha }_2=0.04$ , ${\mu }_1=-0.0843$ and ${\mu }_2=-0.0695$ can be calculated from Eq. 56, while their positions are denoted by "*". It is easy to see that they belong to the region sandwiched by two blue lines (see Figure 6). Hence, the equilibrium position of system (1) is stable, and its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are shown in Figure 7.

FIGURE 6

FIGURE 6. The specific region to which ${\mu }_1,{\mu }_2$ belong at the time of ${\alpha }_1=0.04,{\alpha }_2=0.04$.

FIGURE 7

FIGURE 7. (A), (B), (C) are the phase trajectories, and (D) is the Poincaré map. The blue lines or points correspond to transient motions, and the red points correspond to steady-state motions (i.e., the final approach to the equilibrium position).

Take the second set of parameters: ${\alpha }_1=-0.03,{\alpha }_2=0.08$, ${\mu }_1=-0.0122$ and ${\mu }_2=0.0175$ can be calculated from Eq. 56, while their positions are denoted by "*". It is easy to see that they belong to the area sandwiched by the blue and red lines (see Figure 8). Therefore, system (1) has a stable periodic motion, and its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are shown in Figure 9.

FIGURE 8

FIGURE 8. The specific region to which ${\mu }_1,{\mu }_2$ belong at the time of ${\alpha }_1=-0.03,{\alpha }_2=0.08$.

FIGURE 9

FIGURE 9. (A), (B), (C) are the phase trajectories, and (D) is the Poincaré map. The blue line or point corresponds to transient motion, and the red line or point corresponds to steady-state motion (i.e., eventually trending towards periodic motion).

Take the third set of parameters: ${\alpha }_1=-0.04,{\alpha }_2=0.04$ , ${\mu }_1=0.0295$ and ${\mu }_2=0.0443$ can be derived from Eq. 56, and their positions are represented by "*". Thus, it can be seen that they belong to the region sandwiched by two red lines (see Figure 10). Hence, system (1) has a stable torus motion, where its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are presented in Figure 11.

FIGURE 10

FIGURE 10. The specific region to which ${\mu }_1,{\mu }_2$ belong at the time of ${\alpha }_1=-0.04,{\alpha }_2=0.04$.

FIGURE 11

FIGURE 11. (A), (B), (C) are the phase trajectories, (D) is the Poincaré map. The blue line or point corresponds to the transient motion, and the red line or point corresponds to the steady-state motion (i.e., the final trend towards torus motion).

Take the fourth set of parameters: ${\alpha }_1=0.03,{\alpha }_2=-0.08$, ${\mu }_1=0.0122$ and ${\mu }_2=-0.0175$ can be derived from Eq. 56, and their positions are represented by "*". Thus, it can be seen that they belong to the region sandwiched by the red and blue lines (see Figure 12). Hence, system (1) has a stable periodic motion, where its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are presented in Figure 13. To be honest, the parameter ${\alpha }_2$ (i.e., damping coefficient c) in this area need to be negative, which is attainable mathematically but is not achievable physically.

FIGURE 12

FIGURE 12. The specific region to which ${\mu }_1,{\mu }_2$ belong at the time of ${\alpha }_1=0.03,{\alpha }_2=-0.08$.

FIGURE 13

FIGURE 13. (A), (B), (C) are the phase trajectories, and (D) is the Poincaré map. The blue line or point corresponds to transient motion, and the red line or point corresponds to steady-state motion (i.e., eventually trending towards periodic motion).

In order to validate further the expressions (28–31, 47–50) and the analysis from Eq. 56 to Figure 5, the relevant numerical calculations for the system shown in Figure 1 are conducted with another set of parameters; $m=1,k=3,k_1=8,k_2=14$. And then one can calculate that; ${\omega }_{10}=4.2711$; ${\omega }_{20}=3.1237$, which satisfies the non-resonance condition; $v_c$ remains $1.0541$. Substituting the relevant values into Eqs. (28-31) and Eqs. (47-50), the corresponding values and relations in the expressions (51–56) can be written as:

$\left\{\begin{array}{l}{p}_{11}=-15.7657;p_{12}=-0.9645;\\ p_{21}=-10.5105;p_{22}=-1.4468.\end{array}\right.(60)$

$\left\{\begin{array}{l}{\mu }_{1{\alpha }_1}^{\prime }=-1.4230;{\mu }_{1{\alpha }_2}^{\prime }=-0.8536;\\ {\mu }_{2{\alpha }_1}^{\prime }=-1.4230;{\mu }_{2{\alpha }_2}^{\prime }=-0.1464.\end{array}\right.(61)$

Let

$\left\{\begin{array}{l}\theta =\frac{p_{12}}{p_{22}}=0.6667\\ \delta =\frac{p_{21}}{p_{11}}=0.6667\end{array}\right.(62)$

$\delta \cdot \theta =0.4444(63)$

From Eq. 61, we have

${\left.\det \left(\frac{\partial \mathbf{\mu }}{\partial \mathbf{\alpha }}\right)\right|}_{\mathbf{\mu }=0}\ne 0(64)$

Substituting Eq. 61 into Eq. 42, one can obtain:

$\left\{\begin{array}{l}{\mu }_1\approx {\mu }_{1{\alpha }_1}^{\prime }\cdot {\alpha }_1+{\mu }_{1{\alpha }_2}^{\prime }\cdot {\alpha }_2=-1.4230\cdot {\alpha }_1-0.8536\cdot {\alpha }_2\\ {\mu }_2\approx {\mu }_{2{\alpha }_1}^{\prime }\cdot {\alpha }_1+{\mu }_{2{\alpha }_2}^{\prime }\cdot {\alpha }_2=-1.4230\cdot {\alpha }_1-0.1464\cdot {\alpha }_2\end{array}\right.(65)$

According to equations (60) and (62), the specific form of Figure 3 can be drawn, as shown in Figure 14. Now, $T_1$ and $T_2$ are described as:

$T_1:{\mu }_1=0.6667{\mu }_2,{\mu }_2>0$

$T_2:{\mu }_2=0.6667{\mu }_1,{\mu }_1>0$

FIGURE 14

FIGURE 14. The specific form of Figure 3, determined by the specific parameters of this paper.

Taking four sets of parameters of $\left({\alpha }_1,{\alpha }_2\right)$, i.e.,

${\alpha }_1=0.03,{\alpha }_2=0.03,{\alpha }_1=-0.04,{\alpha }_2=0.06,$

${\alpha }_1=-0.06,{\alpha }_2=0.03\mathrm{a}\mathrm{n}\mathrm{d}{\alpha }_1=0.02,{\alpha }_2=-0.06,$

$\left({\mu }_1,{\mu }_2\right)$ can be derived from Eq. 65, which are

${\mu }_1=-0.0683,{\mu }_2=-0.0471,{\mu }_1=0.0057,{\mu }_2=0.0481,$

${\mu }_1=0.0598,{\mu }_2=0.0810\mathrm{a}\mathrm{n}\mathrm{d}{\mu }_1=0.0228,{\mu }_2=-0.0197,$

respectively, and their positions are represented by "*" in Figure 15.

FIGURE 15

FIGURE 15. The specific region to which ${\mu }_1,{\mu }_2$ belong at the time of: (A) ${\alpha }_1=0.03,{\alpha }_2=0.03$; (B) ${\alpha }_1=-0.04,{\alpha }_2=0.06$; (C) ${\alpha }_1=-0.06,{\alpha }_2=0.03$; (D) ${\alpha }_1=0.02,{\alpha }_2=-0.06$.

According to the analysis from Eq. 56 to Figure 5, it can be seen from Figure 15A that the equilibrium position of system (1) is stable, and its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are shown in Figure 16.

FIGURE 16

FIGURE 16. (A), (B), (C) are the phase trajectories, and (D) is the Poincaré map. The blue lines or points correspond to transient motions, and the red points correspond to steady-state motions (i.e., the final approach to the equilibrium position).

It is noted from Figure 15B that "*" belong to the area sandwiched by the blue and red lines. Then, system (1) has a stable periodic motion, and its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are shown in Figure 17.

FIGURE 17

FIGURE 17. (A), (B), (C) are the phase trajectories, and (D) is the Poincaré map. The blue line or point corresponds to transient motion, and the red line or point corresponds to steady-state motion (i.e., eventually trending towards periodic motion).

The location of "*" in Figure 15C means that system (1) has a stable torus motion, where its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are presented in Figure 18.

FIGURE 18

FIGURE 18. (A), (B), (C) are the phase trajectories, (D) is the Poincaré map. The blue line or point corresponds to the transient motion, and the red line or point corresponds to the steady-state motion (i.e., the final trend towards torus motion).

As shown in Figure 15D, "*" belong to the region sandwiched by the red and blue lines, which indicates that system (1) has a stable periodic motion, where its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are presented in Figure 19.

FIGURE 19

FIGURE 19. (A), (B), (C) are the phase trajectories, and (D) is the Poincaré map. The blue line or point corresponds to transient motion, and the red line or point corresponds to steady-state motion (i.e., eventually trending towards periodic motion).

5 Conclusions and discussions

5.1 Discussions

Comparing this paper with the previous works, the non-linear resonant terms are given analytically for a system without symmetry. In the previous studies: the non-linear resonant terms were not given and only the eigenvalues of linearization matrix and their derivatives were computed [18]; the non-linear resonant terms were obtained by numerical methods [19]; the non-linear resonant terms were derived analytically for a system with O (2) symmetry [24], where the non-linear resonant terms were simplified tremendously by the O (2) symmetry of the system. In the previous Ref. [18], the Hopf-Hopf bifurcation is “Subcritical” and the periodic and torus motions are unstable, where an inversion of time (i.e., t→-t) was needed for the visualization of the periodic and torus motions. However, the current paper gives an example of dynamical system in which “Supercritical” Hopf-Hopf bifurcation could occur; this complements the previous study.

In fact, the methods used in this paper can also be used to analyze the periodic and torus motions of high-dimensional systems. Nevertheless, before calculating the non-linear resonant terms, a longitudinal simplification is needed to reduce the high-dimensional systems to a set of four-dimensional ordinary differential equations, where the center manifold theory will be applied to obtain the reduced equations and the whole computational process could be tedious and numerical. For infinite dynamical systems, e.g., the analysis of vibration characteristics of structures with geometry non-linear caused by large deformation, the current methods are also applicable after reducing the non-linear partial differential equations to a set of four-dimensional ordinary differential equations. However, this process will be more difficult than that of high-dimensional system, which involved the two-point boundary values problem of linear partial differential equation and the projection of function spaces (see Ref. [24] as a simple example). Incidentally, the methods applied in this paper can accurately predict the qualitative dynamical behavior for systems having a pair of purely imaginary eigenvalues for some critical values of the relevant parameters. Nevertheless, the advantages are restricted to some neighborhood of the critical parameters values. Indeed, for the analysis of global bifurcations of mechanical systems, only using the current methods are inadequate.

In lemma 1, there is a non-resonant condition kω₁₀≠lω₂₀, (k, l > 0, k + l ≤ 5). If the eigenvalues ω₁₀, ω₂₀ fail to meet the non-resonant condition, i.e., kω₁₀ equate or equate approximately to lω₂₀ (k, l > 0, k + l ≤ 5), richer dynamical phenomena will happen. Meanwhile, the normal form (9) will become more difficulty and the differential equations in regard to argument variables ϕ₁, ϕ₂ in Eq. 14 should be considered. For this, the detailed analysis remains to be further studied in the future. Some numerical investigations for two kinds of resonant cases are given in what follows.

Case 1. ω₁₀ = 2ω₂₀. Taking $m=1,k=1,k_1=2.23,k_2=11.38$, it can be calculated that:

${\omega }_{10}=3.5338,{\omega }_{20}=1.7669$

They do not satisfy the non-resonance condition because ω₁₀ = 2ω₂₀ (k, l > 0, k + l = 3 ≤ 5). Take the set of parameters: ${\alpha }_1=-0.08,{\alpha }_2=0.005$. Numerical simulations indicate that system (1) has a stable periodic motion whose frequency is near ω₁₀/2 = ω₂₀, and its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are shown in Figure 20.

FIGURE 20

FIGURE 20. (A), (B), (C) are the phase trajectories, and (D) is the Poincaré map. The blue line or point corresponds to transient motion, and the red line or point corresponds to steady-state motion (i.e., eventually trending towards periodic motion).

Case 2. 2ω₁₀ ≈ 3ω₂₀. Taking $m=1,k=0.24,k_1=4.05,k_2=9.06$, it can be calculated that:

${\omega }_{10}=3.0515,{\omega }_{20}=2.0685$

They do not satisfy the non-resonance condition because 2ω₁₀ ≈ 3ω₂₀ (k, l > 0, k + l = 5 ≤ 5). Take the set of parameters: ${\alpha }_1=-0.05,{\alpha }_2=0.01$. Numerical simulations indicate that system (1) has a stable periodic motion whose frequency is near ω₁₀/3 = ω₂₀/2, and its phase trajectories and Poincaré map (with ${\dot{z}}_2=0$ as the Poincaré section) are shown in Figure 21.

FIGURE 21

FIGURE 21. (A), (B), (C) are the phase trajectories, and (D) is the Poincaré map. The blue line or point corresponds to transient motion, and the red line or point corresponds to steady-state motion (i.e., eventually trending towards periodic motion).

5.2 Conclusions

This paper studied the periodic and torus motions of a two-degree-of-freedom dynamical system with dry friction by theoretical analysis and numerical simulations in the non-resonance case. The Poincaré-Birkhoff normal form was calculated analytically, and the simplified equations were obtained. The analysis of the simplified equations reflected stable periodic and torus motions for suitable parameter regions. Numerical simulations were compatible with the theoretical results. The novelty of this paper is focused on the analytical calculations of non-linear resonant terms for the current asymmetric mechanical system and the discovery of “Supercritical” Hopf-Hopf bifurcation.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation. All authors contributed to the article and approved the submitted version.

Author contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Funding

This work was supported by the National Natural Science Foundation of China [12002096] and the 2022 Doctoral Foundation of Anshun University (asxybsjj202201).

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Leine RI, Nijmeijer H. Dynamics and bifurcations of non-smooth mechanical systems. New York: Springer-Verlag (2004).

Google Scholar

2. Luo ACJ. Discontinuous dynamical systems on time-varying domains. Beijing: Higher Education Press and Springer-Verlag (2009).

Google Scholar

3. Galvanetto U. Dynamics of a three DOF mechanical system with dry friction. Phys Lett A (1998) 248:57–66. doi:10.1016/s0375-9601(98)00644-6

CrossRef Full Text | Google Scholar

4. Stefański A, Wojewoda J, Wiercigroch M, Kapitaniak T. Chaos caused by non-reversible dry friction. Chaos, Solitons and Fractals (2003) 16:661–4. doi:10.1016/s0960-0779(02)00451-4

CrossRef Full Text | Google Scholar

5. Bellido F, Ramírez-Malo JB. Periodic and chaotic dynamics of a sliding driven oscillator with dry friction. Int J Non-Linear Mech (2006) 41:860–71. doi:10.1016/j.ijnonlinmec.2006.05.004

CrossRef Full Text | Google Scholar

6. Guo Y, Xie JH. Neimark-Sacker (N-S) bifurcation of oscillator with dry friction in 1:4 strong resonance. Appl Math Mech (2012) 34:27–36. doi:10.1007/s10483-013-1650-9

CrossRef Full Text | Google Scholar

7. Guo Y. Torus and subharmonic motions of a forced vibration system in 1 : 5 weak resonance. Adv Math Phys (2020) 2020:1–9. doi:10.1155/2020/5017893

CrossRef Full Text | Google Scholar

8. Karami H, Mobayen S, Lashkari M, Bayat F, Chang A. LMI-Observer-Based stabilizer for chaotic systems in the existence of a nonlinear function and perturbation. Mathematics (2021) 9:1128. doi:10.3390/math9101128

CrossRef Full Text | Google Scholar

9. Gritli H, Belghith S. Diversity in the nonlinear dynamic behavior of a one-degree-of-freedom impact mechanical oscillator under OGY-based state-feedback control law: Order, chaos and exhibition of the border-collision bifurcation. Mechanism Machine Theor (2018) 124:1–41. doi:10.1016/j.mechmachtheory.2018.02.001

CrossRef Full Text | Google Scholar

10. Li C, Fan J, Yang Z, Xue S. On discontinuous dynamical behaviors of a 2-DOF impact oscillator with friction and a periodically forced excitation. Mechanism Machine Theor (2019) 135:81–108. doi:10.1016/j.mechmachtheory.2019.01.020

CrossRef Full Text | Google Scholar

11. Stefani G, De Angelis M, Andreaus U. Scenarios in the experimental response of a vibro-impact single-degree-of-freedom system and numerical simulations. Nonlinear Dyn (2020) 103:3465–88. doi:10.1007/s11071-020-05791-4

CrossRef Full Text | Google Scholar

12. Stefani G, De Angelis M, Andreaus U. Influence of the gap size on the response of a single-degree-of-freedom vibro-impact system with two-sided constraints: Experimental tests and numerical modeling. Int J Mech Sci (2021) 206:106617. doi:10.1016/j.ijmecsci.2021.106617

CrossRef Full Text | Google Scholar

13. Stefani G, De Angelis M, Andreaus U. Numerical study on the response scenarios in a vibro-impact single-degree-of-freedom oscillator with two unilateral dissipative and deformable constraints. Commun Nonlinear Sci Numer Simulation (2021) 99:105818. doi:10.1016/j.cnsns.2021.105818

CrossRef Full Text | Google Scholar

14. Peng Y, Fan J. Discontinuous dynamics of a class of 3-DOF friction impact oscillatory systems with rigid frame and moving jaws. Mechanism Machine Theor (2022) 175:104931. doi:10.1016/j.mechmachtheory.2022.104931

CrossRef Full Text | Google Scholar

15. Wang J, Shan Za., Chen S. Bifurcation and chaos analysis of gear system with clearance under different load conditions. Front Phys (2022) 10. doi:10.3389/fphy.2022.838008

CrossRef Full Text | Google Scholar

16. Han X, Bi X, Sun B, Ren L, Xiong L. Dynamical analysis of two-dimensional memristor cosine map. Front Phys (2022) 10. doi:10.3389/fphy.2022.911144

CrossRef Full Text | Google Scholar

17. Akram S, Nawaz A, Yasmin N, Ghaffar A, Baleanu D, Nisar KS. Periodic solutions of some classes of one dimensional non-autonomous equation. Front Phys (2020) 8. doi:10.3389/fphy.2020.00264

CrossRef Full Text | Google Scholar

18. Wen G, Xu H, Xie J. Controlling Hopf–Hopf interaction bifurcations of a two-degree-of-freedom self-excited system with dry friction. Nonlinear Dyn (2010) 64:49–57. doi:10.1007/s11071-010-9844-x

CrossRef Full Text | Google Scholar

19. Guo P, Huang C, Zeng J, Cao H. Hopf–Hopf bifurcation analysis based on resonance and non-resonance in a simplified railway wheelset model. Nonlinear Dyn (2022) 108:1197–215. doi:10.1007/s11071-022-07274-0

CrossRef Full Text | Google Scholar

20. Dhooge A, Govaerts W, Kuznetsov YA. Matcont: A MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans Math Softw (Toms) (2003) 37:141–64. doi:10.1145/779359.779362

CrossRef Full Text | Google Scholar

21. Revel G, Alonso DM, Moiola JL. Numerical semi-global analysis of a 1:2 resonant Hopf–Hopf bifurcation. Physica D: Nonlinear Phenomena (2013) 247:40–53. doi:10.1016/j.physd.2012.12.007

CrossRef Full Text | Google Scholar

22. Xie J, Ding W. Hopf–Hopf bifurcation and invariant torus of a vibro-impact system. Int J Non-Linear Mech (2005) 40:531–43. doi:10.1016/j.ijnonlinmec.2004.07.015

CrossRef Full Text | Google Scholar

23. Wen G-L, Xu D. Feedback control of Hopf–Hopf interaction bifurcation with development of torus solutions in high-dimensional maps. Phys Lett A (2004) 321:24–33. doi:10.1016/j.physleta.2003.12.005

CrossRef Full Text | Google Scholar

24. Guo Y, Xie J, Wang L. Three-dimensional vibration of cantilevered fluid-conveying micropipes—types of periodic motions and small-scale effect. Int J Non-Linear Mech (2018) 102:112–35. doi:10.1016/j.ijnonlinmec.2018.04.001

CrossRef Full Text | Google Scholar

25. Jianhua X, Yuan Y, Denghui L. Nonlinear dynamics. Beijing: Science Press (2018).

Google Scholar

26. Kuznetsov YA. Elements of applied bifurcation theory. Third Edition ed. New York: Springer-Verlag (2004).

Google Scholar