Nonsmooth modal analysis of an elastic bar subject to a unilateral contact constraint

This contribution proposes a numerical procedure capable of performing nonsmooth modal analysis (mode shapes and corresponding frequencies) of the autonomous wave equation defined on a finite one-dimensional domain with one end subject to a Dirichlet condition and the other end subject to a frictionless time-independent unilateral contact condition. Nonsmooth modes of vibration are defined as one-parameter continuous families of nonsmooth periodic orbits satisfying the local equation together with the linear and nonlinear boundary conditions. The analysis is performed using the wave finite element method, which is a shock-capturing finite volume method. As opposed to the traditional finite element method with time-stepping schemes, potentially discontinuous deformation, stress and velocity wave fronts induced by the unilateral contact condition are here accurately described, which is critical for seeking periodic orbits. Additionally, the proposed strategy introduces neither numerical dispersion nor artificial dissipation of energy, as required for modal analysis. As a consequence of the mixed time–space discretization, no impact law is needed for the well-posedness of the problem in line with the continuous framework. The frequency–energy dependency nonlinear spectrum of vibration, shown in the form of backbone curves, provides valuable insight on the dynamics. In contrast to the linear system (without the unilateral contact condition) whose modes of vibration are standing harmonic waves, the nonsmooth modes of vibrations are traveling waves stemming from the unilateral contact condition. It is also shown that the vibratory resonances of the periodically driven system with light structural damping are well predicted by nonsmooth modal analysis. Furthermore, the initially unstressed and prestressed configurations exhibit stiffening and softening behaviors, respectively, as expected.

general configurations, the solution should be approximated numerically. However, common numerical methods, such as Finite Element Methods (FEMs) combined with time-stepping schemes, exhibit limitations that might be unacceptable [15,1,31]: spurious oscillations commonly known as Gibb's phenomenon, dispersion and dissipation errors or chattering when constitutive impact laws are implemented. Also, most energypreserving numerical schemes dedicated to contact dynamics feature numerical dispersion issues [23].
In the framework of undamped linear continuous systems, natural modes of vibration are characterized by their natural frequencies (eigenvalues) and normalized shapes (eigenfunctions) [30]. In the phase space, linear modes can be seen as continuous families of periodic orbits describing two-dimensional planes. For nonlinear dynamical systems, the natural modes of vibration are known to form two-dimensional manifolds in the phase space tangent to the previously mentioned planes at equilibrium points [19]. These manifolds are invariant, i.e. for any initial condition on the manifold, the solution will remain on it as time advances [33]. Additionally, the nonlinear natural frequency of an orbit depends on the total energy of the latter.
Most investigations and developments on nonlinear modal analysis, conventionally within a finite dimensional framework or equivalent, have dealt with smooth nonlinearities [19], i.e. nonlinear terms differentiable with respect to the unknown of the problem, mostly polynomial functions of the displacement and velocity. The continuous mechanical system investigated in this contribution involves a unilateral contact condition stricto sensu. As a consequence, displacements, velocities and accelerations are not differentiable with respect to time and space in the usual sense. Such systems are named nonsmooth systems [2]. Classical tools of nonlinear modal analysis rely on smoothness and thus no longer apply. Accordingly, "nonlinear normal modes" (NNMs) of vibration for nonsmooth systems will be referred to as nonsmooth modes (NSMs), defined as one-parameter continuous families of periodic nonsmooth orbits.
The common method to approximate solutions of continuous systems described by partial differential equations is to first obtain, through a semi-discretization in space, a finitedimensional nonlinear Ordinary Differential Equation solely depending on time. Unilateral constraints can then be incorporated via regularization to enable the use of classical nonlinear modal analysis [7,29,6,35,9,18,21]; the main issue lies in the high computational effort induced by the contact stiffness. Another strategy consists in preserving the nonsmooth nature of the system. It can be handled via appropriate time-stepping nonsmooth numerical schemes [24] or semi-analytic derivations [25,37] involving energy-preserving impact laws. A few closed-form solutions in the form of periodic shock-waves of an elastic string vibrating against a point obstacle are re-2 ported using the method of characteristic lines [11,12]. This string system shares similarities with the bar considered in the present work for which a single continuous family of periodic solutions was exhibited through analytical considerations [5]. In the remainder, this family will be shown to form the first NSM of the system. The current work proposes an approach to numerically approximate any NSM.
By definition of NSMs, their computation relies on periodicity: the state of the system at time t 0 must be recovered at some t 0 C T where T is the period, requiring numerical methods that preserve total energy and are not prone to numerical dispersion. A considerable amount of research has been devoted to this topic with limited success [23,20,3]. In the present work, the Wave Finite Element Method (WFEM) [34] is implemented in order to numerically approximate the modes of vibration of a one-dimensional finite elastic bar subject to frictionless unilateral contact constraints. The WFEM is a shock-capturing method-similar to the Godunov method [16] commonly used in computation fluid dynamicswhich consists in finding and tracking waves propagating in a mechanical system. Such shock waves are expected because of the unilateral contact condition. The proposed method does not introduce artificial dissipation of energy nor numerical dispersion [41,34]. Multidimensional problems incorporating contact conditions are not targeted in this work due to the expected intricate interaction of compressional and shear waves.
The paper is organized as follows. The one-dimensional system of interest is outlined in section 2, where three different configurations are detailed: unstressed, prestressed and statically grazing. The capability of the WFEM to handle unilateral contact dynamics is illustrated in section 3 through a benchmark problem offering an exact solution. The overall strategy capable of performing nonsmooth modal analysis is described in sections 4 and 5. The approximate nonsmooth modes are then characterized in section 6. Comparison with periodically forced responses is undertaken in section 7. The implementation of WFEM is thoroughly detailed in Appendix A for a one-dimensional problem. The enforcement of contact constraints within WFEM is explained in Appendix B.

System of interest
The system of interest is a homogeneous elastic bar of length L and constant cross-sectional area S rigidly fixed to the ground at its left end and whose right end is subject to a conservative unilateral constraint as shown in of infinitely small transformations, the unknown displacement, velocity, strain and stress fields are denoted by u.x; t /, v.x; t /, ".x; t / and .x; t / respectively where x 2 OE0 I L is the coordinate of a point of the bar along its longitudinal axis in the initial configuration and t denotes time. The quantity r.t / is the unilateral contact force emerging at x D L. The mass per unit volume is denoted by > 0 and E > 0 stands for Young's modulus which are both, by assumption, space and time independent. Any elastic wave traveling within the bar thus propagates at constant velocity c D p E= . The signed distance between the right end of the bar and the obstacle, or gap function, is defined as g.u.L; t // D g 0 u.L; t /, where g 0 is the signed distance between the unconstrained resting position and the obstacle. In linear elasticity, the stresses read D E" where the axial strains " D u ;x should be physically admissible, that is u ;x > 1. Stresses are related to the contact force by .L; t / D Eu ;x .L; t / D r.t /=S . Unless stated otherwise, there is no other external force acting on the system, either per unit length or pointwise at the boundary. The full formulation reads: -Local equation: where . / ;t t stands for the second derivative in time and . / ;xx , for the second derivative in space. -Homogeneous fixed (Dirichlet) boundary condition (BC) at x D 0 and 8t > 0: -Signorini complementarity condition at x D L and 8t > 0: g.u.L; t // 0; r.t / Ä 0; r.t /g.u.L; t // D 0 -Initial conditions 8x 2 0 I LOE: where u 0 and v 0 are prescribed functions. This problem has a unique solution which conserves the total energy [32]. The local equation is the wave equation defined on a finite domain. D'Alembert's exact solution is the sum of backward and forward travelling waves propagating with velocity c and reads 8x 2 OE0 I L and 8t > 0 [17]: (5) where u ? 0 and v ? 0 are periodic extensions on the real axis of u 0 and v 0 defined on OE0 I L only. The period of u ? 0 and v ? 0 depends on the boundary conditions considered. Non-trivial periodic solutions of the contact problem described by Eqs. (1) to (4) are successions of free phases (open gap) and contact phases (closed gap). They can therefore be understood as the combination of solutions to an hyperbolic Partial Differential Equation with a switching boundary condition at x D L between vanishing stress when the gap is open to prescribed displacement when the gap is closed. The nonlinearity in the formulation comes from the fact that the time of switch depends on the solution. There is also a subtlety in the BC switch at x D L: an open gap implies u ;x .L; / D 0 while a closed gap implies a non-homogeneous prescribed Nonsmooth modal analysis of an elastic bar subject to a unilateral contact constraint 3 displacement of the form u.L; / D g 0 . Accordingly, when the gap is open, BCs at x D 0 and x D L are homogeneous and will be referred to as fixed-free BCs in the remainder. For a closed gap, the BC at x D 0 is homogeneous while the BC at x D L is non-homogeneous: they will be referred as fixed-fixed BCs.
The question is now to find the two functions u 0 and v 0 generating a time-periodic state of the system which satisfies the local equation (1) and properly switches between fixedfree and fixed-fixed BCs to comply with the unilateral contact conditions (3). This is achieved using the Wave Finite Element Method.

Illustration of the capabilities of WFEM
The WFEM is detailed, including a readily-implementable algorithm, in Appendices A (the WFEM method itself) and B (unilateral contact conditions in WFEM). In the sequel, the governing Eqs. (1) to (4) are equivalently recast in terms of stresses and velocities stacked in a state vector q. The system is discretized simultaneously in space and time with mesh sizes x D L=N and t D x=c respectively, where N is the number of cells used for space discretization. The overall dynamics is approximated by Q .n/ D AQ .n 1/ where Q .n/ D OEQ .n/ 1 : : : Q .n/ N and Q .n/ i stands for the averaged state q in cell C i at time t n D nt. The matrix A encompasses stiffness and inertial terms as well as the appropriate boundary conditions. Two matrices are constructed: A D A f for free-free BCs (open gap) and A D A c for fixed-free BCs (closed gap).
In this section, the capabilities of the WFEM are succinctly illustrated on a unilaterally constrained linearly elastic finite bar benchmark, for which analytical solutions are known [15]. It can be compared to other time-integration schemes based on a semi-discretization in space through the commonly used FEM [41]. The problem consists in an unclamped homogeneous elastic rod of length L and constant cross-sectional area S bouncing against an obstacle due to a distributed internal force, as shown in Fig. 2. Parameters reported in [15] are considered and listed in Tab. 1. The displacement (retrieved from stress via an integration in space), velocity, contact force an total energy is depicted for an isolated periodic solution in Fig. 3, for N D 100 cells. Total energy and periodicity are accurately preserved and the nonsmooth behaviour is well recovered. Curves in Fig. 3 and the exact curves are undistinguishable [41]. In contrast, traditional numerical schemes dedicated to contact dynamics and commonly based on the standard FEM show spurious oscillations on the contact displacement and stress [31]. Moreover, these oscillations do not  disappear when the time step decreases and might increase instead. Various solutions have been proposed to alleviate these difficulties, mainly by adapting the time-domain discretization. They all fail by either adding non-physical damping, or not strictly respecting the contact constraint. This is partly explained by the inability of standard finite elements to properly propagate information. Also, the implementation of an impact law yields unavoidable drawbacks, such as chattering if energy conservation is targeted [1]. None of these issues are present with WFEM that produces accurate approximations with a reasonable number of cells. This can be partly explained by 4 the fact that shock wave propagation is accurately captured by this scheme.

Periodic solutions
NSMs are regarded as continuous families of periodic solutions satisfying the local equations and the boundary conditions including the unilateral contact conditions. Finding a periodic solution translates into finding initial conditions u 0 .x/ and v 0 .x/ and a period T which generate a solution satisfying Eqs. (1) to (4) in conjunction with the periodicity conditions Without loss of generality, it is assumed that a period starts with a free phase at t D 0 and ends with a contact phase at t D T . In the general case, a pattern of free and contact phases will arise within one period. The k successive transition instants between free and contact phases are denoted by T i with 0 < T 1 < : : : < T k 1 < T k D T and are unknown. The sought solution is then an unknown combination of functions of the form (5) where u ? 0 and v ? 0 switch between fixed-free BCs over OE0 I T 1 , OET 2 I T 3 , etc. and fixed-fixed BCs over OET 1 I T 2 , OET 3 I T 4 , etc., which must be periodic in time. Connecting these portions of D'Alembert solutions in order to form a periodic solution is quite a formidable task.
The WFEM approximately solves the above problem by building two matrices, A f for the fixed-free BC and A c for the fixed-fixed BC, so that mapping an initial state to the current state after a succession of free and contact phases is straightforward via Eq. (40). In contrast to more commonlyused time-integration schemes based on FEM, the WFEM perfectly preserves energy (in this one-dimensional configuration at least), which is crucial for the computation of periodic solutions. Additionally, for the problem at hand, this scheme is not prone to numerical dispersion and does not suffer from the well-known Gibb's phenomenon commonly observed in unilateral contact dynamics [15].
The FEM framework without regularization of the complementarity conditions in continuous time is used for the calculation of nonsmooth modes in [25,37,38]. However, the formulation relies on a space-discretization and hence requires the incorporation of an impact law-for instance Newton's impact law-to ensure the well-posedness of the problem. Since energy conservation is necessary to find nonsmooth modes, a perfectly elastic impact law should be considered, from which would emanate impulsive dynamics and chattering. This excludes periodic motions with "lasting" non-impulsive contact phases similar to those of the periodically bouncing elastic bar exposed earlier. "Lasting" contact phases can only be obtained with a dissipative inelastic impact law. By simultaneously discretizing the governing equations in space and time in accordance with the characteristic lines, the WFEM accurately propagates shock waves. This probably explains why no impact law is needed in this framework; however, the authors are not aware of a formal proof of this assertion.

Nonsmooth modal analysis
Nonsmooth modal analysis characterizes vibratory mechanical systems with nonsmooth nonlinearities [24,25,38] by searching for one-parameter continuous families of periodic trajectories forming manifolds in the phase space. In this work, the targeted solutions are assumed to comprise free phases (open gap) as well as contact phases (closed gap), as described in section 4. An example of an admissible periodic solution with two free phases and two contact phases is depicted in Fig. 4.

Problem formulation
As stated previously, the motion is uniquely determined by the initial conditions. The aim is to find such initial conditions which generate a periodic motion in the spirit of shooting methods [4]. If such periodic solutions exist, the formulation reduces to finding their period T and the appropriate initial conditions satisfying q.x; T / D q.x; 0/; 8x 2 OE0 I L: In the discrete setting of the WFEM, this reduces to finding n T and the vector of initial conditions Q .0/ such that Q .n T / D Q .0/ where n T is the number of steps required to cover the period T . In their simplest incarnation, the sought periodic solutions are only composed of one free phase and one contact phase per period 1 even though more complicated patterns are expected to exist. Accordingly, the solution is assumed to be in open contact (free phase) for m consecutive steps and in closed contact (contact phase) for p consecutive steps. Invoking Eq. (40), the state of the system after n T D m C p time steps emanating from initial state Q .0/ reads with the origin of time taken at the beginning of the free phase. Both matrices A c and A f are known, see Appendix A. The duration of free and contact phases corresponding to periodic motions are not known a priori, so the integers m and p leading to acceptable solutions are unknowns of the problem. By enforcing the periodicity conditions, Eq. (8) simplifies to An initial condition Q .0/ satisfying Eq. (9) is called a potential solution for given p and m. Indeed, potential solutions may not be actual solutions of the initial problem: nothing prevents them from penetrating the rigid foundation during free phase and/or the corresponding contact force could be non-negative: these conditions cannot be enforced in Eq. (9). Consequently, potential solutions are called admissible solutions if they satisfy the following additional conditions: COND1 Free phase: u .n/ N C1=2 6 g 0 and r .n/ D 0 for n D 0; 1; : : : ; m. COND2 Contact phase: u .n/ N C1=2 D g 0 and r .n/ 6 0 for n D m C 1; m C 2; : : : ; m C p. where u .n/ N C1=2 u.L; T / and " .n/ i is the discretized strain field, both retrieved from discretized stresses. To summarize, periodic solutions satisfying contact conditions are obtained by solving the problem: Find m 2 N , p 2 N , and initial condition Q .0/ 2 R 2N such that Eq. (9) together with COND1, COND2 and COND3 are simultaneously satisfied.

Solution procedure
For specified m and p, a periodic motion necessarily corresponds to an initial condition which is in the kernel of the operator S T D A p c A m f I as highlighted in Eq. (9), that is Q .0/ 2 ker S T . The dimension h of ker S T depends on the doublet .m; p/. A non-trivial solution may exist only if h > 0. If h > 0 and fe 1 ; e 2 ; : : : ; e h g is a basis of ker S T , then Q .0/ D˛1e 1 C˛2e 2 C : : : C˛he h for some˛1;˛2; : : : ;˛h. The state Q .n/ is completely determined by Q .0/ via Eq. (46). Accordingly, Q .n/ can be expressed in terms of the coefficients 1 ;˛2; : : : ;˛h only and it suffices to find appropriate values such that COND1, COND2 and COND3 are satisfied.
Finding m, p such that h > 0 is achieved by systematically computing h for every combination of m and p within a given range. Once an admissible solution is known, other admissible solutions can be straightforwardly found in its vicinity. As soon as a family of periodic orbits-involving one contact phase and one free phase-emerges, it defines a nonsmooth mode of vibration.

Results and discussion
Nonsmooth modes of vibration of the system described in section 2 are now constructed. The parameters E, and L are arbitrarily chosen to be unity and units are discarded. The results are obtained for 1000 cells and time step is calculated accordingly, i.e. t D x=c D 1=1000.
In the linear framework (that is without the contact conditions), linear natural frequencies of vibration denoted ! k for fixed-free boundary conditions and k for fixed-fixed boundary conditions are It is expected that the slightly damped linear system under consideration periodically forced in a neighborhood of these frequencies will resonate. Conversely, in the nonlinear framework (that is with contact conditions), such frequencies do not have the same physical interpretation anymore: the nonlinear system will resonate close to other unknown frequencies which are calculated by the nonsmooth modal analysis, for a given level of energy. The nonlinear system can also feature internal resonances where two or more NSMs interact [22]; for instance, driving the system in the neighborhood of a highfrequency NSM may activate a large amplitude low-frequency NSM [19]. Note that the system of interest satisfies a complete internal resonance condition in the sense that ! k and k , k D 2; : : : ; 1 are all multiples of ! 1 and 1 , respectively. Moreover, a nonlinear system may exhibit subharmonic and superharmonic resonances in the vicinity of .p! k /=q, for p, q positive integers such that 0 < p=q < 1 and 1 < p=q, respectively [19].
In this work, Frequency-Energy plots are preferably used, quoting [19]: "a nonlinear modal motion is represented by a point in such plots, which are drawn at a frequency corresponding to the minimal period of the periodic motion and at an energy equal to the conserved total energy during the periodic motion. A branch, represented by a solid line, is a family of nonlinear modal motions possessing the same qualitative features". The terminology "branch" and "backbone curve" is used interchangeably in the remainder. The reported frequencies are normalized with respect to ! 1 and the energy is normalized with respect to the energy of the first linear mode grazing orbit satisfying max t 2R u.L; t / D g 0 : When the calculated periodic solutions have one free phase and one contact phase, admissible solutions are only found when h D 1 for various .m; p/. For h > 2, potential solutions failed to satisfy COND1 and COND2 simultaneously, hence the distinction between potential and admissible. Admissible solutions were found only when h equals the number of contact phases. Additionally, recent developments, out of the scope of the present paper, have shown that if a periodic trajectory belongs to a continuum parameterized by T , then the displacement must be piecewise-linear with respect to time and space. Thus, if a periodic solution does not satisfy this property, it is either isolated or part of a continuum of orbits sharing the same period. A previous investigation on the first main NSM [5] showed that the duration of the contact phase pt is related linearly to the duration of free phase mt when a solution belongs to a family of periodic orbits. Numerical experiments suggest that other branches could present similar relations, which are employed to bound the range where m and p are iterated for each potential backbone curve. Once a solution is found, a point in a NSM branch is defined by the total energy of the admissible periodic motion and the respective frequency ! D 2 =T , see Fig. 5.
If the initial gap is positive g 0 > 0 (unstressed bar at rest), NSM branches arise in the vicinity of the linear natural frequencies ! k and subharmonics p! k =q for p; q 2 N , and present a hardening behavior. For a "negative" gap g 0 < 0 (prestressed bar at rest), the NSM branches start in the vicinity of the linear natural frequencies k and display a softening behavior. Finally, if the bar is statically grazing with the rigid foundation (grazing bar at rest, that is g 0 D 0), then the NSM branches are the vertical asymptotes common to the previous backbone curves, as shown in Figs. 5 and 6 for the first and second NSM, respectively. Internal resonance branches, explored in subsection 6.3, are omitted in these figures. Identical results were already computed for the first main NSM only [5]. Also, a similar behavior has been observed in the frequency-energy plot of simplified discrete models for the unstressed case with penalized contact constraints [7,29] and with a nonsmooth formulation [24].
The frequency-energy dependence exhibited by NSM is a typical feature of nonlinear systems [19]. Interestingly, when g 0 D 0 the system mimics a linear behavior where the frequency and the shape of the orbits do not depend on their energy. However, the behaviour remains globally nonlinear because of the BC switches.
Though the behavior and ranges of frequencies of the NSM branches depend on the sign of g 0 , they "converge" to the branches of g 0 D 0 when g 0 ! 0, as seen in Fig. 7 for the main NSM branch in OE! 1 ; 1 . In other words, there is no difference between NSMs with very small positive initial gap, zero initial gap or very small initial negative gap. Moreover, the normalized shape of the backbone curves depend only on the sign of g 0 . Without loss of generality and for illustration purposes, g 0 D˙10 3 or g 0 D 0 in the remainder.
In the following sections, a detailed study of the NSM branches is presented. It should be noted that all the results were obtained numerically. The proposed results have to be understood as conjectures.

Main backbone curves
The main backbone curves are defined as the NSM branches, in the Frequency-Energy plot, starting in the vicinity of: -! k , the k-th linear natural frequency for fixed-free BC, for positive initial gap g 0 > 0, k , the k-th linear natural frequency for fixed-fixed BC, for negative initial gap g 0 < 0. -For g 0 D 0, the main backbone curve is a vertical line asymptotic to its two above counterparts. Such asymptotes are located at ! 1 .! k =! 1 C 1/ 2 =.! k =! 1 C 2/. Periodic motions corresponding to the main backbone curves in the range OE! 1 ; 1 are depicted in Fig. 8. Such piecewiselinear displacements were already observed as isolated periodic solutions, without describing a continuum [11,12]. However, previous numerical investigations to locate families of periodic orbits were inaccurate because of numerical dispersion stemming from the space semi-discretization [29,24].
A generic space-time plot of the displacement of the bar is provided in Fig. 9. The solution is the combination of Every periodic solution corresponding to the main backbone curve comprises exactly one free phase and one contact phase. For a positive as well as negative initial gap, the contact phase duration increases while the free phase duration decreases as the frequency of vibration increases. When the initial gap is zero, periodic motions have constant period equal to T D 3L=c and the duration of contact closure is L=c while the free phase time is 2L=c. It should also be noted that the contact force associated to all the depicted motions is piecewise constant on a period.
The invariant manifold in the unstressed case is plotted in the cross-section .u N C1=2 ; v N C1=2 ; u 3=2 /-i.e. (displacement of the contacting end, velocity of the contacting end, displacement of the right-hand side interface of cell C 1 )-of the state-space in Fig. 11. The linear portion of the first mode corresponds to an ellipse (like essentially all linear modes). The linear and nonlinear portions of the manifold are not con- 8 tinuously connected because of the internal resonance property. This is visualized in Fig. 12 which shows the grazing periodic displacements for the first linear mode and first main NSM, both having frequency ! 1 . It should be understood that any autonomous solution of frequency ! 1 shall be expressed as an infinite Fourier sequence separated in space and time of the form where ! k D .2k 1/ c=.2L/ D .2k 1/! 1 , k 2 N , are the linear natural frequencies; the coefficients a k and b k are computed from the initial conditions [17] and are not of interest here. In other words, the full internal resonance condition enjoyed by the system is such that the latter exhibits infinitely many solutions of frequency ! 1 , the piecewise-linear grazing solution being one of them. Accordingly, a discontinuity in the manifold is made possible. The above internal resonance condition can be annihilated by changing the BC at x D 0 from Dirichlet type to Robin type saying that the bar is attached to a spring of stiffness ı at x D 0: ES u ;x .0; t / ıu.0; t / D 0. The solution can still be expressed as an infinite Fourier sequence similar to Eq. (12) with k .x/ D p k cos.! k x=c/ C sin.! k x=c/ where p k D cot.! k L=c/. The ! k are now solutions to the transcendental equation ES! D cı cot.!L=c/ and are no longer commensurate [30]. Assuming that the grazing solution of the first NSM has frequency ! 1 (this should be proven!), it will necessarily see the sole contribution of the first linear mode, which is the only motion of frequency ! 1 in sequence (13), that is: This necessarily induces a continuous connection between linear and nonlinear grazing trajectories. To summarize, the above argument stipulates that the internal resonance condition is necessary to see a discontinuity in the modal manifold between the linear and the nonlinear portions. However, this is not a sufficient condition. Figure 13 depicts the generic space-time plot (irrespective of g 0 ) of the second main backbone curve in the range OE! 2 ; 2 showing that the periodic solutions have one "node of vibration" along x-similar to the second linear mode of the fixed-free bar [17], depicted in Fig. 14. The displacement u.L; / also features one contact phase per period, by assumption. Furthermore, for g 0 D 0, the period is 5L=.4c/ comprising one contact phase of duration L=.4c/ and one free phase of duration L=c, irrespective of the energy of the system.
For illustration purposes, only the first two main NSMs are shown in this section. However, the presented approach can be also used to find higher frequency NSM which requires a smaller time-step t and therefore, higher computational effort. Periodic solutions corresponding to the main backbone curves for higher frequencies display a similar structure: one free phase, one contact phase, and nodes of vibration along x. ! D 2k.n C 1/ 1 n C 1 ! 1 ; n; k 2 N : For g 0 < 0, a subharmonic backbone curve starts in the left vicinity of k and shares an asymptote with the subharmonic backbone curve for positive initial gap. Furthermore, for g 0 D 0, a subharmonic backbone curve is asymptotic to its above counterparts. Several numerical experiments suggest that the asymptotes (vertical brown lines in Figs. 5 and 6) are located at The number of computable subharmonics is determined by the discretization time-step; for example, Figs. 5 and 6 display only six branches while the time step was t D 10 3 . Four periodic solutions corresponding to subharmonic backbone curves are depicted in Fig. 15 for g 0 > 0. These plots suggests that the periodic motions present n grazing instants together with one contact phase, when the NSM branch starts in the vicinity of the n-th subharmonic of ! k , as inferred from Eq. (15). The displacement of the contacting end for points "d " and "g" (see Fig. 6) are similar, however the latter has higher frequency. Such feature is also expected in periodic motions around subharmonics of higher NSM.
Moreover, periodic solutions belonging to distinct NSM branches coexist when the bar is prestressed (g 0 < 0), see point "h" in Fig. 6. The corresponding periodic displacements of the contacting end are plotted in Fig. 16 where the period of vibration is the same. Here, the solutions have longer free phases and increasing number of grazing instants with higher energy. These grazing instants coincide with the BC switching time of lower energy solutions.
Interestingly, extensive numerical experiments show that NSM branches do not appear in the range of frequencies OE k ; ! kC1 , k 2 N ; the reason for this is unknown.

Internal resonance branches
The intrinsic discrete nature of the employed numerical algorithm is such that WFEM cannot capture an actual continuum of solutions (see section 5). However, the periodic solutions previously presented all belong to numerically "apparent" continua parametrized by the energy where the solutions possess similar qualitative and quantitative features. The proposed approach also leads to sparse periodic solutions, see clouds of points in Fig. 17. Presumably, these points are also part of a complicated network of backbone curves stemming from main and subharmonic branches, however they are much more challenging to capture numerically because of the high-frequency content in the corresponding mode shapes, see Fig. 18. These points seem to be organized on internal resonances backbone curves [19] emanating from the main one in the vicinity of the frequencies ! D 2k! 1 =.2k 1/; k D 2; 3; : : : Fig. 15: Periodic motions in the vicinity of corresponding respectively to points "d ", "e", "f " and "g" in Figs. 5 and 6.
Moreover, it is observed in Fig. 17 that there are no internal resonance branches in the frequency ranges OE.2k 1/! 1 =.2k 2/I 2k! 1 =.2k 1/, k D 3; 4; : : : This is not further explored in this work. Periodic motions featuring an internal resonance in the vicinity of .! 3 C ! 1 /=5 D 6! 1 =5 are displayed in Fig. 18. These solutions present two grazing instants and appear to correspond to a modal interaction between the first and third NSM. In contrast to solutions on the main and subhar-  It is difficult to compute the occurrence of internal resonances in systems involving unilateral contact. However, the data obtained by WFEM give suggestions about the existence of these NSM branches. The characterization of the periodic motions associated to internal resonances is helpful to predict the possible sudden resonance of real life applications, when vibrating around frequencies defined in Eq. (17).

Periodic solutions with two contact phases per period
The previous sections were devoted to solutions with one contact phase and one free phase per period, as expressed in Eq. (8). In this section, solutions with more than one contact phase per period are explored. The methodology is unchanged, but two additional unknowns (duration of second contact phase`t and duration of second free phase mt ,à nd m 2 N ) arise in the new equation to be compared with Eq. The solution procedure explained in subsection 5.2 leads to a one-dimensional continuum of periodic orbits emerging in the vicinity of ! 2 =5 D 3! 1 =5, where ! 2 D 3! 1 , defining a subharmonic backbone curve with two contact phases, as illustrated in Fig. 19. This family of periodic orbits appears only for g 0 > 0 and presents hardening behavior. The corresponding space-time displacement field is shown in Fig. 20 for an arbitrary frequency. The asymptote of this NSM branch appears to be located at 2! 1 =3.
The contacting end, shown in Fig. 20 features one free phase with two grazing instants similar to the first subharmonic of the third NSM at ! 3 =3. Also, as in the case with one contact phase per period, the amplitude of the displacements and the duration of the contact phases get larger with the frequency of vibration. Interestingly, it resembles the two impact-per-period trajectories of a serial spring-mass system constrained by an obstacle with a purely elastic contact law [37]. In particular, the solutions have two axes of symmetry per period along the time axis, which are located in the middle of each free phase. A finer discretization is required to obtain additional NSM branches with higher number of contact phases. These branches might be of interest for the prediction of vibratory responses where the elastic bar resonates for excitation frequencies much lower than the first linear natural frequency ! 1 .

Stability analysis of periodic solutions
Due to the space-time coupled discretization and the conditional switching of boundary conditions, a rigorous stability analysis of the periodic orbits is a challenging task and remains an open problem left to future investigations. However, as a preliminary insight, the effects of a small sinusoidal perturbation on the initial conditions corresponding to the first main NSM were investigated by comparing the unperturbed and perturbed solutions time-integrated via WFEM over 1000 periods of the unperturbed solution. We found that the perturbed motion remains in the vicinity of the periodic motion, thus suggesting that a fraction of first main NSM motions are orbitally stable. This is illustrated in Fig. 21 which compares the unperturbed and perturbed motions, over the last period of integration t 2 OE999T a I 1000T a , corresponding to point "a" and g 0 > 0 in Fig. 5. 7 Forced response of a mechanical system subject to contact constraints One important purpose of performing nonsmooth modal analysis of a mechanical system is to predict its behavior when periodically forced [38,29]. Accordingly, the frequency response of the previously-investigated elastic bar with periodic external excitation is now compared to the NSM spectrum. The excitation consists in an external distributed harmonic force or a harmonically moving rigid wall that compresses the bar, see Fig 22. The force acting on the bar is f .x; t / D f 0 .x/ sin.!t / and the displacement of the moving wall is w.t / D w 0 sin.!t / where ! is the frequency of excitation and w 0 > g 0 . The gap function is now defined as g.u.L; t /; w.t // D g 0 C w.t / u.L; t /. The system has been slightly damped by adding a velocity-dependent term with a small viscous damping coefficient in the left-hand side of Eq. (1). The forced response of this system is obtained for various amounts of damping and was computed using the WFEM version detailed in Algorithm 2 using a time-stepping approach not specifically targeting periodic motions. Periodically moving wall The first tested configuration is f 0 D 0 and w 0 ¤ 0. The total energy of the steady-state solution averaged over one forcing period for increasing frequencies of excitation and various amounts of damping is shown in Figs. 23, 24 and 25 (top) for a positive, zero and negative initial gap, respectively. In the frequency ranges where NSM branches do not exist, WFEM could not find externally forced periodic steady-states. Instead, quasiperiodic or chaotic forced responses were detected. Calculating steady-state for each frequency requires long computational times which complicates the construction of a detailed forced response.
Periodic distributed force The second tested configuration is f 0 ¤ 0 and w 0 D 0, corresponding to the elastic bar excited by a distributed force. The results for a positive, zero and negative initial gap are depicted in Figs. 23, 24 and 25 (bottom) respectively.
For a positive gap and high frequencies, large damping confines the system to linear operating conditions. This is one of the reasons why the corresponding plots differ greatly. Moreover, similarly to the first configuration, quasiperiodic and chaotic solutions are observed in the frequency intervals where NSM branches do not exist. Presumably, there exist infinitely many backbones curves and corresponding "inbetween" regions with no periodic solutions; only a few are plotted in Figs. 23 to 25.
More importantly, the backbone curves obtained with the nonsmooth modal analysis provide an excellent approximation of the response resonances. In addition, the internal resonances of the system produce small protuberances in the forced response that coincide with the main resonance and subharmonic backbone curves.
The main advantage of the nonsmooth modal analysis is the characterization of the vibratory response without relying on very expensive numerical time-integration procedures: via NSM as presented, the prediction of frequencies of excitation at which the system will vibrate with high energy is straightforward.

Conclusions
Families of periodic orbits (known as nonsmooth modes of vibration) of an autonomous finite elastic bar subject to frictionless unilateral contact are investigated in this work. Three cases were explored: unstressed (g 0 > 0), prestressed (g 0 < 0) and zero initial gap (g 0 D 0). The computation of periodic solutions was achieved using the Wave Finite Element Method (WFEM), chosen because it preserves energy and avoids numerical dispersion. This method consists in discretizing simultaneously in time and space the governing dynamic equations, resulting in a simple matrix form. Then, the problem of finding periodic solutions was formulated as finding a vector in the kernel of a matrix supplemented by complementarity conditions enforcing unilateral contact constraints. The presented methodology can be adapted to multiple contact phases per period or to systems coupled by unilateral contact conditions. However, it is presently limited to simplified onedimensional problems with a single point of contact on the contact boundary.
It is shown that the elastic bar with unilateral contact has a rich dynamical behavior involving subharmonic resonances and internal resonances. Similar results are already reported in the literature [7,29]. However, the proposed methodology does not necessitate regularized contact conditions. In this work, the unilateral contact conditions are treated as a switch between Dirichlet and Neumann-type boundary conditions when the gap opens or closes: from fixed-free BCs (no contact) to fixed-fixed BCs (contact) and vice-versa. The found nonlinear periodic solutions lying on a NSM are combinations of travelling waves with discontinuous wave fronts, as opposed to their linear counterpart (without contact conditions) which are standing harmonic waves. Also, the behavior of the NSMs depends on the gap at rest: hardening branches for g 0 > 0, softening branches for g 0 < 0 and discrete spectrum for g 0 D 0. Such behavior was already observed only for the first NSM [5]. Moreover, the frequency ranges at which NSM branches exist are conjectured. Quasi-closed form solutions can be extracted from the provided results. They could act as benchmark solutions for researchers designing advanced numerical schemes in unilateral contact dynamics. NSMs stability and the role of coexisting periodic solutions (same energy and frequency) needs be further explored. The relevance of nonsmooth modal analysis was illustrated by the accurate prediction of periodically forced responses, irrespective of the way the external force is applied.
Extension of the proposed approach to multidimensional settings with several points of contact is a remarkable challenge as the non-dispersive and energy-preserving properties of the WFEM in the one-dimensional context will then be lost. However, such extension seems feasible with Godunov-type discretization and shooting techniques by slightly relaxing the imposed periodicity conditions.
Acknowledgments CY gratefully acknowledges the financial support of SENESCYT-Government of Ecuador and MEDA-McGill University fellowships. AT and ML acknowledge the financial support of the grants: 421542-2012 NSERC Discovery Grant and 2014-NC-173113 FRQNT Établissement de Nouveaux Chercheurs Universitaires. ML would like to thank Cyril Touzé for mentioning references [11] and [12].

A Description of the Wave Finite Element Method
In this section, the WFEM, introduced by Shorr for the simulation of shock wave propagation in solids [34], is thoroughly described.

A.1 Hyperbolic system of conservation laws
The local equation (1) can be equivalently written as a system of two first order partial differential equations in terms of velocities v.x; t/ and stresses .x; t/ ;t where . / ;t is the derivation in time and . / ;x is the derivation in space of quantity . / [10]. Recall that axial strains " D u ;x are bounded by non physical inter-penetration, which translates into u ;x > 1, see section 2. Displacements can be straightforwardly recovered by spaceintegration of the strains u.x; / D R x 0 u ;s .s; / ds u.0; /, where u.0; / D 0. By posing q D OE v > , Eq. (19) can be recast as The eigenvalues of matrix B are 1 D p E= and 2 D p E= , coinciding with the algebraic propagation velocity of the elastic wave: positive and negative for the two waves propagating in opposite directions. Since both eigenvalues are distinct and real, Eq. (20) is also referred to as a hyperbolic system of conservation laws [10].
Equation (20) involves time and space derivatives of q. However, observing that q ;t C Bq ;x D 0 is a local form of the conservation law of q (implying q. ; t/ can only change due to fluxes at the boundaries) corresponding to the following integral form it appears that the condition on the smoothness of q is no longer required. Therefore, q is allowed to exhibit discontinuities in time and space [14].

A.2 Discretization
The WFEM consists in dividing the spatial and temporal domain into grid cells of equal size and keeping track of an approximation to the integral of q within every single cell. As depicted in Fig. 26, the bar is discretized using a uniform grid of N cells. Each i th cell, denoted by C i , is delimited by the interval .x i 1=2 ; x i C1=2 /. Similarly, time is discretized into intervals of equal length t D t nC1 t n . The average of q. ; t/, over the ith cell at time t n is This expression is now employed to develop an explicit time-stepping algorithm where Q .nC1/ i is approximated by a function of Q .n/ i . Equation (23) is integrated between t n and t nC1 yielding Rearranging and dividing by x leads to This equation describes how the cell average should be updated within a time step in order to satisfy the conservation of q. In general, the two integrals involving Bq on the right-hand side of the equation cannot be evaluated exactly. Following [26], we pose and Eq. (25) simply becomes The next subsection is dedicated to the computation of F .n/ i˙1=2 , which are the time-averaged fluxes at Fig. 27: Reconstruction of q.x; t n / from the average fluxes Q .n/ i .

A.3 Approximation of the time-averaged fluxes
To approximate the fluxes at the interfaces defined by Eq. (26), the state q.x; t n / at time t n is assumed to be a piecewise constant function defined for all x, constructed from the cell averages Q .n/ i as depicted in Fig. 27. This piecewise reconstruction of the function q.x; t n / is identical to the Godunov's approach widely employed in computational fluid dynamics [27]. A suitable approximation of the flux F .n/ i C1=2 can be obtained by solving the problem, either numerically or exactly, of the conservation law Eq. (20) together with the following discontinuous conditions at time t n [27]: which constitutes a Riemann problem centered at x i C1=2 between cells C i and C i C1 [27]. The solution to this Riemann problem consists of two shock waves propagating along the characteristic lines x D˙ct, one moving to the left into cell C i and one moving to the right into cell C i C1 as depicted in a space-time plot in Fig. 28. The shock wave traveling to the left, indicated by W 1 , propagates at velocity s 1 and connects the state Q .n/ i and the interior state Q i generated by such shock wave. Moreover, the solution t 7 ! q.x iC1=2 ; t/ is constant over the time interval OEt n ; t nC1 . The Rankine-Hugoniot jump condition is proven to hold across any propagating discontinuity [36] which can be written for the left wave W 1 propagating at velocity s 1 and the right wave W 2 propagating at velocity s 2 : Nonsmooth modal analysis of an elastic bar subject to a unilateral contact constraint 15 Because of material continuity, cells C i and C i C1 cannot separate. This requires that the interior states must be equal across the material interface, Q i C1=2 D Q i C1 D Q i . By knowing that the shock speeds s 1 D s 2 D p E= D c are known and constants, the intermediate state is approximated with q.x iC1=2 ; t / Q i C1=2 such that t n 6 t 6 t nC1 and can be calculated from Eq. (29): Equation (30) is regarded as the exact solution of the Riemann problem involving linear elastodynamics [26,8]. The flux approximation in Eq. (26) can be calculated with the solution of a Riemann problem at the cell interface as In a nonlinear framework for the local equation, a solution to the Riemann problem should be approximated numerically using Riemann solvers [39].

A.4 Formulation for inner grid cells
Inserting Eq. (31) into Eq. (27) produces the iterative scheme Equation (32) describes the evolution in time of the states of the grid cells C i . This subsection provides the formulation for the inner cells, where i D 2; : : : ; N 1. The boundary cells C 1 and C N require a different treatment explained in the next subsection. Expressing the flux approximation, employing Eq. (31) on the right side of an inner cell yields Performing the same operation on the left side of the cell reads Accordingly, the total flux within an inner cell is the quantity BQ i C1=2 BQ i 1=2 which, when substituted into Eq. (32) yields Since the exact solution of a Riemann problem is being used, WFEM incorporates an appropriate time-step t D x=c. Then, the stress and velocity of a grid cell C i at time t nC1 are calculated as The above strong assumption is suitable only for 1D elastodynamics problems, since the wave velocity and the direction of the propagation is known. Also, such assumption enforces energy conservation and eliminates numerical dissipation [34]. In the multidimensional framework, even though the waves velocities are known, the waves could propagate in various direction throughout the physical domain. Equation (36) provides the main equation of the WFEM and characterizes how the average value Q .n/ i of q in an inner cell C i is updated at each time step. As required by the local conservation law (21) resulting from the absence of body forces, the evolution of the state of the inner cells depends only on the values of the adjacent cells. WFEM can be seen as the transference of the whole information embedded in cell C i to its adjacent cells at each time step. Employing the latter approach to obtain the evolution of cell states, involving discontinuities such as shock and rarefaction waves, is well known by the Fluid Mechanics community employing Finite Volume Methods [27].

A.5 Formulation for boundary grid cells
To compute the state of the boundary grid cells, the computational domain is extended by including additional cells on both boundaries, known as ghost cells [26], whose average values depend on the boundary conditions. This concept is taken from the Finite Volume Methods. Figure 29 depicts ghost cells for a system discretized using N cells. at ghost cell C 0 : at ghost cell C N C1 : Such average values coincide with the theory of reflection of elastic waves from fixed and free boundaries [17], which states that stress waves reflect from a fixed boundary with the same sign and from a free boundary with the changed sign; similarly, velocity waves reflect from a fixed boundary with an opposite sign and from a free boundary with the same sign. The evolution of the average values of the boundary cells, C 1 and C N , can be calculated by introducing the average values of Eq. (37) into Eq. (36), yielding for C 1 : and for C N : Compute stress and velocity in cell C N via Eq. (39) end Output: stresses, velocities at instants t 0 ; : : : ; t nT framework of linear elastodynamics, the propagation at finite speed c of a wave and accounts for the reflection conditions at the boundaries. By definition of t , the CFL condition t Ä x=c is always satisfied, so the method is always stable [26,28]. Additionally, because the global error is proportional to the discretization steps, the WFEM is first-order accurate both in space and time [34].

A.7 Matrix formulation
Similar to other numerical methods applied on linear systems, the WFEM can be rewritten in a convenient matrix form which facilitates the process of finding nonsmooth modes of vibration. More specifically, the state vector of the system Q .n/ D Q .n/ 1 : : : Q .n/ N > 2 R 2N at time t n satisfies the identity 2 Q .n/ D AQ .n 1/ ; 8n 1 and Q .0/ is the initial state. The matrix A gathers stiffness and inertial terms as well as the type of boundary conditions. It is now derived for the fixed-free BC. In a matrix form, Eq. (36), which governs the evolution of inner cells, reads Another matrix A can be constructed in the same way for the fixed-fixed BC. Finally, the unknown Q .n/ can be directly expressed in terms of the initial conditions Q .0/ from Eq. (40) by where A n is known for each type of BC.

B Treatment of unilateral contact in WFEM
The unilateral contact constraints involved in the formulation are enforced using the concept of floating boundary conditions [34] which can be regarded as a conditional switch between fixed-free and fixed-fixed boundary conditions [13] when a penetration is detected during a time iteration, as illustrated in Fig. 30. In the continuous framework, these two boundary conditions are where u .n/ N C1=2 is the displacement of the contacting interface which can be calculated by numerically integrating strains of the bar u .n/ ;x at time t n . This equality is used at each time step to check whether contact is active during the next iteration: if g .n/ > 0, a free boundary condition is enforced while if g .n/ 6 0, a fixed boundary condition is considered via the change of matrix A.  Based on the theory of reflection of elastic waves from boundaries [17], the state of the ghost cell C N C1 is updated as follows -Active contact [g .n/ 6 0] -Inactive contact [g .n/ > 0] Equation (36) is subsequently used to calculate the evolution of the average values inside the boundary cell C N , as detailed in section A.5. Additionally, the contact force r .n/ is calculated by employing Eq. (30): r .n/ D S .n/ N cv .n/ N : The sign of this quantity is tracked to locate the time of release, when the bar returns to free condition at x D L. Algorithm 1 is modified to include the floating boundary conditions as described in Algorithm 2. Keep free boundary condition at x N C1=2 Compute gap g .nC1/ and stress and velocity at cell C N via Eq. (39) end -end of floating boundary conditionsend Output: stresses, velocities, contact force at instants t 0 ; : : : ; t nT From the matrix formulation of the WFEM in Eq. (40), two matrices A shall then be distinguished: A f for the fixed-free condition (no contact) and A c for the fixed-fixed condition (contact). Both matrices embed the same stiffness and inertial terms of the system of interest; the only unshared information are the boundary conditions. To summarize, the developed WFEM with floating boundary conditions is a numerically conservative and stable scheme able to properly propagate shock waves induced by a switch in the boundary conditions, the latter being governed by complementarity constraints.