Periodic solutions of n-dofs autonomous vibro-impact oscillators with one lasting contact phase

In the continuum structural mechanics framework, a unilateral contact condition between two flexible bodies does not generate impulsive contact forces. However, finite-dimensional systems, derived from a finite element semi-discretization in space for instance, and undergoing a unilateral contact condition, require an additional impact law: Unilateral contact occurrences then become impacts of zero duration unless (i) the impact law is purely inelastic, or (ii) the pre-impact velocity is zero. This contribution explores autonomous periodic solutions with one contact phase per period and zero pre-impact velocity [case (ii)], for any n-dof mechanical systems involving linear free-flight dynamics together with a linear unilateral contact constraint. A recent work has shown that such solutions seem to be limits of periodic trajectories with kimpacts per period as k increases. Minimal analytic equations governing the existence of such solutions are proposed, and it is proven that, generically, they occur only for discrete values of the period. It is also shown that the graphs of such periodic solutions have two axes of symmetry in time. Results are illustrated on a spring–mass system and on a 4-dof two-dimensional system made of 1D finite elements. Animations of SPPs with up to 30 dofs are provided.


Introduction
When a continuous medium collides with a rigid obstacle, (i) travelling waves emerge and (ii) the gap remains closed for some time until the two bodies separate. In the spatially semi-discretized framework, (i) travelling waves are not properly captured and this is commonly overcome via an additional impact law [3]. The impact law makes the finite-dimensional dynamics deterministic, meaning that any initial conditions determines a unique time evolution. In general, (ii) lasting contact phases are only possible for inelastic impact laws, resulting in energy dissipation. However, some very specific initial conditions generate solutions which have the property of being uniquely determined without requiring an impact law and featuring lasting contact phases. Such non-impulsive solutions "avoid" impacts by smoothly closing the gap with zero velocity and zero acceleration every time contact is activated. Among this class of trajectories, this work investigates periodic ones with one contact phase per period. They are referred to as 1 sticking phase per period (1 SPP or simply SPP) solutions of ndof autonomous systems with piecewise-linear dynamics and one unilateral contact constraint. While the investigation of such trajectories was initially motivated by their similarity with the continuous framework (no requirement for an impact law and dissipationfree lasting contact phases), they also seem to be limits of periodic solutions with k impacts per period as k → ∞ [14,16].
The dynamics of finite-dimensional systems subject to unilateral contact constraints has been widely investigated for several decades. Forced response of 1-dof vibro-impact oscillators was investigated [20] with particular interest in bifurcations [9,12]. Multidimensional systems were later studied, see, e.g., [11] and further references. The well-known phenomenon of chattering, consisting in an infinite number of gap closures in finite time, was also investigated for 1dof systems [2] and n-dof systems [1,4]. When chattering is said to be complete, it is followed by a phase of lasting closed contact, referred to as sticking, 1 where the external force maintains the gap closed for a finite amount of time. Such sticking phases have also been widely investigated [7,8,[17][18][19], including stability analysis [10]. The major difference in this work lies in the fact that autonomous periodic solutions with sticking phases are targeted: There is no external force. Such trajectories were initially investigated by Le Thi et al. [6] for a 2-dof spring-mass system with a unilateral constraint on one mass, where the terminology SPP was first introduced. The present contribution extends the results to n-degree-of-freedom systems, general mass and stiffness matrices such as those encountered in the finite element method, as well as any affine unilateral constraint.
More precisely, we consider a system with mass matrix M and stiffness matrix K , both symmetric positive-definite, undergoing a unilateral contact condition g(x) ≤ 0 where x is the vector of generalized coordinates. The gap function g is assumed affine, i.e., there exists a unit vector w such that g(x) = w x + g 0 where g 0 is the signed gap from the unconstrained resting position. Contact is closed if g(x) = 0. We now assume that M, K and w are fixed once for all and in the remaining, the term generically is to be understood in the sense of for almost every M, K and w.
Without loss of generality, the initial time t = 0 is fixed such that the free phase corresponds to time t ∈ [−2t 1 , 0] while the contact phase corresponds to t ∈ [0, 2t 2 ] (the coefficient 2 is introduced to simplify further derivations). The unknown period is T = 2t 1 + 2t 2 . The governing equations on [−2t 1 , 0] read complemented on [0, 2t 2 ], by the Signorini conditions [5]: Since periodic solutions are targeted, the unknowns of the problem are the initial conditions x(0),ẋ(0) ∈ R n , the durations t 1 , t 2 ∈ R and the normal contact force λ ∈ R. Inequality (2b) implies that a zero contact velocity w ẋ(0) = 0 may occur only if w ẍ(0) > 0 or if w ẍ(0) = 0. The first case corresponds to grazing contact and cannot lead to sticking because of inequality (refeq:nonsticking). The second case means that if the obstacle were to be removed, t = 0 would correspond to an inflection point of the normal acceleration w ẍ. A sticking phase can emerge with a zero contact velocity through this mechanism only. The free phase of the trajectory is connected to the contact phase via continuity conditions on x andẋ, together with the periodicity condition, the contact condition at 0 and 2t 2 and the zero contact velocity and acceleration equations: The schematic gap and normal contact force of a SPP are represented in Fig. 1.
Instead of solving these equations directly, the present approach consists in mimicking the contact phase by a free phase with a different stiffness matrix K . The new stiffness matrix should not affect the solutions of the original problem. That point is addressed in Sect. 2, leading to the explicit characterization ofK . Necessary conditions on t 1 , t 2 and x(0),ẋ(0) for the existence of solutions with 1 SPP are given in Sect. 3. The derivations are provided in Sect. 4. Summary of the results and illustrations are given in Sect. 5.

Unilateral contact condition seen as a stiffness matrix switch
We consider a solution of the free system Eq. (1) meeting the obstacle at t = 0, that is g(x(0)) = w x(0) + g 0 = 0, with a zero contact velocity w ẋ = 0. The aim is to describe the contact phase by a differential equation Mẍ +K x = 0 (S ) such that:

Consequences of energy conservation
The second condition means that w ẋ = 0 implies a constant energy. Equivalently, its time derivative must vanish: This scalar number must be zero for all x soẋ (K − K ) = 0, which must hold for anyẋ such that w ẋ = 0. Hence, each column of K −K must be orthogonal to the hyperplane orthogonal to w, that is, each column must be proportional to w. Accordingly, there exists a vector v such that K −K = wv . We are now going to use the first condition to characterize the vector v such that: We have just shown the following result.
Theorem 1 There exists a unique matrixK satisfying conditions C1 and C2, given by:

Stiffness matrix switch illustrated on a simple example
This abstract result can be easily understood with a simple system. Consider the 2-dof spring-mass system illustrated in Fig. 2 (left) as in [6]. The gap function is The mass and stiffness matrices are expressed as: From Theorem 1, the equivalent stiffness matrixK during the contact phase isK This is interpreted as follows: The first mass does not see any change as its dynamics is still dictated by the action of both springs. The dynamics of the second mass is not restricted by a unilateral contact condition since Mẍ +K x = 0 but instead the connected spring induces no action: As a consequence, the second row ofK vanishes. This can be understood as a "surgery" on the stiffness matrix, resulting in mimicking the contact phase by a free phase and thus eliminating the contact force as an unknown-which can be recovered from: It is noteworthy thatK is not symmetric, meaning it violates Maxwell-Betti reciprocal work theorem: It is not possible to build a linear (unconstrained!) mechanical system of stiffness matrixK .

Solution method
The problem is now reduced to finding one solution satisfying during the free phase and during the contact phase, together with the continuity conditions (3a), periodicity conditions (3b), zero contact velocity and acceleration (3c) and gap closure (3d) at t = 0. For the discussion to come, we introduce a terminology to distinguish actual solutions (admissible solutions) and solutions satisfying only the equalities, but not necessarily the two inequalities (potential solutions). (12a) is referred to as potential solution (or potential SPP). An admissible solution (or admissible SPP) is a solution of Eqs. (3), (11) and (12).

Definition 1 A solution of Eqs. (3), (11a) and
By construction of the solution, Eq. (2c) is always satisfied during the free phase since λ = 0 and Eq. (2b) is always satisfied during the contact phase, since g(x) = 0. This implies that Eq. (2d) is also verified. It is worth mentioning that Eq. (12b) follows from Eq. (10), We are now showing that the dynamics of both phases can be written in the formÿ = −L 2 y for some diagonalizable matrix L and y = M 1/2 x.

Free-flight dynamics
The matrix M has a unique symmetric and positivedefinite root M 1/2 ; let y = M 1/2 x. Equation (11a) becomes The matrix 3.2 Dynamics of the contacting phase Expressing Eq. (12) in terms of y gives: From Theorem 1 and introducing w 1 = M −1/2 K , it comes that so that finally: J 1 is the matrix of the orthogonal projection on the hyperplane w ⊥ 1 , so there is an orthogonal change of basis of matrix P such that P J 1 P is the orthogonal projection on the hyperplane H = e ⊥ n , of diagonal matrix J = diag(1, . . . , 1, 0). This change of basis transforms K 1 into K 2 = P K 1 P where K 2 is symmetric positive-definite. The restriction of J K 2 to H is symmetric positive-definite, hence diagonalizable, invertible and has positive eigenvalues. The direct sum of H and the kernel of J K 2 is the whole vector space; hence, J K 2 is diagonalizable and so is J 1 K 1 . Additionally, it follows that J K 2 has nonnegative eigenvalues (1 zero eigenvalue and n −1 strictly positive eigenvalues), implying that J 1 K 1 has a square root L 2 . In the end, we have shown that the dynamics during the contact phase can be described by the following ODE: In practice, L 1 and L 2 are found numerically.
3.3 Solutions ofÿ = −L 2 y the dynamics for some initial conditions y(0) andẏ (0) The dynamics of both the free and contact phases is captured by an ODE of the formÿ = −L 2 y (see Eqs. 14,18). This simple ODE can easily be solved as follows. Let .
3.4 Junction at t = 0 of the solutions for t < 0 and t > 0, and symmetry We consider one solution y ∈ C 2 (R) satisfyingÿ = −L 2 1 y over R − andÿ = −L 2 2 y over R + , both sharing the same initial conditions y(0),ẏ(0). Finding a periodic solution reduces to finding t 1 , t 2 > 0 such that there exists a nonzero solution y defined over [−2t 1 , 2t 2 ], with y(−2t 1 ) = y(2t 2 ) andẏ(−2t 1 ) = y(2t 2 ). Introducing A 1 , A 2 corresponding to L 1 , L 2 in Eq. (19), this condition can also be written as: . Nonzero solutions may only exist for t 1 , t 2 such that is singular. We are going to show the following result. Corollary 1 directly results from Theorem 2: if g 0 = 0, at the middle of the contact phase, that is g(x(t 2 )) = 0, then w x(t 2 ) = −g 0 which is in contradiction with y having a punctual symmetry over [0, 2t 2 ] (we would have y(t 2 ) = 0 so x(t 2 ) = 0). It follows from the theorem that any solution has two axes of symmetry. If g 0 = 0 and det(T 2 ) = 0, by central symmetry which is in contradiction with the unilateral contact condition w x ≥ 0.
To prove Theorem 2, we start with the following simple observation.
This ends the proof of Theorem 2.
Remark 1 For some very specific values of M and K , it is possible to find t 1 , t 2 such that det(T 2 ) = 0, leading to potential solutions with punctual symmetries instead of axial symmetries. Such an example is provided in Sect. 6.2, and however, these cases are never admissible: Corresponding λ has a point of symmetry, so cannot satisfy Eq. (12b).

Methodology to find SPPs
Put together, the previous derivations lead to a methodology to compute SPPs. We assume g 0 = 0. (t 1 , t 2 ) The condition of periodicity is governed by det( ) = 0. However, from Eq. 27 in the proof of Theorem 2, det( ) = 2 n det(R 2 ) det(T 2 ). From Corollary 1, det(T 2 ) = 0 does not lead to admissible solutions; hence, it suffices to find (t 1 , t 2 ) solution to det(R 2 ) = 0. This provides a first equation.

Periodicity as a first condition on
When the periodicity condition det(R 2 ) = 0 is satisfied, all the columns of adj(R 2 ) are colinear and are in the kernel of R 2 , hence in the kernel of , see Eq. (25). Choosing z(0) as a nonzero column of adj(R 2 ) ensures that det( )z(0) = 0, meaning that z(0) is proportional to a set of initial conditions leading to a potential SPP.
We insist that y(0) stores the n first rows of a nonzero column of adj(R 2 ) and is known: It is completely determined by (t 1 , t 2 ). Equation (33) hence provides a second equation on (t 1 , t 2 ), which is generically independent of the first one. Unknowns (t 1 , t 2 ) are thus solutions to two independent equations (except maybe for some peculiar M and K ), which yields the following important result. (t 1 , t 2 ), and thus isolated period T . If g 0 = 0 and there is no SPP.

Theorem 3 For generic M and K , if g 0 = 0 then SPPs may only occur for isolated couples
In particular, there is no continuum of SPP, contrary to trajectories with k impacts per period [14].
The case g 0 = 0 can be shown by observing that w x(0) = 0 yields an additional scalar condition on x(0), a priori independent of the other ones, resulting in an over-determined system, generically with no solution. A subgeneric example where M and K have been specially designed to induce an SPP with g 0 is provided in Sect. 6.1. Until then, it is assumed that g 0 = 0.
Gap closure at t = 0 We have just seen that an appropriate doublet (t 1 , t 2 ) determines z(0) = [y(0),ẏ(0)] chosen as a nonzero column of adj(R 2 ) , up to a multiplicative constant. This constant is fixed by the condition g(x(0)) = 0: w x(0) = −g 0 , that is: This sets the initial velocities tȯ The value g 0 = 0 does not affect the potential existence of SPPs, only dictated by (t 1 , t 2 ). It shows that generically, there is no non-trivial SPP with g 0 = 0. Recall that admissible SPPs must also satisfy the two following conditions, corresponding to Eqs. 11b and 12b. Remark 2 Whether the two above inequalities are satisfied depends on the sign of g 0 , but not on its magnitude. In particular, if there is a SPP for some (t 1 , t 2 ) and a g 0 = 0, there will be no SPP for the same (t 1 , t 2 ) and the opposite gap −g 0 .

Algorithm for finding SPPs
To find an admissible SPP: 1. ComputeK from (7), the positive-definite square roots L 1 of M −1 K and a square root L 2 of M −1K . 2. Find (t 1 , t 2 ) such that: where y(0) is a nonzero column of adj(R 2 ) . 3. For such (t 1 , t 2 ) and y(0), define the potential SPP via its initial conditions: 4. Check the admissibility of the potential SPP determined by (x 0 ,ẋ 0 ) via Eq. (22) and definition y = M 1/2 x: 5. If not admissible, return to step 2.
Animated SPPs are provided in ref. [15], for all the following examples as well as larger systems with up to 30 dofs. The mathematical results are first illustrated using the spring-mass system depicted in Fig. 3. For simplicity, all masses and stiffnesses are chosen equal to one, so that: The resting positions of the masses are zero so that positions and displacements are equal. The gap writes g 0 ≥ x 5 or g(x) ≥ 0 with g(x) = w x + g 0 and w = 0 0 0 0 −1 , and g 0 is arbitrarily chosen equal to one. Equation (7) yields which physically corresponds to disconnecting the last mass, only from the point of view of the last mass (see 2.3). Square roots L 1 and L 2 can be easily computed, and we can proceed with finding numerically (t 1 , t 2 ) solutions of Eq. 36, see Fig. 4. Instead of plotting the roots of det(R 2 ), in this figure each curve corresponds to (t 1 , t 2 ) such that a column c i of adj(R 2 ) is orthogonal to v = w M −1 K M −1/2 . This "1trick" reduces the computational cost and originates from the fact that when R 2 is singular, adj(R 2 ) is generally of rank one; the conditions vc i = 0 for i ∈ {1, . . . , n} become identical, meaning that solution points (t 1 , t 2 ) are located at the intersections of all the root curves

SPPs with slip and non-diagonal mass matrix
The presented results also apply to more complex geometries and general mass or stiffness matrices. In this subsection, an application is presented on the 4dof system depicted in Fig. 6. Contrary to the previous system with concentrated masses, here the mass and stiffness matrices are chosen as with x = x 1 x 2 y 2 y 3 corresponding to Fig. 6. We arbitrarily place the obstacle at position 4y = 1−x so that w = 0 −1 −4 0 / √ 17 and g 0 = 1/17. As previously, computingK , L 1 and L 2 is straightforward. Solutions of system Eq. (36) have to be found numerically, and then, inequalities 38 have to be verified (steps 3 and 4). Solutions of an admissible SPP are represented in Fig. 7 (left) together with the nonnegative gap and the nonnegative normal contact force (right). Again, the trajectories feature two axes of symmetry, at t = −t 1 and t = t 2 (Theorem 3).

Special behaviors for specific parameters
In this section, subgeneric cases are investigated, for completeness.
6.1 Closed gap at rest (g 0 = 0) Theorem 3 assumes g 0 = 0. We now investigate the case when g 0 = 0. The gap closure at t = 0 implies w x(0) = −g 0 = 0: non-trivial solutions can be found only if x(0) is orthogonal to w. This adds one independent scalar equation to the system (36), requiring an additional unknown. To illustrate SPPs with zero initial gap, we consider a 2-dof spring-mass system, with an unknown stiffness k 2 : Introducing an additional unknown k 2 (compared to 5.2, for example) compensate the additional equation w x(0) and solutions can be found numerically by solving: as illustrated in Fig. 8. The system of the previous subsection is considered. Solving for det(T 2 ) = 0 instead of det(R 2 ) = 0 provides (t 1 , t 2 ) leading to SPPs with points of symmetry in position, see Fig. 9. However, since on [0, 2t 2 ], λw and w M −1 K x have the same sign (see Eq. 10), and because x has a point of symmetry at t = t 2 , w M −1 K x is either zero on [0, 2t 2 ] or its sign  Like all SPPs with punctual symmetry, it is not admissible because λ is not nonnegative over [−2t 1 , 2t 2 ]. See [15] for the animated version changes, breaking condition 2c. In other words, there is no admissible SPP with axial symmetry. Moreover, because of the punctual symmetry of the positions, the obstacle cannot be still (except if g 0 = 0) and may only have a linear velocity during the contact phase. In Fig. 9, the motion of the obstacle has been smoothly extrapolated on [−2t 1 , 0] to satisfy periodicity.

Conclusion
Previous works have shown that autonomous periodic motions with one sticking phase per period (SPP) play a particular role in the frequency-energy plot of piecewise-linear dynamical systems with one unilateral constraint. Indeed, they seem to be the limit of autonomous periodic trajectories with k impacts per period (IPP) as k goes to infinity. Additionally, SPPs are uniquely determined by a set of initial conditions even though no impact law is specified. The minimal set of equations governing 1 SPPs was derived in the general n-dof framework with a single unilateral constraint. The first step is to solve a system of two scalar equations (one nonlinear and one linear) for the duration 2t 1 of the free-flight phase and the duration 2t 2 of the contact phase. As with k IPP solutions for which the sequence of free-flight durations determine appropriate initial conditions, these two durations then determine 1 SPPs. It was shown that if g 0 = 0, 1 SPP may only be isolated, irrespective of the sign of g 0 . The displacements of SPPs trajectories were shown to exhibit two axis of symmetry in time: one in the middle of the free flight (t = −t 1 ) and one in the middle of the contact phase (t = t 2 ). It was proven that 1 SPP occur only for discrete values of (t 1 , t 2 ). Results were illustrated using various geometries, unilateral constraints and mass matrices.
While it is clear that potential SPPs exist, there is no result on the existence of admissible SPPs.
This result is known [13], however it can be briefly proven.
Proof Let's first assume A is invertible. Then holds at least for any x which is not eigenvalue of A, hence for all x (it's a polynomial equality), in particular for x = 0. Similarly, if C is invertible, and it comes that which also holds for singular C.