Nonlinear modal analysis of a one-dimensional bar undergoing unilateral contact via the time-domain boundary element method

Finite elements in space with time-stepping numerical schemes, even though versatile, face theoretical and numerical difficulties when dealing with unilateral contact conditions. In most cases, an impact law has to be introduced to ensure the uniqueness of the solution: total energy is either not preserved or spurious high-frequency oscillations arise. In this work, the Time Domain Boundary Element Method (TD-BEM) is shown to overcome these issues on a one-dimensional system undergoing a unilateral Signorini contact condition. Unilateral contact is implemented by switching between free boundary conditions (open gap) and fixed boundary conditions (closed gap). The solution method does not numerically dissipate energy unlike the Finite Element Method and properly captures wave fronts, allowing for the search of periodic solutions. Indeed, TD-BEM relies on fundamental solutions which are travelling Heaviside functions in the considered one-dimensional setting. The proposed formulation is capable of capturing main, subharmonic as well as internal resonance backbone curves useful to the vibration analyst. For the system of interest, the nonlinear modeshapes are piecewise-linear unseparated functions of space and time, as opposed to the linear modeshapes that are separated half sine waves in space and full sine waves in time. Introduction In structural dynamics, autonomous conservative systems most commonly exhibit continuous families of periodic orbits in the phase space, usually referred to as modes of vibration. A major task of modal analysis is to accurately compute natural frequencies and corresponding mode shapes as they are known to properly predict the frequencies under which the associated periodically forced systems will resonate, at least in linear and smooth nonlinear frameworks. Characterizing the modes of vibration of smooth nonlinear mechanical systems (systems governed by PDEs that are smooth with respect to the unknown displacement and velocity) is a current topic of interest in the academic and industrial spheres. It relies on computing periodic solutions which are sensitive to numerical accuracy. For example, a dissipative time-marching numerical scheme cannot describe autonomous periodic solutions. The dynamics of two impacting bodies is characterized by the travelling waves emanating from the contact interface. In the one-dimensional setting chosen in this work, these waves couple time and space, in the sense that they are functions of the form f .x  ̇ ct/ where c is the wave velocity. Uncoupling time t and space x leads to numerical and theoretical issues. In the Finite Element Method (FEM), the displacement commonly takes the form u.x; t/ D P i i .x/ui .t/, where ui .t/ is the i-th displacement participation and i .x/, the corresponding shape function. This leads to spurious oscillations, dispersion, and energy dissipation, for most numerical schemes dealing with unilateral contact conditions [11]. Additionally, an impact law is required to uniquely describe the time-evolution of a space semi-discretized formulation [5]. The impact law should be purely elastic to preserve energy, making it difficult to describe lasting contact phases which are expected to exist in the continuous framework. The Wave Finite Element Method (WFEM), which appropriately combines space and time, has shown promising results for one-dimensional systems undergoing contact conditions [28]. In this work, a variant of the Boundary Element Method (BEM), called the Time Domain Boundary Element Method (TDBEM) [19], is used to solve for the nonlinear modes [14] of the one-dimensional bar with fixed boundary at one end and a unilateral contact condition on the other. Both BEM and FEM are weighted residual methods but the weighting function in BEM is defined as the fundamental solution [27] of the considered PDE. The fundamental solution is defined as the response of a body subjected to a Dirac delta input, irrespective of the boundary conditions. When boundary conditions are included, the fundamental solution transforms to the classical Green’s function [12]. It is then straightforward to compute the response of a linear system to any distributed body forces and any boundary conditions through the principle of superposition reflected by a convolution operation [9]. However, a major limitation of BEM is that fundamental solutions are known exactly only for simple PDEs. In general, they can only be approximated thus reducing the accuracy of BEM. Various types of BEM are available in literature such as Domain-BEM (D-BEM) [8], Time Domain-BEM (TD-BEM), Dual Reciprocity-BEM (DR-BEM) [2, 17], Frequency DomainBEM (FD-BEM) [10], Convolution Quadrature-BEM (CQBEM) [1, 24]. D-BEM and DR-BEM both discretize time and space separately, hence are not of interest here. In contrast, CQBEM and TD-BEM provide a formulation in space and time allowing to precisely capture traveling waves [19], at least for one-dimensional problems in space. CQ-BEM, first introduced by Lubich [18] and later used for transient analysis [1, 24] differs from TD-BEM in the way the integrals are computed, i.e. using the Convolution Quadrature Method (CQM). Application of CQM to the TD-BEM improves the numerical stability of the


Introduction
In structural dynamics, autonomous conservative systems most commonly exhibit continuous families of periodic orbits in the phase space, usually referred to as modes of vibration. A major task of modal analysis is to accurately compute natural frequencies and corresponding mode shapes as they are known to properly predict the frequencies under which the associated periodically forced systems will resonate, at least in linear and smooth nonlinear frameworks.
Characterizing the modes of vibration of smooth nonlinear mechanical systems (systems governed by PDEs that are smooth with respect to the unknown displacement and velocity) is a current topic of interest in the academic and industrial spheres. It relies on computing periodic solutions which are sensitive to numerical accuracy. For example, a dissipative time-marching numerical scheme cannot describe autonomous periodic solutions.
The dynamics of two impacting bodies is characterized by the travelling waves emanating from the contact interface. In the one-dimensional setting chosen in this work, these waves couple time and space, in the sense that they are functions of the form f .x˙ct / where c is the wave velocity. Uncoupling time t and space x leads to numerical and theoretical issues. In the Finite Element Method (FEM), the displacement commonly takes the form u.x; t/ D P i i .x/u i .t /, where u i .t / is the i -th displacement participation and i .x/, the corresponding shape function. This leads to spurious oscillations, dispersion, and energy dissipation, for most numerical schemes dealing with unilateral contact conditions [11]. Additionally, an impact law is required to uniquely describe the time-evolution of a space semi-discretized formulation [5]. The impact law should be purely elastic to preserve energy, making it difficult to describe lasting contact phases which are expected to exist in the continuous framework. The Wave Finite Element Method (WFEM), which appropriately combines space and time, has shown promising results for one-dimensional systems undergoing contact conditions [28].
In this work, a variant of the Boundary Element Method (BEM), called the Time Domain Boundary Element Method (TD-BEM) [19], is used to solve for the nonlinear modes [14] of the one-dimensional bar with fixed boundary at one end and a unilateral contact condition on the other.
Both BEM and FEM are weighted residual methods but the weighting function in BEM is defined as the fundamental solution [27] of the considered PDE. The fundamental solution is defined as the response of a body subjected to a Dirac delta input, irrespective of the boundary conditions. When boundary conditions are included, the fundamental solution transforms to the classical Green's function [12]. It is then straightforward to compute the response of a linear system to any distributed body forces and any boundary conditions through the principle of superposition reflected by a convolution operation [9]. However, a major limitation of BEM is that fundamental solutions are known exactly only for simple PDEs. In general, they can only be approximated thus reducing the accuracy of BEM.
Various types of BEM are available in literature such as Domain-BEM (D-BEM) [8], Time Domain-BEM (TD-BEM), Dual Reciprocity-BEM (DR-BEM) [2,17], Frequency Domain-BEM (FD-BEM) [10], Convolution Quadrature-BEM (CQ-BEM) [1,24]. D-BEM and DR-BEM both discretize time and space separately, hence are not of interest here. In contrast, CQ-BEM and TD-BEM provide a formulation in space and time allowing to precisely capture traveling waves [19], at least for one-dimensional problems in space. CQ-BEM, first introduced by Lubich [18] and later used for transient analysis [1,24] differs from TD-BEM in the way the integrals are computed, i.e. using the Convolution Quadrature Method (CQM). Application of CQM to the TD-BEM improves the numerical stability of the solution [23]. However, it has a low computational efficiency for large scale problems. For these reasons, TD-BEM is chosen in this work.
First, the investigated mechanical system is described. Then, the simulation methods used to compute the time evolution of the system are detailed, and a comparison with a benchmark problem is provided to validate the methodology and illustrate its accuracy. This is followed by a brief explanation of the shooting technique used to find the periodic solutions. The frequencyenergy plot and the mode shapes of main vibratory response, subharmonic response as well as the internal resonances of the system of interest are presented and discussed.

Problem of interest
where g 0 is the gap at the resting position. When contact occurs (g D 0), an elastic wave propagates inside the bar at velocity c D p E= . The local equation which dictates the displacement u.x; t/ of the one-dimensional bar is A@ 2 t u.x; t / EA@ 2 x u.x; t/ D 0; 8x 2 0 I LOE; 8t 0 (2) with the boundary condition u.0; t/ D 0; 8t 0: (3) Contact is described using Signorini's conditions These inequalities are responsible for the nonlinear behavior of the dynamics. The objective is to find the nonlinear modes of the abovedescribed system, defined as continuous families of periodic orbits. Formally, the goal is to find functions u satisfying (2), (3) and (4) together with real numbers T > 0, such that 8t 0 and 8x 2 OE0 I L, u.x; t C T / D u.x; t/.

Simulation methods
This section introduces the background of the one-dimensional TD-BEM, including the algorithm used to implement unilateral contact conditions. The methodology is then validated using a benchmark problem [11].

Formulation of TD-BEM
In TD-BEM, a time-dependent fundamental solution of the PDE (2) is used. The fundamental solution u captures, at the field point x and time t , the effect of a unit impulse ı applied at the source point and time , that is the solution of Solving (5) leads to the fundamental solution for this problem [13] u .
where H is the Heaviside function and jx j is the distance between the field and source points. Noting that and x are any point in the interval OE0 I L, the variables and x are interchangeably used when required [9] and the same applies to t and in the time interval. The method of weighted residuals can be applied to Eqn. (2) using u .x; t; ; / as the weighting function. From here on, u.x; / and u .x; t; ; / are written as u and u respectively, when required by compactness. The weighted residual statement takes the form Substituting Eqn. (6) in (7), and taking the second weak form, i.e. integrating by parts twice, yields The fundamental solution features the following properties [9]: 1. Causality: u .x; t; ; / D 0 if c.t / < jx j; 2. Invariance by time translation: u .x; t; ; / D u .x; t C t 1 ; ; C t 1 /; 3. Reciprocity: u .x; t; ; / D u . ; ; x; t/. The distributional derivative of the displacement u with respect to x and and the causality property of the fundamental solution (8) yield the internal point equation where u.x; / is extended for negative t by u.x; t/ D 0. The last integral in Eqn. (10) should be understood in the following distributional way using Eqn. (9): The integrals over OE0 I L are dealt with by discretizing this interval into sub-intervals onto which are defined piecewise-linear polynomials [7] to approximate u 0 .x/ and v 0 .x/. This discretization is used to find the unknown functions u 0 and v 0 . In a similar fashion, time integrals over OE0 I t are dealt with by considering a time discretization scheme with n time steps of length ; between two successive time steps, @ x u.0; / and @ x u.L; / are approximated by constant interpolation functions. Equation (10) shows that u.x; t/ is formulated as a linear combination of the boundary conditions u.0; t/, u.L; t /, @ x u.0; t /, @ x u.L; t /, u 0 .x/ and v 0 .x/. Exactly half of these boundary conditions are unknown and need to be calculated. This is done by evaluating Eqn. (10) at D 0 and D L, leading to two linear equations at every instant t i D it which can be gathered in the matrix form where H is a square matrix of dimension 2 2, G is a rectangular 2 2n matrix and b is the vector computed from the two last terms of Eqn. (10). Quantity u is the vector (2 1) of boundary displacements at instant t i and @ x u is the vector (2n 1) of boundary tractions computed from the time integration over OE0I t i . Equation (11) can be solved for the two unknown boundary conditions at t i stacked in either u or @ x u or both, and then inserted back into Eqn. (10) to recover the solution.

Switching boundary conditions
The complementarity conditions (4)  When g.t i / < 0 (penetration occurs), the displacement of the contacting node is adjusted to satisfy the Signorini conditions: u.L; t i / D g 0 , as shown in Fig. 2. In both cases, the contact force can be computed from the reaction force exerted by the wall. Two cases are to be considered depending on the sign of this contact force.
Lasting Contact If g.t i / D 0 and the contacting force is positive, i.e. reaction force @ x u.L; t i / is negative, the contact will remain at the next time step. This is modelled by a fixed boundary condition. Release When the contacting force becomes negative, i.e. @ x u.L; t i / > 0, the contacting node is released and the gap will be open at the next time step t i C1 . This is modelled by a free boundary condition.

Validation of the proposed algorithm
The numerical properties of TD-BEM are illustrated on a onedimensional bar, bouncing on a rigid foundation and subjected to constant external body force. One end of the bar is free and the other undergoes unilateral contact conditions. For some specific parameters, the bar bounces periodically against the rigid foundation [11]. The displacement of the contacting node of the bar is compared with the analytical solution [11] and FEM with forward Lagrange multipliers with an explicit time-marching technique [6] in Fig. 3. Time steps are chosen such that both computation time are comparable.

Autonomous periodic solutions
Nonlinear modal analysis helps understand the vibratory signature of nonlinear dynamical systems [26]. Various techniques and tools exist in the literature to compute the nonlinear modes, such as asymptotic-numerical methods [4], invariant manifold techniques [22], Fourier methods [16] and shooting [15,21]. To characterize nonlinear modes, we compute periodic solutions. Because the system is deterministic, it suffices to verify that the initial displacement u 0 and the initial velocity v 0 repeat themselves after a period T to be found, that is u 0 .x/ D u.x; T / and v 0 .x/ D v.x; T /: (12) In this work, initial velocity is assumed to be zero. This implies the existence of an axis of symmetry in the solution explaining the mode shapes observed in the next section. Shooting and TD-BEM are employed to find the sought families of periodic solutions. Since v 0 D 0, Eqn. (12) reduces to just solving u 0 .x/ D u.x; T /. The space domain is discretized into N 1 cells with N nodes. The initial displacement is then approximated as u 0 .x i / u 0i , i D 1; : : : ; N , denoted u 0 . Similarly, the displacement at T is approximated by its values u i , i D 1; : : : ; N , denoted u. This last quantity is computed from the unknowns u 0 and T using the above-described TD-BEM. Periodicity with zero initial velocity is enforced by solving where f W R N C1 ! R N since N independent equations are generally provided for N C1 unknowns. Accordingly, the solution space is expected to be a one-dimensional manifold [3]. However, it was observed that in the subharmonic case, the N equations provide N 1 independent equations, yielding a two-dimensional manifold.
Eqn. (13) is solved using a Newton's solver that shoots for values of initial displacement u 0i , i D 1; : : : ; N . Since the system is underdetermined, outputs of the Newton's solver are elements of a continuum of solutions. Parametric continuation is employed to recover this continuum of solutions, starting from a known solution, the limit case linear grazing mode. When parametric continuation misses the solution as frequency increases, a more sophisticated arc-length continuation is used instead. The TD-BEM solver as well are arc-length continuation are implemented using MATLAB R 2015.

Results
Young's Modulus E, mass density and length of the bar L are arbitrarily chosen equal to one and the initial gap is chosen as g 0 D 0:001 so that g 0 L. The resonant frequencies of undamped linear systems are independent of the vibratory energy: this corresponds to vertical backbone curves in the frequencyenergy diagram. This no longer holds for nonlinear systems, as illustrated by the first two nonlinear modes of vibration of the system of interest, see  They are briefly discussed in the sequel. The backbone curves in Fig. 4 have a vertical part and a curved part. The vertical part corresponds to the linear mode and denotes that contact is not activated. The energy is frequency-independent until contact is initiated which gives rise to a non-straight backbone curve. Linear mode shapes of a fixed-free bar are standing sine waves, but this no longer holds when a contact nonlinearity is introduced: instead, travelling waves are observed because contact induces shock waves. Figure 5 shows the displacement of the contacting node and contact force over one period of the first nonlinear mode. The contacting node first travels freely (fixed-free bar), then hits the rigid foundation and remains on it (fixed-fixed bar); eventually, the contact force vanishes and the contacting end is released. Figures 6 and 7 show the motion corresponding to first and second modes over one period. Mode shapes are no longer separated half sine waves in space and full sine waves in time as in the linear case. Instead, they are unseparated piecewise-linear space Also, the surface plots show that the two modes are different, even though the displacement of the contacting node has a similar pattern. In particular, one point of the bar is a node for the second mode: it has a constant zero displacement over time, see the blue line in Fig. 7. In this respect, it is similar to the second mode of the linear system.  Fig. 4

Subharmonic response
Subharmonic resonances, defined as the resonances at special frequencies equal to an integer sub-multiple of the natural frequencies, exist only in nonlinear systems [14]. The second mode, considered over two periods, defines a new periodic trajectory of frequency ! 2 =2. In the energy-frequency graph, this corresponds to a new backbone curve similar to the second mode, but of half frequency (and same energy), called subharmonic backbone curve. A phenomenon, new in the continuous framework, is observed along this subharmonic backbone curve. As illustrated in Fig. 8, a continuum of periodic orbits is observed at every single point in the subharmonic curve with a minimal period 2T 2 . For a given frequency, co-existing solutions are found with identical energy but with distinct shapes. A similar property named bridge is reported for a two-dof vibro-impact spring-mass system [25].
A trajectory is said to graze when the contacting end reaches the rigid foundation with zero velocity and recedes away without lasting contact: it is the limit case between no contact and contact. The shaded portion in Fig. 8 shows the previously mentioned continuum delimited by two solutions with one contact phase and one grazing instant per period, which are actually the same but shifted by a duration of 2 =! 2 D T 2 . Every other solution in the continuum has two contact phases per period. Fig. 9 shows the motion corresponding to grazing solution over one period.  Fig. 4

Internal resonance
Another phenomenon existing only in nonlinear systems is the internal resonance. In some experiments with nonlinear systems, the excitation of a mode at a frequency actuates the response of a distinct higher frequency mode. This interaction property has been used to design vibration absorbers [20] for instance. Figure 10 shows an internal resonance of a high-frequency mode interacting with a lower frequency mode: the mode shape exhibits a large similarity with the first mode shape but also includes high-frequency content. Figure 11 displays the displacement of the contacting end for the first mode and the internal resonant mode with the same frequency of vibration. In the internal resonant case, a mode with A lowmagnitude high-frequency solution combines with the first nonlinear mode of vibration frequency about 12 times (see Fig. 11, the wave has 12 nodes in one period) the frequency of the first non-smooth mode seems to interact with the first nonlinear mode.

Conclusions
The periodic autonomous dynamics of a one dimensional bar fixed on one end and subject to unilateral contact conditions on the other was investigated. Periodic solutions were targeted in order to build the nonlinear modes of vibration. Unilateral contact conditions give rise to travelling waves which cannot be accurately captured using FEM. In contrast, TD-BEM formulated in space-time domain showed promising numerical characteristics in capturing travelling wave phenomenon.
First, TD-BEM with boundary conditions depending on the contact state was shown to simulate the time-evolution of a bouncing bar with high accuracy, opening doors to the search of periodic solutions of unilateral contact problems. Such periodic solutions were computed via an implementation of TD-BEM within a shooting method, and continuation techniques were used to recover the whole modes.
Backbone curves in the energy-frequency diagram were presented for the first two modes. One of such curves is a subharmonic curve of the second mode. The backbone curve of the subharmonic mode was shown to correspond to a two-dimensional continuum of periodic solutions, delimited by two grazing solutions and centered around the second mode. Vertical branches emanating from the first mode backbone curve were found to correspond to internal resonances. The periodic motions associated to such internal resonances were computed and these are helpful in predicting the possibility of sudden resonances in real life applications, when vibrating in the vicinity of these frequencies.
The next step will consist in extending the presented methodology to higher dimensions in space [10]. Future works also include stability analysis of the computed modes.