Near-resonant dynamics, period doubling and chaos of a 3-DOF vibro-impact system

A mechanical system composed of two weakly coupled oscillators under harmonic excitation is considered. Its main part is a vibro-impact unit composed of a linear oscillator with an internally colliding small block. This block is coupled with the secondary part being a damped linear oscillator. The mathematical model of the system has been presented in a non-dimensional form. The analytical studies are restricted to the case of a periodic steady-state motion with two symmetric impacts per cycle near 1:1 resonance. The multiple scales method combined with the sawtooth-function-based modelling of the non-smooth dynamics is employed. A conception of the stability analysis of the periodic motions suited for this theoretical approach is presented. The frequency–response curves and force–response curves with stable and unstable branches are determined, and the interplay between various model parameters is investigated. The theoretical predictions related to the motion amplitude and the range of stability of the periodic steady-state response are verified via a series of numerical experiments and computation of Lyapunov exponents. Finally, the limitations and extensibility of the approach are discussed.


Introduction
Vibro-impact systems can be found in many areas of science and engineering. Even with a simple structure, they are strongly nonlinear. Therefore, there are no general analytical solutions or strategies for dynamic problems of this type. Many approximate analytical methods designed for nonlinear systems are based on the classical perturbation approach, and they can be applied only to weakly nonlinear problems [18,24]. In recent decades, some analytical techniques have been developed to cope with vibro-impact systems by imposing certain conditions on their motion (e.g., periodicity at the resonant state): the power-law phenomenological modeling, the non-smooth temporal transformation (NSTT), the concept of impact modes, the Ivanov non-smooth coordinate transformation, and few others [12,27,28].
There is also a solution strategy that particularly fits in with the asymptotic approach to nonlinear problems and near-resonant dynamics. It consists in combining the multiple time scales method (MTSM) with a use of a sawtooth function which describes the non-smooth nature of the given system. In recent years the technique has been applied to the systems including vibroimpact nonlinear energy sink (VI NES) [9][10][11][13][14][15][16].
Problems of this kind are strictly connected with the increasingly extensive studies on targeted energy transfer (TET) and energy harvesting [32].
Approximate analytical solutions and the related stability analysis can be useful, among others, in the control and improvement of periodic orbits of vibroimpact systems [4]. However, in the framework of the 'MTSM-sawtooth function' approach, there is a lack of comprehensive works on stability studies of the periodic steady-state near-resonant solutions. Perhaps the most commonly cited work related to stability of impact dampers is the PhD thesis by Masri [19]. The proposed conception and results were also discussed in later papers [20,21]. Similar method was presented, for example, by Popplewell et al. [29]. Unfortunately, in all these works, motion of the primary system between impacts is described by the general solution for a piecewise linear oscillator with viscous damping and periodic excitation. Such a model is usually inconsistent with the approximate solution formulated with the method of multiple scales under the assumption of small forcing and/or damping.
Results of the stability analysis focused on an approximate analytical solution can be verified numerically, by computation of Lyapunov exponents associated with the given system. For nonlinear discontinuous problems (with impacts and/or dry friction), the standard, well-established technique for determination of Lyapunov spectrum [23,25,26] need to be modified in order to handle the discontinuities (jumps). A theoretically based and straightforward generalization of the method was proposed by Müller [22]. This approach has been successfully applied to various multi-degreeof-freedom systems with impacts (e.g., see [1][2][3]).
In practical applications, machine components that undergo impacts are often interconnected with other, non-colliding parts (e.g., in percussive machines/tools such as hammer drills). Mostly, these couplings are strong. The mechanical system under consideration consists of two coupled oscillators, of which one includes an internally colliding small body. Due to the assumption of weak coupling, the paper can be treated as a preliminary work towards more realistic and practical problems. The main aim of this paper is to adapt and discuss the conception of the stability studies suited for the 'MTSM -sawtooth function' approach.
The paper is divided into six sections. A mathematical formulation of the problem is presented in Sect. 2. Section 3 includes an outline of the approximate ana- Fig. 1 The analyzed vibro-impact system lytical procedure of solution, and admissibility of the steady-state solutions is discussed from the viewpoint of different criteria. Analytical results and their verification via numerical simulations are reported in Sect. 4. Section 5 includes some remarks on limitations and extensibility of the presented approach. Conclusions and final remarks are given in Sect. 6.

Physical and mathematical model
Let us consider the system illustrated schematically in Fig. 1. The primary part is a vibro-impact unit comprised of a damped linear oscillator (box) of mass m 1 and an embedded block of mass m 2 . The latter one can move freely between two rigid stops (L is the symmetric gap size). The unit is driven by a harmonic force F 1 , and the components undergo mutual impact interactions. The secondary part of the system, in turn, is an ordinary damped linear oscillator of mass m 3 , elastically coupled with the block. The stiffness constants associated with the linear springs are denoted by k 1 , k 3 and k 23 ; the viscous damping coefficients are c 1 , c 3 .
It is assumed that the masses of the block and the secondary oscillator are identical (m 2 = m 3 ) and relatively small (m 2 m 1 and m 3 m 1 ). Impacts of the block against the main body are characterized by the restitution coefficient κ. The case of imperfectly elastic collisions is considered, that is 0 < κ < 1. Moreover, the spring coupling between the left and right subsystem is assumed to be weak (k 23 k 1 , k 3 ). The external force is harmonic: F 1 (t) = F 10 sin(ω 1 t).
Let x 1 , x 2 and x 3 be the displacements of the components from equilibrium (−L ≤ x 2 −x 1 ≤ L). The equations of motion of the three-degree-of-freedom system between impacts (for |x 2 − x 1 | < L) take the form: Classically, it is assumed that the duration of the collisions is very short, and any external forces are negligible compared to the impact force. When Newton's restitution rule is used together with the momentum conservation principle [31], the change in momentum of the main oscillator due to an impact is given by where μ = m 1 m 2 m 1 + m 2 andẋ − i ,ẋ + i denote the velocities of the ith body just before the impact and just after the impact, respectively. Including the effect of consecutive collisions, the dynamic equations (for where δ(•) is the Dirac delta function, and t j is the time instant of the jth impact. For convenience of analytical studies, we introduce the dimensionless time and displacements: Now, Eq. (3) can be rewritten in their counterpart nondimensional form (cf. [9][10][11]): where the overdots denote differentiation with respect to τ , and the dimensionless parameters are as follows: In order to cope with the impulsive right-hand sides of the equations, we define new variables, namely the coordinates describing (approximately) the displacement of the center of mass of the vibro-impact unit, and the relative displacement of the main body: Introducing these variables into Eq. (4) leads tö Although the left-hand sides of the equations are more complicated than before, only one differential equation (corresponding to the relative displacement) includes the impact-related sum.
3 Analytical approach to the problem

Analytical approximate solution
Let the mass ratio ε play a role of the small parameter (ε 1). Some factors within the mechanical system, such as damping, elastic coupling and dynamic excitation, are assumed to be weak, and the rescaled parametersγ 1 ,γ 3 ,α 23 ,f 10 are formally introduced: Since the construction of a general solution of the nonsmooth dynamical problem is essentially impossible, certain restrictions on the system's motion are indispensable. Our considerations are focused on a specific vibro-impact regime, i.e., the periodic motion with two symmetric impacts per cycle, near 1:1 resonance. As it was indicated by experimental investigations, this particular type of behavior predominates in dynamics of such vibro-impact systems, and is the most efficient from the perspective of targeted energy transfer [16,19]. Accordingly, we assume that whereσ 1 andσ 30 are the detuning parameters.
Applying the method of multiple scales [18,24], we express the approximate solution in terms of two different time scales: where τ k = ε k τ for k = 0, 1. Substituting parameters (8) and (7) as well as expansions (9) into Eqs. (6), and then equating coefficients of order ε 0 leads to where D n k (•) = ∂ n ∂τ n k (•). Obviously, the solutions of Eqs. (10) and (12) can be written as while for Eq. (11) we postulate where Z is the sawtooth function Its role is to reflect the non-smooth nature of the solution as an effect of the cyclic symmetric impacts that occur at τ 0, j = jπ + θ where j = 0, 1, 2, . . .. The sawtooth function is illustrated in Fig. 2a for constant C and θ . Expanding the inverses of Eqs. (5) in a Taylor series with respect to ε, and keeping terms up to second order, we obtain Thus, in the near-resonant case, the block is bounced against the inner wall of the box (from W 0 = −1 to W 0 = 1 and vice versa) whose harmonic motion is slightly disturbed by the impacts. What is more, B as well as C multiplied by the ε-related factor play the role of the amplitudes. Due to the geometric constraints, B and C are mutually dependent. As discussed in Refs. [9,10], analysis of the impact conditions (W 0 = 1 at τ 0 = θ ) described by Eq. (11) together with solutions (14) and (13) produces Combining these two equations leads to the relationship between the slow-time-scale variables which describes the so called slow invariant manifold (SIM) of the problem. Moreover, the phase angles are defined by Geometrically speaking, the curve B(C) is a set of points (C, B) corresponding to potential steady-states of the dynamical system. An exemplary curve B(C) is shown in Fig. 2b. The minimal possible value of B and the related value of C are Moreover, the relation between ρ and κ is graphically presented in Fig. 2c. Parameter ρ changes from 2/π to 0 over the whole range of the restitution coefficient. For the studies of evolution of the system on the SIM, the equations corresponding to U and X 3 at the ε 1 order of approximation are necessary: Function Z involved in W 0 can be expanded into a Fourier series with respect to τ 0 . Eliminating secular terms from U 1 and X 31 , we obtain the solvability conditions which form a system of differential equations for the amplitudes and phases: For the sake of brevity, functions g B , g φ , g A3 , g φ3 are not presented in a full detailed form. System (23) can be transformed into an autonomous one by using formulae (20), and introducing the modified phases The full form of the autonomous system is The steady-state near-resonant motions correspond to the fixed points of (26). Thus, the set of algebraic equations must be solved: Eliminating ψ and ψ 3 , we obtain where Δσ =σ 30 − 2σ 1 . Now, recalling the relation B(C) given by (19), we determine two solutions for C: where The dependence of B on C follows from (19), while amplitude A 3 can be written as

Classification of periodic solutions
Numerical experiments presented in Sect. 5 are conducted by treating some model coefficients as control parameters. However, not all of the theoretically possible responses (14) are admissible. Assessment of the steady-state periodic solutions will be based on the three factors: -geometric constraints, -stability of fixed points, -stability of periodic motions.
They will be the source of certain restrictions on values of amplitude C. In further considerations, the inverse of relation (19) is used where the curve B(C) is divided into the left (L) and right (R) branches. Moreover, the following set of dimensionless parameters is treated as the basic one: Firstly, degenerate solutions should be excluded. By degenerate solutions we mean such functions W that violate the unilateral geometric constraints |W | ≤ 1, i.e., the no-penetration condition. As illustrated in Fig. 3a, negative values of C are not allowed, although the condition W (τ 0, j ) = ±1 is satisfied. On the other hand, high positive values lead to local extrema which may exceed the constraints. The maximum in the interval (θ, θ + π) can be easily determined: and, taking into account relation (32) for the right branch, the equation W max = 1 can be solved numerically for the critical value of B. The geometrically acceptable range of C is given by where the upper limit as a function of κ is shown in Fig. 3b. For example, when κ = 0.75, C G max ≈ 6.07 for which B ≈ 5.1.
In the next step, linear stability of the steady-states corresponding to fixed points P s should be verified. The autonomous system (26) in vector notation can be expressed as follows: where prime denotes the derivative with respect to τ 1 , and Eliminating C through the relations (32), the Jacobian matrix can be evaluated analytically for both the branches of SIM: For a given fixed point P s , the eigenvalues of Classically, the corresponding steadystate periodic solution is asymptotically stable if all the eigenvalues have negative real parts [24,30]. The analytical form of the Jacobian matrix has been reported in Appendix 1. Unfortunately, its complexity does not allow to present concise stability conditions in terms of the model parameters.
Perhaps the most important restrictions imposed on C can be specified by stability studies of periodic motions of the vibro-impact unit. Such analysis is not an easy task due to the non-smooth character of the solution. It may be conducted by means of the method proposed by Masri [19][20][21]. Similar approaches to this kind of systems were presented for example in [6,7,29]. The general idea is to analyze the propagation of small initial perturbations in the solution (see Fig. 4a). However, it should be emphasized that in all of the abovementioned works, the motion of the main body between impacts was described by the general solution for a periodically excited linear oscillator with viscous damping. Then, periodicity conditions were used together with Newton's restitution rule and the principle of momentum conservation. In the next subsection, we adapt this method to the form of the presented approximate analytical solution.

Stability analysis of the periodic motion
Let us rewrite the approximate periodic solutions (17) in an equivalent form: so that Φ = φ +θ is the phase shift between the motion of the box and the internal block. For the sake of brevity, the zero subscript of τ is dropped, and τ j denotes the time instant of the jth impact (τ j = τ 0, j ). The displacement and velocity of the box can be expressed in a piecewise manner for the intervals between the impacts of number j and ( j + 1) as well as between the impacts of number ( j + 1) and ( j + 2): where v 2 =Ẋ 2 denotes the steady-state absolute velocity of the block. In the case of the original (unper- Similarly to the stability analysis discussed in [19], we initially take as the state vector of the subsystem. Using formulae (40) for the first interval (with positive v 2 , see Fig. 4b), we can describe position and velocity of the main body immediately after the jth impact in both the original and disturbed case: where d 1 and v + 1 are the given displacement and velocity of the box at the impact; the tilde denotes the steadystate values perturbed by small variations at τ = τ + j : Equations (42), and all the subsequent ones describing the disturbed state, can be linearized near the original motion (the first-order approximation). Hence, a comparison of the unperturbed and perturbed behavior, i.e., Eqs. (41) with (42), provides the relationships between the deviations: For the time instant just before the ( j + 1)th impact, we have where Comparing the second equations of (45) and (46) leads to The time of the block motion between the stops is given by where l 2 = 2d 2 , From Eqs. (49) we obtain Moreover, using solutions (40) for the second interval (with negative v 2 ), we can describe the motion immediately after the ( j + 1)th impact: where the perturbed values at τ = τ + j+1 are Confronting (42) and (53) yields Now, the discontinuities in the velocities should be described. From Newton's rule and the momentum conservation it follows thaṫ In the unperturbed case, one can puṫ where v − 1 is described by (45). In the disturbed state, in turn, the velocities should be replaced with their perturbed values: where v − 1 is given by (48). Again, by comparing the original and perturbed motion, we get the new relations: As can be seen from Eqs. (44) and (55), the perturbations of the components of z * are not mutually independent. Thus, it is possible to reduce the state vector to Substituting (44) and (55) into Eq. (51) gives Again, using (44) and (55) combined with Eqs. (59) leads to one and the same relation: Finally, Eqs. (60) and (61) form the system of relationships between the perturbations occurring after the consecutive impacts: where and D j is a kind of transition matrix. Its elements are dependent on C, κ and ε in the considered problem. By analogy, this whole procedure can be applied to the second half period of the analyzed cyclic motion, which gives However, the two-impact motion is assumed to be symmetric, therefore D j+1 = D j and the final transformation has the form Stability of the periodic motion can be determined on the basis of the eigenvalues r 1 and r 2 (real or complex) of constant matrix D. The motion is asymptotically stable if |r i | < 1 for every i [19,21]. The resulting stability regions on the planes (κ, C) and (ε, C) are shown in Fig. 5. The acceptable range of C can be written as Although formally the lower boundary depends on both κ and ε, the dependence on ε is negligible; C S min turns out to be a monotonically increasing function of κ, and it practically coincides with C min (κ) given by formula (21). The upper stability boundary, in turn, is weakly dependent on both the restitution coefficient and small parameter. When κ = 0.75 and ε = 0.1, the critical values are C S min ≈ 0.99 and C S max ≈ 1.68. For comparison with literature results, the stability region has been plotted on another plane (see Fig. 6),  10 .
Such a graphical representation of stability boundaries was shown, for example, by Masri [19,21]. His results were reported for the mass ratio ε = 0.05 and damping ratio that is ζ = εγ 1 /2 = 0.01. These values have been taken into account in Fig. 6. A visual comparison of the results reveals a close similarity of the boundaries shapes and, above all, a high agreement of the characteristic ξ values at κ = 0 and κ = 1. The analytical expression for matrix D j as well as an approximation of C S min and C S max have been presented in Appendix 1. Unfortunately, the resulting complicated form of the eigenvalues r 1 and r 2 does not allow to formulate concise stability conditions. Figure 7 presents different parts of the SIM related to the three considered criteria. The geometrically allowable solutions can be found along the branch of stable fixed points (SFP) within the range of stable periodic motions (SPM, see Fig. 7a). The final stability region (S) is illustrated in Fig. 7b. For the basic data set (33), only one fixed point belongs to the stable branch. It is worth noting that the location and length of the stable segment of the SIM are comparable with the results reported in Refs. [13][14][15][16].

Analytical and numerical results
Let us start with the results of an analytical nature. The time histories of U , W and X 1 , X 2 , X 3 corresponding to the stable fixed point P 2 (Fig. 7b) are plotted in Fig. 8a. As can be seen, the effect of collision on the box motion is hardly visible. Due to the very weak coupling between the block and the secondary oscillator, the latter one moves harmonically. The smooth and non-smooth types of the components behavior can be better seen in the phase diagrams (Fig. 8b).
The size of the stability zone and the location of fixed points on the SIM are strongly affected by the model parameters. For instance, when treating the excitation amplitudef 10 as a control parameter, we can observe various changes in the existence and stability of the periodic solutions (see Fig. 9). Asf 10 is decreased to 1.2, points P 1 and P 2 get closer to each other, and P 2 still remains stable. Forf 10 = 0.6, P 2 is shifted to the left branch and becomes unstable. However, increasing the forcing amplitude fromf 10 = 1.5 tof 10 = 2.5 also lead to instability of P 2 . In the case of two unstable fixed points, the assumed periodic solution with two symmetric impacts per cycle is not valid. Therefore, the actual motion of the mechanical system should be determined via numerical simulations. Figure 10 contains the plots of the frequencyresponse curves for various values off 10 andγ 1 . The result corresponds only to the right branch of the SIM (C > C S min ). The characteristics resemble the case of linear system, but they consist of stable and unstable segments. Functions C(σ 1 ) and A 3 (σ 1 ) are mostly single-valued, since there is no spring-related nonlinearity bending the response curves. Only for low forcing amplitudes and small damping, one stable and one unstable branch occur simultaneously forσ 1 < 0. Although the two subsystems are weakly coupled, amplitude A 3 noticeably increases withf 10 . However, it should be emphasized that A 3 is a function of the parameters that are crucial for stability of the system as well as the other parameters (α 23 ,σ 30 ,γ 3 ).
The variation of amplitude C with the excitation amplitude for several values of the detuning parameter is presented in Fig. 11. Naturally, the curves can also be divided into stable and unstable segments. Functions C(f 10 ) seem to be linear forσ 1 > 0, and they take more complicated shapes forσ 1 < 0. In the two cases, the length of the stable branches increases with increasing |σ 1 |.
In order to get a deeper insight into the effect of various parameters on the occurrence of the steady-state periodic response of the system, more systematic calculations have been conducted. Figure 12 shows several parameter planes with the regions of existence and stability of the solutions. The variation ranges and steps of the parameters are given in Table 1. The domains of one stable solution are filled with blue. In the gray-shaded areas the two existing solutions are unstable. The white areas, in turn, relates to the lack of periodic solutions. As it could be expected, low amplitudes of the excitation are insufficient to ensure a steady state in the form of two-impact symmetric motion. Such regions and the stability zones are usually separated by relatively narrow instability strips. Moreover, unstable regions arise near the axisσ 1 = 0. Apart from amplitudef 10 , the restitution coefficient is crucial for the existence and stability of the solutions. On the plane (κ,f 10 ), for example, the stable region is bounded from above and below. Similar effect can be observed for the pair (γ 1 ,f 10 ), but the upper critical force level increases withγ 1 . When it comes to the interplay between κ and γ 1 , high damping values together with low restitution coefficient may exclude the occurrence of the periodic response. Now let us turn to purely numerical results. It should be noted that simulation studies of vibro-impact dynamics can be a challenging task from the computational viewpoint. Firstly, long-term motion of the system should be analyzed to omit the transient response, which becomes cumbersome in parametric investigations. Secondly, a special attention must be paid to collision detection, i.e., a very accurate determination of the impact time instants. In what follows, Eq. (1) in a nondimensional form are integrated numerically between impacts by using the Runge-Kutta-Fehlberg method (RKF45). As soon as the unilateral constraints are violated, the algorithm switches to the classical fourthorder Runge-Kutta method (RK4) modified with the concept of time-step bisection.
Numerical solutions for the basic data set as well as forf 10 = 0.6 are presented in Fig. 13. In the first case, function W has a similar shape as the sinusoidalsawtooth response given in Fig. 8. The numerically obtained amplitudes of X 1 and X 3 are very close to the analytical ones. Time histories for the unstable case are shown for a longer time interval. As can be seen, time history of the relative displacement W is composed of resonant subintervals (two-impact motion) and shorter episodes of irregular motion. The main oscillator exhibits the so-called strongly modulated response (SMR). Such a behavior was reported in many papers devoted to VI NES of this kind (e.g., see [9,11,15,16]). What is more, the secondary oscillator coupled with the colliding block undergoes a similar beating-like motion, but with a slighter amplitude modulation. Completely different complexity of the system behavior in the two cases is also reflected in the phase portraits.
As a result of parametric studies, the bifurcation diagrams of the displacement X 2 are shown in Fig. 14a. The excitation amplitude is the control parameter. The first diagram was calculated for 0.5 ≤f 10 ≤ 5 with a step Δf 10 = 0.01. The initial n t = 400 excitation periods were discarded, and then motion of the system was analyzed for n s = 800 periods. At the beginning of a new step, the last solution (for the previous value of f 10 ) was used as the initial conditions. Instead of using the Poincare map, local maxima of X 2 were found in each period [26].
The system dynamics turns out to be regular in a large part of the analyzed range. The motion is mostly periodic in the subintervals 0.9 ≤f 10 ≤ 2.8, 2.9 ≤ f 10 ≤ 3.4 and 3.8 ≤f 10 ≤ 4. The periodic windows are surrounded and separated by regions of chaotic behavior. The relatively narrow smears of points at (a) (b) Fig. 10 Frequency-response curves of the system: a effect off 10 , b effect ofγ 1 ; stable (bold solid) and unstable (dashed) branches the transition between these two regimes may indicate quasi-periodic motion. More details of the qualitative changes in the system dynamics can be seen in the second diagram computed with the same technique for 1.9 ≤f 10 ≤ 4 with higher resolution: Δf 10 = 0.002.
In order to more precisely specify the character of vibro-impact motion, spectra of the Lyapunov exponents (LEs) {λ 1 , λ 2 , . . . , λ 7 } of the system were determined within the same ranges off 10 . In the procedure of numerical integration of the variational equations along with the equations of motion, we have applied the Müller algorithm for non-smooth dynamical systems [22] and the Gram-Schmidt reorthonormalization [25,26]. The calculations were performed over n λ = 1500 excitation periods, and were terminated if the maximal value of standard deviation of the set of LEs within the last m λ steps decreased below a given threshold.  The diagrams of the three largest Lyapunov exponents are presented in Fig. 14b. As can be seen, the ranges of the control parameter where the maximal exponent λ 1 = 0 correspond to the periodic windows. In the regions of chaotic behavior, in turn, λ 1 becomes positive and λ 2 = 0. Moreover, it seems that λ 1 = λ 2 = 0 at some single points. The third LE is always negative. Table 2 contains the LE spectra in a numeric form for selected values of the excitation amplitude. Now it is clear that whenf 10 increases and the periodic motion with two symmetric impacts loses its stability, the system exhibits a sequence of period-doubling bifurcations. Next, a relatively narrow range of chaotic motion arises. Within the next two periodic windows, this classical route to chaos is not so evident, but it can be observed that 2T -and 4T -periodic motions predominate. Phase portraits of the periodic and chaotic solutions forf 10 > 2 are shown in Fig. 15. It should be noticed that the bifurcation diagram based on local maxima can be confronted directly with the amplitude of the periodic steady-state motion investigated analytically. According to Eq. (17) the amplitude of X 2 is given by A 2 = (1 − ε + ε 2 )C. By analogy to Fig. 11, we can plot A 2 as a function off 10 when the other parameters are kept constant. The resulting force-response curve overlaid on the bifurcation diagram of X 2 is depicted in Fig. 16. As can be seen, there is a quite good agreement between the results. The amplitude based on the approximate analytical solution is slightly lower than the numerical values. Moreover, although the length of the stable periodic range may seem to be overestimated, it can be observed that the numerical solutions atf 10 ≈ 1 andf 10 ≈ 2.1 are mostly 1T -periodic.

Limitations and extensibility of the approach
The approximate analytical approach can be regarded as a specific perturbation technique, where the small parameter, equivalent to the mass ratio, plays the central role. Usually, in the investigations of various VI NESs the parameter value is assumed to be very small, i.e., ε ≤ 0.02. On the contrary, in our studies, we have taken ε = 0.1 as the limiting value: the masses m 2 and m 3 are at least one order of magnitude smaller than m 1 . So, a question about sufficiency of this assumption and accuracy of the solutions arises.
A graphical comparison between the analytical and numerical time histories of U and W for three values of ε is shown in Fig. 17a. For the sake of clarity, the plots are focused on one period of the system motion. Although the vibro-impact behavior reflected by W (always limited by the values A W = 1 and −1) is constantly well approximated, the differences for U increase with ε considerably. Obviously, it has its consequences for solutions X 1 , X 2 and X 3 (it is visible, for example, in Fig. 16). The broader results included in Fig. 17b indicate that the discrepancies in the amplitudes B, A 1 , A 2 and A 3 grow linearly with the small parameter. When ε = 0.1, the amplitude differences for B and A 1 are about 20% (with reference to the numerical values), while the differences in the other amplitudes do not exceed 8%.  It should be noticed that the small parameter is involved in many stages of the procedure, which leads to various simplifications. For instance, the assumption m 2 /m 1 = m 3 /m 1 = ε (i.e., m 2 = m 3 ) together with α 23 = ε 2α 23 results in a very specific idealization of the system. Indeed, the coupling between the colliding block and the secondary oscillator is weak (of order ε 1 ), and their interaction is neglected at the first approximation. Consequently, the stability analysis of the periodic motion can be restricted to the vibro-impact unit, and dynamics of the body of mass m 3 is considered only in the stability analysis of the steady states.
Moreover, the analytical approach has been applied just to the simplest type of the system behavior: the periodic motion with two symmetric impacts per cycle. As it can be concluded from numerous simulations conducted for this and similar systems, the 2T -, 4T -or 8T -periodic motions do not include simple impacts patterns (e.g., equal time intervals between impacts, only full-span travels between the stops), which practically Fig. 16 Comparison of the analytical and numerical results: the force-amplitude curve A 2 (f 10 ) overlaid on the bifurcation diagram makes it impossible to construct analytical solutions to such vibro-impact problems.
As a potential extension of our considerations, we will outline the analysis suited to the case when the assumption m 3 /m 1 = ε is eliminated. So, let m 2 /m 1 = ε while m 3 /m 1 ≈ 1. Now, the equations of motion are given by [cf. Eq. (4)] where .
Let Ω 2 30 = 1 + εσ 30 . If α 23 = εβ 2 and α 32 = εμ 13β 2 where μ 13 = m 1 /m 3 , at the level of the first approximation the vibro-impact motion of the block will be affected by vibrations of the secondary oscillator (but not vice versa). After all, the resulting model will be much more realistic. For the case of weak damping, Eqs. (10) and (12) The solution of this equation can be written as where Z * is the modified sawtooth function [5,8]: For computational purposes, Z * can be expressed as a continuous function, by analogy to formula (15). It should be noticed that Z * reduces to Z in the limiting case forβ = 0 (see Fig. 18). Due to the extended form of W 0 , the relationship between the slow-time-scale variables is much more complicated than before, and involves not only B and C, but also A 3 , φ and φ 3 [cf. Eq. (19)]: with Therefore, searching for the steady states becomes hindered. Moreover, the stability analysis of the periodic motions is more difficult, since the spring described byβ participates in the travel of the block between the stops and this motion is not uniform.
As can be seen, despite the high value of the small parameter, the discrepancies in the amplitudes obtained analytically and numerically are relatively small.

Conclusions
The three-degree-of-freedom system under a harmonic force has been considered. The analytical studies have been restricted to a periodic steady-state motion with two symmetric impacts per cycle near 1:1 resonance. The multiple scales method combined with the sawtooth-function-based modelling of the non-smooth dynamics has been successfully applied to the problem.
Taking into account geometric constraints and stability criteria, certain restrictions on the amplitudes have been specified. The central role belongs to the stability studies of periodic motions. The presented method is something between the approach involving the general solution for a periodically excited, viscously damped piecewise linear oscillator [19,29] and the studies of the Poincaré sections of periodic solutions for discrete dynamical systems (e.g., see [17]).
The approximate analytical approach has allowed us to determine both the frequency-response and forceresponse curves having stable and unstable branches. Moreover, the interplay between the model parameters and their effect on the existence and stability of the periodic solutions have been determined on various parameter planes. When treating the excitation amplitude as the control parameter, one can make the system leave the stability zone. Consequently, the response of the primary oscillator is strongly modulated. This kind of motion is largely followed by the secondary oscillator, in spite of a weak coupling between the subsystems.
The bifurcation diagrams and the Lyapunov exponents spectra have shown that the window of stable 1T -periodic solution is surrounded by the regions of chaotic motion and period-doubling route to chaos. The To sum up, in the area of near-resonant dynamics, the analytical approach is a useful tool in stability studies, estimation of the motion amplitudes, and systematic parametric studies. However, an extension of the applicability of the method to more complex problems with weaker assumptions (e.g., strong and/or nonlinear coupling, mass ratio of order 1) is highly desirable.

Data Availability Statement
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

A Appendix: details of stability analysis
In this appendix, the analytical expressions for the Jacobian matrix J and the transition matrix D, as well as the approximation of the stability boundaries are reported. Let