Nonsmooth modal analysis: from the discrete to the continuous settings

This chapter addresses the prediction of vibratory resonances in nonsmooth structural systems via Nonsmooth Modal Analysis. Nonsmoothness in the trajectories is induced by unilateral contact conditions in the governing (in)equations. Semi-analytical and numerical state-of-the-art solution methods are detailed. The significance of nonsmooth modal analysis is illustrated in simplified one-dimensional space semi-discrete and continuous frameworks whose theoretical and numerical discrepancies are explained. This contribution establishes clear evidence of correlation between periodically forced and autonomous unilaterally constrained oscillators. It is also shown that strategies using semi-discretization in space are not suitable for nonsmooth modal analysis. The spectrum of vibration exhibits an intricate network of backbone curves with no parallel in nonlinear smooth systems. Terminology and acronyms The purpose of this chapter is to provide a general picture of the state-of-the-art vibratory analysis of nonsmooth systems. This topic lies at the interface between modal analysis of smooth nonlinear systems and nonsmooth contact dynamics dedicated to the time-evolution of nonsmooth systems, undergoing impact or dry friction, for instance. Some elementary concepts are succinctly recalled for the purpose of completeness. Unless otherwise stated, the epithet discrete (as in “discrete systems” or “discrete setting”) designates semidiscretization in space, while continuous refers to everything else. BEM Boundary Element Method FEM Finite Element Method FEP Frequency–Energy Plot FVM Finite Volume Method HBM Harmonic Balance Method IPP Impact Per Period NSM Nonsmooth Mode ODE Ordinary Differential Equation PDE Partial Differential Equation SM Shooting Method SPP Sticking Per Period WFEM Wave Finite Element Method


Introduction to nonlinear modal analysis
Mechanical systems, from those of large scale (buildings) to those of small scale (MEMS switches), commonly undergo forced vibrations. The efficient and accurate characterization of the response of such systems to external periodic loading is essential for ensuring safe design. It also has various other applications, such as retrofit, damage detection and model reduction, to name a few. In this context, frequency-response curves play a key role for the dynamicist: they indicate the energy level of a periodic solution produced by an external periodic forcing of (angular) frequency !, as a function of !. For nonlinear systems, computing these frequency-response curves is not a straightforward task. Actually, they are known to depend, in a possibly intricate manner, on the forcing amplitude, the forcing frequency, and the forcing shape [50]. A brute-force time-domain approach consisting in solving the governing equations for various external forces and initial conditions is, in practice, not conceivable for large-scale systems. Instead, modal analysis provides a means of computing, for a much more reasonable cost, the so-called backbone curves that shape the forced response curves. Such backbone curves correspond to the underlying autonomous (i.e., unforced) and conservative (i.e., undamped) periodic solutions of the governing differential equation. 1 Autonomous periodic solutions of conservative systems may seem "unrealistic" in the sense that no undamped systems are observable in the physical world. Their investigation can yet provide germane information on periodically-forced and slightly damped systems. Essentially, they extend the concept of spectrum, defined for linear systems, to the nonlinear framework. In particular, they show the energy-dependence of vibration frequencies. The above statements are illustrated by considering a finite-dimensional system relevant to structural dynamics and governed by a linear Ordinary Differential Equation (ODE) of the form where u is the vector of generalized displacements, M is its positive-definite mass matrix, C its damping matrix, K its positive-definite stiffness matrix and f ext its vector of external loadings. The backbone curves are trivially obtained by considering the autonomous conservative counterpart M R u C Ku D 0 and its periodic solutions, yielding vertical lines in the energy-frequency diagram at the eigenfrequencies of the system defined as the square roots of the eigenvalues of M 1 K: in linear dynamics, the frequency of vibration is independent of the magnitude of vibration. This can be seen in Figure 1 where smooth refers to the smoothness of f with respect to u and P u. Its dynamics is unsurprisingly more subtle than the previous linear case, and systematic solution methods for characterizing the vibrations globally are not available [41, 93, 1.3]. However, it is known that fixed points, periodic and quasi-periodic limit cycles may exist, in the vicinity of which it may be possible to approximate the nonlinear dynamical response. In particular, the centre manifold theorem [41,49]; together with Lyapunov's centre theorem [11, p. 5], show that under sufficient regularity 2 and no internal resonance conditions, two-dimensional invariant manifolds exist locally in the phase space and are tangent to the linear modes of the system linearized at the fixed points. Such two-dimensional invariant manifolds were later defined as nonlinear normal modes of vibration [96,114] in the vibration community. They can be understood as curved extensions of linear modes that correspond to flat two-dimensional invariant manifolds defined by one-parameter continuous families of elliptic trajectories. However, the nonlinear framework encompasses many phenomena that are not observed in linear systems, such as internal resonance, frequency-energy dependence, emergence of subharmonics or chaos and existence of isolated loops in the FEP [51]. Again, the relevance of nonlinear modal analysis is illustrated in Figure 1 (bottom), depicting the forced response of a two-dof Duffing oscillator. The response curves warp around the backbone curves. As opposed to the linear spectrum, the nonlinear backbone curves are frequency-dependent and stiffening is exhibited here. The kink of the forced response in the neighborhood of ! 2 =3 corresponds to a subharmonic resonance. The loop near ! 1 corresponds to an internal resonance, where the first nonlinear mode and the third subharmonic of the second nonlinear mode interact. Among all nonlinearities found in mechanics, unilateral and frictional contact nonlinearities form a specific class in which nonsmoothness arises in the dynamics [97]. Typically, the impact between two bodies induces velocity discontinuities and acceleration impulses [2]. The present chapter focuses on the frictionless framework. The governing equation can no longer be written in the form (2), where f is a smooth function of u and P u. However, classical analytical techniques available for computing nonlinear modes [51] require smoothness of the governing equation. Indeed, the invariant manifold approach is based on the Taylor series of the solution written as a function of a pair of master coordinates [96]; the method of multiple scales [77], as a subclass of perturbation methods, requires asymptotic expansions; normal forms rely on the nonlinearity being an analytic function [46]. When it comes to nonsmoothness, such strategies no longer apply.
Nonsmooth modal analysis is the extension of nonlinear modal analysis to nonsmooth systems. This is accomplished by computing nonsmooth modes, that is, families of nonsmooth periodic solutions of the autonomous and conservative dynamics. Even simplistic nonsmooth oscillators exhibit intricate responses [24,107,110]. The regularizing approach, consisting in replacing nonsmoothness with smooth strong nonlinearities [8,18,45,60,73,75,95,117], has the adverse effect of introducing issues such as numerical stiffness [80,81,104] and is not further discussed in this work. Another approach is to include nonsmoothness as such. Many investigations on the dynamics of forced vibro-impact oscillators [36,86,88,120,122] and grazing bifurcations [24,32,78,84] or stability issues [62] are available. The specific targeting of families of periodic solutions of a conservative nonsmooth problem has emerged recently for space-discretized systems [61,108,111] or continuous ones [42,126].
Multiple applications that could benefit from nonsmooth modal analysis can be listed: rotor-stator contact interactions in rotating machinery involving unilateral contact occurrences between blades and casings [39], boiler tube dynamics with a loose support [73,79], grid-to-rod fretting [43], percussive drilling systems [82,83], cutting tools [123] or, on a smaller scale, capsule systems (capsubots) [65,66], and electrostatically-driven and piezoelectric actuators [33,71]. The sensitivity of an atomic force microscope, in tapping mode, can be improved through understanding of the response of impact oscillators [116]. Additional examples include impact dampers implemented to reduce vibrations [64,94], or fret-string contact interactions within musical instruments [14,19,44]. More applications can be found in [8]. All applications have in common the need to properly characterize nonsmooth vibratory resonances.
The purpose of this chapter is to give a picture of the state-of-the-art nonsmooth modal analysis. While the standard procedure in mechanical engineering is to approximate continuous systems by n-dof systems, complications arise when contact is involved. Nonsmooth modes of a continuous system have intricate relationships with that of their semi-discretized counterparts, which raises open-ended questions. The available analytical and numerical methods for nonsmooth modal analysis are first presented for finite-dimensional systems (Section 2) and continuous systems (Section 3). The relationships between forced response and Nonsmooth Modes (NSMs) are then illustrated in Section 4. The comparison between modal analysis of continuous systems and semi-discretized counterparts is addressed in Section 5, which concludes the chapter.

Nonsmooth modal analysis of discrete oscillators
Consider the dynamics governed by a differential equation of the form (2), where f ext D 0. A contact condition, which prevents penetration between two colliding bodies, is commonly expressed as a unilateral constraint g.u; t / 0, where g stands for gap, that is, the distance between the bodies. This constraint is incorporated into the dynamics via a Lagrange multiplier corresponding to the reaction force in the outward normal direction of the contact surface. The non-sticking condition implies that 0 and can be non-zero only if the gap is closed: g.u; t / .t / D 0. These three conditions, known as the Signorini conditions [2], are commonly written in the synthetic form 0 Ä ? g.u; t / 0. In the case of multiple unilateral constraints, each gap function and its corresponding Lagrange multiplier can be stacked in vectors g and , respectively; the inequalities and the orthogonality operator ? are then defined component-wise. Altogether, the autonomous dynamics now writes and nonsmooth modal analysis consists in finding continuous families of periodic solutions to this problem. Equation (3a) should be read in a weak sense, since u is only of regularity C 0 because of the complementarity condition (3b). Various other formalisms are available to describe the dynamics [2].

Necessity of an impact law
An aspect that does not always seem to be understood is that Problem (3), together with some initial conditions u.0/ and P u.0/, does not uniquely determine a solution. For instance, consider a punctual ball of mass m located above a rigid ground and subjected to gravity. When dropped from a given height, the ball first undergoes a free flight uniquely determined by its initial position and initial velocity, together with an ODE of the form m R u C mg D 0 (Cauchy problem). It then reaches the ground: from there, infinitely many solutions are possible, all satisfying Eq. (3) adapted to the problem at hand. The ball could remain on the ground: P u D 0 and D mg. It could also bounce with the same kinetic energy: P u C D P u and D 2m P u at the impact time, where P u (respectively P u C ) denotes the normal pre-impact (respectively post-impact) contact velocity. These two acceptable solutions are depicted in Fig. 2. This non-uniqueness indicates that information is missing. 3 To ensure well-posedness, Eq. complemented with a constitutive impact law. If the latter does not lead to an increase of kinetic energy, uniqueness is guaranteed as soon as the unilateral constraints, the smooth nonlinear terms and the smooth external forces are analytic functions [9,90]. Nevertheless, even with impact laws, the continuity of the solutions with respect to the initial conditions is not guaranteed in the case of multiple unilateral constraints [9].
The necessity of an impact law holds for any unilateral constraint arising in systems semi-discretized in space, unless special treatment is enforced [53]. Numerical strategies which do not explicitly include an impact law, such as [16,124], produce only one among infinitely many possible solutions.
Among possible impact laws, only conservative ones should be considered in the framework of nonsmooth modal analysis, since autonomous periodic solutions are sought. The most common choice 4 is Newton's purely elastic impact law, P u C D P u at impact times. This choice, dictated by the periodicity condition, is incompatible with lasting contact phases observed during collisions in the continuous framework. This can be illustrated by considering the position of the contacting end of a one-dimensional bar colliding with a rigid obstacle, as depicted in Fig. 3. In the continuous framework, contact phases last a finite amount of time, while the energy is preserved (left plot). When the bar is discretized, the conservative impact law implies instantaneous bounces (right plot). For  In the continuous framework, no impact law is needed for well-posedness and the contact is lasting, even for energy-preserving motions. The discretized bar with a conservative impact law exhibits chattering instead.
n-dof systems, lasting contact phases necessitate a purely inelastic impact law of the form P u C D 0, leading to a loss of kinetic energy incompatible with the conservative framework of modal analysis. Also, it is worth mentioning that when subjected to an external load, a unilaterally-constrained system can exhibit lasting contact phases after a countable infinity of impacts occurring in finite time, for non-purely elastic impact laws [13,68]. This phenomenon is called chattering and is illustrated in Sec. 5.
For very specific initial conditions, systems governed by (3), together with a purely elastic impact law, may have solutions with lasting contact phases, also called sticking phases despite the non-sticking condition on the contact force. Such trajectories can be seen as one specific type of contact, as impact or grazing (see Fig. 4). They were investigated in [58] for a linear two-dof spring-mass system. An extension to n degrees of freedom, general mass matrices and a single unilateral constraint is proposed in [109]. In both cases, T -periodic trajectories with one lasting contact phase were shown to exist only for isolated values of T . While they may seem of purely theoretical interest, it was recently demonstrated that such trajectories play an important role in the response spectrum of piecewise-linear impact oscillators [110,Fig. 4]. No systematic results are presently available in the literature on periodic motions with lasting contact phases of systems with additional smooth nonlinearities or multiple unilateral constraints.

Quasi-analytical techniques in simple cases
The systematic analytic derivation of NSM for n-dof systems has recently been provided for a piecewise-linear spring-mass system with one Impact Per Period (IPP) [61], as well as for any piecewise-linear system with a single linear unilateral constraint and an arbitrary number of IPPs [108]. Preliminary investigations show that there exist strong relationships between the forced response of piecewise-linear impact oscillators and backbone curves obtained using NSM, as detailed in Sec. 4. Additional weak smooth nonlinearities do not seem to change the overall picture, as succinctly discussed in Sec. 4.1. Extension to multiple unilateral constraints quickly becomes tedious, because of the combinatorial nature of the sequence of unilateral constraint activations. We now derive the main ideas on how to carry out nonsmooth modal analysis on n-dof piecewise-linear impact oscillators [108]. The generalized displacements and velocities are denoted by u and P u; the state x is such that x > D OEu; P u > 2 R 2n . The unilateral constraint is assumed to be a linear function of the u. As a consequence, there exists a vector w 2 R n and a constant g 0 2 R such that g.u/ D w > u C g 0 . The elastic impact law can be written as [108, 4.2] g.u/ D 0 H) x C D Nx ; (4) where N is similar to a reflection matrix with respect to a hyperplane of R 2n (also known as a Householder matrix), which depends only on w and the mass matrix M. This describes the impact as a simple relationship in terms of the system state x. In the same spirit, let S. /x denote the state after a free flight of duration from a state x. A k IPP motion (k 2 N ) is the succession of one free flight of duration 1 > 0, one impact, one free flight of duration 2 > 0, one impact, and so on, k times. Such a motion is depicted in Fig. 5. Starting from a post-impact state, the where T D 1 C C k . This condition comes with the k gap closure conditions at impact times, that is, g.x.0// D 0; g.x. 1 // D 0; g.x. 1 C 2 // D 0; : : : ; g.x. 1 C C k 1 // D 0: The initial conditions x 0 , determining a motion x that satisfies conditions (5) and (6) for some s D . 1 ; : : : ; k /, define an autonomous periodic motion, provided the gap remains non-negative, in line with (3b). Finding such x 0 reduces to determining a vector 2 R k that satisfies [108] ….s/ D 0 and †.s/ D g 0 j; with j D OE1; : : : ; 1 > 2 R k and where … and † are two k k matrices, whose expressions are known explicitly [108, 3.1] and depend on the parameters M, K and w. The physical interpretation of vector is that it is proportional to the pre-impact contact velocities. Several major consequences follow from (7): it suffices to find the k components of instead of the 2n unknown components of x 0 ; † is invertible almost everywhere in R k , so can be eliminated by combining (5) and (6). As a result, all the periodic solutions are governed by the equation ….s/ †.s/ 1 j D 0. The first step is to solve for s. Then, the corresponding initial state is recovered via x 0 .s/ D '.g 0 †.s/ 1 j/, where ' is a known function (see [108], not recalled here for conciseness); the skew-symmetry of … is such that ….s/ †.s/ 1 j D 0 generically leads to k 1 independent equations. As a result, the set of solutions is a curve in R k and periodic orbits with k IPPs belong to a one-parameter continuous family, corresponding to a two-dimensional manifold in the phase space (see Fig. 6). This feature is shared by smooth nonlinear systems away from internal resonances.   [108]). This NSM is a continuum of periodic nonsmooth trajectories with 1 IPP continuously connected to a linear grazing mode (green ellipse). This two-dimensional manifold is invariant: if a motion starts on it, it will remain on it as time unfolds. In particular, this manifold cannot be intersected by other trajectories in the phase space.
The above methodology is summarized in Fig. 7. Each NSM corresponds to a backbone curve in terms of FEP. An example of such FEP is provided in Fig. 8 for a two-dof spring-mass system (see Fig. 14), with up to seven IPPs. For one to three IPPs, the spectrum was computed with the quasi-analytical method described above. A multiple shooting method was used for four to seven IPPs (see Sec. 2.3.2). Figure 8 displays no isolated branches. Indeed, all backbone curves can either [110]: diverge to unbounded energy, which corresponds to a singularity of †.s/; be connected to a linear grazing mode (this is true in the case for 1 IPP); be connected to another backbone curve, with the junction then corresponding to a nonsmooth trajectory with impacts and grazing; Fig. 7: Summary of analytical nonsmooth modal analysis in the generic case for g 0 ¤ 0. The dependency of x 0 to s is highlighted by the notation x 0 .s/. converge to a motion with one Sticking Per Period (SPP).
In the neighborhood of a 1 SPP, backbone curves seem to converge to the SPP as the number of IPPs increases. This phenomenon is illustrated in Fig. 9. Convergence to trajectories combining 1 IPP and 1 SPP have also been  observed. While very likely to be true, there is no formal proof of such convergence.
As already reported [108,111], seemingly independent backbone curves might be connected through a vertical backbone curve: this additional non-generic feature was referred to as a bridge. This occurs for isolated s, making †.s/ singular. However, such s and those leading to unbounded energy are distinct.
Stability analysis of k IPP motions is carried out in a straightforward fashion by linearizing the kth return map on the hyperplane g.u/ D 0. A perturbation of an initial condition x 0 propagates through the mapping where ı i is an unknown yet small change of duration of the i th free flight. The first-order Taylor expansion of this assumed smooth mapping yields an equation of the form NS. k / : : : NS 0 . `/ N : The unknowns ı 1 ; : : : ; ı k are found by solving the linearized system g.u.
Ultimately, there exists a linear mapping between ıx 0 and ıx through a matrix A.x 0 / such that The eigenvalues of A.x 0 / determine the spectral stability of the periodic solutions emanating from x 0 [93, Summary 7.5].

Numerical techniques
The above (semi-)analytical developments provide essential insight in understanding nonsmooth modes. They are inevitable for proving mathematical results, but are limited to piecewise-linear systems. Numerical techniques take over for more challenging vibro-impact systems, for instance, with multiple unilateral constraints or polynomial nonlinearities.
In the following, we restrict ourselves to two well-known procedures devoted to periodic solutions: Harmonic Balance Method (HBM) and Shooting Method (SM). HBM enforces periodicity exactly by construction, while contact conditions are only approximated. In contrast, SM handles contact conditions accurately, to the detriment of periodicity. Other methods, such as multiple scales, invariant manifold approach and alike, are not considered, as they essentially apply to smooth nonlinearities.

Harmonic Balance Method and its variants
For n-dof systems, setting the unilateral constraints apart, smooth dynamics is described by ODEs in the form where f is a nonlinear function of the displacements u and velocities P u. The unknown displacement u is approximated by u h , which is defined as a linear combination of N chosen shape functions stacked in a vector ' so that where A is a n N matrix of unknown coefficients. Equation (11) is approximately solved by making the residual f.A'; A P '; A R '; t / orthogonal to a well-chosen set of M test functions for the usual inner product Such integrals collectively form a system of nonlinear equations and can be evaluated numerically if the integrand does not easily simplify. Choosing M D N , the nN coefficients in A are then found using a root-finding algorithm such as Newton-Raphson to solve the nN equations (13). This method can be used to compute periodic responses to either forced or autonomous ODEs. Equation (12) shows that the periodicity condition is transferred to a condition on ', which must therefore be periodic. In the case of a periodic external force of angular frequency , periodic solutions are expected to have a frequency multiple of , so T D 2 = can be chosen, or T D 2p = , p 2 N to accommodate possible subharmonics [15,20,55]. In the autonomous case, T is unknown and continuation procedures must be used [3].
The HBM is a well-established technique [56,106] for finding approximations of periodic solutions to (13). It is obtained by specifying ' D D 1 e i!t : : : e iN !t > ; (14) with ! D 2 =T . While commonly producing accurate results for weak nonlinearities, HBM is mostly used heuristically, and there is no proof that a truncated series is a valid approximation of the exact solution [34]. Other shape and test functions ' and shall be adopted. Another well-known method is the collocation method, which corresponds to a low-order piecewise periodized polynomial for ' and 8t 2 OE0; T ; .t / D ı.t t 1 / : : : where t 1 ; : : : ; t N are the collocations points. The Dirac deltas have the property to transform the computation of the inner product (13) into the simple evaluation of u h at the collocation points. The derivatives of u h are computed from the shape functions if they are differentiable, through a finite difference scheme, for instance, or via a conservative Simo scheme [6,101]. When orthogonal polynomials are chosen as shape functions and the collocation points are the roots of one of the orthogonal polynomials, the method is called orthogonal collocation or pseudospectral [12,35] and is reported to be efficient for dealing with sharp fronts [25].
For unilateral contact problems, HBM has mostly been implemented in conjunction with regularizing techniques [15,20,39,55,73] and the contact forces are directly included in the governing ODE (11). A variant of HBM in which the truncated Fourier series is replaced by wavelets has been proposed to compute periodic solutions of a turbine blade with regularized contact conditions [48]. HBM with regularized friction has been investigated in [47].
The unilateral contact conditions can also be treated without regularization, and the problem reads as Above, no impact law is specified. It is instead replaced by the periodicity conditions. We are not aware of any formal proof of this supposed equivalence. However, within HBM, the impact law with e D 1 is implied by the conservation of the total energy in an autonomous problem with no simultaneous impacts, but it is unclear which solutions are picked by the numerical procedure in other cases, such as in the presence of external forces.
The Signorini conditions are transformed by means of a max operator, observing that, for any˛> 0 [91], in the component-wise sense, and can be readily included in (11) at the cost of reducing the regularity of f [99]. The inner product (13) is computed numerically and the solution is found through a semi-smooth Newton solver. An alternative is to implement HBM together with an augmented Lagrangian in a case of unilateral conditions only [57] or a variation of the augmented Lagrangian in the case of friction [74]. Another possibility is to approximate u and with adapted periodic shape functions and satisfy the Signorini conditions at discrete times (collocation points), leading to a Linear or Nonlinear Complementary Problem [70].
Another possible strategy could consist in adding a chosen nonsmooth function with the same regularity as the expected solution (C 0 in the case of impacts) as a shape function; a faster convergence would then be expected, as in the dry frictional case [54]. Irrespective of the chosen discretization, contact-induced nonlinearities require a large number of harmonics (see Fig. 20).

Shooting Method
The Shooting Method is a well-known procedure capable of tracking periodic solutions of ODEs in the form (11) [7,72]. It consists in finding initial conditions .u 0 ; P u 0 / such that they are recovered after a time integration over some interval OE0; T for some T > 0. The analytical method presented in Sec. 2.2 can be understood for one IPP as a SM in which exact integration is performed through matrix exponentials. In more general cases, time integration can be carried out either by event-driven schemes or time-stepping methods [2]. Enforcing periodicity conditions reduces to finding the roots of a vector function z.u 0 ; P u 0 ; T /, bearing in mind that z might be nonsmooth. In modal analysis, T is unknown and there are a priori 2n C 1 unknowns for 2n independent equations: the solution space is a curve, which can be found, as in the HBM, via numerical continuation. The multiple shooting method enlarges the domain of attraction of the root-finding algorithm by splitting the integration domain, increasing the robustness of the numerical procedure [93,105].
This approach was applied to contact problems with regularized nonsmoothness [87,115]. It was also used to locate grazing [26]. The merits of SM for nonsmooth modal analysis rely on the fact that efficient numerical schemes dedicated to nonsmoothness, such as the Moreau-Jean scheme [2,92], can be employed. Convergence proofs exist for a few schemes [29].
For solutions with multiple impacts per period, period T can be replaced by the succession of unknown free flight durations 1 , : : : , k (see 2.2). Complementing the set of equations with k 1 additional conditions of gap closure (g. 1 / D 0, : : : , g. 1 C C k / D 0) imposes prescribed times of impact, which has the advantage of eliminating the nonsmoothness without regularizing the contact conditions. Again, continuation can be used to recover a backbone curve with a given number of IPPs. The robust features of Manlab could be explored in this context [21,113]. One drawback of the shooting method is that it hardly captures unstable parts of backbone curves, because of the time integration [25].
HBM and SM have been combined in the context of forced nonsmooth dynamics in a hybrid method [92]. The linear part of the dynamics is captured by HBM, while SM deals with the nonlinearities.

Gauss' principle
Another possibility would be to consider Gauss' principle for translating the problem of finding a periodic solution into an optimization problem [112]. This principle is known to be equivalent to d'Alembert's or Jourdain's in the nonsmooth dynamics framework [37]. The acceleration field R u solution to an ODE of the form (3) obeys Gauss' principle with the unilateral constraints The idea is to seek periodic solutions by replacing u with a truncated Fourier series u h , as in Eq. (12), and to express Gauss' principle in a weak sense in which the cost function is This yields a problem of the form: find A solution to where S is a chosen set of discrete times in the interval OE0; T . This approach has been adopted for a one-dof system in [4].

Nonsmooth modal analysis of continuous systems
Contact between two linear elastic media generates shock waves featuring discontinuous stress and velocity fronts. For example, when a bar hits the rigid ground, a shock wave emanates at the contact interface, propagates to the free surface of the bar and reflects. The bar departs from the ground when the reflection of the shock wave comes back to the contact interface. Mathematically, the dynamics is described by a Partial Differential Equation (PDE), a solution of which is completely determined by the initial displacement and velocity fields, even in the presence of unilateral constraints [59]: in contrast to the semi-discretized framework (see 2.1), no impact law is needed for well-posedness. The situation is already quite sophisticated for three-dimensional isotropic homogeneous linear elastic materials, where uncoupled longitudinal and transverse waves propagate at distinct velocities. When a nonlinear constitutive law is considered instead, the governing equations are still hyperbolic, but the longitudinal and transverse waves are coupled [30,Chap. 4]. Here, we focus on one-dimensional homogeneous linear elastodynamics and explore solution methods that do not rely on space semi-discretization techniques exposed in Sec. 2.

One-dimensional problem of interest
The system of interest is a fixed-free bar, whose free end is subjected to a unilateral constraint, as illustrated in Fig. 10. The displacement u is assumed to be small compared to the length L of the bar. The dynamics is governed by where g 0 denotes the gap at rest and c D p E= is the wave propagation speed, defined from the Young modulus E and the density of the material. It is worth mentioning that the eigenfrequencies of the linear fixed-free bar are all multiples of the first one, that is, ! k D k! 1 , k 2 N : any initial condition generates a periodic motion and all linear frequencies satisfy an internal resonance condition.

Analytical solution
A few analytical solutions of (20) are available for colliding bars [38] or vibrating strings with an obstacle [14,42] which share similar governing equations. New ingredients are introduced below.
The general solution to (20a) is of the form u.x; t / D f .ct C x/ C h.ct x/, for x 2 OE0; L and t 2 R. In the weak sense, it suffices to require continuity and piecewise C 1 -regularity for f and h. Condition (20b) implies that f D h. Let ' denote the derivative of f . It follows that Condition (20c) implies that @ x u.L; t / D 0 when the gap is open, in other words, '.x C L/ D '.x L/, which means that ' is a 2L-antiperiodic function on R. When the gap closes, it remains closed as long as @ x u Ä 0. In particular, @ t u.L; t / D 0, which is equivalent to '.ct C L/ '.ct L/ D 0, or ' is 2L-periodic. Consider a free phase over OE0; t 1 . On this interval, the displacement field is associated with a 2L-antiperiodic function '. Assume the gap is then closed over OEt 1 ; t 1 C t 2 . The displacement field is then associated with a 2L-periodic function ' 1 .
Introducing the function , defined over R by 2L-antiperiodicity and the value 1 over OE L; L/, it can be shown that the periodicity condition reduces to the following condition on ' 5 : The problem of finding (potential) periodic solutions with one contact phase per period for the unilaterally constrained bar hence reduces to finding ' solutions of (22). The period is given by T WD t 1 C t 2 . Three additional conditions apply, which can be understood as admissibility conditions [108,109]: the contacting end of the bar must not penetrate the obstacle during the free flight: at x D L, the bar must remain in compression during the contact phase: the gap must be closed at t 1 : Eqs. (22) and (23) can either be solved collectively to find periodic solutions or be used to check the correctness of a candidate periodic solution identified from numerical methods. An interesting direct consequence follows from the absolute value of (22): 8x 2 R, j'.x/j D j'.x C cT /j. Recall that ' is 2L-antiperiodic, so j'j is 2L-periodic, and also cT -periodic. This is possible only if cT =L is a rational, or if j'j is constant. Continua parametrized by T are hence possible only if j'j is constant, meaning that all backbone curves, which are not vertical lines, correspond to piecewise-linear displacement fields, that is, piecewise-constant velocity fields.
Not only does it provide a sound mathematical basis, this approach was proven successful for rediscovering the nonsmooth modes previously conjectured [126], see Fig. 11. The main backbone curves emanate from the linear eigenfrequencies of the fixed-free bar. The additional curves correspond to subharmonics of higher frequency modes. The functions ' labeled a , b and c in Fig. 11 are plotted in Fig. 12.
Among the solutions to (22) are the two main NSMs determined by ' 1 and ' 2 . Each of these functions is defined by its value over OE L; L and its 2L-antiperiodicity over R. For the first one, where the duration of the contact phase t 2 relates to T through T D 4L=c t 2 . For the second mode, where t 2 satisfies T D .4L=c t 2 /=3. In both cases, the mode is parametrized by t 2 , or equivalently, T or !. The coefficient˛is such that u.L; 0/ D g 0 and is not explained for the sake of conciseness. The displacement 0 a 1st NSM field u, calculated from ' using Eq. (21), is depicted for the first two linear and nonsmooth modes in Fig. 13 with appropriate labels. The first nonsmooth mode and its linear counterpart show similar features: this also holds for the second mode, where both exhibit nodes of vibration. However, standing waves in the linear setting become travelling waves in the unilateral setting, where the characteristic lines are clearly identified. The analytical approach developed above is limited to simple systems such as the one considered. Numerical techniques capable of handle more general systems are now exposed.

Finite volumes and the Wave Finite Element Method
Finite Volume Methods (FVMs) form a family of numerical methods widely used in fluid mechanics [63] to solve PDEs. By construction, they are designed to enforce conservation laws. They consist in discretizing the space domain into cells. As opposed to other well-known numerical techniques such as the Finite Element Method (FEM), the strong form of the PDE is considered and the unknown field is averaged in every cell through volume integrals. Time evolutions are calculated via fluxes on the cell boundaries. In the one-dimensional case, the wave equation (20a) is recast into a system of two hyperbolic conservation laws: where v and are the velocity and stress fields, respectively. The Wave Finite Element Method (WFEM) is a shock-capturing FVM, where the time-discretization is coupled to space in such a way that waves propagate along the characteristics lines of the hyperbolic PDE [100]. The Dirichlet-type fixed boundary at x D 0 can be dealt with straightforwardly using ghost cells [63]. The treatment of the unilateral contact condition is more challenging: one possibility is to use the floating boundary condition technique [100], which can be understood as a conditional switch between free and fixed boundary conditions. Finding periodic solutions of the colliding bar reduces to finding the initial stress and velocity fields, in the form of constant averaged values in every cell, which propagate along the characteristic lines, satisfy the clamped boundary condition at x D 0 and the switches between fixed and free boundary at x D L such that the initial state is recovered at time T after a prescribed number k of contact phases per period. The analytical backbone curves in Fig. 11 are retrieved with this approach [126]. A more complicated configuration in which, at x D 0, the Dirichlet condition is replaced with a Robin condition of the form @ x u.0; t / D˛u.0; t / is also of interest, since the internal resonance condition previously mentioned no longer holds. No analytical results could be derived, but nonsmooth modes can be numerically computed.
Also, WFEM implies a projection step when penetration is predicted. This should not be confused with an impact law, since the exact solution of a bouncing bar [27] and the exact solutions in Fig. 13 are retrieved. 6 The main drawbacks of semi-discretization in space are not observed: in particular, there is no chattering, the velocity of the contacting end undergoes a jump at gap openings, and the energy is accurately preserved. Forced responses can be computed as well [126]. However, extension to higher dimensions looks challenging. Indeed, the description of how a discontinuity (between two finite volumes) propagates, the so-called Riemann problem, can no longer be solved exactly. Moreover, conservation laws in the multidimensional framework raise a number of issues that are not well-understood [67].

Boundary Element Method
Problem (20) where H is the Heaviside function. Using u as the weighting function in the space-time integral form of Eq. (20a) yields where the last integral stands in the distributional sense. This is the principle of the TD-BEM in one dimension. The general solution is a linear combination of u defined in (29), which is a progressive wave. The Convolution Quadrature-BEM (CQ-BEM) [1,89] computes the integrals in (31) via the Convolution Quadrature Method. They can also be computed with piecewise-constant or piecewise-linear polynomials [17]. After discretizing space and time integrals, the sought solution u becomes a linear combination of the boundary conditions u.0; /, @ x u.0; /, u.L; /, @ x u.L; / and the initial conditions u 0 and v 0 . Due to clamping at x D 0, u.0; / is known and @ x u.0; / is unknown. The contact condition at x D L corresponds to either a free or a fixed boundary condition, and the switch is triggered by monitoring the gap and the normal contact force. In either case, exactly half of the boundary conditions are known and half are unknown. The latter can be deduced from the evaluation of (31) at D 0 and D L, providing two equations in two unknowns at each prescribed time step. The displacement of internal prescribed nodes can then be recovered through (31).
When targeting periodic solutions, the shooting method (see 2.3.2) can be used, together with the discretized governing equations obtained from (31), providing 2n equations where n is the number of discretized space nodes, for 2n C 1 unknowns (the initial conditions at the n nodes plus the period T ). Again, numerical continuation techniques involving a semi-smooth Newton solver are employed to find nonsmooth modes of vibration.
This approach was successful in computing the first two main backbone curves, some subharmonics and internal resonance backbone curves, see Fig. 11. The main challenge for the extension to the multi-dimensional framework is that the fundamental solutions are only known exactly for simple geometries. In such cases, Green's functions (which are fundamental solutions with specified boundary conditions) can, however, be approximated numerically [28,Chap. 7].

Space-time finite differences
Many other discretization schemes relying on finite differences are available for hyperbolic PDEs [69]. We focus on numerical methods that simultaneously discretize space and time. When discontinuous solutions are expected, common methods include Lax-Friederich, Lax-Wendroff, MacCormack Upwind, Forward-Time-Centered-Space (FTCS), and Leapfrog. These schemes all stem from truncated Taylor series, and differ mostly according to their order in space and in time. For example, the FTCS method is second-order in space and first-order in time. Another method can be derived as follows. Writing the second-order Taylor series of u in time then replacing @ t u with c@ x u (and thus @ 2 t t u with c 2 @ 2 xx u) [85], and applying a first-order central difference for @ x u and second-order central difference for @ 2 xx u yields which is the well-known Lax-Wendroff scheme. Stability is governed by the Courant-Friederich-Lewy (CFL) condition, which provides a necessary condition (sometimes sufficient) on the time step t , given the wave celerity c i in direction i and the space discretization step x i , taking the following form in N dimensions: where C CFL depends on the finite difference scheme. The interpretation of this condition is that numerical "information" should not propagate slower than physical "information". It is not a sufficient condition for stability. The main issue with these finite difference schemes for propagating discontinuous fields is that they are either first-order accurate, thus numerical viscosity 7 smoothens the solution, or second-order accurate, in which case they are dispersive, 8 leading to numerical oscillations known as Gibbs phenomenon. Apart from Glimm's method, which suffers from inaccuracy during smooth phases [22], this is clearly illustrated by several examples in [103].
A common strategy for reducing spurious oscillations is to add numerical diffusion tuned to the Gibbs phenomenon. This approach is problem-dependent, and may therefore be tedious to accomplish, and is hardly compatible with periodic solutions. Possibly more promising are limiters [10,31]. Signorini conditions and impact laws have yet to be incorporated into this formalism. As of now, it is unclear whether these approaches would be suitable for nonsmooth modal analysis. Other potentially relevant methods are listed in [31].
Mixed space-time HBM [121] or a time-space FEM with a discretization along the characteristics 9 for onedimensional systems might also be useful for nonsmooth modal analysis, and are hence worthy of further investigation. 7 Numerical viscosity, or diffusion, arises when the numerical scheme introduces a velocity term with a positive prefactor. 8 Numerical dispersion occurs when the numerically approximated propagation celerity of a wave depends on its frequency. Note that dispersion and numerical dispersion are two distinct concepts. 9 The unknown displacement field would be expanded as where the i could be the usual hat functions, for instance.

Relationships between forced response and nonsmooth modes
Various analytical and numerical methods capable of performing nonsmooth modal analysis have been reviewed in the preceding sections, both in the discrete and continuous frameworks. Some are sufficiently mature for nonsmooth modal analysis, while others have yet to be thoroughly explored, as their usability has not been comprehensively assessed. In the following, the nonsmooth modal analysis of a FEM-like semi-discretization of the colliding bar is explored by means of the analytical and the multiple shooting methods. In the continuous framework, nonsmooth modal analysis of the same system is carried out with WFEM and analytical techniques. The fact that peaks of resonance in the forced response emerge along the backbone curves in the FEP demonstrates the main purpose of nonsmooth modal analysis.

Discrete oscillators
Recall that all numerical methods detailed in Sec. 2.3 are capable of computing periodic solutions of a forced system. The brute-force approach is another possibility, which does not work in the autonomous case. It consists in time-integrating Eq. (2), where f ext is periodic in time, until a periodic response is obtained or a stopping criterion is reached [76]. This simple technique is CPU-intensive when damping is light and the detection of the steady-state may be delicate. Nevertheless, it was implemented to compute the forced response of the two-dof spring-mass system in Fig. 14 with n D 2. Results are presented as a function of the forcing period in Fig. 15 (top), where colors indicate the number of impacts per period. For clarity, only the five lowest IPPs are shown, even though solutions with as high as 24 IPPs were found, an example of which is depicted in Fig. 16 Fig. 15: FEP of a two-dof impact oscillator. Responses with 6 IPPs or more are excluded for clarity. Colors correspond to IPPs with labels used in Fig. 17.
of the forcing. For instance, the period of the 24 IPP-response in Fig. 16 is 8T 0 , where T 0 is the forcing period. Accordingly, the results in Fig. 15 (top) can also be plotted as a function of the response period (see Fig. 15 (bottom)). This results in a correlation between the number of IPP and the response period: IPP curves are clustered. It also shows that identical nonsmooth resonances can be obtained for distinct forcing periods: the two 1 IPP resonance peaks in the top plot seem to correspond to the same resonance in the bottom plot. The purpose of modal analysis is to predict vibratory resonances. Using the analytical method described in Sec. 2.2, it appears that resonance peaks in the forced response mostly emanate in the vicinity of NSM backbone the branches in Fig. 15 look like they were connected to NSMs. The way in which nonsmooth modes relate to forced responses is not restricted to peaks in the FEP, but rather extends to shapes. This is illustrated in Fig. 18 for a two-dof system in which the forced response trajectories of the masses are compared to the nonsmooth modal shape of the same period. Though no longer symmetric, the forced response is strikingly similar to the periodic solution of the autonomous problem. The above observations extend, in part, to Duffing impact oscillators, as depicted in Fig. 14 for n D 2. The corresponding autonomous dynamics between impacts is governed by where is user-defined. The previous piecewise-linear system corresponds to D 0, but no analytical techniques exist to compute the periodic solutions when ¤ 0, except for n D 1. Using the shooting method between time instants 0 C and T , as described in Sec. 2.3.2, both the backbone curves and the forced response curves can be computed for several values of . They are exposed in Fig. 19, in the neighborhood of a backbone curve with 1 IPP. In this figure, the thick backbone curve is the one in Fig. 17 (left), plotted in terms of frequency. It continuously deforms as the cubic nonlinearity increases. The forced response changes accordingly, so that even in the piecewise-nonlinear case, nonsmooth modal analysis seems to provide backbone curves that perfectly support the forced response curves. We now proceed with the illustration of HBM, as described in Sec. 2.3.1 for a one-dof impact oscillator [98]. Figure 20 shows the approximated backbone curves for an increasing number of harmonics: the backbone curve converges to the exact one. Also, the approximated forced response is seen to be perfectly organized around the backbone curves. The time evolution of position (right plot) shows that the residual penetration gets smaller as N increases. This very simple example establishes numerical evidence that when periodicity is enforced, constitutive impact laws are unnecessary.

Continuous oscillators
This subsection succinctly extends the previous results to the continuous framework by exploring the autonomous and forced dynamics of a one-dimensional bar colliding with a rigid wall (see Fig. 10). As explained previously, backbone curves can be obtained analytically (Sec. 3.2), via WFEM (Sec. 3.3) or using TD-BEM (Sec. 3.4). The first four main backbone curves are depicted in Fig. 21 together with the periodically forced response at various energy levels. The top plot corresponds to an excitation induced by a harmonically moving obstacle, while the bottom plot considers an external periodic and distributed force along the bar. As in the discrete case, the main peaks of the forced response align, in both cases, with the main backbone curves. The additional minor peaks on the top plot might correspond to internal resonances. However, this point requires further work.
Similarities between autonomous and forced responses also emerge in terms of frequency and modal shapes. For instance, Fig. 22 compares one periodic solution belonging to the first nonsmooth mode to a forced response arising in its vicinity. It is remarkable that the forced response is dominated by the resonant response, that is, the first mode shape (see Fig. 22, left), which is only slightly altered by the type of external forcing.

From discrete to continuous NSM: similarities and differences
We have seen that space-continuous and space-discrete models fall under two different paradigms. In the first category, contact is simply a constraint from which emanate shock waves propagating in the continuous solid. The second category introduces a number of pitfalls and difficulties. An impact law is required, propagation of a wave is difficult to approximate accurately, and lasting contact phases are hardly compatible with the conservation of energy required by the periodicity condition. Additionally, the regularity of the generalized positions is higher than in the continuous case, characterized by discontinuous velocity waves and not just the degree-of-freedom involved in the unilateral constraint.
This last section attempts to highlight the similarities and differences between the two "worlds" within the unidimensional framework presented in Section 3.1.

Without unilateral contact constraints
Unilateral contact conditions are temporarily set aside. In structural dynamics, the Finite Element Method is widely used to discretize PDE (20a). Loosely speaking, the weak form of (20a) consists in finding u such that for all v in an appropriate space x/v i for some chosen shape functions 1 ; : : : ; n , approximating u and v by u h and v h in (37), respectively, leads to a system of ODEs standard in structural dynamics: where M and K are calculated from (37). In the sequel, we consider, for simplicity, the space semi-discretization of the clamped-free bar with punctual masses (see Fig. 14). Accordingly, M D mI n and K D k Space-discretization formulations are not able to capture the progressive nature of shock waves properly and may lead to non-causal spurious oscillations in space [40]. In order to explain this, let us compare the modal properties of the continuous bar with those of the spring-mass system. The eigenfrequencies and corresponding modeshapes of the continuous bar are given by [38] In contrast, the eigenfrequencies of the discrete system are When n p, using (40), the result is that the eigenfrequencies of the discrete and the continous bar are equivalent: Relating the node j of the discrete system to the position x in the bar via x D L.j 1/=.n 1/, an analogous consequence holds for the mode shapes Q p and p .x/. This is shown in Fig. 23 where both the eigenfrequencies and the eigenmodes are in good agreement in the low-frequency range. However, when the index p is no longer negligible compared to n, the approximation becomes inaccurate. By injecting a progressive monochromatic wave of the form u.x; t / D e i.!t Äx/ into the wave equation (here, i stands for the imaginary unit), it results that Ä D c=!, which constitutes a linear dispersion relation: the phase and group velocities coincide and there is no dispersion. Now, let x denote the space discretization step such that x D L=n. A progressive monochromatic wave of the form u p .t / D e i.!t p Q Äx/ in (38) propagates by satisfying When x Ä D c=!, then Q Ä Ä, translating the fact that low-frequency waves propagate at the same velocity as in the continuous bar. Nevertheless, dispersion appears for higher frequencies, as illustrated in Fig Fig. 24: Time evolution of a spring-mass system with n D 100. All initial displacements and velocities are zero, except a unit initial velocity for the free node n. The main wave propagates at the velocity c, but spurious oscillations become visible due to dispersion. Index j goes from 1 to n. The black curves correspond to the trajectory of every fifth degree of freedom. Trajectories are merged in a surface to facilitate visualization.
shows the time histories for zero initial displacements and velocities, except a unit initial velocity on the free node n. Even with a relatively high number of degrees of freedom (n D 100), the solution displays spurious oscillations due to the dispersion of high frequency waves. This questions the relevance of the space semi-discretization formalism when shock waves are sought. A comparison between numerous different schemes is proposed in [27]. Even the most accurate of them yields significant discrepancies with the exact solution, even after only one (pseudo-)period of motion [5,23].

With unilateral contact constraints
The relationships between nonsmooth modes and forced response curves have been presented in Sec. 4 for discrete and continuous systems, separately. The relationships between discrete NSMs and continuous NSMs is now examined in an exploratory and qualitative manner.
As discussed in Sec. 2.1, the space semi-discretization of a PDE brings in the necessity of an impact law: modal analysis requires e D 1 for energy conservation, while e D 0 is needed if sticking phases are of interest. Sticking phases are meaningful as they emerge naturally in the continuous framework (see Fig. 3). Some authors have proposed the mass redistribution method. It removes the mass of the contacting node and redistributes it to other nodes [23,52], so that kinetic energy is not affected. However, it is not clear how it differs from a penalization approach. In the same vein, a recent exploratory work that incorporates an elastic law e D 0 proposed to redistribute the kinetic energy of the non-massless contacting node to the neighboring mass [125]. Let us now analyze the sensitivity of the responses to e with n fixed, and to n with e fixed, respectively.
It is observed that the sensitivity of the solution to e reduces when n increases. Figure 25 displays the periodic forced responses for various e and n, obtained using a Moreau-Jean scheme together with a Â-method (Â D 1=2) [2]. For n as small as 20, displacements of the masses are not much affected by e, meaning the forced response curves computed for various e are very similar. Chattering obtained for e > 0 seems to have a negligible effect on the overall dynamics [81].
Interestingly, when scaled with respect to the length L, the local behavior of the contact node for large values of n is indistinguishable from that of the continuous bar. This is illustrated with e D 1 in Fig. 26 where the periodic solution with n 2 f5; 20; 1000g is compared with the continuous periodic solutions produced by WFEM (see Sec. 3.3). In particular, no elastic bounces are visible and the contact behaves very much like the lasting contact experienced in the continuous framework. Indeed, the solutions seem to converge with n, irrespective of e, to the solution of the continuous bar. Overall, this paradigmatic difference between the continuous and the discrete systems with forcing and damping vanishes as n becomes large. The chattering phenomenon appears to be the pivot between the discrete and the continuous frameworks. Damping is likely to play an important role as well, since it acts like a low-pass filter, and thus reduces the discrepancies between continuous and discrete models mentioned in Sec. 5.1.
Naturally, one may wonder, in the autonomous and conservative framework, how the backbone curves of the continuous bar compare to the ones of the semi-discretized bar. More explicitly, we would like to approach the backbone curves in Fig. 11 according to the ones exhibited in its n-dof counterpart, as in Fig. 8, for a sufficiently large n. The challenge comes from the fact that when n becomes large, the spectrum is extremely dense and numerically demanding and currently not accessible. Nonetheless, we provide a few clarifications. In Fig. 27, the energy averaged over six forcing periods for n D 5 and n D 20 is plotted, for two levels of forcing and several levels of damping. For n D 5, the resonance peaks roughly correspond to the main backbone curves of the continuous bar. For n D 20, the agreement is clear, and thus, irrespective of the level of damping, for the first two modes. Fig. 27 also shows that the forced response curve is very jagged for low levels of damping (dark curves) and becomes a smooth function of the forcing period as damping increases. This can be understood by plotting the position as a function of time for distinct damping levels, as illustrated in Fig. 28. The forcing magnitude is tuned to approximately maintain the magnitude for the position u n , to compensate for increasing damping levels. The three following types of forced response curves can be distinguished: For low damping, the forced response curve is governed by k IPP nonsmooth modes, as stated in Sec. 4.1. The backbone curves feature a number of small branches 10 (see Fig. 15 for n D 2). It follows that a forced response is very sensitive to the forcing frequency, as witnessed by the numerous irregularities in the forced response curves in Fig. 27. This situation corresponds to the top left plot in Fig. 28 where a 6T -periodic response with 3 IPPs is observed. For moderate damping, the forced response curve is smoother and the trajectory is simpler. This corresponds to the 2T -periodic response with 4 IPPs (top right) and T -periodic response with 2 IPPs (bottom left) in Fig. 28. For large damping, the response curves are very smooth, as shown in Fig. 27. Contact settles through chattering mechanisms, and the macroscopic coefficient of restitution, i.e., seen from the scale of the whole system, is e D 0, even though the computations were performed with e D 1.
Given that the motion of the discrete system converges to that of the continuous bar for sufficient damping, it is not surprising that for medium or high levels of damping, the resonance peaks are close to those of the bar, for n D 20. More surprising is the fact that nonsmooth resonances for low damping also match the continuous backbone curves for large n. In other words, for low damping, as shown in Sec. 4, the forced response of a n-dof system appears to be driven by its (discrete) NSMs, at least for small n. Figure 27 shows that, for large n, this forced response resonates along the backbone curves corresponding to the NSMs of the continuous system. Accordingly, there must be a relationship between the backbone curves of the discrete system and those of the continuous one. This is presently to be clarified, even in the one-dimensional framework, because computing the FEP for the autonomous case with large n is challenging.
We close this chapter with two observations, which tend to confirm some degree of correlation between backbone curves in the continuous and discrete settings. The first one is concerned with the non-existence of nonsmooth modes for the continuous bar in certain frequency ranges. The continuous bar does not seem to possess any backbone curves within the range !=! 1 2 OE2; 3 in Fig. 11 and the same applies to the discretized bar for large n: nearly no periodic solutions are detected in this range, at least for 1 IPP, which is encouraging.
The second point relates to similarities in grazing motions: In the vicinity of ! 1 , the continuous bar features two grazing modes, as illustrated in Fig. 29: the linear grazing mode of the clamped-free bar, which is a sinusoidal function in time, and the nonsmooth grazing mode, which is a triangular function in time (see Fig. 30 (left)). This triangular shape corresponds to the limit case when the mode shape shown in Fig. 13 (bottom left) has a contact duration approaching 0, and can be found exactly from ' given in (24), that is, '.x/ 2 D 1 for x 2 OE L; L and 2L-antiperiodic, by evaluating integral (21). For the discrete bar, the corresponding linear grazing mode is a sine of frequency Q ! 1 as well. There is a priori no equivalent for the triangular function found for the continuous bar, since the modal manifold is known to be continuous for any fixed n [61].
However, the triangular shape can be recovered in the discrete setting for large n, as a 1 IPP trajectory. Let us focus on the contacting end of the first nonsmooth mode, for a grazing amplitude. From [108, Eq. (93a)], the position of the nth mass with 1 IPP is cos. Q ! j .t T =2// Q ! j sin. Q ! j T =2/ Q 2 j;n ; where Q ! j and Q j;n are given by (42) and (43). The value of is such that u n .0/ D g 0 (closed gap); to simplify, only the time-domain shape of u n .t / is studied, its magnitude being dropped. When T approaches 2 = Q ! 1 , the first In the continuous framework, the first linear grazing (a sine function in time) and first nonsmooth grazing (a triangular function in time) trajectories share the same frequency ! 1 . For the discrete system, the linear grazing mode is also a sine. When n is sufficiently large, the triangular shape of the continuous nonsmooth grazing mode is retrieved. [ ] Displacement of the contacting end u.L; t / and u n .t /, respectively, along the first linear grazing mode. Other positions within the bar/discrete oscillator are also indicated from x 0 to x L (white to black).
term of the sum dominates and the shape converges to cos..t T =2/ =2/. This situation corresponds to the first linear grazing mode. For the triangular shape, it should be observed that the sum is dominated by the first terms such that for some n 0 n and j Ä n 0 , Q ! j .2j 1/=2 and Q j;n sin..2j 1/ =2/ D . 1/ j C1 , u n .t / can be approximated by Now, concerning the period of the first grazing nonsmooth mode T D 2 =! 1 D 4, and since sin.2 Q ! j / .2j 1/ =.2n/ when j n, the shape is also similar to which is the truncated Fourier series of a triangular wave. As n ! 1, and Q ! 1 approaches the neighborhood of ! 1 (or its multiples), both the exact grazing sine and the approximated grazing triangular function are found (see Fig. 30 (right)). The n-dof system thus mimics the continuous bar's nonsmooth grazing behaviour. The corresponding energies, for the discrete and continuous systems are also found to be comparable.
The (nonsmooth grazing) triangular displacement reported above emerges because it can be expressed as a combination of the linear modes of the clamped-free bar, whose time-domain participations in the nonsmooth periodic solution follow a Fourier sequence of fundamental frequency ! 1 : this unique attribute stems from the full internal resonance condition enjoyed by the continuous bar considered and no longer holds when this condition is not satisfied.

Conclusion
In the literature, nonlinear modal analysis is recognized as a matured tool for smooth nonlinear vibratory systems of small to moderate size. However, new methods of analysis are needed when vibro-impact dynamics and unilateral contact conditions are involved. Nonsmooth modal analysis is one such tool. It consists in finding continuous families of periodic solutions of unforced nonsmooth systems, as specified by the definition of modal analysis. Existing solution methods serving that purpose, including very recent developments, were presented in this chapter for simplified systems in the form of a one-dimensional continuous bar and a corresponding n-dof discrete spring-mass oscillator. Conceptual dissimilarities between these two frameworks are summarized as follows: For modal analysis purposes, the discrete setting necessitates an energy-preserving impact law with restitution coefficient e, while the continuous setting does not. The discrete setting with an energy-preserving impact law generates chattering, which manifests itself as k-IPP trajectories that are challenging to capture numerically when k and n grow. Chattering was found to be the pivot between the discrete and continuous worlds. It is not clear whether the backbone curves (which define the nonlinear spectrum of vibration) of the discrete oscillator converge towards the backbone curves of the continuous system as n increases.
The sensitivity to the restitution coefficient e of the periodically forced displacement of the discrete oscillator with low damping decreases with n. The backbone curves calculated in the continuous setting accurately predict the vibratory resonances of the discrete oscillator for a sufficiently large n, irrespective of e. By virtue of the above comment, vibratory resonances of the continuous bar and discrete oscillator are in good agreement as soon as n is sufficiently large. The peaks of resonance are not much affected by the type of forcing (distributed, concentrated at the contacting end, or far from the contact zone).
In the long run, the aim is to settle Nonsmooth Modal Analysis as an attractive and standard engineering tool aiding in the the efficient prediction and comprehension of nonsmooth vibratory signatures, in replacement of tedious time-domain computations. Among the various possible avenues to be explored in the future, the following are pressing issues: In the finite element framework, removing the problematic chattering could be overcome by taking advantage of the vanishing influence of the impact laws for large n and choosing a purely inelastic impact law, that is, e D 0. The very small loss of energy should be compensated for in some way. Nonsmooth Modal Analysis of multi-dimensional systems should be tentatively performed employing Finite Volumes and the Time-domain Boundary Element Method. In the context of continuum mechanics with the assumptions of large displacements and strains, smooth nonlinearities emerge. The resulting dynamics involving unilateral contact constraints should be addressed.