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

HAL
open science

# Harmonic Balance based Nonsmooth Modal Analysis of unilaterally constrained discrete systems

Tianzheng Lu, Mathias Legrand

► To cite this version:

Tianzheng Lu, Mathias Legrand. Harmonic Balance based Nonsmooth Modal Analysis of unilaterally constrained discrete systems. Nonlinear Dynamics, 2023, 112, pp.1619-1640. 10.1007/s11071-023-09014-4. hal-03970830v2

HAL Id: hal-03970830

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

Submitted on 28 Sep 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Harmonic Balance based Nonsmooth Modal Analysis of unilaterally constrained discrete systems

Tianzheng Lu & Mathias Legrand

Structural Dynamics and Vibration Laboratory, Department of Mechanical Engineering, McGill University, Canada

Nonsmooth Modal Analysis of a unilaterally constrained one-dimensional bar with constant cross-sectional area was recently proposed. The corresponding formulation took advantage of the d'Alembert solution available for such systems and does not require any space semi-discretization of the governing equations. However, it is unable to cope with non-constant cross-sectional area bars, for instance. The present work suggests a formulation relying on various space semi-discretization methodologies (such as finite elements, Rayleigh-Ritz techniques, component mode synthesis, modal superposition and other reduced-order models) where the complementarity Signorini condition, reflecting the unilateral contact constraint, is enforced in a weighted residual sense in time through the Harmonic Balance Method. Importantly, an impact law, such as Newton's impact law for instance, classically required for uniqueness purposes in a space semi-discrete dynamical framework, is here explicitly ignored in the proposed formulation and is, instead, implicitly satisfied in a weighted residual sense. As required for the existence of periodic solutions, the predicted vibratory responses would then correspond to an energy-preserving impact law in case the latter had to be explicitly implemented. Periodic responses are investigated in the form of classical energy-frequency backbone curves along with the associated displacement fields. It is found that for the constant cross-section benchmark system, the results compare well with existing works and the proposed methodology stands as a viable option in the field of interest when semi-discretization in space is required.

## Contents

1 Introduction

2 System of interest

3 Solution method

3.1 Space semi-discretization and shape functions .................................... 3

3.2 Time-domain discretization for periodic solutions .................................... 4

3.3 Unilateral contact condition ................................................................. 5

4 Convergence analysis

4.1 Transfer function .................................................................................... 6

4.2 Unilateral contact residual ........................................................................ 8

4.3 Unilateral contact condition ........................................................................ 9

5 Nonsmooth modal analysis

5.1 Benchmark case: constant cross-section .................................................... 11

5.2 Non-constant cross-section ........................................................................ 12

6 Conclusion

A Periodic solutions of one-dof systems

A.1 Forced one-dof system: the bouncing ball .................................................... 15

A.2 Autonomous one-dof system ........................................................................ 16

B Linear mode expansion convergence analysis

## References

## 1 Introduction

Nonsmooth modal analysis (NSA) is one incarnation of nonlinear modal analysis dedicated to nonsmooth systems. A few numerical schemes have been proposed to perform NSA of continuous systems. For instance, the Wave Finite Element Method (WFEM) with a switch in the boundary conditions [34] could partially solve the case of a one-dimensional bar system. Another solution based on the Time-Domain Boundary Element Method (TD-BEM) is also available [31]. However, for the problem at hand, it requires challenging computations involving the initial conditions and the attendant space semi-discretization of the domain of interest: this is not optimal since it heavily reduces the computational efficiency of TD-BEM. In [15], the authors suggested a numerical scheme which combines the boundary element method (BEM) to the Harmonic Balance Method (HBM), equivalent, for the system of interest, to numerically approximating the periodic solutions to the governing equations expressed in the frequency domain. Similarly to [34], the solution strategy is able to perform NSAs for bar systems with constant cross-section by taking advantage of the existence of the d'Alembert solution expressed in the frequency domain. However, for one-dimensional systems with space-dependent parameters, the methodology would fail. Accordingly,
the present work suggests a family of methods combining classical space semi-discretization along with HBM in time. Such space semi-discretization has two major common consequences on the governing equations: (1) the need of an impact law at the contact interface where “discrete” mass is generated by the discretization scheme and (2) the detrimental perturbation of the spectrum of the system, all the more in the high-frequency range. Although common, this second aspect is a challenge for the application of interest since it is known that the targeted nonlinear modal motions are not everywhere differentiable progressive waves [34].

The present paper is organized as follows. In Section 2, the one-dimensional non-homogeneous system and attendant governing equations are briefly introduced. Section 3 explains the proposed solution method, which includes FEM-like semi-discretization in space, HBM in time and a weighted residual enforcement of the Signorini condition at the contact boundary. Thorough convergence analysis is detailed in Section 4. Nonsmooth modal analysis is conducted in Section 5 to illustrate the capabilities of the developed tool.

## 2 System of interest

As an extension of [15], the present work explores the modal response of one-dimensional academic systems with unilateral conditions on the boundary, such as the one illustrated in Figure 1. As already said, the main novelty

<figure><img src="image_3.png" /><figcaption>Figure 1: System of interest: a one-dimensional bar with non-constant cross-sectional area.</figcaption></figure>

compared to [15] is on the required semi-discretization of the governing equations in space, in addition to the already existing semi-discretization in time implemented via HBM. Young’s modulus E and mass density ρ are assumed space-independent but the cross-section area A(x) is space-dependent. As already said, the novelty compared to [15] is on the needed semi-discretization of the governing equations in space, in addition to the already existing semi-discretization in time. The length of the bar system is denoted by L. Young’s modulus E and mass density ρ are assumed space-independent but the cross-section area A(x) is space-dependent.

A brief non-dimensional analysis is introduced to facilitate the exposition of the work and analysis. The non-dimensional quantities come with an overbar notation such that $\bar{\bullet}$ is the non-dimensional version of $\bullet$. Are thus defined non-dimensional space $\bar{x} = x/L$, time $\bar{t} = t/\tau$ and displacement $\bar{u} = u/L_1$ with the characteristic time $\tau = L/c$ where $c = \sqrt{E/\rho}$. Derivatives of u are found using the chains rule: $u_x = L\bar{u}_{\bar{x}}\bar{x}_x = \bar{u}_{\bar{x}}$ and $u_t = L\bar{u}_{\bar{t}}\bar{t}_t = L\bar{u}_{\bar{t}}/\tau$. Higher derivatives can be found in a similar way: $u_{xx} = \bar{u}_{\bar{x}\bar{x}}/L$ and $u_{tt} = L\bar{u}_{\bar{t}\bar{t}}/\tau^2$. Meanwhile, non-dimensional $\bar{A}(x)$ is introduced via $\bar{A}(x) = A(x)/A_0$ where $A_0$ is the average cross-section area of the bar. In the remainder, the upper bar notation is dropped and all considered quantities are non-dimensional.

With the introduced notation, the governing equation of the system illustrated in Figure 1 reads

$$ (A(x)u_x(x,t))_x - A(x)u_{tt}(x,t) = 0, \quad x \in ]0; 1[ \qquad (1) $$

The system is not solicited by any other forces. The boundary conditions are as follows:

**Dirichlet** System clamped at $x = 0$ so that a homogeneous Dirichlet boundary condition $u(0, t) = 0$ applies.

**Signorini** Unilateral contact at $x = 1$. Defining the gap $g(t) = g_0 - u(1, t)$ separating the bar tip to the rigid foundation where $g_0$ is the initial gap distance, the Signorini condition takes the classical form $g(t) \ge 0$, $p(t) \le 0$, and $g(t)p(t) = 0$ where $p(t)$ is the contact force acting on the bar at $x = 1$. This boundary condition could equivalently be expressed as the equality [1, 28]

$$ r(u(1,t), p(t)) = p(t) + \max[\alpha(u(1,t) - g_0) - p(t), 0] = 0. \qquad (2) $$

among other possible equality-based expressions. In Equation (2), $\alpha$ is an arbitrary strictly positive number that can be tuned to numerically ensure that the various terms in the expression have similar magnitudes.

To the authors’ knowledge, little if anything is known concerning the existence and uniqueness of solutions for the non-constant cross-section bar with a unilateral contact condition at one end, in the context of continuum mechanics. This is in contrast with the constant cross-section bar, for which results are available. However, the solution scheme proposed in the present work implements a discretization in space and time (or frequency) and taking the limits in space and time to recover results in the continuous framework seems quite challenging.

## 3 Solution method

The present work targets the modal response of the system detailed above, as an extension of previous works investigating bars with constant cross-sections [15, 30, 34]. Space semi-discrete systems were already explored in
[13, 29] where the proposed formulation had to include an energy-preserving impact law for the contacting mass. This resulted in solutions exhibiting chattering-like responses with the additional difficulty of counting the number of impacts per period in the solution method. Methodologies with regularization of the unilateral contact conditions also exist, as discussed in [15]. In the present contribution, we adapt the solution strategy introduced in [15] to discrete systems where the unilateral contact conditions are neither regularized nor complemented with an impact law but satisfied in a weak sense only over a period of the energy-preserving modal motion.

## 3.1 Space semi-discretization and shape functions

Since periodic responses are targeted, a frequency-domain formulation has appealing features and a Fourier transform of the unknown displacement and contact force is used [15], which yields the Helmholtz equation

$$ (A(x)\hat{u}_x(x, \omega))_x + \omega^2 A(x)\hat{u}(x, \omega) = 0, \quad x \in ]0; 1[. \quad (3) $$

Note that this transformation is not really needed since the classical HBM applied on an Ordinary Differential Equation would result in the same formulation<fnref n="1" />. However, it simplifies the exposition of the derivations. Equation (3) is first expressed in its weak form in space: Find $\hat{u}$ such that $\forall v$ with $v(0) = 0$,

$$ \int_{0}^{1} v_{x}(x) A(x) \hat{u}_{x}(x, \omega) \mathrm{d}x - \omega^{2} \int_{0}^{1} v(x) A(x) \hat{u}(x, \omega) \mathrm{d}x = v(1) A(1) \hat{u}_{x}(1, \omega). \quad (4) $$

The Fourier transform of the contact force $p(t)$ is $\hat{p}(\omega)$ satisfying the identity

$$ \hat{p}(\omega) = A(1)\hat{u}_x(1, \omega) \quad (5) $$

used in the coming developments. The two quantities $\hat{u}$ and $\hat{p}$ are linearly independent and considered as two unknowns in the sequel. Space discretization of the trial and test functions reads

$$ \hat{u}(x, \omega) \approx \hat{u}_N(x, \omega) = \sum_{i=1}^{N} \phi_i(x) \hat{u}_i(\omega) \quad \text{and} \quad v(x) \approx v_N(x) = \sum_{i=1}^{N} \phi_i(x) v_i \quad (6) $$

where $\phi_i(x)$, $i = 1, \dots, N$ are the chosen shape functions stored in the vector $\phi(x) = [\phi_1(x), \dots, \phi_N(x)]^\top$. Inserting $\hat{u}_N(x, \omega)$ and $v_N(x)$ from Equation (6) into Equation (4) leads to the system of linear equations

$$ \int_{0}^{1} \left( A(x) \phi_{x}(x)^{\top} \phi_{x}(x) - \omega^{2} A(x) \phi_{x}(x)^{\top} \phi_{x}(x) \right) \mathrm{d}x \hat{u}(\omega) = \hat{p}(\omega) \quad (7) $$

with the frequency-domain displacement vector<fnref n="2" /> $\hat{u}(\omega) = [\hat{u}_1(\omega), \dots, \hat{u}_N(\omega)]^\top$ and the contact force vector $\hat{p}(\omega) = \phi(1)\hat{p}(\omega)$. Evaluating the integral leads to

$$ (\mathbf{K} - \omega^2 \mathbf{M})\hat{u}(\omega) = \hat{p}(\omega). \quad (8) $$

Equation (8) linearly connects the unknown discretized displacements to the unknown unilateral contact force. Let us define $\mathbf{G}(\omega) = \mathbf{K} - \omega^2 \mathbf{M}$ so that Equation (8) becomes

$$ \mathbf{G}(\omega)\hat{u}(\omega) = \hat{p}(\omega) \quad (9) $$

that is $\hat{u}(\omega) = \mathbf{G}(\omega)^{-1}\hat{p}(\omega)$ as soon as $\mathbf{G}(\omega)$ is invertible, which is true here away from its eigenvalues, which are the natural frequencies of the clamped-free bar. Accordingly, the internal displacements of the bar can be condensed to the boundary $x = 1$ through the expression

$$ \hat{u}_N(1, \omega) = \phi(1)^{\top} \hat{u}(\omega) = \phi(1)^{\top} \mathbf{G}(\omega)^{-1} \phi(1) \hat{p}(\omega) = G_N(\omega) \hat{p}(\omega) \quad (10) $$

with $G_N(\omega) = \phi(1)^{\top}\mathbf{G}(\omega)^{-1}\phi(1)$, approximation of the exact ratio $G(\omega) = \hat{u}(1, \omega)/\hat{p}(\omega)$ induced by the space semi-discretization with the property $\lim_{N \to \infty} G_N(\omega) = G(\omega)$ in some sense to be specified. For a constant cross-section bar system, it is known that $G(\omega) = \tan \omega/\omega$ [15, Eq. (23)]<fnref n="3" />. However, there is an important difference in the construction of the functions $G(\omega)$ in [15] and $G_N(\omega)$ in the present paper. The former is established on the strong solution to the problem while the latter is established on the weak solution. This results in a distinct definition of the strain term, ie $u_x$ and $\hat{u}_x$, either strong or weak $x$-derivative of $u$ and $\hat{u}$, respectively, which has consequences on the convergence of the proposed scheme, as explained later.

<footnote n="1">¹The following operations are performed starting from the Equation (1): (1) a Fourier transform in time, (2) a semi-discretization in space and (3) a prescribed periodicity in time. It is straightforward to show that this is strictly equivalent to (1) a semi-discretization in space, (2) a prescribed periodicity in time and (3) a Fourier transform. These last two operations form what is commonly known as the Harmonic Balance Method, applied here on a system of Ordinary Differential Equations. This second approach is used in [26, Eq. (3.24)]. The above three operations commute and there is no preferred order.</footnote>

<footnote n="2">²Note that the participation $\hat{u}_N(\omega)$ of shape function $\phi_N(x)$ should not be confused with the discretized displacement $\hat{u}_N(x, \omega)$.</footnote>

<footnote n="3">³Note that the scalar $\mathbf{G}(\omega)$ should not be confused with the matrix $\mathbf{G}(\omega)$.</footnote>
### 3.1.1 Lagrange linear finite elements

Equation (4) is discretized into N classical Lagrange P1 finite elements (labelled as 'FE' in the sequel) with N + 1 nodes $x_i$, $i = 0, \dots, N$ and attendant hat shape functions

$$ \phi_i(x) = \frac{1}{h} \begin{cases} x - x_{i-1} & x \in [x_{i-1}, x_i], \\ x_{i+1} - x & x \in [x_i, x_{i+1}], \\ 0 & \text{otherwise.} \end{cases} \qquad (11) $$

where the mesh size is $h = 1/N$. Accordingly, $\hat{u}_i(\omega) \approx \hat{u}(ih, \omega)$, $i = 0, \dots, N$. Finding the expression of $G_N(\omega)$ in Equation (10) requires an inversion of the matrix $G(\omega)$ for every $\omega$, which might be computationally expensive. This can be accelerated since $G(\omega)$ has a symmetric tridiagonal structure of the form

$$ G(\omega) = \begin{bmatrix} d_1 & c_1 & & 0 \\ c_1 & \ddots & \ddots & \\ & \ddots & c_{N-1} & d_N \\ 0 & c_{N-1} & & d_N \end{bmatrix} \qquad (12) $$

for the considered finite element family, where $c_i$ and $d_i$, $i = 1, \dots, N$ are abbreviations for $c_i(\omega)$ and $d_i(\omega)$, respectively. A formal Gaussian elimination leads to $G_N(\omega) = 1/(d_N - c_{N-1}a_{N-1})$ with $a_1 = c_1/d_1$ and $a_i = c_i/(d_i - c_{i-1}a_{i-1})$ for $i = 2, \dots, N-1$.

### 3.1.2 Component mode synthesis: Craig-Bampton method

Component mode synthesis (CMS) techniques are widely employed to reduce the computational cost associated to the analysis of complex engineering structures. Among many incarnations, an attractive option for the space semi-discretization is the Craig-Bampton formulation [6], labelled as 'CB' in the remainder. To generate the CB model, both fixed-interface normal modes and constraint modes of the considered structure are calculated. For the problem at hand, the interface is located at $x = 1$ and the family consists in $N - 1$ natural modeshapes of the clamped-clamped bar complemented by one static mode, which can be seen as a lifting function, induced by a prescribed unit displacement at the interface. Collectively, for a constant cross-section bar, this reads

$$ \phi_i(x) = \sin(i\pi x), \quad i = 1, 2, \dots, N-1 \quad \text{and} \quad \phi_N(x) = x, \quad x \in [0, 1]. \qquad (13) $$

In case of a non-constant cross-section bar, the family is computed numerically instead, based on a FE discretization of the formulation using a sufficiently fine mesh.

### 3.1.3 Linear normal modes

The last considered discretization lies on a truncated family of natural modes of the clamped-free bar, referred to as 'LM' in this work. For a constant cross-sectional area, the normalized natural modes read

$$ \phi_i(x) = \sqrt{2} \sin(\omega_i x) \quad \text{with} \quad \omega_i = (2i - 1)\pi/2, \quad i = 1, 2, \dots, N, \quad x \in [0, 1]. \qquad (14) $$

Such shape functions are inserted in Equation (6) and Equation (8) is updated accordingly: $M = I$ and $K$ is a strictly diagonal matrix with entries $\omega_i^2$, $i = 1, 2, \dots, N$.

For a space-dependent cross-sections, the discretized modeshapes of the system are used instead, as exact solutions are only known for limited cases [20]. Such modes are computed from the above FE approximation with $\bar{N}$ finite elements by solving the classical eigenvalue problem $(K - \omega^2 M)u = 0$. Only the first $N \ll \bar{N}$ eigenpairs $(\omega_i, u_i)$ are kept in the (modal) discretization of the equations.

The expansion family suggested in Equation (14), or its numerical counterpart, seems awkward for the problem of interest, because it implies that the stress $A(1)\hat{u}_{N,x}(1, \omega)$ of the approximate solution vanishes. In other words, the family satisfies an homogeneous Neumann boundary condition at $x = 1$ and is thus not appropriate to sustain the unilateral contact force arising in the Signorini boundary condition. A classical approach is to introduce an additional static mode induced by a unit prescribed force at $x = 1$ into Equation (14), similar to the CB static mode. This is known as the Craig-Chang method [17, 22], labelled 'CC' below. For the considered bar system, it improves the accuracy of the approximation compared to LM since it can cope with a non-vanishing stress at $x = 1$ thus capturing the participation of the unilateral contact term. However, the diagonal structure of the reduced matrices is lost.

## 3.2 Time-domain discretization for periodic solutions

Given Equations (9) and (10), the only unknowns of the problem are (1) the nodal displacement $\hat{u}_N(\omega)$ and (2) the corresponding nodal contact force $\hat{p}(\omega)$. Since periodic solutions are of interest, it is natural to expand both of them
as truncated Fourier series with a common fundamental frequency $\Omega$, that is

$$p(t) \approx p_M(t) = \sum_{k=0}^{M} a_k \cos(k\Omega t) = \mathbf{F}(t)^{\top}\mathbf{a} \quad \text{and} \quad u_N(1,t) \approx u_{NM}(1,t) = \sum_{k=0}^{M} b_k \cos(k\Omega t) = \mathbf{F}(t)^{\top}\mathbf{b} \quad (15)$$

with the notation $\mathbf{F}(t) = [1, \cos(\Omega t), \dots, \cos(M\Omega t)]$. Only the cosine terms are considered in the Fourier series as suggested in [15]. The corresponding Fourier transforms read

$$\hat{p}_M(\omega) = \sum_{k=0}^{M} a_k \delta_{k\Omega} \quad \text{and} \quad \hat{u}_{NM}(\omega) = \sum_{k=0}^{M} b_k \delta_{k\Omega}. \qquad (16)$$

Accordingly, the real coefficients $(a_k, b_k)$, $k = 0, \dots, M$ are the unknowns of the problem now discretized in time, along with $\Omega$. Inserting Equation (16) into Equation (10) leads to

$$b_k = G_N(k\Omega)a_k, \quad k = 0, \dots, M \qquad (17)$$

which can also be read as a Harmonic Balance version of Equation (10). Other time-domain families such as wavelets [11, 12] could also be implemented in order to replace Fourier functions. However, Equation (17) would have to be updated accordingly with a more complicated expression.

## 3.3 Unilateral contact condition

Through Equation (17), the coefficients $b_k$ are expressed in terms of the coefficients $a_k$ which are thus the only remaining unknowns along with $\Omega$. As suggested in [15], they can be solved for via a purely numerical version of the Harmonic Balance Method performed on Equation (2) (unilateral contact conditions) and resulting in

$$g_k(\Omega, \mathbf{a}) = \frac{1}{T} \int_0^T \cos(k\Omega t) (\max[\alpha(u_{NM}(1,t) - g_0) + p_M(t), 0] - p_M(t)) dt = 0, \quad k = 0, \dots, M \quad (18)$$

where $T = 2\pi/\Omega$ is the period and where the expansions (15) are first inserted in the integrand. In the right hand-side of Equation (18), it should be noted that both $u_{NM}(1,t)$ and $p_M(t)$ solely depend on $\Omega$ and $\mathbf{a}$ through Equations (15) to (17), hence the notation in the left hand-side of the equation. In order to solve Equation (18) which bears $M + 2$ unknowns for $M + 1$ equations, sequential continuation is conducted with $\Omega$ chosen as the continuation parameter [15]. The system of nonlinear equations (18) in the coefficients $a_k$, $k = 0, \dots, M$, is solved numerically using the trust-region dogleg solver available in Matlab® [15, 21]. The integrals in Equation (18) are approximated numerically using classical quadrature schemes. The gradient of $g(\Omega, \mathbf{a})$ is also expressed analytically in a piecewise manner and evaluated numerically in the nonlinear solver. It reads

$$\partial_{\mathbf{a}} g_k(\Omega, \mathbf{a}) = \frac{1}{T} \int_{0}^{T} \cos(k\Omega t) (\partial_{\mathbf{a}} \max[\alpha(u_{NM}(1,t) - g_0) + p_M(t), 0] - \partial_{\mathbf{a}} p_M(t)) dt, \quad k = 0, 1, \dots, M \quad (19)$$

where the partial differentiation of the max operator is evaluated following Clarke's subdifferential rule [5]:

$$\partial_{\mathbf{a}} \max[\alpha(u_{NM}(1,t) - g_0) + p_M(t), 0] = \begin{cases} \alpha \partial_{\mathbf{a}} u_{NM}(1,t) + \partial_{\mathbf{a}} p_M(t) & \text{if } \alpha(u_{NM}(1,t) - g_0) + p_M(t) > 0 \\ 0 & \text{otherwise.} \end{cases} \quad (20)$$

The partial derivatives in Equations (19) and (20) have explicit expressions: $\partial_{\mathbf{a}} p_M(t) = \mathbf{F}(t)$ from Equation (15) and $\partial_{\mathbf{a}} u_{NM}(1,t) = \mathbf{F}_G(t)$ with $\mathbf{F}_G(t) = [G_N(0), \cos(\Omega t)G_N(\Omega), \dots, \cos(M\Omega t)G_N(M\Omega)]$ according to Equations (15) and (17). They are explicitly implemented in the nonlinear solver, which greatly reduces the computing cost with a gain of about 50% compared to a purely numerical procedure where the gradient of the nonlinear function is approximated via numerical differentiation.

After the space semi-discretization exposed in Section 3.1, the space-semi discretized counterparts of Equations (1) and (2), would normally be supplemented with an impact law to guarantee uniqueness. This aspect is intentionally ignored in the present work. Our conjecture is that by satisfying the Signorini conditions in a weighted residual sense only within a space of periodic functions in time, the purely elastic Newton's impact law with $e = 1$ is recovered as the number $N$ of shape functions tends to infinity. We are not aware of any mathematical proof of this statement. However, in addition to the remainder of the paper, Appendix A shows some numerical evidence for a one-degree-of-freedom academic system. It should also be noted that the $k$-impact-per-period solutions exposed in [13, 29], where an energy-preserving Newton impact law is explicitly implemented, could also be approximated by the proposed methodology for discrete $N$-degree-of-freedom spring-mass systems.

Solution strategies combining frequency and time domains for nonsmooth systems have already been proved successful, such as the Dynamic Lagrangian Frequency-Time (DLFT) methodology [18] or the hybrid Shooting and Harmonic Balance method for unilaterally constrained systems [25]. The DLFT was also shown to compare well with experimental investigations [4] and extended to turbomachinery applications with fretting wear [23, 24]. The
solution method developed above for discrete systems with unilateral contact condition can be extended to systems with friction since the latter can also be expressed in terms of nonsmooth equalities [9, 14, 27]. Even though the above approach and the DLFT differ in the way they handle the contact term in the governing equation and the solution algorithm, they share the common feature of generating solutions satisfying the unilateral constraint in an integral rather than a pointwise sense. Another difference is that the DLFT explicitly requires the AFT procedure [3] in order to navigate between the frequency domain where the equation of motion is solved and the time domain (over a single period) where the unilateral contact and friction conditions are enforced. In the proposed method, this operation is performed implicitly in the sense that the inverse (Fast) Fourier Transform of the AFT is equivalent to Equation (15) while the direct (Fast) Fourier Transform is equivalent to Equation (18). Finally, the AFT is tightly connected to Fourier functions, also chosen in the present work. However, in Equations (15) and (18), it would be straightforward to implement other periodic families such as periodized wavelets for instance [11].

# 4 Convergence analysis

Errors introduced in the approximated solution are twofold: space discretization through the number of shape functions $N$ and time/frequency discretization through the number of harmonics $M$. They are *not* independent and at least related via Expression (17) which acts as a necessary yet not sufficient condition to ensure an accurate approximation. Once a frequency $\Omega$ and a number of harmonics $M$ are specified, it is advised to accurately compute $G_N(k\Omega)$ for $k = 0, \dots, M$ to properly satisfy Equation (17). In other words, $N$ should increase with $M$.

## 4.1 Transfer function

The first part of the convergence analysis is on the “transfer function” $G(\omega)$, which captures the linear dynamics within the clamped-free bar system, and its (various) approximation(s) $G_N(\omega)$ defined in Equation (10). The frequency range on which $G(\omega)$ is investigated is $\omega \in ]0; 200[$. This is dictated by the fact that the frequency of the nonsmooth modes computed later lies near $\Omega = 5$, along with $M = 40$ so that the highest harmonic is around $\omega = M\Omega = 200$, $M$ being selected based on the full convergence analysis detailed in Section 4.2. Covering such frequency range in the convergence analysis ensures the accuracy of the results. However, this range can be reconsidered depending on the modal motion of interest. Also, as already said, convergence on the transfer function is a necessary yet insufficient ingredient for the convergence of the proposed approximation.

**Constant cross-section** First, a constant cross-section bar is considered, for which $G(\omega) = \tan \omega / \omega$. Figure 2 shows a comparison for various $N$ with FE discretization. The approximation with $N = 128$ shows a severely growing phase shift with $\omega$. However, the FE approximation with $N = 4096$ compares very well with the exact transfer function in the considered frequency interval.

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Constant cross-section bar transfer function: (a) values of transfer function: exact solution [---], FE with $N = 128$ [---] and $N = 4096$ [---]; (b) relative error of transfer function: FE with $N = 128$ [---] and $N = 4096$ [---].</figcaption></figure>

Figure 3 shows a comparison where the LM approximation embeds $N$ clamp-free linear modes computed numerically via FE with $\bar{N} = 20N$ finite elements. For $N = 50$, the increasing phase shift between the exact and approximate solutions remains limited in the considered frequency range. At frequencies higher than 150 rad/s, the approximation is no longer meaningful and the number of modes in the reduction should be increased to approximately $N = 400$ to recover accuracy.

Figure 4 provides a comparison with the CB approximation. The Craig-Bampton modes are computed numerically via FE with $\bar{N} = 20N$ finite elements. At $N = 50$, the approximate function $G_N(\omega)$ matches the exact transfer
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Constant cross-section bar transfer function: (a) values of transfer function: exact solution [---], LM for $N = 50$ [—] and $N = 400$ [—]. With $N = 50$, $|G_N(\omega)|$ goes out of bounds for $\omega > 155$; (b) relative error of transfer function: LM for $N = 50$ [—] and $N = 400$ [—].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Constant cross-section bar transfer function: (a) values of transfer function: exact solution [---], CB for $N = 50$ [—] and $N = 400$ [—]. With $N = 50$, $|G_N(\omega)|$ goes out of bounds for $\omega > 190$; (b) relative error of transfer function: CB for $N = 50$ [—] and $N = 400$ [—].</figcaption></figure>

function well in the range $\omega \in ]0; 110[ \text{ rad/s}$. The accuracy quickly reduces in the range $\omega \in ]110; 150[ \text{ rad/s}$. At frequency higher than 150 rad/s, the approximation diverges. Again, accuracy is recovered with $N = 400$ in the range $\omega \in ]0; 200[ \text{ rad/s}$.

The relative error

$$E_N(\omega) = \left| \frac{G_N(\omega) - G(\omega)}{G(\omega)} \right| \qquad (21)$$

is now introduced to better understand how the error depends on $N$ and $\omega$. In Figure 5(a), the relative error at $\omega = 200$ (the upper bound frequency in the considered interval where error is the largest) for increasing $N$ is shown for all three methods in addition to CC. Both LM and CB reduced-order models start to converge around $N = 100$ while the FE starts to converge around $N = 1000$. CB and FE converge at rate close to $O(1/N^2)$ whereas LM converges at a rate close to $O(1/N)$. The lower convergence rate of LM is caused by the homogeneous Neumann boundary condition used at $x = 1$ to construct the expansion family, as discussed above and investigated further in Appendix B. As expected and shown in Figure 5(a), CC converges faster, at a rate close to $O(1/N^2)$. However at $\omega = 200$, $N = 400$ CC's relative error is not significantly better than that of LM with $N = 400$.

The relative error versus $\omega$ is shown in Figure 5(b) for a sufficiently large $N$, again for each of the proposed methods, see Figures 2 to 4. All methods exhibit relative error spikes (vertical asymptotes) located at the clamped-free and clamped-clamped system's natural frequencies where $G(\omega)$ diverges or vanishes, respectively. The error for FE, CB and CC increases at the rate $O(\omega^3)$ whereas the error of LM increases at the rate $O(\omega^2)$ which is better.

The relative error in the constant cross-section case is bounded by $O(\omega^3/N^2)$ for FE, which matches known results [2, 10] and CB behaves in the same way. LM error is bounded by $O(\omega^2/N)$ which indicates a slower convergence versus $N$. CC improves over LM and is bounded by $O(\omega^3/N^2)$. Overall, FE and CB converge the
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Relative error on the transfer function for a constant cross-section bar: (a) in terms of N at $\omega = 200$: FE [—], LM [—], CC [—] and CB [—]; (b) in terms of $\omega$: FE with $N = 4096$ [—], LM with $N = 400$ [—], CC with $N = 400$ [—] and CB with $N = 400$ [—]. Vertical asymptotes exist on every max and min (green curve) value but are not shown for clarity purposes.</figcaption></figure>

fastest but CB exhibits a lower error at specified N and is thus considered the best option in the estimation of the transfer function.

**Non-constant cross-section** A non-constant cross-section case is also explored with $A(x) = 1.5 - x$. The exact expression for $G(\omega)$ is no longer available. Instead, it is approximated and defined as the “true” solution via FE with $N = 25600$. The relative error is provided in Figure 6. The transfer function plots are skipped in this paragraph since they are very similar to that of the constant cross-section case, see Figures 2 to 4. Overall, the convergence is similar to that of the constant cross-section case. CB is still the optimal method. CC shares the same convergence rate as CB but has lower accuracy at a specified N while its computational cost is almost the same as CB. It is no longer considered in the remainder of this work.

<figure><img src="image_8.png" /><figcaption>**Figure 6:** Relative error on the transfer function for the non-constant cross-section bar: (a) in terms of N at $\omega = 200$: FE [—], LM [—], and CB [—]; (b) in terms of $\omega$: FE with $N = 4096$ [—], LM with $N = 400$ [—] and CB with $N = 400$ [—].</figcaption></figure>

## 4.2 Unilateral contact residual

The second convergence analysis focuses on the nonlinear Signorini condition introduced in Equation (2). A periodic solution at a specified frequency is targeted for each configuration. Convergence is explored via the $L^2$ error norm

$$R_{NM} = \int_{0}^{T} \left( \max[\alpha(u_{NM}(1, t) - g_{0}) + p_{M}(t), 0] - p_{M}(t) \right)^{2} dt \quad (22)$$

which depends on both N and M. A similar convergence analysis in terms of M only is provided in [15].

**Constant cross-section** This section considers the initial gap $g_0 = 0.001$. The convergence analysis is conducted on a periodic solution located on the branch of the first nonsmooth mode in Figure 16 at $\Omega = 1.788 \approx 1.138\omega_1$. Also, in the sequel, the i th natural frequency of the linear clamped-free bar is labelled $\omega_i$. Computed boundary displacement and contact force are shown in Figure 7. The displacement plots via all three methods are similar, so do the contact force plots. The corresponding error norm for each method is shown in Figure 8. As explained in [15], the error is affected by Gibbs’ phenomenon which arises along time in the contact force (and even in space for LM), and could theoretically be reduced, in $L^2$ sense, by increasing M. For a large specified N, the approximate
<figure><img src="image_9.png" /><figcaption>**Figure 7:** Boundary displacement $u_{NM}(1,t)$ [—] and contact force $p_M(t)$ [—] along with the exact contact force solution $p(t)$ [—] for a periodic solution with active contact at $\Omega = 1.138\omega_1$, $M = 40$. Constant cross-section. (a) FE with $N = 4096$, (b) LM with $N = 400$ LM model and (c) CB with $N = 400$.</figcaption></figure>

transfer function $G_N(\omega)$ is sufficiently accurate and the error norm reduces with respect to $M$. For a specified $M$, the error norm reduces only for a sufficiently large $N$, above which the transfer function is well approximated, but reaches a plateau, mainly because Gibbs' phenomenon is dominant and dictated by $M$. Overall, convergence is achieved with the convergence rate of $O(1/M)$ for a sufficiently high $N$.

<figure><img src="image_10.png" /><figcaption>**Figure 8:** Error norm $R_{NM}$. Constant cross-section bar at $\Omega = 1.138\omega_1$. (a) FE, (b) LM and (c) CB.</figcaption></figure>

**Non-constant cross-section** This section considers the non-constant cross-section $A(x) = 1.5 - x$. Again, one specific solution located on the branch of the first nonsmooth mode in Figure 18 at $\Omega = 1.102\omega_1$ is investigated. The displacement and contact force at $x = 1$ are shown in Figure 9. All three methods converge to the same solution. The overall convergence behavior shown in Figure 10 is similar to that of the constant cross-section case. However,

<figure><img src="image_11.png" /><figcaption>**Figure 9:** Boundary displacement $u_{NM}(1,t)$ [—] and contact force $p_M(t)$ [—] for a periodic solution with active contact at $\Omega = 1.102\omega_1$, $M = 40$. Non-constant cross-section $A(x) = 1.5 - x$. (a) FE with $N = 4096$, (b) LM with $N = 400$ and (c) CB with $N = 400$.</figcaption></figure>

it is worth noting that the error reduces faster than for the constant cross-section counterpart. The convergence rate reaches approximately $O(1/M^2)$ for a sufficiently high $N$. This can be explained by the difference in the “nature” of the contact force between the two solutions. For the constant cross-section configuration, the exact contact force is known to be a step function, with attendant numerical difficulties induced by Gibbs’ phenomenon. For the non-constant cross-section configuration, the contact force is instead a continuous function of time with limited Gibbs’ effect. However, this is only for the chosen solution and various others are compared in Section 5.2.

Overall, the suggested methodology exhibits appropriate convergence properties for both constant and non-constant cross-sections along with the proposed space discretizations.

## 4.3 Unilateral contact condition

In this section, an accurate $G_N$ is assumed. The question to be addressed is on how the relationship between the approximated contact force $p_M(t)$ and the approximated stress $A(1)u_{NM,x}(1,t)$ is enforced. The two quantities are
<figure><img src="image_12.png" /><figcaption>**Figure 10:** Convergence of error norm $R_{NM}$. Non-constant cross-section system with $A(x) = 1.5 - x$ at $\Omega = 1.102\omega_1$. (a) FE, (b) LM and (c) CB.</figcaption></figure>

theoretically equal to each other, see in Equation (5) used in the weak formulation. However, the contact force is an unknown of the problem while the stress term is retrieved from the displacement field through space differentiation of the shape functions. The various space semi-discretizations have implications on how these quantities are actually computed.

Constant cross-section and gap $g_0 = 0.001$ are considered. The convergence analysis is conducted on a solution located on the branch of the first nonsmooth mode at $\Omega = 1.138\omega_1$. In Figure 11, stress fields are displayed and the focus is on comparing the solutions at $x = 1$. All plots exhibit oscillatory overshoot dominated by Gibbs’

<figure><img src="image_13.png" /><figcaption>**Figure 11:** Stress field of a periodic solution with unilateral contact at $\Omega = 1.788$. Constant cross-section. (a) FE with $N = 4096$, (b) LM with $N = 400$ and (c) CB with $N = 400$. Axis labels detailed in [15, Fig 3(a)].</figcaption></figure>

phenomenon initiated at the contact interface by the Fourier approximation of the step-like contact force which then propagates to the whole domain by wave propagation mechanisms. The approximation of the spectrum of the wave equation differential operator might play a role as well but the differences between the plots being minimal, it does not seem to be critical here. Gibbs’ pollution in the strain field can be reduced via filtering techniques for instance [8]. This aspect is not discussed in the present work. Also, in contrast to the WFEM-based solution [34], the same Gibbs’ phenomenon arises in the stress fields computed using the methodology exposed in [15] (again, with a Fourier discretization in time only and not in space) even though they were not shown in the paper. FE and CB with their respective discretization show almost no difference between $p_M(t)$ and $EAu_{NM,x}(1, t)$. However, with LM, the two quantities can be clearly distinguished. Such a difference is due to the nature of the clamp-free modes used in the LM discretization: each expansion mode satisfies $\phi_{i,x}(1) = 0$ such that $u_{NM,x}(1, t) = 0$ holds. However, in the proposed plot, each spatial derivative is estimated numerically via finite difference with the consequence that the approximated strain is approximately 0 only, as shown in Figure 12(b). Clarifying information is provided in Appendix B. Away from the contact interface, for $x < 1$, the stress is approximated properly if we ignore the discussed pollution effect.

<figure><img src="image_14.png" /><figcaption>**Figure 12:** Boundary displacement $u_{NM}(1, t)$ [---] and contact force $p_M(t)$ [---] along with approximated contact force $A(1)u_{NM,x}(1, t)$ [---] for a periodic solution with active contact at $\Omega = 1.138\omega_1$, $M = 40$. Constant cross-section. (a) FE with $N = 4096$, (b) LM with $N = 400$ and (c) CB with $N = 400$.</figcaption></figure>
Convergence analysis is conducted through the L² error norm

$$R_N = \frac{1}{T} \int_0^T (p_M(t) - A(1)u_{NM,x}(1, t))^2 dt \quad (23)$$

which depends on N and M and measures how well Equation (5) is enforced. Since convergence with respect to N is of interest, the number of harmonics is set to M = 40. The convergence plot is shown in Figure 13(a). The error of both FE and CB behaves in O(1/N) with CB error below. LM does not converge due to the considered shape functions. The L²-norm plot in Figure 13(b) for the non-constant cross-section case contact force convergence analysis also leads to the same conclusion.

<figure><img src="image_15.png" /><figcaption>Figure 13: Contact force error $R_N$ defined in Equation (23) for a NSM solution: FE [—], LM [—] and CB [—]. (a) Constant cross-section, M = 40, Ω = 1.138ω₁, (b) non-constant cross-section A(x) = 1.5 - x, M = 40, Ω = 1.102ω₁.</figcaption></figure>

# 5 Nonsmooth modal analysis

In this section, the above solution strategy is used to perform nonsmooth modal analysis on a few bar systems mainly to illustrate its capabilities. The classical sequential continuation technique [19] is implemented where Ω is successively increased by a small given increment on a given interval of interest. Accordingly, the nonlinear system is solved for the Fourier coefficients only.

Attendant subharmonics of the natural frequencies are also labelled to help identify potential nonlinear modal behaviors. The backbone curve is normalized with respect to the energy of the first linear grazing mode of the system. Finally, $g_0 = 0.001$ and M = 40 are considered in all simulations.

## 5.1 Benchmark case: constant cross-section

A constant cross-section is considered as a benchmark. This configuration was investigated in previous works [15, 16, 30–34] using a different solution strategy. The FE discretizations N = 4096 and N = 128 are used to capture high and low resolution solutions, and compared to FD-BEM/HBM results with the exact transfer function G(ω) [15]. The purpose of NSA of the benchmark system is to show how the accuracy of G_N(ω) at high frequency affects the NSA.

Figure 14 shows the three shown backbone curves computed via sequential continuation using FE in the current investigation. For a sufficiently large N, the results in Figure 14(a) agree well with the FD-BEM/HBM results

<figure><img src="image_16.png" /><figcaption>Figure 14: First NSM backbone curve. Constant cross-section. (a) FE with N = 4096, (b) FE with N = 128 and (c) FD-BEM/HBM.</figcaption></figure>

exposed in Figure 14(c). However, the approximation with N = 128 cannot accurately follow the main backbone curve, as shown in Figure 14(b) which features a cloud of points. This phenomenon is induced by the error in G_N with a low N: high-frequency content in the solution is missed. Such statement is further supported by the plots of the displacement and contact force in Figure 15. It should be noted that a number of internal resonances are known to exist along this main backbone curve and are numerically challenging to track.
<figure><img src="image_17.png" /><figcaption>**Figure 15:** Boundary displacement $u_{NM}(1, t)$ [—] and strain $p_M(t)$ [—] along with Signorini residual $r$ [—] for a periodic solution with active contact at $\Omega = 1.146\omega_1$ (a) FE with $N = 4096$, (b) FE with $N = 128$ and (c) FD-BEM.</figcaption></figure>

The backbone curves computed via FE, LM and CB are compared in Figure 16. In the vicinity of $\omega_1$, all three methods were unable to capture the grazing-like solutions similar in nature to the linear modes. In this region, the nonlinear solver tends to converge to the trivial solution because the contact duration is short which would require higher harmonics in the approximation to be well captured. Various versions of a nonlinear grazing solution in the vicinity of $\omega_1$ are plotted in Figure 17. Examples of such a grazing mode can also be found in [34]. Around $4\omega_1/3$, all three methods perform similarly. FE and CB still exhibit thin clouds of points just above the main backbone curve as shown in Figure 16(a) and (c). They actually represent many short vertical backbone branches induced by internal resonances in the system. LM instead captures a few such vertical branches, on a larger energy range however, as shown in Figure 16(b). All three methods with appropriate choice of $N$ and $M$ can correctly predict

<figure><img src="image_18.png" /><figcaption>**Figure 16:** Backbone curve of first NSM. Constant cross-section bar: (a) FE with $N = 4096$, (b) LM with $N = 400$ and (c) CB with $N = 400$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17:** Boundary displacement $u_{NM}(1, t)$ [—] and strain $p_M(t)$ [—] along with Signorini residual $r$ [—] for the grazing-like solutions in the vicinity of $\omega_1$. Constant cross-section bar. In Figure 16: (a) FE with $N = 4096$ at ‘a’, (b) LM with $N = 400$ at ‘b’ land (c) CB with $N = 400$ at ‘c’.</figcaption></figure>

the topology of the low-energy first backbone curve, and some of the subsequent ones as shown in the remainder, even though it is clear that FE and even more CB perform better. The latter is thus selected in the nonsmooth modal analysis of the non-constant cross-section bar detailed in Section 5.2.

## 5.2 Non-constant cross-section

In this section, the modal response of the system with $A(x) = 1.5 - x$ is explored via CB with $N = 400$. It should be noted that the dynamics observed in the bar responses is very rich and a thorough examination is out-of-scope of this paper. Instead, this part of the work plays the role of a proof-of-concept of the developed methodology. Accordingly, the analysis is focused on the similarities shared with the constant cross-section system. To provide an idea of the computational cost associated to these calculations, a personal laptop with a 3.8GHz CPU and 12 threads was used and the computation of a single backbone curve with 1000 points required about one minute.

In Figure 18, the first [green] and second [red] nonsmooth modes along with subharmonic 2 of the second mode [orange], in the vicinity of $\omega_2/2$, are shown. In contrast to the constant cross-section bar, a major discontinuity in the first backbone curve arises between points labelled ‘a’ and ‘b’. The corresponding boundary displacement,
<figure><img src="image_20.png" /><figcaption>**Figure 18:** Backbone curves of the non-constant cross-section case via CB, *N* = 400. First mode [●], second mode [●], and subharmonic 2 of the second mode [●].</figcaption></figure>

contact force and Signorini residual are compared in Figure 19(a) and (b). The displacement fields within the bar are indicated in Figure 20(a) and (b) which also shows the displacement field of the second mode along with subharmonic 2 of the second mode.

<figure><img src="image_21.png" /><figcaption>**Figure 19:** Boundary displacement $u_{NM}(1, t)$ [—] and contact force $p_M(t)$ [—] along with Signorini residual $r$ [—] of different modes: (a) Mode ‘a’ in Figure 18, (b) Mode ‘b’ in Figure 18.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20:** Displacement field of the second mode: in Figure 18, (a) Mode ‘a’, (b) Mode ‘b’, (c) Mode ‘c’, (d) Mode ‘d’, (e) Mode ‘e’, (f) Mode ‘f’. Axis labels detailed in [15, Fig 3(a)].</figcaption></figure>

The backbone curve of a constant cross-section bar system is almost continuous, as shown in Figure 17. However,
the backbone curve for the varying cross-section system shown in Figure 18 exhibits a jump between point ‘a’ and point ‘b’. The corresponding modal motions feature distinct patterns in displacement and contact force. Since this phenomenon is not observed in the constant cross-Section case, we decided to project the solution found at points ‘a’ and ‘b’ on the first linear clamped-free modes of the system. This projection, explained in [32], shows the participation of each linear mode in the nonsmooth mode (for a given frequency and energy). The projection on each linear mode is defined to be

$$a_k(t) = \frac{\int_0^L U_k(x)u_{NM}(x,t) \, dx}{\int_0^L U_k(x)U_k(x) \, dx} \qquad (24)$$

where $U_k(x), k \in 1, 2, \dots$ are the linear modes of the clamped-free system. The first seven linear modes are considered in this work. The average participation $P_k$ of the k-th linear mode in the NSM over one period T of the nonlinear motion is defined as

$$P_k = \sqrt{\frac{1}{T} \int_0^T a_k^2(t) \, dt} \qquad (25)$$

Another metric is the computation of the harmonics $\bar{a}_{k\ell}$ of the periodic functions $a_k(t)$, such that

$$a_k(t) = \sum_{\ell=0}^{\infty} \hat{a}_{k\ell} \cos(\ell\Omega t). \qquad (26)$$

Figure 21(a) shows the various modal participations $P_k, k = 1, \dots, 7$, for the modal motions on the first backbone curve. The first linear mode clearly dominates in the participation of the fundamental harmonic in the nonlinear solution. This indicates that all solutions on the so-called first backbone curve actually share the same spatial features: they have no nodal point along $x$ with zero displacement, see Figures 20(a) and (b) as opposed to Figures 20(c) and (d) (around $x = 1/2$) and even (e) and (f). Also, the following values provide an idea of the

<figure><img src="image_23.png" /><figcaption>**Figure 21:** Average linear mode participations along the first backbone curve: mode 1 [•], mode 2 [•], mode 3 [•], mode 4 [•], mode 5 [•], mode 6 [•], mode 7 [•]. (a) Participations $P_k$ defined in Equation (25) of each linear mode $k, k = 1, \dots, 7$. Fourier harmonics $\hat{a}_{k\ell}$ of the corresponding LNM projections defined in Equation (26): (b) Mode ‘a’ in Figure 18 and (c) Mode ‘b’ in Figure 18.</figcaption></figure>

frequencies involved in the dynamics: $\omega_1 \approx 1.92$, $\omega_3/4 \approx 1.98$, $\omega \approx 2.07$ and $\omega_b \approx 2.14$. No subharmonics other than $\omega_3/4$ were found in the vicinity of the sudden change in the response, between $\omega_a$ and $\omega_b$. It is also clear that for this system, $\omega_1 \approx \omega_3/4$. Figures 21(b) and (c) display the Fourier harmonics in every linear mode projection, at the points ‘a’ and ‘b’ displayed in Figure 18, normalized with respect to the maximum amplitude. By comparing Figures 21(b) and (c), a sudden drop of the participation of linear mode 3/harmonic 4 (LNM index = 3, harmonic index = 4) and an increase of linear mode 2/harmonic 2 (LNM index = 2, harmonic index = 2) is observed, which agrees with Figure 21(a). However, a detailed conclusion on the mechanism in place is too early but an internal resonance between modes 1 and 3 is suspected.

From the above brief modal analysis we can say with confidence that CB is successful in finding NSM solutions for a non-constant cross-section bar system which requires a space semi-discretization in the solution procedure. Nonlinear phenomenons including frequency-energy dependency, internal resonance solution and subharmonic solution are well captured even though more detailed investigations would be required on specific frequency intervals depending on the need of the analyst.

# 6 Conclusion

In this paper, a solution methodology relying on a space semi-discretization scheme combined to the Harmonic Balance Method (HBM) is introduced to perform Nonsmooth Modal Analysis of one-dimensional bar systems satisfying a Signorini condition at one end. The non-constant cross-section is introduced here to motivate the need of a space semi-discretization which might rely on any acceptable linear combination of well-chosen shape functions,
as classically done in the Finite Element method or the Component Mode Synthesis techniques (Craig-Bampton, Craig-Chang and alike).

A systematic convergence analysis shows that the investigated methodologies exhibit a rate of convergence $O(1/M)$, where $M$ is the number of harmonics, as soon as the space semi-discretization is sufficiently fine. By sufficiently fine, we mean that the integer $N(M)$, function of $M$, defining the number of shape functions in space is sufficiently large. In other words, the space discretization is dictated by the time discretization (the number of harmonics here). A fine time discretization with a coarse space discretization leads to a divergence whereas a coarse time discretization with a fine space discretization generates an 'acceptable' approximation. It should be also noted that in the proposed methodology, the Signorini condition mathematically reflecting a unilateral contact condition is satisfied in a weighted-residual sense only within the HBM. As a notable consequence is that Newton's law, commonly introduced to retrieve uniqueness for such discrete mechanical systems, is here implicitly satisfied in a distributional sense, for a restitution coefficient preserving energy, without being explicitly implemented in the solution algorithm.

A brief nonsmooth modal analysis of the non-constant cross-section bar shows that the nonlinear spectrum can be approximated through the construction of the well-known frequency-energy backbone curves. The extension of the suggested solution technique to systems involving elaborate contact interfaces not limited to a unique point is certainly possible at the cost of more CPU-expensive numerical simulations.

# A Periodic solutions of one-dof systems

This appendix briefly discusses the proposed weak enforcement of the Signorini boundary condition for two academic one-dof systems featuring exact solutions in the sense of distributions.

## A.1 Forced one-dof system: the bouncing ball

A rigid mass $m$ subject to gravity, illustrated in Figure A.1(a), is considered. The (partial) governing equations are:

$$m\ddot{u}(t) = -mg + p(t), \quad u(t) \ge 0, \quad p(t) \ge 0, \quad u(t)p(t) = 0, \quad u(0) = u_0, \quad v(0) = 0 \quad (A.1)$$

where $u(t)$ is the displacement of the ball; $g$ stands for gravity and $p(t)$ is the contact force. The unilateral contact condition is supplemented with an impact law reading: for all $t$ such that $u(t) = 0$ with $\dot{u}(t-) < 0$, then $\dot{u}(t+) = -e\dot{u}(t-)$ to guarantee solution uniqueness [7]. The restitution coefficient is set to $e = 1$ to preserve energy and thus existence of periodic solution. Normalization $\bar{u} = u/2u_0$, $\bar{t} = t/\sqrt{2u_0/g}$ and $\bar{p} = p/mg$ yields $\ddot{u}(t) = -g + p(t)$, $u(0) = 0.5$, $v(0) = 0$ and $g = 1$. The exact solution, shown in Figure A.2(a), is periodic with period $T = 2$ and the solution in the range of $t \in [0, T]$ reads

<figure><img src="image_24.png" /><figcaption>Figure A.1: System of interest: bouncing ball.</figcaption></figure>

$$u(t) = \begin{cases} -0.5t^2 + 0.5 & 0 \le t \le 1 \\ -0.5t^2 + 2t - 1.5 & 1 < t \le 2 \end{cases} \quad (A.2)$$

along with $p(t) = \delta(t - 1)$ where $\delta(\cdot)$ is the Dirac delta distribution. The corresponding Fourier series are

$$u(t) = \frac{1}{3} - \frac{2}{\pi^2} \sum_{m=1}^{\infty} \frac{(-1)^m}{m^2} \cos(m\pi t) \quad \text{and} \quad p(t) = 1 + 2 \sum_{m=1}^{\infty} (-1)^m \cos(m\pi t) \quad (A.3)$$

where $p(t)$ is also the time derivative of $u(t)$ twice, termwise. Pointwise convergence is achieved for the displacement and distributional convergence for the contact force even away from the delta.

The methodology detailed in Sections 3.2 and 3.3 is now used with

$$u(t) \approx u_M(t) = \sum_{m=0}^{M} \hat{u}_m \cos(m\Omega t), \quad p(t) \approx p_M(t) = \sum_{m=0}^{M} \hat{p}_m \cos(m\Omega t). \quad (A.4)$$

The unknowns of the problem are the $2M + 3$ quantities $\{\hat{u}_m, \hat{p}_m, \Omega\}, m = 0, \dots, M$. A Fourier transform on Equation (A.1) leads to the $M + 1$ equations $\hat{p}_m = -(m\Omega)^2\hat{u}_m, m = 1, \dots, M$ along with $\hat{p}_0 = g$. The weak enforcement of the Signorini condition as developed in the paper leads to the $M + 1$ additional equations

$$g_m(\hat{u}, \hat{p}, \Omega) \equiv \frac{1}{T} \int_0^T \cos(m\Omega t)(-p_M(t) + \max[-\rho u_M(t) + p_M(t), 0])dt = 0, \quad m = 0, 1, \dots, M, \quad (A.5)$$

where $\hat{u} = (\hat{u}_0, \hat{u}_1, \dots, \hat{u}_M)$ and $\hat{p} = (\hat{p}_0, \hat{p}_1, \dots, \hat{p}_M)$. The last equation is $\sum_{m=0}^{M} \hat{u}_m = u_0$ since the initial displacement of the bouncing ball is prescribed and the frequency of the response is unknown. The initial velocity $v_0 = 0$ is implicitly enforced since Fourier series are limited to cos terms only. The above equations are collectively
solved using a nonlinear solver. In practice {$\Omega, \hat{u}$} are the sole independent unknowns and Equation (A.5) only should now be solved. We set $\rho = 10$.

The resulting numerically approximated contact force and displacement are shown in Figure A.2. For a low $M \approx 10$, the approximation exhibits a slightly longer period compared to the exact solution and the instantaneous impact is not accurately captured. For a higher $M \approx 40$, the approximate contact force converges to Dirac delta function (in a distributional sense). Convergence on the displacement and period is also achieved. Frequency-domain convergence analysis is illustrated in Figure A.3. It is self-explanatory and it is concluded that the proposed solution strategy to enforce a Signorini boundary condition via HBM is able to reproduces, in a weak sense, an energy-preserving impact law.

<figure><img src="image_25.png" /><figcaption>**Figure A.2:** Time-domain system response: (a) displacement, (b) contact force. Exact solution [—] and approximation with $M = 10$ [—], 20 [—] and 40 [—].</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure A.3:** Frequency-domain system response: (a) displacement, (b) contact force. Exact solution [—] and approximation with $M = 10$ [—], 20 [—] and 40 [—].</figcaption></figure>

## A.2 Autonomous one-dof system

Another simple spring-mass example, illustrated in Figure A.4, closer in spirit to the one considered in the present paper is briefly explored. A mass $m$ is attached to a spring $k$ and subject to a unilateral condition with initial gap $g_0$. Normalized quantities are: time $\bar{t} = \sqrt{k/m} t$, displacement $\bar{u} = u/g_0$ and impact force $\bar{p} = p/(g_0 k)$. The corresponding normalized governing equation reads $u(t) + \ddot{u}(t) = p(t)$ where the upper bar notation is dropped. Unilateral contact is governed by the complementarity Signorini condition along with the Newton impact law with $e = 1$. The main difference with the previous system is the absence of external forcing. Accordingly, the exact solution feature a continuum of (nonlinear) natural frequency $\Omega \in ]1, 2[$ and reads

<figure><img src="image_27.png" /><figcaption>**Figure A.4:** 1-dof spring-mass system with unilateral contact.</figcaption></figure>

$$u(t) = \frac{1}{\cos(T/2)} \begin{cases} \cos(t) & 0 \le n \le T/2 \\ \cos(t-T) & T/2 < n \le T \end{cases} \qquad (A.6)$$

with $T = 2\pi/\Omega$ and $p(t) = 2\tan(T/2)\delta(T/2)$. The exact solution for $T = 4.5$ is shown in Figure A.5. The corresponding Fourier series read

$$u(t) = \frac{2\tan(T/2)}{T} \left( 1 + \sum_{m=1}^{\infty} \frac{2(-1)^m}{1-m^2\Omega^2} \cos(m\Omega t) \right) \quad \text{and} \quad p(t) = \frac{2\tan(T/2)}{T} \left( 1 + \sum_{m=1}^{\infty} 2(-1)^m \cos(m\Omega t) \right). \quad (A.7)$$

The governing equations are again solved using the proposed solution strategy with $T = 4.5$. A total of $2M + 2$ unknowns {$\hat{u}, \hat{p}$} are to be found by solving the $M + 1$ equations $((1 - (m\Omega)^2)\hat{u}_m = \hat{p}_m, m = 0, 1, 2, \dots, M$ and the $M + 1$ additional Signorini conditions in the form (A.5) where $\Omega$, and thus $T$, is specified.

The displacement and contact force are shown in Figure A.5 for various $M$ (assuming a discretization of the type (A.4)). Pointwise convergence is achieved for the displacement and distributional convergence for the contact force with a rate slightly slower than $O(1/M)$ (not shown). Frequency-domain convergence analysis is illustrated in Figure A.6. Again, the proposed scheme finds approximate periodic solutions satisfying the impact law $e = 1$ in a weak sense, in agreement with results exposed in [26].
<figure><img src="image_28.png" /><figcaption>**Figure A.5:** Time-domain system response: (a) displacement, (b) contact force. Exact solution [—] and approximation with $M = 10$ [—], 20 [—], and 40 [—].</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure A.6:** Frequency-domain system response: (a) displacement, (b) contact force. Exact solution [—] and approximation with $M = 10$ [—], 20 [—], and 40 [—].</figcaption></figure>

## B Linear mode expansion convergence analysis

In this paper, the LM method was developed by considering a family of clamped-free (homogeneous Dirichlet-homogeneous Neumann) modeshapes. It is here shown how the method can generate an approximate solution to a bar problem with homogeneous Dirichlet-non-homogeneous Neumann boundary condition. To this end, we reduce the Helmholtz equation to a time-independent configuration. It becomes $u_{,xx}(x) = 0, x \in [0, 1]$, with $u(0) = 0$ and $u_x(1) = 1$ (non-homogeneous Neumann BC). The exact solution is $u(x) = x, x \in [0, 1]$. The clamp-free modes of the bar are listed in Equation (14). Using them to expand the approximation leads to the weak form

$$ \int_{0}^{1} (\phi_{x}(x)^{\top} \phi_{x}(x)) dx \mathbf{u} = \phi(1) \qquad (B.1) $$

with<fnref n="4" /> $\mathbf{u} \equiv [a_1, \dots, a_N]^\top$ storing the contributions of the shape functions (ie, the linear modes) collected in $\phi(x) = [\phi_1(x), \dots, \phi_N(x)]^\top$. The system of linear equations simplifies to $\omega_i^2 a_i = (-1)^{i-1} \sqrt{2}, i = 1, \dots, N$ so that $a_i = (-1)^{i-1} \frac{\sqrt{2}}{\omega_i^2}$ and the approximated displacement solution of the system reads

$$ u_N(x) = \sum_{i=1}^{N} (-1)^{i-1} \frac{2}{\omega_i^2} \sin(\omega_i x), \quad \omega_i = \frac{2i-1}{2}\pi, \qquad (B.2) $$

which converges pointwise to the exact solution $u(x) = x$ as $N \to \infty$.

The approximate transfer function defined in the paper reduces here to $G_N = u_N(1)$ since it is assumed that $p = 1$. Accordingly, the convergence of $u_N(1)$ as $N \to \infty$ is clarified. Equation (B.2) yields

$$ \lim_{N \to \infty} u_N(1) = \frac{8}{\pi^2} \sum_{i=1}^{\infty} \frac{1}{(2i-1)^2} = 1 = u(1). \qquad (B.3) $$

The error caused by truncation at order $N$ is

$$ E_N = |u(1) - u_N(1)| = \frac{8}{\pi^2} \sum_{i=N+1}^{\infty} \frac{1}{(2i-1)^2} \qquad (B.4) $$

with the following bounds:

$$ \int_{N+1}^{\infty} \frac{1}{(2x-1)^2} dx \le \sum_{i=N+1}^{\infty} \frac{1}{(2i-1)^2} \le \int_{N}^{\infty} \frac{1}{(2x-1)^2} dx. \qquad (B.5) $$

Further evaluation of the upper and lower bounds leads to

$$ \frac{8}{\pi^2} \frac{1}{4N+2} < E_N < \frac{8}{\pi^2} \frac{1}{4N-2} \qquad (B.6) $$

<footnote n="4">⁴The notation $a_i, i = 1, \dots, N$ is used instead of $u_i$ to avoid any possible confusion with the term $u_N(x)$ in Equation (B.2). Also the terms $a_i$ in this section should not be confused with coefficients $a_k$ in Section 3.2.</footnote>
indicating that the sequence {$u_N(1)$} for $N$ sufficiently large has the rate of convergence $O(1/N)$ and so has $G_N$, which matches the derivations reported in Section 4.1.

Although $u_N(x)$ converges pointwise to the exact solution, see Figure B.7(a), the approximated strain

$$u_{N,x}(x) = 2 \sum_{i=1}^{N} \frac{\cos(\omega_i x)}{\omega_i}, \quad \omega_i = \frac{2i-1}{2} \pi \qquad (B.7)$$

only converges in $L^2$ sense to $u_x$. Obviously, the approximated strain vanishes at $x = 1$, that is $u_{N,x}(1) = 0$, which is independent of $N$, and does not converge to the exact strain $u_x(1) = 1$. A Gibbs phenomenon is thus generated (here in space) around $x = 1$. The conclusion is the following: in the approximate solution, the LM method generates a Gibbs phenomenon in space which adds to the Gibbs phenomenon in time induced by HBM. The associated expansion family is clearly not optimal for the problem at hand but convergence is still achieved in a weak sense.

<figure><img src="image_30.png" /><figcaption>**Figure B.7:** Convergence analysis: (a) displacement, (b) strain. Exact solution [—] and LM solution with $N = 5$ [—] and $N = 20$ [—].</figcaption></figure>

**Acknowledgments** The authors gratefully acknowledge the financial support by the Natural Sciences and Engineering Research Council of Canada through the Discovery Grants program (421542-2018).

**Conflict of Interest** The authors declare that they have no conflict of interest.

**Data availability** The Matlab scripts used to generate the results exposed in this contribution are available at the permalink 10.5281/zenodo.7577381.

## References

[1] Acary, V. and Brogliato, B. *Numerical Methods for Nonsmooth Dynamical Systems. Applications in Mechanics and Electronics*. 2008. [DOI].

[2] Bayliss, A., Goldstein, C. I., and Turkel, E. "On accuracy conditions for the numerical computation of waves". 59 (1985), pp. 396–404. [DOI].

[3] Cameron, T. M. and Griffin, J. "An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems". *Journal of Applied Mechanics* (1989). [DOI], [OAI] a.

[4] Charleux, D., Gibert, C., Thouverez, F., and Dupeux, J. "Numerical and Experimental Study of Friction Damping Blade Attachments of Rotating Bladed Disks". *International Journal of Rotating Machinery* 2006 (2006), pp. 1–13. [DOI].

[5] Clarke, F. H. "Generalized gradients and applications". *Transactions of the American Mathematical Society* 205 (1975), pp. 247–262. [DOI].

[6] Craig Jr, R. R. and Bampton, M. C. "Coupling of substructures for dynamic analyses". *AIAA journal* 6.7 (1968), pp. 1313–1319. [DOI], [OAI] a.

[7] Doyen, D., Ern, A., and Piperno, S. "Time-Integration Schemes for the Finite Element Dynamic Signorini Problem". *SIAM Journal on Scientific Computing* 33.1 (2011), pp. 223–249. [DOI], [OAI] a.

[8] Duchon, C. E. "Lanczos filtering in one and two dimensions". *Journal of Applied Meteorology and Climatology* 18.8 (1979), pp. 1016–1022. [DOI].

[9] Hüeber, S., Stadler, G., and Wohlmuth, B. I. "A primal-dual active set algorithm for three-dimensional contact problems with coulomb friction". *SIAM Journal on Scientific Computing* 30.2 (2008), pp. 572–596. [DOI], [OAI] a.

[10] Ihlenburg, F. and Babuška, I. "Finite element solution of the Helmholtz equation with high wave number Part I: The h-version of the FEM". *Computers & Mathematics with Applications* 30.9 (1995), pp. 9–37. [DOI].

[11] Jones, S. and Legrand, M. "Forced vibrations of a turbine blade undergoing regularized unilateral contact conditions through the wavelet balance method". *International Journal for Numerical Methods in Engineering* 101.5 (2014), pp. 351–374. [DOI], [OAI] a.

[12] Jones, S. and Legrand, M. "On Solving One-Dimensional Partial Differential Equations With Spatially Dependent Variables Using the Wavelet-Galerkin Method". *Journal of Applied Mechanics* (2013). [DOI], [OAI] a.

[13] Legrand, M., Junca, S., and Heng, S. "Nonsmooth modal analysis of a N-degree-of-freedom system undergoing a purely elastic impact law". *Communications in Nonlinear Science and Numerical Simulation* 45 (2017), pp. 190–219. [DOI], [OAI] a.
[14] Legrand, M. and Pierre, C. “Compact weighted residual formulation for periodic solutions of systems undergoing unilateral contact and frictional occurrences”. *10th European Nonlinear Oscillations Conference (ENOC 2022)*. Lyon, France, 2022. [OAI] a.

[15] Lu, T. and Legrand, M. “Nonsmooth modal analysis via the boundary element method for one-dimensional bar systems”. *Nonlinear Dynamics* 107.1 (2022), pp. 227–246. [DOI], [OAI] a.

[16] Lu, T. and Legrand, M. “Nonsmooth modal analysis with boundary element method”. *XI International Conference on Structural Dynamics*. Greece, 2020, pp. 205–212. [DOI], [OAI] a.

[17] MacNeal, R. H. “A hybrid method of component mode synthesis”. *Computers & Structures* 1.4 (1971). Special Issue on Structural Dynamics, pp. 581–601. [DOI], [OAI] a.

[18] Nacivet, S., Pierre, C., Thouverez, F., and Jézéquel, L. “A dynamic Lagrangian frequency–time method for the vibration of dry-friction-damped systems”. *Journal of Sound and Vibration* 265.1 (2003), pp. 201–219. [DOI], [OAI] a.

[19] Nayfeh, A. H. and Balachandran, B. *Applied nonlinear dynamics: analytical, computational, and experimental methods*. Wiley, 2008. [DOI].

[20] Nunes, A. W., da Silva, S., and Ruiz, A. “Exact general solutions for the mode shapes of longitudinally vibrating non-uniform rods via Lie symmetries”. *Journal of Sound and Vibration* 538 (2022), p. 117216. [DOI].

[21] Powell, M. J. *A Fortran subroutine for solving systems of nonlinear algebraic equations*. 1968. [URL].

[22] Rubin, S. “Improved component-mode representation for structural dynamic analysis”. *AIAA journal* 13.8 (1975), pp. 995–1006. [DOI].

[23] Salles, L., Blanc, L., Thouverez, F., and Gouskov, A. M. “Dynamic analysis of fretting-wear in friction contact interfaces”. *International Journal of Solids and Structures* 48.10 (2011), pp. 1513–1524. [DOI], [OAI] a.

[24] Salles, L., Blanc, L., Thouverez, F., and Gouskov, A. M. “Dynamic Analysis of Fretting-Wear in Friction Contact Interfaces”. *Journal of Engineering for Gas Turbines and Power* 132.1 (2010). [DOI], [OAI] a.

[25] Schreyer, F. and Leine, R. I. “A mixed shooting–harmonic balance method for unilaterally constrained mechanical systems”. *English. Archive of Mechanical Engineering* LXIII.2 (2016), pp. 297–314. [DOI], [OAI] a.

[26] Shi, Y. “Computation of nonlinear modes of vibration of systems undergoing unilateral contact through the semi-smooth Newton approach”. MA thesis. Canada: McGill University, 2016. [URL] a.

[27] Stadler, G. “Semismooth Newton and augmented Lagrangian methods for a simplified friction problem”. *SIAM Journal on Optimization* 15.1 (2004), pp. 39–62. [DOI].

[28] Stewart, D. E. *Dynamics with Inequalities: impacts and hard constraints*. Vol. 59. SIAM, 2011. [DOI].

[29] Thorin, A., Delezoide, P., and Legrand, M. “Non-smooth modal analysis of piecewise-linear impact oscillators”. *SIAM Journal on Applied Dynamical Systems* 16.3 (2017), pp. 1710–1747. [DOI], [OAI] a.

[30] Urman, D., Legrand, M., and Junca, S. “D’Alembert function for exact non-smooth modal analysis of the bar in unilateral contact”. *Nonlinear Analysis: Hybrid Systems* 43 (2021). [DOI], [OAI] a.

[31] Venkatesh, J., Thorin, A., and Legrand, M. “Nonlinear modal analysis of a one-dimensional bar undergoing unilateral contact via the time-domain boundary element method”. *ASME 2017 International Design Engineering Technical Conferences*. Cleveland, United States, 2017. [DOI], [OAI] a.

[32] Yoong, C. “Nonsmooth modal analysis of a finite elastic bar subject to a unilateral contact constraint”. PhD thesis. Canada: McGill University, 2019. [URL] a.

[33] Yoong, C. and Legrand, M. “Nonsmooth modal analysis of a non-internally resonant finite bar subject to a unilateral contact constraint”. *37th IMAC: A Conference and Exposition on Structural Dynamics*. Vol. 1. Nonlinear Structures and Systems. USA, 2019, pp. 1–10. [DOI], [OAI] a.

[34] Yoong, C., Thorin, A., and Legrand, M. “Nonsmooth modal analysis of an elastic bar subject to a unilateral contact constraint”. *Nonlinear Dynamics* (2018), pp. 1–24. [DOI], [OAI] a.