Diversity and Regularity of Periodic Impact Motions of a Mechanical Vibration System with Multiple Rigid Stops

The mechanical model of a two-degree-of-freedom vibration system with multiple rigid stops was established, and the eﬀects of the multiple rigid stops to dynamic characteristics of two mass blocks of the system were studied. The judgment conditions and diﬀerential equations of motion of the system masses impacting rigid stops were analyzed. Based on the multiparameter and multiobjective collaborative simulation analysis, the correlation between the dynamic characteristics of the vibration system and the model parameters is studied. The basic periodic and subharmonic impact motions are analyzed with emphasis on the inﬂuences of dynamical parameters on the mode diversity and the distribution characteristics, and the law of emergence and competition of various periodic impact motions on the parametric plane is revealed. The singular points, the hysteresis transition domains, and the accompanying codimension-two bifurcations, caused by the irreversibility of the transition between adjacent basic periodic impact motions in the low-frequency domain, are analyzed. The reasonable parameter matching range, associated with dynamic characteristic optimization of the system, is determined.


Introduction
e typical nonlinear factors of mechanical systems in engineering practice, such as the geometric relationships, stop conditions, topology, and mode coupling, are considered in this work. Examining nonlinear vibrations helps to describe the dynamic characteristics of a mechanical system more accurately. For mechanical systems with multiple nonsmooth mechanical factors, the interactions and influence of nonsmooth switching points lead to complex dynamic behaviours of the system. is increases the difficulty of performing dynamics research, and the description of the motion of the system becomes difficult with the increase in the number of nonsmooth mechanical factors. For the clearance and stop characteristics in the dynamic study of piecewise smooth mechanical systems, more attention has been paid to unilateral or symmetric constraints, but less attention is paid to the effects of multiple rigid stops contained in multiple mass. However, mechanical systems with multiple rigid stops are common in engineering practice, and their dynamic behaviours become extremely complex due to the influence of multiple rigid stops. Shaw [1] laid a foundation based on impact Poincaré mapping for the study of the dynamics of mechanical systems with clearances and constraints through qualitative analysis and numerical simulation. Whiston [2] used singularity to study the fragmentation mechanism of globally stable manifolds caused by impact mapping discontinuities. Several studies [3][4][5][6] analyzed the dynamic stability of vibration systems with clearances and constraints. Li and Wen [7] studied the bifurcation characteristics and control of a mass-spring-belt friction self-excited vibration system. Scholars around the world have conducted considerable research on the codimension-two bifurcations of vibrating systems, including double codimension-2 Hopf bifurcations [8], Hopf-flip bifurcations [9], codimension-two grazing bifurcations [10,11], flip-grazing bifurcations [12], and fold-grazing bifurcations [12], and they have begun to study codimension-three bifurcations with residual dimensions [13][14][15]. Nordmark [16] found that a special bifurcation occurs during a stable periodic vibration with grazing incidence in an impact oscillator under single parameter control. Hu [17] detected numerous grazing phenomena and the accompanying bifurcation characteristics caused by a forced continuous, piecewise-linear oscillator. Chillingworth [18] studied the nondegraded grazing bifurcation corresponding to the nondegraded chatter vibrations and the microdegraded grazing bifurcation corresponding to degraded chatter vibrations for the discontinuity during low-speed rubbing of a single-degree-of-freedom impact oscillator. Humphries and Piiroinen [19] described the relationship between saddle node and grazing bifurcations from a discontinuous geometric perspective. Budd and Dux [20] determined the region where the chattering vibration of a single-degree-of-freedom impact oscillator exists. Toulemonde and Gontier [21] analyzed the sticking motion responses and characteristics of impact oscillators under harmonic excitations based on the predictor-corrector method. Wagg [22] discussed various nonsmooth events, such as sliding bifurcations, which may occur in the lowfrequency domain during the chatter and sticking of twodegree-of-freedom impact oscillators under low forcing frequencies. Pavlovskaia and Wiercigroch [23] conducted experimental and semianalytical studies on the dynamic characteristics of linear oscillators supported by secondary elasticity under scratching conditions. Smooth and nonsmooth bifurcations observed in the experiment were explained by simulations based on mapping solutions of local smooth subspaces. e characteristics and transition process of chattering-impact and complete chattering-impact vibrations were further studied by numerical simulations [24][25][26]. Luo et al. considered the low-frequency vibration characteristics of a periodic forced vibration system with unilateral [27] and symmetric rigid stops [28] and found that the irreversibility of the mutual transfer of adjacent fundamental periodic impact vibrations caused the transition domain node to become a double grazing codimension-two bifurcation. Wang [29] experimented with a vibration system with symmetric elastic constraints and verified the fact that the test results of the two electronic circuits were basically the same and were consistent with the numerical simulation results of the system. Along with the development of engineering applications and the in-depth study of nonsmooth dynamics, impact vibration dynamics have been widely used in mechanical systems. Relevant research results have revealed the basic dynamic characteristics of various mechanical vibration systems and broadened the knowledge of the nonsmooth characteristics of such systems, such as the three-axle locomotive bogie system [30,31], ground moling [32], the fly-wheel model of a bouncing ball [33], ultrasonic percussive drilling [34], vibroimpact capsule systems [35], pressure relief valves [36], Jeffcott rotors with bearing clearance [37], impact dampers [38], and gear drives [39].
A two-degree-of-freedom vibratory system with multiple rigid stops is considered in this paper, which is similar to the mechanical model with mutual impact stops between the bolster side frame system and the vehicle body, as well as the axle of the heavy-duty truck. e characteristics and occurrence mechanism of the basic periodic motion groups and chattering-impact vibration in the low-frequency domain are analyzed. Based on the multiparameter and multiobjective cosimulation analysis, the incidence relation between dynamic characteristics of the vibration system and the model parameters is studied. e hysteresis transition domains, singular points, and codimension-two bifurcations induced by the irreversibility of the transition between adjacent basic periodic motions in each group are discussed. e reasonable parameter matching range for dynamic characteristic optimization of the system is determined.

The Mechanical Model
e mechanical model of a mechanical vibration system with multiple rigid stops is presented in Figure 1. e two vibrating masses in the system are denoted by M i (i � 1,2), and the corresponding displacement is recorded as X i . e two masses (M i ) are connected to the stops by a linear spring with a stiffness of K i and a linear damper with a damping coefficient of C i . e harmonic exciting forces (P i sin(ΩT + τ)) of amplitude P i , frequency Ω, and phase angle τ act on the corresponding mass (M i ). When the harmonic excitation force is small, the system shows a noimpact forced vibration. e mass will impact on the corresponding rigid stops when the displacement of the mass moving towards the stops increases to the clearance threshold with the increase of the harmonic excitation force. e impact recovery coefficient of the system is R. e dimensionless parameters, variables, and time are defined as follows: 2 Shock and Vibration e dimensionless differential equation of motion is as follows: According to the law of conservation of momentum and the coefficient of restitution of impact, the following can be deduced: (1) When the displacement of the mass M 1 is equal to the clearance threshold − δ 1 (i.e., x 1 � − δ 1 ), the mass M 1 impacts the bottom stop A 1 , and the velocities of the mass M 1 immediately before and after the impact are written by the impact law where the subscript signs "-" and "+" denote the states just before and after the impact, respectively. If impact velocity _ x 1 � 0 and the external force acting on mass M 1 is (F 1 < 0), then us, mass M 1 completes the chattering-impact vibration at stop A 1 . When external force F 1 increases to zero, then mass M 1 is separated from bottom stop A 1 .
(2) When the displacement of mass M 2 is equal to the clearance threshold − δ 2 (i.e.x 2 � − δ 2 ), then the mass M 2 impacts the bottom stop A 2 , and the impact equation of the mass M 2 is given by If impact velocity _ x 2 � 0, and the external force acting on mass M 2 is F 2 < 0, then us, mass M 2 completes the chattering-impact vibration at stop A 2 . When external force F 2 increases to zero, mass M 2 is separated from stop A 2 .
(3) When the relative displacement of the masses M 1 and M 2 is equal to the clearance threshold δ 12 (i.e.x 1 − x 2 � δ 12 ), then the velocities of these two masses are changed according to the conservation law of momentum, and the impact equation is written as If relative impact velocity _ x 1 − _ x 2 � 0 and the combined external force acting on the two masses is F 1 − F 2 > 0, then e two masses exhibit a complete chattering-impact vibration at rigid stopA 12 , and when the combined external force (F 1 − F 2 ) decreases to zero, the two masses separate.
In formulas (3), (5), and (7), _ x i− and _ x i+ are the respective instantaneous velocities before and after the impact of the mass M i . e value ranges of some parameters in dimensionless parameter (1) are μ m ∈ (0, 1), μ k ∈ (0, 1), μ c ∈ (0, 1), and f 20 ∈ [0, 1]. In this way, the correlation and matching law between the dynamic characteristics and model parameters of the system are discussed globally. To systematically analyse the mode types of periodic impact vibration and the bifurcation characteristics, the periodic impact motions of mass M i are represented as n-p-q, where n is the ratio of the vibration period (T n � 2nπ/ω) to the period of the excitation force (T 0 � (2π/ω)) (n � 1, 2, 3, . . .), and p and q represent the number of impacts of mass M 1 or M 2 and the rigid stops, respectively (p, q � 0, 1, 2, 3, . . .). To calculate the values of n, p, and q, Poincaré section σ k I (k � 1, 12, 2) is selected as the system state variable at the moment of impact at rigid stopsA 1 , A 12 , and A 2 in the system, and Poincaré section σ II is selected as the system state variable at the minimum displacement moment of mass M 1 in the period of excitation force. Four kinds of Poincaré section are defined as follows: For periodic impact vibration, the number of fixed points on Poincaré section σ II is the number of exciting force periods n, the numbers of fixed points on Poincaré sections σ 12 I and σ 1 I , respectively, represent the impact times (p and q) of mass M 1 at rigid stops A 12 and A 1 , and the fixed points on Poincaré sections σ 2 I and σ 12 I represent the impact times (p and q) of mass M 2 at rigid stops A 2 and A 12 . erefore, according to the number of branches in the bifurcation map corresponding to Poincaré sections σ II and σ k I , the mode type (n-p-q) of the periodic impact vibration of mass M 1 or M 2 can be determined. e impact Poincaré mapping of the system is expressed as 4 , and υ are the parameters; i.e., υ ∈ R 8 .

Periodic Motions and Bifurcations of the Vibration System with Multiple Clearances
A two-degree-of-freedom forced vibration system with multiple rigid stops was considered, and the correlations between the dynamic characteristics and model parameters of the system were revealed. e effects of key parameters, such as the clearance threshold (δ) and excitation frequency (ω), on the dynamic characteristics of the system are discussed in this section. Dimensionless parameters μ m � 0.5, 8, and f 20 � 0 are selected as the reference parameters. rough the multiparameter and multiobjective cosimulation analysis, the mode types, parameter domains, and bifurcation characteristics of the periodic impact motions of masses M 1 and M 2 in the (ω, δ)-parameter plane are calculated, as shown in Figure 2. In the (ω, δ)-parameter plane obtained in this study, different types of periodic impact motions are distinguished by different colours and the symbols (n-p-q). Furthermore, 70% of the black domain marked with "C" in the two-parameter figure indicates chaos or long-period multi-impact motions. is finding can be attributed to the small domains of these kinds of long-period multi-impact motions, which are consistent with the dynamic characteristics of the system and therefore are not accurately identified further. e 1 − p − q(p ≥ 6 or q ≥ 6) symbol indicates chattering-impact vibration, and the overbars above the symbols 1 − p − q, 1 − p − q, and 1 − p − q indicate that the mass exhibits complete chattering-impact vibration at this rigid stop. Figure 2 shows that, in the large clearance threshold (δ) domain, mass M 1 mainly undergoes 1-0-0 and 1-1-1 periodic motions, while mass M 2 is mainly characterised by 1-0-0 and 1-0-1 periodic motions and tongue-shaped regions on the boundary of their adjacent periodic motion transitions. In the small clearance (δ) domain, the system exhibits abundant basic periodic impact motions (such as 1- chattering-impact vibration, complete chattering-impact vibration, and subharmonic motions in the nonhysteresis transition domains (tongue-shaped regions). In the lowfrequency (ω) domain with a small clearance (δ), a series of grazing-contacts occur between each mass of the system and its corresponding rigid stops with the decrease in ω or δ. e 1 − p − q basic periodic impact motion undergoes a grazing bifurcation transition to 1 − (p + 1) − q (or 1 − p − (q + 1)) periodic impact motion, resulting in an increase in the number of impacts (p (or q)) of each mass during this period. Furthermore, the system exhibits chattering-impact characteristics when the number of impacts p (or q) is large enough. As this dynamic characteristic persists, if (4), (6), or (8) is satisfied at this time, then the two masses of the system exhibit corresponding complete chattering-impact vibration. Under the reference parameters, when the clearance threshold (δ) is large, periodic impact motions such as 1-1- , and 4-3-3 occur in the mass M 1 , while mass M 2 is accompanied by periodic impact motions such as 2-0-0, 2-0-1, 2-0-2, 3-0-1, 4-0-2, and 4-0-3. In the parameter domain of low frequency (ω) and small clearance (δ), mass M 1 not only exhibits basic periodic motions such as 1-1-1, 1-2-2, 1-3-3, 1-4-4, and 1-5-5 of symmetrical impact number, but also results in basic periodic motions such as 1-2-1, 1-3-2, 1-4-3, 1-5-3, and 1-5-4 of asymmetrical impact number. At the same time, the mass M 2 exhibits basic periodic motions such as 1-0-1, 1-0-2, 1-1-2, 1-0-3, 1-1-3, 1-0-4, 1-1-4, 1-0-5, and 1-1-5. According to the calculation results, many periodic impact motions (n-p-q) occur with the asymmetric impacts in the periodic impact vibration modes of the system due to the influence of multiple rigid stops. In the low-frequency (ω) domain with a small clearance (δ), the system exhibits not only the basic periodic motions of symmetric impact times (1-p-p) but also a large number of competitive basic periodic motion groups with asymmetric impact times (1-pq) (p ≠ q). Figure 3 illustrates the local details of Figure 2(a). In the small clearance (δ) parameter domain, the 1 − p − p(p > 0) periodic motion of mass M 1 transitions to 1 − (p + 1) − (p + 1) periodic motion via a real-grazing bifurcation or to subharmonic motion or chaos via a bare-grazing bifurcation with decreasing frequency ω. As frequency ω or clearance δ increases, the 1 − (p + 1) − (p + 1) motion changes from a saddle-node bifurcation to 1 − p − p motion, or the perioddoubling bifurcation causes 2 − 2(p + 1) − 2(p + 1) motion. On the grazing bifurcation boundary G 1− p− p , mass M 1 undergoes grazing contact with rigid stops at both ends (A 12 and A 1 ), and the (p + 1) impact velocity is zero ( _ x 1 � 0). At this time, the number of impacts of mass M 1 increases by one, and the vibration period is unchanged, nearly moving into stable 1 − (p + 1) − (p + 1) motion. e grazing bifurcation of this feature is defined as a real-grazing bifurcation. Naturally, another kind of bifurcation exists, i.e., the bare-grazing bifurcation [40], which is characterised by the unstable 1 − (p + 1) − (p + 1)motion induced by the bifurcation, and the system then moves into subharmonic motion or chaos. As the frequency ω decreases continuously, the system experiences a series of real-grazing bifurcations, increasing the number of impacts p and decreasing the impact velocity _ x i . When the number of impacts (p) is large enough, the system exhibits chattering-impact vibration and undergoes a further transition from a sliding bifurcation to a complete chattering-impact vibration [41,42]. eoretically, the complete chattering-impact vibration is characterised by infinite impact sequences and successive impact velocity attenuation in a one period vibration, resulting in a sticking phenomenon of mass M 1 or M 2 at the corresponding rigid stops. It can be seen from Figure 3 that when the clearance threshold is small and excitation frequency is extremely small, the system presents basic periodic motion groups, chattering-impact vibration, and complete chattering-impact vibration. In Figure 3, the clearance threshold is set at δ � 0.6, and, as the excitation frequency decreases, the time history diagram of the evolution of mass M 1 from basic periodic motions (1 − p − q) to complete chattering-impact vibration (1 − p − q) is shown in Figure 4. e clearance threshold and excitation frequency of the system are δ � 0.6 and ω � 0.33, respectively, and mass M 1 exhibits a 1-4-4 periodic impact motion. Mass M 1 undergoes an increase in the number of impacts caused by the real-grazing bifurcation and exhibits a 1-5-6 periodic impact motion with the decrease in the excitation frequency to ω � 0.295. When the excitation frequency is further reduced to ω � 0.275, mass     bifurcation boundary (SN 1− (p+1)− (p+1) ) and the vibrating real-grazing bifurcation boundary (G 1− p− p ); see the 40% of black domains marked as HR 1 , HR 2 , and HR 3 in Figure 3(c). e parameter domain marked as HR 3 represents a hysteresis transition domain formed by the combination of saddle bifurcation boundary SN 1-3-4 for 1-3-4 period motion and real-grazing bifurcation boundary G 1-3-3 for 1-3-3 period motion. Due to the different initial conditions of the system, two adjacent basic periodic motions (1-3-3 and 1-3-4) coexist in this parameter domain. e phase portrait of the mass M 1 shows that the 1-3-3 periodic motion in the hysteresis transition domain HR 3 is transferred to 1-3-4 periodic motion by real-grazing bifurcation, as shown in Figure 5. In order to clearly observe the evolution of realgrazing bifurcation, local details of mass M 1 at the A 1 rigid stop are superimposed in Figure 5. When the clearance threshold is δ = 0.7 and the excitation frequency is ω = 0.3149, the mass M 1 exhibits a 1-3-3 periodic motion. When the frequency is ω = 0.3147, the 1-3-3 periodic motion of mass M 1 makes grazing contact with the A 1 rigid stop. As the frequency ω continues to decrease, it evolves into a stable 1-3-4 period motion. Similar phenomena occur in adjacent periodic motions within various basic periodic impact groups. e nonhysteresis transition domain is a tongue-shaped region enclosed by the baregrazing bifurcation boundary of G b 1− p− p and the perioddoubling bifurcation boundary (PD 1− (p+1)− (p+1) ) in the (ω, δ)-parameter plane, and the subharmonic motions in the nonhysteresis transition domain change with the excitation frequency ω and clearance threshold δ, exhibiting either a period-doubling or grazing bifurcation. e results in Figure 3 show the mode types of the subharmonic motions within the nonhysteresis transition domain on the boundary of the 1-4-4 and 1-5-5 periodic motions. Based on the calculation results in Figure 3, the mode types and bifurcation characteristics of subharmonic motions in the nonhysteresis transition domain surrounded by periodic doubling bifurcation boundary PD 1-5-5 of 1-5-5 motion and bare-grazing bifurcation boundary G b 1− 4− 4 of 1-4-4 motion are further observed. erefore, Figure 6 shows the local description bifurcation diagram of mass M 1 with excitation frequency ω � 0.272 and clearance threshold δ ∈ [0.77, 0.82]. Figure 7 shows the local description bifurcation diagram of mass M 1 with the clearance threshold δ � 0.61 and excitation frequency ω ∈ [0.296, 0.326]. Figure 6(a) corresponds to Poincaré section σ II , and Yaxis(x 1mp ) represents the minimum displacement of mass M 1 during the period of the excitation force T � (2π/ω). Figure 6(b) corresponds to Poincaré section σ 12 I , and Yaxis( _ x 1 − _ x 2 ) represents the relative impact velocities of the two masses. Figure 6(c) corresponds to Poincaré section σ 1 I , Y-axis( _ x 1 ) represents the impact velocities of mass M 1 and its bottom stop A 1 , and X-axis represents the clearance threshold (δ ∈ [0.77, 0.82]). e number of branches in Figures 6(a), 6(b), and 6(c) can determine the values of n, p, and q in the periodic or subharmonic vibration type (n-p-q) of the system at the corresponding clearance threshold (δ). e coordinate meaning of the bifurcation diagram in Figure 7 is the same as that in Figure 6. Figures 6(d) and 6(g) are local descriptive periodic bifurcation diagrams of Figure 6(a); Figures 6(e), 6(h), 6(f ), and 6(i) are local descriptive impact number bifurcation diagrams of Figures 6(b) and 6(c), respectively. From Figure 6, it can be seen that when the clearance threshold δ decreases, the 1-4-4 periodic motion migrates into 3-13-12, 2-9-8, and 4-19-16 subharmonic motions and chaos via periodic doubling-grazing bifurcation and then evolves into 2-10-8 and 1-5-4 motions via inverse periodic doubling bifurcation. When the clearance threshold δ decreases continuously, the 1-5-4 period motion evolves into the paroxysmal chaos, which induces subharmonic periodic motions such as 6-31-24⟶C⟶5-26-20⟶C⟶4-21-16⟶C⟶3-16-12⟶C⟶2-11-8⟶C. Alternatively, as the clearance threshold δ increases, the 1-5-5 periodic motion undergoes periodic doubling or saddle bifurcation and then evolves into 2-10-10, 2-10-9, and 4-20-18 subharmonic periodic motions and chaos. erefore, there are periodic impact motions, and subharmonic motions n − (np + 1) − np and 2n − 2(np + 1) − 2np exist in the nonhysteresis transition domain with the change in the excitation frequency ω or clearance threshold δ; i.e., the period-doubling bifurcation sequence evolves into chaos or the bare-grazing bifurcation transitions to chaos. When the clearance threshold value is δ � 0.61 and with the increase in the excitation frequency ω (Figure 7), the mass M 1 evolves from stable 1-5-5 periodic motion to 2-9-10 periodic motion through doubling-saddle codimensiontwo bifurcations. With the increase in the excitation frequency ω, the mass M 1 transforms into 1-4-5 periodic motion via a saddle-inverse doubling bifurcation. When the excitation frequency increases to ω � 0.3166, the mass M 1 again evolves into 2-8-9 periodic motion through a doubling-saddle bifurcation. e excitation frequency ω increases continuously, and the mass M 1 moves into stable 1-4-4 periodic motion through saddle-inverse doubling codimension-two bifurcations. e simulation results show that the irreversibility of the transition between 1 − p − p and 1 − (p + 1) − (p + 1) modes of two adjacent basic periodic motions results in an alternating distribution of the hysteresis transition domains and nonhysteresis transition domains along their boundary lines and the singularity is the node of the two types of transition domains. e singularity is the twofold-grazing bifurcation point of the 1 − p − p periodic vibration, i.e., the codimension-two bifurcation point of the saddle-node and the period-doubling bifurcation of the 1 − (p + 1) − (p + 1) periodic motion (Figures 3, 6, and 7).
Impact vibration experiments of flexible beams were designed previously [43,44], and the vibration system of a mass's relative collision and unilateral constraint was studied.
e experimental results verified the types and existence of basic periodic impact vibration groups in the low-frequency domain, and the typical nonsmooth properties of the experimental system, such as grazing bifurcations, mutation, and hysteresis, were tested. However, for vibration systems with multiple rigid stops, basic periodic impact vibration groups with various mode types exist in the low-frequency domain.

e Influence of the Mass Ratio on Dynamic Characteristics of the Vibration
System. e mass distribution ratio (μ m ) of the two masses is discretised at equal intervals within the range of its values, and the correlation between the dynamic characteristics and mass ratio μ m of the vibration system in the (ω, δ)-parameter plane is obtained by numerical calculation. Figure 8 shows four representative sets of mass ratios (μ m � 0.05, μ m � 0.25, μ m � 0.75, and μ m � 0.95), and other dimensionless parameters are selected as references: μ c � 0.5, μ k � 0.5, ζ � 0.1, R � 0.8, and f 20 � 0. e damping distribution ratio (μ c ), stiffness distribution ratio (μ k ), and damping coefficient (ζ) in the dimensionless parameters of the system are simulated with the same method as mass ratio μ m to reveal the relationship between the system dynamic characteristics and the model parameters systematically. Figures 8(a) and 8(e) show that when mass ratio (μ m �0.05) is small, the system shows chatteringimpact vibration and complete chattering-impact vibration (1 − p − q) in the small clearance (δ) domain. In the large clearance (δ) domain, masses M 1 and M 2 of the system exhibit a 1-0-0 no-impact free vibration. e system shows abundant mode types of periodic impact vibration in the domain of ω ∈ [0.5, 1.2]. When the clearance threshold (δ) decreases, the periodic motions of mass M 1 exhibit 3-1-1, 2-0-1 motions which move into 2-1-1, 2-2-1 motions via grazing bifurcation sequence and chaos. Subsequently, the mass M 1 exits chaos and evolves into 2 periodic motions via saddle bifurcation sequence and then undergoes inverse doubling bifurcation to 1-0-1. Subsequently, it experiences realgrazing bifurcation sequence and then shows basic periodic impact motion groups such as    low-frequency and small clearance domain. With the further increase of mass ratio μ m , the basic periodic impact motion groups of the system show diversification and competition in the low-frequency and small clearance domain (Figures 2 and  7). During the increase of mass ratio μ m , the occurrence domains of the two types of periodic impact motions (1-1-1 and 1-2-1) of mass M 1 and the two types of periodic impact motions (1-0-1 and 1-0-2) of mass M 2 increase first and then decrease, and the parameter domain of chaos increases significantly in the (ω, δ)-parameter plane.

e Correlation between Dynamic Characteristics and
Damping Ratio. Figures 2 and 9 show the mode types, existence domains, and bifurcation characteristics of the impact vibration of the system corresponding to the typical damping ratio μ c values in the (ω, δ)-parameter plane. e results show that, with the increase of damping ratio μ c , the occurrence domain of the basic periodic impact motion groups (1-p-q) of masses M 1 and M 2 in the (ω, δ)-parameter plane increases obviously and the existence domain of the chatter-impact vibration and complete chattering-impact vibration of the system expands to the direction of increasing excitation frequency ω. However, with the decrease of damping ratio μ c , the system shows that the parameter domains of chaos or periodic motions, such as 2-1-1, 2-1-2, and 2-2-2, gradually increase. In the low-frequency ω and small clearance δ domain, the existence domains of the basic  periodic impact motion groups (1-p-q) exhibit a laminar flow distribution. With the increase of damping ratio μ c , the boundary of the existence domains of the basic periodic impact motion groups (1-p-q) becomes smooth gradually and the chaos decreases. When the damping ratio is μ c � 0.95 and the excitation frequency is ω ∈ [1.05, 2.70], the basic periodic impact motion groups (1-p-q) of the system present a series of stromatolite distributions existing. And the boundary of the adjacent basic periodic impact motion in the parameter domain is also composed of a series of singularities and alternately distributed hysteresis transition and nonhysteresis transition domains, as shown in Figures 9(d) and 9(h). e calculation results show that small damping ratio μ c induces various types of subharmonic motions and chaos, which results in more complex dynamic system characteristics. However, a large damping ratio μ c shrinks the subharmonic motions and chaos domains, and the existence domain of various basic periodic impact motion groups of the system expands significantly. In addition, a complex competition law is observed in the small clearance δ domain, and the mode types of the basic periodic impact motion groups of mass M 1 have shown diverse characteristics.

e Influence of the Stiffness Ratio on Dynamic Characteristics of the Vibration System.
rough the multiobjective coordination, multiparameter, and numerical simulation of the discretisation value, four groups of representative stiffness distribution ratio μ k corresponding to the system impact vibration mode types and distribution domains in the (ω, δ)-parameter plane are selected, as shown in Figure 10. A smaller stiffness ratio μ k indicates that the spring stiffness (K 2 ) of mass M 2 is smaller than that (K 1 ) of mass M 1 to the stop connection. When the stiffness ratio is small enough (μ k �0.05), the parameter domains of the chatteringimpact vibration and complete chattering-impact vibration of mass M 1 are much larger than those of mass M 2 . e system shows many subharmonic motions in the (ω, δ)-parameter plane. Subharmonic motions, such as 2-1-2, 3-1-1, 3-1-2, 3-1-3, 4-1-4, 5-1-4, and 7-1-4, and small domains 2-0-1, 2-0-2, 2-2-2, 2-1-3, 2-1-4, and 4-2-4 are induced by mass M 1 , and mass M 2 exhibits subharmonic motions of 2-0-1, 3-0-1, 4-0-1, 5-0-1, and 7-0-1 and small domains 2-0-0, 2-0-2, and 4-0-2. As stiffness ratio (μ k >0.25) increases gradually, masses M 1 and M 2 show the distribution domains of chattering-impact and complete chattering-impact vibration which expand to the high-frequency ω direction in the (ω, δ)-parameter plane. Furthermore, the existence domain of the 1-1-q (q ≥ 2) basic periodic impact group of masses M 1 and M 2 decreases gradually until they disappear. Meanwhile, the basic periodic impact motion groups 1-p-q (p ≥ 0q ≥ 1) of mass M 1 exhibit diversification and the competition law, and the domains of the 1-0-q (q ≥ 1) basic periodic motion group of mass M 2 expand and then shrink. When the stiffness ratio is μ k � 0.95 and, in the small clearance δ domain, as shown in Figures 10(d) and 10(h), the boundary of the adjacent basic periodic impact motion distribution domain of the system is sharpened and the chaos is obviously increased.
e system shows complex bifurcation characteristics. e system mainly exhibits a noimpact free vibration of 1-0-0 in the high-frequency ω and large clearance δ domain. e calculation results show that the number of impacts at mass M 1 and bottom stop A 1 increases with stiffness ratio μ k and the number of impacts at mass M 2 and bottom stop A 2 decreases.

e Incidence Relation between Dynamic Characteristics and the Damping Coefficient.
Considering the dimensionless damping coefficient in the range of ζ ∈ (0, 0.3], the numerical calculation is carried out by using the value of equal separation dispersion. Figure 11 shows the mode types and occurrence domain of the periodic impact motions in the (ω, δ)-parameter plane of the mechanical system obtained with four different damping coefficients ζ. When the damping coefficient is ζ � 0.25, in the large clearance δ domain, the system mainly exhibits 1-0-0 no-impact periodic vibration, and, in the small clearance δ domain, the mass M 1 shows various basic periodic impact motion groups (1 − p − q) and the mass M 2 shows basic periodic impact motion group (1 − 0 − q). e nonhysteresis transition domain appears in the local parameter domain during the transition between the 1-0-0 and 1-0-1 motions (transitions between the 1-1-1 and 1-2-1 motions or 1-2-1 and 1-3-1 motions are also possible) of mass M 1 . Owing to the linear relationship between the geometry of the mass M 2 and the mass M 1 (Figure 1), a nonhysteresis transition domain also appears in the corresponding parameter domain. As the damping coefficient decreases, the peak values of the 1-1-1 periodic motion of system mass M 1 and the 1-0-1 periodic motion of the mass M 2 gradually increase, along with the enlargement of the hysteresis transition and nonhysteresis transition domains of the system. As the damping coefficient ζ decreases continuously, in the small clearance (δ) threshold and low excitation frequency (ω) region, the domains in which chattering-impact vibration and complete chatteringimpact vibration of the mass occur shift in the low-frequency direction, and the mode types, existing domains, and bifurcation characteristics of various periodic impact motions of the system are also extremely complex. In particular, when the damping coefficient is small (ζ� 0.06), the system has a large clearance threshold δ and the excitation frequency range is ω ∈ (0.5, 1.55), mass M 1 exhibits a large number of subharmonic periodic motions, such as 2-0-1, 2-1-0, 2-1-1, 2-1-2, 2-2-1, 2-2-2, 3-1-1, 3-2-2, 3-1-2, 3-1-3, 4-2-3, 4-3-3, and 4-4-2 motions, and mass M 2 undergoes 2-0-0,     motions, which results in complex dynamic characteristics of the system as a whole. e calculation results show that increasing the damping coefficient ζ can restrain the typical nonsmooth characteristics, such as grazing bifurcations, catastrophes, and hysteresis of the system's periodic impact vibration.

Conclusions
In this study, a model of a mechanical vibration system with multiple rigid stops is established. Using multiparameter and multiobjective cosimulation analysis, the influences of the model parameters and the excitation frequency on dynamic characteristics of the mechanical vibration system are revealed: (1) In the parameter domain of low frequency and small clearance, the irreversibility of the dynamic evolution of the adjacent basic periodic impact motion results in the occurrence of hysteresis transition domain, nonhysteresis transition domains, and singular points. Given the different bifurcation characteristics, two adjacent basic periodic motions coexist in the hysteresis transition domain, and the nonhysteresis transition domain contains n − (np + 1) − np and 2n − 2(np + 1) − 2np subharmonic motions. Owing to rigid stops, the system exhibits complex bifurcation characteristics, such as codimension-twobifurcations (flip-saddle nodebifurcation, doubling grazing bifurcation, and doubling saddle node bifurcation) and sliding bifurcation; the former are caused by the irreversibility of the transition between adjacent basic periodic impact motions, and, the latter leads to chattering-impact sticking events. erefore, when optimizing the functional objectives of a mechanical system, it is necessary to avoid the parameter domains that induce complex dynamic characteristics.
(2) e dimensionless mass ratio has a significant influence on the pattern types and bifurcation characteristics of the periodic impact vibration of the system. e results show that the smaller the mass ratio is, the larger the parameter domains of the chattering-impact vibration and complete chattering-impact vibration of the system are. e larger the mass ratio is, the more diverse the basic periodic impact motion groups and the more complex the dynamic characteristics of the system are in the parameter domain related to low frequency and small clearance threshold. e large parameter domains of the dimensionless mass ratio should be avoided as much as possible in the optimization design of the mechanical system. e basic periodic impact motions of system mass M 1 exhibit a competitive excitation with the increase in the dimensionless damping ratio. Although mass M 2 does not induce new periodic impact motions, the parameter domain of the basic periodic impact motions increases toward the highfrequency domain. e small damping ratio makes the basic periodic impact vibration groups, and pattern diversity of subharmonic impact motions become more significant; that is to say, smaller damping leads to an increase in the nonhysteresis transition domain. As a result, the damping ratio is larger and the nonlinear vibration of the system in low frequency is simpler. Large damping ratio is beneficial to suppress subharmonic impact motion and basic periodic motions with multiple impacts in the low-frequency domain, and the energy consumption of the system increases accordingly.
(3) e stiffness ratio has a significant effect on the mode type, distribution domain, and bifurcation characteristics of the various basic periodic impact groups  and subharmonic motions of the system in the (ω, δ)-parameter plane. e calculation results show that the stiffness ratio in the range of μ k ∈ [0.25, 0.75] is most beneficial to the optimum design of mechanical system. With the increase in the dimensionless parameter damping coefficient, the complex periodic motions of the system, such as subharmonic motions and chaos, are significantly reduced, and the parameter domains of the hysteresis transition domain and nonhysteresis transition domain are reduced accordingly. Furthermore, the existing domains of the mass chattering-impact vibration and the complete chattering-impact vibration shift in the high-frequency direction and the typical nonsmooth characteristics of the system are effectively suppressed. erefore, the dimensionless damping coefficient should be set slightly larger when optimizing the design of a mechanical system. e numerical results reveal the essential relationships and matching behaviours between the dynamic characteristics and model parameters of a vibrating system with multiple rigid stops. e reasonable matching ranges of the parameters of the mechanical vibrating systems with different functions can be determined, and the collaborative optimization of their dynamic characteristics and functional objectives can be achieved.

Data Availability
Relevant research data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.