M Attar

A Karrech
email: ali.karrech@uwa.edu.au

K Regenauer-Lieb

Non-linear modal analysis of structural components subjected to unilateral constraints

Keywords: Unilateral constraint, Non-linear modal analysis, Piecewise-smooth system, Non-linear normal modes

In this paper, we present a detailed numerical study of the non-linear dynamics in structural components under unilateral contact constraints. Here, the unilateral term characterises the constitutive law of the restoring force in the constraints as they only sustain elastic reactions in one direction, either compressive or tensile.

Thus, the non-differentiability of the contact law at the discontinuity point is the only source of non-linearity. In our approach, the discrete lattice method (DLM) is used to treat the continuous system as a piecewise linear model. Thus, the trajectory of each node in the discrete model would be a sequence of smooth solutions with the switching times between them. The application of the one-step integration scheme allows us to detect the occurrence of contact (i.e. the instants that the lattice nodes cross the discontinuity boundary) and consequently update the active

Introduction

During the past few decades, significant research efforts have been focused towards the theory and application of analytical and numerical models to predict the dynamic behaviour of mechanical systems. Various substructures of these systems can be presented as thin members such as prismatic plates, slender rods or beams that usually operate in contact with other flexible members.

Conventionally, the structural supporting components and their mechanical characteristics can be modelled by spring and damper elements with a smooth constitutive law for their restoring forces. This constitutive relation is based on the assumption that the interconnected body and deformable supports interact through a perfect connection at their interface. Therefore, the resulting reaction is proportional to the relative displacement at the contact interface with the same stiffness coefficient in compression and tension (i.e. bilateral constraint).

The study of the dynamic properties of such a system has widespread applications, such as predicting the effect of scour on the frequency response of a pile partially embedded in soil [START_REF] Prendergast | An investigation of the changes in the natural frequency of a pile affected by scour[END_REF], dynamic analysis of damaged engineering structures [START_REF] Attar | A transfer matrix method for free vibration analysis and crack identification of stepped beams with multiple edge cracks and different boundary conditions[END_REF][START_REF] Attar | Free vibration analysis of a cracked shear deformable beam on a two-parameter elastic foundation using a lattice spring model[END_REF][START_REF] Attar | Dynamic response of cracked timoshenko beams on elastic foundations under moving har-51 monic loads[END_REF], detection of resonance speed of railway tracks in high-speed transportation systems [START_REF] Tran | Vertical dynamic response of nonuniform motion of high-speed rails[END_REF], stability analysis of fluid-conveying pipes [START_REF] Doaré | Dissipation effect on local and global stability of fluid-conveying pipes[END_REF] and post-buckling behaviours of functionally graded plates resting on elastic foundations [START_REF] Duc | Nonlinear postbuckling of symmetric S-FGM plates resting on elastic foundations using higher order shear deformation plate theory in thermal environments[END_REF]. For these linear dynamic settings, the free and forced vibratory responses are always combinations of the system's linear normal modes (LNMs).

However, in practice, the global dynamic behaviour of a dynamic setting can be strongly influenced by different sources of non-linearity, such as non-linear material laws or large deformations [START_REF] Mullapudi | Nonlinear finite element modeling of beams on two-parameter foundations[END_REF][START_REF] Hendou | Transient analysis of nonlinear eulerbernoulli micro-beam with thermoelastic damping, via nonlinear normal modes[END_REF]. Non-linearities may also arise from the connecting elements in the aforementioned bilateral constraints where they exhibit non-linear mechanical properties to describe the supporting reaction [START_REF] Vakakis | Dynamics of linear discrete systems connected to local, essentially non-linear attachments[END_REF]. As a result, the conventional linear techniques no longer present a valid mathematical model of the real physics and they cannot be used anymore to express the general transient response as a linear superposition of normal mode-shapes. Moreover, any attempt to overlook the non-linearity or to seek a simplified model that linearises the system at the forcing level will yield inaccurate results, or at best, a suboptimal design.

The non-linearities due to the interactions between the solid body and deformable supports can be generally classified into two broad categories: smooth non-linearity where the local contact force is only a non-linear function of the local relative displacement at the interface but with equal magnitudes in compression and tension [START_REF] Vakakis | Dynamics of linear discrete systems connected to local, essentially non-linear attachments[END_REF] (i.e. the force-displacement constitutive contact model is always differentiable with respect to the displacement component), and non-smooth non-linearity where the magnitude of local restoring force is also a function of the local interpenetration direction between two objects [START_REF] Attar | A novel technique for dynamic analysis of beam-like structures on tensionless elastic foundations subjected to moving loads[END_REF] (i.e. there is a discontinuity boundary in the vicinity of the interface where the contact model is different at two sides of this border and this makes the constitutive equation non-differentiable at the break-point).

Dynamic analysis of systems with smooth non-linearities has been an active research domain for many years [START_REF] Vakakis | Dynamics of linear discrete systems connected to local, essentially non-linear attachments[END_REF][START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF][START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF][START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF][START_REF] Vakakis | Non-linear normal modes and their applications in vibration theory: An overview[END_REF][START_REF] Chati | Modal analysis of a cracked beam[END_REF][START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF]. Among the first scholars who studied such problems, Rosenberg [START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF] extended the concept of modal analysis theory to non-linear systems and firstly introduced the definition of non-linear normal modes (NNMs) in the 1960s as an essential and robust tool to interpret dynamic properties of non-linear models. Since then, various definitions of NNMs for conservative and non-conservative models have been proposed over the years to gain more insight into the behaviour of non-linear systems under different sets of operating conditions [START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF][START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF][START_REF] Vakakis | Non-linear normal modes and their applications in vibration theory: An overview[END_REF]. Among them, we use the extended version of Rosenberg's definition in the present study where NNMs are characterised as not necessarily synchronous periodic motions of the undamped free non-linear system in which all physical coordinates vibrate equiperiodically [START_REF] Chati | Modal analysis of a cracked beam[END_REF][START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF]. In contrast to LNMs, although NNMs are not orthogonally independent (i.e. superposition principle is not valid for them), it is clear that they have significant potential to describe salient features of non-linear systems which have no counterparts in the underlying linear model [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF]. For instance, non-linear modes can describe the resonance frequencies of the smooth non-linear system as a function of their conserved energy where the number of NNMs may even exceed the number of degrees of freedom (DOFs) in the space-discretised system. From an engineering standpoint, the importance of resonance frequencies is obvious as the system forced response can lead to large persistent oscillations in the neighbourhood of them and the structure is at the highest risk of failure. Moreover, NNMs may be used to explain other distinct phenomena such as bifurcations, modal interactions (internal resonances), subharmonic resonances and superharmonic resonances that cannot be captured by the linear theory.

Recently, significant progress has been achieved in the development of analytical and numerical techniques for determining NNMs, such as the method of multiple scales [START_REF] Lacarbonara | Resonant non-linear normal modes. part i: analytical treatment for structural one-dimensional systems[END_REF], modal representation approach [START_REF] Arquier | Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes[END_REF] and optimisation technique [START_REF] Slater | A numerical method for determining nonlinear normal modes[END_REF]. While in practice the analytical methods are limited to the systems with only a small number of DOFs as they need the explicit form of the equations of motion, the numerical techniques have become more common due to their applicability for large-scale non-linear models. The current most popular methodology is based on the shooting scheme [START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF][START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF] which iteratively improves the initial estimates for the non-linear mode-shape and respective frequency by enforcing the free response of the conservative system to satisfy the periodicity condition. The pseudo-arclength continuation method [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF] can also be utilised to elucidate the evolution of the oscillation frequency with the variation of the system's energy in the frequency-energy plots and capture the internal resonance phenomenon (i.e. modal interactions) as emanating tongues from the NNM branch. It should be noted that each iteration of the shooting scheme relies on the time integration of the non-linear dynamic setting for the prescribed initial condition over the corresponding period and consequently, refining the initial conditions based on the Newton-Raphson method.

Thus, the accuracy of the tools that are used for structural dynamic analysis directly affects the efficiency of the NNM computation process. For models with smooth non-linearities, a great deal of research has been carried out on treating the second-order equations of motion in the time integration step and a number of reliable numerical schemes are available (e.g. Newmark family of methods [START_REF] Géradin | Mechanical vibrations: theory and application to structural dynamics[END_REF]).

On the other hand, despite the importance of the second family of nonlinearities, only a few notable contributions have been devoted to the detailed understanding of the corresponding dynamics. Furthermore, when the nonlinearities involved are non-smooth, even small systems can lead to complex behaviour with no analogy in the linear theory [START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF]. This could be the main reason for the limited number of studies dedicated to the non-linear modal analysis of the classical structural components under unilateral constraints.

However, non-smooth constraints arise in many practical situations. Numerous examples have been already identified, such as critical velocity of the moving load applied on the beam-like components resting on the unilateral substrate [START_REF] Attar | A novel technique for dynamic analysis of beam-like structures on tensionless elastic foundations subjected to moving loads[END_REF], moored floating systems where the stiffness of the mooring component is highly influenced by the tether operating in the taut-slack state [START_REF] Lu | Nonlinear dynamics of a submerged floating moored structure by incremental harmonic balance method with FFT[END_REF], suspension bridges which are subjected to the slackness response of hangers, granular media with the non-smooth dynamics due to the contact between the grains (beads) [START_REF] Jayaprakash | Nonlinear normal modes and band zones in granular chains with no precompression[END_REF], delamination of thin coating films adhered to elastic substrates in flexible electronic devices [START_REF] Mei | Buckling modes of elastic thin films on elastic substrates[END_REF] and complex unilateral behaviour of soil in the offshore riser-seabed interactions [START_REF] Neto | Structural stability of flexible lines in catenary configuration under torsion[END_REF]. The primary objective of the present study is to develop a numerically efficient scheme for dynamic analysis of the structural members subjected to the unilateral elastic constraints.

Herein, the unilateral (also described as damaged, one-sided, bimodular or bilinear) term is used to characterise the constitutive law of the restoring force model in the constraining elements as they only sustain contact forces in one direction, either compressive or tensile. For example, they only sustain compression (i.e. repulsive forces) if they are considered as a substrate and tension (i.e. attractive forces) if they are considered as hangers. We focus on the unilateral constraints with the linear elastic contact law at two regions, before and after contact. A detailed understanding of the free and forced vibrations in such a system is a necessary initial step towards the characterisation of nonsmooth dynamics in many physical systems, such as prismatic members (i.e. bars, beams and rectangular plates) subjected to the unilateral constraints.

A number of studies exist for the static and dynamic analysis of structuresupport problem with unilateral connection (e.g. [START_REF] Coşkun | The response of a finite beam on a tensionless pasternak foundation subjected to a harmonic load[END_REF][START_REF] Zhang | Response of a finite beam in contact with a tensionless foundation under symmetric and asymmetric loading[END_REF][START_REF] Silveira | Nonlinear analysis of structural elements under unilateral contact constraints by a ritz type approach[END_REF][START_REF] Ma | Response of an infinite beam resting on a tensionless elastic foundation subjected to arbitrarily complex transverse loads[END_REF][START_REF] Giannini | Experimental evidence of bifurcating nonlinear normal modes in piecewise linear systems[END_REF][START_REF] Casini | Persistent and ghost nonlinear normal modes in the forced response of non-smooth systems[END_REF]), but, mostly, they assume that the location of the contact and non-contact zones are given as prior information, or they are limited to the systems with limited DOFs. An indepth understanding is still incomplete, as the ongoing challenge is to perform the time integration procedure in a predictive and efficient way. The difficulty stems from the inherent non-smoothness of the time evolution that leads to accumulative error. In this work, we use the idea that was first introduced in [START_REF] Attar | A novel technique for dynamic analysis of beam-like structures on tensionless elastic foundations subjected to moving loads[END_REF] and develop a stable integration method by reducing the continuous non-smooth system to a piecewise-smooth multi-degree of freedom (MDOF) model. Therefore, the trajectory of each node in the space-discretised model will be a sequence of smooth solutions with switching time instants between them. In order to minimise the accumulative error, the transition instants that the nodes cross the boundary between the subdomains must be accurately captured. For this reason, a bisection technique can be embedded in the generalised α-method to detect and localise these moments during the time integration procedure. Although the NNM frequency of the systems with no initial gap between the solid body and unilateral constraint is energy independent [START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF], we utilise the developed dynamic solution and the shooting technique at a fixed level of energy to follow the NNM branch and present the evolution of the periodic motions as a function of the contact stiffness. These plots can provide useful information about the system non-linear resonance frequencies and modal properties which cannot be predicted using the linear and bilinear theories. This scenario is very similar to the frequency-energy plots in the smooth non-linear systems and it would help us to study salient features of the computed NNMs. The local stability of these non-linear modes can also be assessed by the Lyapunov theory. Our analysis suggests that the loss of stability on the computed branch can be associated with a bifurcation point. This may happen due to the modal interaction with the higher modes (internal resonance phenomenon), or evolution of new non-linear modes from the main branch with periods equal to the integer multiples of the fundamental period.

The outline of this paper is as follows. In Section 2, we first present the class of mechanical systems under study and the form of corresponding governing equations of motion. It is also elaborated how to reduce the non-smooth dynamic setting into a piecewise-smooth discrete model. Then, we present the basic steps of the unconditionally stable time integration scheme and the so-lution algorithm. In Section 3, we describe how to implement the developed numerical method in the shooting technique and obtain the non-linear vibrational modes of the system. We also apply the path following method and stability analysis to study non-smooth dynamics of the chain of oscillators with an elastic unilateral obstacle. Section 4 is devoted to demonstrating the computational efficiency of the presented approach on typical numerical examples. Finally, the concluding remarks are presented in Section 5.

2 Non-smooth dynamics of the structures subjected to the unilateral elastic contact

Continuous model

Several structural components in dynamic systems can be represented as one-dimensional (1D) elastic members, such as bars and beams, or two-dimensional (2D) thin structures, such as plates and shells. These components may undergo longitudinal, torsional, and lateral vibrations while they are in contact with other interconnecting deformable element(s) as the unilateral constraint(s). The non-smooth connection between the structure and support(s) permissible in such configurations can give rise to non-negligible nonlinearities, even if the material properties remain linear and the assumption of the elastic behaviour within the range of small deformation theory remains valid. In such a case, the general form of the equations of motion for the interaction of the deformable body and elastic support(s) can be obtained by the use of classical theories and the Heaviside step function to describe the non-smooth nature of the contact constraint [START_REF] Attar | A novel technique for dynamic analysis of beam-like structures on tensionless elastic foundations subjected to moving loads[END_REF].

Consider a straight homogenous elastic bar where its right end is subjected to the unilateral constraint due to the presence of an elastic wall, as shown in Fig. 1(a). This bar with the length of L, cross-sectional area of A and density of ρ is characterised by its centroidal axis. Hence, the displacement field at any point on the centreline is composed of the longitudinal displacement u (x, t) where x is the axial coordinate and t is time. The unilateral support is also characterised with the stiffness k w and damping coefficient c w (in SI units, k w and c w for the bar-support system are expressed by N/m and N.s/m, respectively) where its unstressed position is shifted by the gap of δ from the right end of the unstressed bar. The magnitude of the contact force can be determined by the penetration of the bar end into the constraint, i.e. u (L, t)-δ (if any penetration occurs), and the time derivative of the end displacement, i.e. ∂u (L, t) /∂t, as

F c (t) = c w ∂u (L, t) ∂t + k w [u (L, t) -δ] • H u (L, t) -δ, ∂u (L, t) ∂t (1)
where H (ε, ∂ε/∂t) is the non-smooth contact function to represent the unilateral nature of the constraint and it is expressed as

H(ε, ∂ε ∂t ) = Θ(ε) • Θ(k w ε + c w ∂ε ∂t ) (2)
where Θ(ε) is the Heaviside step function which is defined as Θ(ε) = 1 when ε ≥ 0 and Θ(ε) = 0 when ε < 0. Note that the contact function H (ε, ∂ε/∂t) prescribes the non-smooth property of the constraint and it verifies that the contact force F c acts axially on the bar at point x = L, if and only if, u (L, t) ≥ δ and k w u (L, t) + c w ∂u (L, t) /∂t ≥ 0 (i.e. the contact force is always compressive). Using the mathematical description of the contact force from Eq. 1 and the definition of the Dirac delta function δ (•) to apply this non-smooth force at the bar end (x = L), the equation of motion for the axial displacement field u (x, t) can be identified as

ρA ∂ 2 u (x, t) ∂t 2 -EA ∂ 2 u (x, t) ∂x 2 = F (x, t) -F c (t) δ (x -L) (3)
where E is the Young's modulus and F (x, t) is the externally applied force on the component. It is important to note that Eqs. 1 and 3 are obtained by assuming that the unilateral constraint (obstacle) is always at rest at the moment that the contact occurs. Therefore, we can ignore the free motion of the obstacle. This assumption is reasonable for small values of damping [START_REF] Leine | Dynamics and bifurcations of non-smooth mechanical systems[END_REF] and we use it in present study. For high levels of damping, two separate sets of differential equations should be constructed to govern the behaviour of the bar and obstacle independently (i.e. a second order equation for the bar axial displacement and a first order one for the obstacle deformation), and the contact and non-contact states will be regulated by the contact load at their interface. There is another approach to describe the non-smooth character of the contact mechanism which is utilised in [START_REF] Andreaus | Numerical simulation of the soft contact dynamics of an impacting bilinear oscillator[END_REF] for a 1-DOF model by considering H(ε, ∂ε/∂t) = Θ(ε) and consequently, assuming that the discontinuity of both of the elastic and damping contact forces are always concentrated at u (L, t) = δ. It is also obvious that the contact function reduces to the Heaviside step function for the undamped constraints (c w = 0). Moreover, by setting H(ε, ∂ε/∂t) = 1, Eq. 3 will be reduced to the bar interacting with the conventional bilateral obstacle (smooth model).

We can pursue the same approach to derive the governing equation of motion for the lateral vibrations of beam-like structures subjected to the unilateral supports, as illustrated in Fig. 1(b). Utilising the Timoshenko theory (i.e.

considering the effect of shear deformation and rotary inertia with no axial deformation), the displacement field at any point on the beam centreline is composed of the lateral displacement w (x, t) and rotation of the cross-section ϕ (x, t). In order to describe the homogenous unilateral support, we can use a one-way viscoelastic Winkler model with the stiffness of k w and damping coefficient of c w (note that for the case of beam-support system, the units of k w and c w are N/m 2 and N.s/m 2 in SI units, respectively) where its unstressed position is shifted by the gap of δ from the unstressed configuration of the beam section that is subjected to the contact. In a similar fashion as the previous system, the magnitude of the contact force per unit length can be determined by the lateral penetration of the beam into the constraint, i.e.

w (x, t)δ, and the time derivative of the lateral displacement ∂w (x, t) /∂t as

P c (x, t) = c w ∂w (x, t) ∂t + k w [w (x, t) -δ] • H w (x, t) -δ, ∂w (x, t) ∂t (4)
The physical mechanism of such a constraint can be illustrated by considering, for instance, a hanger connected to the beam at point x where the slackness property of the hanger switches on the contact force as soon as the hanger becomes under tension. Therefore, the governing differential equations at point

x of the centreline can be expressed as (with the damping of the constraint as the only source of dissipation)

ρA ∂ 2 w ∂t 2 -κ s GA ∂ 2 w ∂x 2 + ∂ϕ ∂x = F (x, t) -P c (x, t) (5) ρI ∂ 2 ϕ ∂t 2 -EI ∂ 2 ϕ ∂x 2 + κ s GA ∂w ∂x + ϕ = µ (x, t) ( 6
)
where G is the shear modulus, κ s is the Timoshenko correction factor, I is the second moment of area about the y-axis, F (x, t) is the externally applied load on the component per unit length, and µ (x, t) is the externally applied moment per unit length. Once more, we note that in the dynamic equations of motion, the contact function is the only source of non-linearity. This function introduces inequalities into the system, which are the multi-valued op-erators and they establish the length and number of the contact regions. In other words, the contact conditions for a particular x are ambiguous as it is unknown beforehand whether the respective displacement and velocity components satisfy H (wδ, ∂w/∂t) = 1 or H (wδ, ∂w/∂t) = 0. This gives rise to the non-smooth dynamics of the system even for the present linear elastic structure and substrate. We can also realise that the parameters of the unilateral constraint (i.e. k w , c w and δ) can be specified locally at point x.

Thus, multiple unilateral constraints may be considered over the length of the component with variable properties with respect to x. For instance, Fig. 1(b) depicts a cantilever component where for x ∈ [0, Lx f ] there is no constraint and for x ∈ [Lx f , L] the beam is subjected to the unilateral support with the parameters of k w , c w and δ.

Reducing the Timoshenko formulation to a system with an Euler-Bernoulli beam as the slender member is straightforward. To this end, the role of shear deformation and rotary inertia are neglected in Eqs. 5 and 6 by making EI/ (κ s GA) → 0 and ρI → 0 (note that we can also obtain the formulation corresponding to the Rayleigh beam theory by keeping the rotary inertia term and only neglecting the shear deformation). We can extend the same approach for non-smooth systems involving 2D structural members such as plates and shells. For instance, we can consider the Kirchhoff plate theory [START_REF] Timoshenko | Theory of plates and shells[END_REF] for dynamic analysis of a rectangular thin plate subjected to the unilateral Winkler constraint (see Fig. 1(c)) where the plate deformation is described by the lateral displacement of its middle surface w (x, t). Within the range of small deformation, the lines vertical to the middle surface of the plate thickness remain straight and vertical to the deformed middle surface. In a similar fashion as the beam component, the contact function can be used again in the classical theory to consider the distributed contact reaction per unit area as a non-smooth unilateral load. Thus, the governing non-smooth equation of motion for the isotropic thin plate of length L, width c and thickness b can be expressed as

Eb 3 12 (1 -ν 2 ) ∂ 4 w ∂x 4 + 2∂ 4 w ∂x 2 ∂y 2 + ∂ 4 w ∂y 4 + ρb ∂ 2 w ∂t 2 = F (x, y, t) -c w ∂w ∂t + k w [w -δ] H(w -δ, ∂w ∂t ) ( 7
)
where ν is the Poisson's ratio, F (x, y, t) is the externally applied force on the plate per unit area and in SI units, the foundation parameters for the platesupport system are expressed by N/m 3 and N.s/m 3 for k w and c w , respectively.

Discretisation

Unlike for linear systems, the superposition principle and the direct application of modal expansion method [START_REF] Attar | Dynamic response of cracked timoshenko beams on elastic foundations under moving har-51 monic loads[END_REF] cannot solve the governing equations of non-smooth systems. Note that multiple/global constraints may exist in the model (e.g. when the lateral motion of the structure is fully constrained by the unilateral substrate) and thus, the combination of the modal expansion method and a non-linear iterative solver such as Newton-Raphson method is also inappropriate. In such a scheme, all the points where the non-smooth support reaction is present should be captured in the solution [START_REF] Papalukopoulos | Dynamics of large scale mechanical models using multilevel substructuring[END_REF] and it makes the number of generalised coordinates and consequently, the number of modes required to reach a converged solution very high. Moreover, it involves the stability and convergence issues where this method may not converge toward the desired solution, or any solution. Therefore, we are looking for a general stable solution method that can be utilised for any length of the non-smooth contact area between the flexible body and support. Subsequently, this method can be employed as the time integration step in the computation of NNMs.

It is important to note the discontinuity boundary for the non-smooth behaviour divides the space into two adjacent regions for point x (or point (x, y) of the plate). In each of these regions, the system's equations of motion are smooth and crossing the border is the source of non-smoothness. This implies that the non-smooth continuous system can be partitioned into several segments where the behaviour of each segment is piecewise linear. Here, we utilise the discrete lattice method (DLM) [START_REF] Attar | Free vibration analysis of a cracked shear deformable beam on a two-parameter elastic foundation using a lattice spring model[END_REF][START_REF] Griffiths | Modelling of elastic continua using a grillage of structural elements based on discrete element concepts[END_REF] to discretise the system into a network of identical particles. Based on the DLM, these particles interact with their neighbours through the connecting elements between them. These connecting elements can be springs [START_REF] Attar | Free vibration analysis of a cracked shear deformable beam on a two-parameter elastic foundation using a lattice spring model[END_REF][START_REF] Pasternak | On the possibility of elastic strain localisation in a fault, pure and applied geophysics[END_REF][START_REF] Pasternak | Generalised homogenisation procedures for granular materials[END_REF], truss elements [START_REF] Beex | A multiscale quasicontinuum method for dissipative lattice models and discrete networks[END_REF], beam elements [START_REF] Nikolic | Rock mechanics model capable of representing initial heterogeneities and full set of 3D failure mechanisms[END_REF] or a combination of them. The interactions between the discrete units at small scale allow us to reproduce the macro-mechanical response of the continuous system. In this study, we use a combination of the bar elements and springs in the DLM for the dynamic analysis of the bar-support system in such a way that each particle has one DOF, i.e. {u i (t)} with u i as the longitudinal displacement in x direction. The lumped masses are connected through bar elements and also, the unilateral elastic obstacle is modelled by the oneway spring-damper element, as shown in Fig. 1(a). The contact element is only active when the corresponding condition is satisfied, i.e. u n (t) ≥ δ and

K N u n (t) + C N ∂u n (t) /∂t ≥ 0,
where n is the number of particles (nodes) and K N and C N are the constraint constants in the DLM.

The same strategy can be extended to study the non-smooth dynamics of beam-like components and 2D thin structures. In the discrete model of the beam-substrate system each particle has two DOFs, i.e. {w i (t) , ϕ i (t)} with w i as the lateral displacement in z direction and ϕ i as the rotation about y-axis.

The lattice nodes are connected through the beam elements and the nonsmooth support(s) is represented by a chain of parallel mutually independent spring-damper elements (see Fig. 1(b)). These contact spring-dampers are switched on when H (w iδ i , ẇi ) = 1 where δ i is the gap length for particle i and the over-dot denotes the time derivative. In a similar way, the DLM of the 2D thin structure involves FEM plate elements and one-way spring-dampers, as illustrated in Figs. 1(c) and (d). It is important to note that particles of the plate-support model have one more DOF, compared to the beam structure, which is the rotation about the x-axis (θ) and it makes the DOFs of each point as {w i,j (t) , ϕ i,j (t) , θ i,j (t)}.

The relationships between the lattice model parameters at the discrete level and continuum properties at the macro scale can be obtained from the homogenisation technique [START_REF] Attar | Free vibration analysis of a cracked shear deformable beam on a two-parameter elastic foundation using a lattice spring model[END_REF]. For instance, the constraint constants in the beamsupport system can be identified as K N = k w h and C N = c w h and the particle concentrated mass and rotational moment of inertia about the y-axis are M = ρAh and J yy = (M/12) (h 2 + ˇ b 2 ), respectively, where h is the length of the unit cells, b is the thickness of the component with rectangular crosssection, ˇ = 1 for the Timoshenko theory and ˇ = 0 for the Euler beam model.

The stiffness matrix of the beam elements which relates the nodal displacements of particle i to i + 1 is also expressed as [START_REF] Reddy | On the dynamic behaviour of the timoshenko beam finite elements[END_REF] k = 2EI

h 3 + 12ǧh        6 -3h -6 -3h -3h 2h 2 + 6ǧ 3h h 2 -6ǧ -6 3h 6 3h -3h h 2 -6ǧ 3h 2h 2 + 6ǧ        (8)
where ǧ = 0 for the Euler beam model and ǧ = EI/ (κ s GA) for the Timoshenko theory. We can now derive the non-linear discrete equations of motion for the particles by applying the Lagrange equation as

d dt ∂L (t) ∂ Ḋq - ∂L (t) ∂D q + ∂R (t) ∂ Ḋq = F q (t) (9)
where L (t) = T (t) -U (t) is the so-called Lagrangian function of the system, with T (t) and U (t) being the kinetic energy and potential energy, respectively, R (t) is the Rayleigh dissipation function, F q (t) is the q-th member of the generalised load vector F (t) ∈ R N applied on the system, N = ñp • n is the total number of DOFs in the discrete system when each node has ñp degrees of freedom (ñ p = 1 for the bar elements, while for the beam-support and plate-support systems ñp = 2 and ñp = 3, respectively), and D q is the q-th member of the general displacement vector D (t) ∈ R N . The kinetic energy, potential energy and Rayleigh functions of the smooth system can be expressed respectively as the following

T (t) = 1 2 ḊT (t) • M • Ḋ (t) (10)
U (t) = 1 2 D T (t) • K • D (t) + 1 2 K N D T (t) -G T • [D (t) -G] (11)
R (t) = 1 2 C N ḊT (t) • Ḋ (t) (12)
herein, (•) T denotes the transpose operator, M ∈ R N ×N is the diagonal mass matrix with the lumped mass and rotational inertia of the particles on the main diagonal, K ∈ R N ×N is the symmetric stiffness matrix containing the contribution of the connecting elements (e.g. for beam-support systems, it can be obtained by assembling Eq. 8 into the global form), and G ∈ R N is a vector containing the length of initial gaps between the unstressed supports and lattice nodes (note that G has n nonzero components where n ≥ 0 is the number of DLM nodes that are subjected to the constraints and it depends on the system configuration). The second-order equation of motion for the smooth discrete model can be obtained in the matrix form by using Eqs. 10-12 in Eq. 9 as

M • D (t) + C N Ḋ (t) + K • D (t) + K N [D (t) -G] = F (t) (13)
In order to include the unilateral contact effects, we introduce the contact constraint matrix H (t) ∈ R N ×N which is a diagonal matrix with n nonzero members on the main diagonal containing the contact function of the interacting lattice nodes (note that, like G, the position of n nonzero components on the main diagonal of H (t) also depend on the system configuration). Accordingly, P (t) ∈ R N , which is the vector of general discrete unilateral forces developed at the contact points, is characterised by the following equation

P (t) = C N H (t) • Ḋ (t) + K N H (t) • [D (t) -G] (14)
Therefore, the second-order equation of motion for the non-smooth discrete model can be obtained as

M • D (t) + C N H (t) • Ḋ (t) + K • D (t) + K N H (t) • [D (t) -G] = F (t) (15)
For example, if a beam is interacting with a unilateral substrate through its full length (i.e. N = 2n and n = n) , M, D (t), G, H (t) and F (t) will have the following forms

M = diag [M, J yy , M, J yy , • • • , M, J yy ] ( 16
)
D T (t) = {w 1 , ϕ 1 , w 2 , ϕ 2 , • • • , w n , ϕ n } (17)
G T = {δ 1 , 0, δ 2 , 0, • • • , δ n , 0} (18)
H (t) = diag H (w 1 -δ 1 , ẇ1 ) , 0, H (w 2 -δ 2 , ẇ2 ) , 0, • • • , H (w n -δ n , ẇn ) , 0 (19)
F T (t) = {f 1 , µ 1 , f 2 , µ 2 , • • • , f n , µ n } (20)
where f i and µ i are the nodal external force and moment applied on point i.

Eq. 15 can help us to treat the non-smooth structural dynamics of the problem by the corresponding discrete piecewise smooth model. This model partitions the state-space R 2 of each lattice node, which is subjected to the unilateral constraint, into two neighbouring subdomains within which the trajectory of the node is smooth. Accordingly, when n is large enough, the mechanical behaviour of the continuous system can be simulated with high accuracy using the DLM. For instance, the equations of motion for particle i in the beamsubstrate system can be obtained from Eq. 15 as

M ẅi + C N ẇi H (w i -δ i , ẇi ) + K N {w i -δ i } H (w i -δ i , ẇi ) + 6EI h 3 + 12ǧh {-2w i-1 + hϕ i-1 + 4w i -2w i+1 -hϕ i+1 } = f i ( 21
)
J yy φi + 2EI

h 3 + 12ǧh { -3hw i-1 + h 2 -6ǧ ϕ i-1 + 4h 2 + 12ǧ ϕ i + 3hw i+1 + h 2 -6ǧ ϕ i+1 } = µ i ( 22
)
Note that for applying Eq. 2 to calculate the contact function at discrete level, we use the well-known property of the Heaviside function as Θ(k

w ε/ ȟ + c w ε/ ȟ) = Θ(k w ε + c w ε) for any ȟ > 0.
In order to establish the relation between the continuum properties and DLM parameters and also, quantify the role of h as the characteristic length of the discrete model, we can utilise the homogenisation technique to obtain the continuous form of Eqs. 21 and 22

and compare them to Eqs. 5 and 6. To this end, the discrete displacement components w i (t) and ϕ i (t) are replaced by the respective continuous functions as w (x, t) and ϕ (x, t). Therefore, the displacements of the neighbouring points can also be substituted by the corresponding truncated Taylor series expansions [START_REF] Attar | Free vibration analysis of a cracked shear deformable beam on a two-parameter elastic foundation using a lattice spring model[END_REF].

Time integration algorithm

Upon knowing the boundary conditions in terms of external forces or kinematics, the smooth counterpart of Eq. 15 could be solved with classic numerical techniques. Moreover, coupled with the well-known non-linear iterative solvers (e.g. Newton-Raphson), one may be able to obtain the solution of a specific problem. However, as the inherent mathematical issues (e.g. convergence problems and/or the accumulative error due to the non-smoothness) may arise in the time integration procedure, the calculation of NNMs can be infeasible when the solution involves theses typical integration instabilities. It is also noteworthy that the application of the reduced model techniques is not computationally effective here, as these methods are efficient only if there is a local contact constraint in the system (i.e. n << n), which is not generally valid for the present study. In this section, we propose an alternative algorithm that is based on tracking the trajectories of the DLM particles, localisation of the events that they cross the constraint boundary, and successively capturing the back-to-back intervals that they spend in the smooth subdomains.

Taking advantage of the piecewise linear nature of the problem, the Jacobian matrix of the shooting method process can also be determined concurrently with the system solution. This outstanding numerical property reduces the computational cost of computing the NNMs compared to the finite difference scheme where the components of Jacobian matrix are evaluated numerically by perturbing each DOF.

We solve the differential equilibrium equation of motion Eq. 15 by combined use of the generalised-α time integration scheme [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method[END_REF] and the bisection method.

To this end, an appropriate time increment is assigned as ∆t = t r+1t r to discretise the time domain where the subscript r is the time step number starting from r = 0. The first step is to compute the initial acceleration vector D0 at t = t 0 by means of the equation of motion and known initial conditions, i.e. initial displacement D 0 , velocity Ḋ0 and external applied load F 0 , as (the initial contact constraint matrix H 0 can be calculated from D 0 )

D0 = M -1 • F 0 -C N H 0 • Ḋ0 -K • D 0 -K N H 0 • [D 0 -G] (23)
We use the generalised-α method [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method[END_REF] to rewrite Eq. 15 in a general midpoint of the time step as (note that in the classical Newmark method, the respective equation is considered at the end-point of the time step)

(1 -α m ) M • Dr+1 + α m M • Dr + (1 -α f ) C N H r+1 • Ḋr+1 + α f C N H r • Ḋr + (1 -α f ) K • D r+1 + α f K • D r + (1 -α f ) K N H r+1 • [D r+1 -G] + α f K N H r • [D r -G] = (1 -α f ) F r+1 + α f F r (24)
Moreover, the generalised-α method uses the Newmark difference approximation equations as

D r+1 = D r + ∆t Ḋr + (∆t) 2 2 (1 -2β) Dr + 2β Dr+1 (25)
Ḋr+1 = Ḋr + ∆t (1 -γ) Dr + ∆tγ Dr+1 (26)
The parameters for this unconditionally stable integration scheme in Eqs. 24-26 are defined as [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method[END_REF]

α m = 2ρ ∞ -1 ρ ∞ + 1 , α f = ρ ∞ ρ ∞ + 1 , γ = 1 2 -α m + α f , β = 1 4 (1 -α m + α f ) 2 (27)
where ρ ∞ ∈ [0, 1] is the spectral radius of the method. We note that choosing 0 ≤ ρ ∞ < 1 introduces numerical dissipation in the integration procedure, while ρ ∞ = 1 makes it conservative. Thanks to the generalised-α one-step scheme, we only require data at t = t r (i.e. one single previous time step) to obtain the system solution at t = t r+1 . Therefore, by using Eqs. 24-26 and after some algebra, the displacement vector at the end of the time step can be obtained from the following equation

K eff r+1 • D r+1 = F eff r+1 ( 28
)
where K eff r+1 and F eff r+1 are identified as

K eff r+1 = (1 -α f ) (K + K N H r+1 ) + (1 -α f ) γ β (∆t) C N H r+1 + 1 -α m β (∆t) 2 M (29) F eff r+1 = (1 -α f ) F r+1 + α f F r + (1 -α f ) K N H r+1 • G + α f K N H r • G + (1 -α m ) M • 1 β (∆t) 2 D r + 1 β (∆t) Ḋr + ( 1 2β -1) Dr + (1 -α f ) C N H r+1 • γ β (∆t) D r + ( γ β -1) Ḋr + ∆t( γ 2β -1) Dr -α m M • Dr -α f C N H r • Ḋr -α f (K + K N H r ) • D r (30)
Then, the acceleration and velocity of the system can be computed as

Dr+1 = 1 β (∆t) 2 (D r+1 -D r ) - 1 β (∆t) Ḋr -( 1 2β -1) Dr (31)
Ḋr+1 = Ḋr + ∆t (1 -γ) Dr + ∆tγ Dr+1 (32)
We can now proceed and obtain the transient step-by-step solution of the system from Eqs. 28, 31 and 32. The application of the Heaviside step function in the contact constraint matrix of the system (H) will enable the solution procedure to automatically detect the position of the lattice points with respect to the unilateral elastic constraint, and accordingly update the system stiffness matrix. Therefore, if we assume that no lattice node crosses the boundary over ∆t (i.e. H r+1 = H r ), the system response can be obtained at t = t r+1 from the input data at t = t r . Once the crossing event occurs and the contact condition is satisfied for any of the lattice nodes at t = t r+1 , the system solution activates the corresponding constraint elements (i.e. the discrete spring-damper components) in the subsequent time step and considers their contribution in the effective stiffness matrix and force vector (Eqs. 29 and 30). It is important to note that the time step must be appropriately small to enable this direct time integration procedure to accurately localise the crossing events and minimise the numerical error due to the approximation of the crossing time instants (the main source of error comes from approximating the boundary-crossing time instants by t = t r+1 , while they occur at t ∈ [t r , t r+1 ]). In practice however, the error due to this approximation would increase progressively at every occurrence so that the accuracy of the following estimate would gradually decrease. As a consequence, it can lead to an accumulative error at the end of the simulation and may significantly affect the resolution of the results. One possible method to resolve this issue is to choose a smaller ∆t and perform convergence study to ensure that the solution is independent of the length of the time step. However, in some cases it may impractically increase the overall computational time.

To minimise the accumulative error, we can also embed the bisection algorithm into the integration procedure to capture the time instants that the lattice nodes impact the constraint or recover the gap length at the end of their contact period. To this end, we proceed with the integration procedure and obtain D r+1 from Eq. 28 with the assumption that no particles cross the boundary over ∆t (i.e. H r+1 = H r ). Then, we calculate H r+1 from the resulting displacement vector to verify this precondition. The solution at t = t r+1 can be accepted as the system dynamic response (i.e. we are allowed to go to the next time step) if the resulting contact constraint matrix verifies our assumption, i.e.

H r+1 = H r . But, if any event is detected (i.e. H r+1 = H r ),
it implies that a number of lattice points have crossed the boundary over ∆t where they can be localised by comparing the position of unequal members on the main diagonal of H r+1 with H r . Therefore, we locally break down the current time step [t r , t r+1 ] into a number of subintervals to compute the time instants that the events occur as

t s m ∈ t s 1 = t r , t s 2 , t s 3 , t s 4 , • • • , t s m-1 , t s m , t s m+1 , • • • , t s η-2 , t s η-1 , t s η = t r+1 (33)
where η -1 is the number of subintervals (the required number of subintervals depends on the bisection module), and the subscript m shows the counter of local time instants t s , starting from t s 1 = t r and ending with t s η = t r+1 . Thus, the bisection module is called to use Eqs. [START_REF] Silveira | Nonlinear analysis of structural elements under unilateral contact constraints by a ritz type approach[END_REF] ), all the localised lattice points have crossed the boundary and we can accept the solution at t = t r+1 = t s η to exit the local subintervals. Otherwise, this procedure continues by calling the bisection module again for m = m + 1, until all the localised discrete segments cross the boundary at t = t r+1 . In order to avoid numerical instabilities, if the crossing lattice nodes at the beginning of the subintervals are in contact with the constraint (i.e. H [w i (t s m )δ i , ẇi (t s m )] = 1), the new subinterval ∆t s m is computed from the bisection search in such a way that

-ε s < K N [wi(t s m+1 )-δi]+CN ẇi( t s m+1 ) K N [D(t s m+1 )-G]+CN Ḋ(t s m+1 )
< 0 where • implies the Euclidean norm and ε s ∈ R + is the specified precision tolerance (note that this parameter is specified to localise the event of crossing the non-smoothness boundary by the lattice nodes up to a given precision. Therefore, the suitability of the selected value should be studied by the convergence analysis to make sure that it is small enough to obtain the system response with sufficient accuracy). On the other hand, when the crossing lattice nodes at the beginning of the subintervals are not in contact with the constraint (i.e. H [w i (t s m )δ i , ẇi (t s m )] = 0), the new subinterval ∆t s m is obtained in such a way that 0 ≤

w i( t s m+1 )-δi D(t s m+1 )-G < ε s if w i (t s m ) < δ i , and 0 ≤ K N [wi(t s m+1 )-δi]+CN ẇi( t s m+1 ) K N [D(t s m+1 )-G]+CN Ḋ(t s m+1 ) < ε s if w i (t s m ) ≥ δ i .

Non-linear normal modes

Consider the unforced conservative form of Eq. 15 for the non-smooth discrete system

M • D (t) + K • D (t) + K N H (t) • [D (t) -G] = 0 ( 34
)
where 0 ∈ R N is the zero vector. Based on Rosenberg's definition, and as highlighted in [START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF], the energy dependency of the NNMs in this system arises from the presence of the gap (clearance) between the structure and unilateral constraint. The physical interpretation of this fact is obvious. For instance, in the presence of the gap, if the energy is not sufficient to make the structure interact with the constraint, the system always remains within one of the smooth subdomains and the solution is linear. Thus, unlike the smooth non-linear systems where the frequency-energy dependence is the inherent property of their oscillations, the non-linear frequencies of the structure-support systems are energy invariant as long as the unilateral interaction is concentrated at the interface passing through the origin (i.e. G = 0). This means that changing the energy level (e.g. initial conditions) of the present non-smooth model with no gap does not lead to a new NNM frequency for the system. Accordingly, for λ ∈ R we have z(t, λz (0) ) = λz(t, z (0) ) where z(t, z (0) ) = D Ḋ ∈ R 2N refers to the solution of Eq. 34 with G = 0 which can be recast into the state space form as

ż = Q • z ( 35
)
where

Q ∈ R 2N ×2N reads Q = 0 I M -1 • (K + K N H) 0 (36)
herein, I ∈ R N ×N is the identity matrix and 0 ∈ R N ×N is the zero matrix. We can now use the proposed time integration algorithm as the dynamic analysis tool of the shooting algorithm to compute NNMs of the system presented by Eq. 35. Although the NNMs of this system are energy invariant, they can exhibit salient properties with no counterparts in the underlying smooth model. For instance, the number of resonance frequencies and corresponding NNMs may be more than the number of system DOFs. Furthermore, the non-linear modes of the system can be excited at frequencies which are fractional or integer multiples of the NNM frequencies (subharmonic and superharmonic modes). Another noticeable feature is the internal resonance phenomenon where the non-linear modes interact and generate new multimodal periodic motions with no analogy in the linear system. The internal resonance may occur once the system non-linear frequencies are commensurate or nearly commensurate [START_REF] Kerschen | Nonlinear normal modes, part I: A useful framework for the structural dynamicist[END_REF]. Therefore, since the deployment of the structural members onto the unilateral supports is commonly found in a variety engineering applications, a systematic study of the NNMs can greatly improve our predictive understanding of their dynamic behaviour. Herein, we briefly outline the application of the shooting methodology for computing NNMs in Appendix A where some modifications are also required to adopt it to the present nonsmooth system. The reader is referred to [START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF][START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF] for more details on this approach for systems with smooth and non-smooth non-linearities.

If we consider the interaction of the component with a homogenous unilateral support with no gap at the interface (e.g. a bar in contact with a unilateral elastic obstacle at one of its boundaries, a beam supported by a series of hangers, or a plate resting on the unilateral substrate), the resulting initial conditions together with its period define a portion of the non-smooth periodic motions branch in terms of the given contact stiffness. For a range of the contact coefficient, the backbone of the NNMs branch can be constructed by sequentially incrementing the contact coefficient and using the result of each step as the initial assumption of the shooting scheme in the following step. This plot can provide useful information about the system nonliear frequencies and corresponding mode-shapes in a similar fashion as the frequency-energy plot in the smooth non-linear systems [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF]. For illustration purposes, consider the chain of oscillators in Fig. 2 with D T (t) = {u 1 , u 2 , u 3 } (i.e. n = 3, ñp = 1 and N = 3). This chain is interacting with the unilateral elastic obstacle with the stiffness of k w = k n . This simple but insightful example can represent the nonsmooth axial vibration of a thin bar (considering the stiffness of the connecting elements as k x = EA/h), such as a drill-string, with non-smooth contact at its end. The corresponding unforced conservative equations of motion can be obtained from Eq. 34 as

m 0 0 0 m 0 0 0 m • ü1 (t) ü2 (t) ü3 (t) + 2kx -kx 0 -kx 2kx -kx 0 -kx kx • u 1 (t) u 2 (t) u 3 (t) + k n 0 0 0 0 0 0 0 0 H(u 3 -δ) • u 1 (t) u 2 (t) u 3 (t)-δ = 0 0 0 (37)
Three natural frequencies of the underlying linear system with k n = 0 can be easily obtained from the eigenvalue analysis as ω u 1 = 0.445 k x /m, ω u 2 = 1.247 k x /m and ω u 3 = 1.802 k x /m. These frequencies and corresponding linear normal modes can be used as the start-points in computing the NNMs for a range of the unilateral contact stiffness k n (note that the superscript u in ω u i denotes the unilateral interaction). We perform the non-smooth modal analysis of this system with δ = 0 by sequentially increasing k n and using the shooting method. The NNMs are depicted for 0 ≤ k n /k x ≤ 10 in Fig. 3. The first NNM is presented in more details where the incremental evolution of the NNMs for points (a)-(l) on its backbone is demonstrated as the time-history of the system free response over the first period T u 1 = 2π/ω u 1 in Fig. 4(a) and the configuration space, i.e. u 1 (t) vs. u 3 (t) and u 2 (t) vs. u 3 (t), in Fig. 4(b). For small values of k n /k x , the NNMs are close to the linear shape, i.e. they exhibit linear symmetric configuration and simultaneously pass through the origin (see Figs. [START_REF] Attar | Dynamic response of cracked timoshenko beams on elastic foundations under moving har-51 monic loads[END_REF](a.a) and (a.b)). By increasing the contact stiffness, the NNMs lose their symmetric linear shape and do not pass through zero at the same instant. There are some other important considerations that must be taken into account. Firstly, the bilinear approximation of the non-linear frequencies loses its accuracy at higher values of the contact stiffness (see the orange line in Fig. 3). Secondly, the stability analysis of the NNMs reveals that additional periodic motions may appear in the non-smooth system once the NNMs on the branch become unstable (see the new periodic motions that are born at points (d) and (i) on the backbone of the first NNM in Fig. 3(a)). The well-known bilinear formula has been used in previous studies [START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF][START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF][START_REF] Chati | Modal analysis of a cracked beam[END_REF] to approximate the non-linear frequency of the system with the zero gap as

ωa = 2ω - a ω + a ω - a + ω + a (38)
herein, ωa is the bilinear frequency of the a-th mode, and ω - a and ω + a are the resonance frequencies before and after making contact with the unilateral constraints (i.e. ω -= K -/M and ω + = K + /M). Zuo and Curnier [START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF] believe that the accuracy of this formula depends on the commutativity condition of stiffness matrix i.e. the bilinear formula leads to accurate results if

K -• K + = K + • K -
where K -and K + are the system stiffness before and after making contact with the support, respectively. The present 3-DOF system is not commutative as the contact stiffness partially contributes to the main diagonal of the global stiffness matrix (see Eq. 37). Thus, it is not surprising that the bilinear relation cannot provide accurate approximations for the NNMs, specially for higher modes, as it is shown in Figs. 3(b) and 3(c).

The stability of the computed non-linear modes can also be studied in the linear sense by computing its characteristic multipliers (also known as Floquet multipliers) [START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF][START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF]. Based on the Lyapunov theory, the periodic solution z(t, z (0) ) is stable if for an arbitrary infinitesimal disturbance of the initial conditions ( d 0 z (0) ) the distance between the perturbed solution and initial conditions, i.e. d(t) = z(t, z (0) + d 0 )z (0) , tends to zero when t → ∞. This means that the non-linear mode is stable when the eigenvalues (characteristic multipliers) of the monodromy matrix have magnitudes less than one. The monodromy matrix can be identified as φ T (z (0) ) = ∂z(T,z (0) )

∂z (0)
.

It is worth mentioning that the components of φ T (z (0) ) have already been evaluated at the last iteration of the Newton-Raphson algorithm during the numerical computation of NNMs as g (t = T ) for constructing the Jacobian matrix (see details of the shooting scheme in Appendix A, in particular Eqs.

A.5 and A.10). Therefore, the stability analysis of the computed NNM can be simply performed by calculating the eigenvalues of g (T ) after the last step of the iterative loop. As long as the eigenvalues are located inside the unit circle in the complex plane, the computed NNM is stable. If any of the characteristic multipliers moves towards the outside of the unit circle, the corresponding NNM becomes unstable. The unstable NNMs in Fig. 3 are shown with the dashed line. Interestingly, the stability analysis reveals that when the nonlinear modes become unstable on the branch, new periodic motions are born with the periods equal to the integer multiples of the fundamental period.

This behaviour is similar to the well-known flip bifurcation that occurs in the non-linear systems to switch to a new dynamic response with twice the period of the original system [START_REF] Seydel | Practical bifurcation and stability analysis[END_REF]. When the instability threshold for the onset of the bifurcation is identified by the system Floquet multipliers, the new periodic motion can be predicted by applying a perturbation technique.

For the first NNM of the present system, the first instability occurs at point (d) on the backbone (see Fig. 3(a)) where the new non-smooth periodic motion is born with the period doubling (T u 1 → 2T u 1 ). The incremental evolution of this newly born NNM with the dimensionless stiffness k n /k x is qualitatively different from the backbone branch, as illustrated in Fig. 5 (both the time response and configuration space). The second bifurcation also occurs at point (i) where the new branch is formed with the period equal to the quintuple of the fundamental period (T u 1 → 5T u 1 ) and the system traces a different trajectory. The corresponding evolution of this branch is also presented in Fig. 6. It is important to note that the new NNM branches can also switch from stable to unstable modes where at the location of the instabilities again new periodic solutions arise. For the sake of clarity and conciseness, other branches of higher levels are not presented here.

Another distinct feature of the NNMs is the internal resonance phenomenon where the non-linear modes may interact with the system's higher modes once their frequencies are commensurate. This interaction can be observed at higher resolutions of the frequency-stiffness plot where new multimodal periodic mo-tions can be generated and the corresponding branch may emanate from the main branch. In the previous strategy, it is implicitly assumed that the nonlinear frequency always increases with the contact coefficient. However, if one is interested in predicting the internal resonances, sequentially increasing the contact stiffness cannot locate the turning points and trace the branch after that. A proper treatment of this issue can be achieved through the application of the pseudo-arclength continuation method in a similar fashion as is proposed in [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF] for the smooth non-linear systems. The details of this continuation algorithm for the present non-smooth systems are provided in Appendix B. Nevertheless, we believe that this is the backbone of the frequency-stiffness plot that contains more useful information about the non-smooth dynamics of the structure in contact with the unilateral support. A higher resolution of the stiffness stepsizes can be used only if additional information is required to investigate the occurrence of the internal resonance in the vicinity of a specific contact stiffness.

For the first NNM of the present 3-DOF system, the internal resonance is identified in the vicinity of k n /k x = 0.018 where 4ω u 1 ≈ ω u 3 . Thus, this modal interaction can be specified as modes 1-3 (4:1) where a smooth transition between the first and third mode occurs. A detailed description of this scenario is shown in Fig. 7 which is the close-up view of the frequency-stiffness plot for 0.015 ≤ k n /k x ≤ 0.022. Following the branch path in Fig. 7 and the evolution of the corresponding trajectories in Fig. 8(a) (the time-history) and Fig. 8(b) (the configuration space) reveals that the contribution of the third mode gradually increases during the interaction from point (a) to (d) where the NNM becomes unstable (point (d) is the first turning point which is identified by the use of the proposed method). From points (d) to (f), the NNM is more like the third mode until it converts totally to the third mode at the intersection of the blue dashed line and orange line (between points (f) and (g) in Fig. 7). Then, from point (g) the first mode again starts to contribute (note the phase shift of π between points (f) and (g) where the first mode again starts to contribute to the interactive mode) and at point (h) the NNM becomes stable again (point (h) is the second turning point). From points (h) to (i) the third more gradually diminishes and finally, the first NNM forms again on the plot backbone. It is obvious that the bilinear equation is not able to reproduce this interaction scenario.

Having introduced, developed and explained our mathematical and computational tools, we can now perform the non-smooth modal analysis of particular systems subjected to the unilateral constraints.

Numerical examples

A computer code has been developed for the non-smooth dynamic analysis and also computing the NNMs of structural members subjected to unilateral constraints. In the following section, we proceed by verifying the time integration algorithm for a MDOF system against the standard FEM package of ANSYS [START_REF]ANSYS Academic Research, Release 13.0, Mechanical APDL and Mechanical Applications Theory Reference[END_REF]. The present approach is then utilised to perform non-linear modal analysis for two case studies.

Non-smooth dynamics of a beam constrained by elastic walls

The transient dynamic response of beam-like members constrained by elastic unilateral walls is of interest to many engineering applications, such as petroleum engineering and drilling industries. For the first dynamic example, we considered a cantilever beam-like component that is partially constrained on both sides by the unilateral supports as shown in Fig 9 . An extended ver-sion of this example can be used to study the vibration of a drill-string with unilateral string-borehole interactions. In this example, the material of the structure is linear elastic with the Young's modulus of E = 200 GPa, density of ρ = 7800 kg/m 3 (no gravity is considered here) and Poisson's ratio of ν = 0.3. The beam length is L = 3.2 m and it has a circular cross-section with the radius of R/L = 0.025 (the cross-section correction factor for the Timoshenko theory is considered here as κ s = 0.925). Note that second moment of area for this system can be obtained as I = πR 4 /4 and the rotational moment of inertia in the mass matrix is J yy = (M/12) (h 2 + 3ˇ R 2 ). Two unilateral elastic walls are constraining the structure on both sides with equal elastic properties as shown in Fig. 9. The dimensionless stiffness of the undamped unilateral support is k w L 4 / (EI) = 500 and its length is x f /L = 0.4.

To perform the non-linear dynamic analysis, one may employ the well-known contact algorithms, i.e. penalty approach (or Lagrange multiplier method in case of contact impenetrability condition) [START_REF]ANSYS Academic Research, Release 13.0, Mechanical APDL and Mechanical Applications Theory Reference[END_REF], and construct a full-scale contact model at the interface of the wall and structure. Consequently, the unilateral constraint mechanism can be captured by using the non-linear iterative solvers (e.g. Newton-Raphson). However, it should be noted that the inherent stability and convergence issues of these solution methods are major concerns for a comprehensive investigation and also the design optimisation.

In the present example, we only use the full-scale FEM analysis to verify the results obtained from the proposed unconditional stable algorithm. The Euler-Bernoulli (BEAM3) and Timoshenko (BEAM188) elements are used in ANSYS [START_REF]ANSYS Academic Research, Release 13.0, Mechanical APDL and Mechanical Applications Theory Reference[END_REF] to construct the FEM mesh. The zero-length contact elements (COMBIN39) are also utilised to specify the non-smooth constraints for the beam length that is constrained by the elastic walls. The unilateral constitutive law is introduced for each element of the walls by using the non-linear generalised capability of the contact elements and considering the gap of δ/R = 0.0125 between the walls and the structure. The transient dynamic response of the system subjected to the harmonic excitation of F (t) = F 0 sin (ωt) at x p /L = 0.5 and with F 0 = 4000 kN and ωL 2 ρA/ (EI) = 3.5 are displayed in Fig. 10. Excellent agreement between the present DLM and the reference FEM results can be observed for both Euler-Bernoulli and Timoshenko beam models. Moreover, no discrepancy exists between the two beam theories until tω/ (2π) = 2. However, it is important to note that the stability of Ansys non-linear solver strongly depends on the convergence values that are set as input parameters, while the proposed method is unconditionally stable. Also, in terms of computation time, the present algorithm is approximately one-third of the Ansys solver for this MDOF example.

To gain further insight into the role of the rotary inertia and shear deformation in the system non-smooth dynamics, we use the present DLM to compare the dynamic response of the Euler-Bernoulli and Timoshenko beams over a longer time span. Since L/R = 40, one may initially expect that the consistency between the two beam models in Fig. 10 continues after tω/ (2π) = 2.

In practice however, it is interesting to notice that as the slender structure progressively interacts with the walls over time, the roles of rotary inertia and shear deformation may become increasingly stronger leading to a large discrepancy between two beam models for a same time span. In Fig. 11, we monitor the dynamic behaviour of the two beam models for the same system as the previous example, but with three different gap lengths and over an extended time period, i.e. 0 ≤ tω/ (2π) ≤ 4. Focusing first on the case of the zero gap (δ = 0), the elastic walls can be interpreted as the smooth bilateral foundation where its reaction is identical in the compression and tension. Therefore, two beam theories lead to consistent dynamic responses, as shown in Fig. 11( When the gap takes nonzero values (δ = 0.0125R in Fig. 11(a.2) and δ = 0.025R in Fig. 11(a.3)), the role of non-smooth unilateral contact interaction becomes important. As the beam-like member makes more contacts with the elastic unilateral supports over time, the error due to the Euler-Bernoulli assumption continuously accumulates. As a result, the Euler-Bernoulli model can lose its accuracy to reproduce the beam's dynamic displacement at the end of the time period. To further quantify this process, the resulting beam's configuration from the Euler-Bernoulli model can be compared with the Timoshenko theory at t = 0.8 and t = 4 in Figs. 11(b.2) and 11(b.3). It is obvious that this discrepancy between the two beam theories depends on the gap length and also, elastic properties of the constraints.

Modal analysis of a beam on the unilateral substrate

To further illustrate the versatility of the present approach, the NNMs of two beam-like structures with the pinned-pinned and fixed-free boundary conditions are computed by using the proposed piecewise smooth model. These beams are interacting with the unilateral elastic substrate throughout their whole lengths while there is no gap at their interface (see Fig. 1(b) with F (t) = 0, δ = 0, c w = 0 and x f = L). Therefore, their NNMs are energy independent and the results of non-linear modal analysis can be presented in the frequency-stiffness plots. The DLM in this example has been constructed using the identical discrete unit cells with the length of h = L/100, resulting in a total of 202 DOFs, before applying the boundary conditions (i.e. n = 101, ñp = 2 and N = 202).

In Fig. 12, we present the first three NNMs of the pinned-pinned model and the corresponding maximum amplitudes w (x, t = 0) at the same level of energy for three scenarios, i.e. the Euler-Bernoulli theory (results are independent of its L/b) and Timoshenko formulation with L/b = 10 and L/b = 20.

Since the fundamental modes are of primary importance, only the dimensionless non-linear frequencies ω u i L 2 ρA/ (EI) on the plot backbone are presented with respect to the normalised contact stiffness kw = k w L 4 / (EI). The internal resonance branches that may emanate from the main branch during the continuation algorithm are not shown here. Each branch of non-smooth modes is initiated on the vertical axis at the linear mode of the corresponding beam with no supporting foundation ( kw = 0). Note that these dimensionless frequency-stiffness plots enable us to obtain the resonance frequencies of the Euler-Bernoulli beams for 0 ≤ kw ≤ 1000 with any geometrical and material properties and in a similar fashion, for the Timoshenko beams with L/b = 10 and L/b = 20. These plots can be utilised for the practical engineering applications, e.g. computing the resonance velocities of the moving loads applied on the beams or calculating the non-smooth critical velocities of the fluidconveying pipes which are resting on the tensionless subgrade. In parallel to the numerical results, the system bilinear frequencies are obtained from Eq. 38 and plotted in Figs. 12(a.1)-(a.3) for the three scenarios.

We observe an excellent agreement between the first NNM frequency and bilinear approximation in Fig. 12(a.1), although the system is not commutative (i.e. K -• K + = K + • K -). This can be explained by the shape of the beam maximum amplitudes for the first NNM. Fig. 12(b.1) demonstrates that the periodic motion of the first NNM always preserves its bilinear (smooth) nature over the whole period, in the sense that no non-smooth effects are experienced during the free motion. Subsequently, increasing the substrate stiffness does not change its initial form which is obtained at kw = 0 for both Timoshenko and Euler-Bernoulli models. In other words, over the first period T u 1 = 2π/ω u 1 , all the DLM nodes vibrate synchronously. They start to move from their maximum amplitudes at t = 0, as depicted in Fig. 12(b.1), and simultaneously reach the contact interface line, i.e. z = 0. At this time instant, the whole energy of the system, which was purely elastic at t = 0, is transformed to purely kinetic. Subsequently, the nodes start to interact with the substrate and reach their maximum amplitudes on the other side of the interface at t = T u 1 /2 (note that the nodal maximum amplitudes in contact with the substrate can be different from what is shown in 12(b.1) and it depends on the value of kw ). From t = T u 1 /2 to t = T u 1 , the lattice nodes travel back through the same path as the first half of the period to reach their initial positions. Therefore, the system dynamics associated with the time span that the beam has no contact with the substrate is decoupled from the interacting span. As a consequence, the bilinear formula can accurately predict the first non-linear frequency, regardless of the commutativity condition of the stiffness matrix.

The inaccurate approximation of the bilinear formula for the second and third NNMs in Figs. 12(a.2) and 12(a.3) can now be understood in light of the computed maximum amplitudes. If we increase kw , the maximum amplitudes of these modes do not keep their initial forms at kw = 0. Instead, Figs. 12(b.2) and 12(b.3) depict how the shape of the maximum amplitudes change in accordance with the contact coefficient. Unlike the first mode, the non-linear periodic motions do not feature bilinear behaviour. As a result, the system experiences the non-smooth effects and the bilinear approximation loses its accuracy by increasing the unilateral substrate stiffness.

Based on the smooth character of the first mode in the modal analysis, we can utilise the bilinear formula to propose an analytical expression for the first NNM frequency of the beam-substrate system with the pinned-pinned boundary conditions. Using the separation of variables approach, it is straightforward to obtain the eigenvalues of the corresponding linear system. Thus, after reducing Eqs. 5 and 6 to the Euler-Bernoulli theory and assuming F (x, t) = 0, c w = 0, δ = 0 and H (w) = 1 to obtain the free vibration equation, the natural frequencies of the Euler-Bernoulli beam on the bilateral foundation can be obtained as

ω i = i 4 π 4 EI ρAL 4 + k w ρA ( 39
)
Using Eq. 39 for the first mode (i = 1), we can compute the natural frequencies of the linearised system with and without the foundation. Substituting them into Eq. 38, the first NNM frequency is obtained as

ω u 1 = 2 π 4 EI/L 4 + k w √ ρA 1 + 1 + k w L 4 / (π 4 EI) (40)
A similar equation can also be derived for the first NNM frequency in the analogous problem of rectangular Kirchhoff plate with the simply supported edges. The linear natural frequency is first obtained from Eq. 7 for the plate on the bilateral foundation. Subsequently, the first eigenvalues of the system in the presence of the foundation and with no foundation are used in Eq.

38 to achieve the following equation for the first non-linear frequency of the plate-substrate system

ω u 11 = π 4 Eb 2 3ρ(1-ν 2 ) 1 c 2 + 1 L 2 12kw(1-ν 2 ) π 4 Eb 3 + 1 c 2 + 1 L 2 2 1 c 2 + 1 L 2 + 12kw(1-ν 2 ) π 4 Eb 3 + 1 c 2 + 1 L 2 2 (41)
For the case of beam with the fixed-free boundary conditions, the NNM frequencies on the plot backbone exhibit qualitatively similar behaviour, as shown in Fig. 13. Once again, the accuracy of the bilinear equation to predict the first frequency is remarkable (see Fig. 13(a.1)), as the system does not experience any non-smooth effects. Thus, the shape of the maximum amplitudes for the first mode does not vary with the substrate stiffness for 0 ≤ kw ≤ 1000 (see Fig. 13(b.1)). However, when kw is increased for the higher modes, the initial shape of the maximum amplitudes evolves smoothly into a new form with a smaller contact area between the beam and substrate, as shown in Figs. 13(b.2) and 13(b.3). Our results confirm that the error due to the bilinear approximation for the second and third NNMs becomes more significant at higher kw , as indicated in Figs. 13(a.2) and 13(a.3).

Granular chain with the non-smooth contacts

The last application example is inspired from the complex non-smooth behaviour of granular materials. Two main reasons can be considered as the source of non-linearity in the granular medium: the geometrical effects of the interacting particles, which can be modelled by the Hertzian contact theory, and the repulsive (tensionless) interactions between the grains, which cause the non-smoothness in the system [START_REF] Jayaprakash | Nonlinear normal modes and band zones in granular chains with no precompression[END_REF]. The former makes the non-linear vibrational modes energy-dependent, while assuming the latter as the only source of non-linearity keeps the system dynamics energy-invariant.

In this example, we consider the basic but important example of 1D granular chain with only the second type of non-linearity, i.e. the non-smooth contact mechanism between the identical particles (see Fig. 14(a)). Note that this discrete system is allowed to exhibit the contact/detachment phenomenon between its components where in the previous examples the constraints are introduced over the boundary of the solid. However, as the constitutive contact law is defined with respect to the relative displacements between the particles, the non-smoothness is still concentrated at the origin (see Fig. 14(b)). Such a system is of significant interest due to its tunability to exhibit from weakly non-linear characteristics (with the static constraining compression) to highly non-linear dynamics (with no constraining load) [START_REF] Starosvetsky | Effective particles and classification of the dynamics of homogeneous granular chains with no precompression[END_REF]. Here, no longitudinal compressive load is assumed to be able to preserve the highest non-smooth conditions. As our objective is to show the applicability of the proposed approach to compute the fundamental NNMs of the granular chain, a detailed bifurcation analysis of the energy-independent periodic motions is not presented in this study. Interested readers may refer to [START_REF] Jayaprakash | Nonlinear normal modes and band zones in granular chains with no precompression[END_REF][START_REF] Starosvetsky | Effective particles and classification of the dynamics of homogeneous granular chains with no precompression[END_REF] for a detailed discussion on the energy-variant NNMs of the same problem, but with the Hertzian contact model.

Following the same strategy as in Section 3 to compute the NNMs in a stepwise manner, we first define the interaction law between the oscillators by the bilinear equation with k c and k t as the contact stiffness in compression and tension, respectively, where 0 ≤ k t ≤ k c (see Fig. 14(b)). Thus, ď = 1k t /k c can be expressed as the damage parameter where ď = 0 (i.e. k t = k c ) corresponds to the linear elastic case (perfect connection), 0 < ď < 1 describes the damaged (imperfect) connection and ď = 1 (i.e. k t = 0) specifies the unilateral contact mechanism. Before proceeding to apply the shooting technique, we ought to modify the second-order equation of motion to express the interactions in terms of the relative displacements between the oscillators.

Adopting the approach that we utilised in Eq. 15, it is convenient to convert the diagonal contact constraint matrix H (t) into two separate sparse matrices for the contact (compressive) state as H c (t) and the detachment (tensional) state as H t (t). The corresponding unforced equations of motion in the matrix notation can be obtained as (no damping is present)

M • D (t) + [k c H c (t) + k t H t (t)] • D (t) = 0 (42)
For the sake of clarity, let us consider the case of n = 3 particles. The sparse matrices of H c (t) and H t (t) can be expressed as the following

H c (t) = H(-u 1 )+ H(u 1 -u 2 ) -H(u 1 -u 2 ) 0 -H(u 1 -u 2 ) H(u 1 -u 2 )+ H(u 2 -u 3 ) -H(u 2 -u 3 ) 0 -H(u 2 -u 3 ) H(u 2 -u 3 )+ H(u 3 ) (43)
H t (t) = 2-H(-u 1 )-H(u 1 -u 2 ) -1+ H(u 1 -u 2 ) 0 -1+ H(u 1 -u 2 ) 2-H(u 1 -u 2 )-H(u 2 -u 3 ) -1+ H(u 2 -u 3 ) 0 -1+ H(u 2 -u 3 ) 2-H(u 2 -u 3 )-H(u 3 ) (44)
The equations of motion for other values of n can be derived by extending Eqs. It can be observed that the present numerical continuation method is able to robustly and efficiently trace the backbone branch and obtain the non-smooth periodic orbits across the range of k t /k c .

Interestingly, we can find a qualitative resemblance between the behaviour of these NNMs and the energy-dependent system with the Hertzian contact interactions, which has been studied in [START_REF] Jayaprakash | Nonlinear normal modes and band zones in granular chains with no precompression[END_REF]. Correspondingly, three types of (fundamental) NNMs can be recognised in Fig. 15. There are non-linear modes that are influenced by the detachment phenomenon between the particles and their neighbours. Thus, the particles become motionless at a distance from their zero positions for a fraction of the whole period, as illustrated in Figs.

15(a.2), (c.2), (f.2) and (g.2). Due to the occurrence of separation, it is obvious that the continuum homogenised theories are not able to capture this type of NNMs. Therefore, it can highlight the importance of lattice dynamics as the role of micro-structure is considered in the response of the granular system.

There also exist periodic orbits where the non-smooth character reveals itself only as the non-synchronous motion (see Figs. The results of Fig. 15 also indicate that the NNMs may become unstable as their Floquet multipliers move towards the outside of the unit circle. This behaviour is analogous to the process seen for the 3-DOF bar model in Fig. 3. Thus, the number of modes in this system exceeds the number of DOFs.

The instability can be triggered by the interaction between two NNMs or the evolution of a new mode with a period equal to the integer multiple of the fundamental period. Only as an illustrative example, the incremental evolution of the new NNM for the chain with n = 2 particles is displayed in Fig. 17.

We already observed that the second NNM of this system manifests itself as a smooth mode (see Fig. 15(b.1) and (b.2) for n = 2). However, when it becomes unstable, a small perturbation of the maximum amplitudes can make the system to form a new non-smooth periodic orbit. This can be clearly observed in Fig. 17 where the second NNM is still stable at ď = 0.7 with the frequency of ωu 2 = 1.379 and the maximum amplitudes of {-1.369, 1.369}. After the onset of instability, the new NNM starts to grow incrementally over a new branch with a change of T u 2 → 27T u 2 for the period. Three representative points on this new branch are depicted in Fig. 17 where the non-linear frequencies change from ωu 2 = 0.051 (at ď = 0.704, right after the bifurcation) to ωu 2 = 0.043 (at ď = 1 which is the chain with the unilateral contact). The maximum amplitudes are also indicated as the insets in Fig. 17(b).

Conclusions

We have presented a numerical framework to study the non-smooth dynamics of mechanical systems subjected to unilateral elastic constraints. In our model, different stiffnesses of the constitutive law at two sides of the contact interface is the only source of non-linearity, under the assumption that the material properties are linear and the system remains in the range of small deformation theory. The discrete lattice method is used to obtain the space-discretised model of the problem. Taking advantage of the piecewise smooth nature of this model, we developed an unconditionally stable inte-gration scheme to solve the corresponding second-order equations of motion.

The proposed method is based on tracking the trajectory of lattice nodes to update the active constraints at each time step, without implementation of any additional contact strategy.

The main contribution of this research lies in the combination of proposed unconditionally stable time integration method with the shooting technique to compute the NNMs of the system. Accordingly, the variation of these nonlinear periodic orbits with respect to the contact stiffness has been presented in the frequency-stiffness plots, which are analogous to the frequency-energy plots in the smooth non-linear systems.

Important observations and conclusions can be made from the non-smooth modal analysis of the representative examples. The stability analysis of the computed modes on the backbone of the frequency-stiffness plot reveals that the NNMs may become unstable, in the sense that a small disturbance of the initial conditions leads to the elimination of the corresponding mode. This oscillatory instability can be triggered due to the modal interaction between the NNMs with the higher modes (the internal resonance phenomenon) or the bifurcation of new non-smooth periodic motions with a period equal to the integer multiple of the fundamental period. The former can be detected once the frequencies of two modes are commensurate by utilising the numerical continuation algorithm. The latter can also be computed by using the perturbation technique. In addition, we observed that the bilinear equation can predict the system resonance frequencies only if the associated modes manifest themselves as smooth modes (e.g. the first mode of the beam-substrate system), regardless of the commutativity conditions of the stiffness matrix. Otherwise, the bilinear formula loses its accuracy. This can illustrate the importance of the proposed numerical algorithm, as the condition of smooth character for the non-linear modes may not be satisfied for many important cases of practical interest. Moreover, the computed NNMs of the tensionless granular chain exhibited qualitatively similar characteristics as the energy-dependent modes in the same model but with the Hertzian contact between the particles. In particular, the 1D tensionless case study has significant implications for modeling the dynamics of granular materials by the discrete lattice methods as the homogenised continuum theories are not able to capture the NNMs associated with the detachment phenomenon between the particles. Future work on this subject could include the modal analysis of mechanical systems in the presence of the gap, which makes the NNMs to vary with the total energy. Moreover, it would be desirable to extend the proposed algorithm to treat non-smooth problems when the source of discontinuity is in the velocity field, e.g. the non-smooth stick-slip vibrations induced by the dry friction or the bit-rock interaction in the horizontal drilling problem [START_REF] Ritto | Drill-string horizontal dynamics with uncertainty on the frictional force[END_REF]. It is also conceivable that the principles of the present method can be utilised for the non-linear buckling analysis of a beam/rod inside a horizontal/cylindrical unilateral constraint.

Acknowledgment

The first author gratefully acknowledges the scholarship support from The University of Western Australia (SIRF and UPIAS) and the Australian Government towards a successful completion of the present work. The first author also would like to thank Professor A. Dyskin and Professor E. Pasternak (both UWA) for insightful discussions before starting this work. to zero. Therefore, the following three conditions have to be verified for the actual solution of the NNMs

       S(T k + ∆T k , z (0) k + ∆z (0) k ) = 0 Ḋ(0) k (z (0) k + ∆z (0) k ) = 0 E(z (0) k + ∆z (0) k ) -ě = 0 (A.3)
where ∆T k ∈ R and ∆z (0) k ∈ R 2N are the corrections at each iteration. Eq. A.3 sets the necessary and sufficient conditions of the present constrained minimisation problem (i.e. minimisation of the shooting function with constraints on the phase condition and system energy). The corrections can be obtained from the Taylor expansion of Eq. A.3, neglecting about the higher order terms and constructing the system of equations with ∆T k and ∆z (0) k as the unknown parameters. Thus, the corrections at each iteration are computed by utilising the Moore-Penrose matrix inverse [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF] to solve the following overdetermined system

  J 11 J 12 J 21 0 J 31 0   • ∆z (0) k ∆T k = -    S(T k ,z (0) k ) Ḋ(0) k E(z (0) k )-ě    (A.4)
where J 11 ∈ R 2N ×2N , J 12 ∈ R 2N , J 21 ∈ R N ×2N and J 31 ∈ R 1×2N are the components of the Jacobian matrix (0 ∈ R N in Eq. A.4 is the zero vector) and they can be identified from the Taylor expansion of Eq. A.3 as

J 11 = ∂S ∂z (0) | (T k ,z (0) k ) = ∂z(t, z (0) ) ∂z (0) | (T k ,z (0) k ) -I 2N (A.5) J 12 = ∂S ∂T | (T k ,z (0) k ) = ż(T k , z (0)
k ) (A.6)

J 21 = ∂ Ḋ(0) ∂z (0) = [0 I] (A.7) J 31 = ∂E ∂z (0) | (z (0) k ) = z (0) k T • B (A.8)
with I 2N ∈ R 2N ×2N and I ∈ R N ×N as the identity matrices and 0 ∈ R N ×N as the zero matrix. Thus, the components J 12 , J 21 and J 31 can be computed from Eqs. A.6-A.8. We only need to evaluate ∂z(t,z (0) )

∂z (0)
∈ R 2N ×2N in Eq. A.5 to obtain J 11 and complete the Jacobian matrix. Considering g = ∂z(t,z (0) )

∂z (0)
, the extended form of g is expressed as

g =                         ∂D 0 ∂D (0) 0 ∂D 0 ∂D (0) 1 ••• ∂D 0 ∂D (0) N ∂D 0 ∂ Ḋ(0) 0 ∂D 0 ∂ Ḋ(0) 1 ••• ∂D 0 ∂ Ḋ(0) N ∂D 1 ∂D (0) 0 ∂D 1 ∂D (0) 1 ••• ∂D 1 ∂D (0) N ∂D 1 ∂ Ḋ(0) 0 ∂D 1 ∂ Ḋ(0) 1 ••• ∂D 1 ∂ Ḋ(0) N . . . . . . . . . . . . . . . . . . ∂D N ∂D (0) 0 ∂D N ∂D (0) 1 ••• ∂D N ∂D (0) N ∂D N ∂ Ḋ(0) 0 ∂D N ∂ Ḋ(0) 1 ••• ∂D N ∂ Ḋ(0) N ∂ Ḋ0 ∂D (0) 0 ∂ Ḋ0 ∂D (0) 1 ••• ∂ Ḋ0 ∂D (0) N ∂ Ḋ0 ∂ Ḋ(0) 0 ∂ Ḋ0 ∂ Ḋ(0) 1 ••• ∂ Ḋ0 ∂ Ḋ(0) N ∂ Ḋ1 ∂D (0) 0 ∂ Ḋ1 ∂D (0) 1 ••• ∂ Ḋ1 ∂D (0) N ∂ Ḋ1 ∂ Ḋ(0) 0 ∂ Ḋ1 ∂ Ḋ(0) 1 ••• ∂ Ḋ1 ∂ Ḋ(0) N . . . . . . . . . . . . . . . . . . ∂ ḊN ∂D (0) 0 ∂ ḊN ∂D (0) 1 ••• ∂ ḊN ∂D (0) N ∂ ḊN ∂ Ḋ(0) 0 ∂ ḊN ∂ Ḋ(0) 1 ••• ∂ ḊN ∂ Ḋ(0) N                         (A.9)
According to the numerical finite-difference method [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF], computing g requires to repeat the time integration of the MDOF system for 2N times where the components of z

k are successively perturbed to numerically capture the members of g. Obviously, this procedure needs far too many computations to generate the Jacobian matrix at every iteration. To circumvent this issue, we follow the approach presented by [START_REF] Zuo | Non-linear real and complex modes of conewise linear systems[END_REF] and differentiate Eq. 35 with respect to z (0) to obtain ġ = Q • g as the governing differential equation for g. Starting with the initial conditions of g 0 = ∂z(0,z (0) ) ∂z (0) = I 2N , this ordinary differential equation can be solved concurrently with the solution of the piecewise linear system and the solution is stored at each step of the time integration procedure to finally obtain g at t = T . For example, using the trapezoidal rule (i.e. α m = α f = 0, γ = 1/2, β = 1/4), the numerical solution of g at the end of each time step is obtained as

g r+1 = I 2N - ∆t 2 Q -1 • I 2N + ∆t 2 Q • g r (A.10)
Therefore, the Jacobian matrix is formed at each iteration of the shooting scheme and the corrections are calculated from Eq. A.4. The procedure exits the iterative loop and returns z

k and T k as the NNM and its period when vector to the coefficient matrix in Eq. A.4 as

  J 11 J 12 J 21 0 J 31 0   • p z,j p T,j = 0 3N +1 (B.2)
with 0 3N +1 ∈ R 3N +1 as the zero vector. Thus, p j can be computed by fixing one of its components, solving Eq. B.2 by the Moore-Penrose matrix inverse, and normalising the resulting vector to the unit length. The stiffness of the system at point j +1 can also be specified as k w,j+1 = k w,j +∆k w,j where ∆k w,j ∈ R is the contact stiffness step-size and it should be carefully controlled in accordance with s j . In the corrector step, the initial guess obtained from Eq.

B.1 and contact stiffness of k w,j+1 are used in the shooting iterative process to find the periodic orbit at point j + 1 by solving Eq. A.3. There is only an extra condition which forces the corrections to converge to the solution in the direction orthogonal to the tangent vector obtained in the predictor step.

Thus, the corrections can be computed from the modified form of Eq. A.4 as

     J 11 J 12 J 21 0 J 31 0 p z,j p T,j      • ∆z (0) k,j+1 ∆T k,j+1 = -        S(T k,j+1 ,z (0) k,j+1 ) Ḋ(0) k,j+1 E(z (0) k,j+1 )-ě 0        (B.3)
herein, the superscript indicates the transpose operator and the components of the Jacobian matrix are computed by the use of Eqs. A.5-A.8 at point j + 1. The last equation of this overdetermined system satisfies the orthogonality condition by prescribing the inner product of the tangent vector and the corrections vector to be zero.

In order to properly trace the branch in the same direction during the present numerical continuation scheme, in addition to the constraints set by Eqs. B.2 and B.3, we also need to change the sign of ∆k w,j at the turning points.

The turning points can be identified with two important properties. First, the tangent to the branch with respect to the period, i.e. p T,j , changes its sign at these points, i.e. p - T,j • p + T,j < 0 where p - T,j and p + T,j refer to the value of p T,j just before and after the turning points (it is obvious that p T,j = 0 at the turning points). Second, once the branch reaches the turning points, following the continuation process with the same sign for ∆k w,j results in failure of the shooting algorithm to converge. Therefore, the path following method requires a decision module to track both sign(p T,j ) and the solution convergence. This module can detect the occurrence of turning points and accordingly, handles sign(∆k w,j ) with an adoptive control over it. It should be noted that in the areas close to the turning points, a higher resolution is used to capture how the tongues emanate from the branch, whereas larger steps can be used once we are far from the turning points.

Frequency-stiffness plots for the NNMs of the 3-DOF chain of oscillators with δ = 0 (see Fig. 2). The vertical axis shows the dimensionless frequencies ω u i m/k x where the superscript u denotes the unilateral interaction and i = 1, 2, 3. The horizontal axis shows the range of normalised contact stiffness k n /k x . The unstable NNM frequencies are shown with the dashed lines. The orange line also depicts the bilinear approximation for the NNM frequencies obtained from Eq. ). (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The vertical axis shows the displacement components of the 3-DOF chain of oscillators with δ = 0 (see Fig. 2). The horizontal axis depicts the dimensionless time t/T u 1 where T u 1 = 2π/ω u 1 is the period. (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The vertical axis shows the displacements of the first and second mass, i.e. u 1 (t) and u 2 (t) respectively. The horizontal axis depicts the displacement of the third mass ). (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The vertical axis shows the displacement components of the 3-DOF chain of oscillators with δ = 0 (see Fig. 2). The horizontal axis depicts the dimensionless time t/T u 1 where T u 1 = 2π/ω u 1 is the period. (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The vertical axis shows the displacements of the first and second mass, i.e. u 1 (t) and u 2 (t) respectively. The horizontal axis depicts the displacement of the third mass u 3 (t). The gray coloured plots show the unstable NNMs. 69 7 Close-up view of the modes 1-3 (4:1) internal resonance that occurs for the first NNM of the 3-DOF chain of oscillators with δ = 0 (see Fig. 2 and the corresponding frequency-stiffness plots in Fig. 3). The first and third modes start to interact once 4ω u 1 ≈ ω u The cantilever beam is constrained by the unilateral elastic walls on both sides with the gap of δ. The structure has a circular cross-section with the radius of R and it is subjected to the harmonic load of F (t) applied at point x p , i.e.

F (x, t) = F (t) δ (xx p ). 72

10 Transient dynamic response of the beam with the circular cross-section constrained by the elastic unilateral walls on both sides (see Fig. 9) with k w L 4 / (EI) = 500, c w = 0,

x f /L = 0.

L

x p F(t)

d i i+1 i+2 i-1 i-2 h (a) (b) x Elastic unilateral support (kw,cw,d) z L x f x p x F(t) d i i+1 i+2 i-1 i-2 KN CN h h x z c b L d x p y p x f y f x y hx (i,j) (i-1,j) (i-2,j) (i+1,j) (i+2,j) (i,j+2) (i-1,j+2) (i-2,j+2) (i+1,j+2) (i+2,j+2) (i,j+1) (i-1,j+1) (i-2,j+1) (i+1,j+1) (i+2,j+1) (i,j-1) (i-1,j-1) (i-2,j-1) (i+1,j-1) (i+2,j-1) (i,j-2) (i-1,j-2) (i-2,j-2) (i+1,j-2) (i+2,j-2)
hy F(t)

E la st ic u n il a te ra l su p p o rt (k The DOFs of each lattice node is identified as {w i,j (t) , ϕ i,j (t) , θ i,j (t)}. The unilateral support is also modelled as a network of parallel one-way spring-damper elements (K N = k w h x h y , C N = c w h x h y where h x and h y are the length of the unit cells in x and y directions, respectively). The incremental evolution of the NNMs are presented by the time-history of the free system over one period T u 1 = 2π/ω u 1 . The vertical axis shows the displacement components of the 3-DOF chain of oscillators with δ = 0 (see Fig. 2). The horizontal axis depicts the dimensionless time t/T u 1 where T u 1 = 2π/ω u 1 is the period. (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The vertical axis shows the displacements of the first and second mass, i.e. u 1 (t) and u 2 (t) respectively. The horizontal axis depicts the displacement of the third mass u 3 (t). The gray coloured plots show the unstable NNMs. where the new periodic motion is born at point (d) with a period equal to the double of the fundamental period (T u 1 → 2T u 1 ). (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The vertical axis shows the displacement components of the 3-DOF chain of oscillators with δ = 0 (see Fig. 2). The horizontal axis depicts the dimensionless time t/T u 1 where T u 1 = 2π/ω u 1 is the period. (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The vertical axis shows the displacements of the first and second mass, i.e. u 1 (t) and u 2 (t) respectively. The horizontal axis depicts the displacement of the third mass u 3 (t). The gray coloured plots show the unstable NNMs. ). (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The vertical axis shows the displacement components of the 3-DOF chain of oscillators with δ = 0 (see Fig. 2). The horizontal axis depicts the dimensionless time t/T u 1 where T u 1 = 2π/ω u 1 is the period. (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The vertical axis shows the displacements of the first and second mass, i.e. u 1 (t) and u 2 (t) respectively. The horizontal axis depicts the displacement of the third mass u 3 (t). The gray coloured plots show the unstable NNMs.

w, c w, d ) (c) (d)
k n /k x ω u 1 • m/k x (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)
(h) (d) u 2 u 1 u 3 (t) u i (t) u i (t) (l)
(d.3) u 2 u 1 u 3 (t) u i (t) u i (t)
(i.3) u 2 u 1 u 3 (t) u i (t) u i (t)

  a.1) by the time-history of the displacements at x = L/2 and x = L and in Fig. 11(b.1) by the deformed shape of the structure at two time instants t = 0.8 and t = 4 where t = tω/ (2π).

43 and 42

 42 for different number of particles in a similar manner. Now, using the linear normal modes and the corresponding frequencies at ď = 1k t /k c = 0 as the initial modal information, we can employ the combination of present time-integration technique and shooting method to compute the fundamental NNMs as a function of ď when k t changes from k t = k c to k t = 0. The results of each step are used as the initial guess for the following step and the nonsmooth modal results are presented over the backbone branch, similar to the previous examples. The fundamental NNMs for the chain with n = 2, n = 3 and n = 4 particles are depicted in Figs. 15(a.1)-(i.1). The maximum amplitudes of the non-linear modes at ď = 0, ď = 0.5 and ď = 1 are also presented in these figures for a fixed level of energy. Moreover, the time-history of the fundamental modes with k t = 0 (i.e. the unilateral constraint) are depicted in Figs. 15(a.2)-(i.2).

  [START_REF] Vakakis | Non-linear normal modes and their applications in vibration theory: An overview[END_REF](e.2), (h.2) and (i.2)). Furthermore, depending on the number of particles, we may observe modes with the smooth character over the whole period, as indicated in Fig.15(b.2) and (d.2) for n = 2 and n = 3, respectively. In analogy with the smooth nature of the first NNM in the beam-substrate system (see Section 4.2), the bilinear formula is able to accurately compute the non-linear frequencies for these modes too.The excitability of the computed fundamental NNMs due to the harmonic load is studied in Fig.16for four representative examples with n = 2, n = 4and n = 6. The transient response of the chain is driven into resonance as the excitation frequency reaches the NNM frequency. We can also observe that the forced response of the system occurs in the vicinity of the corresponding NNM (for example, compare Fig.16(a) with Fig.15(a.2)).

  38. (a) The first NNM is presented where the backbone branch starts from ω u 1 m/k x = 0.445 on the vertical axis. More details about the first NNM are presented here. The incremental evolution of the NNMs for the representative points (a)-(l) on its backbone is illustrated in Fig. 4. New periodic motions are also born at points (d) and (i) where the NNMs on the backbone branch become unstable. The corresponding incremental evolution of these newly born NNMs for representative points (d.1)-(d.11) and (i.1)-(i.11) are depicted in Figs. 5 and 6, respectively. (b) The second NNM is presented where its backbone branch starts from ω u 2 m/k x = 1.247 on the vertical axis. (c) The third NNM is presented where its backbone branch starts from ω u 3 m/k x = 1.802 on the vertical axis. The start-points on the vertical axis are the three natural frequencies of the 3-DOF system with k n = 0 which can be obtained from the linear eigenvalue analysis of Eq. 37. 66 4 Trajectories of NNMs for points (a)-(l) on the backbone of the first NNM in Fig 3(a). (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period T u 1 = 2π/ω u 1 . The vertical axis shows the displacement components of the 3-DOF chain of oscillators with δ = 0 (see Fig. 2). The horizontal axis depicts the dimensionless time t/T u 1 where T u 1 = 2π/ω u 1 is the period. (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The vertical axis shows the displacements of the first and second mass, i.e. u 1 (t) and u 2 (t) respectively. The horizontal axis depicts the displacement of the third mass u 3 (t). The gray coloured plots show the unstable NNMs. 67 58 Trajectories of NNMs for points (d)-(d.11) of the first NNM in Fig 3(a) where the new periodic motion is born at point (d) with a period equal to the double of the fundamental period (T u 1 → 2T u 1

u 3

 3 (t). The gray coloured plots show the unstable NNMs. 68 6 Trajectories of NNMs for points (i)-(i.11) of the first NNM in Fig 3(a) where the new periodic motion is born at point (i) with a period equal to quintuple of the fundamental period (T u 1 → 5T u 1

3 .

 3 The vertical axis shows the first dimensionless frequency ω u 1 m/k x . The horizontal axis shows the range of normalised contact stiffness k n /k x . The unstable NNM frequencies are shown with the dashed line. The orange line also depicts the (ω u 3 /4) m/k x obtained from the shooting method. 70 59 Trajectories of NNMs for points (a)-(l) of Fig 7 where the first and third modes start to interact once 4ω u 1 ≈ ω u 3 and accordingly, the modes 1-3 (4:1) internal resonance occurs. (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The displacement components of the 3-DOF chain of oscillators (see Fig. 2 with δ = 0) are shown with u 1 (t), u 2 (t) and u 3 (t). (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The gray coloured plots show the unstable NNMs. 71 9

  [START_REF] Attar | Dynamic response of cracked timoshenko beams on elastic foundations under moving har-51 monic loads[END_REF] and δ/R = 0.0125. The structure is subjected to the harmonic excitation of F (t) = F 0 sin (ωt) at x p /L = 0.5 with F 0 = 4000 kN and ωL 2 ρA/ (EI) = 3.5. The vertical axis shows the time-history of the dimensionless displacement w = w/d s in z-direction at x = L/2 (thin solid line for the present discrete model and thin dashed line for the ANSYS results) and x = L (thick solid line for the present discrete model and thick dashed line for the ANSYS results) where d s = F 0 L 3 / (3EI). The horizontal axis depicts the dimensionless time t = tω/ (2π) which is normalised with respect to the period of the applied harmonic load. (a) Present discrete lattice model (DLM) versus the FEM results (ANSYS) for the Euler-Bernoulli theory (ǧ = 0). (b) DLM versus the FEM results for the Timoshenko theory (ǧ = 1). 73 60 11 Transient dynamic response of the beam with the circular cross-section constrained by the elastic unilateral walls on both sides (see Fig. 9) with k w L 4 / (EI) = 500, c w = 0, x f /L = 0.4 and three different gap lengths. The structure is subjected to the harmonic excitation of F (t) = F 0 sin (ωt) at x p /L = 0.5 with F 0 = 4000 kN and ωL 2 ρA/ (EI) = 3.5. The results obtained from the Timoshenko beam model (ǧ = 1) are compared with the Euler-Bernoulli theory (ǧ = 0). (a) The vertical axis shows the time-history of the dimensionless displacement w = w/d s in z-direction at x = L/2 and x = L where d s = F 0 L 3 / (3EI). The horizontal axis depicts the dimensionless time t = tω/ (2π) which is normalised with respect to the period of the applied harmonic load. (b) The vertical axis shows the dimensionless lateral displacement throughout the beam length (i.e. the deformed shape of the structure) at two different time instants t = 0.8 and t = 4. The horizontal axis depicts the dimensionless axial coordinate x/L. Note that the gap length for Figs. (a.1) and (b.1) is δ = 0, Figs. (a.2) and (b.2) is δ = 0.0125R, and Figs. (a.3) and (b.3) is δ = 0.025R. 74 12 (a.1)-(a.3) Frequency-stiffness plots for the first three NNMs of the pinned-pinned beam fully supported by the unilateral substrate (see Fig. 1(b) with the pinned-pinned boundary conditions and F (t) = 0, δ = 0, c w = 0 and x f = L). The vertical axis shows the dimensionless frequency ω u i L 2 ρA/ (EI) where the superscript u denotes the unilateral interaction. The horizontal axis shows the range of dimensionless contact stiffness kw = k w L 4 / (EI). Only the fundamental non-linear frequencies are shown here and the internal resonance branches that may emanate from the plot backbone are not included. The dot-dashed lines also depict the bilinear approximation for the NNM frequencies obtained from Eq. 38. (b.1)-(b.3) The evolution of the maximum amplitudes as a function of kw for the first three NNMs. Theses maximum amplitudes are shown for the Timoshenko beam with L/b = 20. Two other scenarios (Euler-Bernoulli beam and Timoshenko beam with L/b = 10) have qualitatively similar shapes. Note that the maximum amplitudes of the first NNM do not change with kw . 75 61 13 (a.1)-(a.3) Frequency-stiffness plots for the first three NNMs of the fixed-free beam fully supported by the unilateral substrate (see Fig. 1 with F (t) = 0, δ = 0, c w = 0 and x f = L). The vertical axis shows the dimensionless frequency ω u i L 2 ρA/ (EI). The horizontal axis shows the range of dimensionless contact stiffness kw = k w L 4 / (EI). Only the main frequencies on the plot backbone are shown here and the internal resonance branches that may emanate from the backbones are not included. The dot-dashed lines also depict the bilinear approximation for the NNM frequencies obtained from Eq. 38. (b.1)-(b.3) The evolution of the maximum amplitudes as a function of kw for the first three NNMs. Theses maximum amplitudes are shown for the Timoshenko beam with L/b = 20. Two other scenarios (Euler-Bernoulli beam and Timoshenko beam with L/b = 10) have qualitatively similar shapes. Note that the maximum amplitudes of the first NNM do not change with kw . 76 14 (a) Granular chain of identical particles with non-smooth constraint between them. (b) The bilinear constitutive model for the restoring contact force between the particles (oscillators); k t = k c represents the linear elastic model (perfect connection) and k t = 0 corresponds to the unilateral constraint (fully damaged connection). Accordingly, ď = 1k t /k c can be considered as the damage parameter which can vary from ď = 0 (perfect or intact connection between particles) to ď = 1 (tensionless contact constraint). 77 15 (a.1)-(i.1) The frequency plots of the granular chain with n = 2, n = 3 and n = 4 particles (see Fig. 14). The dimensionless frequency ωu i = ω u i m/k c is shown as a function of the damage parameter ď = 1k t /k c . The maximum amplitudes of the NNMs at ď = 0, ď = 0.5 and ď = 1 are also presented. (a.2)-(i.2) The time-history of the fundamental modes with k t = 0 (i.e. ď = 1). The displacement of particles are indicated versus dimensionless time t/T u i where T u i = 2π/ω u i . 78 62 (k t = 0) due to the harmonic excitation. The vertical axis shows the dimensionless displacement of the particles ūi (t) = u i (t) /u s where u s = 10F 0 /k c . The horizontal axis represents dimensionless time tω/ (2π) where ω is the frequency of the applied force(s). The inset figures show the number and the direction of the applied harmonic loads. (a) The number of particles is n = 2 and the frequency of the applied load is ω = ωu 1 k c /m = 0.739 k c /m. (b) n = 4 and ω = ωu 2 k c /m = 0.829 k c /m. (c) n = 4 and ω = ωu 4 k c /m = 1.281 k c /m. (d) n = 6 and ω = ωu 6 k c /m = 1.323 k c /m. 79 17 Trajectories of NNMs for the second mode of the granular chain with n = 2 particles (see Fig. 14). The second mode becomes unstable at ď = 0.704 (see Fig. 15(b.1) where the dashed line shows the unstable modes). As a consequence, the new periodic motion is born with a change of T u 2 → 27T u 2 for the period. (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The displacement of two particles are shown versus dimensionless time t/T u 2 where T u 2 = 2π/ω u 2 is the second NNM period. (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The displacement of the first particle u 1 (t) is shown versus the displacement of the second particle u 2 (t). The maximum amplitudes of each NNM are displayed on the right-top corner of each plot. 80 Elastic unilateral support (kw,cw,d)

Fig. 1 .

 1 Fig. 1. Structural components subjected to unilateral supports. (a) The longitudinal vibration of the bar is constrained by the elastic unilateral obstacle at its right end. The contact element is active when u (L, t) ≥ δ where u (x, t) is the bar axial deformation. (b) The beam-like component is interacting with the elastic unilateral substrate. The contact occurs once w (x, t) ≥ δ for Lx f ≤ x ≤ L where w (x, t) is the beam lateral deflection. (c) Rectangular thin plate subjected to the unilateral Winkler constraint. The contact elements are only active once w (x, y, t) ≥ δ for Lx f ≤ x ≤ L where w (x, t) denotes the transverse displacement of plate middle surface. (d) The top view of the discrete model of the plate.The DOFs of each lattice node is identified as {w i,j (t) , ϕ i,j (t) , θ i,j (t)}. The unilateral support is also modelled as a network of parallel one-way spring-damper elements (K N = k w h x h y , C N = c w h x h y where h x and h y are the length of the unit cells in x and y directions, respectively).

Fig. 2 .

 2 Fig. 2. Sketch of a 3-DOF chain of oscillators in contact with the elastic unilateral constraint. The unilateral obstacle k n is only active when u 3 ≥ δ.

Fig. 3 .

 3 Fig.3. Frequency-stiffness plots for the NNMs of the 3-DOF chain of oscillators with δ = 0 (see Fig.2). The vertical axis shows the dimensionless frequencies ω u i m/k x where the superscript u denotes the unilateral interaction and i = 1, 2, 3. The horizontal axis shows the range of normalised contact stiffness k n /k x . The unstable NNM frequencies are shown with the dashed lines. The orange line also depicts the bilinear approximation for the NNM frequencies obtained from Eq. 38. (a) The first NNM is presented where the backbone branch starts from ω u 1 m/k x = 0.445 on the vertical axis. More details about the first NNM are presented here. The incremental evolution of the NNMs for the representative points (a)-(l) on its backbone is illustrated in Fig.4. New periodic motions are also born at points (d) and (i) where the NNMs on the backbone branch become unstable. The corresponding incremental evolution of these newly born NNMs for representative points (d.1)-(d.11) and (i.1)-(i.11) are depicted in Figs.5 and 6, respectively. (b) The second NNM is presented where its backbone branch starts from ω u 2 m/k x = 1.247 on the vertical axis. (c) The third NNM is presented where its backbone branch starts from ω u 3 m/k x = 1.802 on the vertical axis. The start-points on the vertical axis are the three natural frequencies of the 3-DOF system with k n = 0 which can be obtained from the linear eigenvalue analysis of Eq. 37.

Fig. 4 .

 4 Fig. 4. Trajectories of NNMs for points (a)-(l) on the backbone of the first NNM in Fig 3(a). (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one periodT u 1 = 2π/ω u 1 .The vertical axis shows the displacement components of the 3-DOF chain of oscillators with δ = 0 (see Fig.2). The horizontal axis depicts the dimensionless time t/T u 1 where T u 1 = 2π/ω u 1 is the period. (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The vertical axis shows the displacements of the first and second mass, i.e. u 1 (t) and u 2 (t) respectively. The horizontal axis depicts the displacement of the third mass u 3 (t). The gray coloured plots show the unstable NNMs.

Fig. 5 .

 5 Fig. 5. Trajectories of NNMs for points (d)-(d.11) of the first NNM in Fig 3(a)where the new periodic motion is born at point (d) with a period equal to the double of the fundamental period (T u 1 → 2T u 1 ). (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The vertical axis shows the displacement components of the 3-DOF chain of oscillators with δ = 0 (see Fig.2). The horizontal axis depicts the dimensionless time t/T u

Fig. 6 .

 6 Fig. 6. Trajectories of NNMs for points (i)-(i.11) of the first NNM Fig 3(a) where the new periodic motion is born at point (i) with a period equal to quintuple of the fundamental period (T u 1 → 5T u 1 ). (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The vertical axis shows the displacement components of the 3-DOF chain of oscillators with

Fig. 7 .≈ ω u 3 .

 73 Fig.7. Close-up view of the modes 1-3 (4:1) internal resonance that occurs for the first NNM of the 3-DOF chain of oscillators with δ = 0 (see Fig.2and the corresponding frequency-stiffness plots in Fig.3). The first and third modes start to interact once 4ω u 1 ≈ ω u 3 . The vertical axis shows the first dimensionless frequency ω u 1 m/k x . The horizontal axis shows the range of normalised contact stiffness k n /k x . The unstable NNM frequencies are shown with the dashed line. The orange line also depicts the (ω u 3 /4) m/k x obtained from the shooting method.

Fig. 8 .

 8 Fig. 8. Trajectories of NNMs for points (a)-(l) of Fig 7 where the first and third modes start to interact once 4ω u 1 ≈ ω u 3 and accordingly, the modes 1-3 (4:1) internal resonance occurs. (a) The incremental evolution of the NNMs are presented by the time-history of the free system over one period. The displacement components of the 3-DOF chain of oscillators (see Fig. 2 with δ = 0) are shown with u 1 (t), u 2 (t) and u 3 (t). (b) The incremental evolution of the NNMs are presented in the configuration space over one period. The gray coloured plots show the unstable NNMs.

Fig. 9 .

 9 Fig.9. The cantilever beam is constrained by the unilateral elastic walls on both sides with the gap of δ. The structure has a circular cross-section with the radius of R and it is subjected to the harmonic load of F (t) applied at point x p , i.e. F (x, t) = F (t) δ (xx p ).

Fig. 10 .

 10 Fig. 10. Transient dynamic response of the beam with the circular cross-section constrained by the elastic unilateral walls on both sides (see Fig. 9) with k w L 4 / (EI) = 500, c w = 0, x f /L = 0.4 and δ/R = 0.0125. The structure is subjected to the harmonic excitation of F (t) = F 0 sin (ωt) at x p /L = 0.5 with F 0 = 4000 kN and ωL 2 ρA/ (EI) = 3.5. The vertical axis shows the time-history of the dimensionless displacement w = w/d s in z-direction at x = L/2 (thin solid line for the present discrete model and thin dashed line for the ANSYS results) and x = L (thick solid line for the present discrete model and thick dashed line for the ANSYS results) where d s = F 0 L 3 / (3EI). The horizontal axis depicts the dimensionless time t = tω/ (2π) which is normalised with respect to the period of the applied harmonic load. (a) Present discrete lattice model (DLM) versus the FEM results (ANSYS) for the Euler-Bernoulli theory (ǧ = 0). (b) DLM versus the FEM results for the Timoshenko theory (ǧ = 1).

Fig. 11 . 20 0Fig. 12 .

 112012 Fig. 11. Transient dynamic response of the beam with the circular cross-section constrained by the elastic unilateral walls on both sides (see Fig. 9) with k w L 4 / (EI) = 500, c w = 0, x f /L = 0.4 and three different gap lengths. The structure is subjected to the harmonic excitation of F (t) = F 0 sin (ωt) at x p /L = 0.5 with F 0 = 4000 kN and ωL 2 ρA/ (EI) = 3.5. The results obtained from the Timoshenko beam model (ǧ = 1) are compared with the Euler-Bernoulli theory (ǧ = 0). (a) The vertical axis shows the time-history of the dimensionless displacement w = w/d s in z-direction at x = L/2 and x = L where d s = F 0 L 3 / (3EI). The horizontal axis depicts the dimensionless time t = tω/ (2π) which is normalised with respect to the period of the applied harmonic load. (b) The vertical axis shows the dimensionless lateral displacement throughout the beam length (i.e. the deformed shape of the structure) at two different time instants t = 0.8 and t = 4. The horizontal axis depicts the dimensionless axial coordinate x/L. Note that the gap length for Figs. (a.1) and (b.1) is δ = 0, Figs. (a.2) and (b.2) is δ = 0.0125R, and Figs. (a.3) and (b.3) is δ = 0.025R.

Fig. 15 .

 15 Fig. 15. (a.1)-(i.1) The frequency plots of the granular chain with n = 2, n = 3 and n = 4 particles (see Fig. 14). The dimensionless frequency ωu i = ω u i m/k c is shown as a function of the damage parameter ď = 1k t /k c . The maximum amplitudes of the NNMs at ď = 0, ď = 0.5 and ď = 1 are also presented. (a.2)-(i.2) The time-history of the fundamental modes with k t = 0 (i.e. ď = 1). The displacement of particles are indicated versus dimensionless time t/T u i where T u i = 2π/ω u i .

Fig. 16 .

 16 Fig. 16. The transient response of the tensionless granular chain (k t = 0) due to the harmonic excitation. The vertical axis shows the dimensionless displacement of the particles ūi (t) = u i (t) /u s where u s = 10F 0 /k c . The horizontal axis represents dimensionless time tω/ (2π) where ω is the frequency of the applied force(s). The inset figures show the number and the direction of the applied harmonic loads. (a) The number of particles is n = 2 and the frequency of the applied load is ω = ωu 1 k c /m = 0.739 k c /m. (b) n = 4 and ω = ωu 2 k c /m = 0.829 k c /m. (c) n = 4 and ω = ωu 4 k c /m = 1.281 k c /m. (d) n = 6 and ω = ωu 6 k c /m = 1.323 k c /m.

  -30 and compute the earliest time instant over t r+1t s m that the lattice nodes cross the boundary. This can be achieved by monitoring the diagonal components of H s m+1 versus H s m during the bisection search to capture the earliest instant that they change from 0 to 1 or vice versa (note that the superscript s shows that the contact constraint matrices are associated with the subintervals). The new subinterval obtained from the bisection module is ∆t s

	m
	to perform the numerical integration over t s m , t s m+1 and accept the solution as
	the system response at t = t s m+1 . Subsequently, we use the system output data
	at t s m+1 to calculate H s m+1 and considering H r+1 = H s m+1 , perform the numer-
	ical integration over t s m+1 , t r+1 . If the resulting contact constraint matrix
	supports the validity of the preceding assumption (i.e. H r+1 = H s m+1

m (note that ∆t s m ≤ t r+1t s m ≤ ∆t) and the new local time instant is t s m+1 = t s m + ∆t s m . Now, we assume that H s m+1 = H s

Appendix A Computation of NNMs for non-smooth systems using the shooting method

The shooting technique is commonly used in the literature to numerically obtain the NNMs of a non-linear MDOF setting and present a family of them in the frequency-energy plots. The resulting NNMs form a branch in the 2N + 1dimensional space of initial conditions and oscillation period at different energy levels. Following the principles of this technique (bearing in mind that we are looking for the system NNMs at a certain level of energy), the combination of the Newton-Raphson algorithm and present solution scheme (Section 2.3) can be utilised in an iterative process to find a particular period T ∈ R + and initial conditions z (0) = D 0 Ḋ0 ∈ R 2n that satisfy the periodicity condition of the free response in the phase plane, i.e. z(t, z (0) ) = z(t + T, z (0) ). In order to apply this method to our case, we consider the system energy at a fixed level, as E(z (0) ) = ě where the energy of the conservative system at any time instant is identified as

Setting k = 1, 2, 3, ... as the iteration index number, it is required to examine the validity of periodicity condition at each step of the Newton-Raphson algorithm for the refined period T k and initial conditions z (0)

k . It can be achieved by introducing the shooting function S(T, z (0)

Since the phase of periodic solution is not fixed, we also need to remove the arbitrariness of the initial conditions [START_REF] Peeters | Nonlinear normal modes, part II: Toward a practical computation using numerical continuation techniques[END_REF]. Here, we fix the phase of each nonlinear periodic motion by setting all initial velocities of the discrete points

k < ε where • implies the Euclidean norm and ε is the precision tolerance.

Appendix B The continuation of periodic solution for the non-

smooth system

The pseudo-arclength continuation approach uses a known periodic motion on the branch, i.e. (z (0) j , T j ), to trace the branch and identify the next solution, i.e. (z (0) j+1 , T j+1 ), where subscript j denotes the step number during the path following process. This method includes two steps, i.e. a predictor step and a corrector step. In the predictor step, the tangent of the branch at the current known point j is utilised to generate the initial guess of the shooting algorithm at the next point j + 1 along the tangent as

is the NNM and its period identified using the shooting scheme at point j with the unilateral contact stiffness of k w,j , z (0) 0,j+1

is the initial guess of the shooting algorithm at point j + 1 with the contact stiffness of k w,j+1 , s j ∈ R is the step-size where its sign is selected in order to trace the branch, and

∈ R 2N +1 is the tangent vector to the (2N + 1)-dimensional space of initial conditions and period at point j. The tangent vector is comprised of two components; p z,j ∈ R 2N denoting the tangent to the branch with respect to z 78