Periodic solutions of n-dof 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-discretisation 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 k impacts 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.

(ii) the gap remains closed for some time until the two bodies separate. In the spatially semi-discretised 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 n-dof 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 dissipation-free lasting contact phases), they also seem to be limits of periodic solutions with k impacts per period as k ! 1 [15,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 were 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 1-dof 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 [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/ D w > x C g 0 where g 0 is the signed gap from the unconstrained resting position. Contact is closed if g.x/ D 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 D 0 is fixed such that the free phase corresponds to time t 2 OE 2t 1 ; 0 while the contact phase corresponds to t 2 OE0; 2t 2 (the coefficient 2 is introduced to simplify further derivations). The unknown period is T D 2t 1 C 2t 2 . The governing equations on OE 2t 1 ; 0 read complemented on OE0; 2t 2 , by the Signorini conditions [5]: Since periodic solutions are targeted, the unknowns of the problem are the initial conditions x.0/; P x.0/ 2 R n , the durations t 1 ; t 2 2 R and the normal contact force 2 R. Ineq. (1.2b) implies that a zero contact velocity w > P x.0/ D 0 may occur only if w > R x.0/ > 0 or if w > R x.0/ D 0. The first case corresponds to grazing contact and cannot lead to sticking because of ineq. (1.2c). The second case means that if the obstacle were to be removed, t D 0 would correspond to an inflection point of the normal acceleration w > R x. 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 P x, together with the periodicity condition, the contact condition at 0 and 2t 2 and the zero contact velocity and acceleration equations:

Unilateral contact condition seen as a stiffness matrix switch
We consider a solution of the free system eq. (1.1) meeting the obstacle at t D 0, that is g.
The aim is to describe the contact phase by a differential equation M R x C Q Kx D 0 (S) such that: C1.
[gap closure] all the solutions of S satisfy the contact condition during the contact phase: w > P x D 0 over OE0; 2t 2 ; C2. [energy conservation] the solutions of S which satisfy condition w > P x D 0 preserve the energy.

Consequences of energy conservation
The second condition means that w > P x D 0 implies a constant energy. Equivalently, its time derivative must vanish: This scalar number must be zero for all x so P x > .K Q K/ D 0, which must hold for any P x such that w > P x D 0. Hence, each column of K Q 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 Q K D wv > . We are now going to use the first condition to characterize the vector v such that:

Consequences of gap closure
We have just shown the following result.
Theorem 2.1 There exists a unique matrix Q K 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 g.
The mass and stiffness matrices are expressed as: From theorem 2.1, the equivalent stiffness matrix Q K during the contact phase is 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 R x C Q Kx D 0 but instead the connected spring induces no action: as a consequence, the second row of Q K vanishes. This can be understood as a "surgery" on the stiffness matrix, resulting in mimicking Contact phase viewed from first mass: unchanged.
Contact phase viewed from second mass: disconnected. the contact phase by a free phase and thus eliminating the contact force as an unknown-which can be recovered from: It is noteworthy that Q K is not symmetric, meaning it violates Maxwell-Betti reciprocal work theorem: it is not possible to build a linear (unconstrained!) mechanical system of stiffness matrix Q K.

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 (1.3a), periodicity conditions (1.3b), zero contact velocity and acceleration (1.3c) and gap closure (1.3d) at t D 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). By construction of the solution, eq. (1.2c) is always satisfied during the free phase since D 0 and eq. (1.2b) is always satisfied during the contact phase, since g.x/ D 0. This implies that eq. (1.2d) is also verified. It is worth mentioning that eq. (3.2b) follows from eq. (2.7), given that w > M 1 w > 0.
We are now showing that the dynamics of both phases can be written in the form R y D L 2 y for some diagonalisable matrix L and y D M 1=2 x.

Free-flight dynamics
The matrix M has a unique symmetric and positive-definite root M 1=2 ; let The matrix K 1 WD M 1=2 KM 1=2 is symmetric hence orthogonally diagonalisable; let K 1 D Q 2 Q 1 for some orthogonal matrix Q and D diag.! 1 ; : : : ; ! n / a diagonal matrix of positive diagonal entries.
3.2. Dynamics of the contacting phase Expressing eq. (3.2) in terms of y gives: From theorem 2.1 and introducing w 1 D 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 D e ? n , of diagonal matrix J D diag.1; : : : ; 1; 0/. This change of basis transforms K 1 into K 2 D P > K 1 P where K 2 is symmetric positive-definite. The restriction of JK 2 to H is symmetric positive-definite, hence diagonalisable, invertible and has positive eigenvalues. The direct sum of H and the kernel of JK 2 is the whole vector space, hence JK 2 is diagonalisable and so is J 1 K 1 . Additionally, it follows that JK 2 has non-negative 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.4.
Junction at t D 0 of the solutions for t < 0 and t > 0, and symmetry We consider one solution y 2 C 2 .R/ satisfying R y D L 2 1 y over R and R y D L 2 2 y over R C , both sharing the same initial conditions y.0/, P y.0/. Finding a periodic solution reduces to finding t 1 ; t 2 > 0 such that there exists a nonzero solution y defined over OE 2t 1 ; 2t 2 , with y. 2t 1 / D y.2t 2 / and P y. 2t 1 / D P y.2t 2 /. Introducing A 1 , A 2 corresponding to L 1 , L 2 in eq. (3.9), this condition can also be written as: with D exp. 2t 1 A 1 / exp.2t 2 A 2 /. Nonzero solutions may only exist for t 1 ; t 2 such that is singular. We are going to show the following result. Theorem 3.1 Instants t 1 , t 2 satisfy det D 0 iff there exists a nonzero solution y such that either y. t 1 / D y.t 2 / D 0; in that case, y has a point of symmetry at t 1 over OE 2t 1 ; 0, and at t 2 over OE0; 2t 2 . or P y. t 1 / D P y.t 2 / D 0; in that case, y has one axis of symmetry at t 1 on the interval OE 2t 1 ; 0, and another axis at t 2 on the interval OE0; 2t 2 . Corollary 3.1 All generic solutions have two axes of symmetry per period: one in the middle of the contact phase and one in the middle of the free flight phase.
Corollary 3.1 directly results from theorem 3.1: if g 0 ¤ 0, at the middle of the contact phase, that is g.x.t 2 // D 0, then w > x.t 2 / D g 0 which is in contradiction with y having a punctual symmetry over OE0; 2t 2 (we would have y.t 2 / D 0 so x.t 2 / D 0). It follows from the theorem that any solution has two axes of symmetry. If g 0 D 0 and det.T 2 / D 0, by central symmetry If for some t 0 2 . 2t 1 ; t 1 /, w > x.t 0 / > 0, then by punctual symmetry w > x. 2t 1 t 0 / < 0 with 2t 1 t 0 2 . t 1 ; 0/, which is in contradiction with the unilateral contact condition w > x 0.
To prove theorem 3.1, we start with the following simple observation.
Then, using matrix trigonometric double-angle identities such as sinc.P / cos.P / D sinc.2P / for any matrix P , can be factorized in two different ways: From the technical lemma A.1, it follows that det.T 1 / D det.R 2 / and det.T 2 / D det.R 1 /. In the end, the following equality holds: det./ D 2 n det.R 2 / det.T 2 /: (3.17) For to be singular, either det.R 2 / D 0: which, through Eq. (3.12), is equivalent to P y.t 1 / D P y.t 2 / D 0, that is y has two axial symmetries; or det.T 2 / D 0: which is equivalent to y.t 1 / D y.t 2 / D 0, that is y has two punctual symmetries. This ends the proof of Theorem 3.1.
Remark 3.1. For some very specific values of M and K, it is possible to find t 1 ; t 2 such that det.T 2 / D 0, leading to potential solutions with punctual symmetries instead of axial symmetries. Such an example is provided in section 6.2, however these cases are never admissible: corresponding has a point of symmetry, so cannot satisfy eq. (3.2b).

Methodology to find SPPs
Put together, the previous derivations lead to a methodology to compute SPPs. We assume g 0 ¤ 0.
Periodicity as a first condition on .t 1 ; t 2 / The condition of periodicity is governed by det./ D 0. However, from eq. (3.17) in the proof of Theorem 3.1, det./ D 2 n det.R 2 / det.T 2 /. From Corollary 3.1, det.T 2 / D 0 does not lead to admissible solutions, hence it suffices to find .t 1 ; t 2 / solution to det.R 2 / D 0. This provides a first equation.
When the periodicity condition det.R 2 / D 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. (3.15). Choosing z.0/ as a nonzero column of adj.R 2 / ensures that det./z.0/ D 0, meaning that z.0/ is proportional to a set of initial conditions leading to a potential SPP.
Smooth normal contact velocity at t D 0 as a second condition on .t 1 ; t 2 / The condition of zero pre-impact normal velocity is satisfied over OE0; t 2 by construction of Q K: on OE0; t 2 , hence w > P x is constant and by symmetry, it vanishes over OE0; t 2 . However eq. (1.3c) also includes w > R x.0/ D 0, which yields a condition on the free flight dynamics, since R x D M 1 Kx on OE 2t 1 ; 0 so w > M 1 Kx.0/ D 0. In terms of y D M 1=2 x, the equality also writes We insist that y.0/ stores the n first rows of a non-zero column of adj.R 2 / > and is known: it is completely determined by .t 1 ; t 2 /. Equation (4.2) 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. In particular, there is no continuum of SPP, contrary to trajectories with k impacts per period [15]. The case g 0 D 0 can be shown by observing that w > x.0/ D 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 section 6.1. Until then, it is assumed that g 0 ¤ 0.
Gap closure at t D 0 We have just seen that an appropriate doublet .t 1 ; t 2 / determines z.0/ > D OEy.0/; P y.0/ chosen as a non-zero column of adj.R 2 / > , up to a multiplicative constant. This constant is fixed by the condition g.x.0// D 0: w > x.0/ D g 0 , that is: This sets the initial velocities to 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 D 0.
Recall that admissible SPPs must also satisfy the two following conditions, corresponding to eqs. (3.1b) and (3.2b).
Remark 4.1. 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 . To find an admissible SPP: 1. Compute Q K from (2.4), the positive definite square roots L 1 of M 1 K and a square root L 2 of M 1 Q K. 2. Find .t 1 ; t 2 / such that:

Algorithm for finding SPPs
where y.0/ is a non-zero column of adj.R 2 / > . 3. For such .t 1 ; t 2 / and y.0/, define the potential SPP via its initial conditions:

SPPs of 5-dof spring-mass system
The mathematical results are first illustrated using the springmass system depicted in fig. 3. For simplicity, all masses and stiffnesses are chosen equal to one, so that:

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 4-dof system depicted in fig. 6. Contrary to the previous system with concentrated masses, here the mass and stiffness matrices are chosen as M D Time t 4y D 1 x so that w > D 0 1 4 0 = p 17 and g 0 D 1=17. As previously, computing Q K, L 1 and L 2 is straightforward. Solutions of system eq. (5.1) have to be found numerically, then ineqs. 5.3 have to be verified (steps 3 and 4). Solutions of an admissible SPP are represented in fig. 7 (left) together with the non-negative gap and the non-negative normal contact force (right). Again, the trajectories feature two axes of symmetry, at t D t 1 and t D t 2 . (theorem 4.1).

Special behaviours for specific parameters
In this section, subgeneric cases are investigated, for completeness.
6.1. Closed gap at rest (g 0 D 0) Theorem 4.1 assumes g 0 ¤ 0. We now investigate the case when g 0 D 0.
The gap closure at t D 0 implies w > x.0/ D g 0 D 0: non-trivial solutions can be found only if x.0/ is orthogonal to w. This adds one independent scalar equation to the system (5.1), 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.
Time t

Non-symmetric potential SPPs with moving obstacle
The system of the previous subsection is considered. Solving for det.T 2 / D 0 instead of det.R 2 / D 0 provides .t 1 ; t 2 / leading to SPPs with points of symmetry in position, see fig. 9. However, since on OE0; 2t 2 , w and w > M 1 Kx have the same sign (see eq. (2.7)), and because x has a point of symmetry at t D t 2 , w > M 1 Kx is either zero on OE0; 2t 2 or its sign changes, breaking condition 1.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 D 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 OE 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 SPPs 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 D t 1 ) and one in the middle of the contact phase (t D 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. Using the proposed methodology, computation of SPPs with 30 dofs was straightforward. No attempt with a larger number of dofs was made, but the equations are indeed expected to be challenging to solve numerically.
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.

A. Technical lemma
Proof. Let's first assume A is invertible. Then which also holds for singular C .