Vibrating systems or devices with clearance or gaps between moving parts are frequently encountered in the aeronautical, civil, and mechanical engineering fields, particularly in mechanisms and machines,^[1,2] such as ships colliding with fenders,^[3] rattling of gear boxes,^[4] structural interaction in piping systems,^[5] and nonlinear soil-structure interaction.^[6] This type of systems or devices with repeated impact is called vibro-impact systems or impacting oscillators. Due to the existence of the impact, a vibro-impact system presents strong nonlinearity and discontinuity, accompanied with rich and complicated dynamic behaviors. Therefore, the analyses of the impact motion are important in the proper design of the corresponding machines and devices.^[7] Notable contributions were made by many researchers,^[8–25] and some interesting phenomena, such as grazing bifurcation,^[9] codimension two bifurcation,^[11] and torus bifurcation,^[14] chatter and sticking motions were observed in the vibro-impact systems.

However, vibro-impact systems, as one kind of non-smooth systems, always have troubles due to their inherent difficulty of analytical study. The main difficulty depends on the additional finite relations between impacts and rebound velocities. On one hand, the impact instants could not predict in advance, but are governed by the equations of motion. Thus, it is reasonable that the impact instants are described by constraints which depend on the state of systems. On the other hand, the impact instants govern the corresponding velocity jumps, i.e., the differential equations of motion on the impact instants are complemented by impact conditions. The previous literatures, which were related to the frequency response of vibro-impact systems, mainly concentrated on traditional frequency responses with single-valued response^[16,18] or three-valued response (with jump). However, for more complicated multi-valued responses,^[24] such as two-valued and four-valued responses, especially the frequency island phenomenon in the vibro-impact systems, which are not mentioned. For this purpose, we consider the same oscillators as that in Ref. [21], and estimate the multi-valued responses and dynamic properties of the system.

This paper is organized as follows. In Section 2, by using the non-smooth transformation the vibro-impact system is transformed a non-impact system. In Section 3, the Krylov-Bogoliubov averaging method is used to obtain the frequency response equation and the backbone curve equation. By using the linearization method, the stability conditions are derived. In Section 4, the directly numerical simulation is carried out to verify the analytical results obtained in Section 3, the effects of system parameters on the dynamic behaviors are revealed. The classification of the steady-state solutions is discussed. Finally, conclusions are presented in Section 5.

A general nonlinear vibro-impact oscillator with a unilateral nonzero offset barrier is considered, and the schematic model of this oscillator is shown in Fig. 1, where x is the displacement of the mass, F(t) is the external excitation, and h is the distance between the oscillator’s static equilibrium position and the rigid barrier. In the following discussion, the non-dimensional equations of the motion are introduced,^[21] which are presented as follows:

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Schematic model of the vibro-impact oscillator with a unilateral rigid barrier.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Phase portraits of system (1) for Ω = 2.08: (a)c1 = 0.1, c3 = −0.1; (b) c1 = 0.1, c3 = 0.0; (c) c1 = 0.1, c3 = 0.1; (d) c1 = 0.2, c3 = −0.1; (e) c1 = 0.2, c3 = 0.0; (f) c1 = 0.2, c3 = 0.1.

By using the non-smooth transformation^[25]

According to Eq. (2), equation (1b) is converted to

Thus, the jump of converted velocity ẏ becomes proportional 1 − r instead of 1 + r for the jump of original velocity ẋ. It is reasonable to use the Dirac function to introduce this jump into the equation of motion as an additional impulsive damping term.

Since y(t_* ±0) = y(t_*) = 0 and ẏ_± = ẏ(t_* ± 0), by introducing the Dirac function δ (y(t_*)) = δ (t − t_*)/|ẏ| and further combining with ẏ(t)δ (t − t_*) = ẏ(t_*)δ (t − t_*), the impulsive term can be obtained

According to Eqs. (3) and (4), equation (1) can be changed into

Therefore, the original system with impact is reduced to one without impact expressed by Eq. (5).

For the necessity of the Krylov-Bogoliubov averaging method, we may assume that the parameter c₁, c₂, c₃, h, and 1 − r are all small and proportional to a small parameter, then the asymptotic method of averaging to one period can be used to analyze Eq. (5). Considering the subharmonic response, then the frequency Ω of external excitation is near the resonance frequency 2ω₀(Ω ≈ 2hω₀),

In order to use the averaging method, the following transformation is used to Eq. (5):

The right sides of Eqs. (8) and (9) are proportional to a small parameter, the derivatives a′ and φ′ are both small under the assumption that c₁, c₂, c₃, h, and 1 − r are all small, then a and φ are slowly varying random processes with respect to time t, and Φ is a fast varying random progress. The average procedure is applied to this system over the fast state variables Φ, and then the following equations can be derived:

The steady-state solutions of Eqs. (10) and (11) are considered. Leta = a₀, φ = φ₀, a′ = 0, and φ′ = 0, then we can obtain the equations that the steady-state solutions of Eqs. (10) and (11) satisfy as follows:

According to Eqs. (12) and (13), the frequency response equation can be derived as

For given parameters of system (1), equation (14) can be solved by a numerical method. We introduce the perturbation terms which can be used to examine the stability of the steady-state solutions,

Substituting Eq. (16) into Eqs. (10) and (11) and neglecting the nonlinear terms, we obtain the linearization modulation of Eqs. (10) and (11) as follows:

According to the characteristic equation of the coefficient matrix of Eqs. (17) and (18), the necessary and sufficient condition of the stability of the steady-state solutions a₀ and φ₀ is that the real parts of the eigenvalues of the coefficient matrix are less than zero. Due to Routh–Hurwitz rule,^[26] it needs to satisfy

Particularly, the stability and classification of the steady-state solutions determined by the sign of Σ₁, Σ₂, and can be found in Ref. [27].

In this section, by combining the numerical simulation, the analytical results can be verified, the frequency response is constructed by using the fourth-order Runge–Kutta algorithm. For the non-smooth characteristic of the impact system, two types of step size to improve the computation precision are chosen. Firstly, the large time step size t₀can be chosen. Once the displacement x is less than or equal to −h, we pull back the oscillator to the former step position and integrate Eq. (1) by a small step size of 0.01t₀. When the displacement x ≤ −h is satisfied under the small step size of 0.01t₀, the velocity transformation (1b) can be used to derive the new velocity and the integration step size is then recovered to t₀. In addition, from Eqs. (2) and (7), we know that a² = y² + (ẏ/ω₀)² = (x + h)² + (ẋ/ω₀)². All results are obtained for the first-order subharmonic response (n = 1), with mean excitation frequency being close to 2ω₀.

The parameters of the system are listed as follows: ω₀ = 1.0, r = 0.8, h = 0.08, c₂ = 0.15, and f = 0.3. Three cases of the nonlinearity intensity (c₃ = −0.1, c₃ = 0.0, and c₃ = 0.1) are considered, respectively. In order to confirm the existence of impact in system (1), the phase portraits of system (1) are depicted corresponding to two different cases of damping coefficient c₁, it can be seen that under these cases the response of the system should be the impact ones. In this paper, only the impact case is considered.

4.1. Response characteristic for different damping parameters

In the first case, damping coefficient c₁ = 0.1 is considered. In Fig. 3, steady-state response curves defined by Eq. (14) are plotted for different values of excitation amplitude corresponding to three different cases of c₃. Obviously, when f is smaller (in Fig. 3 for f = 0.1), the system response is double-valued: one is stable while the other is unstable. This kind of response only occurs in the lower frequency region, since in the higher frequency region equation (14) only has complex solutions, which means that the system does not have steady-state solutions. In Table 1, the maximum frequencies for double-valued responses are presented corresponding to three different values of nonlinear intensity c₃ as f = 0.1. However, it is worth noting that when excitation amplitude f is larger (in Fig. 3 for f = 0.3), the system response is the traditional frequency response and with jump phenomenon in the case of c₃ = −0.1 and c₃ = 0.1. There exists a critical value f = 0.24 between the traditional frequency response and the double-valued response. As shown in Fig. 3, when f = 0.24 there is a branch of the response curve is close to zero. For the softening type (in Fig. 3(b) for c₃ = −0.1), this branch only appears in the higher frequency branch (Ω > 2); for the hardening type (in Fig. 3(c) for c₃ = 0.1), it exists in the lower frequency branch (Ω < 2); however, for the linear type (in Fig. 3(a) for c₃ = 0.0), it distributes in the whole frequency axis. Additionally, the steady-state numerical results confirm a very good accuracy of the analytical ones.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Steady-state responses of system (1): (a) c1 = 0.1, c3 = 0.0; (b) c1 = 0.1, c3 = −0.1; (c) c1 = 0.1, c3 = 0.1. Solid lines: analytical results; the marks of ◼, *, and •: numerical results; dashed line: stability boundaries; and dash-dotted lines: Backbone curves.

Table 1.

Table 1.

Table 1. Critical values of excitation frequency Ω for the boundary of two steady-state solutions. . c3 –0.1 0.0 0.1 Ω 1.888 1.932 2.008

Table 1. Critical values of excitation frequency Ω for the boundary of two steady-state solutions. .

The effectiveness discussed above is also verified by the stability boundaries, the dash lines as shown in Fig. 3, which are derived by the inequality (19). In the case of c₁ = 0.1, Σ₁ is negative in the whole Ω– a₀ plane. Hence, only the steady-state solutions in the interior regions of the dash lines, namely, the shaded regions Σ₂ < 0 are unattainable. Moreover, the backbone curve given by Eq. (15) is also plotted in Fig. 3 by the dash-dotted lines, where the change of its shape and position can be seen.

As a comparison, figure 4 depicts the response characteristic of the vibro-impact system when damping coefficient c₁ = 0.2. It can be seen that when the excitation amplitude is smaller (such as f = 0.05), the frequency response curve is separated into two parts, which is the frequency island phenomenon. Actually, this also leads to the occurrence of the four-valued response as shown in Fig. 4(d) (the local amplification of Fig. 4(b)) in the interval of Ω ∈ [1.745,1.833] for f = 0.05; when Ω is smaller than 1.745, the response keeps two-valued. As indicated in Table 2, the response characteristic for f = 0.05 is denoted as 2S–4S–2S–0S (nS denotes that the numbers of steady-state solutions are n). It can be seen that in the case of f = 0.2, there also exists a smaller interval Ω ∈ [1.642,1.763] (Fig. 4(d)), in which the response is four-valued. However, for the linear like type (c₃ = 0.0) and the hardening type (c₃ = 0.1), there only exists the two-valued response for f = 0.05 as well as f = 0.2, as shown in Figs. 4(a) and 4(c).

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Steady-state responses of system (1): (a) c1 = 0.2, c3 = 0.0; (b) c1 = 0.2, c3 = −0.1; (c) c1 = 0.2, c3 = 0.1. Solid lines: analytical results; the marks of ◼, *, and •: numerical results; dashed line: stability boundaries; and dash-dotted lines: Backbone curves.

Table 2.

Table 2.

Table 2. Multi-valued response characteristic of the system for c1 = 0.2. . c3 f Interval of excitation frequency Ω for corresponding numbers of the steady-state solutions Response characteristic one or two zero or two or three or four one or two zero −0.1 0.05 [0.0,1.745] [1.745,1.833] [1.833,1.842] [1.842, +∞) 2S–4S–2S–0S 0.2 [0.0,1.642] [1.642,1.763] [1.763,1.954] [1.954, +∞) 2S–4S–2S–0S 0.3 [0.0,1.584] [1.584,1.725] [1.725, +∞) 1S–3S–1S 0.0 0.05 [0.0,1.856] [1.856,1.904] [1.904,1.949] [1.949, +∞) 2S–0S–2S–0S 0.2 [0.0,1.990] [1.990, +∞) 2S–0S 0.3 [0.0, +∞) 1S 0.1 0.05 [0.0,1.882] [1.882,1.975] [1.975,2.137] [2.137, +∞) 2S-0S–2S–0S 0.2 [0.0,2.264] [2.264, +∞) 2S–0S 0.3 [0.0,2.057] [2.057,2.329] [2.329, +∞) 1S–3S–1S

Table 2. Multi-valued response characteristic of the system for c1 = 0.2. .

As a result, the stability characteristic is also very different from the case of c₁ = 0.1. Only the steady-state solutions above the dash line Σ₁ = 0 and outside the region Σ₂ < 0 are attainable, which is similar to the results derived in Ref. [28]. As can be seen that for the frequency island response only the branch with the largest amplitude is stable and the steady-state solutions on the mainland are unstable, which is different from the case shown in Ref. [29]. These findings are checked numerically by carrying out direct integration of system (1) and a good agreement can be found between the numerical results and the analytical ones.

4.2. Analysis of the classification of the steady-state solutions in Section 4.1

As discussed in Section 3, the responses are quite different for different damping coefficient c₁ and nonlinear stiffness c₃; therefore, the classification of the steady-state solutions is also very different. In the following part, it can be analyzed specifically.

The stability classification^[27] is illustrated in Fig. 5, where the critical stability boundaries divide the plane into stable and unstable regions. Comparing Figs. 5(a1), 5(b1), and 5(c1) with Figs. 5(a2), 5(b2), and 5(c2) indicates that the damping coefficient can cause a drastic change in the response characteristic. Through the boundary Σ₂ = 0, the steady-state solutions lose stability via a saddle node bifurcation, and they are saddle points interior the boundary Σ₂ = 0. Above the dash line Σ₁ = 0 as well as outside the saddle points region, the steady-state solutions are stable, which include stable nodes and stable foci; and the critical boundary is the division line for nodes and foci. However, the critical boundary Σ₁ corresponding to the occurrence of Hopf bifurcation keeps unchanged in spite of different values of nonlinearity intensity c₃. The non-trivial solutions lose stability via a Hopf bifurcation yielding aperiodic response.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Classification of the steady-state solutions: (a1) c1 = 0.1, c3 = 0.0; (b1) c1 = 0.1, c3 = −0.1; (c1) c1 = 0.1, c3 = 0.1; (a2) c1 = 0.2, c3 = 0.0; (b2) c1 = 0.2, c3 = −0.1; (c2) c1 = 0.2, c3 = 0.1.

4.3. Response analysis according to the jump avoidance condition

In previous sections, the existence of the saddle node bifurcation which results in the jump phenomenon is found. A jump will occur as long as there are two points on the steady-state response curve with a vertical slope,^[30–32] corresponding to dΩ/da₀ = 0. The jump can be avoided if it can be ensured that the frequency response of the system is unique. According to Eq. (14), we have

where

Then, the excitation frequency can be derived as follows:

Differentiating both sides of Eq. (21) with respect to a₀and substituting dΩ/da₀ = 0 yield

where

According to Eq. (22), figure 6 depicts the critical amplitude a₀ as a function of excitation amplitude f for c₁ = 0.1 and c₁ = 0.2, respectively. These results also verify the characteristic shown in Figs. 2 and 3. It can be seen that for c₃ = 0.0 and c₁ = 0.1 (Fig. 6(a1)), as f is smaller than 0.24, there is one point satisfying dΩ/da₀ = 0, which is due to the two-valued response characteristic. Then, as f is greater than 0.24, there is no critical amplitude (Fig. 6(a1)), and therefore no jump. Hence, f = 0.24 is the maximum and acceptable value of f in practical design to avoid the jump for c₃ = 0.0 and c₁ = 0.1.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Response analysis according to jump avoidance condition: (a1) c1 = 0.1, c3 = 0.0; (b1) c1 = 0.1, c3 = −0.1; (c)c1 = 0.1, c3 = 0.1; (a2) c1 = 0.2, c3 = 0.0; (b2) c1 = 0.2, c3 = −0.1; (c2) c1 = 0.2, c3 = 0.1.

Next, as shown in Fig. 6(a2) for c₃ = 0.0 and c₁ = 0.2, in the interval of f ∈ [0,0.054], there are three points satisfying dΩ/da₀ = 0, which is induced by the two-valued response with frequency island phenomenon. Then, in the interval of f ∈ [0.054,0.24], there is one point satisfying dΩ/da₀ = 0 corresponding to the complete two-valued response without frequency island. However, in the interval of f ∈ [0.24,0.26], there still exist two points satisfying dΩ/da₀ = 0. This is due to the existence of traditional frequency response with jump phenomenon, such as the example shown in Fig. 7(a) for f = 0.25. When f is greater than 0.26, there is no critical amplitude for c₃ = 0.0 and c₁ = 0.2 (Fig. 6(a2)), and therefore no jump.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Frequency response of system (1) for verifying the case in Fig. 6.

For c₃ = −0.1 and c₁ = 0.1, as shown in Fig. 6(b1), there exist two branches, the points in the lower amplitude correspond to the two-valued response, while the branch in the higher amplitude is relative to the traditional frequency response with jump phenomenon. Particularly, in the interval f ∈ [0.152,0.24], there exist three points satisfying dΩ/da₀ = 0 for every value of f. Actually, corresponding to certain value of f, the response is four-valued, such as the example shown in Fig. 7(b) for f = 0.2. For c₃ = −0.1 and c₁ = 0.2, as shown in Fig. 6(b2), when f is smaller than 0.24 there always exists the four-valued response. Therefore, it has three points satisfying dΩ/da₀ = 0. When f is greater than 0.24, there exist two points satisfying dΩ/da₀ = 0, which corresponds to the response with jump phenomenon.

For c₃ = 0.1 and c₁ = 0.1, as shown in Fig. 6(c1), in the interval f ∈ [0.0,0.24], it is pure two-valued response. For c₃ = 0.1 and c₁ = 0.2, as shown in Fig. 6(c2), in the interval f ∈ [0.0,0.055], the response is two-valued with frequency island phenomenon (Fig. 4(c)). Thus, there exist three points which satisfy dΩ/da₀ = 0. In the interval f ∈ [0.055,0.24], it is pure two-valued response. Both in Fig. 6(c1) and Fig. 6(c2), f is greater than 0.24, which is a traditional frequency response.

4.4. Effects of the distance between the system’s static equilibrium position and the barrier on the frequency response

When the excitation frequency is fixed at Ω = 2.08, the distance h is varied. The steady-state response of a₀ is shown in Fig. 8.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Effects of the distance h on the steady-state response: (a) c3 = −0.1; (b) c3 = 0.1.

As can be seen that, when the nonlinear stiffness is of softening type, with the increase of the distance h, the steady-state response from the maximum value keeps decreasing to zero in spite of different values of excitation amplitude. However, when the nonlinear stiffness is of hardening type, with larger excitation amplitude, three-valued and two-valued responses can be induced. Obviously, when f = 0.1, with the increase of h there exists a small interval in which the response is three-valued. For f = 0.2 and f = 0.3, as described in Table 3, three-valued response and two-valued response coexist. Hence, the distance is also a very important factor to influence the response of the vibro-impact system.

Table 3.

Table 3.

Table 3. Multi-valued response characteristic of the system with hardening type of stiffness. . f The effective response interval of the distance h The response characteristic One Three Two or one   0.05 [0.0,0.0168] 1S 0.1 [0.0,0.0063] [0.0063,0.0188] [0.0188,0.0335] 1S–3S–1S 0.2 [0.0,0.0382] [0.0382,0.0667] [0.0667,0.1061] 1S–3S–2S 0.3 [0.0,0.0696] [0.0696,0.1000] [0.1000,0.2048] 1S–3S–2S

Table 3. Multi-valued response characteristic of the system with hardening type of stiffness. .

The multi-valued responses and dynamic stability of a nonlinear vibro-impact system with a unilateral non-zero offset barrier are analyzed in this paper. Using the Zhuravlev non-smooth transformation and Krylov–Bogoliubov averaging method, the frequency response, the critical equation for the stability boundary are derived. Meanwhile, the backbone curves are also obtained. Then, by combining numerical simulation, we further investigate the effects of the parameters on the responses of the vibro-impact systems. It is found that multi-valued response, especially two-valued and four-valued responses, exist in the vibro-impact system. Meanwhile, the damping coefficient can induce the frequency island phenomenon, which is verified by the direct simulation results. In addition, the condition for suppressing the jump phenomenon is derived, and by using this condition the characteristic of multi-valued response is also confirmed. The classification of the steady-state solutions is discussed. Finally, the distance between the system’s static equilibrium position and the barrier can induce two-valued and three-valued responses when the nonlinear stiffness is of hardening type.