Response Analysis of the Three-degree-of-freedom Vibroimpact System with an Uncertain Parameter

1. Introduction

Collisions and impacts are common physical phenomena in complex engineering structures and systems, which can significantly affect their stability, reliability, and even cause irreversible damage to mechanical equipment. Such systems, referred to as vibroimpact systems, are typically described by differential equations with constraints due to the non-smooth nature of collisions and impacts [1-3]. This inherent non-smoothness introduces strong nonlinearity into the system dynamics, leading to complex behaviors and high sensitivity to parameter changes [4,5]. Moreover, engineering structures are inevitably influenced by random factors, such as air temperature variations, humidity changes, and ground vibrations [6-9]. These factors can introduce randomization characteristics to the system parameters, making the analysis of dynamic properties even more challenging. Thus, studying the stochastic dynamics of vibroimpact systems is not only interesting but also crucial, albeit facing significant challenges.

Due to the non-smoothness and nonlinearity, closed-form solutions for vibroimpact systems are generally unavailable [10]. Numerical simulations play a crucial role in understanding and analyzing the dynamics of these systems [11-13]. By employing a mapping that relates conditions at subsequent impacts，Shaw and Holmes[14-16] discussed the dynamical behaviors of harmonically excited impact systems, including periodic motions, harmonic and subharmonic motions, chaotic motions, and global bifurcations. Mason [17] analyzed codimension-one, -two and -three bifurcations in a periodically-forced impact oscillator with two discontinuity surfaces using discontinuity-geometry methodology. Li [18] investigated the effects of various system parameters on the stochastic P-bifurcation phenomenon in a vibro-impact system using path integration and the generalized cell mapping method. Liu [19] studied the crises in the Duffing vibro-impact oscillator with non-viscously damping by the composite cell coordinate system method. Qian [20] utilized the radial basis function neural networks method to analyze typical randomly excited Vibroimpact systems. Ding [21] established a six-dimensional Poincaré map to explore the double Neimark–Sacker bifurcation, torus T2 of a three-degree-of-freedom vibro-impact system. As research on vibroimpact systems deepens, unique dynamic phenomena have been discovered, including grazing bifurcations [22-24], stick-slip motion [25,26], chattering [27,28], and corner point bifurcation [29,30].

As the study of vibroimpact systems progresses, it becomes evident that these systems exhibit strong parameter sensitivity dependencies [31-34]. This implies that even small variations or adjustments in the system parameters can have a significant impact on the system’s dynamic behavior. Unfortunately, uncertainties and errors are inevitable in various stages such as modeling, measurement, and manufacturing. Even minor deviations can result in fundamental alterations in the dynamic characteristics of the system. Consequently, studying vibro-impact systems with stochastic parameters becomes highly essential. Feng [35] studied the period doubling phenomenon of a single-degree-of-freedom vibroimpact system with an uncertain parameter by means of Chebyshev polynomial approximation method. However, research on vibroimpact systems with multiple degrees of freedom and uncertain parameters is limited. Although extensive studies have been conducted on single-degree-of-freedom vibroimpact systems, extending the analysis to multiple degrees of freedom introduces additional complexities and challenges. The presence of uncertain parameters further exacerbates the complexity of system dynamics.

In this paper, we aim to analyze the dynamic response of a three-degree-of-freedom vibroimpact system with an uncertain parameter using polynomial approximation methods. Section 2 presents the model of a three-degree-of-freedom vibroimpact system with an uncertain parameter. In Section 3, the Chebyshev polynomial approximation is utilized to transform the vibroimpact system with an uncertain parameter into an equivalent deterministic system. Section 4 provides numerical results to demonstrate the effectiveness of our method and analyzes the influence of system parameters on the dynamic behavior of the system. Finally, the paper concludes with key findings.

2. Model of the Three-degree-of-freedom Vibroimpact System with an Uncertain Parameter

Vibroimpact systems often involve multiple impacts with rigid barriers, and the model represented by differential equations with constraints is referred to as the classical impact model. A schematic of the three-degree-of-freedom vibroimpact system with unilateral rigid constrain [36] is shown in Figure 1. The stiffness coefficients of the three springs are $K_1,K_2,K_3$. Three masses $M_1,M_2,M_3$ are attached to the wall by these three springs. Simultaneously, three linear dampers with damping coefficients $C_1,C_2,C_3$ are connected to the masses at the same time. Above every mass, there is a harmonic excitation force$P_i\sin (\omega t+\tau )\{i=1,2,3\}$, respectively. When the first mass $M_1$ moves the distance of $B$, it will impact with the rigid barrier $A$. The energy loss during the impact is determined by the restitution coefficient.

The dimensionless form of the system motion differential equation can be written as

$$M\ddot{X}+2C\dot{X}+KX=P\sin (\omega t+\tau ),(x_1<b)$$

(1)

$${\dot{x}}_{1A+}=-r{\dot{x}}_{1A-},(x_1=b)$$

(2)

where, $\overrightarrow{X}={(x_1,x_2,x_3)}^T$,$\overrightarrow{P}={(p_1,p_2,p_3)}^T$,$\overrightarrow{M}=\left[\begin{array}{ccc}{m}_1& 0& 0\\ 0& m_2& 0\\ 0& 0& m_3\end{array}\right]$,

$$\overrightarrow{K}=\left[\begin{array}{ccc}{k}_1& -k_1& 0\\ -k_1& k_1+k_2& -k_2\\ 0& -k_2& k_2+k_3\end{array}\right],\overrightarrow{C}=\left[\begin{array}{ccc}{\zeta }_1& -{\zeta }_1& 0\\ -{\zeta }_1& {\zeta }_1+{\zeta }_2& -{\zeta }_2\\ 0& -{\zeta }_2& {\zeta }_2+{\zeta }_3\end{array}\right]$$

During the dimensionless form transformation, some assumptions have been done such as $M_1\ne 0$, $K_1\ne 0$ and $P_0=\sqrt{P_1{}^2+P_2{}^2+P_3{}^2}$, and the dimensionless quantities are $m_i=\frac{M_i}{M_1}$,$k_i=\frac{K_i}{K_1}$,$p_i=\frac{P_i}{P_1}$,${\zeta }_i=\frac{C_i}{2\sqrt{K_1{M}_1}}$,$\omega =\mathrm{\Omega }\sqrt{\frac{M_1}{K_1}}$,$t=T\sqrt{\frac{K_1}{M_1}}$,$b=\frac{B{K}_1}{P_0}$,$x_i=\frac{X_i{K}_1}{P_0}$,$i=1,2,3$.

Note that the linear damper $C_i$ is considered to be the random parameter of the system, thus ${\zeta }_i$ can be reduced to the form of ${\zeta }_i={\overline{\zeta }}_i+{\delta }_{i}u$, in which ${\overline{\zeta }}_i$ and $\frac{{\delta }_i}{2}$ are the mean value and variance of $m_i$ respectively, and $u$ is the random variable whose probability density function defined on $[-1,1]$ following the arch distribution. $x_i$ represents the dimensionless displacement of each mass block relative to its initial position. When the mass $M_1$ arrives at the position of the rigid barrier $A$, its velocity will change abruptly, and $r$ is Newton’s restitution coefficient used to characterize the magnitude of energy loss during the impact instant. ${\dot{x}}_{1A+}$ and ${\dot{x}}_{1A-}$ represent the velocity before and after the impact, respectively.

3. The Approximation of the Vibroimpact System with an Uncertain Parameter

3.1. Chebyshev polynomial Approximation

Orthogonal polynomial approximation is an effective method to solve complex nonlinear equations. However, the selection of a polynomial basis must be consistent with the probability density function satisfied by the random parameters in the equation. Though the normal distribution corresponding to Hermite polynomials will provide us useful procedure to deal with the case that system parameter obeys Gaussian distribution, the value range of normal distribution random variable is from negative infinity to positive infinity, which is contrary to the actual situation. Therefore, the random parameter is considered as the random variable subject to arch distribution, which not only meets the implementation situation, but also will not cause the phenomenon of system instability. The probability density function of random variables satisfying the arch distribution is as follow:

$$p(u)=\{\begin{array}{ll}\frac{2}{\pi }\sqrt{1-u^2}& \left|u\right|\le 1,\\ 0& \left|u\right|>1.\end{array}$$

(3)

The curve of $p(u)$ is shown in Figure 2. It is easy to see from the figure that the random parameter is bounded, which is more reasonable in mechanical engineering.

The corresponding polynomial of arch distribution is the Chebyshev polynomial of the second kind, and its general expression is

$$U_n(u)={\displaystyle \sum _{i=0}^{\frac{n}{2}}{(-1)}^i\frac{(n-i)!}{i!(n-2i)!}}{(2u)}^{n-2i}.$$

(4)

Applying Eq. (4), we have

$$\begin{array}{c}{U}_0(u)=1\\ U_1(u)=2u\\ U_2(u)=4u^2-1\\ U_3(u)=8u^3-4u\\ U_4(u)=16u^4-12u^2+1\\ \mathrm{...}\end{array}$$

(5)

Besides, the second kind of Chebyshev polynomials also have a recursive formula

$$u{U}_i(u)=\frac{1}{2}[U_{i-1}(u)+U_{i+1}(u)].$$

(6)

The second type of Chebyshev polynomials has orthogonality, whose weighted orthogonality can be expressed as

$${\displaystyle {\int }_{-1}^1\frac{2}{\pi }\sqrt{1-u^2}{U}_i(u)}{U}_j(u)du=\{\begin{array}{l}1(i=j),\\ 0(i\ne j).\end{array}$$

(7)

Eq. (7) represents the weighted orthogonal relation of polynomials. Although $u$ is a random variable between$-1$ and $+1$, since the weight function is the same as the probability density function of the random variable $u$, the left side of Eq. (7) can be considered as the expectation of the random function $U_i(u)U_j(u)$. According to the nature of function of random variable, the randomness does not affect the orthogonality of polynomial. And because of the orthogonality of the second kind of Chebyshev polynomials, any measurable function $f(u)$ of random variable $u$ can be expanded into the following series:

$$f(u)={\displaystyle \sum _{i=0}^{\infty }{c}_i{U}_i(u)},$$

(8)

in which,

$$c_i=\frac{{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}p(u)f(u)U_i(u)du}}{{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}p(u)U_i(u)U_i(u)du}}={\displaystyle \underset{-1}{\overset{1}{\int }}p(u)f(u)U_i(u)du}.$$

This expansion method based on the generalized Fourier series is called the orthogonal decomposition of random function $f(u)$, and this approximation method is the best mean square approximation.

3.2. Equivalent Deterministic System

In this section, Eq. (1) will be approximated by orthogonal polynomial. Under unconstrained conditions, the response of the system can be expressed as a function of time $t$ and random variable $u$, and can be expanded into the following series form:

$$\overrightarrow{X}(t,u)={\displaystyle \sum _{t=0}^{\infty }{\overrightarrow{X}}_i(t){\overrightarrow{U}}_i(u)}.$$

(9)

When $i$ only takes a finite number $N$, Eq. (9) can be approximately expressed as

$$\overrightarrow{X}(t,u)\approx {\displaystyle \sum _{t=0}^N{\overrightarrow{X}}_i(t){\overrightarrow{U}}_i(u)}.$$

(10)

Where, ${\overrightarrow{U}}_i(u)$ represents the three-dimensional vector composed of the $i-th$ Chebyshev polynomial. We also have$X_i(t)={\displaystyle \underset{-1}{\overset{1}{\int }}p(u)X(t,u)U_i(u)}du$, in which $X_i(t)={(X_{1i},X_{2i},X_{3i})}^T$.

When the impact factor is not taken into account, substituting Eq. (10) is into Eq. (1) yields the following result

$$\overrightarrow{M}\frac{d^2{\displaystyle \sum _{i=0}^N{\overrightarrow{X}}_i(t){\overrightarrow{U}}_i(u})}{d{t}^2}+2\overrightarrow{C}\frac{d{\displaystyle \sum _{i=0}^N{\overrightarrow{X}}_i(t){\overrightarrow{U}}_i(u})}{dt}+\overrightarrow{K}{\displaystyle \sum _{i=0}^N{\overrightarrow{X}}_i(t){\overrightarrow{U}}_i(u})=\overrightarrow{P}\sin (\omega t+\tau ).$$

(11)

By dividing the vector form of Eq. (11) into scalar form, three equations can be obtained

$$\begin{array}{c}(m_1\frac{d^2}{d{t}^2}+k_1){\displaystyle \sum _{i=0}^N{X}_{1i}(t)U_i(u})-k_1{\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})\\ -2\overline{{\zeta }_1}\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{1i}(t)U_i(u})+2{\delta }_{1}u\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{1i}(t)U_i(u})\\ -2\overline{{\zeta }_1}\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})-2{\delta }_{1}u\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})=p_1\sin (\omega t+\tau )\end{array}$$

(12-1)

$$\begin{array}{c}(m_2\frac{d^2}{d{t}^2}+k_1+k_2){\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})-k_1{\displaystyle \sum _{i=0}^N{X}_{1i}(t)U_i(u})-k_2{\displaystyle \sum _{i=0}^N{X}_{3i}(t)U_i(u})\\ +2(\overline{{\zeta }_1}+\overline{{\zeta }_2})\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})-2\overline{{\zeta }_1}\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{1i}(t)U_i(u})\\ -2\overline{{\zeta }_2}\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{3i}(t)U_i(u})+2({\delta }_1+{\delta }_2)u\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})\\ -2{\delta }_{1}u\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{1i}(t)U_i(u})-2{\delta }_{2}u\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{3i}(t)U_i(u})=p_2\sin (\omega t+\tau )\end{array}$$

(12-2)

$$\begin{array}{c}(m_3\frac{d^2}{d{t}^2}+k_2+k_3){\displaystyle \sum _{i=0}^N{X}_{3i}(t)U_i(u})-k_2{\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})-2\overline{{\zeta }_2}\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})\\ +2(\overline{{\zeta }_2}+\overline{{\zeta }_3})\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{3i}(t)U_i(u})+2({\delta }_2+{\delta }_3)u\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{3i}(t)U_i(u})\\ -2{\delta }_{2}u\frac{d}{dt}{\displaystyle \sum _{i=0}^N{X}_{2i}(t)U_i(u})=p_3\sin (\omega t+\tau )\end{array}$$

(12-3)

The recursive relation Eq. (6) of Chebyshev polynomial can be obtained as

$$\begin{array}{c}u{\displaystyle \sum _{i=0}^N{X}_i(t)U_i(u})=\frac{1}{2}{\displaystyle \sum _{i=0}^N{X}_i(t)}[U_{i-1}(u)+U_{i+1}(u)]\\ =\frac{1}{2}{\displaystyle \sum _{i=0}^N{U}_i(u)}[X_{i-1}(t)+X_{i+1}(t)].\end{array}$$

(13)

Due to the approximation of Eq. (10), $X_{-1}$ and$X_{N+1}$ are set to zero. By substituting Eq. (13) into each of the expressions in Eq. (12), and then multiplying both ends of each of the resulting equations by$U_i(u)(i=0,\mathrm{...},N)$, and finally taking the expectation with regard to the random variable $u$, the following result is obtained.

$$\{\begin{array}{c}(m_1\frac{d^2}{d{t}^2}+k_1)X_{10}(t)-k_1{X}_{20}(t)+2\overline{{\zeta }_1}\frac{d}{dt}{X}_{10}(t)+{\delta }_2\frac{d}{dt}{X}_{11}(t)=p_1\sin (\omega t+\tau )\\ (m_2\frac{d^2}{d{t}^2}+k_1+k_2)X_{20}(t)-k_1{X}_{10}(t)-k_3{X}_{30}(t)+2(\overline{{\zeta }_1}+\overline{{\zeta }_2})\frac{d}{dt}{X}_{20}(t)\\ -2\overline{{\zeta }_1}\frac{d}{dt}{X}_{10}(t)-2\overline{{\zeta }_2}\frac{d}{dt}{X}_{30}(t)+({\delta }_1+{\delta }_2)\frac{d}{dt}{X}_{21}(t)\\ -{\delta }_1\frac{d}{dt}{X}_{11}(t)-{\delta }_2\frac{d}{dt}{X}_{31}(t)=p_2\sin (\omega t+\tau )\\ (m_3\frac{d^2}{d{t}^2}+k_1+k_2)X_{30}(t)-k_2{X}_{20}(t)+2(\overline{{\zeta }_3}+\overline{{\zeta }_2})\frac{d}{dt}{X}_{30}(t)-2\overline{{\zeta }_2}\frac{d}{dt}{X}_{20}(t)\\ +({\delta }_3+{\delta }_2)\frac{d}{dt}{X}_{31}(t)-{\delta }_2\frac{d}{dt}{X}_{21}(t)=p_3\sin (\omega t+\tau )\\ (m_1\frac{d^2}{d{t}^2}+k_1)X_{1i}(t)U_i(u)-k_1{X}_{2i}(t)+2\overline{{\zeta }_1}\frac{d}{dt}{X}_{1i}(t)-2\overline{{\zeta }_1}\frac{d}{dt}{X}_{2i}(t)\\ +{\delta }_1\frac{d}{dt}(X_{1i-1}(t)+X_{1i+1}(t))-{\delta }_1\frac{d}{dt}(X_{2i-1}(t)+X_{2i+1}(t))=0\\ (m_2\frac{d^2}{d{t}^2}+k_1+k_2)X_{2i}(t)-k_1{X}_{1i}(t)-k_2{X}_{3i}(t)\\ +2(\overline{{\zeta }_1}+\overline{{\zeta }_2})\frac{d}{dt}{X}_{2i}(t)U_i(u)-2\overline{{\zeta }_1}\frac{d}{dt}{X}_{1i}(t)-2\overline{{\zeta }_2}\frac{d}{dt}{X}_{3i}(t)U_i(u)\\ +({\delta }_1+{\delta }_2)\frac{d}{dt}(X_{2i-1}(t)+X_{2i+1}(t))-{\delta }_1\frac{d}{dt}(X_{1i-1}(t)+X_{1i+1}(t))\\ -{\delta }_2\frac{d}{dt}(X_{3i-1}(t)+X_{3i+1}(t))=0\\ (m_3\frac{d^2}{d{t}^2}+k_2+k_3)X_{3i}(t)U_i(u)-k_2{X}_{2i}(t)-2\overline{{\zeta }_2}\frac{d}{dt}{X}_{2i}(t)\\ +2(\overline{{\zeta }_3}+\overline{{\zeta }_2})\frac{d}{dt}{X}_{3i}(t)+({\delta }_3+{\delta }_2)\frac{d}{dt}(X_{3i-1}(t)+X_{3i+1}(t))U_i(u)\\ -{\delta }_2\frac{d}{dt}(X_{2i-1}(t)+X_{2i+1}(t))=0(i=1,\mathrm{...},N)\end{array}$$

(14)

Therefore, the unconstrained three-degree-of-freedom unilateral rigid constraint vibroimpact system Eq. (1) is reduced to an equivalent deterministic system Eq. (14).

According to the approximation of Eq. (10), the ensemble mean response of the system can be approximated as

$$E[\overrightarrow{X}(t,u)]=E[{\displaystyle \sum _{i=0}^N{\overrightarrow{X}}_i(t){\overrightarrow{U}}_i(u)]={\displaystyle \sum _{i=0}^N{\overrightarrow{X}}_i(t)E[{\overrightarrow{U}}_i(u)]=X_0}},$$

(15)

where,$\overrightarrow{X_0}={(x_{10},x_{20},x_{30})}^T$.

Note that the random variable $u$ takes a random value of$\left\{u_i\right\}$ on the interval$[-1,1]$, there exists a certain time $t$ in the constrained case in which some orbits of the sample system have impacted with the constraint plane $\sum (b)$ (the virtual constraint plane here), while some have not reached the constraint plane $\sum (b)$ due to the influence of random factors. Therefore, it is necessary to average the constraint conditions.

According to Eq. (15), the average response can be obtained. The average constraint surface is:

$$\overline{\sum }=\{(\overrightarrow{X},\overrightarrow{\dot{X}})|\overrightarrow{X_0}=\overrightarrow{b}\}.$$

(16)

The average constraint condition is:

$$E[\overrightarrow{X}(t,u)]=\overrightarrow{X_0}<\overrightarrow{b}.$$

(17)

The average jump equation:

$$\overrightarrow{{\dot{X}}_{0+}}=-r\overrightarrow{{\dot{X}}_{0-}}.$$

(18)

Based on the above Chebyshev polynomial approximation, the introduced average constraint conditions Eq. (17) and the average jump Eq. (18), the stochastic three-degree-of-freedom vibroimpact system is transformed into an equivalent deterministic system. When $N\to \infty$, the solution ${\displaystyle \sum _{i=0}^N{\overrightarrow{X}}_i(t){\overrightarrow{U}}_i(u)}$ approximates $\overrightarrow{X}(t,u)$ in the sense of the best square error. The value of $N$ determines the accuracy of the solution. In this paper, $N=2$ is taken to obtain an equivalent deterministic unilateral constraint system.

4. Reponses of the Three-degree-of-freedom Vibroimpact System

4.1. Period-doubling Bifurcation

Through the analysis of Eq. (1) and Eq. (10), three systems can be obtained: one is the deterministic three-degree-of-freedom vibroimpact system in which the amplitude of the random response$\delta =0$; one is the original stochastic three-degree-of-freedom vibroimpact system; and the last one is the equivalent deterministic system. The corresponding responses of the three systems are deterministic system response (DSR), stochastic system response (SSR) and equivalent system response (ESR), respectively. The response of the original stochastic system is obtained by random simulation using ARM algorithm and Monte Carlo method. While responses of DSR and ESR are obtained by using numerical simulation methods of deterministic systems. The phase orbit diagrams of the responses of the three systems are compared by numerical simulation to study the dynamical behaviors of the stochastic three-degree-of-freedom vibroimpact system. The maximum Lyapunov exponent is a state quantity that can effectively judge the chaotic state of the system. At present, the Lyapunov exponent method for smooth system has been relatively mature. For vibroimpact systems, the Jacobian matrix discontinuity caused by impact factors exists in the model, the calculation method of smooth system cannot be directly applied to the system. In this paper, a Lyapunov exponent calculation method for non-smooth systems based on discontinuous mapping is used to calculate the Lyapunov exponent.

In order to investigate the effect of the frequency of resonant forces, some parameters of the system are taken as follows: the amplitude of random parameter$\delta =0.01$, the restitution coefficient$r=0.7$, and the impact boundary$r=0.7$, resonant forces$P={(1,0,0)}^T$,$m_2=m_3=1.5$,$k_2=k_3=4$. Since $\delta$ is a small quantity, the initial values of both the stochastic system and the deterministic system can be taken as $\overrightarrow{X}={(1,5,0)}^T$and $\overrightarrow{\dot{X}}={(1,0,0)}^T$. For the parameter$\omega =2.25$, the system has a steady-state periodic response, and period is $1T$ for this case, where $T=\frac{2\pi }{\omega }$. This is the so-called $1-1$ periodic motion in which $n-p$ represents periodic motion, $n$ represents the number of cycles and $p$ represents the number of impacts with the constraint surface. The time history diagram and phase orbit diagram are shown in Figure 3(a) and 3(b), respectively. The corresponding results for the parameter $\omega =2.2$ are shown in Figure 3(c) and 3(d), respectively. It is easy to see that $1-1$ motion in Figure 3(b) becomes a stable period $2-2$ motion with the frequency change. That is the period-doubling bifurcation has occurred.

To further observe the dynamical behaviors of the three systems, Figure 4 show the responses when$\omega =1.65$ and $\omega =1.787$. As shown in Figure 4(a) and (b), the system exhibits a stable $1-1$ symmetric motion at $\omega =1.65$, which is close to the motion form of simple harmonic motion. when $\omega =1.787$ , $1-1$ motion has period-doubling bifurcation and becomes a stable period $2-2$ motion, as shown in Figure 4(c) and (d).

In Figure 5, the curve of the maximum Lyapunov exponent is given for different frequencies in Figure 3 and Figure 4. As can be seen from the figure, the maximum Lyapunov exponent at all four positions is stable and less than 0, which means that the system is indeed in a stable state, matching with the phase orbit diagrams. However, with the appearance of period-doubling bifurcation, the maximum Lyapunov exponent also becomes larger in response, which means that the degree of chaos of the system increases.

4.2. From Period-doubling Bifurcation to Chaos

Considering the system parameter$P={(3,0,0)}^T$, $m_2=m_3=2$. The influence on the system responses for different frequency values is investigated. The time history diagram and phase orbit diagram of the system responses are shown in Figure 6. When $\omega =2.45$, the system has periodic response, and the system motion mode is $1-1$, as shown in Figure 6(a) and 6(b). As the angular frequency of the harmonic excitation $P$ decreases to $1-1$, the system motion mode changes to$2-2$ after period-doubling bifurcation, as shown in Figure 6(c) and 6(d). The successive occurrence of such period-doubling bifurcations indicates the arrival of chaos. It can be seen from the maximum Lyapunov exponent shown in Figure 7 that the maximum Lyapunov exponent changes as the parameters change. When the angular frequency of harmonic excitation $P$ decreases to$\omega =2.21$. The maximum Lyapunov exponent becomes a positive number, which means the system appears chaotic, as shown in Figure 8.

As can be seen from Figure 7 and Figure 8, when other parameters have been determined, as the angular frequency keeps decreasing, the system responses produce abundant period-doubling bifurcation phenomenon in the process from$\omega =2.45$ to $\omega =2.21$, and finally lead to chaos. The numerical simulation results show that the phase diagrams of the three systems agree well with each other in the path from period-doubling bifurcation to chaos, which reflects the evolution process of the system.

4.3. Influence of the Restitution Coefficient

In the stochastic three-degree-of-freedom vibroimpact system Eq. (1) and Eq. (2), The restitution coefficient $r$ represents the energy loss in the impact process, which will further affect the periodic phenomenon of the system and lead to bifurcation. The parameters are considered as:$P={(1,0,0)}^T$,$m_2=m_3=1.5$, and $k_2=k_3=4$. Bifurcation diagram with respect to restitution coefficient $r$ for deterministic three-degree-of-freedom vibroimpact system is shown in Figure 9. It is obvious that the system response experienced two period doubling bifurcations, and then leading to chaos.

To clearly show the influence on the system responses for different restitution coefficients, the phase diagrams are presented in Figure 10. It is easy to see that the larger the restitution coefficient $r$ is, the smaller the energy loss during the impact of the system is. That is the greater the energy available when the system is stable. Thus, different stable phase orbital diagrams may be generated. As shown in Figure 10(a), when $r<0.295$, the energy left by the system is no longer enough to reach the impact surface, so a motion state similar to simple harmonic motion will be formed. When $0.295<r<0.613$, the first block has the energy to reach the constraint surface, so the impact phenomenon occurs shown in Figure 10(b). When $0.613<r<0.776$, the system energy further increases and the period-doubling bifurcation occurs, as shown in Figure 10(c). When $r>0.776$, the system enters a chaotic state, as shown in Figure 10(d).

4.4. Influence of the Uncertain Parameter

For the stochastic three-degree-of-freedom vibroimpact system, When the random variable $u$ randomly takes a series of values $\left\{u_i\right\}$ in the interval $[-1,1]$, each value of $u_i$ corresponds to a certain non-smooth sample system ${\displaystyle \prod u_i}$. If there is a period-doubling bifurcation in the system, a deterministic bifurcation must occur at a critical point. However, a series of sample systems of a stochastic system may not present bifurcation at the same critical point. There must be an interval in which the period-doubling bifurcation may or may not have occurred in the samples of the random system. When the bifurcation parameters pass through this transition interval, almost all samples have bifurcations, and the period-doubling bifurcation of the random system is considered to be completed.

Take system parameters as$P={(1,0,0)}^T$,$m_2=m_3=1.5$,$k_2=k_3=4$. When bifurcation parameters change from$\omega =2.14$ to$\omega =2.28$, the deterministic three-degree-of-freedom vibroimpact system presents reverse period-doubling bifurcation, as depicted in Figure 11(a). Meanwhile, when bifurcation parameters change from$\omega =1.6$ to $\omega =1.85$, period-doubling bifurcation to chaos from Figure 11(b) will also occur.

It can be seen from Figure 12 that, when $\omega =2.135$, the phase orbit diagrams of the deterministic system and the random system will be greatly different. The deterministic system has not yet had bifurcation, but the random system has already had bifurcation. Therefore, under the same initial values and parameters, the bifurcation interval of the stochastic system is indeed different from the bifurcation point of the deterministic system. In this bifurcation interval, the phase orbit diagrams of the stochastic vibroimapct system and the deterministic system are quite different. The period doubling bifurcation has occurred in the former, but not yet in the latter.

5. Conclusion

This paper establishes a model of three-degree-of-freedom unilateral rigid restraint vibroimpact system with an uncertain parameter. The damping coefficient is considered as the uncertain parameter and approximated by the arch distribution. The vibroimpact system with the uncertain parameter is firstly transformed to a deterministic system by Chebyshev polynomial approximation method. Then the numerical method of the deterministic system is used to study the bifurcation behavior of the system under asymmetric constraint. The result from Monte Carlo method and the response of the deterministic system are compared, the period-doubling bifurcation behavior of the system is studied, and the bifurcation behavior of the system is analyzed. The results show that Chebyshev polynomial approximation is an effective method for solving multi-degree-of-freedom systems with uncertain parameter. At the same time, the results of numerical simulations show that the asymmetry of the constraint conditions makes the period-doubling bifurcation of the stochastic constraint system more sensitive to the change of the bifurcation coefficient. The general period response will not change much under the influence of weak random factors. However, at the nodes where period doubling bifurcation occurs, even small random vibrations can have a significant impact. Besides, the method may be further generalized to more general vibroimpact systems with uncertain parameters in future research.