# Periodic solutions of a one-dimensional elastic bar subject to a unilateral constraint

Charlélie Bertrand

► To cite this version:

Charlélie Bertrand. Periodic solutions of a one-dimensional elastic bar subject to a unilateral constraint. [Technical Report] ENTPE, Vaulx en velin. 2020. hal-02897846

HAL Id: hal-02897846

https://hal.science/hal-02897846

Submitted on 15 Jul 2020

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
A technical note on:

Periodic solutions of a one-dimensional
elastic bar subject to a unilateral
constraint

Charlélie BERTRAND

Mathias LEGRAND

July 15, 2020
# Contents

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Analytical framework for periodic solutions of one-dimensional unilateral contact problems</td><td>4</td></tr><tr><td>1.1</td><td>Integral form of d'Alembert solution</td><td>4</td></tr><tr><td>1.1.1</td><td>System of interest</td><td>4</td></tr><tr><td>1.1.2</td><td>Analytical solution</td><td>5</td></tr><tr><td>1.1.3</td><td>Enforcement of boundary conditions</td><td>5</td></tr><tr><td>1.1.4</td><td>Admissibility conditions</td><td>6</td></tr><tr><td>1.2</td><td>Motions with one closing contact per period</td><td>7</td></tr><tr><td>1.2.1</td><td>Periodicity condition</td><td>7</td></tr><tr><td>1.2.2</td><td>Necessary conditions for existence of a continuum of solutions</td><td>7</td></tr><tr><td>1.2.3</td><td>First step to solve Equation (1.24)</td><td>8</td></tr><tr><td>1.2.4</td><td>Reformulation</td><td>8</td></tr><tr><td>1.2.5</td><td>Sequences matching of φ function</td><td>9</td></tr><tr><td>1.2.6</td><td>Admissibility</td><td>10</td></tr><tr><td>1.2.7</td><td>Results and discussion</td><td>12</td></tr><tr><td>1.3</td><td>Possible developments</td><td>14</td></tr><tr><td>1.3.1</td><td>Towards exhaustive continuum of solutions</td><td>14</td></tr><tr><td>1.3.2</td><td>Framework for n boundary switches per period</td><td>14</td></tr><tr><td>1.3.3</td><td>Checking periodicity for inferred solutions</td><td>15</td></tr><tr><td>1.3.4</td><td>Extension to Robin boundary condition at x = 0</td><td>16</td></tr><tr><td>2</td><td>On a stability analysis</td><td>18</td></tr><tr><td>2.1</td><td>A stress-velocity formulation</td><td>18</td></tr><tr><td>2.1.1</td><td>Hyperbolic system of conservation law</td><td>18</td></tr><tr><td>2.1.2</td><td>Switching mechanism</td><td>19</td></tr><tr><td>2.1.3</td><td>Mapping for a given initial condition</td><td>19</td></tr><tr><td>2.2</td><td>Stability under small perturbations</td><td>21</td></tr><tr><td>2.2.1</td><td>Definition</td><td>21</td></tr><tr><td>2.2.2</td><td>Assumptions</td><td>21</td></tr><tr><td>2.2.3</td><td>Linearized mapping for an initial perturbation</td><td>22</td></tr><tr><td>2.2.4</td><td>Linearized "First Return Map" operator</td><td>24</td></tr><tr><td>2.3</td><td>Numerical routine for eigenvalues computation</td><td>24</td></tr><tr><td>2.3.1</td><td>Discretized mapping for an initial perturbation</td><td>24</td></tr><tr><td>2.3.2</td><td>Results and discussions</td><td>25</td></tr><tr><td colspan="2">Appendices</td><td>28</td></tr><tr><td colspan="2">A Notations</td><td>29</td></tr><tr><td colspan="2">B Computation of the linearized First Return Map L</td><td>31</td></tr></tbody></table>
<table><tr><td>C</td><td>Discretization of the linearized First Return Map L</td><td>32</td></tr><tr><td>C.1</td><td>Discretized first return map</td><td>32</td></tr><tr><td>C.1.1</td><td>Space discretization</td><td>32</td></tr><tr><td>C.1.2</td><td>Discretized equivalents of C and F</td><td>32</td></tr><tr><td>C.1.3</td><td>Expression of the B matrix</td><td>33</td></tr><tr><td>C.1.4</td><td>Discrete Derivation</td><td>33</td></tr><tr><td>C.1.5</td><td>Discretization of integral parts of the operator</td><td>34</td></tr><tr><td>C.1.6</td><td>Discretized operator</td><td>35</td></tr><tr><td>C.2</td><td>Convergence of the discretized operator</td><td>36</td></tr><tr><td>C.2.1</td><td>Consistency of WFEM</td><td>37</td></tr><tr><td>C.2.2</td><td>Stability of WFEM</td><td>38</td></tr></table>
# Introduction

In the framework of non-smooth modal analysis (mode shapes and corresponding frequencies) of structural systems, earlier investigations by Carlos Yoong [12] allowed to spot autonomous periodic motions of a one-dimensional elastic bar with a unilateral contact constraint on the boundary. This analysis was achieved using the Wave Finite Element Method (WFEM) which is a shock capturing numerical method. This approach allows to capture accurately discontinuous stress and velocity wavefronts without any numerical dispersion [11]. In order to further investigate the problem, the same problem is solved in an analytical manner in order to show the existence of continuum of solutions. New developments within this framework lead to necessary conditions for a solution to be a periodic solution. Moreover, solving this problem analytically is the only way to ensure that we have an exhaustive description of the solutions.

This report also displays a numerical procedure to investigate the orbital stability in the sense of Poincaré of these periodic motions. We attempted to exhibit a linear First Return mapping in a continuous framework. We also introduce some details in the formulation of the WFEM which are specific to a stress-velocity formulation. In appendix are detailed developments, notations and explanations about scripts that have been developed.
# Chapter 1

## Analytical framework for periodic solutions of one-dimensional unilateral contact problems

### 1.1 Integral form of d'Alembert solution

#### 1.1.1 System of interest

The system of interest, depicted in Figure 1.1, is a homogeneous elastic bar of length $L > 0$ and Young modulus $E > 0$. Working in the framework of linear elasticity, we assume that both the cross section of the bar $S > 0$ and its mass per unit volume $\rho > 0$ are constant. At rest, the bar is clamped rigidly at its left extremity and its right extremity is at a distance $g_0$ from a rigid wall. We assume that the boundary conditions at $x = L$ can switch during the motion from a homogeneous Neumann boundary condition to inhomogeneous Dirichlet boundary condition. The boundary switching are ruled by Signorini boundary conditions.

<figure><img src="image_2.png" /><figcaption>Figure 1.1: Elastic bar subject to a unilateral constraint</figcaption></figure>

The full formulation of the problem reads:

* Equilibrium of the bar:

$$ \forall x \in (0, L), \quad \forall t \ge 0, \quad \partial_t^2 u - c^2 \partial_x^2 u = 0, \quad c = \sqrt{E/\rho}; \qquad (1.1) $$

* Boundary conditions at $x = 0$:

$$ \forall t \ge 0, \quad u(0, t) = 0; \qquad (1.2) $$

* Initial conditions:

$$ \forall x \in (0, L), \quad u(x, 0) = u_0(x), \quad v(x, 0) = v_0(x); \qquad (1.3) $$

* Complementary conditions (Signorini conditions) at $x = L$:

$$ \forall t \ge 0, \quad g(t) = u(L, t) - g_0 \ge 0, \quad \sigma(L, t) \le 0, \quad \sigma(L, t)g(t) = 0. \qquad (1.4) $$
### 1.1.2 Analytical solution

The solution to Equation (1.1) is known to be a superposition of a forward and a backward traveling wave [8, 7] of the form

$$ \forall x \in \mathbb{R}, \quad \forall t \ge 0, \quad u(x, t) = f(ct + x) + h(ct - x). \qquad (1.5) $$

The Dirichlet boundary condition at $x = 0$ implies $h = -f$ and thus

$$ \forall x \in \mathbb{R}, \quad \forall t \ge 0, \quad u(x, t) = f(ct + x) - f(ct - x). \qquad (1.6) $$

We assume that the function $f$ is piecewise continuous on $\mathbb{R}$. It can thus be recast as an integral of a function $\phi$ which is defined almost everywhere in $\mathbb{R}$

$$ \forall x \in \mathbb{R}, \quad \forall t \ge 0, \quad u(x, t) = \int_{ct-x}^{ct+x} \phi(s) \, ds. \qquad (1.7) $$

It implies that $\forall x \in \mathbb{R}$, $\lim_{s \to x^{-}} \phi(s)$ and $\lim_{s \to x^{+}} \phi(s)$ exist. Accordingly, stress and velocity can be derived as

$$ v(x, t) = c\phi(ct + x) - c\phi(ct - x), \qquad (1.8) $$

$$ \sigma(x, t) = E\phi(ct + x) + E\phi(ct - x). \qquad (1.9) $$

Hence, we define $\phi$ through the provided boundary conditions. It guarantees that the motion is determined almost everywhere if we know $\phi$.

### 1.1.3 Enforcement of boundary conditions

The developments exposed in this section would not have been possible without the wise advice by Pierre DELEZOIDE who provided Equation (1.24).

We now consider that the right extremity of the bar, $x = L$, switches from a Neumann boundary condition to Dirichlet boundary condition:

**Neumann boundary condition at $x = L$** When the right extremity of the bar is free, the Neumann boundary condition reads $\phi(s + L) = -\phi(s - L)$, $\forall s \in \mathbb{R}$, which implies that $\phi$ is $2L$-antiperiodic during free-flight phases;

**Dirichlet boundary condition at $x = L$** When contact is closed, the right extremity of the bar is essentially clamped. Accordingly, this is reflected on $\phi$ as a homogeneous Dirichlet boundary condition of the form $\phi(s + L) = \phi(s - L)$, $\forall s \in \mathbb{R}$. This implies that $\phi$ is $2L$-periodic during contact phases.

**Switching mechanism**

We assume that the bar leaves the wall at $t = 0$. During this first phase, the motion is fully determined by $\phi_0$. The bar then closes the gap with the wall at $t = t_f$ and rests on the wall: the motion is now described by $\phi_1$. Since both $\phi_0$ and $\phi_1$ stem from the same initial condition, then

$$ \forall s \in [-L, L], \quad \phi_0(s + ct_f^-) = \phi_1(s + 0^+). \qquad (1.10) $$

However, $\phi_0(\bullet + t_f^-)$ is $2L$-antiperiodic and $\phi_1(\bullet + 0^+)$ is $2L$-periodic. To extend this relationship between $\phi_0$ and $\phi_1$ over $\mathbb{R}$, we consider the function $\epsilon(s) = 1$, $\forall s \in [-L, L]$ and $\epsilon$ is $2L$-antiperiodic. This function is shown in Figure 1.2. From Equation (1.10) and the $\epsilon$ function, the following is deduced:

$$ \forall s \in \mathbb{R}, \quad \phi_1(s) = \epsilon(s)\phi_0(s + ct_f). \qquad (1.11) $$

This relation still holds for the other kind of boundary switching.
<figure><img src="image_3.png" /><figcaption>Figure 1.2: $\epsilon$ function on $[-5L, 5L]$</figcaption></figure>

### 1.1.4 Admissibility conditions

We assume that the Signorini complementarity conditions are satisfied during the motion, that is $\forall t \ge 0$, $g(t) = u(L, t) - g_0 \ge 0$, $\sigma(L, t) \le 0$, $\sigma(L, t)g(t) = 0$. This can be solely enforced through $\phi$.

#### Free-phase

If $g(t) \ge 0$, then by using Equation (1.7), $\forall t \in [0, t_f]$

$$ 
\begin{aligned}
g(t) &= g_0 - \int_{ct-L}^{ct+L} \phi_0(s) \, ds = g_0 - \int_{ct-L}^{-L} \phi_0(s) \, ds - \int_{-L}^{L} \phi_0(s) \, ds - \int_{L}^{ct+L} \phi_0(s) \, ds \\
&= -\int_{ct-L}^{-L} \phi_0(s) \, ds - \int_{L}^{ct+L} \phi_0(s) \, ds = \int_{ct-L}^{-L} \phi_0(s+2L) \, ds - \int_{L}^{ct+L} \phi_0(s) \, ds \\
&= \int_{ct+L}^{L} \phi_0(s) \, ds - \int_{L}^{ct+L} \phi_0(s) \, ds = 2 \int_{ct+L}^{L} \phi_0(s) \, ds
\end{aligned} 
\quad (1.12) $$

To be admissible, $\phi_0$ must satisfy the condition

$$ \forall t \in [0, t_f], \quad \int_{L}^{ct+L} \phi_0(s) \, ds \le 0. \quad (1.13) $$

Moreover, to ensure a closing contact at $t = t_f$, $\phi_0$ must satisfy

$$ \int_{L}^{ct_f+L} \phi_0(s) \, ds = 0. \quad (1.14) $$

#### Contact phase

If $\sigma(L, t) \le 0$, then via Equations (1.9) and (1.11) and Figure 1.2 together with the fact that $\phi_0$ is 2L-antiperiodic

$$ 
\begin{aligned}
\sigma(L,t) &= E\phi_1(ct+L) + E\phi_1(ct-L) \\
&= E\epsilon(ct+L)(\phi_0(ct_f+ct+L) - \phi_0(ct_f+ct-L)) \\
&= 2E\epsilon(ct+L)\phi_0(ct_f+ct+L).
\end{aligned} 
\quad (1.15) $$

To be admissible, $\phi_0$ must also satisfy

$$ \forall t \in [0, t_{dc}], \quad \epsilon(ct+L)\phi_0(ct_f+ct+L) \le 0. \quad (1.16) $$
Moreover, to ensure the contact ends at $t = t_{dc}$, $\phi_0$ must satisfy

$$ \epsilon(ct_{dc} + L)\phi_0(ct_f + ct_{dc} + L) = 0. \qquad (1.17) $$

According to Equation (1.16), since $\epsilon(ct_{dc} + L) = -1$, the following equivalence holds:

$$ \begin{gathered} \forall t \in [0, t_{dc}], \quad \phi_0(c(T - t_{dc}) + ct + L) \ge 0 \\ \iff \forall x \in [L + c(T - t_{dc}), L + cT], \quad \phi_0(x) \ge 0. \end{gathered} \qquad (1.18) $$

## 1.2 Motions with one closing contact per period

We are now interested in finding periodic motions for the problem described in the previous section. We draw our focus on motions with one free phase and one closed phase per period.

### 1.2.1 Periodicity condition

We assume that the motion starts at an instant of release, so that $\phi$ has to be $2L$-antiperiodic. We will refer to the extension on the real axis of $\phi$, defined on $[-L, L]$ only, by writing $\phi_0$. We assume contact to close at $t = t_f$, ie we want the boundary condition to switch at this time. We use the switching tool defined earlier via the $\epsilon$ function, and get $\phi_1$ such that

$$ \forall x \in \mathbb{R}, \quad \forall t \le 0, \quad \phi_1(x) = \epsilon(x)\phi_0(ct_f + x). \qquad (1.19) $$

When release occurs, ie the boundary condition switches again, we get $\phi_2$ such that

$$ \forall x \in \mathbb{R}, \quad \phi_2(x) = \epsilon(x)\phi_1(ct_{dc} + x). \qquad (1.20) $$

Equations (1.19) and (1.20) are combined and T-periodicity is enforced on $\phi_0(x)$ as follows:

$$ \phi_0(x) = \phi_2(x), \qquad (1.21) $$

$$ \phi_0(x) = \epsilon(x)\phi_1(ct_{dc} + x), \qquad (1.22) $$

$$ = \epsilon(x)\epsilon(ct_{dc} + x)\phi_0(ct_f + ct_{dc} + x), \qquad (1.23) $$

...

$$ \phi_0(x) = \epsilon(x)\epsilon(ct_{dc} + x)\phi_0(cT + x). \qquad (1.24) $$

Equation (1.24) is of utmost importance in the remainder.

### 1.2.2 Necessary conditions for existence of a continuum of solutions

Squaring Equation (1.24) yields $\phi_0^2(x) = \phi_0^2(ct + x)$, $\forall x \in \mathbb{R}$. We want to investigate on the existence of continua of periodic solutions. Since $\phi_0^2$ is $2L$-periodic and also T-periodic, then $\phi_0^2$ should be constant, to allow for possible cases with $T/(2L) \in \mathbb{R} - \mathbb{Q}$ [4]. Indeed, a non constant $\phi_0^2$ would lead to Q-parametrized solutions only. However, we seek for $T \in \mathbb{R}$ then $\phi^2$ is constant. This leads to the fact that every displacement will be piecewise linear and that stress and velocity will be piecewise constant. Without any loss of generality, we assume that $\phi_0^2 = 1$, that is

$$ \forall x \in \mathbb{R}, \quad \phi_0(x) = \pm 1. \qquad (1.25) $$
The function $\phi_0$ will have to be multiplied by a factor, in order to be compatible with the existing gap with the bar at rest, $g_0$. We denote this multiplicative constant by $\alpha_0$. From Equation (1.25), this constant can be derived from the identity

$$g_0 = u(L, 0) = \alpha_0 \int_{-L}^{L} \phi_0(s) \, ds. \qquad (1.26)$$

Interestingly, the case $g_0 = 0$ could lead to infinitely many possible solutions with various levels of energy. The reader has to keep in mind that complex behaviors can arise when this case is reached. We could expect non “piecewise linear” displacement for instance.

### 1.2.3 First step to solve Equation (1.24)

We consider a solution $\phi_0$ to Equation (1.24) and we assume that it has finitely many discontinuities over $[-L, L]$. The full set of discontinuities of a function $f$ over $\mathbb{R}$ is denoted as $\mathcal{D}(f)$. In physical terms, a discontinuity can be seen as a wave front. We assume that $T$ is irrational and that $\phi_{0,[-L,L]} \neq \pm\epsilon \times \epsilon(\bullet + ct_{dc})$. This implies

If $x \in \mathcal{D}(\phi_0) - \{\mathcal{D}(\epsilon) \cup (\mathcal{D}(\epsilon) - \{ct_{dc}\})\}$ and $\phi_0$ solves (1.24) then $(cT+x) \in \mathcal{D}(\phi_0)$. (1.27)

Using Equation (1.24), we can forward and backward iterate this operation a finitely many times until $x + \kappa cT$ belongs to $\{\mathcal{D}(\epsilon) \cup (\mathcal{D}(\epsilon) - \{ct_{dc}\})\} - \mathcal{D}(\phi)$. If we can iterate to infinity, it leads to a contradiction because $cT$ is irrational and thus

$$\exists(n^+, n^-) \in \mathbb{N}^{*2}, \quad \begin{cases} x + n^+ cT \in \{\mathcal{D}(\epsilon) \cup (\mathcal{D}(\epsilon) - c\{t_{dc}\})\} - \mathcal{D}(\phi_0), \\ x - n^- cT \in \{\mathcal{D}(\epsilon) \cup (\mathcal{D}(\epsilon) - c\{t_{dc}\})\} - \mathcal{D}(\phi_0). \end{cases} \qquad (1.28)$$

Since $\mathcal{D}(\epsilon) \cup (\mathcal{D}(\epsilon) - \{ct_{dc}\}) = (2\mathbb{Z} + 1)L \cup (2\mathbb{Z} + 1)L - ct_{dc}$, it leads to two possible cases

$$\exists(p, q) \in \mathbb{N}^* \times \mathbb{Z}, \quad \begin{cases} pcT = 2qL, \\ pcT = 2qL \pm ct_{dc}. \end{cases} \qquad (1.29)$$

Since $T$ is irrational, only Equation (1.29) is possible and

$$\exists(p, q) \in \mathbb{N}^* \times \mathbb{Z}, \quad T = \frac{q}{p} \frac{2L}{c} \pm \frac{t_{dc}}{p}. \qquad (1.30)$$

If we assume that $T \ge t_{dc} \ge 0$, then $q \in \mathbb{N}$.

### 1.2.4 Reformulation

All the previous results on $\phi$ can be recast as

$$\phi_0(x) = \epsilon(x)\epsilon(ct_{dc} + x)\phi_0(cT + x), \quad \forall x \in \mathbb{R} \qquad \text{periodicity condition} \qquad (1.31a)$$

$$\int_{-L}^{ct-L} \phi_0(s) \, ds \ge 0, \quad \forall t \in [0, t_f] \qquad \text{no penetration} \qquad (1.31b)$$

$$\phi_0(x) \le 0, \quad \forall x \in [L - ct_{dc}, L] \qquad \text{compression condition} \qquad (1.31c)$$

$$\int_{-L}^{ct_f-L} \phi_0(s) \, ds = 0 \qquad \text{closure at } t = t_f \qquad (1.31d)$$

$$T = t_f + t_{dc} \qquad (1.31e)$$

$$\exists(p, q) \in \mathbb{N}^* \times \mathbb{N}, \quad T = (q/p)(2L/c) \pm (t_{dc}/p) \qquad (1.31f)$$

together which a wise use of the $2L$-antiperiodicity of $\phi_0$.
### 1.2.5 Sequences matching of φ function

Through Equation (1.24), $\phi_0$ can be expressed as

$$ \forall x \in \mathbb{R}, \quad \phi_0(x) = \begin{cases} A(x) & x \in [-L, L - ct_{\mathrm{dc}}], \\ B(x) & x \in [L - ct_{\mathrm{dc}}, L]. \end{cases} \tag{1.32a} $$

Cases where $cT < ct_{\mathrm{dc}}$ were excluded because physically non-admissible. Others cases are redundant because of the 4L-periodicity of $\phi_0$. Without any loss of generality, we work on A and B defined respectively on $[0, 2L - ct_{\mathrm{dc}}]$ and $[0, ct_{\mathrm{dc}}]$ only.

**Case** $ct_{\mathrm{dc}} \le 2L - ct_{\mathrm{dc}}$

Using Equation (1.24), we can derive a sequence matching for the $\phi_0$ function which is

$$ \bullet \quad ct_{\mathrm{dc}} \le cT \le 2L - ct_{\mathrm{dc}} $$

$$ \begin{aligned} A(x) &= A(x + cT) & x \in [0, 2L - ct_{\mathrm{dc}} - cT] \\ A(x) &= B(x - 2L + ct_{\mathrm{dc}} + cT) & x \in [2L - ct_{\mathrm{dc}} - cT, 2L - cT] \\ A(x) &= -A(x - 2L + cT) & x \in [2L - cT, 2L - ct_{\mathrm{dc}}] \\ B(x) &= A(x + cT - ct_{\mathrm{dc}}) & x \in [cT - ct_{\mathrm{dc}}, cT] \end{aligned} \tag{1.33} $$

$$ \bullet \quad 2L - ct_{\mathrm{dc}} \le cT < 2L $$

$$ \begin{aligned} A(x) &= B(x + cT - 2L + ct_{\mathrm{dc}}) & x \in [0, 2L - cT] \\ A(x) &= -A(x + 2L - cT) & x \in [2L - cT, 2L - ct_{\mathrm{dc}}] \\ B(x) &= A(x + cT - ct_{\mathrm{dc}}) & x \in [0, 2L - cT] \\ B(x) &= B(x - 2L + cT) & x \in [2L - cT, ct_{\mathrm{dc}}] \end{aligned} \tag{1.34} $$

$$ \bullet \quad 2L < cT \le 2L + ct_{\mathrm{dc}} $$

$$ \begin{aligned} A(x) &= -A(x - 2L + cT) & x \in [0, 4L - cT - ct_{\mathrm{dc}}] \\ A(x) &= -B(x - 4L + ct_{\mathrm{dc}} + cT) & x \in [4L - cT - ct_{\mathrm{dc}}, 2L - ct_{\mathrm{dc}}] \\ B(x) &= B(x + cT - 2L) & x \in [0, 2L + ct_{\mathrm{dc}} - cT] \\ B(x) &= -A(x - 2L - ct_{\mathrm{dc}} + cT) & x \in [2L + ct_{\mathrm{dc}} - cT, ct_{\mathrm{dc}}] \end{aligned} \tag{1.35} $$

$$ \bullet \quad 2L + ct_{\mathrm{dc}} \le cT \le 4L - ct_{\mathrm{dc}} $$

$$ \begin{aligned} A(x) &= -A(x - 2L + cT) & x \in [0, 4L - ct_{\mathrm{dc}} - cT] \\ A(x) &= -B(x - 4L + ct_{\mathrm{dc}} + cT) & x \in [4L - ct_{\mathrm{dc}} - cT, 4L - cT] \\ A(x) &= A(x - 4L + cT) & x \in [4L - cT, 2L - ct_{\mathrm{dc}}] \\ B(x) &= A(x + cT - 2L - ct_{\mathrm{dc}}) & x \in [0, ct_{\mathrm{dc}}] \end{aligned} \tag{1.36} $$

$$ \bullet \quad 4L - ct_{\mathrm{dc}} \le cT \le 4L $$

$$ \begin{aligned} A(x) &= -B(x - 4L + ct_{\mathrm{dc}} + cT) & x \in [0, 4L - cT] \\ A(x) &= A(x - 4L + cT) & x \in [4L - cT, 2L - ct_{\mathrm{dc}}] \\ B(x) &= -A(x - 2L + cT - ct_{\mathrm{dc}}) & x \in [0, 4L - cT] \\ B(x) &= -B(x + cT - 4L) & x \in [4L - cT, ct_{\mathrm{dc}}] \end{aligned} \tag{1.37} $$
**Case** $ct_{dc} \ge 2L - ct_{dc}$

Using Equation (1.24), we can write a sequence matching for the $\phi_0$ function which is

• $ct_{dc} \le cT < 2L$

$$
\begin{align}
A(x) &= B(x - 2L + cT + ct_{dc}) & x \in [0, 2L - cT] \nonumber \\
A(x) &= -A(x - 2L + cT) & x \in [2L - cT, 2L - ct_{dc}] \nonumber \\
B(x) &= A(x + cT - ct_{dc}) & x \in [0, 2L - cT] \nonumber \\
B(x) &= B(x - 2L + cT) & x \in [2L - cT, ct_{dc}] \tag{1.38}
\end{align}
$$

• $2L < cT \le 4L - ct_{dc}$

$$
\begin{align}
A(x) &= -A(x + cT - 2L) & x \in [0, 4L - cT - ct_{dc}] \nonumber \\
A(x) &= -B(x - 4L + ct_{dc} + cT) & x \in [4L - ct_{dc} - cT, 2L - ct_{dc}] \nonumber \\
B(x) &= B(x + cT - 2L) & x \in [0, 2L - cT + ct_{dc}] \nonumber \\
B(x) &= -A(x - 2L + cT - ct_{dc}) & x \in [2L - cT + ct_{dc}, ct_{dc}] \tag{1.39}
\end{align}
$$

• $4L - ct_{dc} \le cT \le 4L$

$$
\begin{align}
A(x) &= -B(x + 4L - cT) & x \in [0, 4L - cT] \nonumber \\
A(x) &= A(x - 4L + cT) & x \in [4L - cT, 2L - ct_{dc}] \nonumber \\
B(x) &= -A(x - 2L + cT - ct_{dc}) & x \in [0, 4L - cT] \nonumber \\
B(x) &= -B(x - 4L + cT) & x \in [4L - cT, ct_{dc}] \tag{1.40}
\end{align}
$$

With this matching, it is possible to compute numerical solution to the problem at hand. Indeed, by dividing the interval $[-L, L]$, the problem is reduced to finding a non trivial kernel of a square matrix which represents each system of equations described above. The reader has to keep in mind that the unknowns are the restrictions of $A$ and $B$ over all subdivisions of $[-L, L]$.

1.2.6 Admissibility

Through Equation (1.31c), $\phi_0$ can be expressed as

$$
\forall x \in \mathbb{R}, \quad \phi_0(x) = \begin{cases} A(x) & x \in [-L, L - ct_{\text{dc}}] \\ -1 & x \in [L - ct_{\text{dc}}, L] \end{cases} \tag{1.41}
$$

Moreover, from Equation (1.31b) we can exclude every case where $A$ is not positive in the vicinity of 0. Indeed, if $\exists \eta > 0$, $A_{[0,\eta]} = -1$ then the gap becomes negative at $t = 0^+$ which is absurd. We remind the reader that $\phi(-L + \eta) = A(\eta)$. These conditions reduce to the previous cases as listed below.

**Case** $ct_{dc} \le 2L - ct_{dc}$

Using Equation (1.24), we can write a sequence matching for the $\phi_0$ function which is

• $ct_{dc} \le cT \le 2L - ct_{dc}$

$$
A(x) = \begin{cases} -1 & x \in [2L - ct_{\text{dc}} - cT, 2L - cT] \cup [cT - ct_{\text{dc}}, cT] \\ A(x + cT) & x \in [0, 2L - ct_{\text{dc}} - cT] \\ -A(x - 2L + cT) & x \in [2L - cT, 2L - ct_{\text{dc}}] \end{cases} \tag{1.42}
$$
• $2L < cT \le 2L + ct_{dc}$

$$A(x) = \begin{cases} 1 & x \in [0, cT - 2L] \cup [4L - ct_{dc} - cT, 2L - ct_{dc}] \\ -A(x - 2L + cT) & x \in [0, 4L - cT - ct_{dc}] \end{cases} \quad (1.43)$$

• $2L + ct_{dc} \le cT \le 4L - ct_{dc}$

$$A(x) = \begin{cases} 1 & x \in [4L - ct_{dc} - cT, 4L - cT] \cup [cT - 2L - ct_{dc}, cT - 2L] \\ A(x - 4L + cT) & x \in [4L - cT, 2L - ct_{dc}] \\ -A(x - 2L + cT) & x \in [0, 4L - ct_{dc} - cT] \end{cases} \quad (1.44)$$

**Case** $ct_{dc} \ge 2L - ct_{dc}$

Using Equation (1.24), we can write a sequence matching for the $\phi_0$ function which is

• $2L < cT \le 4L - ct_{dc}$

$$A(x) = \begin{cases} 1 & x \in [0, cT - 2L] \cup [4L - cT - ct_{dc}, 2L - ct_{dc}] \\ -A(x + cT - 2L) & x \in [0, 4L - cT - ct_{dc}] \end{cases} \quad (1.45)$$

## A first insight on admissible solutions

Via Section 1.2.6, we can be more precise about Equation (1.30). We will show the methodology to derive some results of Equation (1.43), the others being obtained in a similar fashion. The idea is to list all the possible cases for the position of the intervals where the function $A$ is defined. Then we show that there exists a smallest interval included in $[0, 2L - ct_{dc}]$ which must divide every interval where $A$ satisfies an equation listed on each possible case of Equations (1.43) to (1.45) and (1.52). Eventually we use this last property over a given interval to exhibit a relation between $T$ and $t_{dc}$. To be efficient, this interval must be independent of $t_{dc}$.

• We assume that $cT - 2L > 4L - ct_{dc} - cT$, then

$$A(x) = 1, \quad \forall x \in [0, 2L - ct_{dc}]. \quad (1.46)$$

However, $A$ satisfies Equation (1.52). So this case is absurd since $A$ should change sign in the vicinity of $4L - ct_{dc} - cT$.

• We further assume that $cT - 2L \le 4L - ct_{dc} - cT$, then we consider $d = cT - 2L$. $A$ is known over intervals $[0, d]$ and $[2L - ct_{dc} - d, 2L - ct_{dc}]$. If $d$ does not divide $D = 2L - ct_{dc} - 2d$ which is the length of $[d, 2L - ct_{dc} - d]$, then

$$\exists(\kappa, \delta) \in \mathbb{N} \times (0, d), \quad D = \kappa d + \delta. \quad (1.47)$$

Two cases now arise with regard to the parity of $\kappa$, we focus on the case $\kappa$ is even (the same methodology works for odd $\kappa$). If we consider $I = [D - \delta, 2L - ct_{dc} - \delta] \cap [D, 2L - ct_{dc}]$, then $\forall x \in I$, $A(x) = 1$ according to Equation (1.52), but also $\forall x \in I$, $A(x) = -1$ according to Equation (1.52), which is absurd.

• We further assume that $D = \kappa d$ with $k \in \mathbb{N}$. Two cases arise with regard to the parity of $\kappa$. If $\kappa$ is even, then Equation (1.43) yields $A(x) = 1$ and $A(x) = -1$, $\forall x \in [D, D + d]$, which is absurd. Then the only possibility is $D = (2\kappa + 1)d$. This leads to

$$T = \frac{2\kappa + 4}{2\kappa + 3} \frac{2L}{c} - \frac{t_{dc}}{2\kappa + 3}. \quad (1.48)$$

Knowing the relation between $T$, $t_{dc}$ and $d$, a function $\phi$ solution of Equation (1.31) is uniquely determined. This methodology does not give an exhaustive list of all possible solutions.
### 1.2.7 Results and discussion

In this subsection are displayed the final results derived from the methodology and equations introduced in Section 1.2.6.

#### Analytical expression of admissible $\phi$

The analytical expression of $\phi$ leading to a periodic motion can be derived as well. We only give expressions over interval $[-L, L]$ because $\phi$ is $2L$-antiperiodic.

**Duration** $ct_{dc} \le cT \le 2L - ct_{dc}$ The admissible $\phi$ function reads

$$ \phi_{1,\kappa} = \alpha_0 \times \begin{cases} (-1)^n & x \in [-L + nd, -L + (n+1)d] \\ -1 & x \in [L - ct_{dc} - cT, -L + cT] \\ (-1)^n & x \in [-L + cT + nd, -L + cT + (n+1)d] \\ -1 & x \in [L - ct_{dc}, L] \end{cases} \quad (1.49) $$

where $d = 2cT - 2L, \kappa \in \mathbb{N}, n \in [0, 2\kappa]$, $t_{dc} \in [0, L/c]$ and

$$ \alpha_0 = \frac{g_0}{8(\kappa + 1)cT - (8\kappa + 10)L}, \quad T = \frac{2\kappa + 2}{4\kappa + 3} \frac{2L}{c} - \frac{t_{dc}}{4\kappa + 3}. $$

**Duration** $2L < cT \le 2L + ct_{dc}$ or $2L < cT \le 4L - ct_{dc}$ The admissible $\phi$ reads

$$ \phi_{2,\kappa} = \alpha_0 \times \begin{cases} 1 & x \in [-L, cT - 3L] \\ (-1)^n & x \in [cT - 3L + nd, cT - 3L + (n+1)d] \\ 1 & x \in [3L - ct_{dc} - cT, L - ct_{dc}] \\ -1 & x \in [L - ct_{dc}, L] \end{cases} \quad (1.50) $$

where $d = cT - 2L, \kappa \in \mathbb{N}, n \in [0, 2\kappa]$, $t_{dc} \in [0, 2L/c]$ and

$$ \alpha_0 = \frac{g_0}{cT - 2L - ct_{dc}}, \quad T = \frac{2\kappa + 4}{2\kappa + 3} \frac{2L}{c} - \frac{t_{dc}}{2\kappa + 3}. $$

**Duration** $2L + ct_{dc} \le cT \le 4L - ct_{dc}$ The admissible $\phi$ reads

$$ \phi_{3,\kappa} = \alpha_0 \times \begin{cases} (-1)^n & x \in [-L + nd, -L + (n+1)d] \\ 1 & x \in [3L - ct_{dc} - cT, -3L + cT] \\ (-1)^{n+1} & x \in [-3L + cT + nd, -3L + cT + (n+1)d] \\ -1 & x \in [L - ct_{dc}, L] \end{cases} \quad (1.51) $$

where $d = 6L - 2cT, \kappa \in \mathbb{N}, n \in [0, 2\kappa - 1]$, $t_{dc} \in [0, L/c]$ and *

$$ \alpha_0 = \frac{g_0}{2cT - 6L}, \quad T = \frac{6\kappa + 2}{4\kappa + 1} \frac{2L}{c} - \frac{t_{dc}}{4\kappa + 1}. $$

For whichever $\kappa \in \mathbb{N}$ considered, a continuum of solutions parametrized by $t_{dc}$ is depicted. Via Equations (1.49) to (1.51), we can also parametrize every continuum by $T$. Moreover, it is possible to show that all these solutions satisfy Equation (1.31). However the solutions found until now are not an exhaustive set of solutions, more sophisticated cases can arise from Section 1.3.1.
<table><caption>Table 1.1: Overview of continua for energies of the system of interest</caption><thead><tr><th></th><th>ω/ω₁ ∈</th><th>√E(ω) (positive)</th><th>Vertical asymptote</th></tr></thead><tbody><tr><td>First NSM and sub harmonics</td><td>[2κ + 1, 2]</td><td>2g₀ / (8κ + 8)/ω - (2κ + 5)2L/c</td><td>ω/ω₁ = 4κ + 8 / (2κ + 5)</td></tr><tr><td>First NSM Internal resonances</td><td>[16κ + 4, 16κ + 4] / [12κ + 4, 12κ + 3]</td><td>2g₀ / (8/ω - 6L/c)</td><td>ω/ω₁ = 16κ + 4 / (12κ + 3)</td></tr><tr><td>Second NSM and sub harmonics</td><td>[4κ + 3, 4]</td><td>2g₀ / (32(κ + 1)/ω - (8κ + 10)L/c)</td><td>ω/ω₁ = 32(κ + 1) / (8κ + 10)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 1.3: Energy frequency diagram: (—) first NSM [$g_0 > 0$], (—) first NSM [$g_0 < 0$], (---) first NSM internal resonance [$g_0 > 0$], (—) second NSM [$g_0 > 0$], (—) second NSM [$g_0 < 0$]</figcaption></figure>

## Energy-frequency diagrams for solutions with one closing contact per period

Through the previous developments, the normalized energy of each motions can be derived as a function of the normalized frequency. We normalize both quantities with regard to the first clamped-free linear grazing mode. The is summarized in Table 1.1. Admissible periodic solutions are displayed in an energy-frequency diagram in Figure 1.3 and period versus time of contact diagram as shown in Figure 1.4. The current framework provides access to continua of solutions for some period intervals only, see Equation (1.31). We can also exhibit coexisting continuum of solutions for $ω ∈ [5ω₁/4, 4ω₁/3]$. These solutions coexist with the same energy and same frequency. Moreover, no continuum can be spotted within the frequency range [2ω₁, 3ω₁]. However, numerical suggest that other continua exist. They could be found by solving the simplified problem in Section 1.3.1. The main backbone curve of the first NSM is already known [5, 1]. Interestingly, we can satisfy admissibility for $φ$ in the case $g_0 = 0$ if we take the convention that $α_0 = +∞$.
<figure><img src="image_5.png" /><figcaption>Figure 1.4: T versus $t_{dc}$: (—) first NSM, (—) first NSM internal resonances, (—) second NSM</figcaption></figure>

## 1.3 Possible developments

### 1.3.1 Towards exhaustive continuum of solutions

It appears that other continua can be exhibited by solving Equation (1.52). It is possible to write this problem in the following simplified format: Find $A(x)$ such that

$$A(x) = -1 \qquad \forall x \in [a, b] \cup [w, d] \qquad (1.52a)$$

$$A(x) = A(x + d) \qquad \forall x \in [0, a] \qquad (1.52b)$$

$$A(x) = -A(x - b) \qquad \forall x \in [b, b + w] \qquad (1.52c)$$

$$A(x) = \pm 1 \qquad \forall x \in [0, b + w] \qquad (1.52d)$$

$$a + d = b + w \qquad (1.52e)$$

$$qd = p(b + d) - (b - a) \qquad (q, p) \in \mathbb{N}^* \times \mathbb{N}^* \qquad (1.52f)$$

where $(a, b, w, d) = (2L - cT - ct_{dc}, 2L - cT, cT - ct_{dc}, cT)$. The reader may enjoy an illustration of this problem in Figure 1.5. A similar formulation arises for Equation (1.44) with

<figure><img src="image_6.png" /><figcaption>Figure 1.5: Functional equation visualization</figcaption></figure>

different values for $(a, b, w, d)$ and a slight change in Equation (1.52).

### 1.3.2 Framework for n boundary switches per period

We can extend this formalism to the case where the bar sticks $n$ times per period. We denote by $t_{dc,i}$ (respectively $t_{f,i}$), the duration of the $i^{th}$ contact phase (respectively the duration of the $i^{th}$
free-flight phase). Then, Equation (1.11) can be extended to the n boundary-switch case

$$ \forall (x, n) \in \mathbb{R} \times \mathbb{N}^*, \quad \phi_0(x) = \beta_n \phi_0 \left( x + \sum_{k=1}^{n-1} (ct_{\mathrm{dc},n-k} + ct_{\mathrm{f},n-k}) \right) \quad (1.53) $$

with

$$ \beta_n = \prod_{i=0}^{n-1} \epsilon \left( x + \sum_{k=1}^{i} (ct_{\mathrm{dc},n-k} + ct_{\mathrm{f},n-k}) \right) \epsilon \left( x + \sum_{k=1}^{i} (ct_{\mathrm{dc},n-k} + ct_{\mathrm{f},n-k}) + t_{\mathrm{dc},n-i-1} \right). $$

We take the convention that if $i = 0$ then $\sum_{k=1}^{i} = 0$. This formulation shows how difficult the problem can be with many boundary switches. However, we still can see that looking for continuum in $\tilde{T} = \sum_{k=1}^{n-1}(ct_{\mathrm{dc},n-k} + ct_{\mathrm{f},n-k})$ leads to $\phi = \pm \tilde{\alpha}_0$. Then the displacement stays piecewise linear and stress and velocity stay piecewise constant. The argument used is the same as earlier Section 1.2.2. This formalism for two sticking phases per period, ie $n = 2$, yields

$$ \forall x \in \mathbb{R}, \quad \phi_0(x) = \beta_2 \phi_0(x + t_{\mathrm{dc},1} + t_{\mathrm{f},1} + t_{\mathrm{dc},0} + t_{\mathrm{f},0}) \quad (1.54) $$

with $\beta_2 = \epsilon(x)\epsilon(x + t_{\mathrm{dc},1})\epsilon(x + t_{\mathrm{dc},1} + t_{\mathrm{f},1})\epsilon(x + t_{\mathrm{dc},1} + t_{\mathrm{f},1} + t_{\mathrm{dc},0})$.

### 1.3.3 Checking periodicity for inferred solutions

#### General case

A candidate solution of the problem for one boundary switching per period is assumed to be known, together with the stress profile $\sigma(x) = h(x)$ at $t = t_f/2$ and $v(x) = 0$. Through Equations (1.8) and (1.9), we can compute the 2L-antiperiodic $\phi_{t_f/2}$ function associated to half of the free-phase duration:

$$ 2E\phi_{t_f/2}(x) = \begin{cases} h(x) & x \in [0, L], \\ h(-x) & x \in [-L, 0]. \end{cases} \quad (1.55) $$

Then we can compute a function $\phi_0$ associated to the motion at the release time which is

$$ \forall x \in \mathbb{R}, \quad \phi_0(x) = \phi_{t_f/2}(x - ct_f/2). \quad (1.56) $$

Equation (1.56) implies that if $D(\phi_{t_f/2})$ is known, then $D(\phi_0)$ is uniquely determined as

$$ D(\phi_0) = D(\phi_{t_f/2}) + \{ct_f/2\}. \quad (1.57) $$

Roughly speaking, if $x^*$ is a discontinuity of $\phi_{t_f/2}$, then $x^* + ct_f/2$ is a discontinuity of $\phi_0$. This discontinuity shifting is illustrated in Figure 1.6 for an internal resonance NSM located at $\omega_5/7$. We can now check the periodicity of this candidate solution via Equation (2.16) and the admissibility conditions from Equation (1.31) can be checked with $\phi_0$.

#### Illustration with the First NSM

For the first NSM, the function $h$ is defined as $h(x) = -\sigma_0, \forall x \in [0, L]$ with $\sigma \in ]g_0E/L; 0[$. Accordingly, the 2L-antiperiodic $\phi_{t_f/2}(x)$ is

$$ \phi_{t_f/2}(x) = -\sigma_0, \quad x \in [-L, L]. \quad (1.58) $$
<figure><img src="image_7.png" /><figcaption>Figure 1.6: Discontinuity shift of an internal resonance NSM located at $\omega_{5/7}$: (—) $\phi_{t_f/2}$ and (—) $\phi_0$</figcaption></figure>

Previous works have shown that $ct_f/2 = L - g_0E/\sigma_0$. Using the methodology in Section 1.3.3, we have to consider every right shift $ct_f/2 \in [0, 2L]$ of the function $\phi_{t_f/2}(x)$. Using previous results, $ct_f/2 = 2L - ct_{dc}$ and we can write the 2L-antiperiodic function $\phi_0$ as

$$ \phi_0(x) = \begin{cases} -\sigma_0 & x \in [L - ct_{dc}, L], \\ +\sigma_0 & x \in [-L, L - ct_{dc}]. \end{cases} \qquad (1.59) $$

Then admissibility holds because

$$ \int_{-L}^{ct-L} \phi_0(s) \, ds = \begin{cases} \sigma_0 ct & t \in [0, 2L/c - t_{dc}] \\ \sigma_0(2L - ct_{dc}) - \sigma_0(ct - 2L + ct_{dc}) & t \in [2L/c - t_{dc}, 4L/c - 2t_{dc}]. \end{cases} $$

When $t = t_f = 4L/c - 2t_{dc}$, the closure happens because $\int_{-L}^{ct_f-L} \phi_0(s) \, ds = 0$. Moreover, the compression condition (1.31c) is satisfied as seen in Equation (1.59). Considering that $T = 4L/c - t_{dc}$, we can now compute $\phi_0(x + cT)$ using the 2L-antiperiodicity of $\phi_0$:

$$ \phi_0(x + cT) = +\sigma_0, \quad x \in [-L, L] \qquad (1.60) $$

and $\phi_0(\bullet + cT)$ is 2L-antiperiodic. Then the periodicity condition is satisfied because

$$ \epsilon(x)\epsilon(ct_{dc} + x)\phi_0(cT + x) = \phi_0(x) = \begin{cases} -\sigma_0 & x \in [L - ct_{dc}, L] \\ +\sigma_0 & x \in [-L, L - ct_{dc}]. \end{cases} \qquad (1.61) $$

### 1.3.4 Extension to Robin boundary condition at $x = 0$

The system of interest, illustrated in Figure 1.7, is almost identical to the previous one. The bar is now attached to the ground at $x = 0$ through a spring of stiffness $k > 0$. The governing Equation (1.1) still holds as well as the general solution (1.5). However, the Robin boundary condition $u_x(0, t) = ku(0, t)$ implies

$$ \forall x \in \mathbb{R}, \quad f'(x) - h'(x) = k(f(x) + h(x)). \qquad (1.62) $$
<figure><img src="image_8.png" /><figcaption>Figure 1.7: One-dimensional bar with Robin and Signorini boundary conditions</figcaption></figure>

Accordingly, the function h can be written as

$$ 
\begin{aligned}
\forall x \in \mathbb{R}, \quad h(x) &= e^{-kx} \left[ h(0) + \int_{0}^{x} e^{ks} (f'(s) - kf(s)) \mathrm{d}s \right] \\
&= e^{-kx} \left[ h(0) - 2k \int_{0}^{x} e^{ks} f(s) \mathrm{d}s + e^{kx} f(x) - f(0) \right] \\
&= f(x) + e^{-kx} \left[ h(0) - f(0) - 2k \int_{0}^{x} e^{ks} f(s) \mathrm{d}s \right]
\end{aligned} 
\quad (1.63) $$

Then $\forall (x, t) \in \mathbb{R} \times \mathbb{R}$

$$ u(x,t) = f(ct + x) + f(ct - x) + e^{k(x-ct)} \left[ h(0) - f(0) - 2k \int_{0}^{ct-x} e^{ks} f(s) \mathrm{d}s \right] \quad (1.64) $$

**Free flight** During free flights we have that $\forall t \in \mathbb{R}$, $u_x(L,t) = 0$, then $\forall x \in \mathbb{R}$

$$ f'(x + 2L) = h'(x) \Leftrightarrow f(x + 2L) = h(x) + f(2L) - h(0) \quad (1.65) $$

Using Equation (1.63), the problem reduces to finding f such that $\forall x \in \mathbb{R}$

$$ f(x + 2L) - f(2L) + h(0) - f(x) = e^{-kx} \left[ h(0) - f(0) - 2k \int_{0}^{x} e^{ks} f(s) \mathrm{d}s \right] \quad (1.66) $$

**Contact phase** During contact phases we have that $\forall t \in \mathbb{R}$, $u(L,t) = g_0$, then $\forall x \in \mathbb{R}$: $f(x + 2L) + h(x) = g_0$. Using Equation (1.63), the problem is reduced to finding f such that $\forall x \in \mathbb{R}$

$$ g_0 - f(x + 2L) - f(x) = e^{-kx} \left[ h(0) - f(0) - 2k \int_{0}^{x} e^{ks} f(s) \mathrm{d}s \right] \quad (1.67) $$

which implies $h(0) = g_0 - f(2L)$ during contact phases.

These equations do not give valuable insight about f. We see here, that analytical solutions of this problem may be way more difficult to find than in the case where the left extremity is clamped to the wall. However we can still observe that f can be described over $\mathbb{R}$ if we know it over an interval of length $2L$.
# Chapter 2

## On a stability analysis

We are now interested in the stability of the periodic solutions reported in the previous chapter. Our main purpose is to check if a naive extension of the first return map used in [9] is possible in this framework.

### 2.1 A stress-velocity formulation

D'Alembert solution expressed in terms of initial displacement $u_0(x)$ and velocity $v_0(x)$ reads

$$2u(x,t) = u_0(x+ct) + u_0(x-ct) + \frac{1}{c} \int_{x-ct}^{x+ct} v_0(s) \, ds. \qquad (2.1)$$

For sufficiently differentiable $u_0$ and in the framework of linear elasticity, the following holds via the chain rule for differentiation

$$2 \begin{pmatrix} \sigma(x,t) \\ v(x,t) \end{pmatrix} = \begin{pmatrix} \sigma_0(x+ct) + \sigma_0(x-ct) + E/c[v_0(x+ct) - v_0(x-ct)] \\ c/E[\sigma_0(x+ct) - \sigma_0(x-ct)] + v_0(x+ct) + v_0(x-ct) \end{pmatrix}. \qquad (2.2)$$

Hence, for a given initial condition in stress and velocity, we can compute the stress and the velocity for any $x \in [0, L]$ and any $t \ge 0$. The system of interest and the governing equations are the same as in Section 1.1.1.

#### 2.1.1 Hyperbolic system of conservation law

The equilibrium equation of the bar Equation (1.1) can be written as a system of two first order partial differential equations in terms of $\sigma(x, t)$ and $v(x, t)$:

$$\partial_t \mathbf{q} + \mathbf{B} \partial_x \mathbf{q} = \mathbf{0} \quad \text{where} \quad \mathbf{B} = \begin{bmatrix} 0 & -E \\ -c^2/E & 0 \end{bmatrix} \quad \text{and} \quad \mathbf{q} = \begin{pmatrix} \sigma \\ v \end{pmatrix}. \qquad (2.3)$$

The eigenvalues of $\mathbf{B}$ are $c$ and $-c$ which correspond to the wave propagation velocity along the bar. That refers to a hyperbolic system of conservation laws [2]. This local equation of conservation law has the integral form

$$\partial_t \left( \int_{x_1}^{x_2} \mathbf{q}(x,t) \, dx \right) = \mathbf{B} (\mathbf{q}(x_1, t) - \mathbf{q}(x_2, t)). \qquad (2.4)$$
### 2.1.2 Switching mechanism

Through the expressions of $\sigma$ and $v$ in terms of $\phi$, we can write the periodic extensions for the stress and velocity profiles on the whole real axis.

**Free flights** $\sigma$ and $v$ are $2L$-antiperiodic: $\forall x \in \mathbb{R}, \forall t \in \mathbb{R}$

$$ 
\begin{aligned}
\sigma(L + x, t) &= -\sigma(L - x, t) \\
\sigma(x, t) &= \sigma(-x, t) \\
v(L + x, t) &= v(L - x, t) \\
v(x, t) &= -v(-x, t)
\end{aligned} 
\qquad (2.5) $$

**Contact phases** $\sigma$ and $v$ are $2L$-periodic: $\forall x \in \mathbb{R}, \forall t \in \mathbb{R}$

$$ 
\begin{aligned}
\sigma(L + x, t) &= \sigma(L - x, t) \\
\sigma(x, t) &= \sigma(-x, t) \\
v(L + x, t) &= -v(L - x, t) \\
v(x, t) &= -v(-x, t)
\end{aligned} 
\qquad (2.6) $$

### 2.1.3 Mapping for a given initial condition

**Free flight operator**

We can now write d'Alembert's solution at the beginning of the free-flight as

$$ 2 \begin{pmatrix} \sigma(x,t) \\ v(x,t) \end{pmatrix} = \begin{pmatrix} \sigma_0^{**}(x+ct) + \sigma_0^{**}(x-ct) + E/c[v_0^*(x+ct) - v_0^*(x-ct)] \\ c/E[\sigma_0^{**}(x+ct) - \sigma_0^{**}(x-ct)] + v_0^*(x+ct) + v_0^*(x-ct) \end{pmatrix} \qquad (2.7) $$

and in a compact manner

$$ \mathbf{q}(x, t) = \mathcal{F}(\mathbf{q}_0)(x, t). \qquad (2.8) $$

Here $\bullet$** and $\bullet$* are the periodic extensions of stress and velocity for free flights.

**Contact phase operator**

We can also write the d'Alembert solution regarding to initial conditions $\sigma_1$ and $v_1$ at the beginning of the contact phase as following:

$$ 2 \begin{pmatrix} \sigma(x,t) \\ v(x,t) \end{pmatrix} = \begin{pmatrix} \hat{\sigma}_1(x + ct) + \hat{\sigma}_1(x - ct) + E/c[\tilde{v}_1(x + ct) - \tilde{v}_1(x - ct)] \\ c/E[\hat{\sigma}_1(x + ct) - \hat{\sigma}_1(x - ct)] + \tilde{v}_1(x + ct) + \tilde{v}_1(x - ct) \end{pmatrix} \qquad (2.9) $$

or

$$ \mathbf{q}(x, t) = \mathcal{C}(\mathbf{q}_1)(x, t). \qquad (2.10) $$

Here $\bullet$ and $\tilde{\bullet}$ are the periodic extensions of stress and velocity for contact phases. Both operators $\mathcal{F}$ and $\mathcal{C}$ have to be seen as operators which convert an initial wave into its image at time $t$.
## Standalone properties

A few important properties of $\mathcal{F}$ and $\mathcal{C}$ are listed below but basic proofs are provided for $\mathcal{F}$ only for conciseness. This relies on Equation (2.2) along with the linearity of function evaluation.

**Linearity with respect to $q_0$** The linearity of $\mathcal{F}$ stems from the linearity of the governing PDE (without the boundary conditions).

**Characteristics curves relations** $\mathcal{F}$ satisfies this property:

$$2\mathcal{F}(\mathbf{q})(x, t + \Delta t) = \mathcal{F}(\mathbf{q})(x + c\Delta t, t) + \mathcal{F}(\mathbf{q})(x - c\Delta t, t) - \mathbf{B}/c[\mathcal{F}(\mathbf{q})(x + c\Delta t, t) - \mathcal{F}(\mathbf{q})(x - c\Delta t, t)] \quad (2.11)$$

for an arbitrary $\mathbf{q}$. This expression holds for every $\Delta t$ considered and $\forall x \in \mathbb{R}$. We can illustrate this equality in terms of characteristics in Figure 2.1. It is similar to properties reported in [8].

<figure><img src="image_9.png" /><figcaption>Figure 2.1: Characteristics curves around $\mathbf{q}(x, t)$</figcaption></figure>

**Spatial periodicity** The periodic extensions of free-flight or contact phase yield $\mathcal{F}(\mathbf{q})(x, t + 4L/c) = \mathcal{F}(\mathbf{q})(x, t)$. Similarly, $\mathcal{C}(\mathbf{q})(x, t + 2L/c) = \mathcal{C}(\mathbf{q})(x, t)$. Moreover, from Equation (2.11) and extensions of $\sigma$ and $v$ with $\Delta t = 2L/c$, $\mathcal{F}(\mathbf{q})(x, t + 2L/c) = -\mathcal{F}(\mathbf{q})(x, t)$.

**Temporal differentiability of $\mathcal{F}$** As a consequence of Equation (2.3), we can compute the derivative in the sense of distributions of $\mathcal{F}$ through

$$\forall \varphi \in C_{\text{comp}}^{\infty}, \quad \int_{\mathbb{R}} \partial_t \mathcal{F}(\mathbf{q}_0)(x, t) \varphi(x) \, dx = \int_{\mathbb{R}} -\mathbf{B} \partial_x \mathcal{F}(\mathbf{q}_0)(x, t) \varphi(x) \, dx. \quad (2.12)$$

In the framework of almost everywhere (a.e.) differentiable functions, the following holds via d'Alembert solution:

$$\forall \varphi \in C_{\text{comp}}^{\infty}, \quad \int_{\mathbb{R}} \partial_t \mathcal{F}(\mathbf{q}_0)(x, t) \varphi(x) \, dx = \int_{\mathbb{R}} -\mathbf{B} \mathcal{F}(\partial_x \mathbf{q}_0)(x, t) \varphi(x) \, dx. \quad (2.13)$$

The previous means that for almost every $x \in \mathbb{R}$ and for almost every $t \ge 0$

$$\partial_t \mathcal{F}(\mathbf{q}_0)(x, t) = -\mathbf{B} \mathcal{F}(\partial_x \mathbf{q}_0)(x, t). \quad (2.14)$$

This also works for operator $\mathcal{C}$.
Mapping

We now use the fact that a switch happens at $t = t_f$ makes the link between the two steps of the motion to write:

$$ q_2(x) = C[\mathcal{F}(q_0)(\bullet, t_f)](x, t_c). \qquad (2.15) $$

When computed with an initial condition generating a periodic motion of period $T = t_f + t_{dc}$, we have

$$ \forall x \in [0, L], \quad C[\mathcal{F}(q_0)(\bullet, t_f)](x, t_c) = q_0(x). \qquad (2.16) $$

## 2.2 Stability under small perturbations

### 2.2.1 Definition

We call a stable solution, every $q$ such that $\forall \epsilon \ge 0, \forall x \in [0, L]$

$$ \exists \delta \ge 0, \quad \|q_0(x) - \delta\psi_0(x)\| \le \delta \quad \Rightarrow \quad \forall t \ge 0, \quad \|q(x, t) - \delta\psi(x, t)\| \le \epsilon \qquad (2.17) $$

where $q_0$ and $\delta\psi_0$ are respectively the initial condition of a periodic solution and a disturbed initial condition. Roughly speaking, this definition means that if a motion is sufficiently close to a periodic motion at $t = 0$, then the motion will stay close to the periodic motion [6]. If the solution $q$ does not match this, it will be called an unstable solution. An illustration for this definition of stability may be seen on Figure 2.2.

<figure><img src="image_10.png" /><figcaption>Figure 2.2: Periodic motion (—), Disturbed motion 1 (---), Disturbed motion 2 (—)</figcaption></figure>

### 2.2.2 Assumptions

We now want to exhibit an operator which maps an initial small perturbation $\delta\psi_0$ to the perturbation at the time of first return $T$ denoted as $\delta\psi_T$. We will denote this operator as $\mathcal{L}_q$. The index $q$ reminds us that this operator is a function of the periodic solution $q$ considered. We draw our focus on particular initial perturbations $\delta\psi_0$ which satisfy the following conditions:

$$ \forall x \in [0, L], \quad \|\delta\psi_0(x)\| \ll \|q_0(x)\| \qquad (2.18a) $$

$$ q_0 + \delta\psi_0 \text{ satisfies the fixed-free boundary condition} \qquad (2.18b) $$

$$ \delta\psi_0 \text{ is Lipschitz-continous along } \mathbb{R} \qquad (2.18c) $$

$$ \forall t \ge 0, \quad \sigma(L, t_r^+) < -\delta\sigma(L, t) = -E\partial_x\delta u(L, t) < \sigma(L, t_r^-) + J_q \qquad (2.18d) $$
<figure><img src="image_11.png" /><figcaption>Figure 2.3: Assumption (2.18d): Undisturbed stress (—), Disturbed stress (—)</figcaption></figure>

Equation (2.18d) is illustrated in Figure 2.3. The reader can see that the initial stress wave front stays large enough to ensure that the release occurs because of the reflection of the compression wave at the contacting end. This assumption is crucial to compute the variation of closure time.

### 2.2.3 Linearized mapping for an initial perturbation

#### Variation of the duration of free-flight

Without any loss of generality, we assume that the motion starts at a release time with an initial condition $q_0 + \delta\psi_0 = (\sigma_0 + \delta\sigma_0, v_0 + \delta v_0)^\top$. The duration of free-flight will be changed by the perturbation. Here, we use $u(x, t)$ and $\delta u(x, t)$ to refer to the displacement of the bar and the perturbation in displacement of the bar. By expansion around $t_f$ at $x = L$ and using that closure happens at $t_f$, we get:

$$ \delta t_f = - \frac{\delta u(L, t_f)}{v(L, t_f)} + o(\delta) \qquad (2.19) $$

Linear elasticity yields

$$ \delta t_f \approx - \int_0^{t_f} \frac{\mathbf{e}_2^\top \delta\psi(L, \tau)}{\mathbf{e}_2^\top \mathbf{q}(L, t_f)} d\tau \qquad (2.20) $$

or equivalently

$$ \delta t_f \approx - \int_0^L \frac{1}{E} \frac{\mathbf{e}_1^\top \delta\psi(s, t_f)}{\mathbf{e}_2^\top \mathbf{q}(L, t_f)} ds \qquad (2.21) $$

Here the assumptions given by Equations (2.18a) to (2.18c) allow for this current development. The corresponding output is shown in Figure 2.4. The quantity $\delta t_f$ has to be understood as the signed variation of free-flight duration due to the perturbation.

#### Variation of the duration of contact-phase

We now take advantage of assumptions Equations (2.18c) and (2.18d). Here it means that the magnitude of the stress perturbation is small enough to ensure that $(\sigma + \delta\sigma)(L, \cdot)$ keeps the same number of discontinuities and the jump still cross the real axis. Physically, this assumption
<figure><img src="image_12.png" /><figcaption>Figure 2.4: Free-flight duration variation: Numericals (—), δf Linearization (—)</figcaption></figure>

is a translation of Equation (2.18a) in terms of stress. By spotting the discontinuities of the stress during the motion, we know the position α ∈ [0, L] at the closure time of the discontinuity in stress which makes the release occur. Thanks to the wave traveling velocity c, we are able to know the exact contact-phase duration tdc = (L − α)/c. Since we have a perturbation, the α is translated to α + δα and the release is still a consequence of the jump in stress thanks to Equation (2.18d). Hence:

$$ct_{dc} + c\delta t_{dc} = L - (\alpha + \delta\alpha) \quad \Rightarrow \quad \delta t_{dc} = -\delta\alpha/c \qquad (2.22)$$

By using the free-flight variation duration δt<sub>f</sub> and wave velocity c, we can express the variation of the position δα, that leads to:

$$c\delta t_{dc} = -\delta\alpha = -c\delta t_{f} \qquad (2.23)$$

The reader can see the relationship (obtained by numerical) between both duration variations in Figure 2.5. The δt<sub>dc</sub> has to be understood as the signed contact-phase duration variation due

<figure><img src="image_13.png" /><figcaption>Figure 2.5: Contact-phase duration versus free-flight duration: Numericals</figcaption></figure>

to the perturbation.
### 2.2.4 Linearized “First Return Map” operator

We are now able to describe the operator which makes the link between an initial small perturbation $\delta\psi_0$ and the due perturbation after one complete free-flight and one complete contact-phase $\delta\psi_T$. In order to do so, we are going to use the previous operators $\mathcal{F}_{\bullet,t_f}$ and $\mathcal{C}_{\bullet,t_{dc}}$ and the duration variation of both phases of the motion $\delta t_f$ and $\delta t_{dc}$. By expansion of the operators evaluated in $\mathbf{q}_0 + \delta\psi_0$ and in $t = T$ we get the following:

$$ \delta\psi(x, T) \approx C[\mathcal{F}(\delta\psi_0)(\bullet, t_f)](x, t_{dc}) + \gamma \int_0^{t_f} \frac{e_2^\top \mathcal{F}(\delta\psi_0)(L, \tau)}{e_2^\top \mathcal{F}(\mathbf{q}_0)(L, t_f)} d\tau \quad (2.24) $$

where $\gamma = C[\mathcal{B}\mathcal{F}(\partial_x\mathbf{q}_0)(\bullet, t_f)](x, t_{dc}) - C[\partial_x\mathcal{F}(\mathbf{q}_0)(\bullet, t_f)](x, t_{dc})$. Referring to Equation (2.21), it can also be written as

$$ \delta\psi(x, T) \approx C[\mathcal{F}(\delta\psi_0)(\bullet, t_f)](x, t_{dc}) + \frac{\gamma}{E} \int_0^L \frac{e_1^\top \mathcal{F}(\delta\psi_0)(s, t_f)}{e_2^\top \mathcal{F}(\mathbf{q}_0)(L, t_f)} ds, \quad (2.25) $$

or in the compact form

$$ \text{For (a.e) } x \in [0, L], \quad \delta\psi(x, T) \approx \mathcal{L}(\delta\psi_0)(\mathbf{q}_0, x, t_{dc}, t_f). \quad (2.26) $$

## 2.3 Numerical routine for eigenvalues computation

The stability of a given periodic solution $\mathbf{q}$ relies on the analysis of the operator $\mathcal{L}_{\mathbf{q}_0,x,t_{dc},t_f}$. As it is challenging to conduct such investigation in the continuous framework, a discretized version is explored instead.

### 2.3.1 Discretized mapping for an initial perturbation

Using a space step $\Delta x = L/N$ and a time step $\Delta t = \Delta x/c$ and the WFEM matrices [12], we can give a matrix expression of a discretized version of $\mathcal{L}$. The details of the discrete formulation are given in appendix. Two choices arise when we compute this discrete version. We can time integrate velocity or space integrate stress in order to get the perturbed displacement.

The purpose of the discretization is to exhibit a linear operator whose input belongs to a finite dimension space. Indeed, the computation of the spectral radius is easier in finite dimension. Then, we can use the same reasoning than [9] to conclude on the stability of some branches.

#### Time integration of the perturbation

$$ \delta\psi_T = A_c^p A_f^m \delta\psi_0 + \frac{1}{2} \frac{\Delta t}{\Delta x} \langle e_{2N} | \mathbf{F} \delta\psi_0 \rangle \mathbf{x}_{\text{NSM}} \quad (2.27) $$

$$ = [A_c^p A_f^m + \frac{1}{2c} \mathbf{x}_{\text{NSM}} e_{2N}^\top \mathbf{F}] \delta\psi_0 = [\mathbf{R}_{\text{NSM}}^{(1)}] \delta\psi_0 \quad (2.28) $$

$R_{\text{NSM}}^{(1)}$ is the discretized linear operator of first return. We may compute its eigenvalues with numerical methods. The vector $\mathbf{x}_{\text{NSM}}$ is a NSM considered dependent vector. We see clearly here that the operator is linear in $\delta\psi_0$ thanks to the matrix expression of it.

Matrices $A_c^p$ and $A_f^m$ are the discretized equivalents of $\mathcal{C}_{\bullet,p\Delta t}$ and $\mathcal{F}_{\bullet,m\Delta t}$. We denote as $\delta\psi_0$ is the discretized shape function of the initial perturbation along the bar. The $\mathbf{F}$ matrix and $\mathbf{x}_{\text{NSM}}$ vector are described in appendix. We denote as $e_{2N}$ the $2N^{\text{th}}$ canonical basis vector of $\mathbb{R}^{2N}$.
## Space integration of the perturbation

$$ \delta\psi_T = A_c^p A_f^m \delta\psi_0 + \frac{1}{E} \sum_{k=1}^{N} \langle e_k | A_f^m \delta\psi_0 \rangle x_{\text{NSM}} \quad (2.29) $$

$$ = [A_c^p A_f^m + \frac{1}{E} x_{\text{NSM}} (\sum_{k=1}^{N} e_k)^{\top} A_f^m] \delta\psi_0 = [R_{\text{NSM}}^{(2)}] \delta\psi_0 \quad (2.30) $$

## 2.3.2 Results and discussions

We focus the results on the main backbone curve of the first NSM. The reader may see the spectral radius of both first return maps $R_{\text{NSM}}^{(1)}$ and $R_{\text{NSM}}^{(2)}$ in Figure 2.6. What we see is that motions should be unstable. When the frequency converges to linear frequency of the system, then we are close to neutral stability ($\rho = 1$). However the results given by both operators are different in terms of amplitude even if the meaning stays the same. This drawback could be linked to a huge number of numerical operations necessary to get the eigenvalues of $R_{\text{NSM}}^{(1)}$. However, when we compute naively a disturbed motion near a periodic orbit with numerical it

<figure><img src="image_14.png" /><figcaption>Figure 2.6: Spectral radius of discretized operators</figcaption></figure>

appears that the motion is stable. Indeed, we can notice that the perturbation appears to behave as an additional superposition of forward and backward traveling waves as shown in Figures 2.7 and 2.8.

For all these reasons, it appears that we can not extend the first return map tool to continuous system by this methodology. As suggested by Figure 2.9, the perturbed motions could be seen as waves superposition.
<figure><img src="image_15.png" /><figcaption>Figure 2.7: Residual perturbation of initial amplitude $\delta$: initial perturbation (—), perturbation at second release (→)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 2.8: Residual perturbation of initial amplitude $\delta$: initial perturbation (—), perturbation at second release (→)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>(a) Undisturbed motion</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) Perturbed motion with a 2.5% amplitude size initial perturbation</figcaption></figure>

Figure 2.9: Periodic motion living on the main backbone curve of 1<sup>NSM</sup> and a disturbed motion
# Appendices
# Appendix A

## Notations

Notations are given in order of appearance:

<table><tr><td>Lower case</td><td>Scalar quantities and functions</td></tr><tr><td>Lower case, bold font</td><td>Vector quantities and functions</td></tr><tr><td>Upper case, bold font</td><td>Matrix quantities and functions</td></tr><tr><td>L</td><td>Length of the bar</td></tr><tr><td>E</td><td>Young modulus of the bar</td></tr><tr><td>S</td><td>Section of the bar</td></tr><tr><td>ρ</td><td>linear density</td></tr><tr><td>g<sub>0</sub></td><td>Initial gap between the right tip and the wall</td></tr><tr><td>x</td><td>Straight curvilinear abscissa of the bar</td></tr><tr><td>u</td><td>Displacement of the section</td></tr><tr><td>∂<sub>•</sub></td><td>Partial differentiation with respect to •</td></tr><tr><td>c</td><td>Elastic waves celerity</td></tr><tr><td>u<sub>0</sub></td><td>Initial displacement field</td></tr><tr><td>v</td><td>Velocity of the section</td></tr><tr><td>v<sub>0</sub></td><td>Initial velocity field</td></tr><tr><td>g</td><td>Current gap between the right tip and the wall</td></tr><tr><td>σ</td><td>Stress applied to the section</td></tr><tr><td>f</td><td>Forward traveling wave</td></tr><tr><td>h</td><td>Backward traveling wave</td></tr><tr><td>φ</td><td>Wave function</td></tr><tr><td>R</td><td>Real numbers set</td></tr><tr><td>φ<sub>0</sub></td><td>Initial wave function</td></tr><tr><td>φ<sub>1</sub></td><td>Wave after the contact</td></tr><tr><td>φ<sub>2</sub></td><td>Wave after the release</td></tr><tr><td>ε</td><td>Function defined in Figure 1.2</td></tr><tr><td>t<sub>f</sub></td><td>Duration of the free flight</td></tr><tr><td>t<sub>dc</sub></td><td>Duration of the contact</td></tr><tr><td>T</td><td>Period of the motion</td></tr><tr><td>Q</td><td>Rational numbers set</td></tr><tr><td>D(•)</td><td>Set of the discontinuity of function •</td></tr><tr><td>N<sup>(*)</sup></td><td>Natural numbers set (zero excluded)</td></tr><tr><td>Z<sup>(*)</sup></td><td>Relative numbers set (zero excluded)</td></tr><tr><td>A(x), B(x)</td><td>Sequence matching functions</td></tr><tr><td>α<sub>0</sub></td><td>Energy constant for each mode</td></tr><tr><td>ω</td><td>Angular frequency</td></tr><tr><td>k</td><td>Stiffness of the Robin boundary condition</td></tr></table>
$\delta\psi_0$ Initial perturbation

$\delta\psi$ Residual perturbation

$x_{NSM}$ Discretized profile of the mode studied starting at $t = \frac{t_f}{2}$

$L$ Operator of first return applying to a discretized profile

$R_{NSM}^{(1)}$ Discretized operator of first return in matrix form (first expression)

$R_{NSM}^{(2)}$ Discretized operator of first return in matrix form (second expression)

$\delta\bullet$ Increment of $\bullet$ due to a small perturbation

$\Delta t, \Delta x$ Time and space segment

$e$ Canonical basis vector of $\mathbb{R}^d$

$\mathcal{F}$ Propagation operator for the free-flight phase

$\mathcal{C}$ Propagation operator for the contact phase

$\mathcal{L}$ First order approximation of the first return map
# Appendix B

## Computation of the linearized First Return Map $\mathcal{L}$

By expansion of the operators evaluated at $\mathbf{q}_0 + \delta\psi_0$ and $t = T$, we get the following $\forall x \in [0, L]$:

$$ (\mathbf{q} + \delta\psi)(x, T) = \mathcal{C}[\mathcal{F}(\mathbf{q}_0 + \delta\psi_0)(\bullet, t_f + \delta t_f)](x, t_{dc} + \delta t_{dc}) \quad (B.1) $$

$$ \begin{aligned} \mathbf{q}(x, T) + \delta\psi(x, T) &= \mathcal{C}[\mathcal{F}(\mathbf{q}_0)(\bullet, t_f + \delta t_f) + \mathcal{F}(\delta\psi_0)(\bullet, t_f + \delta t_f)](x, t_{dc} + \delta t_{dc}) \\ &= \mathcal{C}[\mathcal{F}(\mathbf{q}_0)(\bullet, t_f) + \delta t_f \partial_t \mathcal{F}(\mathbf{q}_0)(\bullet, t_f) + \mathcal{F}(\delta\psi_0)(\bullet, t_f)](x, t_{dc} + \delta t_{dc}) + \circ(\delta) \\ &= \mathcal{C}(\tilde{\mathbf{q}}_1)(x, t_{dc} + \delta t_{dc}) \\ &= \mathcal{C}(\tilde{\mathbf{q}}_1)(x, t_{dc}) + \delta t_{dc} \partial_t \mathcal{C}(\tilde{\mathbf{q}}_1)(x, t) + \circ(\delta) \end{aligned} \quad (B.2) $$

From the linearity of the operator $\mathcal{C}$, we get $\forall x \in [0, L]$:

$$ \begin{aligned} \mathbf{q}(x, T) + \delta\psi(x, T) &= \mathcal{C}(\mathcal{F}(\mathbf{q}_0)(\bullet, t_f))(x, t_{dc}) \\ &\quad + \delta t_f \mathcal{C}[\partial_t \mathcal{F}(\mathbf{q}_0)(\bullet, t_f)](x, t_{dc}) + \mathcal{C}[\mathcal{F}(\delta\psi_0)(\bullet, t_f)](x, t_{dc}) \\ &\quad + \delta t_{dc} \partial_t \mathcal{C}[\mathcal{F}(\mathbf{q}_0)(\bullet, t_f)](x, t) + \circ(\delta) \end{aligned} \quad (B.3) $$

The previous closure conditions along with Equations (2.20) and (2.23) yields

$$ \delta\psi(x, T) = \mathcal{C}[\mathcal{F}(\delta\psi_0)(\bullet, t_f)](x, t_{dc}) - \left[ \int_0^{t_f} \frac{\mathbf{e}_2^\top \delta\psi(L, \tau)}{\mathbf{e}_2^\top \mathbf{q}(L, t_f)} d\tau \right] \mathbf{a} + \circ(\delta) \quad (B.4) $$

with $\mathbf{a} = \mathcal{C}[\partial_t \mathcal{F}(\mathbf{q}_0)(\bullet, t_f)](x, t_{dc}) - \partial_t \mathcal{C}[\mathcal{F}(\mathbf{q}_0)(\bullet, t_f)](x, t_{dc})$. Now we use the fact that $\mathcal{F}$ relates $\delta\psi_0$ to $\delta\psi(L, t_f)$ (the same result holds for $\mathbf{q}(L, t_f)$). This quantity has to be understood as the modified position of the contacting end due to the perturbation at $t = t_f$, which gives:

$$ \delta\psi(x, T) = \mathcal{C}[\mathcal{F}(\delta\psi_0)(\bullet, t_f)](x, t_{dc}) - \left[ \int_0^{t_f} \frac{\mathbf{e}_2^\top \mathcal{F}(\delta\psi_0)(L, \tau)}{\mathbf{e}_2^\top \mathcal{F}(\mathbf{q}_0)(L, t_f)} d\tau \right] \mathbf{a} \quad (B.5) $$

where terms of order greater than one are discarded and $T = t_f + t_{dc}$. Using the derivatives properties of $\mathcal{F}$ and $\mathcal{C}$ shown in Equation (2.13) we have

$$ \delta\psi(x, T) = \mathcal{C}[\mathcal{F}(\delta\psi_0)(\bullet, t_f)](x, t_{dc}) - \left[ \int_0^{t_f} \frac{\mathbf{e}_2^\top \mathcal{F}(\delta\psi_0)(L, \tau)}{\mathbf{e}_2^\top \mathcal{F}(\mathbf{q}_0)(L, t_f)} d\tau \right] \mathbf{b} \quad (B.6) $$

with $\mathbf{b} = \mathcal{C}[\mathbf{B}\mathcal{F}(\partial_x\mathbf{q}_0)(\bullet, t_f)](x, t_{dc}) - \mathbf{B}\mathcal{C}[\partial_x\mathcal{F}(\mathbf{q}_0)(\bullet, t_f)](x, t_{dc})$. From Equation (2.21), it can also be written as

$$ \delta\psi(x, T) = \mathcal{C}[\mathcal{F}(\delta\psi_0)(\bullet, t_f)](x, t_{dc}) - \frac{1}{E} \left[ \int_0^{t_f} \frac{\mathbf{e}_1^\top \mathcal{F}(\delta\psi_0)(s, \tau)}{\mathbf{e}_2^\top \mathcal{F}(\mathbf{q}_0)(L, t_f)} ds \right] \mathbf{b}. \quad (B.7) $$
# Appendix C

## Discretization of the linearized First Return Map $\mathcal{L}$

### C.1 Discretized first return map

#### C.1.1 Space discretization

Considering the space discretization of the WFEM, we consider a constant space step $\Delta x$ along the bar of length $L$, ie $\exists N \in \mathbb{N}$ such that $\Delta x = L/N$ and assume $\Delta t = \Delta x/c$. An initial vector shape is seen as $(\boldsymbol{\sigma}^{\top}, \mathbf{v}^{\top})^{\top} = (\sigma_1 \dots \sigma_N v_1 \dots v_N)$.

#### C.1.2 Discretized equivalents of $\mathcal{C}$ and $\mathcal{F}$

Regarding the meaning of these operators, the discretized equivalent would be an iteration of the WFEM matrices. Indeed, if we consider $t_{dc} = p\Delta t$ and $t_f = m\Delta t$, the following holds: $\mathcal{C} \approx A_c^p$ and $\mathcal{F} \approx A_f^m$ where these matrices are the same than those described in Appendix of [12]. We are going to explain how we can find the discretized formulation with an other way than in previous works.

**Proof for $\mathcal{F} \approx A_f^m$**

It is shown in Equation (2.11) that the time and space evolutions are linked explicitly. This property allows to show that using $\mathcal{F}$ leads to the discrete formulation introduced by the WFEM for wave propagation along an elastic bar as described in [12] by considering an infinitesimal $\Delta t$ and taking $\Delta x = c\Delta t$ [3]:

$$2 \begin{pmatrix} \sigma_i^{(n+1)} \\ v_i^{(n+1)} \end{pmatrix} = \begin{pmatrix} \sigma_{i+1}^{(n)} + \sigma_{i-1}^{(n)} + E/c(v_{i+1}^{(n)} - v_{i-1}^{(n)}) \\ c/E(\sigma_{i+1}^{(n)} - \sigma_{i-1}^{(n)}) + v_{i+1}^{(n)} + v_{i-1}^{(n)} \end{pmatrix} \quad (C.1)$$

where $\bullet_i^{(n)}$ holds for the value of $\bullet$ in the $i^{th}$ cell of the bar at time step $(n)$. From the parity properties of $\sigma$ and $v$, we can also find the state of the ghost cells:

$$2 \begin{pmatrix} \sigma_1^{(n+1)} \\ v_1^{(n+1)} \end{pmatrix} = \begin{pmatrix} \sigma_2^{(n)} + \sigma_1^{(n)} + E/c(v_2^{(n)} + v_1^{(n)}) \\ c/E(\sigma_2^{(n)} - \sigma_1^{(n)}) + v_2^{(n)} - v_1^{(n)} \end{pmatrix} \quad (C.2)$$

and

$$2 \begin{pmatrix} \sigma_N^{(n+1)} \\ v_N^{(n+1)} \end{pmatrix} = \begin{pmatrix} -\sigma_N^{(n)} + \sigma_{N-1}^{(n)} + E/c(v_N^{(n)} - v_{N-1}^{(n)}) \\ c/E(-\sigma_N^{(n)} - \sigma_{N-1}^{(n)}) + v_N^{(n)} + v_{N-1}^{(n)} \end{pmatrix} \quad (C.3)$$
WFEM free-flight matrix

From previous developments, it is known that

$$ \left( \begin{array}{c} \boldsymbol{\sigma}^{(n+1)} \\ \mathbf{v}^{(n+1)} \end{array} \right) = \frac{1}{2} \left[ \begin{array}{cc} \mathbf{A}_1 & \mathbf{A}_2 \\ \mathbf{A}_3 & \mathbf{A}_4 \end{array} \right] \left( \begin{array}{c} \boldsymbol{\sigma}^{(n)} \\ \mathbf{v}^{(n)} \end{array} \right) \Leftrightarrow \mathbf{q}^{n+1} = \mathbf{A}_f \mathbf{q}^n \qquad (C.4) $$

where $\mathbf{A}_1 = \mathbf{A}_g(1, 1, 1, -1)$, $\mathbf{A}_2 = E/c\mathbf{A}_g(1, 1, -1, 1)$, $\mathbf{A}_3 = c/E\mathbf{A}_g(-1, 1, -1, -1)$ and $\mathbf{A}_4 = \mathbf{A}_g(-1, 1, 1, 1)$ with

$$ \mathbf{A}_g(a, b, c, d) = \begin{bmatrix} a & b & 0 & \cdots & 0 & 0 & 0 \\ c & 0 & b & & 0 & 0 & 0 \\ 0 & c & 0 & & 0 & 0 & 0 \\ \vdots & \ddots & \ddots & \ddots & & & & \\ 0 & 0 & 0 & & 0 & b & 0 \\ 0 & 0 & 0 & & c & 0 & b \\ 0 & 0 & 0 & \cdots & 0 & c & d \end{bmatrix}. \qquad (C.5) $$

We can perform a similar computation to show that $C_{x,t_{dc}}$ leads to $\mathbf{A}_c^p$ with

$$ \mathbf{A}_c = \frac{1}{2} \begin{bmatrix} \mathbf{A}_5 & \mathbf{A}_6 \\ \mathbf{A}_7 & \mathbf{A}_8 \end{bmatrix} \qquad (C.6) $$

and $\mathbf{A}_5 = \mathbf{A}_g(1, 1, 1, 1)$, $\mathbf{A}_6 = E/c\mathbf{A}_g(1, 1, -1, -1)$, $\mathbf{A}_7 = c/E\mathbf{A}_g(-1, 1, -1, 1)$ and $\mathbf{A}_8 = \mathbf{A}_g(-1, 1, 1, -1)$.

### C.1.3 Expression of the B matrix

With regards to the hyperbolic system of conservation laws we got B. We can write an discretized equivalent for B which is

$$ \hat{\mathbf{B}} = - \begin{bmatrix} \mathbf{0} & E\mathbf{I}_N \\ c^2/E\mathbf{I}_N & \mathbf{0} \end{bmatrix}. \qquad (C.7) $$

### C.1.4 Discrete Derivation

Regarding the periodic extension of the functions, we have to consider different kind of discrete derivations (as done usually in finite different schemes) to differentiate stress and velocity at the Nth component. Indeed,

$$ f \in F([0, L], \mathbb{R}), \quad \forall n \in [|1, N-1|], \quad f'(x_n) \approx \frac{1}{\Delta x}(f(x_{n+1}) - f(x_n)) \qquad (C.8) $$

But for the Nth values of f, we have to use the periodic extensions.

### Free-flight

If we consider the free-flight and the corresponding periodic extensions for $\sigma$ and $v$, we can write

$$ q'_{ff} = \frac{1}{\Delta x} \begin{bmatrix} \mathbf{DFF}_1 & 0 \\ 0 & \mathbf{DFF}_2 \end{bmatrix} \begin{pmatrix} \boldsymbol{\sigma} \\ \mathbf{v} \end{pmatrix} = \frac{1}{\Delta x} \mathbf{D}_{ff} q_{ff} \qquad (C.9) $$
where $D_{ff}$ is a $2N \times 2N$ matrix such that

$$ \mathbf{DFF}_1 = \begin{bmatrix} -1 & 1 & 0 & \cdots \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \cdots & -1 & 1 \\ 0 & \cdots & -1 & -1 \end{bmatrix} \quad \text{and} \quad \mathbf{DFF}_2 = \begin{bmatrix} -1 & 1 & 0 & \cdots \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \cdots & -1 & 1 \\ 0 & \cdots & 1 & -1 \end{bmatrix}. \qquad (C.10) $$

The last row of $\mathbf{DFF}_1$ and $\mathbf{DFF}_2$ are explained by the fact that during free-flight $\sigma(L - x) = -\sigma(L + x)$ and $v(L - x) = v(L + x)$.

### Contact phase

If we consider the contact phase and the corresponding periodic extensions for $\sigma$ and $v$, we can write the following:

$$ q'_{c} = \frac{1}{\Delta x} \begin{bmatrix} \mathbf{DFF}_{2} & 0 \\ 0 & \mathbf{DFF}_{1} \end{bmatrix} \begin{pmatrix} \boldsymbol{\sigma} \\ \mathbf{v} \end{pmatrix} = \frac{1}{\Delta x} \mathbf{D}_{c} \mathbf{q}_{c} \qquad (C.11) $$

where $D_c$ is also a $2N \times 2N$ and the use of $\mathbf{DFF}_2$ and $\mathbf{DFF}_1$ are explained by the fact that during contact $\sigma(L - x) = \sigma(L + x)$ and $v(L + x) = -v(L - x)$.

### C.1.5 Discretization of integral parts of the operator

First, notice that the following holds:

$$ \int_{0}^{t_f} \frac{\mathbf{e}_2^{\top} \mathcal{F}(\delta \psi_0)(L, \tau)}{\mathbf{e}_2^{\top} \mathcal{F}(\mathbf{q}_0)(L, t_f)} d\tau = \frac{\mathbf{e}_2^{\top}}{\mathbf{e}_2^{\top} \mathcal{F}(\mathbf{q}_0)(L, t_f)} \int_{0}^{t_f} \mathcal{F}(\delta \psi_0)(L, \tau) d\tau \qquad (C.12) $$

**Discretization of** $\int_{0}^{t_f} \mathcal{F}(\delta \psi_0)(L, \tau) d\tau$

Using the property of integration and the fact that $t_f = m\Delta t$, we can write

$$ \int_{0}^{t_f} \mathcal{F}(\delta \psi_0)(L, \tau) d\tau = \sum_{k=0}^{m-1} \int_{k\Delta t}^{(k+1)\Delta t} \mathcal{F}(\delta \psi_0)(L, \tau) d\tau \qquad (C.13) $$

$$ \frac{\Delta t}{2} \approx \sum_{k=0}^{m-1} \mathcal{F}(\delta \psi_0)(\bullet, k\Delta t) + \mathcal{F}(\delta \psi_0)(\bullet, (k+1)\Delta t) \qquad (C.14) $$

$$ \approx \frac{\Delta t}{2} \left[ \sum_{k=0}^{m-1} \mathbf{A}_{f}^{k} + \mathbf{A}_{f}^{k+1} \right] \delta \psi_{0} = \frac{\Delta t}{2} \mathbf{F} \delta \psi_{0} \qquad (C.15) $$

where $\delta \psi_0$ should be discretized along the interval $[0, L]$. Here we use trapezoidal rule in Equation (C.14) (illustration in Figure C.1) for the integral in order to keep as much information about the motion as possible. Indeed, this method allows to keep informations given by $A_f^0$ and $A_f^m$. It is well-known that this type of integration rule converges for Lipschitz continuous functions, which is the case here for $\delta\sigma$, see Equation (2.18c).

### Discrete projections on state vector components

Let us remind that $\mathbf{e}_2^\top \mathcal{F}$ selects the component $v$. In the discretized framework, the equivalent will be the projection on the $2N^{th}$ coordinate of the discretized state vectors $\delta \psi_T$ or $\mathbf{q}_t$. We
<figure><img src="image_19.png" /><figcaption>Figure C.1: Trapezoidal quadrature rule function: (—) f and (—) approximation of f.</figcaption></figure>

will denote $e_{2N}^\top$ the projection in the discretized framework. $e_i$ is the $i^{th}$ canonical basis vector of $\mathbb{R}^{2N}$. Using the previous results, we can write

$$ \int_{0}^{t_f} \frac{e_2^\top \mathcal{F}(\delta\psi_0)(L, \tau)}{e_2^\top \mathcal{F}(\mathbf{q}_0)(L, t_f)} d\tau \approx \frac{\Delta t}{2} \frac{\langle e_{2N} | \mathbf{F} \delta\psi_0 \rangle}{\langle e_{2N} | \mathbf{A}_f^m \mathbf{q}_0 \rangle} \quad (C.16) $$

where $\langle \cdot | \cdot \rangle$ holds for canonical inner-product of $\mathbb{R}^{2N}$. We could also use the known result for the mode considered : $\langle e_{2N} | \mathbf{A}_f^m \mathbf{q}_0 \rangle = v_{NSM}(L, t_f)$ to fasten computation.

**Discretization of** $\int_0^L \mathcal{F}(\delta\psi_0)(s, t_f) ds$

We saw in Equation (2.21) that two possible expressions of $u(L, t_f)$ exist. Here it is depicted how to discretize the second expression. The faster computation will be considered:

$$ \int_{0}^{L} \mathcal{F}(\delta\psi_0)(s, t_f) ds = \sum_{k=1}^{N} \int_{(k-1)\Delta x}^{k\Delta x} \mathcal{F}(\delta\psi_0)(s, t_f) ds \approx \Delta x \sum_{k=1}^{N} \mathcal{F}(\delta\psi_0)(k\Delta x, t_f) \quad (C.17) $$

$$ \approx \Delta x \sum_{k=1}^{N} (\mathbf{A}_f^m \delta\psi_0)_k = \Delta x \sum_{k=1}^{N} \langle e_k | \mathbf{A}_f^m \delta\psi_0 \rangle. \quad (C.18) $$

Previous results yield

$$ \int_{0}^{L} \frac{1}{E} \frac{e_{1}^{\top} \mathcal{F}_{s, t_f}(\delta\psi_{0})}{e_{2}^{\top} \mathcal{F}_{L, t_f}(\mathbf{q}_{0})} ds \approx \frac{\Delta x}{E} \sum_{k=1}^{N} \frac{\langle e_{k} | \mathbf{A}_{f}^{m} \delta\psi_{0} \rangle}{\langle e_{2N} | \mathbf{A}_{f}^{m} \mathbf{q}_{0} \rangle}. \quad (C.19) $$

### C.1.6 Discretized operator

#### A first version of the discretized operator

Using all the previous subsections we can write

$$ \begin{aligned} \delta\psi_T &= \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{\Delta t}{2\Delta x} \frac{\langle e_{2N} | \mathbf{F} \delta\psi_0 \rangle}{\langle e_{2N} | \mathbf{A}_f^m \mathbf{q}_0 \rangle} (\mathbf{A}_c^p \hat{\mathbf{B}} \mathbf{A}_f^m \mathbf{D}_{ff} \mathbf{q}_0 - \hat{\mathbf{B}} \mathbf{A}_c^p \mathbf{D}_c \mathbf{A}_f^m \mathbf{q}_0) \\ &= \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{\Delta t}{2\Delta x} \frac{\langle e_{2N} | \mathbf{F} \delta\psi_0 \rangle}{\langle e_{2N} | \mathbf{A}_f^m \mathbf{q}_0 \rangle} (\mathbf{A}_c^p \hat{\mathbf{B}} \mathbf{A}_f^m \mathbf{D}_{ff} \mathbf{q}_0 - \hat{\mathbf{B}} \mathbf{A}_c^p \mathbf{D}_c \mathbf{A}_f^m \mathbf{q}_0) \\ &= \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{2} \frac{\Delta t}{\Delta x} \langle e_{2N} | \mathbf{F} \delta\psi_0 \rangle \frac{(\mathbf{A}_c^p \hat{\mathbf{B}} \mathbf{A}_f^m \mathbf{D}_{ff} - \hat{\mathbf{B}} \mathbf{A}_c^p \mathbf{D}_c \mathbf{A}_f^m)}{\langle e_{2N} | \mathbf{A}_f^m \mathbf{q}_0 \rangle} \mathbf{q}_0. \end{aligned} \quad (C.20) $$
In order to write it shorter and more clearly, we define

$$ \mathbf{x}_{\text{NSM}} = \frac{(\mathbf{A}_c^p \hat{\mathbf{B}} \mathbf{A}_f^m \mathbf{D}_{\text{ff}} - \hat{\mathbf{B}} \mathbf{A}_c^p \mathbf{D}_c \mathbf{A}_f^m) \mathbf{q}_0}{\langle \mathbf{e}_{2N} | \mathbf{A}_f^m \mathbf{q}_0 \rangle} \quad (C.21) $$

which is a vector of $\mathbb{R}^{2N}$. This vector carries all the information about the mode considered (free phase, closure time but also switching time), and

$$ \begin{aligned} \delta\psi_T &= \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{2} \frac{\Delta t}{\Delta x} \langle \mathbf{e}_{2N} | \mathbf{F} \delta\psi_0 \rangle \mathbf{x}_{\text{NSM}} \\ &= \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{2c} \langle \mathbf{e}_{2N} | \mathbf{F} \delta\psi_0 \rangle \mathbf{x}_{\text{NSM}} = \mathbf{L}(\delta\psi_0) \end{aligned} \quad (C.22) $$

where $\mathbf{L}$ is a notation which embed the discretized linear operator of first return map. We may compute its eigenvalues with numerical methods as soon as we can express a matrix expression of it which will be used to derive a stability analysis.

Using the fact that $\langle \mathbf{x} | \mathbf{y} \rangle = \mathbf{x}^\top \mathbf{y}$ and that left and right scalar multiplication are the same, we get:

$$ \begin{aligned} \delta\psi_T &= \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{2c} \langle \mathbf{e}_{2N} | \mathbf{F} \delta\psi_0 \rangle \mathbf{x}_{\text{NSM}} = \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{2c} \mathbf{x}_{\text{NSM}} \langle \mathbf{e}_{2N} | \mathbf{F} \delta\psi_0 \rangle \\ &= \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{2c} \mathbf{x}_{\text{NSM}} (\mathbf{e}_{2N}^\top \mathbf{F} \delta\psi_0) = \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{2c} (\mathbf{x}_{\text{NSM}} \mathbf{e}_{2N}^\top \mathbf{F}) \delta\psi_0 \\ &= (\mathbf{A}_c^p \mathbf{A}_f^m + \frac{1}{2c} \mathbf{x}_{\text{NSM}} \mathbf{e}_{2N}^\top \mathbf{F}) \delta\psi_0 = \mathbf{R}_1(m, p, \mathbf{q}_0) \delta\psi_0. \end{aligned} \quad (C.23) $$

## A second version of the discretized operator

Using Equation (C.19) instead of Equation (C.16) we can have an other way to write the discretized operator

$$ \delta\psi_T = \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{E} \sum_{k=1}^{N} \langle \mathbf{e}_k | \mathbf{A}_f^m \delta\psi_0 \rangle \mathbf{x}_{\text{NSM}} = \mathbf{L}(\delta\psi_0) \quad (C.24) $$

Moreover, having two different expressions of this operator could be helpful to check results.

Via the property of the inner product, we can show that

$$ \begin{aligned} \delta\psi_T &= \mathbf{A}_c^p \mathbf{A}_f^m \delta\psi_0 + \frac{1}{E} \sum_{k=1}^{N} \langle \mathbf{e}_k | \mathbf{A}_f^m \delta\psi_0 \rangle \mathbf{x}_{\text{NSM}} \\ &= (\mathbf{A}_c^p \mathbf{A}_f^m + \frac{1}{E} \mathbf{x}_{\text{NSM}} (\sum_{k=1}^{N} \mathbf{e}_k)^\top \mathbf{A}_f^m) \delta\psi_0 = \mathbf{R}_2(m, p, \mathbf{q}_0) \delta\psi_0. \end{aligned} \quad (C.25) $$

## C.2 Convergence of the discretized operator

Since the quadrature rule for integration of Lipschitz continuous functions and discrete derivation are usual tools in finite difference methods, the convergence of this discretized operator towards the explicit time/space one only relies on the convergence of the WFEM's solutions towards the exact solutions of the problem introduced earlier in Section 1.2.4. So we need the WFEM scheme to be proven stable and consistent. For any further proofs or deeper explanation, we refer to [3]. Moreover, we use the known CFL condition ($c\Delta t \le \Delta x$) about time step and space step in order to avoid numerical dispersion as explain in [10].
### C.2.1 Consistency of WFEM

We consider a hyperbolic system of conservative laws of the form

$$ \forall t \ge 0, \quad \forall x \in [0, L], \quad \partial_t \mathbf{q} + \partial_x f(\mathbf{q}) = 0 \qquad (C.26) $$

From [3], the WFEM is consistent if the numerical flux $\tilde{f}$ is consistent with Equation (C.26) and the flux $f$ is Lipschitz continuous.

* We call numerical flux $\tilde{f}$ the function of $\mathbb{R}^2$ such that the WFEM numerical scheme can be written $\forall i \in [1, N]$, $\forall n \in \mathbb{N}$:

$$ \mathbf{q}_i^{(n+1)} - \mathbf{q}_i^{(n)} + \frac{\Delta t}{\Delta x} (\tilde{f}(\mathbf{q}_i^{(n)}, \mathbf{q}_i^{(n+1)}) - \tilde{f}(\mathbf{q}_i^{(n-1)}, \mathbf{q}_i^{(n)})) = 0. \qquad (C.27) $$

* We say that a function $h$ of $\mathbb{R}^2$ is Lipschitz continuous if

$$ \exists K \in \mathbb{R}^{+}, \quad \forall (\mathbf{x}, \mathbf{y}) \in \mathbb{R}^{2} \times \mathbb{R}^{2}, \quad \|h(\mathbf{x}) - h(\mathbf{y})\|_{\mathbb{R}^{2}} \le K \|\mathbf{x} - \mathbf{y}\|_{\mathbb{R}^{2}}. \qquad (C.28) $$

* We say that a numerical flux $\tilde{f}$ is consistent with Equation (C.26) if

$$ \forall \mathbf{q} \in \mathbb{R}^2, \quad \tilde{f}(\mathbf{q}, \mathbf{q}) = f(\mathbf{q}). \qquad (C.29) $$

#### Numerical flux $\tilde{f}$ of WFEM

If we refer to Equation (C.1) the following holds $\forall i \in [1, N]$, $\forall n \in \mathbb{N}$

$$ \mathbf{q}_i^{(n+1)} - \mathbf{q}_i^{(n)} + \frac{\Delta t}{\Delta x} [\tilde{f}(\mathbf{q}_{i+1}^{(n)}, \mathbf{q}_i^{(n)}) - \tilde{f}(\mathbf{q}_i^{(n)}, \mathbf{q}_{i-1}^{(n)})] = 0 \qquad (C.30) $$

where

$$ \tilde{f} : (\mathbf{x}_1, \mathbf{x}_2) \rightarrow -\frac{1}{2} \begin{bmatrix} c & E \\ c^2/E & c \end{bmatrix} \mathbf{x}_1 - \frac{1}{2} \begin{bmatrix} -c & E \\ c^2/E & -c \end{bmatrix} \mathbf{x}_2 = \mathbf{F}_1 \mathbf{x}_1 + \mathbf{F}_2 \mathbf{x}_2. \qquad (C.31) $$

Since we have just shown the shape of $\tilde{f}$, we can now show that the WFEM scheme is consistent.

#### Consistency of the WFEM

We consider a given $\mathbf{q} \in \mathbb{R}^2$ such that

$$ \begin{aligned} \tilde{f}(\mathbf{q}, \mathbf{q}) &= -\frac{1}{2} \begin{bmatrix} c & E \\ c^2/E & c \end{bmatrix} \mathbf{q} - \frac{1}{2} \begin{bmatrix} -c & E \\ c^2/E & -c \end{bmatrix} \mathbf{q} \\ &= - \begin{bmatrix} 0 & E \\ c^2/E & 0 \end{bmatrix} \mathbf{q} = \mathbf{B}\mathbf{q} = f(\mathbf{q}). \end{aligned} \qquad (C.32) $$

We have shown that the numerical flux is consistent with Equation (2.3) because in our case the function $f$ is the matrix product of $\mathbf{q}$ by matrix $\mathbf{B}$.

We consider a given $(\mathbf{q}_1, \mathbf{q}_2) \in \mathbb{R}^2 \times \mathbb{R}^2$. If we consider a sub multiplicative matrix norm (such as $||\mathbf{M}||_{\mathcal{M}_{2,2}} = \sup_{(i,j)}|\mathbf{M}_{i,j}|$), it is easy to check that the flux $f$ is Lipschitz continuous:

$$ \begin{aligned} ||f(\mathbf{q}_1) - f(\mathbf{q}_2)||_{\mathbb{R}^2} &= ||\mathbf{B}\mathbf{q}_1 - \mathbf{B}\mathbf{q}_2||_{\mathbb{R}^2} \\ &\le ||\mathbf{B}||_{\mathcal{M}_{2,2}} ||\mathbf{q}_1 - \mathbf{q}_2||_{\mathbb{R}^2} \le \sup_{E,c^2/E} ||\mathbf{q}_1 - \mathbf{q}_2||_{\mathbb{R}^2}. \end{aligned} \qquad (C.33) $$

With the previous development we proved that the WFEM is consistent. Although we have a switching at the boundary, it does not matter since the consistency only relies on Equation (C.26).
### C.2.2 Stability of WFEM

#### Definitions

We call the total variation of a grid function q

$$TV(\mathbf{q}^{(n)}) = \sum_{i=-\infty}^{+\infty} |\mathbf{q}_i^{(n)} - \mathbf{q}_{i-1}^{(n)}| = \sum_{i=-\infty}^{+\infty} |\mathbf{v}_i^{(n)}|. \qquad (C.34)$$

We call total variation diminishing (TVD), a two-level method which satisfy for any set of data

$$TV(\mathbf{q}^{(n+1)}) \leq TV(\mathbf{q}^{(n)}). \qquad (C.35)$$

According to [10], a TVD method which is consistent with a Lipschitz continuous flux f is convergent.

#### Sufficient condition for TVD method

Here we are going to prove that the flux allow the variation across the i-th cell to be conserved during time

$$\mathbf{v}_i^{(n)} = \mathbf{q}_i^{(n)} - \mathbf{q}_{i-1}^{(n)} \qquad (C.36)$$

$$= \mathbf{q}_i^{(n+1)} + \frac{\Delta t}{\Delta x} [\tilde{\mathbf{f}}(\mathbf{q}_{i+1}^{(n)}, \mathbf{q}_i^{(n)}) - \tilde{\mathbf{f}}(\mathbf{q}_i^{(n)}, \mathbf{q}_{i-1}^{(n)})] \qquad (C.37)$$

$$= \mathbf{v}_i^{(n+1)} + \frac{\Delta t}{\Delta x} [\mathbf{F}_1(\mathbf{q}_{i+1}^{(n)} - 2\mathbf{q}_i^{(n)} + \mathbf{q}_{i-1}^{(n)}) + \mathbf{F}_2(\mathbf{q}_i^{(n)} - 2\mathbf{q}_{i-1}^{(n)} + \mathbf{q}_{i-2}^{(n)})] \qquad (C.38)$$

$$= \mathbf{v}_i^{(n+1)} + \frac{\Delta t}{\Delta x} [\mathbf{F}_1(\mathbf{v}_{i+1}^{(n)} - \mathbf{v}_i^{(n)}) + \mathbf{F}_2(\mathbf{v}_i^{(n)} - \mathbf{v}_{i-1}^{(n)})] \qquad (C.39)$$

$$\Leftrightarrow \mathbf{v}_i^{(n+1)} - \mathbf{v}_i^{(n)} + \frac{\Delta t}{\Delta x} [\tilde{\mathbf{f}}(\mathbf{v}_{i+1}^{(n)}, \mathbf{v}_i^{(n)}) - \tilde{\mathbf{f}}(\mathbf{v}_i^{(n)}, \mathbf{v}_{i-1}^{(n)})] = 0. \qquad (C.40)$$

Using Equation (C.40) and the fact that $\mathbf{F}_2 - \mathbf{F}_1 = c\mathbf{I}$ we can write

$$\forall (i, n) \in \mathbb{Z} \times \mathbb{N}, \quad \mathbf{v}_i^{(n+1)} = \frac{\Delta t}{\Delta x} \mathbf{F}_2 \mathbf{v}_{i-1}^{(n)} - \frac{\Delta t}{\Delta x} \mathbf{F}_1 \mathbf{v}_{i+1}^{(n)} \qquad (C.41)$$

which implies

$$|\mathbf{v}_i^{(n+1)}| \leq \frac{\Delta t}{2\Delta x} \sup(E, c^2/E, c)(|\mathbf{v}_{i-1}^{(n)}| + |\mathbf{v}_{i+1}^{(n)}|). \qquad (C.42)$$

By summing this inequality and rearranging the right side we obtain

$$\sum_{i \in \mathcal{I}} |\mathbf{v}_i^{(n+1)}| \leq \frac{\Delta t}{\Delta x} \sup(E, c^2/E, c) \sum_{i \in \mathcal{I}} |\mathbf{v}_i^{(n)}| \qquad (C.43)$$

where $\mathcal{I}$ is the finite set of $\mathbb{Z}$ where $\mathbf{v}_i^{(\bullet)}$ is non zero.

If the quantity denoted by $\frac{\Delta t}{\Delta x} \sup(E, c^2/E, c)$ is less or equal to unity, then the scheme is a TVD method, so the WFEM is convergent in the weak sense [3]. This can be understood by the fact that the total variation of each cell is consequently bounded by its initial value.
Bibliography

[1] V. Astashev and V. Krupenin. Longitudinal vibrations of a thin rod interacting with an immobile limiter. *Journal of Machinery Manufacture and Reliability*, 36(6):535–541, 2007.

[2] F. Bouchut. Quasilinear systems and conservation laws. In *Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws*, Frontiers in Mathematics, chapter 1, pages 1–11. Birkhäuser Basel, 2004.

[3] F. Bouchut. Conservative schemes. In *Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws*, Frontiers in Mathematics, chapter 2, pages 13–64. Birkhäuser Basel, 2004.

[4] C. Briot and J.-C. Bouquet. *Théorie des fonctions doublement périodiques et, en particulier, des fonctions elliptiques*. Bibliothèque nationale de France, 1859.

[5] H. Cabannes. Mouvement d'une corde vibrante en présence d'un obstacle ponctuel fixé. *Comptes-rendus des séances de l'Académie des sciences*, 298(15):613–616, 1984.

[6] S. de Bièvre, F. Genoud, and S.R. Nodari. Orbital Stability: Analysis Meets Geometry. In *Nonlinear Optical and Atomic Systems*, number 2146 in Lecture Notes in Mathematics, pages 147–273. Springer, 2015.

[7] L. Evans. *Partial Differential Equations*. AMS, 2010.

[8] F. John. Second Order Linear Equations with Constant Coefficients. In *Partial Differential Equations*, chapter 3, pages 87–185. Springer, 1975.

[9] 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.

[10] R. LeVeque. *Nonlinear Stability*. In *Numerical Methods for Conservation Laws*, chapter 15, pages 158–172. Birkhäuser, 1992.

[11] C. Yoong, A. Thorin, and M. Legrand. The Wave Finite Element Method applied to a one-dimensional linear elastodynamic problem with unilateral constraints. In *ASME IDETC & CIE Conference*, Boston, United States, 2015.

[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.