A Solution Procedure Combining Analytical and Numerical Approaches to Investigate a Two-Degree-of-Freedom Vibro-Impact Oscillator

In this paper, a new approach is proposed to analyze the behavior of a nonlinear two-degree-of-freedom vibro-impact oscillator subject to a harmonic perturbing force, based on a combination of analytical and numerical approaches. The nonlinear governing equations are analytically solved by means of a new analytical technique, namely the Optimal Auxiliary Functions Method (OAFM), which provided highly accurate explicit analytical solutions. Benefiting from these results, the application of Schur principle made it possible to analyze the stability conditions for the considered system. Various types of possible motions were emphasized, taking into account possible initial conditions and different parameters, and the explicit analytical solutions were found to be very useful to analyze the kinetic energy loss, the contact force, and the stability of periodic motions.


Introduction
Vibro-impact processes are widely used in mechanical-engineering applications and devices such as hammer-like devices, rotor-casing dynamical systems, wheel-rail interaction of high-speed railway couches, heat exchangers, fuel elements of nuclear reactors, gears, piping systems [1], stiction and electric short circuits of MEMS [2], noise and wearproducing processes, grinding, foraging, drilling, punching, tamping, pile cutting a variety of objects, and vibrating machinery or structures with stops or clearance [3].
Dynamics of vibro-impact systems has received great attention in the last decades. Some numerical or experimental results were obtained, but adequate analytical predictions need to be further studied. For the first time, Masri [4] analyzed the stability of the impact damper by using the perturbation method. The finite-difference method was used by Kobrinski [5], Babitsky [6], Brîndeu [7], and Okolewska et al. [8]. The nonlinear equation of a single-degree-of-freedom oscillator with an impact damper was analyzed approximately using the method of Kryloff and Bogoliuboff by Cronin and Van [9]. Shaw [10] studied the symmetric double-impact motion, both harmonic and subharmonic, by means of a mapping that related conditions at subsequent impacts. A mechanical system with one degree of freedom was investigated by Ivanov [11] using standard Poincare-Bendixon theory, Lyapunov's second method, and Zhuravlev transformations. The exact solutions of the vibro-impact oscillator with two degrees of freedom, including localized states and their bifurcation structure, were obtained by Aziz et al. [12]. Stability and bifurcation for the unsymmetrical, periodic motion of a horizontal impact oscillator under a periodic excitation were employed through four mappings based on two switch planes by Luo [13]. De Souza and Caldas [14] implemented the Ott-Grebogi-Yorke (OGY) method to stabilize desired unstable periodic orbits by applying a small perturbation on an available control parameter and introducing an impact map for the dynamical variables. Chaos and periodic motion of a cantilever beam system with impacts was found readily by Emans et al. [15]. In addition, nonlinear behavior such as coexistence of attractors was found even at modest oscillation levels during investigations with realistic parameters.
The dynamics of an impact oscillator with viscoelastic and Hertzian contact was examined by Mann et al. [16] using the presence of multiple periodic attractors, subharmonics, quasi-periodic motions, and chaotic oscillations. Avramov [17] considered the forced vibrations of an impact Duffing oscillator using Zhuravlev's transformation and the multiple-scale method. The stability and bifurcations of periodic vibrations were explored. The control of vibro-impact dynamics of a single-sided Hertzian contact forced oscillator was considered by Bichri et al. [18]. The multiple-scale technique was applied to study the slow dynamic and the frequency-response curves near the primary resonance. Based on Zhuravlev and Ivanov transformations, Grace et al. [19] developed an analytical model of a ship's roll motion interacting with ice. Extensive numerical simulations were carried out for all initial conditions covered by the ship's grazing orbit for different values of excitation amplitude and frequency of the external wave roll moment.
Jovic et al. [20] considered the motion trajectory of a vibro-impact system based on the oscillator moving along a rough parabolic line in the vertical plane, under the action of extreme force with nonlinearity of the bond originates of sliding Coulomb's type friction force. The study of Askari and Tahani [2] was focused on the effect of mechanical shock on dynamic pull-in instability of electrically actuated microbeams through an alternative reduced-order model and using the fourth-order Runge-Kutta method. By setting up a Poincaré map and using a shooting method, Lin et al. [21] obtained a fixed point of periodic motion in the system, period doubling bifurcation, and Hopf bifurcation. The existing and stability conditions of period-1 motion in a single-degree-of-freedom oscillator with double-side constraints were studied by Wang et al. [22]. A smaller shock gap than impact gap could make the periodic motion more stable.
Reboucas et al. [23] examined limitations of the coefficient of restitution model using experimental observations from a simple vibro-impact setup, and the effect of the magnitude response of gap deviation during an experimental sweep. Numerical simulation using Zhuravlev and Ivanov transformations was obtained. A triboelectric energy harvester on a three-degree-of-freedom vibro-impact oscillator was investigated by Fu et al. [24]. The symmetric mass configurations of the oscillator were more beneficial to energy harvesting than the asymmetric cases. Zukovic et al. [25] explored a vibro-impact system consisting of a crank-slider mechanism with one attached oscillator. The cases with an ideal and a non-ideal excitation were analyzed, and analytical and numerical solutions were obtained for the vibrating system with the ideal excitation. For the system with non-ideal excitation, the average value of the frequency was obtained.
The present paper analyzes the motion of a horizontal vibro-impact system with two degrees of freedom in the case when the external coercive force and viscous damping force are known and a periodic vibro-impact model is realized in the system. A variety of possible cases of motion of the vibro-impact system is presented. This process is very complex, and a nonlinear differential equation for one of the masses and a linear differential equation for the other mass is needed in the analysis. It is very difficult to obtain an exact solution for these types of equations, but an analytical approximate solution with a high accuracy will be given by means of OAFM [26][27][28][29][30][31]. The condition for periodic motion is obtained and the stability of periodic motion is studied.

Dynamical Model of the Vibro-Impact System
The considered dynamical model of the vibro-impact system under investigation is depicted in Figure 1, which details the two-degree-of-freedom system with gaps. The system under investigation is composed of a horizontal direction periodical force F = Acosωt and two oscillators of mass M1 and M2, which also move in the horizontal direction. These two masses are initially separated by three gaps: δ1 between the wall W1 and the mass M1, δ2 between the masses M1 and M2, and δ3 between the mass M2 and the wall W2. The mass M1 is connected to a linear spring K1, nonlinear spring K*, and the damper c1, and the mass M2 is connected to linear spring K2 and damper c2. The oscillator M2 collides with the rigid wall W2 and with the mass M1. The coefficient of restitution between M1 and M2 is R, and R2 is the coefficient of restitution between M1 and W1 and M2 and W2. The displacement between the mass M1 and vertical wall W1 is x1, and the displacement between the mass M2 and the wall W1 is x2.
The equations of motion between any two successive collisions are given as: where the prime denotes the derivative with respect to time t. By means of the following transformations: Equations (1) and (2) can be rewritten in the dimensionless variables, after some simple manipulations: where the dot denotes the derivative with respect to variable τ. Between two successive collisions, the initial conditions of Equations (4) and (5) are: where V1 and V2 are known parameters. Equations (4)-(6) are second-order nonlinear differential equations, and therefore are very difficult to be analytically solved. In the following, we present a relatively new approach, namely the Optimal Auxiliary Functions Method, to obtain an approximate analytical solution in an efficient manner.

Basics of the Optimal Auxiliary Functions Method (OAFM)
Following basic concepts of OAFM presented in [26][27][28][29][30][31], in this work, we consider a general nonlinear equation in the form: The system under investigation is composed of a horizontal direction periodical force F = Acosωt and two oscillators of mass M 1 and M 2 , which also move in the horizontal direction. These two masses are initially separated by three gaps: δ 1 between the wall W 1 and the mass M 1 , δ 2 between the masses M 1 and M 2 , and δ 3 between the mass M 2 and the wall W 2 . The mass M 1 is connected to a linear spring K 1 , nonlinear spring K*, and the damper c 1 , and the mass M 2 is connected to linear spring K 2 and damper c 2 . The oscillator M 2 collides with the rigid wall W 2 and with the mass M 1 . The coefficient of restitution between M 1 and M 2 is R, and R 2 is the coefficient of restitution between M 1 and W 1 and M 2 and W 2 . The displacement between the mass M 1 and vertical wall W 1 is x 1 , and the displacement between the mass M 2 and the wall W 1 is x 2 .
where the dot denotes the derivative with respect to variable τ. Between two successive collisions, the initial conditions of Equations (4) and (5) are: where V 1 and V 2 are known parameters. Equations (4)-(6) are second-order nonlinear differential equations, and therefore are very difficult to be analytically solved. In the following, we present a relatively new approach, namely the Optimal Auxiliary Functions Method, to obtain an approximate analytical solution in an efficient manner.

Basics of the Optimal Auxiliary Functions Method (OAFM)
Following basic concepts of OAFM presented in [26][27][28][29][30][31], in this work, we consider a general nonlinear equation in the form: where L and N are the linear differential and nonlinear differential operators, τ is the independent variable, and U(τ) is an unknown function. The initial or boundary conditions are: In order to obtain an approximate solution U(τ), we consider that this approximate solution can be expressed in the form: in which U 0 (initial approximation) and U 1 (the first approximation) will be determined as described below, and with C 1 , C 2 , . . . ,C n being the unknown parameters whose values will be optimally determined. Substituting Equation (9) into Equation (7), one obtains: The initial approximation U 0 (τ) can be obtained from the linear equation: and therefore, the first approximation is obtained from Equations (10) and (11): The nonlinear term from Equation (12) can be expanded in the form: To accelerate the convergence of the first approximation U 1 (τ,C i ) and therefore of the approximate solution as well, and also to avoid the difficulties that can appear in solving the nonlinear differential Equation (12), we propose another expression, such that Equation (12) can be written as follows: where C i are arbitrary unknown parameters and f i (τ) are auxiliary functions depending on the initial approximation U 0 (τ), on the functions which appear in N[U(τ)] or N[U 0 (τ)] or are combinations of such expressions. These auxiliary functions f i (τ) are very important and are not unique. It should be emphasized that we have much freedom to choose such auxiliary functions, because it is known that the nonlinear differential equation does not have a unique solution. The parameters C i , which appear into Equations (9) and (14), can be optimally identified via rigorous methods such as the Galerkin method, the least square method, the Ritz method, the collocation method, or by minimizing the square residual error. In this last case, if: where U(τ) is given by Equations (9), (11), and (14), then The values of the parameters C i are obtained from the following system: It is clear that in this way, the approximate solution U(τ) is well determined after the identification of the optimal values of the initially unknown convergence-control parameters C i .

Application of the Optimal Auxiliary Functions Method to Equations (4)-(6)
The approximate solutions and exact solutions of Equations (4)-(6) are of the following forms: For the nonlinear differential Equation (4) and for the linear Equation (5), the linear operators are, respectively: The initial approximations X 0 (τ) and the exact solution Y 0 (τ) are obtained from the linear equations: .. ..
and therefore: The nonlinear operator for Equation (4) is obtained from Equations (7) and (20): and the corresponding nonlinear operator for Equation (5) is obtained from Equations (7) and (29) as: Taking into account the initial conditions (6) for the first variable X(τ) and the initial conditions for (22) for the initial approximation X 0 (τ), one can write the initial conditions for the first approximation X 1 (τ): Mathematics 2021, 9, 1374 6 of 17 The initial conditions for Y(τ) are: From the Equations (26) and (28) and from the Equations (27) and (29), we can choose the first approximation in the forms: The approximate solution of Equations (4) and (6) and exact solution of Equations (5) and (6) are respectively: The parameters Q 1 and Q 2 can be determined by the Galerkin method, Q 3 , Q 4 can be determined by identification of coefficients, and C 3 , C 4 can be determined from Equation (29): The approximate solution of Equations (4)-(6) and exact solution are obtained using Equations (32)-(34): Now, to prove the accuracy and the effectiveness of our method, we consider a set of numerical values for the physical parameters involved in the governing equations; i.e., M 1 = 0.072, M 2 = 0.086, K 1 = 110, K 2 = 141, a = 3. Figure 2 presents a comparison between the approximate solution (35) and corresponding numerical integration results obtained by means of a fourth-order Runge-Kutta method.
This comparison emphasizes the high accuracy of the proposed analytical solution obtained through OAFM. Figure 3 presents the corresponding phase portraits. Now, to prove the accuracy and the effectiveness of our method, we consider a set of numerical values for the physical parameters involved in the governing equations; i.e., M1 = 0.072, M2 = 0.086, K1 = 110, K2 = 141, a = 3. Figure 2 presents a comparison between the approximate solution (35) and corresponding numerical integration results obtained by means of a fourth-order Runge-Kutta method.
(a) (b) This comparison emphasizes the high accuracy of the proposed analytical solution obtained through OAFM. Figure 3 presents the corresponding phase portraits.

Analysis of Non-Periodic Motion
In this section, we present some possible situations, depending on the motions of the masses M1 and M2 of the dynamical system considered in Section 2. Depending on the perturbing force F and initial conditions, we present the following possible cases of nonperiodic motions.
, then the oscillator M1 stops at the wall W1. The velocity of the mass M1 after the impact will be numerical values for the physical parameters involved in the governing equations; i.e., M1 = 0.072, M2 = 0.086, K1 = 110, K2 = 141, a = 3. Figure 2 presents a comparison between the approximate solution (35) and corresponding numerical integration results obtained by means of a fourth-order Runge-Kutta method. This comparison emphasizes the high accuracy of the proposed analytical solution obtained through OAFM. Figure 3 presents the corresponding phase portraits.

Analysis of Non-Periodic Motion
In this section, we present some possible situations, depending on the motions of the masses M1 and M2 of the dynamical system considered in Section 2. Depending on the perturbing force F and initial conditions, we present the following possible cases of nonperiodic motions.
, then the oscillator M1 stops at the wall W1. The velocity of the mass M1 after the impact will be

Analysis of Non-Periodic Motion
In this section, we present some possible situations, depending on the motions of the masses M 1 and M 2 of the dynamical system considered in Section 2. Depending on the perturbing force F and initial conditions, we present the following possible cases of non-periodic motions.
Case 5.1. The oscillator of mass M 1 is moving from the static equilibrium position to the position of the impact with the wall W 1 . If the moment of impact M 1 with W 1is τ 1 , then this case takes place if X(τ) ≤ 1, X(τ 1 ) = 0 and . X(τ 1 ) ≤ 0. In particular, if . X(τ 1 ) = 0, then the oscillator M 1 stops at the wall W 1 . The velocity of the mass M 1 after the impact will be .
, where the symbols "+" and "−" indicate the moments after and before the impact, respectively, and R 2 is the coefficient of restitution. For the oscillator of mass M 2 , the following situations are possible (types of movement): (a) The oscillator M 2 is moving between the oscillator of mass M 1 and the wall W 2 without impact and without sticking motion. The condition is equivalent to the existence of the inequalities X(τ) < Y(τ) < 1 + δ 2 /δ 1 + δ 3 /δ 1 and τ < τ 1. It follows that there exists a moment of time τ * 1 when M 2 stops, . Y(τ * 1 ) = 0. Y(τ 2 ). After this collision, the velocities can be written as: (c) The oscillator of mass M 1 is sticking to M 2 (sticking motion). The condition is written as X(τ 2 ) = Y(τ 2 ) and . X(τ 2 ) = .
Y(τ 2 ). The contact force corresponding to masses M 1 and M 2 is obtained from Equations (4) and (5): Case 5.2. The oscillator of mass M 1 is moving without impact from the position of static equilibrium between the wall W 1 and the oscillator of mass M 2 . The oscillator of mass M 1 stops at the moment τ = τ 1 before reaching the wall W 1 , then stops at the moment τ = τ 1 , before reaching the oscillator of mass M 2 , so that . X(τ 1 ) = . X(τ 1 ) = 0. For the oscillator of mass M 2 , one can identify two situations: (a) The oscillator of mass M 2 is moving without impact between the oscillator of mass M 1 and the wall W 2 , so that one can write X(τ) < Y(τ) < 1 + δ 2 /δ 1 + δ 3 /δ 1 . It results that there exist two moments: τ 2 , when the oscillator of mass M 2 stops before the oscillator M 2 ; and τ 2 , when the oscillator M 2 stops before the wall W 2 , so that . Y(τ 2 ) = .
(c) The oscillator of mass M 2 is moving from the position of static equilibrium toward the wall W 2 without impact with the oscillator of mass M 1 or with the wall W 2 . This subcase takes place if X(τ) ≥ 1, Y(τ) ≥ 1 + δ 2 /δ 1 , and X(τ) < Y(τ). It is clear that Y(τ 7 ) = 0, which means the oscillator of mass M 2 stops. (d) The oscillator of mass M 2 is moving from the position of static equilibrium toward the wall W 2 , but with impact between the two oscillators. The impact takes place before the oscillator of mass M 2 to reach at the wall W 2 and without sticking. The conditions will be X(τ) ≥ 1, Y(τ) ≥ 1 + δ 2 /δ 1 , and there exists a moment τ 7 when X( . If the impact of the oscillators takes place with sticking, then the condition is similar to (38), F 12 (τ) > 0, τ > τ 7 . (e) The oscillator of mass M 2 is moving toward the wall W 2 , and the impact with W 2 takes place at the moment τ 8 , but there is no impact between the two oscillators. Therefore, the following conditions should be satisfied: .
After the impact, the velocity of the oscillator of mass M 2 is .
. It should be stated that after the two oscillators undergo one of the movements studied before, the movement of each oscillator should be studied by taking into account the conditions mentioned for the corresponding case. For example, after Case 5.1a takes place, therefore in the condition . X(τ 1 ) < 0 (the oscillator of mass M 1 does not stop), Y(τ * 1 ) = 0, then the approximate solutions (35) and (36) will be written as: where τ * 1 is obtained from the equation is given by Equation (36); and Q 3 and Q 4 are given by Equation (34). Figure 4 shows the contact force of (38). where * 1 τ is obtained from the equation is given by Equation (36); and Q3 and Q4 are given by Equation (34). Figure 4 shows the contact force of (38). It can be seen that the contact force has two maximum points at τ = 0.8 and τ = 2.4, and is null for four points: τ = 0.4, τ = 1.3, τ = 1.9, and τ = 2.9.
Another specific element for the impact phenomenon is the kinetic energy loss. At the moment of impact, the kinetic energy is transformed into work, caloric energy, etc. The kinetic energy loss is given by the difference between the kinetic energy at the beginning and at the end of impact: It can be seen that the contact force has two maximum points at τ = 0.8 and τ = 2.4, and is null for four points: τ = 0.4, τ = 1.3, τ = 1.9, and τ = 2.9.
Another specific element for the impact phenomenon is the kinetic energy loss. At the moment of impact, the kinetic energy is transformed into work, caloric energy, etc. The kinetic energy loss is given by the difference between the kinetic energy at the beginning and at the end of impact: In Case 5.1b, the kinetic energy loss becomes: (44) Figure 5 presents the kinetic energy loss for Case 5.1b. It can be seen that the contact force has two maximum points at τ = 0.8 and τ = 2.4, and is null for four points: τ = 0.4, τ = 1.3, τ = 1.9, and τ = 2.9.
Another specific element for the impact phenomenon is the kinetic energy loss. At the moment of impact, the kinetic energy is transformed into work, caloric energy, etc. The kinetic energy loss is given by the difference between the kinetic energy at the beginning and at the end of impact: In Case 5.1b, the kinetic energy loss becomes: (44) Figure 5 presents the kinetic energy loss for Case 5.1b.

Analysis of Periodic Motion
The solutions (35) and (36) are written for the initial conditions in (6), in the absence of impact. After some time, the free vibration described by (35) and (36) disappear, so that the movement becomes periodic with the same period of the perturbing force F, which is T = 2π/Ω.
The problem of the stability of periodic motion is of interest only in the case of the impact between the two oscillators, since the impact with the walls W 1 and W 2 is assumed to be elastic. In the general case, instead of the conditions in (6), we consider that the approximate solution X(τ) and the exact solution Y(τ) are: where Q 3 and Q 4 are given by Equation (34), and C 1 -C 4 and φ are determined from initial conditions. Solutions (45) and (46) are stable if a negligible change of a small perturbation takes place in the initial conditions and delay at any time interval. If the change is significant, then the solutions become unstable.
According to the Schur criterion, the conditions in (69) are satisfied if the following inequalities take place: In order to study the stability conditions in (7), the parameters C 1 -C 4 and the delay φ should be determined. In this respect, the periodicity conditions will be used: where A and B were defined in Equation (65). From Equations (45), (46), and (71), one obtains the system of equations: . .
A graphical analysis of the stability conditions is presented in Figures 6-9, taking into account the influence of different parameters; the yellow color denotes a zone of stability.

Conclusions
In this paper, we presented a study on the dynamical model of a two-degree-of-freedom vibro-impact oscillator with damping, subject to a perturbing force. Within this system, each of the masses may produce a combination of continuous movements or vibroimpact movements that are quite complicated. The system under study was composed of a linear oscillator and nonlinear one whose movements were determined by means of an

Conclusions
In this paper, we presented a study on the dynamical model of a two-degree-of-freedom vibro-impact oscillator with damping, subject to a perturbing force. Within this system, each of the masses may produce a combination of continuous movements or vibroimpact movements that are quite complicated. The system under study was composed of a linear oscillator and nonlinear one whose movements were determined by means of an

Conclusions
In this paper, we presented a study on the dynamical model of a two-degree-of-freedom vibro-impact oscillator with damping, subject to a perturbing force. Within this system, each of the masses may produce a combination of continuous movements or vibroimpact movements that are quite complicated. The system under study was composed of a linear oscillator and nonlinear one whose movements were determined by means of an

Conclusions
In this paper, we presented a study on the dynamical model of a two-degree-offreedom vibro-impact oscillator with damping, subject to a perturbing force. Within this system, each of the masses may produce a combination of continuous movements or vibroimpact movements that are quite complicated. The system under study was composed of a linear oscillator and nonlinear one whose movements were determined by means of an approximate analytical solution, very close to the numerical solution. This approximate analytical solution was optimally determined using a new method, namely the Optimal Auxiliary Functions Method (OAFM), which was very efficient in practice. This solution was further used in analyzing the kinetic energy loss, the contact force, and the stability of periodic motions. Several cases were emphasized for possible movements of the system depending on the initial conditions and movement parameters, considering various impact possibilities between the two oscillators or between the oscillators and the walls.
In the case of periodic movements, it was considered that the period of the system was proportional with the period of the perturbing force. The boundary conditions of the system between two periodic subsequent arbitrary impacts were considered. After some analytical and numerical computations and by using the Schur principle, we established the stability conditions of the system. For different values of the movement parameters, a graphical analysis of stability conditions was presented. Different expressions for the velocities of the two oscillators were considered, establishing recurrence relations for periodic solutions.
Limitations of this study include the fact that it did not cover bifurcation-and chaosrelated problems, but these aspects will be the subject of future developments. In addition, more research will be necessary to refine and further elaborate the proposed approach to make it applicable to three-degree-of-freedom vibro-impact systems, which would provide a more complete understanding of more complex systems. Moreover, as an alternative to the Schur criterion, use of the Zhukovsky-type stability [33] could be convenient, and the Sommerfeld effect [34,35] could be investigated by means of OAFM.