Bifurcation Characteristics of Fundamental and Subharmonic Impact Motions of a Mechanical Vibration System with Motion Limiting Constraints on a Two-Parameter Plane

A two-degree-of-freedom periodically forced system with multiple gaps and rigid constraints is studied. Multiple types of impact vibrations occur at each rigid constraint and interact with each other, which results in the emergence of some complex transitions in the system. Through the cosimulation of the key parameters gap value δ between the two masses and the excitation force frequency ω , the types, existence areas, and bifurcation regularities of the periodic and subharmonic motions can be obtained on the ( ω , δ )-parameter plane. In the corresponding three-dimensional surface diagram of the maximum impact velocity, the distribution law of the maximum impact velocity at each constraint can be obtained. The transition laws of fundamental impact motions in the low-frequency parameter domain are studied, and two types of transition regions in the transitions of adjacent fundamental impact motions are found: tongue-like regions and hysteresis regions. Moreover, these two types of transition regions show some atypical partitioning and deformation due to the combined eﬀects of impact vibrations at each constraint. By combining the two-parameter plane diagram and the three-dimensional surface diagram, the eﬀect of changing the gap values between each mass and the ﬁxed constraint and the damping coeﬃcient ζ on the dynamic characteristics of the system is studied. Combining the existence areas of periodic motions and the distribution of maximum impact velocity can provide guidance for the reasonable selection of system parameters.


Introduction
In the mechanical system, there are inevitably gaps and constraints between components and between the components and the fixed boundaries due to the needs of installation and the limitations of machining accuracy. e existence of gaps and constraints makes the mechanical system witness impact vibration when working, which will increase the noise level of the mechanical system and accelerate the wear of components. Due to the existence of constraints and gaps, the system becomes a typical nonsmooth system showing rich dynamic characteristics.
People want to reveal the inherent laws of impact vibration systems and use this to guide the effective use of impact vibration characteristics and regularize their potential hazards. A large number of researches in this field can be found in available monographs. Shaw and Holmes [1] considered a one-degree-of-freedom nonlinear oscillator and found period-doubling bifurcation, through which the system can lead to chaos. Aidanpää and Gupta [2] found the characteristics of period doubling, jumping, chaos, and so forth in the two-degree-of-freedom forced vibration system. Peterka and Vacik [3] discussed the law of the transition from periodic motion to chaos and characterized the periodicimpact motion by p/n. As a strong nonlinear factor, dry friction often exists in shock vibration systems, and the coupling of such strong nonlinear factors as gap, constraint, and dry friction was considered in [4][5][6][7]. Multi-degree-offreedom piecewise smooth mechanical system is a multiparameter high-dimensional nonlinear system; the codimension-two bifurcations have been studied in [8][9][10][11][12][13]. Some experiments were also used to verify the results of numerical simulation of impact vibration systems in [14][15][16].
Due to the existence of constraints and gaps, the impact vibration system has bifurcation types that are not available in smooth systems, such as grazing bifurcation. Nordmark [17] studied the singularities caused by grazing impact and found that the motion after the grazing bifurcation may be nonperiodic. Whiston [18] used the singularity theory to study the shredding mechanism of the globally stable manifold caused by the discontinuity of the impact map. Foale and Bishop [19,20] defined different types of grazing bifurcation. e discontinuous geometric method of impact mapping was used to further study the singularity and grazing bifurcation of piecewise smooth dynamic systems, as reported in [13,[21][22][23][24][25][26][27].
Complex chattering behavior has appeared in many impact vibration systems, which has attracted widespread attention from scholars. Chattering is a phenomenon in which there are an infinite number of impacts within one excitation force cycle [28]. Toulemonde and Gontier [29] analyzed the phenomenon of sticking periodic motions of impact oscillators and identified rising bifurcation. Wagg [30] discussed the chattering and sticking behavior that occurs under low-frequency excitation force conditions. Nordmark and Piiroinen [31] introduced the local discontinuous map based on the impact chatting accumulation and used the map to analyze the stability of the limit cycle with chatting.
In the previous studies of nonlinear dynamics, the periodic motion mode and bifurcation characteristics of the system were mainly studied through single-parameter bifurcation diagrams. en, most dynamic systems belong to multiparameter systems, and many parameters affect the dynamic characteristics of the system together.
erefore, how to get the existence areas of the periodic and subharmonic motions and the bifurcation characteristics through the double-parameter cosimulation has become a hot spot for scholars. In [12,22,32,33], the double-parameter bifurcations of one-degree-of-freedom impact vibration system were studied. Luo [34,35] studied the mode types and distribution areas of the periodic motions of the two-degree-of-freedom impact vibration system on the double-parameter plane. e aim of this paper is to get the types, existence areas, and bifurcation regularities of the periodic and subharmonic motions of the two-degree-of-freedom impact vibration system with multiple constraints and gaps through double-parameter cosimulation. Discover the reasonable matching relationship between the vibroimpact dynamics of the system and the parameters. e paper is organized as follows. In Section 2, the physical model of a two-degree-of-freedom periodically forced system with multiple gaps and multiple rigid constraints is introduced. In Section 3, through the double-parameter numerical simulation, the types, existence areas, and bifurcation regularities of periodic motions are obtained on the (ω, δ)-parameter plane, and the distribution of the maximum impact velocity at each constraint is studied. In Section 4, the effects of changing the gaps δ 1 � δ 2 between the two masses and the foundation supports and the damping coefficient ζ on the system dynamics are studied. Conclusions are obtained in Section 5.

Mechanical Model
A two-degree-of-freedom periodically forced system with multiple gaps and rigid constraints is depicted in Figure 1, where the masses M i (i � 1, 2) displacements are represented by X i (i � 1, 2). e masses M i (i � 1, 2) are connected to the supporting foundation by linear springs K i (i � 1, 2) and linear dampers C i (i � 1, 2). e masses M i (i � 1, 2) are excited by harmonic force of amplitude P i (i � 1, 2), frequency Ω, and phase angle τ. e clearance between the mass M 1 and the impact surface A 1 is B 1 . e clearance between the mass M 2 and the impact surface A 2 is B 2 . e clearance between M 1 and M 2 is B 3 . e stop is assumed to be rigid and the coefficient of recovery after a rigid collision is R. e speed before and after the impact of mass M i (i � 1, 2) is represented by _ X i− and _ X i+ (i � 1, 2), respectively. In the study, we introduced the following dimensionless parameters: It can be found that some parameters in the system have a certain value range through dimensionless parameters, such as μ m ∈ (0, 1), μ k ∈ (0, 1), μ c ∈ (0, 1), f ∈ [0, 1]. When the system has no impact on each rigid constrains, the motion of the impact system could be described by the nondimensional equations: When the mass M i (i � 1, 2) has an impact with the fixed constraints, the velocity of mass M i (i � 1, 2) is changed immediately.

2
Shock and Vibration e instantaneous velocity of mass M 1 and mass M 2 before and after collision can be determined by momentum conservation.
where _ x i− and _ x i+ (i � 1, 2) represent the velocities before and after the impact, respectively.
For the periodically forced impact system, the mass M i (i � 1, 2) may stay on the constraint surfaces with a speed of zero for a period of time under the resultant force and this phenomenon is called sticking. e impact vibration system in Figure 1 contains multiple gaps and rigid constraints, so the system may exhibit sticking phenomenon on multiple rigid constraint surfaces in a certain motion time section. From the conditions of the occurrence of the sticking phenomenon, the sticking conditions of the mass M i (i � 1, 2) on each rigid constraint surfaces can be analyzed: (1) e conditions for mass M 1 sticking to the surface of constraint A 1 are x 1 � −δ 1 and _ x 1 � 0, and the resultant force of the harmonic excitation force, spring, and damping force on the mass M 1 is satisfied: Until F 1 (t) changes direction, the sticking motion is ended. (2) e conditions for mass M 2 sticking to the surface of constraint A 2 are x 2 � δ 2 and _ x 2 � 0, and the resultant force of the harmonic excitation force, spring, and damping force on the mass M 2 is satisfied: Until F 2 (t) changes direction, the sticking motion is ended. (3) e conditions for mass M 1 and mass M 2 sticking to each other at constraint A 12 are Whether there is sticking between the two masses depends on the interaction force between the two masses; here, the interaction force is set to N(t). At this time, the resultant forces of mass M 1 and mass M 2 can be expressed, respectively, as According to f 1 (t) � (1 − μ m /μ m )f 2 (t), N(t) can be obtained: When N(t) > 0, the two masses stick together and move synchronously. Until the interaction force disappears, that is, N(t) � 0, the sticking motion is ended.
Here, the symbol p/n is used to describe the types of periodic-impact motions at each rigid constraint, where n represents the ratio of the vibration period of the system to the period of the excitation force (n � 1, 2, 3, ...), and p represents the number of impacts at each rigid constraint (p � 0, 1, 2, 3, ...). When n � 1, that is, when the vibration period of the system and the excitation force period are equal, there are fundamental impact motions p/1 (p ≥ 1). When n ≥ 2 and p ≥ 1, there are subharmonic impact motions p/n (n ≥ 2, p ≥ 1).
e Poincaré section σ n is selected by selecting the minimum displacement of the mass M 1 in the period of the excitation force and it could determine the number. e Poincaré sections σ A1 , σ A12 , and σ A2 could determine the number of p on each of the rigid constraints, respectively. Figure 1: Physical model of a two-degree-of-freedom periodically forced system with multiple rigid constraints.
Generally, the impact mapping of the system corresponding to the sections σ A1 , σ A12 , and σ A2 can be expressed as Where, corresponding to sections σ A1 and σ A12 , , τ (i+1) ) T ; X ∈ R 4 and μ are the system parameters, where μ ∈ R m , m � 8.

Dynamical Characteristics of the Two-Degree-of-Freedom Periodically Forced System with Multiple Rigid Constraints
In this section, we mainly study the mode types, existence areas, and bifurcation regularities of the periodic motions on the two-parameter plane composed of the exciting force frequency ω and the gap value δ, as well as the corresponding distribution law of the maximum impact velocity at each constraint. In model 1, the gap can be divided into two types, one is the gap between the mass M i (i � 1, 2) and the fixed constraints and the other is the gap between M 1 and M 2 . Consider here that the gaps δ 1 and δ 2 are taken to have fixed values, and δ 3 is taken as a variable parameter and recorded as δ. Here μ m � 0.5, μ k � 0.5, μ c � 0.5, ζ � 0.1, R � 0.8, and δ 1 � δ 2 � 1.0 are selected as the basis parameters. e exciting force frequency ω and the gap δ between the two masses are set as bifurcation parameters. rough doubleparameter numerical simulation, the mode types and occurrence areas of periodic motions can be get on the (ω, δ)-parameter plane as depicted in Figure 2. e maximum impact velocity of the system on each constraint surface can be obtained in the (ω, δ, _ x max )-three-dimensional surface program as depicted in Figure 3. On the (ω, δ)-parameter plane, different types of periodic motions are represented by different colors and marked with p/n. CIS/1 represents the sticking motions. In addition, the unmarked dark gray areas are mainly quasi-periodic motions, chaotic motions, and unidentified periodic-impact motions (mainly including n or p is too large and existence areas are extremely small). Moreover, GB, G r , G b , SN, S, and PB represent the grazing, real grazing, bare grazing, saddlenode, sliding, and period doubling bifurcation, respectively. Figure 2 shows the existence types and distribution areas of the fundamental impact motions and subharmonic motions on the (ω, δ)-parameter plane associated with each constraint. e corresponding three-dimensional surface diagrams of the maximum impact velocity can be observed in Figure 3. As shown in Figure 2(a), on the (ω, δ)-parameter plane associated with constraint A 1 , in the frequency domain of ω ≤ 1.4 the system mainly exhibits 1/1 motion, and there is a small amount of 2/1 motion in the lower frequency domain. In the frequency domain of ω > 1.4 the system mainly exhibits impactless motion 0/1. Subharmonic motions mainly exist in the middle frequency domain, mainly 2/2, 1/3, 3/3, 4/4, and 4/ 4 motions. e corresponding three-dimensional surface diagrams of the maximum impact velocity are shown in Figure 3(a); the peaks of the maximum impact velocity mainly exist in the middle frequency domain where the motion form is mainly 1/1. In the high-frequency domain, the peaks of the maximum impact velocity also exist in the existence domain which mainly contains subharmonic motions 1/3 and 2/4. In the low-frequency domain, the maximum impact velocity is relatively small. As shown in Figure 2(c), on the (ω, δ)-parameter plane associated with the constraint A 2 , periodicimpact motions mainly exist in the middle frequency domain. e fundamental motion is mainly 1/1 motion, and the rest are mainly subharmonic motions, including 1/2, 1/2, 1/3, 2/3, and 2/4 motions. e corresponding three-dimensional surface diagrams of the maximum impact velocity are shown in Figure 3(c); the peak values of the maximum impact velocity exist in the regions of 1/1, 1/2, and 1/3 motions. Figure 2(b) shows the (ω, δ)-parameter plane associated with constraint A 12 , which presents richer dynamic characteristics on the entire parameter plane. In the domain of large δ and high ω, it is mainly impactless motion 0/1. e dark gray areas and subharmonic motions (mainly 3/3, 4/3, and 4/4 motions) mainly exist in the domain of middle ω. In the parameter domain located at the lower left corner of Figure 2(b), the types of periodic motions are more complicated. It mainly presented as the fundamental impact motions p/1 (p ≥1), the incomplete chatting motions p/1, and the sticking motions (p/1), as observed in Figure 4. e fundamental impact motions p/1 (p ≥1) present band-like distribution in the domain of low ω, and as the number of impacts p increases, the width of the band-like region gradually narrows. As the exciting force frequency ω decreases, p/1 motion will have a grazing bifurcation, and the number of impacts p increases by one, and it transitions to (p+1)/1 motion. When p keeps increasing, the system will enter incomplete chatting motions p/1. At this time, the number of impacts p is very large, but the minimum impact speed is still not zero, and the masses M 1 and M 2 do not stick together and move synchronously. When the number of impacts p increases sufficiently and the minimum impact velocity decreases to zero, the system will enter the sticking motions p/1. At this time, the interaction force between the two masses N(t) > 0, and the two masses are sticking together for synchronous movement. e three-dimensional surface diagrams of the maximum impact velocity corresponding to Figure 2(b) are shown in Figure 3(b); the fundamental impact motions regions near ω � 0.7 and ω � 1.2 have larger maximum impact velocity; in the middle frequency domain, the peaks of maximum impact velocity appear in the subharmonic motions 3/3 and 4/4 regions.
In order to further study the bifurcation characteristics of the system, a single-parameter bifurcation diagram that crosses the (ω, δ)-parameter plane along a horizontal scan can be taken. e global bifurcation diagrams of the system when the gap value δ � 0.1 is taken can be seen in Figure 5, and the relation between single-parameter bifurcation diagrams and periodic motions on the (ω, δ)-parameter plane can be observed. In the middle frequency domain, impact motions occur in all three constraints. For example, at the point (ω, δ) � (0.85, 0.1), the number of cycles n is two, the number of impacts p at constraint A 1 is two, the number of impacts p at constraint A 12 is six, and the number of impacts p at constraint A 2 is one, as shown in Figure 6. In the lowfrequency range, impact motions mainly occur at constraints A 1 and A 12 . e impact motions at constraint A 1 are mainly fundamental 1/1 and 2/1 motions, while the impact motions at constraint A 12 are more complicated, which include p/1 (p ≥1), (p/1), and (p/1) motions. Figure 7 shows a partial enlarged bifurcation diagram of the low-frequency part of Figure 5(b), in which the red part is the sticking motions. It can be seen that, with the decrease of ω, the 2/1 motion passes through the bare-grazing bifurcation point G b (2/1) (ω � 0.668) and enters the region containing the subharmonic motion 10/ 4 and then transfers to the 3/1 motion. e phase diagram at G b (2/1) is shown in Figure 8(a); it can be seen that the impact velocity at grazing bifurcation point is zero. e phase diagram at G r (3/1) (ω � 0.543) is shown in Figure 8(b), and 3/1 motion transfers to 4/1 motion. e phase diagram at G r Figure 8(c), and 4/1 motion transfers to 5/1 motion. So far, it can be found that, with the decrease of ω in the low-frequency region, the p/1 motion passes through real grazing bifurcation, and the number of impacts p increases once and directly enters (p + 1)/1 motion. When passing through bare-grazing bifurcation, the system will enter an unstable region that may include subharmonic motions, quasi-periodic motions, and chaotic motions, and then it will enter (p + 1)/1 motion. When ω is further reduced, the number of impacts p increases one by one, and the system enters (p/1). e time series of the system when ω � 0.3619 is shown in Figure 9. At this time, ω is very close to the sliding bifurcation point, and the system is in the incomplete chatting motion. It can be seen that the relative displacement   When ω continues to decrease, the system enters the sticking motions through the sliding bifurcation. Time series of sticking motions with ω � 0.35 are shown in Figure 10. It can be seen that, after the system has undergone a chatteringimpact sequence in which the relative displacement x 2 x · 2 (c) and the relative velocity _ x 1 − _ x 2 are gradually reduced, _ x 1 − _ x 2 is reduced to zero, and x 1 − x 2 is equal to the gap value δ � 0.1. e interaction force between the two masses N(t) > 0, and the two masses are sticking together and moving synchronously. In the low-frequency domain, as ω decreases, the transition law of the system on the (ω, δ)-parameter plane associated with constraint A 12 can be summarized as     Shock and Vibration 7 ω↓: From the (ω, δ)-parameter plane associated with constraint A 12 , it can be found that the fundamental impact motions p/1 (p ≥1) are in bands distribution. It can also be found that there are transition regions in the adjacent fundamental impact motions in Figures 2(b), 4, 11(a), and 11(b).
is transition regions look like the shape of the tongue, so these regions can be named as tongue-like regions, marked as TR (p/1)∩((p+1)/1) . When the numerical simulation was performed with decreasing ω and increasing ω, respectively, and depicted on the same (ω, δ)-parameter plane, another type of transition regions, the hysteresis regions, is found, marked as HR p . Due to the difference of initial values, there is coexistence of p/1 and (p+1)/1 in HR p , as shown in the dark gray areas in Figure 11(c).

The Correlation between Dynamics and System Parameters
is section mainly studies the relationship between system dynamics and system parameters. e effect of changing the fixed gap value δ 1 � δ 2 and damping coefficient ζ on dynamics is studied, respectively. In the following analysis, only the parameters that need to be studied are changed; other parameters are the same as the basis parameters.

e Correlation between Dynamic Performance and Fixed
Gap Value δ 1 � δ 2 . In the analysis of basis parameters in the previous chapter, the gaps δ 1 � δ 2 � 1.0 are considered. Here, the larger gap value δ 1 � δ 2 � 1.5 and the smaller gap value δ 1 � δ 2 � 0.5 are used for analysis.
When taking a larger gap value δ 1 � δ 2 � 1.5, Figure 17 shows the existence types and distribution areas of periodic motions on the (ω, δ)-parameter plane associated with each constraint. e corresponding three-dimensional surface diagrams of the maximum impact velocity can be observed in Figure 18. As for low ω, there is no impact vibration at constraint A 1 and constraint A 2 ; the impact vibration of the system is mainly represented as the impact between mass M 1 and mass M 2 at constraint A 12 , and its mode types are fundamental impact motions, the incomplete chatting motions, and sticking motions. Figure 19 is a partial enlarged view of Figure 17(b) in the low ω and small δ domain; impact vibration only occurs at constraint A 12 . e band-like region of fundamental p/1 motion is more complete, and the boundaries of tongue-like region TR p/1∩(p+1)/1 are also more obvious. For example, the upper boundary of TR 2/1∩3/1 is the curve G b 2/1 , and the lower boundary is the curve SN 3/1 ; the subharmonic motions in TR 2/1∩3/1 mainly appear as 5/2 motion. e upper boundary of TR 3/1∩4/1 is the curve G b 3/1 , and the lower boundary is the curve PD 4/1 ; the subharmonic motions in TR 3/1∩4/1 mainly appear as 7/2 motion. It can be found that the upper boundary of TR p/1∩(p+1)/1 is the curve G b p/1 , and the lower boundary is the curve SN (1+p)/1 or the possible PD (1+p)/1 , where the main form of subharmonic motions is (np + 1)/n (p �1, 2, . . .; n � 2, 3, . . .). On the (ω, δ)-parameter plane associated with constraint A 1 , impact vibrations occur only in the vertical-strip region in the middle frequency domain, which mainly presents 1/1 motion, as well as a small amount of subharmonic motions (mainly 2/2, 3/3, 4/4, and 5/5 motions); see Figure 17(b). On the (ω, δ)-parameter plane associated with constraint A 2 , there is a small amount of 1/1 motion in the middle frequency domain.
As shown in Figure 19, compared with the basis parameters, the maximum impact velocity at constraints A 1 and A 2 is reduced, and the maximum impact velocity at constraint A 12 is increased. At constraint A 1 , the peak of the maximum impact velocity appears near the vertical-strip region with ω � 1.0. At constraint A 12 , the peak of the maximum impact velocity exists in the domain of 1/1 motion. At constraint A 2 , the peak of the maximum impact velocity appears in the existence area of 1/1 motion and the existence area of subharmonic motion 1/2.
When taking a smaller gap value δ 1 � δ 2 � 0.5, Figure 20 shows the existence types and distribution areas of periodicimpact motions on the (ω, δ)-parameter plane associated with each constraint. As for low ω, there are impact vibrations at each constraint. At constraint A 1 , the system mainly presents fundamental impact motions, the incomplete chatting motions, and the sticking motions which are mainly distributed like several vertical strips. At constraint A 2 , the system mainly presented as 1/1 motion, as well as a small amount of 2/1 motion. On the (ω, δ)-parameter plane associated with constraint A 12 , due to the complex impact vibrations occurring at constraint A 1 , the band-like region of fundamental p/1 motion is deformed and segmented in lowfrequency domain. Figure 21 is a partial enlarged view of Figure 20(b) in the low ω and small δ domain, the band-like region of fundamental p/1 motion is disrupted, the existence areas and boundaries of the tongue-like regions are no longer clear, and more unrecognized gray areas are embedded in the domain of small δ. For example, the subharmonic motion 7/2 that appears during the transition of 3/ 1 motion and 4/1 motion appears as an interval, and its shape no longer appears as a tongue; the subharmonic motion 9/2 that appears during the transition of 4/1 motion and 5/1 motion also showed a similar distribution law. It can be seen that although the band-like region of fundamental p/ 1 motion is divided, the subharmonic motions between the adjacent p/1 and (p + 1)/1 motions still conform to the law of (np + 1)/n (p �1, 2, . . .; n � 2, 3, . . .). More unrecognized gray areas are embedded in the high-frequency domain of the system. e three-dimensional surface diagrams of the maximum impact velocity can be observed in Figure 22. e maximum impact velocity at constraint A 1 increases compared to the basis parameters, and its distribution area is roughly the same as that of the basis parameters. e peaks of the maximum impact velocity appear in the existence areas of 1/1 and 2/1 motions in the large δ domain. At constraint A 12 , the maximum impact velocity reduces compared to the basis parameters, and the peak appears in the subharmonic 4/2 motion region distributed in the domain of small δ and low ω. At constraint A 2 , the maximum impact velocity is increased compared with the basis parameters, and the peak appears in the fundamental 1/1 motion region distributed in the domain of middle ω.
By comparison with the basis parameters, for large δ 1 � δ 2 , the system has no impact vibration at constraint A 1 and constraint A 2 at low frequency, and impact vibration only occurs at constraint A 12 . On the (ω, δ)-parameter plane associated with constraint A 12 , the band-like regions of fundamental p/1 motions and the tongue-like regions TR p/1∩(p+1)/1 have complete boundaries. e maximum impact velocity increases at constraint A 12 and decreases at Shock and Vibration constraint A 1 and constraint A 2 . With the decrease of δ 1 � δ 2 , the impact vibration at constraint A 1 and constraint A 2 gradually intensified. At constraint A 1 , there is a transition from p/1 motion to p/1 motion. At constraint A 12 , the mode types and distribution areas of the periodicimpact vibrations have not changed much. Due to the interaction of the impact vibrations of each constraint, the band-like regions of fundamental p/1 motions are cut by the gray unidentified motions; the boundaries of the tongue-like regions TR p/1∩(p+1)/1 are also destroyed, but the laws of the main subharmonic motions mode types remain the same. e system embeds a large amount of gray unrecognized motions in the middle frequency domain, which reduces the existence areas of periodic motions. e distribution domains of the peaks of the maximum impact velocity have little change. e maximum impact velocity increases at constraints A 1 and A 2 and decreases at constraint A 12 .   e damping coefficient ζ is the precoefficient of the velocity term in the differential equation of the system, which affects the impact absorption characteristics of the dynamic system. e damping coefficient is ζ � 0.1 in the basis parameters; the larger value ζ � 0.2 and the smaller value ζ � 0.05 are selected to analyze its influence on the dynamics.
When a larger damping coefficient ζ � 0.2 is taken, the mode types and existence regions of the periodic motions on the (ω, δ)-parameter plane associated with each constraint are depicted in Figure 23.
e impact vibrations of the system only occur at constraint A 1 and constraint A 12 , and there is no impact vibration at constraint A 2 . Figure 23(a) is the (ω, δ)-parameter plane associated with constraint A 1 ; similar to the basis parameter, it is mainly 1/1 motion in the domain of ω ≤ 1.4, and the domains of 2/1, 2/2, 3/3, and 4/4 motions are significantly reduced.
ere is subharmonic motion 4/2 in the middle frequency domain. e parameter domain of ω > 1.4 is impactless motion 0/1. e (ω, δ)-parameter plane associated with constraint A 12 is shown in Figure 23(b), the existence domain of the fundamental impact motions p/1 (p ≥1) increases, and the embedded subharmonic motions and unrecognized gray areas decrease significantly. e existence domain of the fundamental impact motions p/1 (p ≥1) extends to the domain of large ω, its bulge forms a semicircular domain, and Figure 24 is a partial enlarged view of it. In the semicircular parameter domain, the transition process with decreasing of ω can be summarized as ω↓: Tongue-like regions appear between adjacent p/1 and (p + 1)/1 motions in the semicircular parameter domain. e tongue-like region TR 3/1∩3/1 contains the main forms of subharmonic motions as 2/5, 7/3, and 8/3 motions. e tongue-like region TR 4/1∩3/1 contains the main forms of subharmonic motions as 10/3 and 11/3 motions. It can be found that the main forms of subharmonic motions in the tongue-like regions TR p/1∩(p+1)/1 are (np + 1)/n and (np + 2)/ n (p �1, 2, . . .; n � 2, 3, . . .). e three-dimensional surface diagrams of the maximum impact velocity can be observed in Figure 25. e maximum impact velocities at constraint A 1 and constraint A 12 are both reduced compared to the basis parameters. At constraint A 1 , the peaks of the maximum impact velocity are similar to those in basis parameters, appearing in the domain of 1/1 motion. At constraint A 12 , the peaks of the maximum impact velocity appear in the region of 1/1 motion in the small δ domain and the region of subharmonic motion 4/4.
When a smaller damping coefficient ζ � 0.05 is taken, periodic motions on the (ω, δ)-parameter plane associated with each constraint are depicted in Figure 26. It can be seen that, on the (ω, δ)-parameter plane associated with each constraint, the types and distribution areas of periodic motions are roughly the same as those in the basis parameters, but the existence regions of subharmonic motions and the unidentified gray areas increase and are embedded to the domains of fundamental impact motions. e distribution and transition laws of the fundamental impact motions p/1 (p ≥1) on (ω, δ)-parameter plane associated with constraint A 12 are mainly analyzed. Figure 27 is a partial enlarged view of Figure 26(b) in the low ω and small δ domain; compared with the basis parameters, the wave-like undulations of the p/1 (p ≥1) motions band-like regions increase. e area of tongue-like regions increases, but, due to the embedding of the gray unrecognized area, the boundaries of the tongue-like regions are destroyed.
In three-dimensional surface diagrams of the maximum impact velocity shown in Figure 28, the maximum impact velocity associated with each constraint has increased compared to the basis parameters. At constraint A 1 , peaks of   By comparison with the basis parameters, changing the damping coefficient ζ has a greater influence on the mode types and existence areas of the periodic-impact motions. e large ζ makes the impact vibrations only occur in constraint A 1 and constraint A 12 . At constraint A 1 , the medium low ω domain is mainly 1/1 motion, and the medium high ω domain is mainly impactless motion 0/1. At constraint A 12 , the band-shaped regions of the fundamental    impact motions extend to the domain of large δ, and the unrecognized gray areas are sharply reduced. e area of the tongue-like regions decreased during the transition process, and its pattern types also changed. e maximum impact velocity at each constraint is smaller, and the number of peaks is reduced. e peaks only appear in the existence domain of 1/1 motion. When the damping coefficient ζ decreases, the impact vibration at constraint A 1 and constraint A 2 is enhanced, and the domain of the fundamental impact motions is embedded with a large number of subharmonic motions and gray unidentified areas. At constraint A 12 , the wave-like undulations of the p/1 (p ≥1) motions band-like regions increase, and the gray area divides the band-like regions, and its integrity is destroyed. e boundaries of the tongue-like regions are no longer complete, but the mode types of the subharmonic motions are consistent with the reference parameters. e maximum impact velocity at each constraint increases; the number of peaks increases, mainly distributed in the occurrence domains of 1/1 motion and subharmonic motions, as well as the unidentified gray areas.

Conclusions
In this paper, we consider the dynamic characteristics of a two-degree-of-freedom periodically forced system with multiple gaps and multiple rigid constraints. On the (ω, δ)-parameter plane associated with each constraint, the mode types, existence areas, and bifurcation regularities of the periodic motions are obtained. e distribution area of the maximum impact velocity and the correlation between the periodic motions distributions are acquired.
(1) Under the basis parameters, on the (ω, δ)-parameter planes associated with constraints A 1 and A 2 , the mode types of impact vibrations are relatively simple. On the (ω, δ)-parameter plane associated with constraint A 12 , the system exhibits rich dynamic characteristics. For low ω, the system presents the fundamental impact motions p/1 (p ≥1), the incomplete chatting motions p/1, and the sticking motions p/1. ere are two types of transition regions between adjacent fundamental impact motions p/1 and (p + 1)/1: tongue-like regions TR p/1∩(p+1)/1 and hysteresis regions HR p . TR p/1∩(p+1)/1 is bounded by the upper boundary G b p/1 and the lower boundary SN (1+p)/1 . Subharmonic motions dominate in TR p/1∩(p+1)/1 and their mode types are (np + 1)/n or (np + 2)/n (p �1, 2, . . .; n � 2, 3, . . .). HR p exists between the upper boundary SN (1+p)/1 and the lower boundary G r p/1 . Due to the different initial values, p/1 and (p + 1)/1 motions coexist in HR p . e peaks of the maximum impact velocity generally appear in the existence domains of 1/1 motion, partial subharmonic motions, and a small part of identified gray areas.
(2) e influence of changing the gap value δ 1 � δ 2 and the damping coefficient ζ on the system dynamics is separately considered. With the decrease of δ 1 � δ 2 , the impact vibrations at each constraint are strengthened and the mutual reaction is increased. At constraint A 1 , the low-frequency region appears as vertical strips of p/1 (p ≥1), p/1, and p/1 motions. On the (ω, δ)-parameter plane associated with constraint A 12 , the band-like domains of p/1 (p ≥1) are divided by the embedded gray areas, and the boundaries of the tongue-like regions TR p/1∩(p+1)/1 are partially destroyed. e maximum impact velocities at constraints A 1 and A 2 increase, and the maximum impact velocity at constraint A 12 decreases. e increase in the number of subharmonic impact vibrations and unidentified gray areas leads to an increase in the peak number of maximum impact velocities. e damping coefficient ζ has a great influence on the types and distribution areas of the periodic motions. With the decrease of the damping coefficient ζ, the domains of fundamental impact motions decrease, and the domains of subharmonic motions and the unrecognized gray areas increase sharply. e wave-like undulations of the p/1 (p ≥1) motions band-like regions increase, and the area of the tongue-like regions increases. At each constraint, the maximum impact velocity increases and the number of peaks increases.
Many mechanical vibration systems can be described by a two-degree-of-freedom periodically forced system with multiple gaps and rigid constraints. In incomplete chatting motions and sticking motions, since a lot of impacts occur in one exciting force cycle, it may accelerate the wear of components and increase noise level. e larger impact velocity mostly appears in the domains of 1/1 motion and partial subharmonic motions, which may aggravate the pitting fatigue in the mechanical system. erefore, multiparameter cosimulation can better provide theoretical guidance for the parameters design of impact vibration systems.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest.

Authors' Contributions
Shijun Wang contributed to conceptualization, formal analysis, writing original draft, writing review, and editing. Guanwei Luo supervised the study.