# Nonsmooth Modal Analysis of a Varying Cross-Sectional Area Bar in Unilateral Contact

David Urman, Mathias Legrand

► To cite this version:

David Urman, Mathias Legrand. Nonsmooth Modal Analysis of a Varying Cross-Sectional Area Bar in Unilateral Contact. Journal of Sound and Vibration, 2023, 545, pp.117385. 10.1016/j.jsv.2022.117385. hal-03619197v2

HAL Id: hal-03619197

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

Submitted on 12 Dec 2022

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 Varying Cross-Sectional Area Bar in Unilateral Contact

David Urman, Mathias Legrand

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

**Abstract:** Nonsmooth modes of vibration allow for identification of resonant behaviours and attendant vibratory frequencies in structures prone to unilateral contact conditions on the boundary. The prominent approach for finding nonsmooth modes of vibration entails finding continua of periodic solutions to the system in question. In this paper, nonsmooth modes of a one-dimensional bar of varying cross-sectional area prone to unilateral contact with a rigid obstacle are determined. While numerical and analytical techniques were previously proposed, they were limited to constant cross section bar and could not be applied on the varying area bar for which the classical d'Alembert solution no longer exists. In this article, nonsmooth modal analysis of the varying area bar is performed via a novel treatment of the Signorini conditions within the finite element framework: the *nodal boundary method*. The nodal boundary method solves the Signorini problem by switching between two sets of shape functions describing either (1) inactive contact motion (motion away from the rigid obstacle) or (2) active contact motion (bar in contact with the rigid obstacle). In the proposed nodal boundary method, the motion of the contacting node does not participate in the resulting governing Ordinary Differential Equation (ODE). Instead, its motion is prescribed by the boundary conditions and is dictated by the motion of internal nodes. The nodal boundary method results in a discontinuous ODE in the internal nodes which can be solved both analytically and via numerical techniques. Solutions obtained by the nodal boundary method exhibit several advantages over existing numerical techniques: no chattering at contact, no penetration of the rigid obstacle, and existence of periodic solutions. Specifically, these periodic solutions are readily detectable via the shooting method with sequential continuation. The nodal boundary method is used successfully for the nonsmooth modal analysis of different models of the varying area bar. Besides, application of the nodal boundary method for nonsmooth modal analysis of the uniform area cantilever bar in Dirichlet or Robin boundary conditions is also demonstrated for comparison with existing literature.

## Contents

1 Introduction

2 Problem statement

3 Finite-Element formulation

4 Nodal Boundary Method

4.1 Comment on application of NBM using other shape functions 6

4.2 Inactive contact motion 6

4.3 Active contact motion 8

4.4 NBM-FEM formulation of Signorini problem 9

4.5 Energy conservation properties of solutions to NBM-ODE 10

4.6 Notes on the NBM-FEM formulation 12

5 Nonsmooth modal analysis

5.1 Time-marching techniques and shooting method 13

5.2 Error estimation 13

5.3 Sequential continuation with correction 14

6 Results

6.1 Convergence of Crank-Nicolson and NBM 14

6.2 Nonsmooth modes 16

6.2.1 Forced-response curves 17

6.2.2 Constant cross-section and internally resonant bar 17

6.2.3 Constant cross-section bar with soft support 17

6.2.4 Varying-area bar 19

7 Scope and limitations of NBM

8 Conclusion

9 Supplementary Material

A Appendix

A.1 Proof of $\phi'_N(1) > 0$ 24

A.2 Conservation of energy away from instant of switch 25

A.3 Derivation of energy jump at switch 26
# 1 Introduction

For vibrations of structures prone to unilateral contact, families of periodic solutions are indicative of resonance frequencies and behaviours and are, in fact, referred to as nonsmooth modes of vibration (NSMs) [1–4]. In essence, NSMs are families of periodic motions of the unforced and undamped investigated structure [5–8]. In the present work, we focus on nonsmooth modal analysis of the bar of varying (cross-section) area in unilateral contact with a rigid structure. The dynamics of the bar are formulated via the wave equation with Signorini conditions for treatment of unilateral contact conditions [9, 10]. While periodic solutions have been obtained analytically for the case of the bar with constant area [1, 11], the same analytical method cannot be generalized to the case of the bar in varying area as a closed-form solution to the wave equation with varying coefficients has yet been obtained. Therefore, periodic solutions to the Signorini problem must be obtained numerically. Recent works [12, 13] showed that application of the Wave Finite Element method (WFEM) to the problem of the bar of constant area in unilateral contact produces energy-conserving schemes allowing for existence of periodic solutions. However, the conservation of energy in WFEM solutions does not persist for the case of the varying-area bar or two-dimensional Signorini problems [13]. Other than via the WFEM, the Signorini problem has been also solved via the finite element method (FEM) [9, 10, 14–16]. In the FEM, application the Signorini boundary conditions leads to an ill-posed problem if not supplied with an impact law or some relationship between the stress and displacement [14, p. 9]. Thus, multiple methods for implementation of Signorini contact have been conceived. Such schemes are those utilizing a Newtonian impact law relating the pre- and post-impact velocities of the contacting node in the form $v(t^+) = -ev(t^-)$ where $e \in [0, 1]$. However, these schemes generate non-physical chattering for energy conserving schemes ($e = 1$) and annihilate periodic solutions with non-zero pre-contact velocity when $e \neq 1$) [14]. Another common treatment of Signorini conditions is to approximate the contact-force using a displacement-dependent model which regularizes the non-penetration condition. For example, penalty methods use a non-linear and possibly smooth function for the contact force univoquely expressed in terms of the displacement [14, 17]. However, since the convergence of penalty methods depend on a penalty parameter reaching infinity, these method inevitably result in stiff ODE and are generally found insufficient in solving the Signorini problem accurately [13, 14]. Next, mass redistribution methods [15, 18, 19] aim to solve the issue of ill-posedness by eliminating the mass of nodes on the contact boundary such that at, every instance, prediction of the location the boundary nodes is obtained by solving the static Signorini problem [15]. Specifically, the mass redistribution method (MRM) proposed in [15] represents a formulation culminating in a an energy-conserving ODE which allows for periodic solutions while eliminating chattering at contact. However, the scheme requires solving a constrained optimization problem to redistribute the mass of the contact nodes. This issue was later readdressed in [18] where the usage of specific quadrature rules has reduced the computational cost of the mass redistribution. Regardless, the mass redistribution method requires usage of a different mass matrix than this used in classical FEM models. At last, it is worth mentioning the Nitsche method in FEM which solves the Signorini boundary condition in a weak sense [16]. While penetration is indeed allowed in the Nitsche method, it has been proven that the Signorini boundary conditions are satisfied as the number of elements increases [16].

In this work, we propose the *Nodal Boundary method* (NBM) for treatment of Signorini boundary conditions in FEM. Applied to the Signorini problem, this method allows for existence of periodic solutions while eliminating chattering. While, similarly to [15], this method requires redistribution of masses, it does not require solution of a constrained optimization problem nor usage of “non-traditional” quadrature rules (as done in [18]). Thus, the proposed method allows for implementation of mass redistribution to classical finite-element models with no significant addition to computational effort. At last, to find periodic solutions of the FEM-NBM reduced model, the shooting method together with the continuation method are used, similar to the procedure suggested in [7]. In particular, nonsmooth modal analysis of the varying area bar, which could not be approached via other techniques such as WFEM and Newton impact law, is performed in this article via the NBM.

In what follows, we detail the problem statement and its finite-element formulation in sections 2 and 3, respectively. Next, we present the nodal boundary method applied to the Signorini conditions in section 4. Application of nonsmooth modal analysis (ie, detection of periodic solutions) the model is presented in section 5. At last, the results section 6 consists of: convergence analysis, presentation of modes and forced motion analysis.
## 2 Problem statement

We consider a bar illustrated in Figure 1. The displacement field of the bar is denoted $u(\bar{x}, \bar{t})$, where $\bar{x}$ and $\bar{t}$ represent the physical position and time, respectively. The bar is clamped to a wall at $\bar{x} = 0$ and is prone to unilateral contact with a rigid obstacle at its other end $\bar{x} = L$, where L denotes the length of the bar. At rest, the rigid obstacle is set at a distance $\bar{g}$ from the bar.

<figure><img src="image_2.png" /><figcaption>Figure 1: Bar of varying area prone to unilateral contact with a rigid wall.</figcaption></figure>

Moreover, we consider the linear-elasticity framework for the deformation of the bar. Hence, the motion of the bar is described by the wave equation

$$u_{tt}(x,t) = (A(x)u_x(x,t))_x, \quad \forall x \in (0, 1), \quad t \in (0, \infty). \tag{1}$$

where subscripts denote partial differentiation with respect to the denoted variable. Furthermore, the following non-dimensional coordinates are introduced: $x = \bar{x}/L$ and $t = \alpha\bar{t}/L$. Here, $\alpha$ is introduced as a normalization factor $\alpha^2 = E\bar{A}_0/m$ where $E > 0$, $m > 0$ and $\bar{A}_0 = \bar{A}(0) > 0$ representing Young's modulus [N/m$^2$], the constant mass per unit length of the bar [kg/m], and the physical area of the bar at $x = 0$ [m$^2$], respectively. In turn, $A(x) > 0$ is a non-dimensional quantity representing the area variation in the bar and $A(0) = 1$ holds. The physical cross-sectional area of the bar abides $\bar{A}(\bar{x}) = \bar{A}(xL) = A(x)\bar{A}_0$. Moreover, we use the following Dirichlet boundary condition at $x = 0$:

$$u(0, t) = 0, \quad \forall t \in (0, \infty). \tag{2}$$

At $x = 1$ (ie, $\bar{x} = L$), we impose the Signorini boundary condition

$$0 \ge u(1, t) - g \perp \bar{\sigma}(t) \le 0 \quad \forall t \in (0, \infty), \quad \bar{\sigma}(t) = A(1)u_x(1, t) \tag{3}$$

to describe contact of the bar with the rigid wall, where the non-dimensional gap distance $g = \bar{g}/L$ and the stress at the contact end of the bar $\bar{\sigma}(t)$ are introduced. Namely, the Signorini conditions (3), in the continuous setting, can be seen as a switching of boundary conditions at $x = 1$:

$$\text{Active contact conditions: } u(1, t) = g \quad \text{and} \quad \bar{\sigma}(t) \le 0 \tag{4}$$

$$\text{Inactive contact conditions: } \bar{\sigma}(t) = 0 \quad \text{and} \quad u(1, t) \le g. \tag{5}$$

For purposes of modal analysis, we are required to find the initial conditions $u_{\text{init}}(x)$ and $v_{\text{init}}(x)$ generating periodic solutions

$$u(x, T) = u(x, 0) = u_{\text{init}}(x), \quad \forall x \in [0, 1] \tag{6}$$

$$u_t(x, T) = u_t(x, 0) = v_{\text{init}}(x), \quad \forall x \in [0, 1] \tag{7}$$

where $T$ denotes the period of motion. Accordingly, to solve the problem numerically, we use the FEM along with the nodal boundary method (NBM) for the treatment of boundary condition on the governing boundary value problem exposed from Equation (1) to Equation (5). Next, the shooting method and continuation are used to depict the continua of solutions answering Equations (6) and (7).

While the problem of modal analysis of the bar in unilateral contact with a constant cross-section, $A'(x) = 0$, has been studied before both analytically and numerically [1, 11, 12, 20], the same techniques could not be implemented for the case of $A'(x) \ne 0$. Namely, analytical techniques have relied on the exact solution to the wave equation to describe both inactive and active contact phases in a closed-form manner [1, 11, 20]. In contrast, the WFEM, in the case of the bar of uniform area, exhibits properties favoring the existence of periodic solutions such as: energy conservation and preservation of characteristic quantities [13]. However, the same
properties are not exhibited for the case for the varying area bar as the upwind-flux, used in the WFEM, does not accurately solve for the transfer of quantities between elements [21, chapter 9]. In contrast, it will be shown in this manuscript that the NBM allows for detection of periodic solutions in the varying area bar by implementation of a Galerkin-Bubnov method and boundary shape functions to satisfy the inactive and active phases of the Signorini conditions.

## 3 Finite-Element formulation

In order to apply FEM on the one-dimensional Signorini problem, the displacement within the bar is approximated by $u^h(x, t)$. In turn, the approximation $u^h(x, t)$ consists of a series of piecewise Lagrange polynomials $\phi_i(x)$, $i = 0, 1, 2, \dots, N$, and corresponding nodal quantities $u_i(t)$ [22] located at the nodes $x_i = i/N$ for $i = 0, 1, 2, \dots, N$, as classically achieved in FEM. The approximation thus reads

$$u(x,t) \approx u^h(x,t) = \sum_{i=0}^{N} \phi_i(x) u_i(t) \equiv \mathbf{P}(x)\mathbf{u}(t), \quad \text{where } \phi_i(x_j) = \delta_{ij} \text{ and } u_i(t) \approx u(x_i,t) \quad (8)$$

with $\mathbf{u}(t)$ storing the time-domain nodal displacements and $\delta_{ij}$ denoting the Kronecker Delta. Furthermore, we introduced the vector quantity $\mathbf{P}(x) = (\phi_0(x) \phi_1(x) \dots \phi_N(x))$ to simplify the representation of the finite-element in matrix form. The finite element applies to the weak form of partial differential equation (PDE) (1). It requires the definition of test functions $\mathbf{w}(t)$ corresponding to the nodal displacements $\mathbf{u}(t)$. In the NBM, the test functions $\mathbf{w}(t)$ will be subject to change through time (according to the phase of contact motion), while it is not the case for common representation in the FEM. This will be clarified in section 4. Hence, application of the finite element approximation on the weak form of PDE (1) with cantilever condition (2) ($u_0(t) = 0$ and such that $\phi_0(x)$ is omitted from $\mathbf{P}(x)$) yields

$$\mathbf{w}^{\top}(t)\mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{w}^{\top}(t)\mathbf{K}\mathbf{u}(t) - w_N\bar{\sigma}(t) = 0, \quad \forall \mathbf{w}(t) \quad (9)$$

where $\mathbf{M}$ is the mass matrix, $\mathbf{K}$ is the stiffness matrix, with respective entries

$$\mathbf{M} = \int_{0}^{1} \mathbf{P}(x)\mathbf{P}(x) \, dx, \quad \mathbf{K} = \int_{0}^{1} A(x)\mathbf{P}_{x}^{\top}(x)\mathbf{P}_{x}(x) \, dx. \quad (10)$$

While this notation is conventional in finite element (FE) analysis, it is reminded here since the NBM formulation will rely on it significantly. Specifically, the NBM modifies the weak form (9) throughout the motion and the test functions $\mathbf{w}(t)$ actively participate in its formulation. We now introduce the NBM to treat the Signorini conditions (4) and (5) in the finite element framework. There, the displacement and the stress at $x = 1$ will be approximated using $u^h(x, t)$

$$u(1,t) \approx u^h(1,t) = u_N(t), \quad \bar{\sigma}(t) \approx A(1)u_x^h(1,t) = A(1) \sum_{i=1}^{N} \phi_i'(1)u_i(t), \quad (11)$$

respectively. While the approximation $u(1,t) \approx u_N(t)$ is equivalent to the one used in the classical FEM, the above stress approximation $\bar{\sigma}(t) \approx A(1)u_x^h(1,t)$ is non-traditional. To clarify, in the classical FEM, implementation of Neumann conditions is done in a weak sense, and the motion of the nodes does not satisfy the Neumann condition for any grid. Instead, the error in answering the Neumann condition is reduced with increasing the number of elements or degree of polynomials. In the NBM, the Neumann conditions are imposed on the shape functions for any choice of number of elements or degree of polynomial such that the stress at the end of the bar is always strictly zero and the Neumann condition is exactly satisfied throughout inactive contact motion. In fact, this is the key to the implementation of the NBM, as described below.

## 4 Nodal Boundary Method

In NBM, the FE approximation is obtained by plugging the approximations (11) into the complementarity conditions (3) such that

$$g \ge u_N(t) \perp A(1) \sum_{i=1}^{N} \phi_i'(1)u_i(t) \le 0. \quad (12)$$
The main proposition in NBM is that condition (12) is solved by constructing shape functions capable of satisfying the inactive and active contact conditions. To do so, the boundary node $u_N(t)$ is isolated in (12) such that

$$g \geq u_N(t), \quad u_N(t) \leq -\sum_{i=1}^{N-1} \frac{\phi_i'(1)}{\phi_N'(1)} u_i(t), \quad (u_N(t) - g) \left( u_N(t) - \sum_{i=1}^{N-1} \frac{\phi_i'(1)}{\phi_N'(1)} u_i(t) \right) = 0. \quad (13)$$

Indeed, if (13) holds then (12) holds as well. By separating $u_N(t)$, we impose that $u_N(t)$ is no longer dictated by the ODE (9) but is dictated exclusively by condition (13). The following points introduce the remaining steps of the derivation (which will be elaborated in the upcoming sections):

1. To solve for the motion during inactive contact, the Signorini conditions (13) require that

$$u_N(t) = -\sum_{i=1}^{N-1} \frac{\phi_i'(1)}{\phi_N'(1)} u_i(t), \quad u_N(t) \leq g. \quad (14)$$

Substitution of the above into the FE approximation (8) effectively creates a family of shape functions that always satisfies the inactive contact conditions

$$u(x, t) \approx \sum_{i=1}^{N-1} \left( \phi_i(x) - \phi_N(x) \frac{\phi_i'(1)}{\phi_N'(1)} \right) u_i(t). \quad (15)$$

In other words, any solution obtained using the above approximation satisfies the homogeneous Neumann condition taking place during inactive contact. In the FEM, all functions $\phi_i(x)$ have local support and are non-zero for elements containing the node $u_i(t)$. In the case of NBM, the principle of local support is followed as well, and only the shape functions at the element including the contact node $u_N(t)$ are affected by the approximation (15). In turn, the inactive contact inequality in (14) reads

$$u_N(t) \leq g \quad \Rightarrow \quad -\sum_{i=1}^{N-1} \frac{\phi_i'(1)}{\phi_N'(1)} u_i(t) \leq g. \quad (16)$$

2. To solve for the motion during active contact, we construct a family of shape functions that always answers

$$u(L,t) \approx u_N(t) = g, \quad u_N(t) \leq -\sum_{i=1}^{N-1} \frac{\phi_i'(1)}{\phi_N'(1)} u_i(t). \quad (17)$$

Evidently, the set of shape functions answering this condition admits

$$u(x, t) \approx \sum_{i=1}^{N-1} \phi_i(x) u_i(t) - \phi_N(x) g \quad (18)$$

Thus, any solution that is obtained using the above approximation satisfies the Dirichlet condition in (17). In turn, the active contact inequality in (17), under this approximation becomes

$$u_N(t) \leq -\sum_{i=1}^{N-1} \frac{\phi_i'(1)}{\phi_N'(1)} u_i(t) \Rightarrow -\sum_{i=1}^{N-1} \frac{\phi_i'(1)}{\phi_N'(1)} u_i(t) \geq g. \quad (19)$$

3. The different sets of shape functions described in (15) and (18) will lead to different sets of ODEs governing inactive and active contact motions, respectively. The Signorini conditions are then satisfied by switching between the two sets of functions according to inequalities (16) and (19), which are mutually exclusive, as expected from the Signorini conditions. At the moment of switch, the internal nodal displacements and velocities (internal nodes are those with indexed $i = 1, 2, \dots, N-1$) are assumed to be continuous in time.

4. The NBM formulation results in an ODE, featuring discontinuous mass and stiffness matrices, which exhibits periodic solutions.
## 4.1 Comment on application of NBM using other shape functions

Although the NBM is derived in this manuscript using the classical FEM piecewise Lagrangian shape functions, this method can be also formulated using other shape functions. However, for other shape functions, precautions must be considered. One such precaution is that the shape functions must admit a stress approximation that is not always vanishing at the contacting end. To clarify, the NBM relies on the approximation of stress for the switching between active and inactive contact phases, as seen in Equation (12). In order to allow for a switching between contact phases, the shape function must be chosen such that the stress approximation at $x = 1$ does not exhibit $\phi_i'(1) = 0$ for all $i = 1, \dots, N$.

Here, the linear modes of the cantilever bar with uniform area

$$ \phi_i(x) = \sin\left(\frac{(2i - 1)\pi x}{2}\right), \quad i = 1, 2, \dots, N. \qquad (20) $$

will be investigated as an example of a set of shape functions that cannot be used in NBM. These mode shapes bear the undesired property $\phi_i'(1) = 0, \forall i$.

$$ g \ge u_N(t) \perp 0 \le g \ge u_N(t). \qquad (21) $$

In simple terms, Equation (21) shows that the linear modes of the cantilever bar exhibit inactive contact motions exclusively and therefore a switch between contact phases cannot occur.

The choice of shape functions is crucial for other methods in contact dynamics and not only the NBM. In fact, for the shape functions (20), both the mass redistribution method [15] and Nitsche method would fail. The mass redistribution method would fail in the same fashion as NBM since it relies on strong enforcement of the Signorini boundary condition. In turn, the Nitsche method would not be able to approximate adequately the Signorini problem. This is illustrated in the remainder of this section. In the Nitsche method, the Signorini boundary conditions are enforced via the following approximation of the stress at the contact boundary

$$ \bar{\sigma}(t) \approx -\max\left(0, \gamma(u_N(t) - g) - A(1)\sum_{i=1}^{N} \phi_i'(1)u_i(t)\right), \quad \gamma > 0 \qquad (22) $$

where $\gamma$ is set to be constant [16]. The Signorini condition is then satisfied as $N \to \infty$

$$ \bar{\sigma}(t) = \lim_{N \to \infty} -\max\left(0, \gamma(u_N(t) - g) - A(1)\sum_{i=1}^{N} \phi_i'(1)u_i(t)\right) \qquad (23) $$

$$ \bar{\sigma}(t) = -\max(0, \gamma(u(1, t) - g) - \bar{\sigma}(t)), \quad \gamma > 0 \qquad (24) $$

where the last term is equivalent to the Signorini condition (3). The convergence of the Nitsche method largely relies on the participation of the stress approximation in the right-hand side of Equation (23). This property is disrupted for the $\phi_i(x)$ in Equation (21). For these shape functions, the Nitsche stress approximation (22) reads

$$ \bar{\sigma}(t) \approx -\gamma \max(0, u_N(t) - g). \qquad (25) $$

Since the stress approximation is omitted from the right hand side of Equation (25), the resulting approximation (25) is equivalent to a penalty force with penalty parameter $\gamma$ [14]. The penalty method does not share the same convergence properties as Nitsche as can be seen that by taking the limit of Equation (25) as $N \to \infty$:

$$ \bar{\sigma}(t) = \lim_{N \to \infty} -\gamma \max(0, u_N(t) - g) = -\gamma \max(0, u(1, t) - g). \qquad (26) $$

This term, in contrast to Equation (24), is not equivalent to the Signorini conditions.

## 4.2 Inactive contact motion

The inactive contact motion condition (5) in the NBM framework is

$$ \sigma(t) \equiv A(1)u_x^h(1, t) = A(1)\sum_{i=1}^{N} \phi_i'(1)u_i(t) = 0, \quad u_N(t) \le g. \qquad (27) $$
where we use $\sigma(t)$ to denote the FE approximation of the stress at $x = 1$. To clarify, the true stress in the bar reads $\bar{\sigma}(t) \approx \sigma(t)$ according to Equation (11), and the homogeneous Neumann condition $\bar{\sigma}(t) = 0$ is satisfied by virtue of Equation (27). For inactive contact motion, we impose that $u_N(t)$ satisfies Equation (27) such that

$$u_N(t) = -\frac{1}{\phi'_{N}(1)} \sum_{i=1}^{N-1} \phi'_{i}(1) u_i(t) \equiv S(\mathbf{u}^{o}(t)), \quad \mathbf{u}^{o}(t) = (u_1(t) \ u_2(t) \ \dots \ u_{N-1}(t))^{\top} \quad (28)$$

Here, $\mathbf{u}^o(t)$ gathers internal nodal displacements whereas the contact node is represented in $u_N(t)$. This distinction is fixed and an internal node cannot turn into a contact node (or vice-versa) throughout the motion of the bar. This distinction between internal nodes and nodes prone to contact is also used in the MRM formulation [15]. Thus, in the NBM, the loci of contact nodes must be known a priori. Although the loci of contact nodes is not always known in advance, for the case of small strains, the loci of contact nodes can be assumed to be known and fixed without loss of generality.

Next, in Equation (28), we introduce the function $S$, acting on $\mathbf{u}^o(t)$, to simplify the notation. Actually, $S$ appears naturally throughout the derivation of both inactive and active contact motions in NBM. Indeed, this function is of integral importance to NBM and serves for multiple purposes: it defines the conditions for the switching between phases, the motion of the contact node during inactive contact phase, and the contact force applied on the bar during active contact phase. These roles of $S$ are illustrated in Figure 2 and are explicitly explained later in Section 4.4. For the remainder of this article, the function $S$ will be referred to as the *switching function*.

<figure><img src="image_3.png" /><figcaption>Figure 2: Switching function and associated quantities. The switching function dictates the phase of motion. For $S(\mathbf{u}^o(t)) \le g$, inactive contact motion takes place, and active contact motion takes place otherwise. Moreover, the displacement at the contact boundary abides $u_N(t) = S(\mathbf{u}^o(t))$ during the inactive contact phase, $S(\mathbf{u}^o(t)) \le g$, and the approximation of the stress at the contact boundary abides $\sigma(t) = A(1)\phi'_{N}(1)(g - S(\mathbf{u}^o(t)))$ during the active contact phase.</figcaption></figure>

Following the substitution $u_N(t)$ in (28), $\mathbf{u}(t)$ can be related to $\mathbf{u}^o(t)$ via a linear operator $\mathbf{B} \in \mathbb{R}^{N \times N-1}$

$$\mathbf{u}(t) = \mathbf{B}\mathbf{u}^o(t) \quad B_{ij} = \begin{cases} \delta_{ij} & i = 1, \dots, N-1 \text{ and } j = 1, \dots, N-1, \\ -\frac{\phi'_{j}(1)}{\phi'_{N}(1)} & i = N \text{ and } j = 1, \dots, N-1. \end{cases} \quad (29)$$

This representation will be helpful when we introduce the ODE governing the inactive contact motion. Specifically, it allows representation of the approximation of inactive contact conditions to be

$$u(x, t) \approx \mathbf{P}(x)\mathbf{B}\mathbf{u}^o(t), \quad u_N \le g \quad (30)$$

which is equivalent to the presentation in (15).

Next, the inequality condition for inactive contact motion (27) (non-penetration) can be put in terms of the nodes $\mathbf{u}^o(t)$ via the switching function

$$u_N \le g \Rightarrow S(\mathbf{u}^o) \le g \quad (31)$$
To derive the ODE governing the inactive contact motion, we insert the inactive motion constraints (29) and (31) into the weak-form of the PDE, Equation (9),

$$ \mathbf{w}^{\top}(t)(\mathbf{M}\mathbf{B}\ddot{\mathbf{u}}^{\circ}(t) + \mathbf{K}\mathbf{B}\mathbf{u}^{\circ}(t)) = 0 \quad (32) $$

$$ u_N(t) = S(\mathbf{u}^o(t)), \quad \ddot{u}_N(t) = S(\ddot{\mathbf{u}}^o(t)), \quad S(\mathbf{u}^o(t)) \le g. \quad (33) $$

In classical FEM, it is generally assumed Equation (32) is true for all $\mathbf{w}(t)$ and the subsequent omission of $\mathbf{w}(t)$ from the equation takes place. Here, however, omission of $\mathbf{w}(t)$ will lead to an over-defined system of ODEs ($N-1$ variables in $\mathbf{u}^o$ for $N$ equations). To remedy this, we use the Galerkin-Bubnov method where we project the residual (the term multiplying $\mathbf{w}$ in Equation (32)) on the same solution space used for $\mathbf{u}(t)$ [23]. This strategy has been proven successful in other applications and is commonly used when shape functions that satisfy the boundary conditions are involved in the approximation [23, 24]; [25, p. 300]. It is further noted that the strategy taken here in deriving the NBM has also been referred to as *basis recombination* [26, p. 112]. Following the Galerkin-Bubnov method, we project the residual resulting from the approximation (30) on the composing trial functions. Under representation (9), this results effectively in modification of the test-function as follows

$$ \mathbf{w}(t) = \mathbf{B}\mathbf{w}^o \quad (34) $$

where $\mathbf{w}^o$ gathers all test function contributions corresponding to the internal nodes $\mathbf{u}^o(t)$. Then, substitution of (34) into (32) and omission of $\mathbf{w}^o$ results in the reduced ODE (of $N-1$ equations)

$$ (\mathbf{M}_N \ddot{\mathbf{u}}^o(t) + \mathbf{K}_N \mathbf{u}^o(t)) = 0, \quad \mathbf{M}_N = \mathbf{B}^\top \mathbf{M} \mathbf{B}, \quad \mathbf{K}_N = \mathbf{B}^\top \mathbf{K} \mathbf{B}. \quad (35) $$

where the subscript $N$ in $\mathbf{M}_N$ (or $\mathbf{K}_N$) is used to denote the NBM coefficients corresponding to answering homogeneous *Neumann* conditions. We remark that, via the NBM, the displacement $u_N(t)$ has been effectively removed from the ODE such that neither the equation nor the inequality constraint in expression (35) includes $u_N(t)$.

## 4.3 Active contact motion

The active contact condition in NBM, derived from Equation (13), reads

$$ u_N(t) = g \quad \text{implying} \quad \dot{u}_N(t) = 0 \quad \text{and} \quad w_N(t) = 0 \quad (36) $$

Effectively, this approximation applies in the active contact phase and the velocity of the contact node is discontinuous at the moment of contact. Here, the test function vanishes on Dirichlet boundaries (ie, $w_N(t) = 0$ during active contact) which is common practice for FEM approximations [22, 23, 27]. Next, the complementarity condition (13) can be rewritten in terms of the switching function

$$ A(1) \sum_{i=1}^{N} \phi'_{i}(L) u_{i}(t) \le 0, \quad (37) $$

$$ -\sum_{i=1}^{N-1} \frac{\phi'_{i}(1)}{\phi'_{N}(1)} u_{i}(t) \ge u_{N}(t), \quad (38) $$

$$ S(\mathbf{u}^o(t)) \ge g. \quad (39) $$

Here, the transition between statements (37) and (38) requires that $\phi'_{N}(1) > 0$ holds. Otherwise, for example in the case $\phi'_{N}(1) \le 0$, both inactive and active contact phases would have to occur for $S(\mathbf{u}^o(t)) < 0$, that is simultaneously, which disagrees with the mutual exclusivity of the inactive and active contact conditions in the complementarity conditions (3). Accordingly, for $\phi'_{N}(1) > 0$, we note that the active contact motion occurs for $S(\mathbf{u}^o(t)) > g$, and the inactive contact motion occurs for $S(\mathbf{u}^o(t)) < g$ such that both are mutually exclusive. Fortunately, for the case of Lagrangian polynomials used here, the statement $\phi'_{N}(1) > 0$ has been proven to hold for any number of elements. The proof can be found in section A.1.

We continue with the substitution of expressions (36) and (39) into the FEM-ODE (9) where (36) is first recast in the matrix format

$$ \mathbf{w}(t) = \mathbf{B}^d \mathbf{w}^o, \quad \mathbf{u}(t) = \mathbf{B}^d \mathbf{u}^o(t) + g \mathbf{b}^d, \quad S(\mathbf{u}^o(t)) \ge g \quad (40) $$
with the notations

$$ B_{ij}^{d} = \begin{cases} \delta_{ij} & i = 1, \dots, N-1; j = 1, \dots, N-1 \\ 0 & i = N; j = 1, \dots, N-1 \end{cases} \quad \text{and} \quad b_{i}^{d} = \begin{cases} 0 & i = 1, \dots, N-1 \\ 1 & i = N. \end{cases} \qquad (41) $$

We then plug (40) into (9) to obtain an ODE in terms of $\mathbf{u}^o$ multiplied by $\mathbf{w}^o$ forming the scalar equation

$$ (\mathbf{w}^o)^{\top}(\mathbf{M}_D \ddot{\mathbf{u}}^o(t) + \mathbf{K}_D \mathbf{u}^o(t) + g\mathbf{f}_D) = 0, \quad S(\mathbf{u}^o(t)) \geq g \qquad (42) $$

where $\mathbf{M}_D = (\mathbf{B}^d)^{\top}\mathbf{M}\mathbf{B}^d$, $\mathbf{K}_D = (\mathbf{B}^d)^{\top}\mathbf{K}\mathbf{B}^d$, and $\mathbf{f}_D = (\mathbf{B}^d)^{\top}\mathbf{K}\mathbf{b}^d$, and the subscript $D$ is used to denote the coefficients corresponding the ODE answering the non-homogeneous Dirichlet conditions. Assuming Equation (42) should be valid for all values $\mathbf{w}^o$, the following ODE formulation is obtained:

$$ \mathbf{M}_D \ddot{\mathbf{u}}^o(t) + \mathbf{K}_D \mathbf{u}^o(t) + g\mathbf{f}_D = \mathbf{0}, \quad S(\mathbf{u}^o(t)) \geq g \qquad (43) $$

Equation (43) is equivalent to the ODE obtained by application of classical Finite Element method on the clamped-clamped bar. Furthermore, the NBM formulation of the ODE for active contact is equivalent to this described by the basis recombination method for non-homogeneous boundary conditions [26, p. 112]. However, the NBM adds the restriction $S(\mathbf{u}^o(t)) > g$ to infer that the bar must be repulsed at all time throughout contact, as required by the active contact condition (4).

In the next section, we combine the ODEs corresponding to both active and inactive motions, Equations (35) and (43) respectively, to construct the ODE approximation for the original Signorini problem.

## 4.4 NBM-FEM formulation of Signorini problem

The switching method [12, 28] is used for the enforcement of the Signorini conditions in NBM. In the switching method, the Signorini complementarity conditions are answered by alternating between the boundary conditions (4) and (5) both in the test and trial functions such that the inequality constraints are satisfied. In NBM, this translates to switching between Equations (35) and (43), and the complete NBM-ODE reads

$$ \begin{cases} \mathbf{M}_D \ddot{\mathbf{u}}^o(t) + \mathbf{K}_D \mathbf{u}^o(t) + g\mathbf{f}_D = \mathbf{0} & S(\mathbf{u}^o(t)) \geq g \\ \mathbf{M}_N \ddot{\mathbf{u}}^o(t) + \mathbf{K}_N \mathbf{u}^o(t) = \mathbf{0} & S(\mathbf{u}^o(t)) \leq g. \end{cases} \qquad (44) $$

We note that at the moment of switch, denoted $t_s$ such that $S(\mathbf{u}^o(t_s)) = g$, the resulting NBM-ODE (44) raises two conflicting definitions to the ODE. This conflict is resolved by extending the active and inactive contact NBM conditions as follows

$$ \text{Active contact NBM: } S(\mathbf{u}^o(t)) > g, \text{ or } S(\mathbf{u}^o(t)) = g \text{ and } S(\dot{\mathbf{u}}^o(t^-)) > 0 \qquad (45) $$

$$ \text{Inactive contact NBM: } S(\mathbf{u}^o(t)) < g, \text{ or } S(\mathbf{u}^o(t)) = g \text{ and } S(\dot{\mathbf{u}}^o(t^-)) < 0 \qquad (46) $$

which is mathematically sound since the velocity at the moment of switch indicates whether an active contact or an inactive contact occurs after the switch. Furthermore, note that the case of zero velocity $S(\dot{\mathbf{u}}^o(t_s^-)) = 0$ before contact is not investigated in definitions (45) and (46). Such solutions are referred to as grazing solutions and propose a challenge that is beyond the scope of this paper (the reader may refer to [10, p. 385] for some of the intricacies involved with determining grazing periodic motions). Since grazing motions are excluded, the NBM is limited in its solutions. Thus, it is important to note that the NBM is not presented as a scheme for generation physically accurate simulations. Rather, the NBM is used for modal analysis and detection of periodic Signorini compliant non-grazing motions.

Next, at the instant of contact, we impose that the internal displacements $\mathbf{u}^o(t)$ and velocities $\dot{\mathbf{u}}^o(t_s)$ are continuous at switching instants

$$ \text{Continuity of internal states: } \mathbf{u}^o(t^+) = \mathbf{u}^o(t^-), \quad \dot{\mathbf{u}}^o(t^+) = \dot{\mathbf{u}}^o(t^-), \quad S(\mathbf{u}^o(t)) = g. \qquad (47) $$

while only the contact node is characterized by discontinuous velocities, as will be shown later in Equation (51). The condition for continuity of the internal states (47) corresponds to common application of the Signorini conditions both numerically and analytically [1, 12, 14–16]. Thus, both internal displacements and internal velocities are assumed to be always continuous and only the acceleration $\ddot{\mathbf{u}}^o(t)$ is discontinuous at instants $t_s$ (as
can be deduced from Equation (44)). The discontinuity of internal accelerations at the moment of switch is a consequence of the NBM and is not expected in the true solution. However, solutions of the NBM still show good agreement with the true solution as $N \to \infty$, as illustrated in section 6.

Given that the acceleration is discontinuous at the moment of switch, it is more convenient to represent the ODE (44) in terms of the acceleration at $t^+$

$$ \ddot{u}^o(t^+) = \begin{cases} -(M_D)^{-1}(K_D u^o(t) + g f_D) & \text{active contact NBM} \\ -(M_N)^{-1}(K_N u^o(t)) & \text{inactive contact NBM} \end{cases} \qquad (48) $$

where $M_D$ and $M_N$ are always invertible (the proof follows from $M$ being non-invertible and is not presented in this manuscript for sake of conciseness).

In turn, the approximation of $u(x, t)$ in NBM is defined via expressions (30) and (40)

$$ u(x, t) \approx u^h(x, t) = \mathbf{P}(x)\mathbf{u}(t) = \begin{cases} \mathbf{P}(x)(\mathbf{B}^d\mathbf{u}^o(t) + g\mathbf{b}^d) & S(\mathbf{u}^o(t)) \ge g \\ \mathbf{P}(x)\mathbf{B}\mathbf{u}^o(t) & S(\mathbf{u}^o(t)) \le g \end{cases} \qquad (49) $$

At last, from expression (49), we obtain nonsmooth expressions for the displacement, velocity and stress at the contacting end (with strict inequality applied on the active contact condition)

$$ u(1, t) \approx u_N(\mathbf{u}^o(t)) = \begin{cases} g & S(\mathbf{u}^o(t)) \ge g \\ S(\mathbf{u}^o(t)) & S(\mathbf{u}^o(t)) \le g \end{cases} \qquad (50) $$

$$ u_t(1, t) \approx \dot{u}_N(\mathbf{u}^o(t), \dot{\mathbf{u}}^o(t)) = \begin{cases} 0 & S(\mathbf{u}^o(t)) > g \\ S(\dot{\mathbf{u}}^o(t)) & S(\mathbf{u}^o(t)) < g \end{cases} \qquad (51) $$

$$ \bar{\sigma}(t) \approx \sigma(\mathbf{u}^o(t)) = \begin{cases} A(1)\phi'_N(1)(g - S(\mathbf{u}^o(t))) & S(\mathbf{u}^o(t)) \ge g \\ 0 & S(\mathbf{u}^o(t)) \le g \end{cases} \qquad (52) $$

Here, we note that $\sigma(t)$ is continuous in the NBM formulation contrarily to the discontinuous behaviour of the true solution to the Signorini problem [1, 12, 20]. However, in Section 6 it will be shown that the method still converges. Furthermore, it is important to note that schemes utilizing Newton's impact law with $e = 0$ are characterized by continuous contact pressure as well and still show convergence to the true solution [14].

The Signorini problem is hence formulated as the nonsmooth ODE (48) in $\mathbf{u}^o(t)$ with unique solutions given initial conditions $\mathbf{u}^o(0)$ and $\dot{\mathbf{u}}^o(0)$. Furthermore, the solutions generated by the NBM are characterized by *sticking phases*. Sticking phases are continuous intervals of time of non-zero measure where active contact motion occurs. To clarify, sticking phases stand in contrast to chattering exhibited by schemes utilizing a Newton impact law with $e = 1$ [1]. This is a noteworthy property of this scheme since sticking phases occur in the true solution to the Signorini problem [1, 11, 12, 14].

On another note, the acceleration $\ddot{u}_N(t)$ obtained by differentiating Equation (51) involves the Dirac-delta distribution at the instance of switch, which may affect the formulation of the NBM-ODE since it participates in its definition, see Equation (33). However, the influence of the Dirac-delta in the NBM formulation was not investigated in this article and is suppressed to simplify the formulation. Nevertheless, numerical experiments show that the NBM-ODE (48) admits solutions that converge to the true motion for large $N$. Such numerical experiment is explored in Section 6.1.

Next, since periodic solutions require energy conservation, we investigate the energy conservation properties of solutions to the NBM-ODE.

## 4.5 Energy conservation properties of solutions to NBM-ODE

Solutions to the NBM-ODE (48) are equipped with the energy metric

$$ 2E(t) = \int_{0}^{1} u_{t}(x, t)^{2} + A(x)u_{x}(x, t)^{2} dx \approx \dot{\mathbf{u}}^{\top}(t)\mathbf{M}\dot{\mathbf{u}}(t) + \mathbf{u}^{\top}(t)\mathbf{K}\mathbf{u}(t) \qquad (53) $$

and exhibit the following properties

1. The ODE preserves energy for $S(\mathbf{u}^o(t)) > g$ and $S(\mathbf{u}^o(t)) < g$, away from instants $t_s$ such that $S(\mathbf{u}^o(t_s)) = g$.
2. At a time instant $t_s$, where a transition between active and inactive contact occurs (namely, at an instant where $S(\mathbf{u}^o(t_s)) = g$ and $S(\dot{\mathbf{u}}^o(t_i)) \neq 0$), an instantaneous change in energy $\Delta E$ occurs

$$ \Delta E = E(t_s^+) - E(t_s^-) = -|S(\dot{\mathbf{u}}^o(t_s))| \left( \frac{1}{2} M_{NN} S(\dot{\mathbf{u}}^o(t_s)) + \sum_{j=1}^{N-1} M_{Nj} \dot{u}_j(t_s) \right). \quad (54) $$

The energy after transition may either decrease ($\Delta E < 0$), increase ($\Delta E > 0$) or be conserved ($\Delta E = 0$). For proofs, Appendix A.2 establishes lemmas regarding Statement 1 and Appendix A.3 details the proof for Statement 2.

An immediate consequence of Statements 1 and 2 is that solutions of the NBM-ODE may exhibit a periodic energy evolution in time. Indeed, this property is favourable for detection of periodic solutions, and it is expected that the NBM can be used for detection periodic solutions to the Signorini problem. From numerical experiments, such periodic solutions to the NBM were found. An example of a periodic solution is shown in Figure 3 from which it is clear that the NBM allows for existence of periodic solutions with sticking phases. In contrast, solutions

<figure><img src="image_4.png" /><figcaption>Figure 3: NBM periodic solution for the varying area bar $A(x) = 1 - x/2$ and $g = 0.1$ of 10 elements and quadratic Lagrangian shape functions ($N = 20$). This model is used to illustrate different aspects of the NBM methodology (a similar model with different $N$ is used in Figure 4). Note that the total energy is dissipated at the moment of contact but is completely regained at the end of the contact phase. Due to this characteristic, the NBM allows for periodic solutions with sticking phases. In Figure 3(a), and all the following displacement field plots, shaded rectangular surfaces and dotted lines highlight active contact phases and gap from wall, respectively.</figcaption></figure>

obtained via Nitsche method exhibit sticking phases only at convergence (ie, for high number of elements) [16] and solutions of scheme with Newton impact law $e = 1$ exhibit chattering [14]. While solutions to Newton's impact law with $e = 0$ exhibit sticking phases, they also exhibit energy dissipation for non-zero impact velocity [14] and therefore do not allow for periodic solutions. The NBM, in comparison to Newton's impact law $e = 0$, allows for regain of energy in transition from active to inactive contact.

Indeed, the existence of periodic solutions in NBM relates largely to the fact that energy can be regained throughout the motion, as seen in Figure 3.

We note here that the NBM does not have an explicit impact law and the energy jumps (whether loss or gain) are instead explained by the change in shape functions at the moment of switch.

In the true solution, at the instance of switch, $u_t(x, t_s^+) = u_t(x, t_s^-)$ holds true for $x \in (0, 1) [1, 11, 20]$. However, in NBM, the approximations of $u_t(x, t_s^+)$ and $u_t(x, t_s^-)$ involve two different sets of shape functions, described in Equations (16) and (19). For example, during closing contact (switch from inactive to active contact), the contact node's velocity $\dot{u}_N$ vanishes as dictated by the shape functions in Equation (19) and a loss in energy ensues at this instance. This is similar to application of $e = 0$ using Newton impact law in FEM. However, in opening contact (switch from active to inactive contact), there is a gain of energy since the shape functions dictating the inactive contact motion (16) allow for non-zero $\dot{u}_N$ after the switch. This is different than application of $e = 0$ using Newton impact law in FEM where the velocity $\dot{u}_N$ is zero at opening contact. In sum, the regain in energy is a consequence of the switching mechanism and shape functions introduced in NBM. To better illustrate this mechanism, we use the example of the NBM model of the bar with linear shape functions.
For the NBM model of the bar with linear shape functions we have $\phi'_{N}(1) = 1$, $\phi'_{N-1}(1) = -1$ and $\phi'_{i}(1) = 0$ for $i = 1, 2, \dots, N-2$. In turn, for the bar with linear shape functions, the displacement and velocity of the contact node during inactive contact, $S(\mathbf{u}^{o}(t)) < g$, are given as follows:

$$u(x = 1, t) \approx u_N(t) = S(\mathbf{u}^o(t)) = u_{N-1}(t) \qquad (55)$$

$$u_t(x = 1, t) \approx \dot{u}_N(t) = S(\dot{\mathbf{u}}^o(t)) = \dot{u}_{N-1}(t). \qquad (56)$$

where Equations (55) and (56) are deduced from Equations (50) and (51) for linear Lagrange shape functions, respectively. From Equations (55) and (56), we see that for linear shape functions, during inactive contact motion, $u_N(t) = u_{N-1}(t)$ holds. Physically, it means that the last element of the bar acts as a rigid body throughout inactive contact motion. During active contact motion, $S(\mathbf{u}^o(t)) > g$, the last element of the bar acts as an elastic body similar to FEM with $e = 0$ since the approximating shape functions used (19) are identical to implementation of a non-homogeneous Dirichlet condition in FEM. Thus, upon closing contact, the last element of the bar transforms from a rigid element to an elastic element and loses energy due to vanishing $\dot{u}_N$. However, in opening contact, the bar transforms from an elastic element to a rigid element, and the bar regains energy as the velocity of the contact node equals the velocity of the previous node $\dot{u}_N(t) = \dot{u}_{N-1}(t)$. For higher order shape functions, a similar mechanism takes effect, however it cannot be illustrated as transformation of the last element between rigid and elastic body. Rather, the jump in energy is explained as a consequence of changing between two families of shape functions of the NBM.

Along the same line, the terms dissipation and gain of energy will be used to describe the evolution of the energy metric in time rather than implying any physical energy transfer in or out of the system. Moreover, in the NBM, it was evident from numerical experiments that the energy jump $\Delta E$ diminishes for large $N$. An example of such experiment is presented in Section 6.1. Nevertheless, an analytical proof for the convergence of $\Delta E$ to zero with higher number of elements is not presented in this article.

To conclude, the behaviour of the energy metric in the NBM should be understood as a consequence of the Galerkin-Bubnov method rather than as a physical imposition on the system (as done by implementing a Newton impact law with $e = 0$, for example).

## 4.6 Notes on the NBM-FEM formulation

In this section, we will explore differences and similarities between the NBM and other existing methods for treatment of the Signorini conditions.

The Signorini conditions is implemented in the NBM-ODE in Equation (48) via switching between systems of ODEs. We note that the Wave-FEM [12] also implements the Signorini conditions via switching between two discrete dynamical systems (one governing inactive contact motion and the other governing active contact motion). Also, the NBM can be considered similar to the MRM [15] since the MRM effectively modifies the stiffness matrix on the switch between active and inactive contact whereas the NBM modifies both the mass and stiffness matrix on the switch.

On another note, we would like to note that the NBM-ODE does not introduce explicitly an impact law. To clarify, in schemes utilizing Newton's impact law, different coefficients of restitution may be used $e \in [0, 1]$ to solve the Signorini problem for a single FE model. In contrast, the NBM formulation of the Signorini conditions is unique with respect to the FE model. Still, we can draw a parallel between the NBM and FEM schemes utilizing Newton's impact law. The NBM coincides with Newton impact law $e = 0$ during the active contact motion and away from switching instants. There, the NBM abides exactly the same ODE abided by the Newton impact law $e = 0$ since the shape functions discretizing the governing PDE are equivalent to those used in classical FEM for a fixed $u_N(t) = g$. However, during inactive contact motions, the ODE dictating the motion of nodes is distinct from the one used in implementation of Newton's impact law since the set of NBM shape functions to discretize the PDE (16) is distinct from this used in classical FEM. For example, for the case of first-order shape functions, the NBM preserves the length of the contact element during inactive contact motion: $u_N(t) = S(\mathbf{u}^o(t)) = u_{N-1}(t)$ (ie, the contact node has the same displacement as the node before it) which does not occur in the classical FEM scheme.
# 5 Nonsmooth modal analysis

In the preceding sections, the NBM and FEM were used to approximate the solution to the initial boundary value problem exposed in Equation (1) to Equation (4). To solve for the remaining conditions, we require that $\mathbf{u}^o(t)$ and $\dot{\mathbf{u}}^o(t)$ are periodic, see Equations (6) and (7),

$$ \mathbf{u}^o(0) = \mathbf{u}^o(T) \qquad (57) $$

$$ \dot{\mathbf{u}}^o(0) = \dot{\mathbf{u}}^o(T). \qquad (58) $$

In this article, we attempt to find such solutions and corresponding period using the shooting method. Moreover, continuation is used for detection of nonsmooth modes, ie families of periodic solutions [7].

## 5.1 Time-marching techniques and shooting method

In the shooting method [7], we find periodic solutions by aiming for a set of initial conditions, $\mathbf{u}_0^o$ and $\mathbf{v}_0^o$, generating a periodic solutions, ie

$$ \begin{aligned} \mathbf{u}^o(T) &= \mathbf{u}^o(0) \equiv \mathbf{u}_0^o \\ \dot{\mathbf{u}}^o(T) &= \dot{\mathbf{u}}^o(0) \equiv \mathbf{v}_0^o \end{aligned} \qquad (59) $$

where $\mathbf{u}^o(t)$ is subject to the ODE (48).

In order to solve for the shooting equations, we solve for $\mathbf{u}^o(T)$ and $\dot{\mathbf{u}}^o(T)$ numerically via a numerical time marching scheme. For example, a Crank-Nicolson (CN) scheme (or Newmark scheme with $\beta = 1/2$ and $\gamma = 1/4$ [29]) may be used. While the CN scheme is considered energy stable for finite element schemes [29], it is an implicit algorithm and requires the implementation of a root solving algorithm. To implement the time-marching techniques, we shall first convert the NBM-ODE (48) into a system of first order ODEs $\dot{\mathbf{q}}(t) = \mathbf{G}(\mathbf{q}(t))$) with

$$ \mathbf{G}(\mathbf{q}(t)) = \begin{cases} \begin{pmatrix} \mathbf{v}^o(t) \\ \mathbf{M}_D^{-1}(g\mathbf{f}_D - \mathbf{K}_D\mathbf{u}^o(t)) \\ -\mathbf{M}_N^{-1}\mathbf{K}_N\mathbf{u}^o(t) \end{pmatrix} & \text{Active contact NBM} \\ \begin{pmatrix} \mathbf{v}^o(t) \\ -\mathbf{M}_N^{-1}\mathbf{K}_N\mathbf{u}^o(t) \end{pmatrix} & \text{Inactive contact NBM} \end{cases} \quad \text{and} \quad \mathbf{q}(t) = \begin{pmatrix} \mathbf{u}^o(t) \\ \mathbf{v}^o(t) \end{pmatrix}. \qquad (60) $$

For application of the time-marching scheme, a discretization of the time span $t \in [0, T]$ into $N_t$ steps is considered with the time-step $\Delta t = T/N_t$. In the time-marching scheme, the quantities $\mathbf{q}_i \approx \mathbf{q}(t_i)$ are solved iteratively for $i = 1, \dots, N_t$. After $N_t$ steps, the state at the end of the period is obtained $\mathbf{q}(T) \approx \mathbf{q}_{N_t}$. Indeed, the approximated state $\mathbf{q}_{N_t}$ is a function of $T$ and $\mathbf{q}_0$. Substitution of the approximation $\mathbf{q}_{N_t}(\mathbf{q}_0, T)$ in Equation (59) reads

$$ \mathbf{F}(\mathbf{q}_0, T) = \mathbf{q}_{N_t}(T, \mathbf{q}_0) - \mathbf{q}_0 = \mathbf{0}. \qquad (61) $$

This problem constitutes an under-defined system of equations with $2(N-1)$ equations for $2(N-1)+1$ variables. This gives rise to a continuum of solutions parameterized in $T$, ie $(\mathbf{q}_0(T), T)$. This continuum of solutions then directly corresponds to families of periodic solutions, otherwise known as nonsmooth modes (NSMs) [1].

## 5.2 Error estimation

For cases where the solution $u(x, t)$ exists, a proper error estimate of the NBM solution is the $L_2$ norm $||u^h(x, t) - u(x, t)||_2$ where $u^h(x, t)$ is defined in Equation (49). However, in the absence of closed-form solution, such as the periodic solutions sought for the varying-area bar, other error metrics are needed. While the error in the residual resulting from the approximation constitutes a common error metric for the Galerkin-Bubnov method [23, 30], the residual of the PDE (1) under the NBM-FEM approximation requires knowledge of $\phi_i''(x)$ for $x \in [0, 1]$. However, since $\phi_i(x)$ is described using the piecewise Lagrangian used in FE framework, the double derivative of $\phi_i(x)$ is not defined on element boundaries. Therefore, to quantify the error, the residual estimator [30, p. 93]

$$ R(t) = h^4 \sum_{j=1}^{N_e} \int_{\mathcal{E}_j} \left( \sum_{i=1}^{N} \phi_i(x) \ddot{u}_i(t) - \phi_i''(x) u_i(t) \right)^2 dx \qquad (62) $$
is used, where h describes the length of the element, $N_e$ denotes the total number of elements and $\mathcal{E}_j$ is the domain of the element j excluding its boundaries such that any given $\phi_i''(x)$ is defined everywhere in $\mathcal{E}_j$, and the boundary of the elements are excluded from the error metric (62). It is noted, that while the error metric effectively excludes points of discontinuity, it evaluates the accuracy of the solution for ranges where the approximation of $u(x, t)$ is clearly defined. Thus, it is considered a proper metric for the evaluation of the solution's accuracy. For the NBM, we must evaluate the integral of $R(t)$ for $t \in [0, T]$. Since the acceleration $\ddot{u}(t)$ is discontinuous at the moment of switch $t_s$, we define the residual error by excluding instances of discontinuity, similarly to the residual estimator,

$$R_{\epsilon} = \int_{0}^{t_s^{1-}} R(t) dt + \int_{t_s^{(N_s-1)+}}^{T} R(t) dt + \sum_{j=1}^{N_s} \int_{t_s^{j+}}^{t_s^{(j+1)-}} R(t) dt \quad (63)$$

where $t_s^j$ denote various distinct instants of switch through the motion and $N_s$ defines the total instances of switch in $[0, T]$. We duly note that the residual error does not take into account the discontinuities in time, and may be an inaccurate error metric for the Signorini problem. However, the metric $R_{\epsilon}$ does define the quality of the approximation of inactive and active phases of motion by their respective shape functions. Thus, $R_{\epsilon}$ is useful in determining the accuracy in the approximation of the active and inactive contact phases. Nevertheless, if the exact solution to the problem is known, the more accurate error norm $||u^h(x, t) - u(x, t)||_2$ will be used instead of $R_{\epsilon}$.

## 5.3 Sequential continuation with correction

As evident from literature on modal analysis of the Signorini problem, NSM can be described by via a continuum of solutions $(\mathbf{q}_0(T), T)$ on a closed interval of periods, referred to as a backbone curve [1, 11, 12, 31]. This continuum of solutions can be found by applying sequential continuation on the system of equations (61) [12, 31]. Applying sequential continuation to the NBM numerically, it has been noticed that the sequential continuation for curves of $N \ge 10$ faces difficulties in obtaining solutions on the backbone curve $(\mathbf{q}_0(T), T)$. To distinguish between curves approximated using different $N$, we will denote a solution continuum as follows: $(\mathbf{q}_0(T), T)^N$.

For large $N$, it is more efficient to obtain solutions by first finding solution with a low $N$ and then apply a shooting algorithm while recursively increasing the number of nodes for each point in the low-$N$ curve. We refer to this method as sequential continuation with correction (SCC). The steps of sequential continuation with correction are:

1. Obtain the nonlinear normal mode with low $N$ (for instance, $N = 4$), stored as the series $(\mathbf{q}_0^k, T^k)^4$, $k = 1, 2, \dots, N_c$.

2. For the point $i = 1$ in the series $(\mathbf{q}_0^1, T^1)^N$, perform shooting on a system with higher number of nodes $N^+$. First, interpolate $(\mathbf{q}_0^1)^N$ using the shape functions $P(x)$ to obtain an initial guess for the desired $N^+$ approximation.

3. Solve the shooting equations with $N^+$ nodes with period $T^1$ and obtain $(\mathbf{q}_0^1, T^1)^{N^+}$.

4. Steps 2 and 3 above may be repeated for higher number of nodes while keeping the period $T^1$ constant.

5. Repeat 2, 3 and 4 for all points in the series discovered in step 1, ie $k = 2, 3, \dots, N_c$. Figures 4 and 5 illustrate the method.

# 6 Results

## 6.1 Convergence of Crank-Nicolson and NBM

In this section, we verify the validity of the NBM for the cantilever bar of uniform area, ie $A(x) = 1, \forall x \in [0, 1]$ and $g = 0.1$. For this model, analytical solutions and nonsmooth modal analysis results are available [1, 11, 13]. We compare the exact solution with its NBM approximation for $T = 3.5$. The corresponding NSM displacement field reads [11]

$$u(x,t) = f(t+x) - f(t-x) \quad \text{with} \quad f(s) = 0.1 \begin{cases} -s & s \in [-1, 1], \\ s-2 & s \in [1, 2.5], \\ 3-s & s \in [2.5, 4.5]. \end{cases} \quad (64)$$
<figure><img src="image_5.png" /><figcaption>Figure 4: Sequential Continuation with Correction: every point represents a periodic solution of specific frequency and energy. The curve for $N = 4$ is found using sequential continuation. For each point on this curve, a shooting algorithm is applied to obtain a periodic solution in higher $N$ and the same frequency. Dotted line relates to Figure 5 depicting solutions along this line. The results in this figure and in Figure 5 were obtained using CN and NBM with $g = 0.1$ and $A(x) = 1 - x/2$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Periodic displacement field corresponding to Figure 4, for period $T = 3.7812$. With higher number of nodes, the obtained displacement field is more accurate as evident by the lower $R_\epsilon$ value.</figcaption></figure>

The NBM model is assigned the initial conditions generating the exact periodic solutions. The exact initial conditions are discretized and their values are taken at loci $x_i$ corresponding to the NBM-FEM nodes $u_i(t)$

$$u_i^o(0) = f(x_i) - f(-x_i) = -0.2x_i \quad \text{and} \quad v_i^o(0) = f'(x_i) - f'(-x_i) = 0 \qquad (65)$$

inserted in the implicit CN time-marching with $N_t = 2000$ steps and $\Delta t = 1.75 \cdot 10^{-3}$. The error used in the convergence analysis is expressed in the $L_2$-norm

$$\| \mathbf{P}(x)\mathbf{u}(t) - u(x, t) \|_2 = \frac{1}{T} \sqrt{ \int_0^1 \int_0^T (u(x, t) - \sum_{i=1}^{N-1} \phi_i(x)u_i^o(t) + \phi_N(x)u_N(\mathbf{u}^o(t)))^2 dt } \quad (66)$$

where $u_N(\mathbf{u}^o(t))$ is defined in Equation (50). A sample of the solution for $N = 200$ is illustrated in Figure 6. The motion plotted in Figure 6 shows very close to piecewise linear surface which is indeed the expected solution for the given initial conditions as presented in Equation (64).

Figure 7 denotes convergence of $L_2$ error and $\Delta E$ for the cantilever bar solution with NBM and C-N algorithm. From Figure 7(a), we note that, for any shape function polynomial degree, convergence is approximately of first order in terms of the number of elements. Indeed, since the exact solution is only once piecewise differentiable as evident from Equation (64), the order of convergence of the FE methods is limited to first order [22, p. 117]. Since the order of convergence in FE based method is expected to be linear, it can be said that the number of elements affects the solution's accuracy more than the order of shape functions. However, it was found that for application of shooting method, dealing with higher shape functions has allowed to reduce the number of variables
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: NBM solution emerging from initial conditions (65) for 100 elements and quadratic shape functions ($N = 200$). We note a small disturbance in the velocity field Figure 6(b). This disturbance in the velocity field seems to travel along characteristic lines and diminishes progressively as higher number of elements and degree of polynomials are used. In Figure 6(b), and all following velocity field plots, a thick line on the $t$ axis highlights active contact phases.</figcaption></figure>

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

<figure><img src="image_10.png" /><figcaption>Figure 7: Error plots for NBM with Crank-Nicolson algorithm for a nonsmooth motion of the cantilever bar. Curves denote different order of Lagrangian shape function in the FE approximation.</figcaption></figure>

while reducing both the residual error and computation time. Thus, in the upcoming nonsmooth modal analysis Section 6.2, the authors deemed a set of 20 elements with quadratic shape functions ideal in terms of accuracy and computation time.

Furthermore, in Figure 7(b), we see that as the number of elements or order of shape functions increases, the energy jump $\Delta E$ (discussed previously in Section 4.5) decreases. Compared to the $L_2$ convergence rate in Figure 7(a), the rate of convergence $\Delta E \to 0$ seems to be affected by the order of shape functions and shows a higher rate with increasing polynomial order. This indicates that, in NBM, a higher order of shape functions or higher number of elements allows for better energy conservation throughout the motion. Although, besides the empirical evidence shown here, no analytical investigation was done to prove this statement.

## 6.2 Nonsmooth modes

In this section, we shall use the nonsmooth modal analysis techniques developed in Section 5 to perform nonsmooth modal analysis for three variations of the bar: the internally resonant cantilever uniform area bar, the uniform area cantilever bar with soft support, and the varying area cantilever bar. To verify the validity of the nonsmooth modes presented here, we will compare each nonsmooth mode with its corresponding forced-response
diagram since it is expected that the backbone branch will align with the frequency and energy at resonance [7, 13]. All results in Section 6.2 are generated for a gap distance g = 0.001 to comply with the models investigated in [13, 31]. Furthermore, all NSMs were generated via FE models of 20 elements and quadratic shape functions for which the backbone branches were depicting resonant points sufficiently.

### 6.2.1 Forced-response curves

The forced response-curves are generated by solving the equation

$$u_{tt}(x,t) + cu_t(x,t) + (A(x)u_x(x,t))_x = \bar{F} \cos(\omega t) \qquad (67)$$

where c, $\bar{F}$ and $\omega$ denote the damping coefficient, forcing amplitude and the forcing frequency, respectively. This governing PDE is complemented with the Signorini boundary conditions (5) and (4) as well as the boundary conditions imposed at $x = 0$ by the model in questions. The resulting Signorini problem is then solved for $c = 0.1, 0.2, \dots, 0.7$ and for frequencies $\omega$ within the range of the detected NSM. For each set of values c and $\omega$, we record the sum of kinetic and potential energies of the structure at steady state to plot the forced-response diagrams.

In practice, it is assumed a steady state is reached as $t \to \infty$. Although we often expect a forced motion to reach a periodic steady state, for some frequencies, quasi-periodic or chaotic solutions take place [12]. Thus, to obtain the forced-response curve, the Signorini problem is solved until a periodic motion is obtained or until the energy's mean value throughout a forcing period is sufficiently stable.

All forced-response curves were obtained using the FEM framework with 20 elements and quadratic shape functions and solved via SICONOS [32]. SICONOS uses a Moreau-Jean scheme to implement the Newton-impact law in the resulting system of ODEs. Here, the Newton-impact law is applied on a classical FEM approximation of the model in question, that is without application of NBM. The coefficient of restitution used to generate the forced responses in Sections 6.2.2 to 6.2.4 is $e = 0$. It is important to note that forced-response diagrams for the models have been also obtained using SICONOS with $e = 1$, Nitsche's method, and NBM. Since all the forced-response curves yielded very similar results, we chose to present just those obtained for SICONOS with $e = 0$ for sake of conciseness.

### 6.2.2 Constant cross-section and internally resonant bar

The internally resonant bar of uniform area $A(x) = 1$ is an example that has been investigated both numerically [12, 31] and analytically [1, 11, 20]. The system is known to manifest an intricate modal space consisting of families of iso-periodic periodic solutions in a dense set of periods and families of periodic solutions of the same frequency and energy.

In Figure 8(a), we compare the curve obtained by NBM to the curve of piecewise-linear analytical solutions obtained in [13]. It is clear that the NSM obtained from NBM lays closely to the exact piecewise-linear NSM obtained in literature [1, 11, 13]. We note that there exist other exact periodic solutions on top of the NSM curve as was concluded in [11]. It is known, from exact solutions and numerical solutions, that the spectrum of the internally resonant bar consists of iso-periodic nonsmooth modes existing as lines above the exact NSM branch in the frequency-energy diagram [11]. Therefore, the group of solutions found by NBM may be considered numerically accurate due to the existence of solutions above the Exact NSM curve. Nevertheless, we note that the NBM backbone curve does cross all points of resonance in the forced-response diagram. This shows that the NBM is useful for the modal analysis of the Signorini problem.

### 6.2.3 Constant cross-section bar with soft support

Similarly to [13, 31], the bar with soft support features a uniform cross-sectional area $A(x) = 1$ where the homogeneous Dirichlet condition at $x = 0$ is replaced with the Robin condition

$$ku(0,t) = u_x(0,t), \quad k = \bar{k}L/(EA_0) \qquad (68)$$

where $\bar{k}$ is the physical stiffness coefficient of the spring. The NBM was used to handle the Signorini boundary condition at $x = 1$ while the soft support condition (68) was treated using the classical finite element technique.
<figure><img src="image_11.png" /><figcaption>(a) Forced response curves</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Velocity field of solution *a*, autonomous motion generated using NBM with CN time-marching</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Displacement field of solution *a*, forced and damped motion generated by SICONOS with $e = 0$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) Displacement field of solution *a*, autonomous motion using NBM with CN time-marching</figcaption></figure>

<caption>Figure 8: NSM of the internally resonant bar detected by the NBM. ($\bar{F} = 0.05$). Exact NSM is the piecewise-linear mode detected for the internally resonant bar in [1, 13]. $\omega_1 = \pi/2$ and $E_1 \approx 5.9 \times 10^{-7}$. We note here that the velocity field Figure 8(c) involves porous oscillations which are common to methods in the FE framework. However, a heat-map representation of the velocity field is useful in demonstrating that the found solution follows (approximately) the characteristic lines, exhibited by the exact solution [13].</caption>

For this experiment, we set $k = 0.5$ to replicate the results in [13, 31]. Corresponding results are shown in Figure 9. Again, as in the case of the internally resonant bar, the alignment between occurrences of resonance and the NSM detected by NBM is clear. The motions obtained in the NBM analysis are similar to those obtained in [31] for the same values. Specifically, the displacement fields depicted in Figures 9(b) and 9(c) are similar to those presented in [31, (a) and (b) in Figure 11].

In Figure 9(a), it seems apparent that both displacements relate to two different branches of the solution. Moreover, the forced response curves of less damped motions have two peaks which may indicate the existence of two distinct NSM branches. Here, the branch corresponding to solution $q_2$ has been detected until a maximal energy point. Backward sequential continuation has then revealed a distinct set of points to which the solution $q_1$ belongs. These points seem to consist of a curve and the origin of this curve coincides with a subharmonic 4 of the second fundamental frequency $\omega_2$. This coincidence with the subharmonic $\omega_2/4$ may suggest the existence of an internal resonance in the proximity of both curves. Further attempts using sequential continuation to reveal the internal resonance between the curves were not successful. Indeed, the use of sequential continuation prevents us from confirming confidently the existence of the two distinct branches since the method does distinguish between branches belonging to different continua [33]. In order to affirm this hypothesis, a continuation method capable of resolving internal resonances is required.

A method generally used for detection of internal resonances is the pseudo-arclength continuation [7]. However, the use of pseudo-arclength continuation relies on the tangent to the backbone curve to formulate the next solution along the curve [33]. Due to the nonsmoothness of the motions in NSMs, such tangent cannot be formulated at every point on the branch. In fact, pseudo-arclength continuation was applied in [31, p.9] for
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>(b) Displacement field of solution $q_1$, autonomous NBM motion</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Displacement field of solution $q_2$, autonomous NBM motion</figcaption></figure>

<caption>Figure 9: NSM analysis of the bar with soft support with $k = 0.5$ detected by the NBM. Forced response curves with $\bar{F} = 0.025$ and $c = 0.1, 0.2, \dots, 0.7$. First natural frequency $\omega_1 \approx 0.65$, second natural frequency $\omega_2 \approx 3.3$ and grazing energy $E_1 \approx 5.9 \times 10^{-7}$.</caption>

detection of nonsmooth modes. While some continuous sections of the NSM were detected, pseudo-arclength continuation has failed to reveal internal resonances and could not reveal backbone curves for long ranges of frequencies [31, p.10].

### 6.2.4 Varying-area bar

While modal analysis of the bar of uniform area has been the subject of both analytical or numerical analysis [11, 13, 31], it required knowledge of the d'Alembert function or Green's function or characteristic lines. Here, the NBM allows for modal analysis of the varying-area bar since it allows detection of periodic solutions in the autonomous varying-area bar prone to contact. For the varying area bar in Figure 1, several area functions $A(x)$ were considered. In order to simplify the discussion for the remainder of this section, the following terminology is introduced to distinguish between the investigated models:

$$ \textbf{heav-bar} \quad A(x) = 1 - 0.5\Theta(x - 0.5) \equiv A_h(x) \qquad (69) $$

$$ \textbf{lin-bar} \quad A(x) = 1 - x/2 \equiv A_l(x) \qquad (70) $$

$$ \textbf{quad-bar} \quad A(x) = 0.5x(2-x) \equiv A_q(x) \qquad (71) $$

where $\Theta(x)$ denotes the Heaviside function, and the heav-bar hence exhibits two cross-sectional areas: $A_0$ for $\bar{x} \in [0, 0.5L]$ and $0.5A_0$ for $\bar{x} \in (0.5L, 1]$. Furthermore, it is worth noting that all models consist of bars with decreasing areas such that $A(0) = 1$ and $A(1) = 0.5$. The corresponding NSM is illustrated in Figure 10(a). The detected NSMs of the lin- and quad-bars are illustrated in Figures 10(b) and 10(c), respectively. Indeed, the NSMs obtained by the NBM coincide the resonant points in the forced-response diagrams as expected from theory of nonlinear modal analysis [1, 7].
<figure><img src="image_18.png" /><figcaption>(a) heav-bar ($A(x) = A_h(x)$). $\omega_1 \approx 1.44$, $\omega_3 \approx 6.54$ and $E_1 \approx 4.43 \times 10^{-7}$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) lin-bar ($A(x) = A_l(x)$). $\omega_1 \approx 1.44$, $\omega_3 \approx 6.72$ and $E_1 \approx 4.98 \times 10^{-7}$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) quad-bar ($A(x) = A_q(x)$). $\omega_1 \approx 1.35$, $\omega_3 \approx 6.32$ and $E_1 \approx 4.31 \times 10^{-7}$</figcaption></figure>

<caption>Figure 10: Backbone and forced response curves ($\bar{F} = 0.05$) for varying area bar. NSMs detected by the NBM with 20 elements of quadratic shape functions. The depicted points in sub-figures (a), (b) and (c) correspond respectively to solutions *b*, *c*, and *d* investigated in Figure 11(b).</caption>

The effect of area variation on the modal space of the bar in unilateral contact is of interest. In contrast to the conclusions from the bar with soft support in Section 6.2.3, no internal resonances were detected for the cases of the varying area bar around the subharmonic $\omega_3/4$ or other sub-harmonics within spectrum of the backbone curve. Next, Figure 11 illustrates the backbone curves of all varying area models and that of the uniform area investigated in Section 6.2.2. In Figure 11(a), the backbones of the varying area bars exhibit higher energies for the same normalized frequencies when compared to the uniform area bar (investigated in Section 6.2.2). In other terms, the behaviour of the varying area bars can be characterized as “softer” in relation to the uniform area bar. Along the same line, it is noted that while the lin- and quad-bars exhibit a similar stiffening pattern, the heav-bar is characterized by the softest stiffening. Thus, it is indicative that the varying area function affects the stiffening behaviour and, in turn, the range size of resonant frequencies. Next, in Figure 11(b), comparison of the varying-area bar NSMs with respect to their true (not normalized) frequencies shows that the area variation causes backward shifts in the backbone curves towards lower frequency ranges. We note that the backward shifts in Figure 11(b) originates in the vertical line portions of the backbone curves below normalized energy 1. These vertical line portions depict linear modes of vibrations where no contact occurs and $u_N(t) < g$ throughout the whole duration of motion. Since the sections of linear vibration modes of the varying area bars are shifted backwards in frequency, we conjure that the backward shifts are a consequence of the area variation's effect on the linear mode shapes of the bar rather than the contact dynamics introduced to the system.

At last, specific motions on the NSMs of the heav-bar, quad-bar and lin-bar are plotted in Figure 12. It is noted that the NSM motions of the varying area bar models, while qualitatively similar to the NSM motions of the uniform area bar depicted in Figure 8(d) and 8(c), exhibit piecewise nonlinear displacement fields in space-time rather than piecewise-linear displacement field of the uniform area bar.

While the results presented here are novel, a more in depth and more conclusive analysis in regard to the effect of different area variation is a subject of future research.
<figure><img src="image_21.png" /><figcaption>Figure 11: Backbone curves of different area cantilever bars corresponding to Figures 8(a) and 10.</figcaption></figure>

# 7 Scope and limitations of NBM

In this section, we provide a summary of all limitations to NBM and include possible spectrum of applications of NBM to other structures.

The limitations of NBM can be summarized as follows:

1. The choice of shape functions to approximate $u(x, t)$, in the NBM, must be able to exhibit non-zero and negative stress $A(1)u_x(1, t)$. This limitation is common to the Nitsche method and MRM as well. See Section 4.1 for further elaboration.

2. The NBM relies on distinction between internal nodes and nodes prone to contact. Furthermore, those are fixed in throughout the motion and must be known a priori. See Section 4.2 for context and further discussion.

3. The NBM formulation does not cover the case of grazing motions as they were beyond the scope of this paper. Currently, the rules for implementing active or inactive contact are (45) and (46). It is suspected that additional rules are required to allow NBM to capture grazing motions. This is further discussed in Section 4.4.

4. The acceleration term $\ddot{u}_N(t)$ exhibits a Dirac-delta function at the switch between inactive and active contact by definition (obtained by differentiating Equation (51)). The Dirac-delta does not appear in the NBM-ODE formulation and its effects were yet explored. Further elaboration on the matter is given at the last paragraph of Section 4.4.

5. We note that application of NBM for problems involving tangential forces (such as friction [9]) is currently under question. The NBM applies only on the normal direction of contact forces and has been investigated for the one dimensional problem of the bar. Although, research on the NBM for two-dimensional problems in space is ongoing. There, the linear complementarity problem resulting from the NBM implementation of the Signorini conditions is more intricate and sets additional restrictions on the choice of shape functions than the restriction presented here in the context of the one-dimensional NBM. Along the same line, it is assumed that conditions on the tangential direction will add more restrictions on the shape functions. Thus, application of NBM to problems involving friction or tangential forces still remains in question and necessitates a thorough examination.

Together with the presented limitations above, the scope of NBM applications remains sufficiently large. Here, we bring forth a few possible applications of NBM that are beyond the scope of this paper:

1. The NBM can be applied for contact problems involving the one-dimensional bar with varying Young's modulus $E$, density $\rho$, or area $A(x)$ and with any boundary conditions on the non-contacting end (so long that these can be implemented via classical finite-element).

2. The NBM can be applied for contact problems in multidimensional structures. Extension of the NBM to the two-dimensional case has been successfully obtained by the authors. Application of NBM to nonsmooth modal analysis of multidimensional structures is currently ongoing.
<figure><img src="image_22.png" /><figcaption>(a) Displacement field solution b</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Velocity field solution b</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(c) Displacement field solution c</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) Velocity field solution c</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(e) Displacement field solution d</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(f) Velocity field solution d</figcaption></figure>

<caption>Figure 12: Comparison of selected autonomous NBM motions from the backbone curves of all varying-area bars. The locations of the selected motions on the FEP is noted in Figures 10 and 11(b).</caption>

# 8 Conclusion

The nodal boundary method, for treatment of Signorini boundary conditions in the framework of FEM, was presented. The method was developed for nonsmooth modal analysis purposes entailing the detection of periodic solutions to the autonomous Signorini problem. Compared to application of Newton impact law in FEM or WFEM, the resulting ODE from the NBM formulation for the varying area bar allows for existence of periodic solutions with a continuous sticking phase at contact.

The NBM assumes different approximations of the contacting nodal displacement $u_N(t)$ during inactive and active contact phases. The state $u_N(t)$ is dictated by (1) boundary conditions and (2) nodes that are not prone to contact (internal nodes). While the treatment of active contact is done similarly to classical FEM (clamped condition at end of bar), in the treatment of inactive contact, the homogeneous Neumann boundary condition is enforced in a strong sense such that the approximation of the contact stress exactly vanishes, that is $\bar{\sigma}(t) = 0$,
throughout the entire inactive contact duration. The two associated approximations of the quantity $u_N(t)$ can be seen as constituting two distinct sets of shape functions. The residuals for the inactive and active motion approximations are then projected onto their respective set of shape functions to form two distinct ODEs. The Signorini problem is then solved by switching between the sets of shape functions both in trial and test functions. At last, the full NBM-ODE consists of a nonsmooth set of equations which can be solved either analytically or numerically.

Moreover, nonsmooth modal analysis via NBM resulted in valid backbone curves aligning with resonances of forced-response diagrams. These results were obtained for three cases: the cantilever bar of uniform area, the cantilever bar with soft support, and the cantilever bar of varying area.

The results presented for the uniform area bar and the bar with soft support have agreed with previous research on the topic. Furthermore, the NBM has allowed for characterization of two distinct NSMs for the bar with soft support. The two distinct curves seem to relate via an internal resonance as one of the curves originates from a sub-harmonic of the motion. However, affirmation of this result could not be achieved with sequential continuation and the detection of internal resonances in nonsmooth modes is subject for future research. Furthermore, application of the NBM to discover the modal space of the bar of varying area in unilateral contact has been proven successful and the results show good agreement with the forced-response curves.

# 9 Supplementary Material

Excerpts of scripts and algorithms used to perform the analysis and generate this manuscript are available on Zenodo [34].

# References

[1] A. Thorin and M. Legrand. Nonsmooth Modal Analysis: From the Discrete to the Continuous Settings. In *Advanced Topics in Nonsmooth Dynamics: Transactions of the European Network for Nonsmooth Dynamics*, pages 191–234. Springer, 2018. [DOI], [HAL] a.

[2] M. Legrand, S. Junca, and S. 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. [DOI], [HAL] a.

[3] E.H. Moussi, S. Bellizzi, B. Cochelin, and I. Nistor. Nonlinear Normal Modes of a Two Degrees-of-Freedom Piecewise Linear System. *Mechanical Systems and Signal Processing*, 64-65:266–281, 2015. [DOI], [HAL] a.

[4] S. Peter, F. Schreyer, and R. Leine. A Method for Numerical and Experimental Nonlinear Modal Analysis of Nonsmooth Systems. *Mechanical Systems and Signal Processing*, 120:793–807, 2019. [DOI].

[5] G. Kerschen. Definition and Fundamental Properties of Nonlinear Normal Modes. In *Modal Analysis of Nonlinear Mechanical Systems*, volume 555, pages 1–46. Springer, 2014. [DOI].

[6] S. Shaw. Invariant Manifold Representations of Nonlinear Modes of Vibration. In G. Kerschen, editor, *Modal Analysis of Nonlinear Mechanical Systems*, volume 555, pages 47–74. Springer, 2014. [DOI].

[7] M. Peeters, R. Viguié, G. Sérandour, G. Kerschen, and J.-C. Golinval. Nonlinear Normal Modes, Part II: Toward a practical computation using numerical continuation techniques. *Mechanical Systems and Signal Processing*, 23(1):195–216, 2009. [DOI], [HAL] a.

[8] G. James, V. Acary, and F. Pérignon. Periodic Motions of Coupled Impact Oscillators. In R.I. Leine, V. Acary, and O. Brüls, editors, *Advanced Topics in Nonsmooth Dynamics: Transactions of the European Network for Nonsmooth Dynamics*, pages 93–134. Springer, 2018. [DOI], [HAL] a.

[9] V. Yastrebov. *Numerical Methods in Contact Mechanics*. John Wiley & Sons, 2013. [DOI].

[10] B. Brogliato. *Nonsmooth Mechanics: Models, Dynamics and Control*. Springer, London, second ed. edition, 1999. [DOI].

[11] D. Urman, M. Legrand, and S. Junca. D'Alembert Function for Exact Non-Smooth Modal Analysis of the Bar in Unilateral Contact. *Nonlinear Analysis: Hybrid Systems*, 43:101115, 2021. [DOI], [HAL] a.

[12] C. Yoong, A. Thorin, and M. Legrand. Nonsmooth Modal Analysis of an Elastic Bar Subject to a Unilateral Contact Constraint. *Nonlinear Dynamics*, 91(4):2453–2476, 2018. [DOI], [HAL] a.

[13] C. Yoong. *Nonsmooth Modal Analysis of a Finite Linear Elastic Bar Subject to Unilateral Contact Constraint*. PhD thesis, McGill University, 2018. [URL] a.

[14] D. Doyen, A. Ern, and S. Piperno. Time-Integration Schemes for the Finite Element Dynamic Signorini Problem. *SIAM Journal on Scientific Computing*, 33(1):223–249, 2011. [DOI], [HAL] a.

[15] H. Khenous, P. Laborde, and Y. Renard. Mass Redistribution Method for Finite Element Contact Problems in Elastodynamics. *European Journal of Mechanics - A/Solids*, 27(5):918–932, 2008. [DOI], [HAL] a.
[16] F. Chouly, M. Fabre, P. Hild, R. Mlika, J. Pousin, and Y. Renard. An Overview of Recent Results on Nitsche's Method for Contact Problems. In *Geometrically Unfitted Finite Element Methods and Applications*, pages 93–141. Springer, 2017. [DOI], [HAL] a.

[17] P. Wriggers and G. Zavarise. *Computational Contact Mechanics*, chapter 6. Wiley, 2004. [DOI].

[18] C. Hager, S. Hüeber, and B. Wohlmuth. A Stable Energy-Conserving Approach for Frictional Contact Problems Based on Quadrature Formulas. *International Journal for Numerical Methods in Engineering*, 73(2):205–225, 2008. [DOI].

[19] C. Hager and B. Wohlmuth. Analysis of a Space-Time Discretization for Dynamic Elasticity Problems Based on Mass-Free Surface Elements. *SIAM Journal on Numerical Analysis*, 47(3):1863–1885, 2009. [DOI], [HAL] a.

[20] C. Bertrand. Periodic Solutions of a One-Dimensional Elastic Bar Subject to a Unilateral Constraint. Technical report, ENTPE, 2020. [HAL] a.

[21] R. J. LeVeque. *Finite Volume Methods for Hyperbolic Problems*, volume 31. Cambridge university press, 2002. [DOI].

[22] T. Hughes. *Finite Element Method - Linear Static and Dynamic Finite Element Analysis*. Dover Publications, 2000.

[23] V. Krys'ko, J. Awrejcewicz, and G. Narkaitis. Nonlinear Vibration and Characteristics of Flexible Plate-Strips with Non-Symmetric Boundary Conditions. *Communications in Nonlinear Science and Numerical Simulation*, 11(1):95–124, 2006. [DOI].

[24] Y. Modarres-Sadeghi and M. Païdoussis. Nonlinear Dynamics of Extensible Fluid-Conveying Pipes, Supported at Both Ends. *Journal of Fluids and Structures*, 25(3):535–543, 2009. [DOI].

[25] M. Païdoussis. *Pipes Conveying Fluid: Nonlinear and Chaotic Dynamics*, volume 1 of *Fluid-Structure Interactions*, chapter 5, pages 277–414. Academic Press, 1998. [DOI].

[26] J. Boyd. *Chebyshev and Fourier Spectral Methods*, volume 49. Springer, 2 edition, 2001.

[27] L. Evans. *Partial Differential Equations*. American Mathematical Society, 2010.

[28] J. Aitchison and M. Poole. A Numerical Algorithm for the Solution of Signorini Problems. *Journal of Computational and Applied Mathematics*, 94(1):55–67, 1998. [DOI].

[29] S. Krenk. Energy Conservation in Newmark Based Time Integration Algorithms. *Computer Methods in Applied Mechanics and Engineering*, 195(44):6110–6124, 2006. [DOI].

[30] E. Süli. *Finite Element Methods for Partial Differential Equations*. Oxford University Computing Laboratory Oxford, 2002.

[31] T. Lu and M. Legrand. Nonsmooth Modal Analysis via the Boundary Element Method for One-Dimensional Bar Systems. *Nonlinear Dynamics*, October 2021. [DOI], [HAL] a.

[32] V. Acary, O. Bonnefon, M. Brémond, O. Huber, F. Pérignon, and S. Sinclair. An Introduction to Siconos. Technical Report RT-0340, INRIA, 2019. [HAL] a.

[33] E. Allgower and K. Georg. *Numerical Continuation Methods: an Introduction*, volume 13. Springer, 2012. [DOI].

[34] D. Urman and Legrand M. MATLAB Scripts for Nonsmooth Modal Analysis of the bar via Nodal Boundary Method, 2022. [DOI] a.

# A Appendix

## A.1 Proof of $\phi'_N(1) > 0$

**Lemma 1.** In the context of classical finite elements for the PDE (1) with Lagrangian shape functions $\phi_i(x)$ for $i = 1, 2, \dots, N$ consisting of uniformly spaced nodes on loci $x_i = i/N$, the value of $\phi'_N(1)$ depends on the order of shape functions only.

*Proof.* The proof follows from the construction of shape functions in the FEM. In the classical finite element, the structure is divided into elements and each element consists of a set of shape functions that are Lagrange polynomials. These shape functions are local to the element and their definition depends only on the order of the polynomial chosen for this specific element [22]. Thus, the value $\phi'_N(1)$ is dependent only on the order of shape functions used. $\square$

This lemma will be necessary in the generalization of the theorem below to any number of elements and any order of shape functions.

**Theorem 2.** For the classical FE approximation for the PDE (1) with Lagrangian shape functions $\phi_i(x)$ $i = 1, 2, \dots, N$ based on uniformly spaced nodes $x_i = i/N$, the statement

$$ \phi'_{N}(1) > 0 \qquad (72) $$

always holds.
*Proof.* The proof consists of first proving Inequality (72) for the case of a single element by inspecting the exact expression of the Lagrangian function. Then, the proof is expanded to any number of elements and/or shape functions by virtue of Lemma 1.

We start by approximating the bar's displacement using a single element. It follows, then, that $N \ge 1$ stands for the order of the shape function such that

$$ \phi_i(x) = \prod_{\substack{0 \le m \le N \\ m \ne i}} \frac{x - x_m}{x_i - x_m}, \quad x \in [0, 1]. \qquad (73) $$

Then, for $\phi_N(x)$, we can simplify the expression using $x_i = i/N$ to simplify the denominator

$$ \phi_N(x) = \prod_{m=0}^{N-1} \frac{x - x_m}{x_N - x_m} = N^N \prod_{m=0}^{N-1} \frac{x - x_m}{N - m}. \qquad (74) $$

Then, we take the derivative of $\phi_N(x)$ to obtain

$$ \phi'_{N}(x) = N^{N} \sum_{j=0}^{N-1} \frac{1}{N-j} \prod_{\substack{0 \le m \le N-1 \\ m \ne j}} \frac{x-x_m}{N-m}. \qquad (75) $$

To show $\phi'_{N}(1) > 0$, we simply evaluate every term in the expression

$$ \begin{aligned} \phi'_{N}(1) &= N^N \sum_{j=0}^{N-1} \frac{1}{N-j} \prod_{\substack{0 \le m \le N-1 \\ m \ne j}} \frac{1-x_m}{N-m} = N^N \sum_{j=0}^{N-1} \frac{1}{N-j} \prod_{\substack{0 \le m \le N-1 \\ m \ne j}} \frac{1}{N} \frac{N-m}{N-m} \\ &= \left(\frac{N^N}{N^{N-1}}\right) \sum_{j=0}^{N-1} \frac{1}{N-j} = N \sum_{j=0}^{N-1} \frac{1}{N-j} = N \sum_{j=1}^{N} \frac{1}{j}. \end{aligned} \qquad (76) $$

We note that all the quantities presented here are exclusively positive and that their sum and product will also be positive such that $\phi'_{N}(1) > 0$ holds. By virtue of Lemma 1, we may conclude that if $\phi'_{N}(1) > 0$ holds for a single element, it will hold for any number of elements. $\square$

## A.2 Conservation of energy away from instant of switch

We introduce two lemmas on the energy conservation during active and inactive contact motions, away from the moment of switch ($t_s$ such that $S(\mathbf{u}^o(t_s)) = g$).

**Lemma 3.** For the FEM-NBM ODE (35) developed for a cantilever bar in inactive contact conditions, the energy (involving boundary nodes and non reduced matrices) term

$$ 2E(t) = \ddot{\mathbf{u}}^{\top}(t)\mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{u}^{\top}(t)\mathbf{K}\mathbf{u}(t) \qquad (77) $$

is conserved, that is $\dot{E}(t) = 0$ for $S(\mathbf{u}^o(t)) < g$.

*Proof.* This lemma is proven by developing the term $\dot{E}$ from Equation (77) and plugging the solution of Equation (35).

First, we differentiate Equation (77) with respect to $t$, and for symmetric **M** and **K** (as they are in the given in the FE formulation [22]) we obtain

$$ 2\dot{E}(t) = \ddot{\mathbf{u}}^{\top}(t)(\mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{K}\mathbf{u}(t)). \qquad (78) $$

Next, we note that for inactive contact $S(\mathbf{u}^o(t)) < g$, the NBM formulation of $u_N$ in (28) admits $\dot{u}_N = S(\dot{\mathbf{u}}^o(t))$. Under this restriction, the relationship between $\dot{\mathbf{u}}(t)$ and $\dot{\mathbf{u}}^o(t)$ can be described via $\dot{\mathbf{u}}(t) = \mathbf{B}\dot{\mathbf{u}}^o(t)$. Next, substitution of the previous identity into expression (78) admits $\dot{E}(t) = \dot{\mathbf{u}}^{\top}(t)(\mathbf{M}_N\ddot{\mathbf{u}}^o(t) + \mathbf{K}_N\mathbf{u}^o(t)) = 0$, for $S(\mathbf{u}^o(t)) < g$ where the last equality completes this proof by virtue of (35). $\square$
**Lemma 4.** For the FEM-NBM ODE (43) developed for a cantilever bar with non-homogeneous Dirichlet conditions, the energy metric (77) (similar to the energy metric used in Lemma (3)) is conserved, that is $\dot{E}(t) = 0$ for $S(\mathbf{u}^o(t)) > g$.

*Proof.* This theorem is proven by developing the term $\dot{E}(t)$ from Equation (77) and plugging the solution of Equation (43). We note that application of Equation (36) implies that

$$ \dot{u}_N(t) = 0. \qquad (79) $$

We can conclude that $\mathbf{u}(t)$ belongs to the same space as $\mathbf{w}$ in (40), ie

$$ \dot{\mathbf{u}}(t) = \mathbf{B}^d \dot{\mathbf{u}}^o(t). \qquad (80) $$

Thus, we plug-in the identities (40) and (80) into (78) such that $\dot{E}(t) = \dot{\mathbf{u}}^o^\top(t)(\mathbf{M}_D\ddot{\mathbf{u}}^o(t)+\mathbf{K}_D\mathbf{u}^o(t)+g\mathbf{f}_D) = 0$, for $S(\mathbf{u}^o(t)) > g$ where the last equality is obtained by virtue of (43). $\square$

### A.3 Derivation of energy jump at switch

In this section, the jump in energy in NBM solutions (see Section 4.5), $\Delta E$, is developed mathematically to obtain the relation in Equation (54). This jump in energy occurs upon switching between active to inactive contact.

In what follows, the energy jump is derived for the case opening contact (transition from active to inactive contact) and is denoted $\Delta E_{A/I}$. This derivation procedure can be then replicated for the case of closing contact (transition from inactive to active contact) denoted $\Delta E_{I/A}$. In this paper, the term $\Delta E_{I/A}$ is presented but not its derivation. Finally, the general term $\Delta E$ describing the energy jump between active and inactive contact (regardless of the order of contact phases) is derived from both $\Delta E_{A/I}$ and $\Delta E_{I/A}$.

To develop the term $\Delta E_{A/I}$, we first denote $t_s$ as the time of switch at opening contact. At opening contact, the stress at the contacting end of the bar, $\sigma(t_s)$, is zero

$$ \sigma(t_s) = A(1)\phi'_N(1)(g - S(\mathbf{u}^o(t_s))) = 0. \qquad (81) $$

Furthermore, the following statements hold

1. The bar is in contact with the wall during the moment of switch (see Equation (50))

$$ u_N(t_s) = g. \qquad (82) $$

2. The velocity of the contact node during active contact (just before the switch) is zero (see Equation (51))

$$ \dot{u}_N(t_s^-) = 0. \qquad (83) $$

3. The derivative of the stress function at the tip of the bar is positive

$$ \dot{\sigma}(t_s^-) = A(1)\phi'_N(1)(g - S(\dot{\mathbf{u}}^o(t_s))) > 0 \qquad (84) $$

since the stress before the switch must be strictly negative during active contact (by Equation (52) and assuming non-grazing motion), zero at the switch (see Equation (81)), and $\sigma(t)$ is continuous at all times (see Section 4.4). From Equation (84) we can also deduce that

$$ S(\dot{\mathbf{u}}^o(t_s)) < 0 \qquad (85) $$

since $A(1) > 0$ by definition, and $\phi'_N(1) > 0$ holds as proven in Appendix A.1.

4. After the switch, the system is in inactive contact, and the velocity of the contact node is negative by virtue of Equations (51) and (85):

$$ \dot{u}_N(t_s^+) = S(\dot{\mathbf{u}}^o(t_s)) < 0. \qquad (86) $$

In fact, it is expected that the velocity after the switch is negative as after the transition from active to inactive contact, the bar should separate from the wall (for non-grazing motions).
To simplify the notations in the upcoming equations, we present the energy of the system (53) using block matrices:

$$2E(t) = [\dot{\mathbf{u}}^o(t)^\top \quad \dot{u}_N(t)] \begin{bmatrix} \mathbf{M}_{oN} & \mathbf{M}_{oN} \\ \mathbf{M}_{oN}^\top & M_{NN} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{u}}^o(t_s) \\ \dot{u}_N(t) \end{bmatrix} + [\mathbf{u}^o(t_s)^\top \quad u_N(t)] \begin{bmatrix} \mathbf{K}_{oN} & \mathbf{K}_{oN} \\ \mathbf{K}_{oN}^\top & K_{NN} \end{bmatrix} \begin{bmatrix} \mathbf{u}^o(t_s) \\ u_N(t) \end{bmatrix}. \quad (87)$$

Then, we can describe the system's energy just before the moment of switch (during active contact) by plugging Equations (82) and (83) into Equation (87)

$$2E(t_s^-) = [\dot{\mathbf{u}}^o(t_s)^\top \quad 0] \begin{bmatrix} \mathbf{M}_{oN} & \mathbf{M}_{oN} \\ \mathbf{M}_{oN}^\top & M_{NN} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{u}}^o(t_s) \\ 0 \end{bmatrix} + [\mathbf{u}^o(t_s)^\top \quad g] \begin{bmatrix} \mathbf{K}_{oN} & \mathbf{K}_{oN} \\ \mathbf{K}_{oN}^\top & K_{NN} \end{bmatrix} \begin{bmatrix} \mathbf{u}^o(t_s) \\ g \end{bmatrix}. \quad (88)$$

Next, the energy of the system at $t_s^+$, presuming that the bar is in inactive contact motion, is given by plugging Equations (82) and (86) into Equation (87):

$$2E(t_s^+) = [\dot{\mathbf{u}}^o(t_s)^\top \quad S(\dot{\mathbf{u}}^o(t_s))] \begin{bmatrix} \mathbf{M}_{oN} & \mathbf{M}_{oN} \\ \mathbf{M}_{oN}^\top & M_{NN} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{u}}^o(t_s) \\ S(\dot{\mathbf{u}}^o(t_s)) \end{bmatrix} + \dots \\ \begin{bmatrix} \mathbf{u}^o(t_s)^\top \quad g \end{bmatrix} \begin{bmatrix} \mathbf{K}_{oN} & \mathbf{K}_{oN} \\ \mathbf{K}_{oN}^\top & K_{NN} \end{bmatrix} \begin{bmatrix} \mathbf{u}^o(t_s) \\ g \end{bmatrix}, \quad S(\dot{\mathbf{u}}^o(t_s)) < 0. \quad (89)$$

Then, the term $\Delta E_{A/I}$ is given as follows (where calculations where omitted for sake of conciseness)

$$\Delta E_{A/I} = E(t_s^+) - E(t_s^-) = S(\dot{\mathbf{u}}^o(t_s)) [\mathbf{M}_{oN}^\top \dot{\mathbf{u}}^o(t_s) + \frac{1}{2} M_{NN} S(\dot{\mathbf{u}}^o(t_s))], \quad S(\dot{\mathbf{u}}^o(t_s)) < 0. \quad (90)$$

Similarly, if we would have started our derivation with the assumption that $t_s$ denotes a instance of closing contact, we would find that the energy jump satisfies

$$\Delta E_{I/A} = -S(\dot{\mathbf{u}}^o(t_s)) [\mathbf{M}_{oN}^\top \dot{\mathbf{u}}^o(t_s) + \frac{1}{2} M_{NN} S(\dot{\mathbf{u}}^o(t_s))], \quad S(\dot{\mathbf{u}}^o(t_s)) > 0. \quad (91)$$

At last, we derive the general term for the energy jump $\Delta E$ by combining Equations (90) and (91)

$$\Delta E = -|S(\dot{\mathbf{u}}^o(t_s))| [\mathbf{M}_{oN}^\top \dot{\mathbf{u}}^o(t_s) + \frac{1}{2} M_{NN} S(\dot{\mathbf{u}}^o(t_s))]. \quad (92)$$

This term is presented using the elements of $\mathbf{M}$ in Equation (54).

**CRediT authorship contribution statement** David Urman: Conceptualization, Methodology, Software, Validation, Formal analysis, Writing – original draft, Writing – review & editing, Visualization. Mathias Legrand: Conceptualization, Formal analysis, Validation, Writing – original draft, Writing – review & editing, Supervision, Project administration, Funding acquisition.

**Declaration of competing interest** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Acknowledgement** The authors gratefully acknowledge the financial support by the Natural Sciences and Engineering Research Council of Canada through the Discovery Grant program, award number RGPIN-2018-04497.