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

Jayantheswar Venkatesh, Anders Thorin, Mathias Legrand

► To cite this version:

Jayantheswar Venkatesh, Anders Thorin, Mathias Legrand. Nonlinear modal analysis of a one-dimensional bar undergoing unilateral contact via the time-domain boundary element method. ASME 2017 International Design Engineering Technical Conferences, Aug 2017, Cleveland, United States. 10.1115/DETC2017-68340. hal-01569943

HAL Id: hal-01569943

https://hal.science/hal-01569943

Submitted on 28 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# NONLINEAR MODAL ANALYSIS OF A ONE-DIMENSIONAL BAR UNDERGOING UNILATERAL CONTACT VIA THE TIME-DOMAIN BOUNDARY ELEMENT METHOD

## Jayantheeswar Venkatesh, Anders Thorin, Mathias Legrand

Structural Dynamics and Vibration Laboratory
Department of Mechanical Engineering
McGill University
Montreal, Quebec, H3A 0C3, Canada

### ABSTRACT

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 \pm 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) = \sum_i \phi_i(x)u_i(t)$, where $u_i(t)$ is the $i$-th displacement participation and $\phi_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 frequency-energy 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.

# 1 Problem of interest

A one-dimensional elastic bar of length L, constant cross-sectional area A, Young's Modulus E and mass density ρ is considered. The bar is fixed at x = 0 and subject to unilateral contact conditions at x = L, as shown in Fig. 1. The initial

<figure><img src="image_2.png" /><figcaption>Figure 1: Bar with unilateral contact condition</figcaption></figure>

displacement field of the bar at time t = 0 is $u_0(x)$ and the corresponding initial velocity is $v_0(x)$ where $x \in [0; L]$. The signed distance, or gap function, between the contact node at $x = L$ and the rigid foundation is

$$g(t) = g_0 - u(L, t), \quad \forall t \quad (1)$$

where $g_0$ is the gap at the resting position. When contact occurs ($g = 0$), an elastic wave propagates inside the bar at velocity $c = \sqrt{E/\rho}$. The local equation which dictates the displacement $u(x,t)$ of the one-dimensional bar is

$$\rho A \partial_t^2 u(x,t) - EA \partial_x^2 u(x,t) = 0, \quad \forall x \in ]0; L[, \quad \forall t \ge 0 \quad (2)$$

with the boundary condition

$$u(0, t) = 0, \quad \forall t \ge 0. \quad (3)$$

Contact is described using Signorini's conditions

$$g(t) \ge 0, \quad \partial_x u(L, t) \le 0, \quad g(t) \cdot \partial_x u(L, t) = 0, \quad \forall t \ge 0. \quad (4)$$

These inequalities are responsible for the nonlinear behavior of the dynamics.

The objective is to find the nonlinear modes of the above-described 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 $\forall t \ge 0$ and $\forall x \in [0; L]$, $u(x, t + T) = u(x, t)$.

# 2 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].

## 2.1 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 $\delta$ applied at the source point $\xi$ and time $\tau$, that is the solution of

$$\partial_x^2 u^*(x, t, \xi, \tau) - \frac{1}{c^2} \partial_t^2 u^*(x, t, \xi, \tau) = \delta_{\xi, \tau}. \quad (5)$$

Solving (5) leads to the fundamental solution for this problem [13]

$$u^*(x, t, \xi, \tau) = -\frac{c}{2} H[c(t - \tau) - |x - \xi|] \quad (6)$$

where $H$ is the Heaviside function and $|x - \xi|$ is the distance between the field and source points. Noting that $\xi$ and $x$ are any point in the interval $[0; L]$, the variables $\xi$ and $x$ are interchangeably used when required [9] and the same applies to $t$ and $\tau$ in the time interval. The method of weighted residuals can be applied to Eqn. (2) using $u^*(x, t, \xi, \tau)$ as the weighting function. From here on, $u(x, \tau)$ and $u^*(x, t, \xi, \tau)$ are written as $u$ and $u^*$ respectively, when required by compactness. The weighted residual statement takes the form

$$\int_0^t \int_0^L \partial_x^2 u u^* dx d\tau - \frac{1}{c^2} \int_0^t \int_0^L \partial_t^2 u u^* dx d\tau = 0. \quad (7)$$

Substituting Eqn. (6) in (7), and taking the second weak form, i.e. integrating by parts twice, yields

$$\begin{aligned} 0 = \int_0^t \int_0^L u \partial_x^2 u^* dx d\tau & - \frac{1}{c^2} \int_0^t \int_0^L u \partial_t^2 u^* dx d\tau \\ & + \int_0^t \left( (\partial_x u u^*) \Big|_0^L \right) d\tau - \int_0^t \left( (u \partial_x u^*) \Big|_0^L \right) d\tau \\ & - \frac{1}{c^2} \int_0^L \left( (\partial_t u u^*) \Big|_0^t \right) dx + \frac{1}{c^2} \int_0^L \left( (u \partial_\tau u^*) \Big|_0^t \right) dx. \end{aligned} \quad (8)$$

The fundamental solution features the following properties [9]:

1. Causality: $u^*(x, t, \xi, \tau) = 0$ if $c(t - \tau) < |x - \xi|$;

2. Invariance by time translation: $u^*(x, t, \xi, \tau) = u^*(x, t + t_1, \xi, \tau + t_1)$;

3. Reciprocity: $u^*(x, t, \xi, \tau) = u^*(\xi, -\tau, x, -t)$.

The distributional derivative of the displacement $u^*$ with respect to $x$ and $\tau$

$$\frac{1}{c} \partial_{\tau} u^*(x, t, \xi, \tau) = \partial_x u^*(x, t, \xi, \tau) = \frac{c}{2} \delta_{|x-\xi|-c(t-\tau)} \quad (9)$$

and the causality property of the fundamental solution (8) yield the internal point equation

$$\begin{aligned} u(\xi, t) = \frac{1}{2} & u(L, t - (L - \xi)/c) + \frac{1}{2} u(0, t - \xi/c) \\ & - \int_0^t \partial_x u(L, \tau) u^*(L, t, \xi, \tau) d\tau \\ & - \int_0^t \partial_x u(0, \tau) u^*(0, t, \xi, \tau) d\tau \\ & + \frac{1}{c^2} \int_0^L v_0(x) u^*(x, t, \xi, 0) dx \\ & - \frac{1}{c^2} \int_0^L u_0(x) \partial_\tau u^*(x, t, \xi, 0) dx \end{aligned} \quad (10)$$

where $u(x, \cdot)$ is extended for negative $t$ by $u(x, t) = 0$. The last integral in Eqn. (10) should be understood in the following
distributional way using Eqn. (9):

$$ \int_{0}^{L} u_{0}(x) \partial_{\tau} u^{*}(x, t, \xi, 0) dx = $$

$$ \begin{cases} u_0(\xi - ct) & \text{if } \xi - ct \ge 0 \text{ and } \xi + ct > L \\ u_0(\xi + ct) & \text{if } \xi - ct < 0 \text{ and } \xi + ct \le L \\ u_0(\xi - ct) + u_0(\xi + ct) & \text{if } \xi - ct \ge 0 \text{ and } \xi + ct \le L \\ 0 & \text{otherwise.} \end{cases} $$

The integrals over $[0; 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 $[0; t]$ are dealt with by considering a time discretization scheme with $n$ time steps of length $\Delta\tau$; between two successive time steps, $\partial_x u(0, \tau)$ and $\partial_x u(L, \tau)$ 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)$, $\partial_x u(0, t)$, $\partial_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 $\xi = 0$ and $\xi = L$, leading to two linear equations at every instant $t_i = i \Delta t$ which can be gathered in the matrix form

$$ \mathbf{H}\mathbf{u} = \mathbf{G}\partial_x\mathbf{u} - \mathbf{b} \qquad (11) $$

where $\mathbf{H}$ is a square matrix of dimension $2 \times 2$, $\mathbf{G}$ is a rectangular $2 \times 2n$ matrix and $\mathbf{b}$ is the vector computed from the two last terms of Eqn. (10). Quantity $\mathbf{u}$ is the vector $(2 \times 1)$ of boundary displacements at instant $t_i$ and $\partial_x \mathbf{u}$ is the vector $(2n \times 1)$ of boundary tractions computed from the time integration over $[0; t_i]$. Equation (11) can be solved for the two unknown boundary conditions at $t_i$ stacked in either $\mathbf{u}$ or $\partial_x \mathbf{u}$ or both, and then inserted back into Eqn. (10) to recover the solution.

## 2.2 Switching boundary conditions

The complementarity conditions (4) at $x = L$ are accounted for by switching the boundary conditions at the time instants when the gap opens or closes. The time discretisation may lead to undesirable residual penetration of the contacting end of the bar into the rigid foundation. Such penetrations are handled by projecting the contacting end of the bar on the rigid foundation. Doing so creates an extra shock wave, reducing the accuracy of the solution; choosing an appropriate time step leading to small penetration is therefore essential. Overall, at any given instant $t_i$, one of the following conditions applies:

1. contact is not activated: $g(t_i) > 0$ and $\partial_x u(L, t_i) = 0$,

2. contact or penetration occurs: $g(t_i) \le 0$ and

(a) contact remains closed: $\partial_x u(L, t_i) \le 0$,

(b) contact is released: $\partial_x u(L, t_i) > 0$.

The darker grey portion in Fig. 2 shows the bar with zero displacement and the lighter grey portion shows the stretched bar at a given instant $t_i$.

**Open Gap** The sign of the gap defined in Eqn. (1) is monitored at every time step $g(t_i)$. If it is positive, the contacting node at $x = L$ is free: $\partial_x u(L, t_i) = 0$.

**Penetration or Contact** If $g(t_i) = 0$, no adjustment is required. When $g(t_i) < 0$ (penetration occurs), the displacement of the contacting node is adjusted to satisfy the Signorini conditions:

$u(L, t_i) = 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) = 0$ and the contacting force is positive, i.e. reaction force $\partial_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. $\partial_x u(L, t_i) > 0$, the contacting node is released and the gap will be open at the next time step $t_{i+1}$. This is modelled by a free boundary condition.

<figure><img src="image_3.png" /><figcaption>Figure 2: Adjustment in case of penetration</figcaption></figure>

**Input:** total number of time steps $n$, initial gap $g_0$, boundary conditions, initial conditions $u_0$ and $v_0$, wave velocity $c = \sqrt{E/\rho}$

1. **for** $i = 1:n$ **do**

2. Compute $g(t_i) = g_0 - u(L, t_i)$ using Eqn. (11)

3. -- Switching Boundary Conditions --

4. **if** $g(t_i) \ge 0$ **then**

5. Set free boundary condition at $x = L$

6. **else**

7. $u(L, t_i) \leftarrow g_0$

8. Set fixed boundary condition at $x = L$

9. Compute $\partial_x u(L, t_i)$ using Eqn. (11)

10. **end**

11. -- End Switching Boundary Conditions --

12. **end**

**Output:** $u(L, t_i)$, $\partial_x u(L, t_i)$, $\partial_x u(0, t_i)$

## Algorithm 1: Unilateral contact in TD-BEM

The proposed algorithm 1 summarizes the approach. It is essentially a time-marching procedure where the gap and the contact force are computed. Once all the boundary values are found, the internal displacements can be computed using Eqn. (10).

## 2.3 Validation of the proposed algorithm

The numerical properties of TD-BEM are illustrated on a one-dimensional 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.
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Bouncing bar solution: analytical solution [—], TD-BEM [---] and FEM [—]</figcaption></figure>

TD-BEM accurately captures the traveling waves propagating at a finite speed and generated by unilateral contact, with no spurious oscillations. Energy is conserved over time, as opposed to the chosen FEM simulation. TD-BEM does not necessitate an impact law to retrieve the exact solution. This allows for both non-impulsive lasting contact and energy preserving solutions.

### 3 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) = u(x, T) \quad \text{and} \quad v_0(x) = v(x, T). \qquad (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 = 0$, Eqn. (12) reduces to just solving $u_0(x) = 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) \approx u_{0i}$, $i = 1, \dots, N$, denoted $\mathbf{u}_0$. Similarly, the displacement at $T$ is approximated by its values $u_i$, $i = 1, \dots, N$, denoted $\mathbf{u}$. This last quantity is computed from the unknowns $\mathbf{u}_0$ and $T$ using the above-described TD-BEM. Periodicity with zero initial velocity is enforced by solving

$$\mathbf{f}(\mathbf{u}_0, T) = \mathbf{u}_0 - \mathbf{u}(\mathbf{u}_0, T) = \mathbf{0} \qquad (13)$$

where $\mathbf{f}: \mathbb{R}^{N+1} \to \mathbb{R}^N$ since $N$ independent equations are generally provided for $N+1$ 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 = 1, \dots, 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® 2015.

### 4 Results

Young's Modulus $E$, mass density $\rho$ and length of the bar $L$ are arbitrarily chosen equal to one and the initial gap is chosen as $g_0 = 0.001$ so that $g_0 \ll L$. The resonant frequencies of undamped linear systems are independent of the vibratory energy: this corresponds to vertical backbone curves in the frequency-energy diagram. This no longer holds for nonlinear systems, as illustrated by the first two nonlinear modes of vibration of the system of interest, see Fig. 4. The backbone curves show the main vibratory resonances in the vicinity of $\omega_1(E)$ and $\omega_2(E)$, subharmonic resonances near $\omega_2(E)/2$ and internal resonances observed along the first nonlinear backbone curve.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Backbone curves via TD-BEM & shooting. linear mode [—], nonlinear main vibratory response [—], subharmonic of the second mode ($\omega_2/2$) [—], internal resonance [---]. labels are referred to in the following figures</figcaption></figure>

Subharmonic resonances and internal resonances are typical of nonlinear dynamics and cannot be observed in linear systems. 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.

#### 4.1 Main vibratory responses

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
<figure><img src="image_6.png" /><figcaption>**Figure 5:** First nonlinear mode trajectory: displacement of contacting node [—] and corresponding unilateral contact force [—]</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6:** First nonlinear mode space-time trajectory: a in Fig. 4</figcaption></figure>

time functions with a clear indication of the characteristic lines. 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.

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Second nonlinear mode space-time trajectory with one node in space: b in Fig. 4</figcaption></figure>

## 4.2 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 $\omega_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.

<figure><img src="image_9.png" /><figcaption>**Figure 8:** Constant energy continuum at every point of the subharmonic backbone curve in the vicinity of $\omega_2/2$: limits of the continuum (grazing) [—] & second mode over two periods [—]</figcaption></figure>

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\pi/\omega_2 = 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.

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Grazing subharmonic nonlinear mode space-time trajectory: c in Fig. 4</figcaption></figure>

## 4.3 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
<figure><img src="image_11.png" /><figcaption>**Figure 10:** Internally resonant space-time trajectory emanating from the first main backbone curve: *d* in Fig. 4</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11:** Displacement of the contacting end of the bar: first nonlinear mode [—] and internal resonant mode [—]. A low-magnitude high-frequency solution combines with the first nonlinear mode of vibration</figcaption></figure>

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.

## Nomenclature

$L, \rho, A$ length, mass density, cross-sectional area of the bar

$c$ wave velocity

$E$ Young's modulus of the bar

$x$ field point in space

$\xi$ source point in space

$t$ time

$\tau$ source point in time

$T$ time period of oscillation

$\Delta t$ time-step

$n$ total number of time-steps

$g$ time dependent gap function

$g_0$ initial gap

$u(x, t)$ space-time displacement field

$\partial_t, \partial_x$ first derivative with respect to $t, x$

$\partial^2_t, \partial^2_x$ second derivative with respect to $t, x$

$\partial_{\tau}$ first derivative with respect to $\tau$

$EA\partial_x u$ internal traction within the bar

$u_0(x)$ initial displacement

$v_0(x)$ initial velocity

$u^*(x, t, \xi, \tau)$ fundamental solution

$\delta$ Dirac distribution

$H(x)$ Heaviside function

## Acknowledgment

This research was sponsored by the Natural Sciences and Engineering Research Council of Canada (NSERC) and Fonds de Recherche du Québec Nature et Technologie (FQRNT) through the “Discovery” and “Établissement de nouveaux chercheurs universitaires” programs.

## References

[1] ABREU, Ana, CARRER, Jose Antonio Marques, and MANSUR, Webe João. “Scalar wave propagation in 2D: a BEM formulation based on the operational quadrature method”. *Engineering Analysis with Boundary Elements*, 2003, 27(2):101–105. [hal-01580922].

[2] AGNANTIARIS, John, POLYZOS, Demosthenes, and BESKOS, Dimitri. “Three-dimensional structural vibration analysis by the dual reciprocity BEM”. *Computational Mechanics*, 1998, 21(4):372–381. [hal-01581502].

[3] ALLGOWER, Eugene and GEORG, Kurt. *Numerical Continuation Methods: an Introduction*. Vol. 13. Springer Science & Business Media, 2012. ISBN: 9783642612572.

[4] ARQUIER, Rémi, BELLIZZI, Sergio, BOUC, Robert, and COCHE-LIN, Bruno. “Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes”. *Computers & Structures*, 2006, 84(24):1565–1576. [hal-00088183].

[5] BALLARD, Patrick and BASSEVILLE, Stéphanie. “Existence and uniqueness for dynamical unilateral contact with Coulomb friction: a model problem”. *ESAIM: Mathematical Modelling and Numerical Analysis*, 2005, 39(1):59–77. [hal-00461828].

[6] CARPENTER, Nicholas, TAYLOR, Robert, and KATONA, Michael. “Lagrange constraints for transient finite element surface contact”. *International Journal for Numerical Methods in Engineering*, 1991, 32(1):103–128. [hal-01389918].
[7] CARRER, Jose Antonio Marques and COSTA, Vivianne Luiza. "Boundary Element Method formulations for the solution of the scalar wave equation in one-dimensional problems". *Journal of the Brazilian Society of Mechanical Sciences and Engineering*, 2015, 37(3):959-971. [hal-01581476].

[8] CARRER, Jose Antonio Marques and MANSUR, Webe João. "Alternative time-marching schemes for elastodynamic analysis with the domain boundary element method formulation". *Computational Mechanics*, 2004, 34(5):387-399. [hal-01581478].

[9] CARTWRIGHT, David. *Underlying Principles of the Boundary Element Method*. WIT, 2001. ISBN: 9781853128394.

[10] DOMINGUEZ, Jose. *Boundary Elements in Dynamics*. WIT Press, 1993. ISBN: 9781853122583.

[11] DOYEN, David, ERN, Alexandre, and PIPERNO, Serge. "Time-integration schemes for the finite element dynamic Signorini problem". *SIAM Journal on Scientific Computing*, 2011, 33(1):223-249. [hal-00440128].

[12] DUFFY, Dean. *Green's Functions with Applications*. CRC Press, 2015. ISBN: 9781138894464.

[13] GRAFF, Karl. *Wave Motion in Elastic Solids*. Dover Publications, 1991. ISBN: 9780486667454.

[14] KERSCHEN, Gaëtan, PEETERS, Maxime, GOLINVAL, Jean Claude, and VAKAKIS, Alexander. "Nonlinear normal modes, Part I: A useful framework for the structural dynamicist". *Mechanical Systems and Signal Processing*, 2009, 23(1):170-194. [hal-01357931].

[15] LAXALDE, Denis and LEGRAND, Mathias. "Nonlinear modal analysis of mechanical systems with frictionless contact interfaces". *Computational Mechanics*, 2011, 47(4):469-478. [hal-00492775].

[16] LAXALDE, Denis and THOUVEREZ, Fabrice. "Complex nonlinear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces". *Journal of Sound and Vibration*, 2009, 322(4):1009-1025. [hal-00343494].

[17] LOEFFLER, Carlos, NEVES, Felipe Patricio, and OLIVEIRA, Patrícia. "Dual reciprocity boundary element formulation applied to non linear darcian diffusive advective problems". *Boundary Element Techniques XI. Proceedings of the 11th International Conference*. EC Ltd, 2010. [hal-01580915].

[18] LUBICH, Christian. "Convolution quadrature and discretized operational calculus". *Numerische Mathematik*, 1988, 52(2):129-145. [eudml-133229].

[19] MANSUR, Webe João. "A time-stepping technique to solve wave propagation problems using the boundary element method". PhD thesis. University of Southampton, 1983.

[20] OUEINI, Shafic, CHIN, Char-Ming, and NAYFEH, Ali. "Dynamics of a cubic nonlinear vibration absorber". *Nonlinear Dynamics*, 1999, 20(3):283-295. [hal-01580924].

[21] PEETERS, Maxime, VIGUIÉ, Régis, SÉRANDOUR, Guillaume, KERSCHEN, Gaëtan, and GOLINVAL, Jean-Claude. "Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques". *Mechanical Systems and Signal Processing*, 2009, 23(1):195-216. [hal-01581480].

[22] PESHECK, Eric, PIERRE, Christophe, and SHAW, Steven. "A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds". *Journal of Sound and Vibration*, 2002, 249(5):971-993. [hal-01580920].

[23] SAITOH, Takahiro, HIROSE, Sohichi, and FUKUI, Takuo. "Convolution quadrature time-domain boundary element method and acceleration by the fast multipole method in 2D viscoelastic wave propagation". *Theoretical and Applied Mechanics Japan*, 2009, 57:385-393. [hal-01581489].

[24] SCHANZ, Martin. *Wave Propagation in Viscoelastic and Poroelastic Continua: A Boundary Element Approach*. Vol. 2. Springer Science & Business Media, 2012.

[25] THORIN, Anders, LEGRAND, Mathias, and JUNCA, Stéphane. "Nonsmooth modal analysis: Investigation of a 2-dof spring-mass system subject to an elastic impact law". *Proceedings of the ASME IDETC/CIE. 11th International Conference on Multibody Systems, Nonlinear Dynamics, and Control*. Vol. 6. 2015. [hal-01185973].

[26] VAKAKIS, Alexander, MANEVITCH, Leonid, MIKHLIN, Yuri, PILIPCHUK, Valery, and ZEVIN, Alexandr. *Normal Modes and Localization in Nonlinear Systems*. Springer, 2001. ISBN: 9789401724524.

[27] VLADIMIROV, Vasiliï Sergeevich. *Generalized Functions in Mathematical Physics*. MIR, 1979. ISBN: 9780714715452.

[28] YOONG, Carlos, THORIN, Anders, and LEGRAND, Mathias. "The Wave Finite Element Method applied to a one-dimensional linear elastodynamic problem with unilateral constraints". *Proceedings of the ASME IDETC/CIE. 11th International Conference on Multibody Systems, Nonlinear Dynamics, and Control*. Vol. 6. 2015. [hal-01194922].