# Nonsmooth modal analysis of a non-internally resonant finite bar subject to a unilateral contact constraint

Carlos Yoong, Mathias Legrand

► To cite this version:

Carlos Yoong, Mathias Legrand. Nonsmooth modal analysis of a non-internally resonant finite bar subject to a unilateral contact constraint. 37th IMAC: A Conference and Exposition on Structural Dynamics, Society for Experimental Mechanics, Jan 2019, Orlando, United States. pp.1-10, 10.1007/978-3-030-12391-8_1. hal-01918511v2

HAL Id: hal-01918511

https://hal.science/hal-01918511v2

Submitted on 29 Jul 2019

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
# Nonsmooth modal analysis of a non-internally resonant finite
bar subject to a unilateral contact constraint

Carlos Yoong<fnref n="1" /> and Mathias Legrand<fnref n="1" />

¹Department of Mechanical Engineering, McGill University, Montréal, Canada

**Abstract:** The present contribution describes a numerical technique devoted to the nonsmooth modal analysis (natural frequencies and mode shapes) of a non-internally resonant elastic bar of length $L$ subject to a Robin condition at $x = 0$ and a frictionless unilateral contact condition at $x = L$. When contact is ignored, the system of interest exhibits non-commensurate linear natural frequencies, which is a critical feature in this study. The nonsmooth modes of vibration are defined as one-parameter continuous families of nonsmooth periodic orbits satisfying the local equation together with the boundary conditions. In order to find a few of the above families, the unknown displacement is first expressed using the well-known d'Alembert's solution incorporating the Robin boundary condition at $x = 0$. The unilateral contact constraint at $x = L$ is reduced to a conditional switch between Neumann (open gap) and Dirichlet (closed gap) boundary conditions. Finally, $T$-periodicity is enforced. It is also assumed that only one contact switch occurs every period. The above system of equations is numerically solved for through a simultaneous discretization of the space and time domains, which yields a set of equations and inequations in terms of discrete displacements and velocities. The proposed approach is non-dispersive, non-dissipative and accurately captures the propagation of waves with discontinuous fronts, which is essential for the computation of periodic motions in this study. Results indicate that in contrast to its linear counterpart (bar without contact constraints) where modal motions are sinusoidal functions “uncoupled” in space and time, the system of interest features nonsmooth periodic displacements that are intricate piecewise sinusoidal functions in space and time. Moreover, the corresponding frequency-energy “nonlinear” spectrum shows backbone curves of the hardening type. It is also shown that nonsmooth modal analysis is capable of efficiently predicting vibratory resonances when the system is periodically forced. The pre-stressed and initially grazing bar configurations are also briefly discussed.

**Keywords:** nonsmooth systems, modal analysis, internal resonance, unilateral contact constraints, wave equation

## 1 Introduction

The concept of linear modes (natural frequencies and mode shapes) is a widely studied subject in the field of structural dynamics [7]. A possible extension of this notion to nonlinear conservative systems sees a mode of vibration as a one-parameter continuous family of periodic orbits displaying similar qualitative features [5]. In the phase space, nonlinear modes emerge as invariant surfaces of periodic trajectories, referred to as *invariant manifolds* [10], where invariant implies that the motion initiated on the manifold stays on it as time unfolds. To some extent, nonlinear modal analysis can be employed for predicting vibratory resonances, computing the nonlinear spectra of vibration or performing model-order reduction. Techniques traditionally employed for nonlinear modal analysis require a certain degree of smoothness in the nonlinearities [12] and thus fail for systems with *nonsmooth nonlinearities* such as unilateral contact constraints. Certainly, an accurate characterization of the vibratory response of these systems is essential to achieving enhanced and safer engineering applications [11]. Modal analysis of nonsmooth mechanical systems, also called *nonsmooth modal analysis*, has been recently proposed for a finite elastic bar of length $L$ subject to a Dirichlet boundary condition at $x = 0$ and a unilateral contact constraint at $x = L$ [13]. This system satisfies a *complete internal resonance* condition, i.e. all linear natural frequencies are commensurate with the first one, which has drastic consequences on the nonlinear modal response. Despite the simplicity of the system, the computed *nonsmooth modes* (NSMs) indicate highly intricate vibratory behaviour. Corresponding periodic displacements were observed to be unseparated piecewise linear functions of space and time, as opposed to their linear counterparts which are sinusoidal functions separated in space and time. Moreover, for certain NSMs such internal resonance generates a discontinuity between the linear and nonlinear portions of the invariant manifold. To further explore the nonlinear dynamics of this one-dimensional contact problem, a non-internally resonant configuration is investigated in the present work. The complete internal resonance condition is annihilated by changing the boundary condition (BC) at $x = 0$ from Dirichlet type $u(0, t) = 0$ to a Robin type $\partial_x u(0, t) - \alpha u(0, t) = 0$ which reflects that the elastic bar is now connected to a rigid support through a simple linear spring<fnref n="1" />. Analytical derivations are first proposed to facilitate the construction of the sought NSMs. Then, a numerical scheme based on the simultaneous discretization of the space and time domains is employed and the nonsmooth modes of vibration are constructed.

<footnote n="1">¹In this document, operators $\partial_{\xi}(\bullet)$ and $\partial_{\xi\xi}^2(\bullet)$ respectively stand for the first and second derivatives of $(\bullet)$ with respect to the argument $\xi$.</footnote>
## 2 Non-internally resonant elastic bar

The system of interest is a homogeneous elastic bar of length $L > 0$ and constant cross-sectional area $S > 0$ subject to a conservative unilateral constraint at its right end. Its left end is connected to a rigid support through a spring of stiffness $\kappa > 0$, as depicted in Fig. 1. The displacement, velocity, strain and stress fields are denoted by $u(x, t)$,

<figure><img src="image_2.png" /><figcaption>Fig. 1: One-dimensional finite elastic bar attached to a spring at its left extremity and subject to unilateral contact constraint on its right tip.</figcaption></figure>

$v(x, t)$, $\epsilon(x, t)$ and $\sigma(x, t)$ respectively where $x$ is the coordinate of a point of the bar in the initial configuration and $t$ denotes time. Young's modulus is denoted by $E > 0$ and $\rho > 0$ stands for the mass per unit volume, which are both, by assumption, space and time independent. Hence, the propagation speed of any longitudinal wave is $\sqrt{E/\rho}$. In the framework of linear elasticity, the stresses read $\sigma = E\epsilon$, where $\epsilon = \partial_x u$ should be infinitesimal and satisfy $|\partial_x u| < 1$ [1, p. 5]. The unilateral contact force $r(t)$ is related to the stresses by $\sigma(L, t) = E\partial_x u(L, t) = r(t)/S$. Further, the signed distance between the right extremity of the bar and the rigid obstacle, termed gap function, is defined as $g(u(L, t)) = g_0 - u(L, t)$, where $g_0$ is the signed distance between the unrestricted resting configuration and the obstacle: it is strictly negative in the pre-stressed configuration, for instance. Unless stated otherwise, there is no external excitation on the system. The full formulation reads:

$$ \text{Local equation} \qquad \partial_t^2 u(x, t) - c^2 \partial_{xx}^2 u(x, t) = 0, \qquad \forall x \in ]0; L[, \quad \forall t > 0, \quad (1) $$

$$ \text{Robin BC} \qquad \partial_x u(0, t) - \alpha u(0, t) = 0, \qquad \forall t > 0, \quad (2) $$

$$ \text{Signorini BC} \qquad g(u(L, t)) \ge 0, \quad r(t) \le 0, \quad r(t)g(u(L, t)) = 0, \qquad \forall t > 0, \quad (3) $$

$$ \text{Initial conditions} \qquad u(x, 0) = u_0(x), \quad v(x, 0) = v_0(x), \qquad \forall x \in ]0; L[. \quad (4) $$

where $\alpha = \kappa/(ES)$. This formulation possesses a unique solution which conserves the total energy [9]. It is worth noting that the local equation (1) is the well-known wave equation (a hyperbolic partial differential equation) defined on a one-dimensional finite domain. The natural frequencies $\omega_k$ and $\Omega_k$ of the underlying linear systems are solutions to the transcendental equations:

$$ \text{spring-free BCs} \qquad \omega_k - \alpha c \cot(\omega_k L/c) = 0, \qquad k \in \mathbb{N}_{>0}, \quad (5) $$

$$ \text{spring-fixed BCs} \qquad \Omega_k + \alpha c \tan(\Omega_k L/c) = 0, \qquad k \in \mathbb{N}_{>0}. \quad (6) $$

The corresponding natural periods are $T_k = 2\pi/\omega_k$ and $P_k = 2\pi/\Omega_k$. In both configurations, the natural frequencies are incommensurate, in the sense that $\omega_k$ and $\Omega_k$ for $k = 2, 3, \dots, \infty$ are not multiples of $\omega_1$ nor $\Omega_1$, respectively [8, p. 245]. Accordingly, the complete internal resonance condition emerging when the bar is clamped at $x = 0$ no longer holds [13].

Non-trivial solutions of the unilateral contact problem described by Eqs. (1) to (4) are successions of free phases (open gap) and contact phases (closed gap) [2]. Hence, these solutions can be perceived as the combination of motions satisfying the wave equation together with a switching boundary condition at $x = L$ between $\partial_x u(L, \cdot) = 0$ when the gap is open, referred to as “spring-free BCs” (or equivalently Robin-Neumann BCs), and prescribed displacement $u(L, \cdot) = g_0$ which implies $v(L, \cdot) = 0$ when the gap is closed, named “spring-fixed BCs” (or equivalently Robin-Dirichlet BCs). The nonlinearity in the formulation arises in the dependence of the solution to the unknown switching time. To further elaborate on this statement, consider the general solution to the local equation (1) comprising the superposition of forward $h$ and backward $f$ travelling waves defined on the real line $\mathbb{R}$ [1, p. 91]

$$ u(x, t) = f(ct + x) + h(ct - x). \quad (7) $$

At $x = 0$, the effect of the attached spring on the reflection of an incident backward wave $f$ is computed by inserting the general solution (7) in (2), yielding $(\partial_\xi f(\xi) - \partial_\xi h(\xi)) = \alpha(f(\xi) + h(\xi))$, $\forall \xi \in \mathbb{R}$, which in turns leads to the identity

$$ h(\xi) = f(\xi) + e^{-\alpha\xi} \left( \beta - 2\alpha \int_{0}^{\xi} e^{\alpha s} f(s) \, ds \right), \quad \forall \xi \in \mathbb{R} \quad (8) $$
Without going into the details, it is straightforward to check that the BC at $x = L$ implies $\beta = 0$. Incorporating (8) in (7) leads to an integral expression of the displacement in terms of $f$ only

$$u(x,t) = f(ct + x) + f(ct - x) - 2\alpha e^{\alpha(x-ct)} \int_{0}^{ct-x} e^{\alpha s} f(s) ds, \quad (9)$$

which incorporates the reflection mechanism of the travelling waves at $x = 0$. At $x = L$, the boundary condition successively switches from Neumann to Dirichlet and vice-versa. The corresponding reflections require additional considerations elaborated later in the text.

Let us first denote by $f|_{[a;b]}$ the set of all ordered pairs $(\xi, f(\xi))$, $\forall \xi \in [a; b]$, which is also referred to as the graph of function $f$ over the interval $[a; b]$. Provided that $f|_{[0;L]}$ is known, the boundary condition at $x = L$ is utilized to specify $f(\xi)$ over the remaining portions of the real axis $\mathbb{R}$, that is for $\xi \in ]-\infty; 0[ \text{ and } \xi \in ]L; \infty[$. During free phases (open gap), a homogeneous Neumann BC is active at $x = L$, that is $\partial_x u(L, \cdot) = 0$. Inserting Eq. (7) into the latter yields $\partial_{\xi} f(\xi + 2L) = \partial_{\xi} h(\xi)$, which in turn results in $f(\xi + 2L) - h(\xi) = f(2L) - h(0)$. Incorporating the influence of the attached spring via Eq. (8), the latter is re-formulated in terms of the backward wave $f$ only, that is

$$f(\xi + 2L) - f(\xi) + 2\alpha e^{-\alpha\xi} \int_{0}^{\xi} e^{\alpha s} f(s) ds = f(2L) - f(0). \quad (10)$$

During contact phases (closed gap), the right end of the bar satisfies a non-homogeneous Dirichlet BC, corresponding to a prescribed displacement at $x = L$, which reads $u(L, \cdot) = g_0$. Again, combining this boundary condition and the general solution (7), the interaction of forward and backward travelling waves shall satisfy $f(\xi + 2L) + h(\xi) = g_0$, which can also be expressed in terms of $f$ only as follows:

$$f(\xi + 2L) + f(\xi) - 2\alpha e^{-\alpha\xi} \int_{0}^{\xi} e^{\alpha s} f(s) ds = g_0, \quad (11)$$

Accordingly, if $f$ is known in any region of length $2L$, it can be expanded over the real axis, via (10) for a free phase (open gap) and (11) for a contact phase (closed gap). It is then employed in (9) to calculate the corresponding displacement.

During each of the above phases, $f|_{[0;2L]}$ is obtained from the associated initial and boundary conditions [4, p. 80]. For the free phase, knowing that the displacement and velocity waves reflect at a free BC without changing signs, $f|_{[0;2L]}$ satisfies

$$f|_{[0;2L]} = \begin{cases} \frac{u_0(\xi)}{2} + \frac{1}{2c} \int_0^\xi v_0(s) ds & \xi \in [0; L], \\ \frac{u_0(2L - \xi)}{2} + \frac{1}{2c} \int_0^L v_0(s) ds + \frac{1}{2c} \int_L^\xi v_0(2L - s) ds & \xi \in ]L; 2L]. \end{cases} \quad (12)$$

However, for the contact phase, knowing that $u_0(L) = g_0$ and that the displacement and velocity waves reflect with opposite sign, $f|_{[0;2L]}$ shall satisfy

$$f|_{[0;2L]} = \begin{cases} \frac{u_0(\xi)}{2} + \frac{1}{2c} \int_0^\xi v_0(s) ds & \xi \in [0; L], \\ g_0 - \frac{u_0(2L - \xi)}{2} + \frac{1}{2c} \int_0^L v_0(s) ds - \frac{1}{2c} \int_L^\xi v_0(2L - s) ds & \xi \in ]L; 2L], \end{cases} \quad (13)$$

Consequently, $f$ can be completely defined over $\mathbb{R}$ for every phase.

To summarize, the displacement of the elastic bar with Robin BC on the left end is obtained through Eq. (9) provided that $f$ is defined everywhere on the real axis. The successive switches in boundary conditions at $x = L$, reflecting the unilateral contact constraint, are incorporated through appropriate extensions: Equation (10) for the free phase or Equation (11) for contact phase. The following section is concerned with analytical derivations for the computation of periodic solutions by employing the expressions (9), (10) and (11).

### 3 Periodic solutions and Nonsmooth Modal Analysis

Nonsmooth modes of vibration of the spring-bar system depicted in Fig. 1 are defined as continuous families of periodic solutions satisfying the formulation (1)-(4) together with periodicity conditions in displacement and velocity: $\exists T > 0$
such that $u(x, t + T) = u(x, t)$ and $v(x, t + T) = v(x, t)$, $\forall x \in [0; L]$ and $\forall t > 0$. Finding such solutions translates into finding corresponding initial conditions $u_0$ and $v_0$ and period $T$ which generate periodic motions. Without loss of generality, it is assumed that within one period, over the interval $t \in [0; T]$, the initial time segment is a free phase that initiates at $t = 0^+$ and the final time segment is a contact phase that ends at $t = T^-$ and switches back to the initial free phase state at $t = T^+$. In general, various successions of free and contact phases might arise within one period. Knowing that the motion of the bar can be uniquely defined by a single function $f$, the targeted trajectory is then an unknown periodic sequence of functions $u$ in the form (9) where $f$ switches between (10) and (11).

Let us consider the simplest combination of one free phase and one contact phase of duration $t_f$ and $t_c = T - t_f$, respectively. When the gap is open, displacement and velocity satisfy the following equalities $\forall x \in [0; L]$ and $\forall t \in [0; t_f]$:

$$u_1(x, t) = f_0(ct + x) + f_0(ct - x) - 2\alpha e^{\alpha(x-ct)} \int_0^{ct-x} e^{\alpha s} f_0(s) ds, \quad (14a)$$

$$\frac{1}{c}v_1(x, t) = \partial_t f_0(ct + x) + \partial_t f_0(ct - x) - 2\alpha f_0(ct - x) + 2\alpha^2 e^{\alpha(x-ct)} \int_0^{ct-x} e^{\alpha s} f_0(s) ds, \quad (14b)$$

where $f_0|_{[0;2L]}$, calculated via Eq. (12) with the initial conditions $u_0$ and $v_0$, is then expanded throughout the real axis $\mathbb{R}$ via (10). During gap closure, the motion is described by the composite function $f_1(f_0(\cdot))$ corresponding to a boundary condition switch. This function arises by defining the graph $f_1|_{[0;2L]}$ with $u_1(\cdot, t_f)$ and $v_1(\cdot, t_f)$ as the “initial conditions” in Eq. (13) and then expanding it on $\mathbb{R}$ via Eq. (11).

During a contact phase, the displacement and velocity read $\forall x \in [0; L]$ and $\forall t \in ]t_f; T]$

$$u_2(x, t) = f_1(c(t - t_f) + x) + f_1(c(t - t_f) - x) - 2\alpha e^{\alpha(x-c(t-t_f))} \int_0^{c(t-t_f)-x} e^{\alpha s} f_1(s) ds, \quad (15a)$$

$$\frac{1}{c}v_2(x, t) = \partial_t f_1(c(t - t_f) + x) + \partial_t f_1(c(t - t_f) - x) \\ - 2\alpha f_1(c(t - t_f) - x) + 2\alpha^2 e^{\alpha(x-c(t-t_f))} \int_0^{c(t-t_f)-x} e^{\alpha s} f_1(s) ds. \quad (15b)$$

Accordingly, admissible $T$-periodic motions involving one lasting contact phase per period are described by functions $f_1 \circ f_0$ and $f_0$ that satisfy the following periodicity condition, arising from $u_0(x) = u_2(x, T)$ and $v_0(x) = v_2(x, T)$, $\forall x \in [0; L]$

$$u_0(x) = f_1(ct_c + x) + f_1(ct_c - x) - 2\alpha e^{\alpha(x-ct_c)} \int_0^{ct_c-x} e^{\alpha s} f_1(s) ds, \quad (16a)$$

$$\frac{1}{c}v_0(x) = \partial_t f_1(ct_c + x) + \partial_t f_1(ct_c - x) - 2\alpha f_1(ct_c - x) + 2\alpha^2 e^{\alpha(x-ct_c)} \int_0^{ct_c-x} e^{\alpha s} f_1(s) ds. \quad (16b)$$

together with the admissibility conditions reflecting Signorini's conditions at $x = L$ in Eq. (3):

* the elastic bar shall not contact the obstacle during a free phase except at gap closure when $t = t_f$; grazing is permitted during free phase

$$g_0 - f_0(ct + L) - f_0(ct - L) + 2\alpha e^{\alpha(L-ct)} \int_0^{ct-L} e^{\alpha s} f_0(s) ds \geq 0, \quad \forall t \in [0; t_f], \quad (17a)$$

$$g_0 - f_0(ct_f + L) - f_0(ct_f - L) + 2\alpha e^{\alpha(L-ct_f)} \int_0^{ct_f-L} e^{\alpha s} f_0(s) ds = 0, \quad (17b)$$

* the contact force should be non-positive during gap closure until contact separation at $t = t_f + t_c = T$

$$\partial_x f_1(ct + L) - \partial_x f_1(ct - L) + 2\alpha f_1(ct - L) - 2\alpha^2 e^{\alpha(L-ct)} \int_0^{ct-L} e^{\alpha s} f_1(s) ds \leq 0, \quad \forall t \in [0; t_c], \quad (18a)$$

$$\partial_x f_1(ct_c + L) - \partial_x f_1(ct_c - L) + 2\alpha f_1(ct_c - L) - 2\alpha^2 e^{\alpha(L-ct_c)} \int_0^{ct_c-L} e^{\alpha s} f_1(s) ds = 0. \quad (18b)$$

From Eq. (16), the equality $c\partial_x u_0(x) + v_0(x) = c\partial_x u_2(x, T) + v_2(x, T)$ provides an additional relationship between $f_0$ and $f_1$ for potential periodic motions to exist: $f_1(ct_c + x) - f_1(ct_c) = f_0(x) - f_0(0)$, $\forall x \in [0; L]$.
Finding solutions to Eqs. (16)-(18) is a noticeably challenging task. The complexity for solving such problems lies in the difficulty of defining a simple relationship between functions $f_0$ and $f_1$ incorporating the BC switching mechanism. Additionally, the effect of the spring, arising as a delay integral term in the displacement (9) further complicates the identification of admissible periodic motions.

## 4 Numerical Scheme

The main objective of the proposed numerical scheme is to simultaneously discretize the space and time domains of the above formulation in order to accurately mimic the propagation of discontinuous waves along the characteristics lines: $x \pm ct$ = constant. The main limitation of the proposed technique is the fact that the *travelling-wave solution* shall be partially known (in the sense that identities such as Eqs. (16)-(18) can be derived) before discretization.

To compute families of periodic orbits, the space and time domains of the integral equations (16) are simultaneously discretized in order to approximate the initial conditions that generate a periodic motion. Space is divided into $N$ intervals of identical length $\Delta x = x_{i+1} - x_i = L/N$ with $i = 0, 1, \dots, N$. Since travelling waves are required to propagate along the characteristic lines, the time-step shall satisfy $\Delta t = t_{n+1} - t_n = \Delta x/c$ for $n = 0, 1, \dots, n_T$, where $n_T$ satisfies $n_T \Delta t = T$. To approximate $f_1$ in Eq. (16), the discretization of the real axis emerges from the discretized space-time coupling argument $x \pm ct$ of the travelling-wave solution, hence the discretized $f_1$ function is denoted as $f_{1j} \approx f_1(\xi_j)$ such that $\xi_j = x_i \pm ct_n$ where $j = i \pm n$ for $i = 0, 1, \dots, N$ and $n = 0, 1, \dots, n_T$. Similar notations also apply to $f_{0j}$. The discretization of the displacement and velocity during free phase produces $u_{1i}^{(n)} \approx u_1(x_i, t_n)$ and $v_{1i}^{(n)} \approx v_1(x_i, t_n)$ respectively, and in a similar fashion during contact phase for $u_2$ and $v_2$. The discretized initial conditions read $u_{0i} \approx u_0(x_i)$ and $v_{0i} \approx v_0(x_i)$.

Since the duration of the free phase is $t_f = n_f \Delta t$ and the duration of the contact phase is $t_c = n_c \Delta t$, the approximations of $f_0|_{[0;2L]}$ and $f_1|_{[0;2L]}$ needed to approximately solve (16) are computed via a trapezoidal rule to compute the integrals:

* for a free phase via discretization of Eq. (12) with discrete initial conditions $u_{0i} \approx u_0(x_i)$ and $v_{0i} \approx v_0(x_i)$

$$f_{0j} = \begin{cases} \frac{u_{0j}}{2} + \frac{1}{4c} \sum_{i=0}^{j} \Delta x(v_{0i} + v_{0(i+1)}) & j = 0, 1, \dots, N, \\ a + \frac{u_{0(2N-j)}}{2} + \frac{1}{4c} \sum_{i=N+1}^{j} \Delta x(v_{0(2N-i)} + v_{0(2N-i+1)}) & j = N+2, N+3, \dots, 2N, \end{cases} \quad (19)$$

* for a contact phase via discretization of Eq. (13) with discrete “initial conditions” $u_{1i}^{(n_f)} \approx u_1(x_i, t_f)$ and $v_{1i}^{(n_f)} \approx v_1(x_i, t_f)$

$$f_{1j} = \begin{cases} \frac{u_{1j}^{(n_f)}}{2} + \frac{1}{4c} \sum_{i=0}^{j} \Delta x(v_{1i}^{(n_f)} + v_{1(i+1)}^{(n_f)}) & j = 0, 1, \dots, N, \\ b - \frac{u_{1(2N-j)}^{(n_f)}}{2} - \frac{1}{4c} \sum_{i=N+1}^{j} \Delta x(v_{1(2N-i)}^{(n_f)} + v_{1(2N-i+1)}^{(n_f)}) & j = N+2, N+3, \dots, 2N, \end{cases} \quad (20)$$

where $4ac = \sum_{i=0}^{N} \Delta x(v_{0i} + v_{0(i+1)})$ and $4bc = \sum_{i=0}^{N} \Delta x(v_{1i}^{(n_f)} + v_{1(i+1)}^{(n_f)}) + 4g_0c$. The above expressions (19) and (20) known for $j = 0, 1, \dots, 2N$ are then expanded on the real axis:

$$f_{0(j+2N)} - f_{0(j)} + \operatorname{sgn}(j) \alpha e^{-\alpha j} \Delta x \sum_{k=0}^{j} (e^{\alpha \xi_k} f_{0(k)} + e^{\alpha \xi_{k+1}} f_{0(k+1)}) = f_{0(2N)} - f_{0(0)}, \quad (21)$$

$$f_{1(j+2N)} + f_{1(j)} + \operatorname{sgn}(j) \alpha e^{-\alpha j} \Delta x \sum_{k=0}^{j} (e^{\alpha \xi_k} f_{1(k)} + e^{\alpha \xi_{k+1}} f_{1(k+1)}) = g_0, \quad (22)$$

where $\operatorname{sgn}(\cdot)$ is the sign function. Since the period $T = t_f + t_c = n_T \Delta t$ is known, where $n_T = n_f + n_c$, the approximated displacement can be calculated for each phase as follows:

* free phase, via discretization of Eq. (14), for $i = 0, 1, \dots, N$ and $n = 0, 1, 2, \dots, n_f$

$$u_{1i}^{(n)} = f_{0(n+i)} + f_{0(n-i)} - \operatorname{sgn}(n-i) \alpha e^{\alpha(x_i - ct_n)} \sum_{k=0}^{n-i} \Delta x (e^{\alpha \xi_k} f_{0(k)} + e^{\alpha \xi_{k+1}} f_{0(k+1)}), \quad (23)$$
• contact phase, via discretization of Eq. (15), for $i = 0, 1, \dots, N$ and $n = n_f + 1, n_f + 2, \dots, n_T$

$$u_{2i}^{(n)} = f_{1(n-n_i+i)} + f_{1(n-n_i-i)} - \operatorname{sgn}(n-n_f-i)\alpha e^{\alpha(x_i-c(t_n-t_i))} \sum_{k=0}^{n-n_f-i} \Delta x (e^{\alpha \xi_k} f_{1(k)} + e^{\alpha \xi_{k+1}} f_{1(k+1)}). \quad (24)$$

The corresponding velocity is calculated through a numerical time-differentiation scheme, such as the forward Euler method $v_i^{(n)} = (u_i^{(n+1)} - u_i^{(n)})/\Delta t$. It is worth remarking that because $f_{1j}$ in Eq. (24) is compounded with $f_{0j}$, and the latter is computed from discrete initial conditions, then all displacements points $u_{1i}^{(n)}$ and $u_{2i}^{(n)}$ can be expressed in terms of $u_{0i}$ and $v_{0i}$.

The discretization of the periodicity condition (16) together with the unilateral contact condition yields a constrained system of linear equations written in terms of the discrete initial conditions $u_{0i}$ and $v_{0i}$, free-phase time-steps $n_f$ and contact-phase time-steps $n_c$:

$$\text{periodicity:} \qquad u_{0i} - u_{2i}^{(n_f+n_c)} = 0 \text{ and } v_{0i} - v_{2i}^{(n_f+n_c)} = 0, \qquad i = 0, 1, \dots, N \quad (25)$$

$$\text{impenetrability:} \qquad g_0 - u_{1N}^{(n)} \ge 0, \qquad n = 0, 2, \dots, n_f \quad (26)$$

$$\text{compressive contact:} \qquad E \epsilon_{2N}^{(n)} \le 0, \qquad n = n_f + 1, n_f + 2, \dots, n_f + n_c \quad (27)$$

where $\epsilon_{2N}^{(n)}$ is the strain at $x = L$ at time $t_n$ calculated from the numerical space-differentiation of discrete displacements $u_{2i}^{(n)}$; for instance, the implementation of the forward Euler scheme yields discrete strains in the form of $\epsilon_{2i}^{(n)} = (u_{2i}^{(n)} - u_{2(i-1)}^{(n)})/\Delta x$. The discretized formulation given by Eqs. (25)-(27) is solved through the following steps:

1. Set numbers of time-steps $n_f$ and $n_c$.

2. Calculate potential initial conditions $u_{0i}$ and $v_{0i}$ that generate a periodic motion satisfying Eq. (25) only.

3. Verify that the generated periodic motion satisfies the conditions of impenetrability (26) and compressive contact (27): the corresponding initial conditions are admissible.

4. If the initial conditions are admissible, compute and store the period of vibration and total energy.

5. Change values of $n_f$ and $n_c$.

For the identification of families of periodic orbits, every feasible pair $(n_f, n_c) \in \mathbb{N}_{>0}^2$ is examined. Nevertheless, verifying among all potential combinations is computationally very inefficient. For that reason, the iteration intervals are bounded for each family of admissible periodic motions by analyzing every possible combination of $n_f$ and $n_c$ on coarse meshes.

# 5 Spectrum of nonsmooth vibration

Three initial configurations of the bar are explored: unstressed ($g_0 > 0$), prestressed ($g_0 < 0$) and initially grazing ($g_0 = 0$). Without loss of generality, the mechanical parameters $E, \rho, S$ and $L$ are arbitrarily chosen to be unity and units are discarded. The autonomous dynamics of the bar is investigated for two spring-bar stiffness ratios $\alpha$. The linear natural frequencies $\omega_k(\alpha)$ and $\Omega_k(\alpha)$ respectively calculated using Eqs. (5) and (6), now only depend on the parameter $\alpha$. The proposed discretization technique with $\Delta t = \Delta x/c = 10^{-3}$ is implemented. Modal responses were obtained with $g_0 = \pm 10^{-3}$ in both unstressed and prestressed configurations. The response of the autonomous elastic bar is depicted in frequency-energy plots (FEPs) where a backbone curve (also known as branch) represents a NSM. Each point pertaining to a backbone curve represents a modal motion whose frequency is indicated on the horizontal axis and constant total energy along the vertical axis. The frequencies in the FEPs are not normalized while the energy is normalized with respect to the energy of the first linear mode grazing orbit.

NSMs computed for $\alpha = 1$ and $\alpha = 1/2$ are now investigated. In contrast to the results exposed in [13], the intricate dynamics caused by the spring complicates the identification of NSM branches and a highly fine discretization is required to discern how admissible periodic solutions organize to form a continuum. The backbone curves, emerging in the frequency range [$\omega_1 : \Omega_1$] are depicted in Fig. 2 as sets of scattered points supposedly belonging to NSMs. Contrary to main NSMs of the internally resonant bar studied in [13] where the energy continuously depends on the frequency, the depicted scattered points indicate more complicated backbone curves. Figure 2 suggests that a “nicely” connected continuum representing a possible “main NSM” does not exist. However, several independent branches emerge around subharmonics frequencies lying in the range for the unstressed configuration. This agrees with NSMs in discrete systems without linear internal resonance [6, 11]. An equivalent phenomenon occurs for the prestressed configuration where independent branches align in a seemingly softening branch. For the initially grazing configuration, the NSMs emerge as vertical lines embedding periodic solutions with identical frequency and increasing energy. It
<figure><img src="image_3.png" /><figcaption>**Fig. 2:** Backbone curves in the frequency interval $[\omega_1; \Omega_1]$ for $\alpha = 1$ and $\alpha = 1/2$ with bar configurations: unstressed [—], prestressed [—] and initially grazing [—]. Subharmonics frequencies [---].</figcaption></figure>

also appears that every given collection of scattered points represents a purely independent subharmonic vibration, as opposed to the internally resonant bar which features subharmonic as well as internally resonant NSMs [13]. This observation coincides with the fact that the system of interest does not satisfy the full internal resonance condition.

Figures 3 and 4 depict NSM periodic displacement fields, for $\alpha = 1$ and $\alpha = 1/2$ respectively, corresponding to points **b** and **c** in Fig. 2. Both points **b** represent periodic motions belonging to backbone curves that emerge

<figure><img src="image_4.png" /><figcaption>**Fig. 3:** Periodic displacement fields for $\alpha = 1$ corresponding to points **b** (left) and **c** (right) in Fig. 2.</figcaption></figure>

in the vicinity of the respective $\omega_1$. On the other hand, points **c** live in apparent NSM branches that arise in the neighborhood of linear subharmonic frequencies. The complicated pattern of each solution is caused by an intricate interplay between various travelling waves embedding the Robin and Signorini boundary conditions. In contrast to linear modes that are purely harmonic functions, the nonsmooth modes of the non-internally resonant system are nonsmooth piecewise-sinusoidal functions. The velocity fields shall present several jump discontinuities per period, which would not be accurately described by traditional semi-discretization strategies [3]. From the reported motions, it also seems that the displacement field for $t \in [0; t_f]$ presents an axis of symmetry in time located at the middle of the free phase: $t = t_f/2$. This observation could facilitate the derivation of closed-form expressions; this is left to future investigations.
<figure><img src="image_5.png" /><figcaption>Fig. 4: Periodic displacement fields for $\alpha = 1/2$ corresponding to points (b) (left) and (c) (right) in Fig. 2.</figcaption></figure>

## 6 Response to periodic external forcing

This section numerically investigates the relationships between NSMs and the system response under periodic excitation. A slight amount of structural viscous damping is introduced in the governing equation (1) and a weak external damper is attached at $x = 0$ so that the system can possibly reach a periodic steady-state with bounded energy. The system is forced via a harmonically moving rigid wall that periodically compresses the bar, as seen in Fig. 5. The formulation for

<figure><img src="image_6.png" /><figcaption>Fig. 5: Spring–bar system excited by a moving rigid wall.</figcaption></figure>

the forced spring–bar system, to be compared to Eqs. (1)-(4), reads

$$ \partial_t^2 u(x, t) + \zeta_1 \partial_t u(x, t) - c^2 \partial_{xx}^2 u(x, t) = 0, \quad \forall x \in ]0; L[, \quad \forall t > 0, \quad (28) $$

$$ \partial_x u(0, t) - \alpha u(0, t) - \zeta_2 \partial_t u(0, t) = 0, \quad \forall t > 0, \quad (29) $$

$$ g(u(L, t), w(t)) \ge 0, \quad r(t) \le 0, \quad r(t)g(u(L, t), w(t)) = 0, \quad \forall t > 0, \quad (30) $$

$$ u(x, 0) = u_0(x), \quad v(x, 0) = v_0(x), \quad \forall x \in ]0; L[, \quad (31) $$

where $\zeta_1$ and $\zeta_2$ are respectively the structural and external damping coefficients. The moving wall excitation corresponds to $w(t) = w_0 \sin \omega t$. For an elastic bar initially at rest, the response of the periodically-forced system is constructed by computing, when possible, the steady-state solution for each frequency of excitation. In this study, such solutions are obtained via the Wave Finite Element Method [13]. The total energy of the steady-state solution, averaged over one forcing period, for increasing frequencies of excitation and various damping coefficients is illustrated in Fig. 6 in the interval $[\omega_1; \Omega_1]$. In the frequency intervals where NSMs possibly do not exist$^2$, highlighted with grayed areas in Fig. 6, the forced system does not seem to exhibit periodic steady-states. When non-periodic forced responses are detected, corresponding portions of the response curves in the frequency diagram consider, for each frequency of excitation, the total energy of the forced motion averaged over ten periods of the external forcing. It is observed that forced responses with high-energy periodic steady-state align well with the NSM backbone curves. Also, Figure 7 illustrates a steady-state displacement of the slightly-damped bar under periodic forcing computed in the vicinity of an NSM. It clearly resembles the corresponding NSM motion, see Fig. 3. Although this forced solution seems to be identical to the autonomous oscillation, its free phase is actually not symmetric with respect to a time axis located in the middle of the free phase.

## 7 Conclusions

This contribution targeted the nonsmooth modal analysis of a non-internally resonant bar through a numerical strategy based on simultaneous space-time discretization of the travelling-wave solution. The system of interest consisted of a

<footnote>$^2$The number of computable NSMs depends on the mesh size, and finer meshes should provide a more detailed spectrum.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6:** Response to periodic forcing over $[\omega_1 : \Omega_1]$ of the spring–bar system under various damping coefficients: $\alpha = 1$ (top) and $\alpha = 1/2$ (bottom). Bar configuration at rest: unstressed $g_0 > 0$ (left), grazing $g_0 = 0$ (center) and prestressed $g_0 < 0$ (right). Bar damping coefficient $\zeta_1$: low [---] to high [---]. Grayed regions where NSMs were not found and forced-responses are not periodic.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7:** Displacement over one steady-state period due to external excitation via moving wall [---]. Result computed with $(\omega, \alpha) = (1.25, 1)$ and $g_0 > 0$. It corresponds to point (d) in Fig. 6 in the vicinity of NSM motion (b) in Fig. 2 top and Fig. 3 left.</figcaption></figure>

finite elastic bar of length $L$ subject to a Robin boundary condition at $x = 0$ and a unilateral contact constraint at $x = L$. Such configuration annihilates the full internal resonance condition featured by its internally resonant counterpart [13]. In addition, the Robin BC, physically corresponding to a simple spring attachment, causes a distortion in the waveform that complicates the numerical construction of nonsmooth modes. A semi-analytical approach is derived from the exact solution to the autonomous wave equation together with a switch in the boundary condition where unilateral contact arises. As nonsmooth periodic solutions in closed form are inaccessible, a discretization strategy is proposed to find families of periodic motions. The periodic nonsmooth motions are piecewise-sinusoidal functions, as opposed to the internally resonant counterparts where modal displacements are piecewise-linear functions. The forced response plots illustrated the capability of the Nonsmooth Modal Analysis to predict the vibratory resonances of the one-dimensional periodically-forced elastic bar. Even though the vibratory characterization is more challenging than for the internally resonant counterpart [13], the computed NSMs forecast most intervals involving nonlinear resonances.

## References

[1] Jan Achenbach. *Wave propagation in elastic solids*. Elsevier, 1973. DOI: 10.1016/C2009-0-08707-8.

[2] Luigi Amerio. Continuous solutions of the problem of a string vibrating against an obstacle. *Rendiconti del Seminario Matematico della Università di Padova*, 59:67–96, 1978. EUDML: 107691.

[3] David Doyen, Alexandre Ern, and Serge Piperno. Time-integration schemes for the finite element dynamic Signorini problem. *SIAM Journal of Scientific Computing*, 33(1):223–249, 2011. OAI: hal-00440128.
[4] Karl Graff. *Wave Motion in Elastic Solids*. Dover, 1975.

[5] Gaëtan Kerschen, Keith Worden, Alexander Vakakis, and Jean-Claude Golinval. Past, present and future of nonlinear system identification in structural dynamics. *Mechanical systems and signal processing*, 20(3):505–592, 2006. OAI: hal-01863233.

[6] Mathias Legrand, Stéphane Junca, and Sokly Heng. Nonsmooth modal analysis of a *n*-degree-of-freedom system undergoing a purely elastic impact law. *Communications in Nonlinear Science and Numerical Simulation*, 45:190–219, 2017. OAI: hal-01185980.

[7] Leonard Meirovitch. *Analytical Methods in Vibrations*. Macmillan New York, 1967.

[8] Singiresu Rao. *Vibration of Continuous Systems*. Wiley, 2007. DOI: 10.1002/9780470117866.

[9] Michelle Schatzman and Michel Bercovier. Numerical approximation of a wave equation with unilateral constraints. *Mathematics of Computation*, 53(187):55–79, 1989. OAI: hal-01295436.

[10] Steven Shaw and Christophe Pierre. Non-linear normal modes and invariant manifolds. *Journal of sound and Vibration*, 150(1):170–173, 1991. OAI: hal-01310674.

[11] Anders Thorin, Pierre Delezoide, and Mathias Legrand. Non-smooth modal analysis of piecewise-linear impact oscillators. *SIAM Journal on Applied Dynamical Systems*, 16(3):1710–1747, 2017. OAI: hal-01298983.

[12] Anders Thorin and Mathias Legrand. Nonsmooth modal analysis: From the discrete to the continuous settings. In Remco Leine, Vincent Acary, and Olivier Brüls, editors, *Advanced Topics in Nonsmooth Dynamics*, Transactions of the European Network for Nonsmooth Dynamics, chapter 5, pages 191–234. Springer, 2018. OAI: hal-01771849.

[13] Carlos Yoong, Anders Thorin, and Mathias Legrand. Nonsmooth modal analysis of an elastic bar subject to a unilateral contact constraint. *Nonlinear Dynamics*, 91(4):2453–2476, 2018. OAI: hal-01471341.