# A compact, equality-based weighted residual formulation for periodic solutions of systems undergoing frictional occurrences

Mathias Legrand, Christophe Pierre

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

## To cite this version:

Mathias Legrand, Christophe Pierre. A compact, equality-based weighted residual formulation for periodic solutions of systems undergoing frictional occurrences. Journal of Structural Dynamics, 2024, 2, pp.144-170. 10.25518/2684-6500.190. hal-04189699v2

HAL Id: hal-04189699

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

Submitted on 29 Mar 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
<figure><img src="image_3.png" /></figure>

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

# A compact, equality-based weighted residual formulation for periodic solutions of systems undergoing frictional occurrences

Mathias Legrand *<fnref n="1" /> and Christophe Pierre †<fnref n="2" />

## Abstract

A very compact weighted residual formulation is proposed for the construction of periodic solutions of oscillators subject to frictional occurrences. Coulomb's friction is commonly expressed as a differential inclusion which can be cast into the complementarity formalism. When targeting periodic solutions, existing algorithms rely on a procedure alternating between the frequency domain, where the dynamics is solved, and the time domain, where friction is dealt with. In contrast, the key idea of the present work is to express all governing equations including friction as equalities, which are then satisfied in a weak integral sense through a weighted residual formulation. The resulting algebraic nonlinear equations are solved numerically using an adapted trust region nonlinear solver and basic integral quadrature schemes. To increase efficiency, the Jacobian of the friction forces is calculated analytically in a piecewise linear fashion. The shape functions considered in this work are the classical Fourier functions. It is shown that periodic solutions with clear multiple sticking and sliding phases can be found with a high degree of accuracy. The equality-based formulation is shown to be effective and efficient, convergence being achieved in all cases considered with low computational cost, including for large numbers of harmonics. Importantly, this new friction formulation does not suffer from the typical limitations or hypotheses of existing frequency-time domain methods for non-smooth systems, such as regularization, penalization, or massless frictional interfaces.

**Keywords:** Non-smooth dynamics; Vibration analysis, Dry friction, Periodic solutions

Received on August 29, 2023, Accepted on February 21, 2024, Published on March 4, 2024

## 1 Introduction

Non-smooth nonlinearities due to dry friction and unilateral contact are ubiquitous in structural engineering systems. Turbomachinery rotors are a prime example of industrial systems that are subject to a wide variety of intermittent contact and frictional occurrences. For example, both friction and contact can occur at the interface of neighboring blade shrouds and as well between the blades and the casing. Friction at the blade root can lead to fretting, and dry friction dampers are commonly used in a number of configurations to mitigate adverse vibrations, for example under-platform or ring dampers [1, 2]. While predicting the dynamic response of non-smooth systems is of great importance and has been the subject of much research over the years, their equations of motion can be remarkably challenging to solve, because Signorini unilateral contact and Coulomb friction conditions are commonly formulated as complementarity problems [3].

The most fundamental and perhaps widely used model of dry friction damping is Coulomb's classical law, though many variations have been developed over the years to model certain phenomena. Coulomb's law nonlinearities are

<footnote n="1">¹Department of Mechanical Engineering, McGill University, Montreal, Canada</footnote>

<footnote n="2">²Department of Mechanical Engineering, Stevens Institute of Technology, Hoboken, USA</footnote>

<footnote>*mathias.legrand@mcgill.ca</footnote>

<footnote>†cpierre@stevens.edu</footnote>
typically expressed as a differential inclusion which can be cast into the complementarity formalism, such that the positive/negative parts of friction forces and the positive/negative parts of relative velocities are the complementary variables. When sliding occurs between two bodies in contact under a normal force $N$, which in this work is assumed to be constant and known, the tangential friction force is constant in magnitude and opposes the relative motion, thus of the form $\pm\mu N$ where $\mu$ is a coefficient of friction. When the relative velocity is zero, then the two bodies are sticking and the friction force is smaller or equal in absolute value to $\mu N$. Thus, in the case of sliding the friction force is piecewise constant in the relative velocity and, importantly, in the case of sticking it is a multivalued function of relative velocity at zero velocity. For a frictional occurrence without separation or impact, the relative displacement and velocity at the friction interface are continuous, but the acceleration and thus the friction force can be discontinuous functions of time [4].

The vibration of dry friction damped systems has been the subject of a very large number of studies over the years, dating back to the early work of Den Hartog [5]. These being too numerous to cite, only representative or sample references are listed in this paper. Many previous works focus on dry friction laws, models and attendant solution methods and algorithms, while others center on applications to structural systems with dry friction damping. Time domain methods are extensively used to analyze frictional occurrences, typically relying on numerical time integration schemes to determine transitions between the different phases of motion and simulate them. Their advantages include the ability to handle multiple transitions between sticking and sliding phases, multiple nonlinear degrees of freedom, and complicated friction damping laws. However, the accurate detection of transitions between phases (sliding, sticking, separation) in the time domain may require advanced time-stepping or event-driven schemes [6]. Contact and friction between deformable bodies which are modeled by the finite element method are typically handled using penalty or Lagrange multipliers methods. The method of augmented Lagrangians is popular for such time simulations of contact and friction, such that the multipliers are the contact forces and the penetration between contacting surfaces is penalized to ensure near zero penetration [7, 8, 9]. Furthermore, time-domain methods are not most suitable for forced vibration problems, as the required time simulation of transients can result in computational costs several orders of magnitude greater than for frequency-domain methods that directly predict steady-state responses. Also, with frequency domain methods, frequency response curves can be generated efficiently by selecting the initial guess at a given frequency as the steady state response at the previous frequency. Finally, unlike frequency domain methods, the time domain formulation is not commonly used for the dynamic condensation of the linear degrees-of-freedom of the system to the nonlinear ones, because of the very large size of the generated algebraic systems. This can result in prohibitively large systems of equations of motion for large-scale structural systems such as dry friction damped bladed disks, even if reduced-order order modeling is performed, e.g., using component mode synthesis. Consequently, for the vibration of systems with frictional occurrences, a frequency domain formulation is often favored, which is the case in the present work.

Methods to calculate the periodic solutions of non-smooth systems include shooting methods [10], which can be effective but suffer from a high sensitivity to the initial guess used in the iterative shooting algorithm and become prohibitively expensive when the number of degrees-of-freedom increases due to the Jacobian computation. The most often used formulation for periodic solutions is the broad class of harmonic balance methods, which are the subject of the present paper [11]. Due to the nature of the dry friction nonlinearity, early frequency-domain studies of dry friction vibration were limited to single harmonic approximations of the response at the frequency of the excitation [12, 13, 14, 15]. Essentially, in the harmonic balance method, the friction force must be expressed in closed form in the frequency domain in terms of the relative velocity over one cycle of oscillations. This is possible only when a single harmonic is retained, as the instants of transition between sticking and sliding states cannot be determined analytically when multiple harmonics are included.

When sliding motion is dominant, the one-harmonic approximation leads to reasonably accurate solutions, but sticking cannot be reproduced accurately since the friction force is then a smooth, harmonic function of time. In order for sticking phases and transitions between sliding and sticking to be captured and to better approximate the variation of the friction force, multiple harmonics of the response must be accounted for. However, except in the case of a smooth, polynomial nonlinearity, the multi-harmonic balance method cannot be fully formulated in the frequency domain, requiring the time domain calculation of the non-smooth nonlinear forces. A number of variants of the multi-harmonic balance method have been developed for nonlinear systems [16], and only those directly relevant to non-smooth systems are reviewed below.

Pierre et al. [17] were first to propose a multi-harmonic formulation of the periodic response of friction damped systems, by extending Lau's incremental harmonic balance method [18] to a non-smooth nonlinearity. They used an approximation of Coulomb's law in which the friction force equals the normal force times the single-valued version of the signum function of the velocity, which can be expanded in terms of a Dirac impulse in the incremental procedure. While the discontinuity in the friction force was captured, its multi-valued nature for zero velocity was not, since for
zero velocity the chosen signum function and thus the friction force are zero. Therefore, here again sticking phases were not adequately modeled. However, circumventing the difficulty posed by the non-uniqueness of the friction force versus the sliding velocity enabled the consideration of multiple harmonics of the excitation in the response and the reasonable approximations of motions with at least two sticking phases per cycle and the associated stick-slip transitions.

The Alternating Frequency/Time (AFT) method subsequently developed by Cameron and Griffin [19] represented a breakthrough for the forced vibration of systems with general nonlinearities, particularly those featuring the hysteresis behavior caused by frictional occurrences. Recognizing that the friction force cannot be expressed in closed form in the frequency domain, in the AFT procedure a Fourier transform of the equations of motion is first taken, resulting in a system of nonlinear algebraic<fnref n="1" /> equations in the Fourier coefficients. At each iteration of the nonlinear Newton-Raphson solution algorithm, the nonlinear friction force is first calculated in the time domain from the displacement and velocity time histories, and then transformed back into the frequency domain via a DFT, until convergence is reached. The AFT method is general and versatile and allows for the inclusion of multiple harmonics in the response. In fact, it can be viewed as a improved variant of the original multi-harmonic balance formulation in which, at each iteration, nonlinear forces are calculated in the time domain and transformed back into the frequency domain. In [19], the AFT was applied successfully to an elastic-perfectly plastic damper without mass at the frictional interface, a modeling simplification that allows for the friction force to be a continuous, piecewise linear function of the relative displacement at the friction interface, and hence a continuous function of time. Such regularization of the friction force in terms of the displacement preserves Coulomb's law but removes the multi-valued nature of the damper force, henceforth greatly facilitating convergence. In the massless frictional interface model, the friction force depends on the history of the motion and features a hysteresis loop in terms of the damper's relative displacement. A comprehensive discussion of the suitability of the AFT method for non-smooth nonlinearities of contact and friction can be found in [20].

Building on the AFT method, the hybrid frequency-time (HFT) domain method was developed for structural systems with dry friction damping, which relied on the accurate calculation of friction forces from the time histories of the system variables based on criteria for precisely detecting sliding and sticking phases and their transitions [21]. More robust nonlinear Broyden and hybrid Powell solution algorithms were used [22]. The HFT method was shown to handle successfully structures with large number of one-dimensional dampers and damper models with variable normal load and separation at the contact point. However, the method was only applied to flexible dampers with a massless interface, for which friction forces are continuous, single-valued functions of relative displacements, and one can expect the HFT to suffer the same convergence problems as the AFT for dampers with mass.

Another major contribution was the development of the Dynamic Lagrangian Frequency/Time (DLFT) method [23]. Like the AFT and HFT, it is a multi-harmonic formulation which requires the calculation of the contact and friction conditions in the time domain at each iteration of the nonlinear solver via a Discrete Fourier Transform. However, the DLFT method is more advanced, as it is an adaptation to the frequency domain framework of the augmented Lagrangians formulation used successfully in the time domain for finite element models of structures in frictional contact. Namely, the dynamic Lagrangians are the nonlinear contact forces obtained from the equations of motion in the frequency domain. Coulomb friction and non-penetration conditions are satisfied with the addition of a penalization on the difference between the interface displacements calculated by the nonlinear solver in the frequency domain and those calculated in the time domain from the nonlinear contact forces. Unlike the AFT and the HFT, the dynamic Lagrangians-based DLFT enables the calculation of friction forces for frictional interfaces with mass. Hence the method is the first to handle Coulomb's law discontinuity of the friction force in the relative velocity, making it the most effective to date.

The DLFT method has been used successfully to calculate the response of complex structures such as bladed disks with under-platform dampers, without requiring either the regularization of the friction law or the assumption of a flexible damper with a massless interface [24]. Also see [25] for a comparative analysis of the DLFT and the linear complementarity problem formulation for a unilateral contact problem. However, the DLFT uses a penalty coefficient, the value of which can influence the rate of convergence and affect the response amplitudes of the converged solutions, such effects being greater when fewer harmonics are retained [26]. Since the adequate range of the penalty coefficient is not known a priori, it is advisable to study its effects in terms of the number of harmonics and to validate sample results with time integrations. Thus, while the DLFT is the most effective and advanced frequency-domain formulation to handle Coulomb's friction law without assumptions, selecting the penalty parameter requires attention.

The great majority of published works on the vibration of dry friction damped systems using multi-harmonic balance methods is based on some kind of a "smoothing" procedure, in order to overcome the difficulty associated with the non-smooth and discontinuous multi- or set-valued nature of Coulomb's friction law. One common approach

<footnote n="1">¹By algebraic, it is meant equations that are not time dependent though not necessarily polynomial functions of the unknowns.</footnote>
in vibrations is to directly approximate Coulomb's discontinuous law with a smoother function which is continuous and single-valued—and in cases, regular—in the velocity at the frictional interface. Examples of regularized functions include the hyperbolic tangent function [27, 28, 29] and the arctangent function [30]. Another common smoothing approximation is to use a piecewise linear function with finite slope in the approximated sticking region [31, 32, 33]. Then, the friction force is a continuous, single-valued function of the damper relative velocity, yielding an approximation of sticking—without true sticking—in a range of small relative velocities. The level of approximation can be adjusted with a parameter of the regularized function, which determines how much the friction force varies versus the velocity in the neighborhood of zero velocity. Note that strictly speaking, this is not a regularization since the friction force is not differentiable in the velocity at the stick-slip transitions. Also see [34, 35, 36, 20] and the references therein for a variety of other regularized friction models, and as well the interesting use of the regularizing sigmoid function in robotics [37, 38].

For the important class of structures of bladed disks with dry friction dampers, a widely used model is to add a contact stiffness at the friction contact point, such that the interface becomes compliant while also being effectively massless, or equivalently to model interfaces with massless compliant finite elements. This model enables the circumvention of the set-valuedness of Coulomb's law. A few of the many works based on flexible, massless frictional interfaces can be found in [10, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]; also see the comprehensive review paper on bladed disks with frictional and contact interfaces [2]. For a frictional interface that is flexible and massless, there is no approximation of Coulomb's friction law and true sticking does occur. However, great simplification results for the solution algorithms, because the friction force is piecewise linear and continuous in the relative displacement of the neighboring degree-of-freedom, featuring a hysteresis loop, hence again eluding the challenge caused by a multi-valued friction force. Note that if the interface model calls for the inclusion of linear damping such as viscous damping in parallel with the interface's stiffness element, then the friction force can no longer be expressed in closed form in terms of the displacement. There are variations of this model, including Dahl's model to represent micro-slip effects [51, 52]. While the assumption of a flexible, massless frictional interface may not seem a priori physical, since friction occurs between bodies which have mass, the introduction of a contact flexibility has been shown to account for the rough and wavy nature of interface surfaces and thus be an appropriate model for dry friction dampers [2]. The validity of the model should be assessed on a case by case basis. For instance, while the effect of damper mass has received very little attention in the literature, Ferri and Heck [53] showed that it can have a significant effect on the response of a blade with an attached friction damper.

To the authors' knowledge, there are very few studies [24, 54] which consider Coulomb's classical law without either regularizing it or assuming a massless friction interface. These are based on the DLFT, which appears to be the only method that produces converged results without requiring explicit smoothing. However, as mentioned above the DLFT introduces a penalization at the contact point to prevent penetration.

A vast body of literature also exists in the area of mechanics of contact and friction for the numerical simulation of non-smooth dynamical systems [55, 6], a significant part of which focuses on rigid bodies in contact. While some of this material is not directly relevant to the present study, a very interesting approach is to convert the complementarity formulations of unilateral contact and friction, which are governed by a system of equalities and inequalities and complementarity conditions, into an equivalent set of non-smooth equalities. These equalities are expressed as the zero-level set of certain functions which have been shown to exist for frictional contact problems [56, 57], and such equivalent non-smooth equalities have been shown to lead to a simpler formulation and effective solution method for the three-dimensional frictional contact problem [58, 59]. This approach, which expresses the non-smooth relationships in an implicit fashion, has also been used effectively for a simplified friction problem [60]. Also see [6, 61] for exhaustive treatments and formulations of non-smooth dynamical systems. Recently, an equivalent equality was successfully formulated for the periodic responses and modal analysis of a one-dimensional bar with a Signorini contact condition [62].

In this paper, the above mentioned equality-based formulations of Coulomb's classical friction law are implemented within the formalism of the harmonic balance method. Namely, the friction law is written as a non-smooth equality and the equations governing the vibration response are satisfied in a weak integral sense through a weighted residual formulation. In Section 2, equivalent equalities for Coulomb's friction derived in the above cited references are presented. To demonstrate the generality of the approach, equivalent equalities are presented in Section 3 for common regularized and approximate friction laws and the single-valued signum function approximation. In Section 4, a two-degree-of-freedom system with a one-dimensional frictional interface is considered, for which the equivalent equality acts as a non-smooth condition in the friction force and velocity at the friction degree-of-freedom, simply augmenting the equations of motion. Periodic responses to a single-frequency excitation are sought by performing a Ritz-Galerkin procedure, whereby all governing equations, including the friction condition as an equality, are satisfied in a weak integral sense in time, over one period of the external force. Contrary to existing formulations, this compact,
weighted residual formulation does not require switching from the frequency to the time domain at each iteration to calculate the nonlinear forces and attendant states via numerical integration. The resulting algebraic nonlinear equations are solved numerically using an adapted trust region nonlinear solver, and at each iteration the integrals of the non-smooth functions arising from the weak formulation are computed using a classical quadrature scheme. To make computations efficient, the Jacobian of the nonlinear forces is expressed in an exact, piecewise linear fashion, corresponding to the various states of the system. Results presented in Section 5.1 include frequency response curves for various specified and constant normal force values and time histories of system displacements, velocities, and friction force. In Section 5.2, a convergence analysis shows an improved accuracy when the number of harmonics is increased, and results are validated with time integration in Section 5.3. To demonstrate the versatility of the approach and assess the common smoothing assumptions of existing techniques, the effect of damper mass on periodic responses is investigated in Section 5.4, including the common case of a massless flexible frictional interface, and an analysis of the piecewise linear approximation of the friction force is conducted in Section 5.5. The key features of the equality-based formulation are discussed in Section 6, and the approach is compared with existing frequency-time domain methods for the forced vibration of friction damped systems.

## 2 Equality-based formulation for Coulomb friction conditions

In order to introduce the proposed formulation for non-smooth systems with frictional occurrences, the basic model of one-dimensional Coulomb friction shown in Figure 1 is considered. The mass $m$ rests on a horizontal surface and its displacement in the $x$ direction is denoted by $u(t)$. The normal force on the mass in the downward vertical direction is denoted by $N$ (such that $N > 0$) and the tangential friction force in the $x$ direction, resisting the motion of the mass, is denoted by $r(t)$.

Coulomb's classical friction model reads

$$ \begin{cases} \dot{u} = 0 & \Rightarrow |r| \le \mu N \\ \dot{u} \neq 0 & \Rightarrow |r| = \mu N \text{ and } \exists \alpha \ge 0 \mid r = -\alpha \dot{u} \end{cases} \qquad (1) $$

where $\mu$ is the coefficient of friction. The admissible set of points in the $(\dot{u}, r)$ plane for the two conditions for sticking and sliding of the mass in Equation (1) can be classically represented as shown in Figure 1. It shows that the relationship between $r$ and $\dot{u}$, loosely denoted $r(\dot{u})$, is not single-valued but set-valued, since the set of acceptable values of $r$ when $\dot{u} = 0$ does not reduce to a single point. Coulomb's conditions in Equation (1) can be recast into an

<figure><img src="image_5.png" /><figcaption>Fig. 1: Coulomb friction conditions.</figcaption></figure>

equivalent non-smooth equality [58, 60, 63, 6, 64, 9]. For instance, consider the two-dimensional level-set function

$$ \Psi(\dot{u}, r) = \dot{u} + \min(0, \rho(r + \mu N) - \dot{u}) + \max(0, \rho(r - \mu N) - \dot{u}) \qquad (2) $$

where $\rho > 0$ is an optional normalization parameter with default value 1, which if needed can be selected so that the quantities $\rho(r + \mu N)$, $\rho(r - \mu N)$, and $\dot{u}$ are of similar magnitudes in order to facilitate computations.<fnref n="2" /> The function in Equation (2) can be expressed as

$$ \Psi(\dot{u}, r) = \begin{cases} \rho(r - \mu N) & \text{if } \rho(r + \mu N) - \dot{u} > 0 \text{ and } \rho(r - \mu N) - \dot{u} > 0, \\ \dot{u} & \text{if } \rho(r + \mu N) - \dot{u} > 0 \text{ and } \rho(r - \mu N) - \dot{u} < 0, \\ \rho(r + \mu N) & \text{if } \rho(r + \mu N) - \dot{u} < 0 \text{ and } \rho(r - \mu N) - \dot{u} < 0. \end{cases} \qquad (3) $$

<footnote n="2">²For all cases in this paper, $\rho = 1$ and there were no convergence or computational issues.</footnote>
Note that the case $\rho(r + \mu N) - \dot{u} < 0$ and $\rho(r - \mu N) - \dot{u} > 0$ is never satisfied, because subtracting the latter quantity from the former yields $2\rho\mu N < 0$, which is not possible. The zero-level set $\Psi_0$ of the function $\Psi(\dot{u}, r)$ is defined as $\Psi_0 \equiv \{(\dot{u}, r) | \Psi(\dot{u}, r) = 0\}$. Accordingly, we have from Equation (3), for the zero-level set:

$$ \begin{array}{l} \Psi(\dot{u}, r) = 0 \Rightarrow \begin{array}{l} r = +\mu N \text{ and } \dot{u} \le 0 \quad (\text{negative sliding}) \\ \text{or } \dot{u} = 0 \text{ and } -\mu N \le r \le \mu N \quad (\text{sticking}) \\ \text{or } r = -\mu N \text{ and } \dot{u} \ge 0 \quad (\text{positive sliding}) \end{array} \end{array} \quad (4) $$

which can be recast into

$$ \Psi_0 = \underbrace{([0, \infty, 0] \times \{\mu N\})}_{\text{negative sliding}} \cup \underbrace{([0, \infty] \times [-\mu N, \mu N])}_{\text{sticking}} \cup \underbrace{([0, +\infty[ \times \{-\mu N\})}_{\text{positive sliding}}. \quad (5) $$

Hence, it is found that the above set $\Psi_0$ exactly describes the Coulomb friction set defined by Equation (1) and shown in Figure 1. In other words, the non-smooth equality $\dot{u} + \min(0, \rho(r + \mu N) - \dot{u}) + \max(0, \rho(r - \mu N) - \dot{u}) = 0$ is another version of the considered Coulomb's friction law. $\Psi(\dot{u}, r)$ and $\Psi_0$ are both plotted in Figure 2(a). For the three cases in Equation (3), $\Psi(\dot{u}, r)$ reduces to two planar surfaces with slope $\rho$ in the $r$ direction and one planar surface with unit slope in the $\dot{u}$ direction.

There are multiple ways to implicitly represent Coulomb's classical law as a level set. Another example of such an equality condition proposed by Hüeber et al. [59] is

$$ \Psi(\dot{u}, r) = r \max(\mu N, |r - \rho \dot{u}|) - \mu N(r - \rho \dot{u}), \quad \rho > 0, \quad (6) $$

a graphical representation of which is depicted in Figure 2(b). Equation (6) can be shown to be equivalent to Equation (1), as follows. If $|r - \rho \dot{u}| < \mu N$, then $\Psi(\dot{u}, r) = \mu N \rho \dot{u}$, corresponding to the portion of the plane with positive slope proportional to $\dot{u}$ in Figure 2(b); hence $\Psi = 0$ yields $\dot{u} = 0$, corresponding to the sticking condition. Similarly, if $|r - \rho \dot{u}| \ge \mu N$, then $\Psi(\dot{u}, r) = (r \mp \mu N)|r - \rho \dot{u}|$, and $\Psi = 0$ yields $r = \pm \mu N$, corresponding to the sliding conditions. Figure 2(b) features more complicated surfaces than Figure 2(a), but the same intersection $\Psi_0$ with the zero-level plane.

<figure><img src="image_6.png" /><figcaption>Fig. 2: Exact implicit description of one-dimensional Coulomb's friction law. The intersection of the surface represented by the level-set function $\Psi(\dot{u}, r)$ [blue] with the zero-level plane [orange] yields the corresponding zero-level set $\Psi_0$ [solid black line] defined in Equation (5), for $\mu N = 1$.</figcaption></figure>

## Proximal points and projectors

A classical equality-based formulation of Coulomb's friction comes from Convex Analysis and the notion of proximal points and projection operators. Much literature exists on this topic [64] and without going into the details, the one-dimensional Coulomb's friction law considered in this work can be expressed as follows:

$$ r = \operatorname{proj}_{D(\mu N)}(r - \rho \dot{u}), \quad \rho > 0 \quad (7) $$
where $D(\mu N) \equiv \{r | |r| \le \mu N\}$ is Coulomb's friction cone (here, the interval $-\mu N \le r \le \mu N$). For the considered configuration, the projection operation explicitly reads³

$$ \mathrm{proj}_{D(\mu N)}(r - \rho \dot{u}) = \begin{cases} + \mu N & \text{if } r - \rho \dot{u} \ge \mu N \\ r - \rho \dot{u} & \text{if } |r - \rho \dot{u}| \le \mu N \\ - \mu N & \text{if } r - \rho \dot{u} \le -\mu N \end{cases} \qquad (8) $$

For $\rho = 1$, the level-set function $\Psi(\dot{u}, r) = r - \mathrm{proj}_{D(\mu N)}(r - \rho \dot{u})$ and the corresponding zero-level set $\Psi_0$ are exactly like their counterparts in Figure 2(a).

# 3 Regularized and approximate laws

In order to avoid the theoretical and computational difficulties associated with the multi-valued nature of the friction force in the case of sticking, it is common practice to use a smoothing procedure for Coulomb's friction law, in which the discontinuous variation of the friction force in Figure 1 is approximated with a single-valued continuous function. These regularized or approximate laws can also be cast in terms of an equivalent equality, which is smooth in some cases but non-smooth in others, as is shown below for commonly used regularizations.

## 3.1 Hyperbolic tangent and arctangent regularizations

Coulomb's law is frequently approximated with the hyperbolic tangent function $r \approx r_\alpha = -\mu N \tanh(\alpha \dot{u})$ where $\alpha > 0$ is a parameter that can be chosen to adjust the rate of variation of $r$ in terms of $\dot{u}$ in the neighborhood of $\dot{u} = 0$ [28]. As $\alpha \to \infty$, the regularized function $r_\alpha$ tends to Coulomb's classical law. Note that because $r_\alpha$ is a single-valued and continuous function of $\dot{u}$, there is no possibility of true sticking motion, since $r_\alpha \neq 0$ except for $\dot{u} = 0$. As a result, the degree of validity of such regularization approaches depends on the appropriate selection of the parameter $\alpha$ to approximate the sticking motions, which might be case dependent. An implicit version is simply given by the zero-level set of the function

$$ \Psi(\dot{u}, r) = r + \mu N \tanh(\alpha \dot{u}) \qquad (9) $$

which in this case is smooth. The surface $\Psi(\dot{u}, r)$ in Equation (9) and its zero-level set are depicted in Figure 3(a) for $\alpha = 20$, which readily shows the regularized approximation of Coulomb's law.

Another frequently used regularization of Coulomb's law is based on the arctangent function [30], in which case the equivalent friction equality takes the form

$$ \Psi(\dot{u}, r) = r + \mu N \frac{2}{\pi} \arctan(\alpha \dot{u}). \qquad (10) $$

Again here the level-set function is smooth. The function and its zero-level set are depicted in Figure 3(b) for $\alpha = 20$, showing a similar yet smoother regularization than that of the hyperbolic tangent function for the same value of $\alpha$.

## 3.2 Piecewise linear approximation

Another approach to make Coulomb's law smoother is to approximate the discontinuity of $r$ versus $\dot{u}$ at $\dot{u} = 0$ in Figure 1 with a line of finite slope [33]:

$$ \left\{ \begin{array}{ll} \dot{u} \le -\delta & \Rightarrow r = \mu N \\ -\delta \le \dot{u} \le \delta & \Rightarrow r = -\mu N \dot{u} / \delta \\ \dot{u} \ge \delta & \Rightarrow r = -\mu N \end{array} \right. \qquad (11) $$

where $\delta > 0$ and as $\delta \to 0$, Equation (11) tends to Coulomb's law (1). This piecewise linear law is single valued, and again it does not capture true sticking since $r = 0$ if and only if $\dot{u} = 0$. An equality that is equivalent to the approximate law in Equation (11) can be readily formulated using the level-set function

$$ \Psi(\dot{u}, r) = (\mu N) \dot{u} / \delta + r + \min(\mu N(1 - \dot{u} / \delta), 0) + \max(-\mu N(1 + \dot{u} / \delta), 0) \qquad (12) $$
<figure><img src="image_7.png" /><figcaption>**Fig. 3:** Regularized one-dimensional Coulomb's friction law: level-set function $\Psi(\dot{u}, r)$ [blue], zero-level plane [orange] and corresponding zero-level set [solid black line].</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 4:** Piecewise linear approximation of Coulomb's one-dimensional friction, see Equation (12): function $\Psi(\dot{u}, r)$ [blue], zero-level plane [orange] and corresponding zero-level set [solid black line] for $\delta = 1$.</figcaption></figure>

which in this case is not smooth since it involves min and max functions. Figure 4 shows that its zero-level set is the piecewise linearly approximated law in Equation (11).

Note that the piecewise linear approximation is not exactly a regularization because the function $r$ in Equation (11) is not smooth as it is not differentiable at $\dot{u} = \pm\delta$, and Equation (12) is not smoother than Equation (2). However, the resulting zero-level set is single-valued when seen in the $(\dot{u}, r)$ plane, rather than multi-valued as is the zero-level set for Coulomb's law in Figure 2(a). It is interesting that the piecewise linear approximation amounts to inserting a linear dashpot between the mass $m_2$ and the ground in Figure 1. In the approximate "sticking" phases of the motion, frictional contact is replaced with a viscous damper of constant inversely proportional to the range of "sticking" velocities $-\delta \le \dot{u} \le \delta$. Clearly, the smaller $\delta$, the better the approximation of sticking. Also, since contact then occurs between the ground and a massless viscous damping element, the piecewise linear approximation has the effect of assuming that the friction interface is massless. While this allows for a smoother, single-valued model of the friction force and leads to better convergence of the nonlinear solvers, it remains an approximation of Coulomb's sticking.

<footnote n="3">³The prox operator could be used as well since the two operators are equivalent here.</footnote>
## 3.3 Single-valued signum function approximation

Another previously used approximation of Coulomb's classical law [17], which does not belong to the realm of regularization, consists of approximating the friction force with the single-valued signum function, such that it vanishes only when the velocity is equal to zero:

$$r = -\mu N \operatorname{sign} \dot{u} = \begin{cases} \mu N & \text{if } \dot{u} < 0 \\ 0 & \text{if } \dot{u} = 0 \\ -\mu N & \text{if } \dot{u} > 0. \end{cases} \qquad (13)$$

Then, the corresponding level-set function is simply $\Psi(\dot{u}, r) = r + \mu N \operatorname{sign}(\dot{u}) = 0$, represented in Figure 5 along with its zero-level set. Here, the friction force is discontinuous in the velocity at $\dot{u} = 0$, but contrary to Coulomb's law r

<figure><img src="image_9.png" /><figcaption>Fig. 5: Single-valued signum function approximation of Coulomb's one-dimensional friction. Function $\Psi(\dot{u}, r)$ [blue], zero-level plane [orange], and corresponding zero-level set [solid black line] showing the single-point intersection at $r = 0$ for $\dot{u} = 0$.</figcaption></figure>

is single-valued in the velocity, such that $r = 0$ for $\dot{u} = 0$ only. Thus the mass is always sliding and sticking is not captured. This is seen in Figure 5, which shows that the zero-level set consists of a single point for $\dot{u} = 0$, jumping to $r = \pm \mu N$ for $\dot{u} \neq 0$. As a result, the sliding phases of the motion are adequately captured but the sticking phases are not, and the mass is found to oscillate with back and forth sliding motions of small amplitudes during "sticking". This model has significant limitations to capture stick-slip behavior, in contrast with the arctangent and hyperbolic tangent regularizations and the piecewise linear approximation, which allow for a better approximation of sticking.

# 4 Weighted residual formulation for periodic solutions

In this section, periodic solutions are sought using the equality-based friction formulation presented in Section 2, using a solution procedure based on a weak formulation and subsequent Ritz-Galerkin projection. While the weighted residual formulation is presented here for a simple two-degree-of-freedom system with one-dimensional friction, it is believed that from a theoretical point of view, the approach can be extended to interfaces with both contact and two-dimensional friction and to multiple contact interfaces.

## 4.1 Two-degree-of-freedom system with one-dimensional friction

The system is shown in Figure 6 and features two degrees-of-freedom of mass $m_1$ and $m_2$, with a friction condition on mass $m_2$. The system is subject to external forces $f_1(t)$ and $f_2(t)$, periodic with natural frequency $\omega$, at the first and second degrees-of-freedom, respectively. Additionally, to limit response amplitudes, viscous damping is added to the system in the form of two dashpots of damping constants $d_1$ and $d_2$, which are parallel to the springs $k_1$ and $k_2$, respectively. Using the zero-level set of the function in Equation (2) at the second degree-of-freedom, and denoting
<figure><img src="image_10.png" /><figcaption>Fig. 6: Two-degree-of-freedom mass-spring system with Coulomb friction at the second degree-of-freedom. Stiffness $k_i$ and viscous damping $d_i$, $i = 1, 2$ are directly included in the springs.</figcaption></figure>

the displacements of the two degrees-of-freedom by $x_1(t)$ and $x_2(t)$, the equations governing the dynamics of the system are<fnref n="4" />

$$ \Psi_1(x_1, x_2, r) = m_1 \ddot{x}_1 + (d_1 + d_2) \dot{x}_1 - d_2 \dot{x}_2 + (k_1 + k_2) x_1 - k_2 \dot{x}_2 - f_1 = 0 \quad (14a) $$

$$ \Psi_2(x_1, x_2, r) = m_2 \ddot{x}_2 + d_2 \dot{x}_2 - d_2 \dot{x}_1 + k_2 x_2 - k_2 x_1 - f_2 - r = 0 \quad (14b) $$

$$ \Psi_r(\dot{x}_2, r) = \dot{x}_2 + \min(0, \rho(r + \mu N) - \dot{x}_2) + \max(0, \rho(r - \mu N) - \dot{x}_2) = 0 \quad (14c) $$

where the parameter $\rho > 0$ is specified by the user. System (14) can be read as a set of Differential Algebraic Equations [64]. The unknowns in Equation (14) are the three functions of time $x_1$, $x_2$, and $r$. In the dry friction damping literature, the system in Figure 6 would be regarded as a model of a one-degree-of-freedom system with an attached friction damper of mass $m_2$, stiffness $k_2$, and viscous dissipation $d_2$. Thus, in the subsequent sections, the second degree-of-freedom may be referred to as the friction damper.

## 4.2 Solution procedure for periodic solutions

A solution procedure based on a weak formulation and a Ritz-Galerkin projection is presented to search for the periodic solutions of the system governed by Equation (14). To recall, the two equations of motion for the two degrees-of-freedom, Equations (14a) and (14b), are augmented by the non-smooth equality for one-dimensional friction, Equation (14c).

Periodic solutions of period $T = 2\pi/\omega$ are sought in a weighted residual sense. Classically, all unknowns of the problem are discretized by expanding them on an appropriate truncated basis of $T$-periodic functions with $N_\phi$ members, as

$$ x_1(t) \approx x_{1h}(t) = \sum_{k=1}^{N_\phi} x_{1k} \phi_k(t), \quad x_2(t) \approx x_{2h}(t) = \sum_{k=1}^{N_\phi} x_{2k} \phi_k(t), \quad r(t) \approx r_h(t) = \sum_{k=1}^{N_\phi} r_k \phi_k(t) \quad (15) $$

where $x_{1h}(t)$, $x_{2h}(t)$ and $r_h(t)$ are the $N_\phi$-term series approximations of the unknowns and $\phi_k(t)$ is the $k$th basis function, which can be viewed as a shape function of time. The $3N_\phi$ coefficients $x_{1k}$, $x_{2k}$ and $r_k$ are the unknowns of the discretized problem, and they can be regarded as the participation of the $k$th shape function $\phi_k$ in the approximate solution.

Commonly, the Fourier functions used in the Harmonic Balance Method are chosen for the truncated basis $\{\phi_k\}_{k=1,...,N_\phi}$. However, this is not a limitation of the formulation and other basis functions can be used, either smooth ones or those featuring some level of non-smoothness, such as wavelets [65]. Depending on the smoothness of the selected basis functions, time derivatives might be obtained either by pointwise differentiation in time or by expanding them on a less smooth basis and relating them to the differentiated quantity in a weak sense.

The next step in the procedure is to substitute the expansions for the unknown approximate functions $x_{1h}$, $x_{2h}$, $r_h$ and possibly their derivatives into the two equations of motion and the non-smooth friction equality condition. The displacements are differentiated term-wise as

$$ \dot{x}_1(t) \approx \dot{x}_{1h}(t) = \sum_{k=1}^{N_\phi} x_{1k} \dot{\phi}_k(t), \quad \dot{x}_2(t) \approx \dot{x}_{2h}(t) = \sum_{k=1}^{N_\phi} x_{2k} \dot{\phi}_k(t), \quad (16a) $$

$$ \ddot{x}_1(t) \approx \ddot{x}_{1h}(t) = \sum_{k=1}^{N_\phi} x_{1k} \ddot{\phi}_k(t), \quad \ddot{x}_2(t) \approx \ddot{x}_{2h}(t) = \sum_{k=1}^{N_\phi} x_{2k} \ddot{\phi}_k(t). \quad (16b) $$

<footnote n="4">⁴Note that the notations $\Psi_1(x_1, x_2, r)$ and $\Psi_2(x_1, x_2, r)$ are a bit loose as the corresponding velocities and accelerations are also involved in the expressions. However, this is sufficiently clear and unambiguous for the upcoming derivations.</footnote>
It is important to note that the friction condition does not feature any time derivatives of the tangential friction force $r$ (or of the normal force $N$, in the case it were not a constant). Hence the procedure does not require the differentiation of these non-smooth functions, which is an important advantage. Substitution leads to the nonzero residuals

$$ \Psi_1(x_{1h}, x_{2h}, r_h) = m_1 \ddot{x}_{1h} + (d_1 + d_2) \dot{x}_{1h} - d_2 \dot{x}_{2h} + (k_1 + k_2) x_{1h} - k_2 x_{2h} - f_1 \quad (17a) $$

$$ \Psi_2(x_{1h}, x_{2h}, r_h) = m_2 \ddot{x}_{2h} + d_2 \dot{x}_{2h} - d_2 \dot{x}_{1h} + k_2 x_{2h} - k_2 x_{1h} - f_2 - r_h \quad (17b) $$

$$ \Psi_r(\dot{x}_{2h}, r_h) = \dot{x}_{2h} + \min(0, \rho(r_h + \mu N) - \dot{x}_{2h}) + \max(0, \rho(r_h - \mu N) - \dot{x}_{2h}) \quad (17c) $$

where $\Psi_1(x_{1h}, x_{2h}, r_h)$, $\Psi_2(x_{1h}, x_{2h}, r_h)$ and $\Psi_r(\dot{x}_{2h}, r_h)$ are the residuals of the two equations of motion and the friction condition, respectively, for the approximation of order $N_\phi$. Recall that either Equation (2) or Equation (6) or any other equivalent friction equality condition could be used for the expression of $\Psi_r$, including the regularized and approximate laws of Section 3. Next, a Ritz-Galerkin projection is performed, by requiring the residuals to be orthogonal to the set of shape functions $\{\phi_k\}_{k=1,...,N_\phi}$. This is written as

$$ \int_{0}^{T} \phi_k(t) \Psi_1(x_{1h}(t), x_{2h}(t), r_h(t)) dt = 0 \qquad k = 1, \dots, N_\phi \quad (18a) $$

$$ \int_{0}^{T} \phi_k(t) \Psi_2(x_{1h}(t), x_{2h}(t), r_h(t)) dt = 0 \qquad k = 1, \dots, N_\phi \quad (18b) $$

$$ \int_{0}^{T} \phi_k(t) \Psi_r(\dot{x}_{2h}(t), r_h(t)) dt = 0 \qquad k = 1, \dots, N_\phi. \quad (18c) $$

Equation (18) consists of a nonlinear system of $3N_\phi$ algebraic equations with $3N_\phi$ unknowns. However, note that Equations (18a) and (18b) govern the dynamics of masses 1 and 2, respectively, and are thus linear in $x_{1h}$, $x_{2h}$ and $r_h$ (or equivalently, in the coefficients $(x_{1k}, x_{2k}, r_k)$), and that the nonlinearity of the system is solely in Equation (18c) which is non-smooth in $\dot{x}_2$ and $r$ (or equivalently, in the coefficients $(x_{2k}, r_k)$).

The accuracy of the approximate solution is solely dictated by the selection of the basis functions $\phi_k$ and by the number $N_\phi$ of these functions. Classically, the solution of Equation (18) becomes more accurate as $N_\phi$ increases, since the residuals become smaller in a weak sense.

## 4.3 Fourier series approximation in time

The general weighted residual procedure presented above is applied to the two-degree-of-freedom system in Figure 6. Here, the basis functions $\phi_k(t)$ are classically chosen to be the Fourier functions. The Fourier series expansions can be simplified a bit, because in the common case of an external excitation which is harmonic, there are no contributions from the even harmonics or from the constant terms due to the symmetry of the Coulomb friction nonlinearity. Thus, an approximation with $N_\phi$ Fourier terms, where $N_\phi$ is even without loss of generality, reduces to an equal number of the odd cosine and sine harmonics, each ranging from harmonic 1 to harmonic $N_\phi - 1$. The unknown functions $x_{1h}$, $x_{2h}$ and $r_h$ are thus expanded into truncated series of Fourier functions with $N_\phi$ coefficients each, as

$$ x_1(t) \approx x_{1h}(t) = \sum_{k=1}^{N_\phi/2} [x_{1_{2k-1}} \cos(2k-1)\omega t + x_{1_{2k}} \sin(2k-1)\omega t], \quad (19a) $$

$$ x_2(t) \approx x_{2h}(t) = \sum_{k=1}^{N_\phi/2} [x_{2_{2k-1}} \cos(2k-1)\omega t + x_{2_{2k}} \sin(2k-1)\omega t], \quad (19b) $$

$$ r(t) \approx r_h(t) = \sum_{k=1}^{N_\phi/2} [r_{2k-1} \cos(2k-1)\omega t + r_{2k} \sin(2k-1)\omega t]. \quad (19c) $$

For the sake of generality, the external periodic excitations on masses 1 and 2 can also be expanded as Fourier series with odd harmonic terms as

$$ f_1(t) = \sum_{k=1}^{N_\phi/2} [f_{1_{2k-1}} \cos(2k-1)\omega t + f_{1_{2k}} \sin(2k-1)\omega t], \quad (20a) $$

$$ f_2(t) = \sum_{k=1}^{N_\phi/2} [f_{2_{2k-1}} \cos(2k-1)\omega t + f_{2_{2k}} \sin(2k-1)\omega t]. \quad (20b) $$
although in the cases considered in this paper, $f_1(t) = f_1 \cos \omega t$ and $f_2(t) = 0$. For convenience, the following vectors of Fourier series coefficients are defined for the first and the second mass and the friction force: $\mathbf{x}_1 = [x_1, \dots, x_{1N_\phi}]^\top$, $\mathbf{x}_2 = [x_{21}, \dots, x_{2N_\phi}]^\top$, and $\mathbf{r} = [r_1, \dots, r_{N_\phi}]^\top$. Similarly, for the external forces:

$$ \mathbf{f}_1 = [f_{12k-1}, f_{12k}]^\top_{k=1,...,N_\phi/2} \quad \mathbf{f}_2 = [f_{2k-1}, f_{2k}]^\top_{k=1,...,N_\phi/2} \qquad (21) $$

The Ritz-Galerkin projection is performed by substituting Equation (19) into Equation (14), multiplying in turn by the $N_\phi$ shape functions $\{\cos(2k-1)\omega t\}_{k=1,...,N_\phi/2}$ and $\{\sin(2k-1)\omega t\}_{k=1,...,N_\phi/2}$, and integrating over one period of the motion $[0, 2\pi/\omega]$ as in Equation (18). The following system of coupled nonlinear equations is obtained:

$$ \mathbf{A}_{11}\mathbf{x}_1 + \mathbf{A}_{12}\mathbf{x}_2 - \mathbf{f}_1 = \mathbf{0} \qquad (22a) $$

$$ \mathbf{A}_{21}\mathbf{x}_1 + \mathbf{A}_{22}\mathbf{x}_2 - \mathbf{r} - \mathbf{f}_2 = \mathbf{0} \qquad (22b) $$

$$ \mathbf{g}_r(\mathbf{x}_2, \mathbf{r}) = \mathbf{0} \qquad (22c) $$

with $\mathbf{A}_{11} = (k_1 + k_2)\mathbf{I} + (d_1 + d_2)\mathbf{D} - m_1\mathbf{\Omega}^2$, $\mathbf{A}_{12} = -k_2\mathbf{I} - d_2\mathbf{D} = \mathbf{A}_{21}$, $\mathbf{A}_{22} = k_2\mathbf{I} + d_2\mathbf{D} - m_2\mathbf{\Omega}^2$, $\mathbf{I}$ being the identity matrix of appropriate size, and the block diagonal matrices

$$ \mathbf{D} = \mathbf{B}\mathbf{D}\mathbf{A}\mathbf{g}_{k=1,...,N_\phi/2} \begin{bmatrix} 0 & (2k-1)\omega \\ -(2k-1)\omega & 0 \end{bmatrix}, \quad \mathbf{\Omega} = \mathbf{B}\mathbf{D}\mathbf{A}\mathbf{g}_{k=1,...,N_\phi/2} \begin{bmatrix} (2k-1)\omega & 0 \\ 0 & (2k-1)\omega \end{bmatrix}. \qquad (23) $$

Also, the components of $\mathbf{g}_r$ defined in Equation (22c) read

$$ g_{r_{2k-1}}(\mathbf{x}_2, \mathbf{r}) = \int_0^T \Psi_r(\dot{x}_{2h}(t), r_h(t)) \cos(2k-1)\omega t \, dt = 0, \qquad k = 1, \dots, N_\phi/2 \qquad (24a) $$

$$ g_{r_{2k}}(\mathbf{x}_2, \mathbf{r}) = \int_0^T \Psi_r(\dot{x}_{2h}(t), r_h(t)) \sin(2k-1)\omega t \, dt = 0, \qquad k = 1, \dots, N_\phi/2. \qquad (24b) $$

where, in the right hand-side, the dependency on $(\mathbf{x}_2, \mathbf{r})$ is implied from Equation (19). This formulation readily holds for a periodic forces of period $\omega$ which are odd functions of time, including in the case forcing is applied at the friction damper, i.e. $\mathbf{f}_2 \neq 0$. Equations (22a) and (22b) can be solved via substitution, eliminating $\mathbf{x}_1$ and expressing $\mathbf{x}_2$ as an affine relationship with $\mathbf{r}$, that is

$$ \mathbf{x}_2(\mathbf{r}) = (\mathbf{A}_{22} - \mathbf{A}_{21}\mathbf{A}_{11}^{-1}\mathbf{A}_{12})^{-1}(\mathbf{r} + \mathbf{f}_2 - \mathbf{A}_{21}\mathbf{A}_{11}^{-1}\mathbf{f}_1). \qquad (25) $$

Substituting $\mathbf{x}_2(\mathbf{r})$ into Equation (22c), System (22) of dimension $3N_\phi$ reduces to $N_\phi$ non-smooth equations

$$ \mathbf{g}_r(\mathbf{x}_2(\mathbf{r}), \mathbf{r}) = \mathbf{0} \qquad (26) $$

in only the unknown $\mathbf{r}$. This reduced system of nonlinear equations can be solved using a non-smooth Newton or alike solver. Here, the fsolve built-in command in MATLAB® is used, which is based on the dogleg trust region procedure. The tolerance is set to $10^{-6}$. The integrals involving non-smooth terms in Equation (24) are numerically computed using a classical quadrature scheme such as a simple Riemann sum for instance. The integral built-in command in MATLAB® is used with an absolute tolerance of $10^{-6}$. Note that in general, the rate of convergence of the procedure might depend on the parameter $\rho$, which has to be assigned a value in the solvers. For all results presented in this paper, the normalization parameter $\rho$ was taken equal to its default value of 1, as no convergence or computational issues were found. Furthermore, a sensitivity study was carried out in the range $\rho \in [10^{-3}, 5]$, and no notable effect of the value of $\rho$ on the convergence rate was observed.

## 4.4 Piecewise Jacobian of the nonlinear terms

Nonlinear solvers require the calculation of the Jacobian of the nonlinear forces at each iteration, which can become computationally intensive and slow convergence, especially for large numbers of harmonics. Computations can be made considerably more efficient by calculating the Jacobian analytically. In the case of friction, the Jacobian can be written as a piecewise linear analytical expression corresponding to the various states of the system, instead of being approximated numerically via finite differences in the nonlinear solver.

The Jacobian of the non-smooth function $\mathbf{g}$ of the harmonic components of the friction force $\mathbf{r}$ in Equation (26) is denoted as $\nabla_r \mathbf{g}_r$ and from Equation (18c), its elements are:

$$ (\nabla_r \mathbf{g}_r)_{ij} = \frac{\partial g_{ri}}{\partial r_j} = \int_0^T \phi_i(t) \frac{\partial}{\partial r_j} \Psi_r(\dot{x}_{2h}(t), r_h(t)) \, dt \quad i, j = 1, \dots, N_\phi \qquad (27) $$
with $\Psi_r$ defined in Equation (2). The chain rule leads to

$$ \frac{\partial \Psi_r}{\partial r_j} = \frac{\partial \Psi_r}{\partial \dot{x}_2} \frac{\partial \dot{x}_2}{\partial r_j} + \frac{\partial \Psi_r}{\partial r} \frac{\partial r}{\partial r_j} \qquad (28) $$

and the piecewise terms are readily obtained by taking the partial derivatives of $\Psi_r(\dot{x}_2, r)$ in Equation (2) with respect to $\dot{x}_2$ and $r$. The weak derivatives of the min and max functions can be written in terms of the Heaviside step functions as

$$ \frac{d}{dx} \min(0, x) = 1 - H(x), \qquad \frac{d}{dx} \max(0, x) = H(x) \qquad (29) $$

with $H(x) = 0$ for $x \le 0$ and $H(x) = 1$ for $x > 0$. This readily yields

$$ \frac{\partial \Psi_r}{\partial \dot{x}_2} = 1 - (1 - H(\rho(r + \mu N) - \dot{x}_2) - H(\rho(r - \mu N) - \dot{x}_2)) \qquad (30) $$

and

$$ \frac{\partial \Psi_r}{\partial r} = \rho(1 - H(\rho(r + \mu N) - \dot{x}_2)) + \rho H(\rho(r - \mu N) - \dot{x}_2) \qquad (31) $$

both of which can be expressed in piecewise linear form as<fnref n="5" />

$$ \frac{\partial \Psi_r}{\partial \dot{x}_2} = \begin{cases} 0 & \text{and } \frac{\partial \Psi_r}{\partial r} = \begin{cases} \rho & \text{if } \rho(r(t) + \mu N) - \dot{x}_2(t) > 0 \text{ and } \rho(r(t) - \mu N) - \dot{x}_2(t) > 0 \\ 0 & \text{if } \rho(r(t) + \mu N) - \dot{x}_2(t) > 0 \text{ and } \rho(r(t) - \mu N) - \dot{x}_2(t) < 0 \\ \rho & \text{if } \rho(r(t) + \mu N) - \dot{x}_2(t) < 0 \text{ and } \rho(r(t) - \mu N) - \dot{x}_2(t) < 0 \end{cases} \end{cases} \qquad (32) $$

Finally, $\partial r/\partial r_j = \phi_j$ and $\partial \dot{x}_2/\partial r_j$ is calculated from the expansion of $\dot{x}_2$ in Equation (16a) after substituting therein the expression of $x_2(r)$ given in Equation (25). The Jacobian (27) is computed numerically with a quadrature scheme where the time dependent conditions in Equation (32) are tested for every discrete time $t_i$ of the quadrature scheme. The Jacobian is provided to the trust-region based nonlinear solver.

# 5 Results and discussion

Periodic responses for the system in Figure 6 are sought for the following values of the system parameters: $m_1 = m_2 = 1$, $k_1 = k_2 = 1$, and $\mu = 0.9$. The damping constants are taken as $d_1 = d_2 = 0.02$, hence the viscous damping matrix is proportional to the stiffness matrix. The system is forced harmonically only at the first degree-of-freedom, such that $f_1 = 20 \cos \omega t$ and $f_2 = 0$. This set of parameters is considered in the rest of this section unless otherwise stated.

Associated linear systems without viscous damping can be considered in the limiting cases of sliding or sticking over an entire period of the motion. In the linear case of no friction, $\mu = 0$, the friction damper is always sliding and the two free vibration undamped natural frequencies of the linear system are $\omega_1 = (\sqrt{5} - 1)/2 \approx 0.618$ for the in-phase vibration mode and $\omega_2 = (\sqrt{5} + 1)/2 \approx 1.618$ for the out-of-phase mode. In the linear case of permanent sticking for the mass $m_2$, $\mu \to \infty$, the associated linear system has a single degree-of-freedom and its free vibration undamped natural frequency is $\omega_1 = \sqrt{2}$.

## 5.1 Frequency responses and typical periodic solutions

The frequency response curves of the system are depicted in Figure 7 for values of the normal load ranging from $N = 8$ to 12, obtained by sequential continuation on the forcing frequency. Here the mean value of the total energy of the system over one period of the motion, defined as the sum of the kinetic and potential energies integrated from 0 to $T = 2\pi/\omega$ and divided by $T$, is shown on a logarithmic scale. Responses are shown from $\omega = 0.1$ to 0.8, a range which encompasses the natural frequency $\omega_1 = 0.618$ of the first vibration mode of the associated "sliding" linear undamped system. The resonance at the second vibration mode at $\omega_2 = 1.618$ is not considered here, because in this frequency range the friction damper features a rather dull behavior of permanent sliding for the considered values

<footnote n="5">⁵The fact that the partial derivatives of $\Psi_r$ are not defined everywhere, notably along the lines separating its three planes, see Figure 2(a), is not an issue since their weak version only is considered in Equation (27). Also, in Equation (32), the conditions on the right hand-side apply to both $\partial \Psi_r/\partial \dot{x}_2$ and $\partial \Psi_r/\partial r$.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 7:** Frequency response curves of the system for various values of *N*, for $N_{\phi} = 60$. The response of the linear system ($\mu = 0$) is indicated in red. The labels A, B, C, D, and E are used in later plots.</figcaption></figure>

of the normal force *N*. In **Figure 7**, the frequency response of the associated linear, viscously damped system without friction (i.e., for $\mu = 0$) is also shown for comparison. The number of harmonics in the Ritz-Galerkin approximation is $N_{\phi} = 60$, such that 30 odd harmonics ranging from the first to the 59th contribute to the response. This number of harmonics was found to be sufficient to achieve accurate converged solutions.

One observes, as expected, that as the normal force *N* increases, the effects of friction become more pronounced and the response amplitudes decrease correspondingly. For *N* = 8, 9 and 10, the frequency response features a primary resonance close to the linear system's first resonant frequency, $\omega_1 = 0.618$. As frequency decreases away from the main peak, a succession of much smaller peaks can be seen in the vicinities of $\omega = 0.54, 0.30, 0.21, 0.16, 0.13$ and 0.11 in **Figure 7**. These are believed to be super-harmonic resonances, similar to those previously observed and studied in detail for a single-degree-of-freedom system with Coulomb friction [66, 17]. These small resonances occur when one of the higher harmonics of the response, which arise due to the friction nonlinearity, drives the system at or near one of its resonance frequencies. Super-harmonic resonances in friction systems are typically associated with the apparition of additional sticking phases of the damper's motion and the greater participation of select higher harmonics of the response in the vicinity of these frequencies – hence with a qualitative change in stick-slip behavior. For the system at hand, both modes of vibration of the linearized frictionless system, namely at $\omega_1 = 0.618$ and $\omega_2 = 1.618$, engender opportunities for super-harmonic resonances. For example, the small resonance around $\omega = 0.54$ seems to correspond to a 3:1 super-harmonic resonance of the second mode (Mode 2), the one near $\omega = 0.30$ to a 5:1 Mode 2 super-harmonic resonance, the one near $\omega = 0.21$ to a 3:1 Mode 1 super-harmonic resonance, and so forth. Examination of the modal participation and harmonic content of time responses at these super-harmonic resonant frequencies confirms the presence of the relevant mode and harmonic. Note that the super-harmonic resonances involve odd harmonics only, since the Fourier series expansion of the displacements and friction force only include odd harmonics due to the symmetry of Coulomb's law. Also observe the presence of anti-resonances in the vicinity of some of the super-harmonic resonances, in most cases to the left of the peak but occasionally to the right, for instance at $\omega = 0.299$ for *N* = 10.5. Finally, observe that as the normal force increases, the primary resonant frequency decreases slightly, with the main resonant peak becoming smaller and disappearing for normal loads greater than 10.5. Correspondingly, the second mode's 3:1 super-harmonic resonant peak emerges in the vicinity of $\omega = 0.56$. It dwarfs the primary resonance peak for *N* = 10.5, and remains the only resonance at $\omega = 0.562$ for *N* = 11. For *N* = 12, resonances have all but nearly disappeared.

At the primary resonance, the friction damper is expected to be sliding back and forth during the entire period, with sticking occurring instantaneously when the mass $m_2$ changes direction, while richer stick-slip motions take place at lower frequencies in the vicinity of the super-harmonic resonances. Representative time responses of the
system are described in **Figures 8 to 11** for several points on the frequency response curves, denoted in **Figure 7** by A, B, C and D, respectively. Each figure shows the time histories over one cycle of motion of the displacement and velocity of the two degrees-of-freedom and of the friction force. As well, the friction force is represented in terms of the velocity of the second mass in order to verify whether Coulomb's classical friction law of **Figure 1** is captured accurately by the approximate Ritz-Galerkin solution. In order to ensure that the solution is converged, $N_{\phi} = 160$ Fourier harmonics are included in the solution expansion (corresponding to 80 odd harmonics, each with 80 cosinus and 80 sinus coefficients).

**Figure 8** shows the system response at point A, corresponding to the primary resonant frequency for the smallest value of normal load, $N = 8$. Both the first degree-of-freedom and the friction damper feature a harmonic-like

<figure><img src="image_12.png" /><figcaption>**Fig. 8:** Response at point A in **Figure 7**, corresponding to the primary resonance condition, for $N = 8$, $\omega = 0.618$, that is $T = 10.17$ s, and $N_{\phi} = 160$. Friction damper degree-of-freedom [solid line] and first degree-of-freedom [dashed line] are shown.</figcaption></figure>

oscillation, such that the damper undergoes one forward phase and one backward phase of sliding motion of equal durations during one period, its velocity vanishing instantaneously twice during each period at the instants of time the mass changes direction. The time history of the friction force $r(t)$ confirms the sliding-only nature of the motion, as it remains nearly constant during each forward and backward phase, undergoing a discontinuity from $r = -\mu N = -7.2$ to $r = \mu N = 7.2$ at the transition from forward to backward motion, and vice versa. In this case of slip-slip resonant response, the square wave nature of the friction force—including the discontinuities— is very well captured by the equivalent equality formulation of Coulomb's friction and associated Ritz-Galerkin solution. Finally, the plot of the friction force $r$ in terms of the damper degree-of-freedom velocity $\dot{x}_2$ in **Figure 8** shows that Coulomb's classical friction law is accurately reproduced.

Notice, in the plots of the friction force versus time and versus $\dot{x}_2$, the small oscillations just before and just after each jump between the two sliding phases. This is an indication of the classical Gibbs phenomenon, which is known to occur for Fourier series at points of discontinuity. In this case of resonant slip-slip motion, the oscillations are exacerbated because the friction force is a square wave function of time, constituting a worst case. Note that taking 80 odd harmonics into account greatly reduces the region of Gibbs oscillations, and these could be further restricted to a smaller neighborhood of the discontinuities by increasing $N_{\phi}$ further. Also, the addition of different functions, such as wavelets, in the series expansion for the solution might mitigate Gibbs oscillations [65].

**Figure 9** shows the system response at point B, also for $N = 8$, at a frequency near the 5:1 Mode 2 super-harmonic resonance. The time responses of the displacements and velocities and the friction force are richer than at the primary

<figure><img src="image_13.png" /><figcaption>**Fig. 9:** Response at point B in **Figure 7**, corresponding to a 5:1 Mode 2 super-harmonic resonance condition, for $N = 8$, $\omega = 0.293$, that is $T = 21.44$ s, and $N_{\phi} = 160$. Friction damper degree-of-freedom [solid line] and first degree-of-freedom [dashed line] are shown.</figcaption></figure>

resonance. While the damper undergoes sliding during most of the period, there are two short sticking phases, which
can be seen on the plot of *r* versus *t*. As time increases from zero to about 0.8 s, notice a short forward-sliding phase during which *r* = -7.2. The friction force then undergoes a jump, after which it rapidly increases to reach *r* = 7.2 at about *t* = 2.5 s—a short phase during which the damper is sticking, such that *x*₂ is constant and *x*₂ = 0, as can be seen in the corresponding plots. Once *r* reaches its maximum value, a long sliding phase takes place for almost half of the period until about *t* = 11.5 s, during which the damper slides in the backward direction. It is interesting to note that at the beginning of this sliding phase, say approximately from *t* = 2.5 s to 5 s, the damper appears to be almost sticking, such that its displacement *x*₂ is nearly constant and its velocity *x*₂̇ nearly zero, while it is in fact slowly drifting in the backward direction before more rapid sliding takes place.

Note the rich harmonic contents and intricate dynamics of the displacements and velocities of both degrees-of-freedom during the period, which the equality-based friction formulation is able to capture with 160 odd harmonics without any convergence issues. Also observe that in the *r* versus *x*₂̇ plot, Coulomb's classical friction law is accurately reproduced. Finally, note that the oscillations in *r* due to the Gibbs phenomenon are less significant than at the primary resonance (point A), due to the fact that *r* undergoes smaller discontinuities during the period.

Next, results for a larger normal load *N* = 10.5 in the vicinity of the 5:1 Mode 2 super-harmonic resonance condition are discussed. Figure 10 is for point C, at the anti-resonance frequency which is slightly smaller than the resonance frequency, while Figure 11 is for point D, which is at that resonant frequency. These two cases illustrate periodic solutions of the damper that are more intricate, with multiple stick-slip transitions.

<figure><img src="image_14.png" /><figcaption>**Fig. 10:** Response at point *C* in Figure 7, corresponding to the anti-resonance associated with the 5:1 Mode 2 super-harmonic resonance, for *N* = 10.5, *ω* = 0.299, that is *T* = 21.01 s, and *N*<sub>φ</sub> = 160. Friction damper degree-of-freedom [solid line] and first degree-of-freedom [dashed line] are shown.</figcaption></figure>

At the anti-resonance, it is seen in Figure 10 that two phases of sticking and sliding share an approximately equal amount of time during one period. Namely, as time increases from zero, the damper slides forward during a short time duration, after which a sticking phase takes place approximately from *t* = 2.5 s to 7.5 s, during which the friction force smoothly increases after a small discontinuity at the transition from sliding to sticking. In each phase of the motion, the solution method captures the response according to Coulomb's law, as in the sticking phases |*r*| ≤ *μN* and *x*₂ is constant and *x*₂ = 0, while in the sliding phases |*r*| = *μN*. Also the residual Gibbs oscillations are very small in this case, in agreement with the relatively small jumps in the friction force *r* at the slide-to-stick transitions.

At the super-harmonic resonant frequency, Figure 11 shows a qualitatively different response, with four sticking and four sliding phases per cycle. As time increases from zero, the damper slides until about *t* = 2 s, after which the

<figure><img src="image_15.png" /><figcaption>**Fig. 11:** Response at point *D* in Figure 7, corresponding to a 5:1 Mode 2 super-harmonic resonance condition, for *N* = 10.5, *ω* = 0.308 that is *T* = 20.4 s, and *N*<sub>φ</sub> = 160. Friction damper degree-of-freedom [solid line] and first degree-of-freedom [dashed line] are shown.</figcaption></figure>

friction force jumps and then gradually increases, corresponding to a long phase of sticking until about *t* = 4.5 s. Note
the constant value of $x_2$ and zero value of $\dot{x}_2$ during this phase. A sliding phase subsequently begins, but interestingly, it abruptly ends around $t = 5.2$ s to yield a short sticking phase until $t = 6.5$ s, after which the damper begins sliding again. Hence, in this case there are four sticking phases per cycle of the motion. Sharp variations of the friction force between the sticking and sliding states can be seen in the graph of $r$ versus $t$. Again, note that in each phase of the motion, however short in duration, the solution method captures the response according to Coulomb's law, as in the sticking phases $|r| \le \mu N$ and $x_2$ is constant and $\dot{x}_2 = 0$, while in the sliding phases $|r| = \mu N$. The residual Gibbs oscillations are slightly larger in this case, due to the more pronounced jumps in the friction force.

In conclusion, the above results show that the equality-based formulation of friction is able to capture accurately complex stick-slip damper motions, with no other approximation of Coulomb's classical law than that introduced by the finite number of Fourier harmonics included in the Ritz-Galerkin projection.

## 5.2 Convergence analysis

Since the only approximation in the procedure stems from the expansion of the unknown displacement and friction force functions in finite series of $T$-periodic functions, it is important to assess convergence in terms of the number of the basis functions. In Figure 7, it was found that increasing the number of harmonics beyond $N_\phi = 60$ led to no visible difference in the frequency responses, hence convergence was deemed to be achieved, at least regarding the system's energy. In Figures 8 to 11, while 20 harmonics (i.e., 10 odd harmonics) were generally sufficient to obtain converged displacements, accurate representation of the velocities and especially the friction force during sliding and sticking phases and the transitions between them required additional harmonics. Therefore, to ensure convergence of the solution, especially at the discontinuity points of the friction force, all these results were derived using $N_\phi = 160$ harmonics, that is, 80 odd harmonics. The nonlinear solver did not experience any convergence issues for such a large number of harmonics.

In order to further explore convergence in terms of the number of harmonics, the case of a rich time history for the friction force is considered, namely the response at point E in Figure 7. At this point, the friction force undergoes four jumps from sliding to sticking during each period and additional sharp variations, making it a good candidate to test convergence. Figure 12 depicts, for numbers of harmonics $N_\phi$ in a doubling sequence from 10 to 160, the time histories of the damper displacement, damper velocity, and friction force over one period of the motion at point E, along with the variation of $r$ as a function of $\dot{x}_2$.

Both the amplitude and the general wave shape of the displacement $x_2$ are seen to be well captured with as few as 10 harmonics, though the accurate representation of the flat areas of the time history in the sticking phases (e.g., after $t = 2$ s) requires at least 40, if not 80 harmonics. Similarly, the amplitude and general wave shape of the damper velocity $\dot{x}_2$ are well approximated for $N_\phi = 10$, but capturing the zero velocity accurately in the four sticking phases requires at least 80 harmonics.

The convergence analysis is more interesting for the friction force, since it features discontinuities and rapid variations. While the amplitude of the friction force is reasonably approximated for as little as 10 harmonics, its variation with time is not. Using less than 80 harmonics leads to poor approximations of friction force discontinuities at the transitions from sliding to sticking, as well as to spurious oscillations and thus a non-constant friction force in the sliding phases. Such oscillations are characteristic of Fourier series expansions of functions with discontinuities, which converge in a mean square sense but not pointwise.

Finally, the plots of $r$ versus $\dot{x}_2$ in Figure 12 go hand in hand with the time histories of $r(t)$ and demonstrate how adding harmonics helps capture Coulomb's classical friction law. For $N_\phi = 40$ or less, the variation of $r$ in terms of $\dot{x}_2$ features oscillations and hysteresis loops, and the discontinuity at $\dot{x}_2 = 0$ is not well approximated. Convergence is nearly achieved for $N_\phi = 80$, and an excellent representation of Coulomb's law can be seen for $N_\phi = 160$.

Regarding the convergence of the Fourier series expansion of the friction force, the worst case corresponds to the discontinuities at the transitions between sliding and sticking phases. At these jumps, the friction force varies like a square wave, and thus one expects the mean square convergence in terms of the number of harmonics to be on the order of $1/N_\phi$.

In all cases considered, the nonlinear solver converged to the solution in a few iterations.

## 5.3 Comparison with time integration

In order to validate the equality-based friction formulation with attendant Ritz-Galerkin procedure, a time-domain integration of the equations of motion (14) is directly performed. An implicit Euler integration scheme is selected, which is known to handle phase transitions of non-smooth systems efficiently. To ensure the accuracy of the numerical integration, the time step is chosen to be very small, $h = 10^{-4}$, with approximately $2 \times 10^5$ time steps per period. Again,
<figure><img src="image_16.png" /><figcaption>**Fig. 12:** Convergence analysis in terms of number of harmonics at point E in Figure 7, at the 5:1 Mode-2 super-harmonic resonance peak, for *N* = 10 and *ω* = 0.308.</figcaption></figure>

the rich response at point E in Figure 7 is considered, which features significant higher harmonics content, friction force discontinuities, and four sticking phases and eight transitions per cycle.

<figure><img src="image_17.png" /><figcaption>**Fig. 13:** Comparison of frequency and time domain solutions. Equality-based Ritz-Galerkin method with $N_{\phi} = 160$ [black] and Implicit Euler ($h = 10^{-4}$) [dashed red] at point E in Figure 7 for $N = 10$ and $\omega = 0.308$. [left] Displacements $x_1$ and $x_2$. [middle] Velocities $\dot{x}_1$ and $\dot{x}_2$. [right] Friction force $r$.</figcaption></figure>
The frequency- and time-domain solutions for the displacement and velocity of the two degrees-of-freedom and the friction force are shown in Figure 13. Observe that the two solutions are virtually undistinguishable in all cases. In particular, the higher harmonics oscillations of the displacement and velocities are well captured, as so are the constant displacement and zero velocity of the friction damper in the sticking phases of the motion. Both the jumps in the friction force at the slip-to-stick transitions and the constant friction force in the sliding phases are also accurately reproduced by both approaches.

Taken together, these results demonstrate the validity of the equality formulation of Coulomb's friction and the convergence of the weighted residual formulation with the Fourier basis functions. For a sufficient number of harmonics, the stick-slip behavior of periodic responses and Coulomb's friction law are accurately replicated, including the transitions between sliding and sticking states, and these results are in excellent agreement with those of time simulations.

Computational times for the equality-based weighted residual formulation are small. Frequency response curves in Figure 7 are obtained by sequential continuation with 700 frequency points. Each point is computed in approximately 3.5 seconds for $N_\phi = 40$ in the nonlinear solver and the recursive quadrature scheme of MATLAB. Euler time integrations, on the other hand, are slower. In order to give an estimate, the solution in Figure 13 was computed by considering the time interval [0,2000]s, required for the periodic steady-state to get established, for a CPU time of approximately 15000s using the same nonlinear solver on the same personal computer (CPU Intel(R) Core(TM) i7-2600 3.4 GHz). Also, while the computational cost of time integration might be manageable for the two-degree-of-freedom system with one friction damper in Figure 6, it can become prohibitive compared to a frequency-domain solution as the number of (linear) degrees-of-freedom increases, for example for a large-scale structure such as a bladed disk. This is because frequency-domain methods lead to algebraic equations, the nonlinear part of which can be expressed in terms of the linear degrees-of-freedom using a condensation technique, reducing the size of the system of equations to the number of nonlinear degrees-of-freedom, which for non-smooth nonlinearities is typically small. While condensation to the nonlinear degrees-of-freedom is possible in time-marching schemes, it is not commonly implemented because of the very large size of the resulting algebraic systems. Reduced-order modeling methods such as component mode synthesis can still be used to decrease the number of linear degrees-of-freedom, but the resulting system of equation will, say for a bladed disk with friction dampers, be considerably larger than the reduced set of frequency-domain nonlinear algebraic equations. This makes frequency-domain methods, such as the equality-based friction formulation presented here, both efficient and effective at obtaining the periodic solutions of non-smooth systems.

## 5.4 Analysis of the effect of damper mass

As mentioned in Section 1, in many studies the assumption of a damper with zero mass is made [19], which allows one to circumvent the difficult problem of the multi-valuedness of the friction force for zero velocity in Coulomb's classical law of friction. In Figure 6, if $m_2 = 0$ then the friction damper has no inertia, but it has flexibility since it is connected to the first degree-of-freedom (or in general, to a multi-degree-of-freedom structure) through a stiffness $k_2$. It is also commonly assumed that $d_2 = 0$, since otherwise the friction force cannot be written in terms of the relative displacement from Equation (14b) as $r = k_2(x_2 - x_1)$ (which of course is not a limitation of the formulation in this paper). This is commonly referred to as an "elastic-perfectly plastic" damper, as for $m_2 = 0$, the friction damper reduces to a spring of stiffness $k_2$ in the stuck position and it is perfectly plastic during sliding. In absolute value the friction force grows as $k_2|x_1|$ during sticking since $x_2$ is then constant, until it reaches $\mu N$ and remains at that value during sliding. Typically, the friction force exhibits a hysteresis loop in terms of the displacement of the neighboring degree-of-freedom, $x_1$. Furthermore, the friction force does not undergo jumps versus time, because for $m_2 = 0$ in Equation (14b), $r$ must be continuous since $x_1$ and $x_2$ are continuous.

Essentially, the assumption of zero damper mass allows one to reformulate the equations of motion by expressing the friction force as a piecewise linear function of the relative displacement, $x_2 - x_1$, rather than directly as a set-valued function of the damper velocity $\dot{x}_2$ per Coulomb's law (1). As a result, the nonlinearity is made smoother and is single-valued for sticking motions, which alleviates the problems associated with non-smooth systems.

While the massless damper assumption results in a more regular expression of the friction nonlinearity, it is different from the Coulomb's law regularizations analyzed in Section 2, such as arctangent, hyperbolic tangent, and piecewise linear. With these approximations, the friction force is expressed as a regular function of the friction damper velocity, while in the massless case, the friction force is a regular function of the relative displacement of the friction damper. In fact, for a massless damper, Coulomb's friction law is correctly accounted for, but indirectly so, and results capture the typical discontinuity of $r$ versus $\dot{x}_2$.

It also needs pointing out that for a massless friction damper, the velocity can be discontinuous since there is
<figure><img src="image_18.png" /><figcaption>**Fig. 14:** System periodic responses for various values of $m_2$ and $N = 6.5$, $\omega = 0.4$, and $N_{\phi} = 160$; friction damper degree-of-freedom [solid line] and first degree-of-freedom [dashed line].</figcaption></figure>

no inertia, which strictly speaking is not physically realistic. Indeed, friction in structural systems occurs between components that have both stiffness and mass. This is a limitation of the elastic-perfectly plastic massless damper model, and how much it affects the validity of the results should be assessed on a case by case basis. Also recall that jumps in the friction force versus time cannot be captured when zero mass is assumed for the friction damper.

Since the equality-based formulation in Equation (14c) is valid for any mass of the friction damper including zero (note that there is no explicit dependency on $m_2$ in the friction equality), a parametric study of the effect of friction damper mass on the periodic response was carried out in a representative case of resonance. Figure 14 shows the system response as the mass of the friction damper decreases from $m_2 = 1$ to $m_2 = 0$, for $N = 6.5$ and $\omega = 0.4$. All other system parameters remain unchanged, except that $d_2 = 0$ in the case $m_2 = 0$ in order the reproduce the results of the common massless damper assumption. Note that in this case the effect of viscous damping at the friction damper degree-of freedom was found to be negligible. Several remarks are in order.

* As the damper mass decreases, the plot of the friction force $r$ versus the displacement of the first degree-of-freedom $x_1$ shows the onset of a hysteresis loop, which for $m_2 = 0$ takes the typical simple pattern found in the literature [19], corroborating this well-known behavior. For $m_2 \neq 0$, the friction force features jumps, both in terms of time and in terms of $x_1$. As explained above, such physical discontinuities cannot be captured if a massless damper is assumed. For $m_2 = 1$, the friction force jumps between $-\mu N$ and $+\mu N$, while these discontinuities become of smaller magnitude as $m_2$ decreases, and vanish for $m_2 = 0$. For the massless damper, the friction force variation is regular, with two "smooth" sliding and sticking phases per cycle and a simple hysteresis loop,
and the more intricate and discontinuous behaviors for $m_2 \neq 0$ are not seen.

* For $m_2 = 1$, the damper undergoes back and forth slip-slip motion during the period, with no sticking phase, though there are two short time intervals of near sticking when the damper changes direction. For $m_2 = 0.5$, two short sticking phases per cycle appear, which as $m_2$ decreases to 0.25 and 0.125 become four sticking phases per cycle. The jumps in the friction force become correspondingly less prevalent as $m_2$ decreases, and again disappear for $m_2 = 0$.

* The plots $r$ versus $x_1$ all show a hysteresis loop, which as $m_2$ increases becomes less regular with more jumps, such that for $m_2 = 1$ the variation of $r$ in terms of $x_1$ is entirely discontinuous and the loop's pattern is square. Hence, the hysteresis loop is not a result of the massless damper assumption.

* Importantly, for $m_2 = 0$ there are two jumps per cycle in the damper velocity, as the damper transitions from sticking to sliding. These discontinuities, which are not physical, are of course not present in the cases for $m_2 \neq 0$.

* Finally, plots of $r$ versus $\dot{x}_2$, which are not shown here, would show that Coulomb's friction law is reproduced in all cases including for $m_2 = 0$.

## 5.5 Analysis of the piecewise linear approximation

Here the equality-based formulation is applied to the piecewise linear approximation of Coulomb's friction described in Section 3.2. The system's dynamics are governed by Equation (14a) and Equation (14b), together with the piecewise linear expression of the friction force in terms of the velocity, described by Equation (12) and governed by the slope $1/\delta$ in the approximated "sticking" phases [33]. Figure 15 shows the time histories of the displacements, velocities and friction force and the corresponding friction law $r$ versus $\dot{x}_2$, for various values of $\delta$, with a comparison to the exact solution without regularization. The system parameters are taken at point $E$ in Figure 7.

As $\delta \to 0$, the slope $1/\delta$ of the approximate friction law becomes infinite, tending to Coulomb's classical friction. While sticking is better approximated as $\delta$ decreases, such that $|r| < \mu N$ over a range of velocity smaller than $2\delta$ in absolute value, no true sticking can occur for $\delta \neq 0$ because the friction damper velocity vanishes only when $r = 0$. As $\delta$ increases, the approximation of sticking becomes poorer and its range greater, and as a result, there is less sliding (defined as $|r| = \mu N$) during one period. For example, for $\delta = 2$, sliding takes place during two short phases only, and both degree-of-freedom displacements are dominated by the first harmonic, while the velocities show relatively small higher harmonic content. Also observe the more gradual variation of the friction force. As the slope $1/\delta$ of the approximate friction law increases, higher harmonics become more prevalent, and the displacements, velocities and friction force vary more abruptly in terms of time. Sliding phases increase in duration, and for $\delta \le 0.5$, there are four sliding phases and four approximated sticking phases per cycle, such that the response is more complex than for the smoother system for $\delta = 2$. As $\delta$ decreases, convergence is first observed for $x_1$ and $x_2$ and then for $\dot{x}_1$ and $\dot{x}_2$ and $r$, which is expected. Results for $\delta = 0.5$ and $0.25$ are very similar, and in the latter case the regularized friction law $r$ versus $\dot{x}_2$ approximates Coulomb's classical friction law closely. Comparing the time histories in Figure 15 for $\delta = 0.25$ and the solution without approximation shows that the displacements, velocities and friction force are all well approximated. However, as expected, the regularized results show that in the sticking phases, the damper displacement is not quite constant and the damper velocity not quite zero, while in the slip-to-stick transitions the friction force jumps are not quite as sharp. Thus, piecewise linear regularization can yield good approximations of periodic responses, but the equality-based formulation of Coulomb's law captures true stick-slip behavior without having to resort to approximate regularizations and without any additional numerical cost.

## 6 Key features and relationship with existing approaches

The primary advantages of the equality-based weighted residual formulation presented herein for non-smooth systems with friction are its compactness, effectiveness, and efficiency. The equations of motion are simply augmented by the friction conditions expressed as non-smooth equalities involving min and max functions. This system of governing equations is cast in a weak form by approximating the unknown displacements and friction forces as truncated series of basis functions, and the time variable is removed by a Ritz-Galerkin projection onto the basis functions, yielding a nonlinear system of algebraic equations which are solved numerically using a Newton nonlinear solver. Since the friction equalities are non-smooth, the integrals derived from them in the Ritz-Galerkin procedure cannot be
<figure><img src="image_19.png" /><figcaption>**Fig. 15:** System periodic response for various values of the inverse slope parameter $\delta$ in the piecewise linear law in Equation (12), at point E in Figure 7 with $N = 10$, $\omega = 0.308$, and $N_{\phi} = 160$; friction degree-of-freedom [solid line] and first degree-of-freedom [dashed line].</figcaption></figure>

expressed analytically, hence they are evaluated numerically using a classical quadrature scheme. The computation of these integrals involves, at each sampled instant of time, the evaluation of the min and max functions of the force and velocity at the friction degree-of-freedom. Information about the sticking and sliding phases and the transitions between them is thus implicitly embedded in a weak sense in the calculation of the integrals.

In contrast, existing harmonic balance-type techniques such as the AFT, HFT, and DLFT methods require a dedicated and somewhat cumbersome treatment of the non-smooth forces, since these are written as complementarity conditions and thus are not expressible or calculable in the frequency domain. At each iteration of the nonlinear solver, an explicit transformation to the time domain is performed, where nonlinear forces are evaluated from the system's displacements and velocities and the states of the system are determined by performing tests to estimate the stick-slip transitions, after which these quantities are transformed back to the frequency domain. In the DLFT method, the passage to the time domain is particularly critical to the correction procedures for the nonlinear forces [23]. For periodic solutions, these transformations are performed using direct and inverse discrete Fourier transforms.

It is interesting to note that in the equality-based weighted residual formulation, the evaluation of the time integrals of the non-smooth functions corresponds to the time-domain tests required in the frequency-time domain methods to determine the sticking and sliding states and their transition when calculating the nonlinear forces. In the equality-based formulation, the system states and the transitions between them are implicitly determined in a weak sense through the evaluation of the integrals. The handling of the non-smooth terms is straightforward and
efficient, necessitating no dedicated treatment different from that of the equation of motions, except for the fact that the non-smooth integrals must be computed numerically. No switching back to the time domain is necessary at each iteration to calculate the nonlinear forces.

The calculation of the integrals of the nonlinear terms, Equations (24a) and (24b) using a Riemann quadrature scheme bears resemblance with taking a Discrete Fourier Transform (DFT). The residue of the Coulomb's equality condition, $\Psi_r(\dot{x}_{2h}, r_h)$, is itself calculated from the functions of time $\dot{x}_{2h}$ and $r_h$, which can be viewed as the inverse Discrete Fourier Transform (iDFT) of the vectors $j\omega x_2$ and $\mathbf{r}$ (using, for simplicity, the complex exponential expression of the Fourier Transform). Hence the process at each iteration of the nonlinear solver can be conceptually viewed as follows: starting from the harmonic representation of the friction force $\mathbf{r}$, (i) take its iDFT, (ii) express the Coulomb's equality function $\Psi_r$ in the time domain, and (iii) take its DFT to obtain the harmonic representation of the residue, which the algorithm seeks to drive to zero. Schematically, one seeks to solve the following system in the $N_\phi$ unknowns $\mathbf{r}$:

$$ \text{DFT}(\Psi_r(\text{iDFT}(j\omega x_2(\mathbf{r}), \text{iDFT}(\mathbf{r}))) = 0. \quad (33) $$

It is interesting that this process is conceptually similar to that of the AFT and HFT methods [19, 21] in which, at each iteration, an inverse DFT of the frequency-domain harmonic displacements and velocities is taken, the corresponding nonlinear forces are calculated via time integration and tests to determine state transitions, the DFT of these is taken to obtain the frequency-domain forces, and finally the updated frequency-domain displacements are solved for. The equality-based formulation of friction, though, is highly compact, applies to any representation or approximation of the Coulomb's law level-set function $\Psi$. It only requires the numerical evaluation (by quadrature or by DFT) of the integral of a non-smooth function, rather than either cumbersome time integrations in which the state of the system must be determined through tests at each time step or the use of approximate friction laws or simplified friction damper models.

Perhaps most importantly, the equality-based weighted residuals formulation is highly effective, at least in all cases reported herein, as it produces converged and accurate solutions at a low computational cost and for large numbers of harmonics. Remarkably, this is achieved without any smoothing, approximation or hypothesis regarding Coulomb's friction law. Existing frequency-time domain techniques also perform well in many cases, but they can suffer from convergence problems, especially so the AFT and the HFT, unless some kind of smoothing of the friction nonlinearity is *de facto* conducted. As covered in Section 1, in order to achieve numerical stability and produce a converged solution, in most studies either the penalization of the friction forces by introducing a finite stiffness or dissipative element or the regularization of the friction force is required. For example, the DLFT method, which is based on a dynamic Lagrangian formulation, relies on penalizing the difference between the relative interface displacements calculated in the frequency and the time domains, an algorithm which introduces artificial energy dissipation and is dependent on a penalty coefficient which influences the rate of convergence. In contrast, an important feature of the equality-based weighted residual formulation is that no approximations of penalization or regularization are made, which goes toward its simplicity and accuracy, also eliminating questionable residual penetrations and the need to select a penalty coefficient. In fact, its only approximation lies in the truncated nature of the basis of functions used in the weak formulation, a feature common to all frequency-domain methods.

Some consequences of smoothing the nonlinear forces may not be fully understood. For example, in most studies of turbomachinery, friction dampers are assumed to be massless. However, while dampers may have small mass, it may not necessarily be negligible, questioning whether the assumption makes physical sense. As seen in Section 5.4, for a massless damper the friction force is a continuous function of displacement and hence of time. This means that back-and-forth slip-slip motions (i.e., with no sticking phase) are not accurately captured, since the friction force and the acceleration are then discontinuous in time (the velocity being always continuous in time). As well, in Section 5.1 the responses which feature jumps in the friction force at the slip-to-stick transitions cannot be reproduced for a massless friction damper. Of course, the massless assumption may be appropriate, but its validity should be assessed on a case by case basis, for example through a parametric study as in Section 5.4. Regularizing the friction force has a similar effect to assuming a massless damper, as it then depends continuously of the velocity, which is a continuous function of time, hence friction force discontinuities cannot be captured. Thus regularization also prevents admissible solutions from being modeled [35].

The equality-based weighted residuals formulation does not suffer from such limitations. It is general and allows, in a natural way, for mass at frictional interfaces and for an exact Coulomb's friction law. As shown, it also applies naturally to the cases of a massless damper and regularized friction laws, if needed.

Another feature of the present approach is that it is not limited to Fourier functions for the truncated basis $\{\phi_k\}_{k=1,...,n}$ in the weighted residual procedure, as the formulation holds for general periodic functions. Functions such as periodized Chebyshev polynomials or wavelets could be used alternatively, in order to better capture discontinuities in
accelerations or friction forces. In contrast, the AFT and DLFT methods are naturally based on Fourier functions, and
while they could conceptually be generalized to other basis functions, this would require a significant redesign.

Finally, note that the present formulation allows for external forcing at the friction degree-of-freedom without any
simplifying assumption.

7 Conclusions

An equality-based weighted residuals formulation for periodic solutions of systems undergoing frictional occurrences was presented. It is highly compact, essentially embedding the non-smooth constraints in the set of governing equations such that they are handled no differently from the equations of motion, except for the numerical evaluation of the integrals arising from the weak formulation. With the method, there is no time-domain integration of nonlinear forces at each iteration and no cumbersome estimations of system states and the transitions between them. The method is tremendously versatile, as in theory it can be applied to any (non-smooth) law formulated in terms of an equality. Unlike existing methods, which require the cumbersome selection of smoothing parameters and are subject to convergence problems, its accuracy is solely dictated by the family of basis functions {$\phi_k$} in the Ritz-Galerkin expansion and the number of these functions.

The equality-based weighted residuals formulation has the potential to be a powerful method for the analysis of vibrating non-smooth systems. It is believed that equivalent equalities can be formulated for any non-smooth system whose dynamics are governed by a complementarity formulation, opening the door to the application of the equality-based weighted residuals formulation to a variety of nonlinearities, for example those with more sophisticated friction laws, two-dimensional friction, multiple frictional interfaces, and combined contact and friction.

References

[1] A. A. Ferri. Friction damping and isolation systems. *Journal of Vibration and Acoustics*, 117(B):196-206, 1995.
doi:10.1115/1.2838663.

[2] M. Mitra and B. Epureanu. Dynamic modeling and projection-based reduction methods for bladed disks with nonlinear frictional and intermittent contact interfaces. *Applied Mechanics Reviews*, 71(5), 2019. doi:10.1115/1.4043083.

[3] P. Flores, R. Leine, and C. Glocker. Modeling and analysis of rigid multibody systems with translational clearance joints based on the nonsmooth dynamics approach. In K. Arczewski, W. Blajer, J. Fraczek, and M. Wojtyra, editors, *Multibody Dynamics: Computational Methods and Applications*, pages 107–130. Springer, 2011. doi:10.1007/978-90-481-9971-6_6.

[4] S. W. Shaw. On the dynamic response of a system with dry friction. *Journal of Sound and Vibration*, 108(2):305-325, 1986.
doi:10.1016/s0022-460x(86)80058-x.

[5] J. P. Den Hartog. Forced vibrations with combined Coulomb and viscous friction. *Journal of Fluids Engineering*, 53(2):107-115, 1931. doi:10.1115/1.4022656.

[6] V. Acary and B. Brogliato. *Numerical Methods for Nonsmooth Dynamical Systems*. Springer, 2008. doi:10.1007/978-3-540-75392-6.

[7] T. A. Laursen and J. C. Simo. A continuum-based finite element formulation for the implicit solution of multibody, large deformation-frictional contact problems. *International Journal for Numerical Methods in Engineering*, 36(20):3451–3485, 1993. doi:10.1002/nme.1620362005.

[8] J.-H. Heegaard and A. Curnier. An augmented Lagrangian method for discrete large-slip contact problems. *International Journal for Numerical Methods in Engineering*, 36(4):569–593, 1993. doi:10.1002/nme.1620360403.

[9] A. Le van and T. T. H. Nguyen. A weighted residual relationship for the contact problem with Coulomb friction. *Computers & Structures*, 87(23-24):1580–1601, 2009. doi:10.1016/j.compstruc.2009.08.013.

[10] L. Charroyer, O. Chiello, and J.-J. Sinou. Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method. *International Journal of Mechanical Sciences*, 144:90–101, 2018. doi:10.1016/j.ijmecsci.2018.05.045.

[11] M. Krack and J. Gross. *Harmonic Balance for Nonlinear Vibration Problems*. Springer, 2019. doi:10.1007/978-3-030-14023-6.

[12] J. H. Griffin. Friction damping of resonant stresses in gas turbine engine airfoils. *Journal of Engineering for Power*, 102(2): 329-333, 1980. doi:10.1115/1.3230256.
[13] A. A. Ferri and E. H. Dowell. The behavior of a linear, damped modal system with a non-linear spring-mass dry friction damper system attached, part II. *Journal of Sound and Vibration*, 101(1):55-74, 1985. doi:10.1016/s0022-460x(85)80038-9.

[14] A. Sinha and J. H. Griffin. Friction damping of flutter in gas turbine engine airfoils. *Journal of Aircraft*, 20(4):372-376, 1983. doi:10.2514/3.44878.

[15] L. Jézéquel. Structural damping by slip in joints. *Journal of Vibration and Acoustics*, 105(4):497-504, 1983. doi:10.1115/1.3269134.

[16] Z. Yan, H. Dai, Q. Wang, and S. N. Atluri. Harmonic balance methods: A review and recent developments. *Computer Modeling in Engineering & Sciences*, 137(2):1419-1459, 2023. doi:10.32604/cmes.2023.028198.

[17] C. Pierre, A. A. Ferri, and E. H. Dowell. Multi-harmonic analysis of dry friction damped systems using an incremental harmonic balance method. *Journal of Applied Mechanics*, 52(4):958-964, 1985. doi:10.1115/1.3169175.

[18] S. L. Lau, Y. K. Cheung, and S. Y. Wu. Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems. *Journal of Applied Mechanics*, 50(4a):871-876, 1983. doi:10.1115/1.3167160.

[19] T. M. Cameron and J. H. Griffin. An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems. *Journal of Applied Mechanics*, 56(1):149-154, 1989. doi:10.1115/1.3176036.

[20] L. Woiwode, N. N. Balaji, J. Kappauf, F. Tubita, L. Guillot, C. Vergez, B. Cochelin, A. Grolet, and M. Krack. Comparison of two algorithms for harmonic balance and path continuation. *Mechanical Systems and Signal Processing*, 136:106503, 2020. doi:10.1016/j.ymssp.2019.106503.

[21] O. Poudou and C. Pierre. Hybrid frequency-time domain methods for the analysis of complex structural systems with dry friction damping. In *44th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference*. AIAA, 2003. doi:10.2514/6.2003-1411.

[22] J. Guillen, C. Pierre, and T. Lagrange. An advanced damper model for the dynamics of dry friction damped systems. In *Volume 7B: 17th Biennial Conference on Mechanical Vibration and Noise*. ASME, 1999. doi:10.1115/detc99/vib-8083.

[23] S. Nacivet, C. Pierre, F. Thouverez, and L. Jézéquel. A dynamic Lagrangian frequency-time method for the vibration of dry-friction-damped systems. *Journal of Sound and Vibration*, 265(1):201-219, 2003. doi:10.1016/s0022-460x(02)01447-5.

[24] S. Quaegebeur, B. Chouvion, and F. Thouverez. Impact of mistuned underplatform dampers on the nonlinear vibration of bladed disks. *Journal of Engineering for Gas Turbines and Power*, 143(12), 2021. doi:10.1115/1.4051868.

[25] T. Vadcard, A. Batailly, and F. Thouverez. On Harmonic Balance Method-based Lagrangian contact formulations for vibro-impact problems. *Journal of Sound and Vibration*, 531:116950, 2022. doi:10.1016/j.jsv.2022.116950.

[26] D. Charleux. *Étude des effets de la friction en pied d'aube sur la dynamique des roues aubagées*. PhD thesis, Ecole Centrale de Lyon, 2006, URL https://theses.hal.science/tel-01130569.

[27] D. Laxalde, L. Salles, L. Blanc, and F. Thouverez. Non-linear modal analysis for bladed disks with friction contact interfaces. In *Turbo Expo: Power for Land, Sea, and Air: Structures and Dynamics, Parts A and B*, volume 5, pages 457-467. ASME, 2008. doi:10.1115/gt2008-50860.

[28] C. Joannin, B. Chouvion, F. Thouverez, M. Mbaye, and J.-P. Ousty. Nonlinear modal analysis of mistuned periodic structures subjected to dry friction. *Journal of Engineering for Gas Turbines and Power*, 138(7), 2016. doi:10.1115/1.4031886.

[29] S. Quaegebeur, B. Chouvion, and F. Thouverez. Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties. *Mechanical Systems and Signal Processing*, 145:106970, 2020. doi:10.1016/j.ymssp.2020.106970.

[30] S. He, S. Cho, and R. Singh. Prediction of dynamic friction forces in spur gears using alternate sliding friction formulations. *Journal of Sound and Vibration*, 309(3-5):843-851, 2008. doi:10.1016/j.jsv.2007.06.077.

[31] A. A. Ferri and E. H. Dowell. Frequency domain solutions to multi-degree-of-freedom, dry friction damped systems. *Journal of Sound and Vibration*, 124(2):207-224, 1988. doi:10.1016/s0022-460x(88)80183-4.

[32] A. Cardona, T. Coune, A. Lerusse, and M. Géradin. A multiharmonic method for non-linear vibration analysis. *International Journal for Numerical Methods in Engineering*, 37(9):1593-1608, 1994. doi:10.1002/nme.1620370911.

[33] M. Berthillier, C. Dupont, R. Mondal, and J.-J. Barrau. Blades forced response analysis with friction dampers. *Journal of Vibration and Acoustics*, 120(2):468-474, 1998. doi:10.1115/1.2893853.

[34] D. Karnopp. Computer simulation of stick-slip friction in mechanical dynamic systems. *Journal of Dynamic Systems, Measurement, and Control*, 107(1):100-103, 1985. doi:10.1115/1.3140698.
[35] D. D. Quinn. A new regularization of Coulomb friction. *Journal of Vibration and Acoustics*, 126(3):391-397, 2004.
doi:10.1115/1.1760564.

[36] P. Vigué, C. Vergez, S. Karkar, and B. Cochelin. Regularized friction and continuation: Comparison with Coulomb's law. *Journal of Sound and Vibration*, 389:350-363, 2017. doi:10.1016/j.jsv.2016.11.002.

[37] D. Schüthe, F. Wenk, and U. Frese. Dynamics calibration of a redundant flexible joint robot based on gyroscopes and encoders. In *13th International Conference on Informatics in Control, Automation and Robotics*. Science and Technology Publications, 2016. doi:10.5220/0005976603350346.

[38] C. Westin and R. Irani. Continuously differentiable stick-slip friction model with applications to cable simulation using nonlinear finite elements. In *2020 IEEE Conference on Control Technology and Applications (CCTA)*. IEEE, 2020. doi:10.1109/ccta41146.2020.9206303.

[39] S. Baek and B. Epureanu. Contact model identification for friction ring dampers in blisks with reduced order modeling. *International Journal of Non-Linear Mechanics*, 120:103374, 2020. doi:10.1016/j.ijnonlinmec.2019.103374.

[40] W. Tang and B. Epureanu. Nonlinear dynamics of mistuned bladed disks with ring dampers. *International Journal of Non-Linear Mechanics*, 97:30-40, 2017. doi:10.1016/j.ijnonlinmec.2017.08.001.

[41] J. J. Chen, B. D. Yang, and C.-H. Menq. Periodic forced response of structures having three-dimensional frictional constraints. *Journal of Sound and Vibration*, 229(4):775-792, 2000. doi:10.1006/jsvi.1999.2397.

[42] E. P. Petrov and D. J. Ewins. Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks. *Journal of Turbomachinery*, 125(2):364-371, 2003. doi:10.1115/1.1539868.

[43] E. P. Petrov and D. J. Ewins. Effects of damping and varying contact area at blade-disk joints in forced response analysis of bladed disk assemblies. *Journal of Turbomachinery*, 128(2):403-410, 2005. doi:10.1115/1.2181998.

[44] J. J. Chen and C.-H. Menq. Periodic response of blades having three-dimensional nonlinear shroud constraints. *Journal of Engineering for Gas Turbines and Power*, 123(4):901-909, 2001. doi:10.1115/1.1385828.

[45] D. Laxalde and F. Thouverez. Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces. *Journal of Sound and Vibration*, 322(4-5):1009-1025, 2009. doi:10.1016/j.jsv.2008.11.044.

[46] S. Zucca and C. M. Firrone. Nonlinear dynamics of mechanical systems with friction contacts: Coupled static and dynamic multi-harmonic balance method and multiple solutions. *Journal of Sound and Vibration*, 333(3):916-926, 2014. doi:10.1016/j.jsv.2013.09.032.

[47] D. SÜß and K. Willner. Investigation of a jointed friction oscillator using the multiharmonic balance method. *Mechanical Systems and Signal Processing*, 52-53:73-87, 2015. doi:10.1016/j.ymssp.2014.08.003.

[48] J. H. Wang and W. K. Chen. Investigation of the vibration of a blade with friction damper by HBM. *Journal of Engineering for Gas Turbines and Power*, 115(2):294-299, 1993. doi:10.1115/1.2906708.

[49] M. Mitra, S. Zucca, and B. Epureanu. Adaptive microslip projection for reduction of frictional and contact nonlinearities in shrouded blisks. *Journal of Computational and Nonlinear Dynamics*, 11(4):1-15, 2016. doi:10.1115/1.4033003.

[50] T. Liu, D. Zhang, and Y. Xie. A nonlinear vibration analysis of forced response for a bladed-disk with dry friction dampers. *Journal of Low Frequency Noise, Vibration and Active Control*, 38(3-4):1522-1539, 2019. doi:10.1177/1461348419834759.

[51] P. R. Dahl. Solid friction damping of mechanical vibrations. *AIAA Journal*, 14(12):1675-1682, 1976. doi:10.2514/3.61511.

[52] C.-H. Menq, J. Bielak, and J. H. Griffin. The influence of microslip on vibratory response, part i: A new microslip model. *Journal of Sound and Vibration*, 107(2):279-293, 1986. doi:10.1016/0022-460x(86)90238-5.

[53] A. A. Ferri and B. S. Heck. Vibration analysis of dry friction damped turbine blades using singular perturbation theory. *Journal of Vibration and Acoustics*, 120(2):588-595, 1998. doi:10.1115/1.2893868.

[54] S. Quaegebeur, B. Chouvion, and F. Thouverez. Nonlinear dynamic analysis of three-dimensional bladed-disks with frictional contact interfaces based on cyclic reduction strategies. *International Journal of Solids and Structures*, 236-237:111277, 2022. doi:10.1016/j.ijsolstr.2021.111277.

[55] J.-J. Moreau. An introduction to unilateral dynamics. In M. Frémond and F. Maceri, editors, *Novel Approaches in Civil Engineering*. Springer, 2004. doi:10.1007/978-3-540-45287-4.

[56] P. Alart and A. Curnier. A mixed formulation for frictional contact problems prone to Newton like solution methods. *Computer Methods in Applied Mechanics and Engineering*, 92(3):353-375, 1991. doi:10.1016/0045-7825(91)90022-x.
[57] M. Jean. The non-smooth contact dynamics method. *Computer Methods in Applied Mechanics and Engineering*, 177(3-4): 235-257, 1999. doi:10.1016/s0045-7825(98)00383-1.

[58] L. Xuewen, A.-K. Soh, and C. Wanji. A new non-smooth model for three dimensional frictional contact problems. *Computational Mechanics*, 26(6):528-535, 2000. doi:10.1007/s004660000202.

[59] S. Hüeber, G. Stadler, and B. I. Wohlmuth. A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction. *SIAM Journal on Scientific Computing*, 30(2):572-596, 2008. doi:10.1137/060671061.

[60] G. Stadler. Semismooth Newton and augmented Lagrangian methods for a simplified friction problem. *SIAM Journal on Optimization*, 15(1):39-62, 2004. doi:10.1137/s1052623403420833.

[61] V. Acary, M. Brémond, and O. Huber. On solving contact problems with Coulomb friction: formulations and numerical comparisons. In R. Leine, V. Acary, and O. Brüls, editors, *Advanced Topics in Nonsmooth Dynamics: Transactions of the European Network for Nonsmooth Dynamics*, pages 375-457. Springer, 2018. doi:10.1007/978-3-319-75972-2_10.

[62] T. Lu and M. Legrand. Nonsmooth modal analysis via the boundary element method for one-dimensional bar systems. *Nonlinear Dynamics*, 107(1):227-246, 2022. doi:10.1007/s11071-021-06994-z.

[63] A. Y. T. Leung, C. Guoqing, and C. Wanji. Smoothing Newton method for solving two- and three-dimensional frictional contact problems. *International Journal for Numerical Methods in Engineering*, 41(6):1001-1027, 1998. doi:10.1002/(SICI)1097-0207(19980330)41:6<1001::AID-NME319>3.0.CO;2-A.

[64] R. Leine and H. Nijmeijer. *Dynamics and Bifurcations of Non-Smooth Mechanical Systems*, volume 18 of *Lecture Notes in Applied and Computational Mechanics*. Springer, 2013. doi:10.1007/978-3-540-44398-8.

[65] S. Jones and M. Legrand. 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):351-374, 2014. doi:10.1002/nme.4807.

[66] G. Csernák, G. Stépán, and S. W. Shaw. Sub-harmonic resonant solutions of a harmonically excited dry friction oscillator. *Nonlinear Dynamics*, 50(1-2):93-109, 2007. doi:10.1007/s11071-006-9145-6.

[67] M. Legrand and C. Pierre. Equal[iseR]: A frequency-domain solver for the vibration analysis of structural mechanical systems subject to unilateral and frictional contact occurrences. 2024, URL https://archive.softwareheritage.org/swh:1:dir:925849db7fa38d37336bc85c4d56ed8b7b2ed2b7.

## Authors' contributions

M.L. and C.P. contributed to the design of the research and to the analysis of the results. M.L. implemented the developed idea and corresponding scripts while C.P. wrote the manuscript. Both authors approved the final version of the work.

## Acknowledgements

The first author gratefully acknowledges the support of the Natural Sciences and Engineering Research Council of Canada through the Discovery Grant program.

## Supplementary material

The numerical code used in this work is shared publicly on a repository [67]. Please cite this paper and the reference to the numerical code for any subsequent publication exploiting it, and comply to the terms of license published on the repository.