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

HAL
open science

# An Equality-Based Weighted Residual Formulation for the Vibration of Systems with Two-Dimensional Friction

Arash Hashemi, Christophe Pierre, Mathias Legrand

► To cite this version:

Arash Hashemi, Christophe Pierre, Mathias Legrand. An Equality-Based Weighted Residual Formulation for the Vibration of Systems with Two-Dimensional Friction. Journal of Structural Dynamics, 2025, 3, pp.127-155. 10.25518/2684-6500.282. hal-04852967v2

HAL Id: hal-04852967

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

Submitted on 14 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# An Equality-Based Weighted Residual Formulation for the Vibration of Systems with Two-Dimensional Friction

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

## Abstract

An equality-based weighted residual formulation is proposed for the periodic responses of vibrating systems subject to two-dimensional dry friction on a plane. Coulomb's law is expressed as two coupled nonsmooth equality conditions which augment the equations of motion, resulting in a mixed displacement-friction force formulation whose periodic solutions are sought using a standard Ritz-Galerkin procedure in time. The shape functions considered are the classical Fourier functions, and a quasi-analytical expression for the Jacobian of the friction terms is derived in a piecewise linear fashion and computed in a weighted residual sense. The method is based on an exact equality representation of Coulomb's law for interfaces with mass, thus avoiding common hypotheses such as regularization, penalization, or massless interfaces. It is entirely carried out in the frequency domain, contrary to existing frequency-time methods which require the calculation of contact forces in the time domain at each iteration of the nonlinear solver. The method is compact and found to be robust and accurate. It is void of convergence or other numerical issues up to very large numbers of harmonics of the response in all cases considered. Periodic responses featuring complex two-dimensional interface motions and multiple stick-slip transitions are calculated accurately at various resonant and sub-resonant excitation frequencies, at a reasonable computational cost. Since the only approximation in the procedure is the finite number of terms in the Ritz-Galerkin expansion, the intricate behavior of the two-dimensional friction force dictated by Coulomb's law can be captured with a high degree of accuracy.

**Keywords:** Non-smooth dynamics; Vibration analysis; Two-dimensional dry friction; Periodic solutions

Received on December 21, 2024, Accepted on July 11, 2025, Published on August 27, 2025

## 1 Introduction

Nonsmooth behavior can be observed across different areas of engineering. In electrical engineering, circuits consisting of nonlinear resistors, switches, and ideal diodes feature nonsmooth current-voltage characteristics [1]. In mechanical engineering, unilateral contact and frictional occurrences are the primary sources of nonsmooth characteristics. Examples of nonsmooth mechanical structures include robot manipulators and cable-driven parallel robots with friction at joints and cable-puley interfaces [2], and railroad car bogies with wedge dampers to reduce force transmission from tracks [3]. Coulomb friction can also be used to improve simulations of contacting objects (hair, cloth, grains...) in computer graphics [4, 5, 6]. It also forms an important class of problems in control theory [7, 8]. In the area of structural dynamics of turbomachinery rotors, blades experience a number of contact and frictional events, including at blade-disk and blade-casing interfaces, shrouds and snubbers, and dry friction dampers that are

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

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

<footnote>*ahashemi1@stevens.edu</footnote>

<footnote>†cpierre@stevens.edu</footnote>

<footnote>‡mathias.legrand@mcgill.ca</footnote>
commonly used to mitigate nefarious vibration levels [9]. These numerous engineering applications require effective and efficient modeling and analysis tools, spurring the attention of many researchers over the years.

Coulomb's classical law is the most common friction model for structures, and a number of more advanced laws, such as Stribeck's model, have been proposed over the years [10, 11]. The equations of motion of systems with friction can be challenging to solve, because the Coulomb conditions feature a multi-valued, or set-valued force-velocity constitutive law, which must be expressed in terms of both equalities and inequalities [12] in order to capture discontinuities.

Building on Den Hartog's seminal work [13], the vibration of dry friction-damped systems has been extensively investigated using two principal approaches: time-domain and frequency-domain methods. Time-domain methods are best suited for simulations of transients and rigid body and quasi-static frictional contact, relying on advanced time-stepping or event-driven schemes to detect state transitions (e.g., from sliding to sticking) [14, 15]. For finite element models of deformable bodies in frictional contact, the Augmented Lagrangian method is favored, where multipliers act as contact forces and penetration between contacting surfaces is penalized [16, 17]. However, time-domain methods can be computationally intensive for small damping and thus are less suitable for vibration or parametric analyses. Shooting methods [18] can be effective for periodic solutions but are highly sensitive to initial guess and computationally expensive for many degrees of freedom (DOF). Frequency-domain methods, by contrast, efficiently provide steady-state solutions and frequency responses through continuation methods. Additionally, linear DOFs can be dynamically condensed to the nonlinear ones, dramatically reducing the problem size for large structures. These advantages make frequency-domain methods the preferred choice for analyzing nonsmooth system vibrations, with the harmonic balance (HB) method being the most widely used for periodic solutions [19].

For one-dimensional (1D) friction, early applications of the HB method were limited to a single harmonic [20, 21, 22, 23]. However, periodic responses inherently involve multi-harmonic content due to the occurrence of sticking and sliding phases. Since the exact instants of state transition cannot be determined analytically, the friction force cannot be expressed in a closed form as a function of relative velocity. A widely used approach to address this challenge is to regularize Coulomb's law into a continuous single-valued function of velocity. This smoothening enables the classical use of HB with multiple harmonics, yielding a system of nonlinear algebraic equations in the Fourier coefficients. Common regularizations or approximations of 1D friction laws include the hyperbolic tangent function [24, 25, 26], arctangent function [27], signum function [28], sigmoid function in robotics [29, 30], piecewise linear function [31, 32, 33], and various adaptations [34, 35, 36, 37]. Drawbacks of regularization are the absence of true sticking, as the friction force vanishes only at zero velocity, and the lack of a priori knowledge of its validity in terms of the regularized function's parameter. Note that regularized or approximate laws have rarely been developed for two-dimensional (2D) friction.

For 1D frictional occurrences, mixed frequency-time domain methods, such as Alternating Frequency/Time (AFT) [38] and Hybrid Frequency-Time (HFT) [39, 40], address the challenge of expressing nonsmooth contact forces in the frequency domain. These methods rely on a Fast Fourier Transform (FFT) to alternate between the frequency and time domains, solving the dynamics in the former and integrating the contact forces in the latter. By assuming massless interfaces, friction forces are simplified into continuous piecewise linear functions of relative interface displacements, preserving Coulomb's law while improving convergence [37, 41]. The Dynamic Lagrangian Frequency/Time (DLFT) method improves upon this by accounting for interfaces with mass and enforcing non-penetration via penalized differences between the frequency-domain displacements and those obtained from the time-integrated contact forces [42]. A significant limitation of mixed frequency-time methods is the time integration required to compute contact forces at each iteration of the nonlinear solver. This requires precise criteria to detect stick-slip transitions at each interface point during one or more periods of motion, which can be computationally challenging and affect convergence.

Two-dimensional frictional contact, such as that of a point mass on a surface, presents additional complexities compared to 1D friction. The direction of the friction force becomes an unknown, and interface points may undergo intricate 2D trajectories on the contact surface. Common assumptions about the sliding angle include elliptical trajectories of the contact point [43, 44, 45] or the decoupling of friction into two orthogonal 1D models [46, 47]. Early studies on modeling 2D friction and 3D flexible contact interfaces were confined to single-harmonic solutions [48, 49, 44, 50, 43]. Later research showed that mixed frequency-time methods, such as the AFT and DLFT, are also applicable to 2D friction and can compute multi-harmonic solutions [51, 48, 52, 53, 54, 55]. In most of these studies, turbomachinery bladed disks are the primary engineering application. However, 2D friction exacerbates the limitations of mixed methods, as the switching between frequency and time domains requires time integration over multiple periods at each nonlinear iteration, until the 2D trajectory of each interface point is compatible with the contact forces [56, 57]. Implementation can be challenging, with convergence issues especially for multiple contact points, necessitating calibration or validation from time-domain simulations [48, 51, 9]. Furthermore, the original DLFT
formulation based on velocity may encounter convergence issues in the presence of a static preload, requiring a displacement-based alternative [52].

In both 1D and 2D friction studies, a flexible massless interface is often assumed to represent the rough and wavy nature of contact surfaces, though identifying accurate contact stiffnesses remains a challenge; for instance, see the models in [58, 59, 60]. In studies of bladed disks, finite element models frequently employ massless contact elements [61, 48]. This approach, especially for 2D friction, helps mitigate convergence issues arising from nonsmooth behavior. Notably, the AFT is not readily applicable when masses are in direct contact, leaving the DLFT as the only frequency-time domain method capable of handling interfaces with mass [62]. However, friction physically occurs between bodies with mass, and the massless assumption can lead to unrealistic jumps in the relative velocity of the frictional interface. For friction dampers, it has also been shown that damper mass can significantly influence blade response [63].

In addition to bladed disk assemblies, wedge dampers are commonly used in railroad car bogies due to their manufacturing simplicity and effectiveness. Two approaches were presented in [11] for 2D friction dampers: a decoupled model which essentially combines two independent 1D models, and a coupled model which accounts for the resultant friction force on the contact surface. Two-dimensional stick-slip motion in train bogies was investigated by time integration in [64] using a hyperbolic secant function to approximate the Stribeck model. A friction direction angle was introduced, determined by the relative velocities during sliding or the applied forces during sticking. In a subsequent study [65], the response of a three-piece freight truck with two wedge dampers in contact with a bolster and side walls was analyzed, capturing 2D frictional interactions as well as separation.

Very few studies of 2D friction made use of regularization of the friction law. The signum function was used to simulate brake and suspension systems for trains [66] and to analyze slipping behavior in piping systems [47]. The hyperbolic tangent function was used to analyze transitions between stick and slip in systems with two-dimensional friction [67]. This regularization method can be highly accurate if the slope-controlling parameter is large enough, with the caveat that the formulation may have convergence issues.

The extensive literature on numerical simulation tools for nonsmooth dynamical systems includes advanced time-stepping and event-tracking schemes for mechanics, electronics, and control systems [68, 14, 69, 70]. In mechanics, the focus has been on unilateral contact with Coulomb friction for rigid or quasi-static deformable bodies. An interesting approach, referred to as the *implicit formulation*, is to convert the complementarity and inequality conditions of unilateral contact and friction into equivalent nonsmooth equalities, expressed as the zero-level set of functions [71, 72, 73, 74, 75]. The resulting formulation is simpler and enables an effective solution method for the three-dimensional frictional contact problem. Several variants of implicit formulations of friction have been developed in contact mechanics [76, 77, 78], including some which introduce an additional variable such as the slip angle at the frictional interface [79]. The implicit formulation was recently adopted successfully for the periodic autonomous responses and modal analysis of a one-dimensional bar with Signorini contact [80], and as well for the periodic forced responses of systems with 1D Coulomb friction [81], demonstrating its high potential for structural dynamics and vibrations problems.

In the present paper, the equality-based weighted residual formulation developed by [81] for the periodic solutions of forced systems with 1D friction is generalized to 2D frictional occurrences. Namely, using a projector, the 2D friction law is expressed as two nonsmooth implicit equalities, which augment the equations governing the vibration response. This yields a primal-dual displacement-force formulation, the governing equations of which are then collectively satisfied in a weak integral sense through a weighted residual formulation over one period of the motion. In this formulation, the direction of sliding motion is an inherent unknown embedded in the two friction equalities. The new method 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. Instead, a quadrature scheme based on a Fast Fourier Transform is implemented to compute the integrals of nonsmooth terms in the weighted residual procedure. In Section 2, various equivalent equality-based formulations for 2D friction, which are expressed as zero-level set functions, are presented for a point mass subject to friction on a surface. These nonsmooth friction conditions are discussed, and the one found most suitable for computational implementation is selected. In Section 3, a Ritz-Galerkin procedure is applied to the equations of motion and the implicit friction conditions to obtain the steady-state periodic solutions of a simple 2-DOF system with 2D friction subject to periodic external excitation. The weak formulation of the Jacobian of the nonlinear system is derived in a piecewise linear fashion to accelerate the numerical computations. Results are presented in Section 4, including frequency response curves for various normal force values and time histories of system displacements, velocities, and friction force. Some of these periodic responses exhibit complex trajectories of the friction point on the plane of contact, with multiple stick-slip phases, necessitating a very large number of harmonics to capture friction forces accurately. A convergence study is performed up to large numbers of harmonics and the frequency-domain results are validated against numerical time integration. To demonstrate
the potential applicability of the new formulation to large-scale structural systems, it is generalized in Section 5 to a two-mass system consisting of a linear component with an attached friction damper component. At the end, a study of the effect of damper mass is presented.

## 2 Equality-Based Formulations of Two-Dimensional Coulomb Friction

Consider the mass depicted in Figure 1, subject to 2D friction on the contact surface engendered by the unit vectors $i_1$ and $i_2$. It is subject to a constant normal force of magnitude $N$ applied in the negative $i_3$ direction, and as a result

<figure><img src="image_4.png" /><figcaption>Fig. 1: Conceptual point mass subject to two-dimensional Coulomb's friction.</figcaption></figure>

experiences tangential frictional forces in the plane of contact. Defining the displacement vector as $\mathbf{u} = [u_1, u_2]^T$ and the tangential friction force vector as $\mathbf{r} = [r_1, r_2]^T$, Coulomb's classical law is expressed by the set of equalities and inequalities in the velocity and friction force vectors:

$$ \left\{ \begin{array}{ll} \|\dot{\mathbf{u}}\| = 0 & \Rightarrow \|\mathbf{r}\| \le \mu N \quad [\text{sticking}] \\ \|\dot{\mathbf{u}}\| \neq 0 & \Rightarrow \mathbf{r} = -\mu N \frac{\dot{\mathbf{u}}}{\|\dot{\mathbf{u}}\|} \quad [\text{sliding}] \end{array} \right. \qquad (1) $$

where $\mu$ is the coefficient of friction. This shows that the mapping between $\mathbf{r}$ and $\dot{\mathbf{u}}$ is not single-valued but set-valued, since the set of admissible values of $\mathbf{r}$ when $\dot{\mathbf{u}} = 0$ does not reduce to a single point. In this study, the normal load $N$ is assumed constant, and Coulomb's model reduces to Tresca's friction model.

### 2.1 Equality-Based Two-Dimensional Friction Conditions

Among several equality-based formulations into which Equation (1) can be recast that have been proposed in the nonsmooth mechanics literature, the one developed in [76] was found to be the most suitable for calculating the periodic solutions of systems subject to 2D friction. In two dimensions, this vector equality takes the form

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

where $\rho$ is a scaling parameter of default value 1. The solution of Equation (2) is expressed as the zero-level set of the function $\Psi$, which can be shown to be identical to the admissible set of points $(\mathbf{r}, \dot{\mathbf{u}})$ of the classical Coulomb formulation in Equation (1). The proof for this is derived by considering the two cases for the max function in Equation (2), as follows.

* **Case** $\|\mathbf{r} - \rho \dot{\mathbf{u}}\| \le \mu N$: Equation (2) readily reduces to $\mu N \rho \dot{\mathbf{u}} = 0$, hence $\dot{\mathbf{u}} = 0$, which in turn yields the sticking condition $\|\mathbf{r}\| \le \mu N$ in Equation (1).

* **Case** $\mu N < \|\mathbf{r} - \rho \dot{\mathbf{u}}\|$: Equation (2) becomes $\mathbf{r}(\|\mathbf{r} - \rho \dot{\mathbf{u}}\| - \mu N) + \mu N \rho \dot{\mathbf{u}} = 0$, from which one infers $\dot{\mathbf{u}} = -\alpha \mathbf{r}$, where $\alpha > 0$. One can then write $\|\mathbf{r} - \rho \dot{\mathbf{u}}\| = (1 + \rho \alpha)\|\mathbf{r}\|$, which, upon substitution into the previous identity, yields $(1 + \rho \alpha)\|\mathbf{r}\| = (1 + \rho \alpha)\mu N$ and thus implies $\|\mathbf{r}\| = \mu N$. Since the norm of the friction force vector is $\mu N$ and it is collinear and opposite in direction to the velocity vector, the friction force vector can be written as $\mathbf{r} = -\mu N \dot{\mathbf{u}}/\|\dot{\mathbf{u}}\|$.
Note that the case $\dot{\mathbf{u}} = \mathbf{0}$ is not possible, since it would yield $\mu N < \|\mathbf{r}\|$, which violates the calculated norm $\mu N$ of the friction force vector. Thus, the sliding condition $\|\dot{\mathbf{u}}\| \neq 0$ implies $\mathbf{r} = -\mu N \dot{\mathbf{u}} / \|\dot{\mathbf{u}}\|$, and Equation (1) is recovered.

This proves that the equality formulation (2) is equivalent to Coulomb's classical law (1). It can be used as a 2D nonsmooth friction condition which augments the equations of motion. This formulation was found to be an excellent choice for the weighted residual scheme, as no convergence issues were observed.

## 2.2 Alternate Implicit Formulations

Several alternative equality-based formulations have been proposed in the field of nonsmooth mechanics to represent Coulomb's 2D law, and some of these are described below to illustrate the versatility of the approach. While these formulations were tested, they exhibited inferior convergence compared to Equation (2) and were therefore not selected.

As a first example, the projective version of 2D friction reads [82]

$$ \mathbf{r} - \operatorname{proj}_{D(\mu N)}(\mathbf{r} - \rho \dot{\mathbf{u}}) = \mathbf{0} \qquad (3) $$

where $\operatorname{proj}_D(\mathbf{x})$ denotes the projection of vector $\mathbf{x}$ on the domain $D$. The domain $D(\mu N)$ is Tresca's friction disk (the cross section of Coulomb's cone) defined as $D(\mu N) = \{\mathbf{r} : \|\mathbf{r}\| \le \mu N\}$ [71] and the projection operator has the explicit form

$$ \operatorname{proj}_{D(\mu N)}(\mathbf{x}) = \begin{cases} \mathbf{x} & \text{if } \mathbf{x} \in D(\mu N), \\ \mu N \mathbf{x} / \|\mathbf{x}\| & \text{if } \mathbf{x} \notin D(\mu N) \end{cases} \qquad (4) $$

with $\mathbf{x} \in D(\mu N) \Leftrightarrow \|\mathbf{x}\| \le \mu N$ and $\mathbf{x} \notin D(\mu N) \Leftrightarrow \|\mathbf{x}\| > \mu N$. Equation (3) thus becomes

$$ \mathbf{r} - \begin{cases} \mathbf{r} - \rho \dot{\mathbf{u}} & \text{if } \|\mathbf{r} - \rho \dot{\mathbf{u}}\| \le \mu N \\ \frac{\mu N(\mathbf{r} - \rho \dot{\mathbf{u}})}{\|\mathbf{r} - \rho \dot{\mathbf{u}}\|} & \text{if } \|\mathbf{r} - \rho \dot{\mathbf{u}}\| > \mu N \end{cases} = \mathbf{0} \qquad (5) $$

which can be rewritten as a nonsmooth vector equality utilizing the max and min operators

$$ \mathbf{r} - \frac{\mathbf{r} - \rho \dot{\mathbf{u}}}{\|\mathbf{r} - \rho \dot{\mathbf{u}}\| - \mu N} \left( \min(\|\mathbf{r} - \rho \dot{\mathbf{u}}\| - \mu N, 0) + \frac{\mu N}{\|\mathbf{r} - \rho \dot{\mathbf{u}}\|} \max(\|\mathbf{r} - \rho \dot{\mathbf{u}}\| - \mu N, 0) \right) = \mathbf{0}. \qquad (6) $$

Note that this formulation (Equation (6)) may face singularities when initial conditions are taken to be zero, due to the denominator in front of the max function. The projective version in Equation (5) can be shown to be equivalent to Equation (2) by simplifying the top expression and multiplying the bottom expression by $\|\mathbf{r} - \rho \dot{\mathbf{u}}\|$, which yields:

$$ \begin{cases} \dot{\mathbf{u}} = \mathbf{0} & \text{if } \|\mathbf{r} - \rho \dot{\mathbf{u}}\| \le \mu N, \\ \|\mathbf{r} - \rho \dot{\mathbf{u}}\| \|\mathbf{r} - \mu N(\mathbf{r} - \rho \dot{\mathbf{u}})\| = \mathbf{0} & \text{if } \|\mathbf{r} - \rho \dot{\mathbf{u}}\| > \mu N. \end{cases} \qquad (7) $$

As a second example, an alternate 2D friction formulation was proposed in [79], where the tangential velocity and friction force are parameterized using polar coordinates $(\dot{u}_s, \theta)$ and $(r_s, \theta)$ as described in [14, p. 431] as

$$ \begin{cases} \dot{u}_1(t) - \dot{u}_s(t) \cos \theta(t) = 0 \\ \dot{u}_2(t) - \dot{u}_s(t) \sin \theta(t) = 0 \end{cases} \quad \text{and} \quad \begin{cases} r_1(t) - r_s(t) \cos \theta(t) = 0 \\ r_2(t) - r_s(t) \sin \theta(t) = 0 \end{cases} \qquad (8) $$

along with the 1D friction equality in $(\dot{u}_s, r_s)$

$$ \dot{u}_s + \min(0, \rho(r_s + \mu N) - \dot{u}_s) + \max(0, \rho(r_s - \mu N) - \dot{u}_s) = 0, \quad \rho > 0. \qquad (9) $$

However, the proposed solution procedure then relies on a regularization of the min and the max operators. Also, while Equation (2) consists of two scalar equations with four unknowns $(\dot{u}_1, \dot{u}_2, r_1, r_2)$, the formulation in [79] involves seven unknowns $(\dot{u}_1, \dot{u}_2, r_1, r_2, \dot{u}_s, r_s, \theta)$ for five equations, collectively in Equations (8) and (9); however, it should be noted that the four equations in Equation (8) can be reduced to the three equations

$$ \begin{aligned} & \dot{u}_1 r_s - \dot{u}_s r_1 = 0 \\ & \dot{u}_2 r_s - \dot{u}_s r_2 = 0 \\ & \dot{u}_1^2 + \dot{u}_2^2 - \dot{u}_s^2 = 0 \quad (\text{or} \quad r_1^2 + r_2^2 - r_s^2 = 0) \end{aligned} \qquad (10) $$
in the six unknowns ($\dot{u}_1, \dot{u}_2, r_1, r_2, \dot{u}_s, r_s$) only. While Equation (9) is a good equality version of Coulomb's friction in 1D, that is when $\theta(t) = 0$ as shown in [81], it does not seem to be ideal in the context of 2D friction as it introduces new unknowns and possibly spurious solutions, notably because the sign of $\dot{u}_s$ and $r_s$ in Equation (8) or Equation (10) is not clearly defined. This formulation was not tested in the present work.

As a third example, the slightly simpler formulation in [77] suppresses the extra unknown $\theta(t)$ and addresses the problem of determining the direction of motion by introducing an equivalent nonsmooth equation, such that the 2D friction equalities read

$$ \min(\|\dot{\mathbf{u}}\|, \mu N - \|\mathbf{r}\|) = 0 \qquad (11a) $$

$$ |\dot{u}_1 r_2 - \dot{u}_2 r_1| + \max(\dot{u}_1 r_1, 0) = 0 \qquad (11b) $$

where Equation (11a) governs the magnitude of the friction force while Equation (11b) governs the directions of the friction force and tangential velocity. This formulation was found to suffer from convergence issues and thus not deemed suitable.

In conclusion, one notes that in the vector Equations (2), (5) and (6), each of the two scalar equations corresponds to one component of the friction force and tangential velocity vectors. Hence these formulations feature a desirable symmetry in the two directions of contact, so that no condition on the relative direction of the friction force and the velocity is required, such as in Equation (9) and Equation (11). In the implicit formulations, the direction of the friction force is inherently embedded in the friction conditions. Furthermore, the formulation in Equation (2) alleviates the singularity problem present in Equation (6). Because of these features, as well as the robust convergence of the nonlinear solution algorithm associated with the weighted residual procedure, Equation (2) was selected in the present study.

## 2.3 Equality-Based Approximations

To avoid challenges associated with the discontinuous slopes in Coulomb's law, it is common to use smoothing functions to regularize the nonsmooth behavior. Single-valued functions are also used to circumvent the difficulties caused by the set-valued nature of Coulomb's law. These approximate methods do not have convergence issues, since there is no jump in the friction force. However, with these approximations, no true sticking can take place and the system will always slip. Also, as the single-valued function approaches Coulomb's law, the equations can become numerically stiff.

The equality-based formulations developed in Sections 2.1 and 2.2 are exact representations of Coulomb's law and thus alleviate the need for any kind of regularization. However, regularized versions of Coulomb's law may still be cast in the form of equivalent equalities, as is described below for two popular regularized functions, namely hyperbolic tangent and piecewise linear functions. Note that in addition to these, the signum function has also been used for quasi-static friction systems [47]. However, this requires considering a 2D system as two uncoupled 1D systems, which often is an unrealistic assumption.

Regularized equalities of 2D friction can be formulated first by writing that the friction force vector is collinear with and in opposite direction to the velocity vector, that is

$$ \mathbf{r} = -\|\mathbf{r}\|\frac{\dot{\mathbf{u}}}{\|\dot{\mathbf{u}}\|}. \qquad (12) $$

In Equation (12) and in the remainder, the convention $0/\|0\| = 0$ is used.

### 2.3.1 Single-Valued Hyperbolic Tangent Regularization

The norm of the friction force $\|\mathbf{r}\|$ is regularized as [67]:

$$ \|\mathbf{r}\| = \mu N \tanh(\alpha \|\dot{\mathbf{u}}\|) \qquad (13) $$

where $\alpha > 0$ is a regularization parameter, such that when $\alpha$ increases, the regularized function better approximates Coulomb's law (an insightful demonstration is provided in [81]). Substituting Equation (13) into Equation (12), one obtains

$$ \mathbf{r} + \mu N \tanh(\alpha \|\dot{\mathbf{u}}\|)(\dot{\mathbf{u}}/\|\dot{\mathbf{u}}\|) = 0 \qquad (14) $$

which is a smooth vector equality representation of the regularized approximation of Coulomb's 2D friction law. It converts the set-valued mapping in Equation (2) into a single-valued mapping such that the corresponding dynamics
is no longer governed by nonsmooth Differential Algebraic Equations but rather by smooth Ordinary Differential Equations. Accordingly, it is not as accurate as the exact equality-based representation of Coulomb's law and results in always-sliding responses where no true sticking motion can be captured. For large $\alpha$, sticking can be very well approximated at the cost of possibly generating stiff Ordinary Differential Equations, as illustrated by the steep shadowed area in Figure 3(a) in [81].

### 2.3.2 Single-Valued Piecewise Linear Approximation

Coulomb's law can also be approximated using a piecewise linear function. In this case the norm of the friction force is expressed as

$$ \|\mathbf{r}\| = \begin{cases} \frac{\mu N \|\dot{\mathbf{u}}\|}{\delta} & \text{if } \|\dot{\mathbf{u}}\| \le \delta \\ \mu N & \text{if } \|\dot{\mathbf{u}}\| > \delta \end{cases} = \mu N \min(\|\dot{\mathbf{u}}\|, \delta)/\delta, \quad \delta > 0 \quad (15) $$

where $\delta$ controls the slope of the single-valued function in the approximated 'sticking' region, such that Coulomb's law is better approximated for small values of $\delta$ [81]. Substituting Equation (15) into Equation (12) results in the single-valued approximated vector equality formulation of Coulomb's law

$$ \mathbf{r} + \mu N(\dot{\mathbf{u}}/\|\dot{\mathbf{u}}\|) \min(\|\dot{\mathbf{u}}\|, \delta)/\delta = \mathbf{0}. \quad (16) $$

Interestingly, Equation (16) is a nonsmooth equality due to the slope discontinuities of the piecewise linear operator $\min$, thus constituting an approximation rather than a regularization. However, the piecewise linear approximation is single-valued, as $\mathbf{r}$ in Equation (15) is a single-valued function of $\dot{\mathbf{u}}$. Similarly to the tanh regularization above, the resulting Ordinary Differential Equations become stiff when $\delta$ tends to 0, see Figure 4 in [81].

### 2.3.3 Set-Valued Regularization

The regularization using the tanh function is very classical (at least in a one-dimensional setting) to suppress the nonsmoothness of the mapping. However, by transforming the set-valued nature of the mapping into a single-valued function $\mathbf{r}(\dot{\mathbf{u}})$, it is unable to capture true sticking phases in the resulting dynamics. Instead, it is possible to regularize the max operator arising in Equation (2) while (partially) keeping the set-valued feature of the mapping. Among others, one option is to proceed as follows [83, 84, 85]:

$$ \max(a, b) = \frac{b + a + |b - a|}{2} \quad \text{and} \quad |b - a| \approx \sqrt{(b - a)^2 + \epsilon} \quad \rightarrow \quad \max_{\epsilon}(a, b) = \frac{b + a + \sqrt{(b - a)^2 + \epsilon}}{2} \quad (17) $$

where $a$ and $b$ are two arbitrary real numbers. In Equation (17), $\epsilon$ dictates the level of regularization in the approximation, as illustrated in Figure 2 for 1D Coulomb friction; a smaller value corresponds to a "better" regularization. This type of regularization does not introduce stiff terms at the cost of not converting an Algebraic Differential Equation into an Ordinary Differential Equation. Also, it should be noted that the selected regularization using the square-root function has a global effect on the max operator, which decreases away from $a \approx b$. However, it would also be possible to regularize locally using a true arc of radius $\epsilon$.

# 3 Weighted Residual Formulation for Periodic Solutions

## 3.1 Two-Degree-of-Freedom System with Two-Dimensional Friction

Consider the one-mass system moving in the plane $O\mathbf{i}_1\mathbf{i}_2$ with displacements $u_1(t)$ and $u_2(t)$, stored in vector $\mathbf{u} = [u_1, u_2]^T$, and connected to the ground via springs $k_1$ and $k_2$ and viscous dampers $c_1$ and $c_2$ along the $\mathbf{i}_1$ and $\mathbf{i}_2$ directions, respectively, as shown in Figure 3. For simplicity, assume that spring and damper restoring forces remain linear and uncoupled in each direction. The mass is subject to external forces $f_1(t)$ and $f_2(t)$, periodic with frequency $\omega = 2\pi/T$. It is also subject to a force $N$ normal to the plane of motion, resulting in 2D frictional contact, and the resulting friction forces are denoted as $r_1(t)$ and $r_2(t)$, stored in vector $\mathbf{r} = [r_1, r_2]^T$. The case of unequal stiffnesses in the two directions of motion, $k_1 \neq k_2$, is referred to as anisotropic stiffness, while $k_1 = k_2$ is the isotropic stiffness case. Note that viscous damping is desirable to avoid the occurrence of singular impedance matrices for the harmonics of the linear terms in the HB procedure.
<figure><img src="image_5.png" /><figcaption>Fig. 2: 1D equality-based version of Coulomb's friction in direction i₁ only (that is r₂ = ṁ₂ = 0 in Equation (2)), for illustration purposes: [blue] Ψ(r₁, ṁ₁) = r max(μN, |r₁ - ρû₁|) - μN(r₁ - ρû₁), [orange] 0 plane, [solid black line] dry friction solution set seen as the zero-level set of the function Ψ, that is the set (r₁, ṁ₁) satisfying Ψ(r₁, ṁ₁) = 0; [left] nonsmooth function Ψ involving the max operator; [right] set-valued regularization of Ψ where max ≈ maxε through Equation (17) with ε = 0.1: stiffening is avoided.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3: Top view of two-DOF mass-spring-damper system with 2D frictional interface.</figcaption></figure>

The equations of motion of the mass and the two equality-based friction conditions derived from the implicit formulation in Equation (2) can be written as:

$$ \psi_{u_1}(u_1, r_1) = m\ddot{u}_1 + c_1\dot{u}_1 + k_1u_1 - r_1 - f_1 = 0 \qquad (18a) $$

$$ \psi_{u_2}(u_2, r_2) = m\ddot{u}_2 + c_2\dot{u}_2 + k_2u_2 - r_2 - f_2 = 0 \qquad (18b) $$

$$ \psi_{r_1}(\dot{u}_1, \dot{u}_2, r_1, r_2) = r_1 \max(\mu N, \|\mathbf{r} - \rho \dot{\mathbf{u}}\|) - \mu N(r_1 - \rho \dot{u}_1) = 0 \qquad (18c) $$

$$ \psi_{r_2}(\dot{u}_1, \dot{u}_2, r_1, r_2) = r_2 \max(\mu N, \|\mathbf{r} - \rho \dot{\mathbf{u}}\|) - \mu N(r_2 - \rho \dot{u}_2) = 0 \qquad (18d) $$

where ρ is a user-defined scaling parameter of default value 1. In Equation (18), the unknowns are u₁(t), u₂(t), r₁(t) and r₂(t). Note that the two nonsmooth friction conditions take symmetric expressions, one along i₁ and the other along i₂, and are coupled through the terms in ||r - ρû||. It is also interesting to note that the two (linear) equations of motion are augmented by two (nonsmooth) friction conditions, resulting in a *primal-dual mixed formulation* of the problem in terms of both displacements (and their derivatives) and friction forces. This contrasts with frequency-time domain methods such as the AFT and DLFT, which are *primal* formulations in the displacements only (or velocities), such that friction and contact forces are calculated via time integration at each iteration of the nonlinear algorithm from the displacements. For the mixed displacement-force formulation in Equation (18), all unknowns including the friction forces are solved for simultaneously in the frequency domain via the nonlinear algorithm.
## 3.2 Ritz-Galerkin Solution Procedure

In order to solve for the steady-state periodic solutions of the system governed by Equation (18), all unknowns are discretized by expanding them on a truncated basis of T-periodic functions with $N_\phi$ members, as

$$u_1(t) \approx u_{1h}(t) = \sum_{k=1}^{N_\phi} u_{1k} \phi_k(t), \quad u_2(t) \approx u_{2h}(t) = \sum_{k=1}^{N_\phi} u_{2k} \phi_k(t) \qquad (19a)$$

$$r_1(t) \approx r_{1h}(t) = \sum_{k=1}^{N_\phi} r_{1k} \phi_k(t), \quad r_2(t) \approx r_{2h}(t) = \sum_{k=1}^{N_\phi} r_{2k} \phi_k(t) \qquad (19b)$$

where $u_{1h}(t)$, $u_{2h}(t)$, $r_{1h}(t)$ and $r_{2h}(t)$ are the series approximations, $\phi_k(t)$ is the kth basis function, and $u_{1k}$, $u_{2k}$, $r_{1k}$ and $r_{2k}$ are the $4N_\phi$ unknowns of the discretized system of equations. Time derivatives of the displacement unknowns should also be expanded as:

$$\dot{u}_1(t) \approx \dot{u}_{1h}(t) = \sum_{k=1}^{N_\phi} u_{1k} \dot{\phi}_k(t), \quad \dot{u}_2(t) \approx \dot{u}_{2h}(t) = \sum_{k=1}^{N_\phi} u_{2k} \dot{\phi}_k(t) \qquad (20a)$$

$$\ddot{u}_1(t) \approx \ddot{u}_{1h}(t) = \sum_{k=1}^{N_\phi} u_{1k} \ddot{\phi}_k(t), \quad \ddot{u}_2(t) \approx \ddot{u}_{2h}(t) = \sum_{k=1}^{N_\phi} u_{2k} \ddot{\phi}_k(t) \qquad (20b)$$

Since the nonsmooth friction equalities do not contain time derivatives of friction forces $r_1$ and $r_2$, the procedure does not require differentiation of these nonsmooth functions. Substituting the approximate series expansions in Equations (19) and (20) into Equation (18), the following non-zero residuals are obtained:

$$\psi_{u_1}(u_{1h}, r_{1h}) = m\ddot{u}_{1h} + c_1\dot{u}_{1h} + k_1u_{1h} - r_{1h} - f_1 \qquad (21a)$$

$$\psi_{u_2}(u_{2h}, r_{2h}) = m\ddot{u}_{2h} + c_2\dot{u}_{2h} + k_2u_{2h} - r_{2h} - f_2 \qquad (21b)$$

$$\psi_{r_1}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) = r_{1h} \max(\mu N, \|\mathbf{r}_h - \rho \dot{\mathbf{u}}_h\|) - \mu N(r_{1h} - \rho \dot{u}_{1h}) \qquad (21c)$$

$$\psi_{r_2}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) = r_{2h} \max(\mu N, \|\mathbf{r}_h - \rho \dot{\mathbf{u}}_h\|) - \mu N(r_{2h} - \rho \dot{u}_{2h}) \qquad (21d)$$

Above, $\psi_{u_1}(u_{1h}, r_{1h})$ and $\psi_{u_2}(u_{2h}, r_{2h})$ are the residuals of the equations of motion, while $\psi_{r_1}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h})$ and $\psi_{r_2}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h})$ are the residuals of the two friction conditions, for the approximation of order $N_\phi$. The residuals are enforced to be orthogonal to the basis functions $\phi_k$ for $k = 1, \dots, N_\phi$ through a Ritz-Galerkin projection:

$$\int_{0}^{T} \phi_k(t) \psi_{u_1}(u_{1h}, r_{1h}) dt = 0, \quad k = 1, \dots, N_\phi \qquad (22a)$$

$$\int_{0}^{T} \phi_k(t) \psi_{u_2}(u_{2h}, r_{2h}) dt = 0, \quad k = 1, \dots, N_\phi \qquad (22b)$$

$$\int_{0}^{T} \phi_k(t) \psi_{r_1}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) dt = 0, \quad k = 1, \dots, N_\phi \qquad (22c)$$

$$\int_{0}^{T} \phi_k(t) \psi_{r_2}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) dt = 0, \quad k = 1, \dots, N_\phi. \qquad (22d)$$

Equation (22) consists of a nonlinear system of $4N_\phi$ equation with $4N_\phi$ unknowns. Equations (22a) and (22b) govern the dynamics of the system and are linear in the unknowns, while Equations (22c) and (22d) are nonlinear due to the presence of the nonsmooth max function. The accuracy of the approximate solution is dictated by the selection of the functions $\phi_k$ and their number $N_\phi$, such that the norm of the residuals decreases as the order $N_\phi$ increases.

## 3.3 Fourier Series Expansion in Time

Here the shape functions $\phi_k(t)$ are chosen to be the classical Fourier functions, although this is not a limitation and other functions (either smooth or nonsmooth, such as periodized wavelets [86]) could be used. The unknowns are
expanded into Fourier series with $N_\phi$ coefficients as

$$u_1(t) \approx u_{1h}(t) = \sum_{k=1}^{N_\phi/2} u_{1_{2k-1}} \cos(2k-1)\omega t + u_{1_{2k}} \sin(2k-1)\omega t \quad (23a)$$

$$u_2(t) \approx u_{2h}(t) = \sum_{k=1}^{N_\phi/2} u_{2_{2k-1}} \cos(2k-1)\omega t + u_{2_{2k}} \sin(2k-1)\omega t \quad (23b)$$

$$r_1(t) \approx r_{1h}(t) = \sum_{k=1}^{N_\phi/2} r_{1_{2k-1}} \cos(2k-1)\omega t + r_{1_{2k}} \sin(2k-1)\omega t \quad (23c)$$

$$r_2(t) \approx r_{2h}(t) = \sum_{k=1}^{N_\phi/2} r_{2_{2k-1}} \cos(2k-1)\omega t + r_{2_{2k}} \sin(2k-1)\omega t \quad (23d)$$

where each odd harmonic of the response contains a cosine and a sine term. The periodic excitation is assumed to contain only odd harmonics of time, such that there are no contributions of the even harmonic or constant terms in the response, due to the odd nature symmetry of the Coulomb friction nonlinearity, which here simplifies the procedure. It is expanded 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 (24a)$$

$$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 (24b)$$

For compactness, Fourier coefficients are stored in the vectors $\mathbf{u}_1 = [u_{11}, \dots, u_{1N_\phi}]^T$, $\mathbf{u}_2 = [u_{21}, \dots, u_{2N_\phi}]^T$, $\mathbf{r}_1 = [r_{11}, \dots, r_{1N_\phi}]^T$, $\mathbf{r}_2 = [r_{21}, \dots, r_{2N_\phi}]^T$, $\mathbf{f}_1 = [f_{11}, \dots, f_{1N_\phi}]^T$, and $\mathbf{f}_2 = [f_{21}, \dots, f_{2N_\phi}]^T$. In the HB method, the expansions in Equations (23) and (24) are substituted into Equation (21), pre-multiplied by the $N_\phi$ Fourier functions, and then integrated over one period of the motion, from 0 to $T = 2\pi/\omega$. The yields the system of equations<fnref n="1" />:

$$\mathbf{A}_{11}\mathbf{u}_1 - \mathbf{r}_1 - \mathbf{f}_1 = \mathbf{0} \quad (25a)$$

$$\mathbf{A}_{22}\mathbf{u}_2 - \mathbf{r}_2 - \mathbf{f}_2 = \mathbf{0} \quad (25b)$$

$$\mathbf{g}_{\mathbf{r}_1}(\mathbf{u}_1, \mathbf{u}_2, \mathbf{r}_1, \mathbf{r}_2) = \mathbf{0} \quad (25c)$$

$$\mathbf{g}_{\mathbf{r}_2}(\mathbf{u}_1, \mathbf{u}_2, \mathbf{r}_1, \mathbf{r}_2) = \mathbf{0} \quad (25d)$$

where $\mathbf{A}_{ii} = k_i \mathbf{I} + c_i \mathbf{D} - m\Omega^2$, $i = 1, 2$, $\mathbf{I}$ is the identity matrix, and the block matrices $\Omega$ and $\mathbf{D}$ are

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

Additionally, the components of the vectors $\mathbf{g}_{\mathbf{r}_1}$ and $\mathbf{g}_{\mathbf{r}_2}$ in Equations (25c) and (25d) take the form

$$g_{r_1 2k-1}(\mathbf{u}_1, \mathbf{u}_2, \mathbf{r}_1, \mathbf{r}_2) = \int_0^T \psi_{r_1}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) \cos(2k-1)\omega t dt = 0, \quad k=1, \dots, N_\phi/2 \quad (27a)$$

$$g_{r_1 2k}(\mathbf{u}_1, \mathbf{u}_2, \mathbf{r}_1, \mathbf{r}_2) = \int_0^T \psi_{r_1}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) \sin(2k-1)\omega t dt = 0, \quad k=1, \dots, N_\phi/2 \quad (27b)$$

$$g_{r_2 2k-1}(\mathbf{u}_1, \mathbf{u}_2, \mathbf{r}_1, \mathbf{r}_2) = \int_0^T \psi_{r_2}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) \cos(2k-1)\omega t dt = 0, \quad k=1, \dots, N_\phi/2 \quad (27c)$$

$$g_{r_2 2k}(\mathbf{u}_1, \mathbf{u}_2, \mathbf{r}_1, \mathbf{r}_2) = \int_0^T \psi_{r_2}(\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) \sin(2k-1)\omega t dt = 0, \quad k=1, \dots, N_\phi/2. \quad (27d)$$

Classical dynamic condensation is then performed by solving Equations (25a) and (25b) for $\mathbf{u}_1$ and $\mathbf{u}_2$, as

$$\mathbf{u}_1(\mathbf{r}_1) = \mathbf{A}_{11}^{-1}(\mathbf{r}_1 + \mathbf{f}_1) \quad (28a)$$

$$\mathbf{u}_2(\mathbf{r}_2) = \mathbf{A}_{22}^{-1}(\mathbf{r}_2 + \mathbf{f}_2) \quad (28b)$$

<footnote n="1">¹A strategy similar to the proposed solution method is very briefly mentioned in [87, Eq. (3.39)] but preference is then given to the DLFT by the author.</footnote>
and then substituting Equation (28) into Equations (25c) and (25d) to yield

$$ \mathbf{g}_{r_1}(\mathbf{u}_1(\mathbf{r}_1), \mathbf{u}_2(\mathbf{r}_2), \mathbf{r}_1, \mathbf{r}_2) = \mathbf{0}, \qquad (29a) $$

$$ \mathbf{g}_{r_2}(\mathbf{u}_1(\mathbf{r}_1), \mathbf{u}_2(\mathbf{r}_2), \mathbf{r}_1, \mathbf{r}_2) = \mathbf{0}. \qquad (29b) $$

Equation (29) consists of $2N_φ$ equations in only the $2N_φ$ unknowns $\mathbf{r}_1$ and $\mathbf{r}_2$. These coupled, nonsmooth equations can be solved using a nonlinear solver such as one based on Powell's dog leg trust region procedure, which is used in this study through the MATLAB® built-in function `fsolve`. The integrals involving nonsmooth terms in Equation (27) are calculated numerically using a FFT. For the results in this study, the number of samples per period varied between $2^9 = 512$ and $2^{14} = 16384$ depending on the number of harmonics considered.

## 3.4 Piecewise Linear Jacobian

An expression for the weak formulation of the Jacobian can be derived and provided as an input to the `fsolve` built-in command. This is achieved by writing the Jacobian as a piecewise linear expression which corresponds to the various states of the system. The Jacobian matrix of the nonsmooth functions $\mathbf{g}_{r_1}$ and $\mathbf{g}_{r_2}$ in Equations (25c) and (25d) is denoted as

$$ \mathbf{J} = \begin{bmatrix} \nabla_{r_1} \mathbf{g}_{r_1} & \nabla_{r_2} \mathbf{g}_{r_1} \\ \nabla_{r_1} \mathbf{g}_{r_2} & \nabla_{r_2} \mathbf{g}_{r_2} \end{bmatrix} \qquad (30) $$

the elements of which are the block matrices

$$ (\nabla_{\mathbf{r}_n} \mathbf{g}_{\mathbf{r}_m})_{ij} = \int_0^T \phi_i(t) \frac{\partial \psi_{r_m}}{\partial r_{n_j}} (\dot{u}_{1h}, \dot{u}_{2h}, r_{1h}, r_{2h}) dt, \quad m, n = 1, 2, \quad i, j = 1, \dots, N_\phi. \qquad (31) $$

Using the chain rule leads to

$$ \frac{\partial \psi_{r_m}}{\partial r_{n_j}} = \frac{\partial \psi_{r_m}}{\partial r_n} \frac{\partial r_n}{\partial r_{n_j}} + \frac{\partial \psi_{r_m}}{\partial \dot{u}_n} \frac{\partial \dot{u}_n}{\partial r_{n_j}}. \qquad (32) $$

The partial derivatives of the nonsmooth functions $\psi_{r_m}$ with respect to $r_n$ and $\dot{u}_n$ in the right hand side of Equation (32) can be expressed using the Heaviside function $H(x)$ as follows:

$$ \frac{\partial \psi_{r_m}}{\partial r_n} = r_m \frac{r_n - \rho \dot{u}_n}{\|\mathbf{r} - \rho \dot{\mathbf{u}}\|} H(\|\mathbf{r} - \rho \dot{\mathbf{u}}\| - \mu N) + (\max(\mu N, \|\mathbf{r} - \rho \dot{\mathbf{u}}\|) - \mu N) \delta_{mn}, \qquad (33a) $$

$$ \frac{\partial \psi_{r_m}}{\partial \dot{u}_n} = -\rho r_m \frac{r_n - \rho \dot{u}_n}{\|\mathbf{r} - \rho \dot{\mathbf{u}}\|} H(\|\mathbf{r} - \rho \dot{\mathbf{u}}\| - \mu N) + \rho \mu N \delta_{mn}, \qquad (33b) $$

where $\delta_{mn}$ is the Kronecker delta. From Equation (19b), one can simply infer that $\partial r_n(t)/\partial r_{n_j} = \phi_j(t)$ for $n = 1, 2$ and $j = 1, \dots, N_\phi$. Also, using Equation (28) to express the velocity vectors $\dot{\mathbf{u}}_n$, while considering Equation (20a) and then differentiating with respect to $r_{n_j}$, results in $\partial \dot{\mathbf{u}}_n/\partial \mathbf{r}_n = (A_{mn}^{-1})^T \dot{\phi}$, with $\dot{\phi} = [\phi_1, \dots, \phi_{N_\phi}]^T$. The Jacobian (Equation (31)) can then be computed by testing the time dependent conditions in Equation (33) for every discrete time $t_i$ of the Heaviside function argument. Then, the integral in Equation (31) is calculated using a FFT.

In the remainder, the implicit formulation of friction and the attendant Ritz Galerkin procedure for periodic solutions, as described in Sections 3.2 to 3.4, are referred to as the Equality-based Weighted Residual (E-WR) method and are summarized in Figure 4. Note that some the integrals at the Ritz-Galerkin projection step, namely those involving nonsmooth terms, do not have a closed form and thus must be approximated numerically using a FFT. The quasi-analytical Jacobian in Equation (31) is called by the nonlinear solver at each iteration. The solution procedure is highly adaptive to other exact or regularized friction formulations, as Equations (21c) and (21d) can simply be replaced by any alternative formulation of interest, the other steps of the procedure remaining exactly the same. Of course the Jacobian should be derived separately for other formulations.

# 4 Results for a One-mass Oscillator with Two-Dimensional Friction

Steady-state periodic responses for the 2D frictional system with two DOFs shown in Figure 3 are calculated for $m = 1$ and $c_1 = c_2 = 0.01$. In most cases below, anisotropic stiffnesses in the two directions of motion are considered
<figure><img src="image_7.png" /><figcaption>Fig. 4: E-WR flowchart for the calculation of periodic responses.</figcaption></figure>

$(k_1 = 1$ and $k_2 = 2)$. Some cases of isotropic stiffnesses ($k_1 = k_2 = 1$) are also treated. Various values of the maximum friction force $\mu N$ are tested, as well as several cases of external excitations $f_1(t)$ and $f_2(t)$. The natural frequencies of the corresponding linear system (i.e. $\mu N = 0$) are $\omega_1^2 = k_1/m$ and $\omega_2^2 = k_2/m$, corresponding to uncoupled modes of vibration in the $\mathbf{i}_1$ and $\mathbf{i}_2$ directions, respectively. Note that the linear equations of motion Equation (21a) and Equation (21b) are uncoupled in the $\mathbf{i}_1$ and $\mathbf{i}_2$ directions so as to highlight the coupling introduced by the friction law.

## 4.1 Frequency Responses

Frequency response curves of the one-mass system are represented in Figure 5 for various values of $\mu N$ over a frequency range encompassing the two modes of the linearized system. Using the notation $\mathbf{u} = [\mathbf{u}_1, \mathbf{u}_2]^T$, the total amplitude $||\mathbf{u}||$ of the displacement is shown. The anisotropic stiffness case $k_1 = 1$ and $k_2 = 2$ is considered, such that the natural frequencies of the linearized system are $\omega_1 = 1$ and $\omega_2 = \sqrt{2}$. The excitation forces are $f_1(t) = f_2(t) = 10 \cos \omega t$. $N_\phi = 60$ harmonics are used in the Fourier expansions to ensure the accuracy of the HB solution.

<figure><img src="image_8.png" /><figcaption>Fig. 5: Frequency response of the one-mass system in the anisotropic case for various values of $\mu N$ and $N_\phi = 60$.</figcaption></figure>

As the normal load increases, as expected the effects of friction become more prevalent and response amplitudes decrease. For $\mu N = 5$ and 8, the two resonant peaks are observed at the resonant frequencies of the linearized system
with viscous damping, but with lower amplitudes. In low frequency ranges (i.e. $\omega < 0.3$), multiple super-harmonic resonances are noticeable as well as their associated anti-resonances, which are characteristic of an increasing number of stick-slip transitions per period of the motion. For large values of the normal load, $\mu N = 10$ and 12, the two primary resonant peaks disappear. Also note the apparition of weak resonance peaks in the range $\omega \in [0.55, 0.75]$ when friction is present.

No convergence problems were observed in Figure 5 for this relatively large number of harmonics, $N_\phi = 60$. The scaling coefficient $\rho$ in Equations (21c) and (21d) was taken equal to 1 in all cases. For each value of the normal force, the frequency response was generated by continuation from $\omega = 0$ to 1.5 in increments of 0.001. The required number of iterations in the nonlinear dog leg algorithm ranged from 0 (for the linear case) to 120 (for $\mu N = 12$) for a tolerance of $10^{-5}$. The initial guess for the nonlinear solver at the first frequency point was taken as the zero friction force, and for the subsequent frequencies it was taken as the solution at the previous frequency. The number of iterations was large only for the first frequency and decreased to less than 10 for the subsequent frequencies. The computational time needed to generate one such frequency response curve on a personal computer (CPU: Intel(R) Core(TM) i7-4750HQ 2.00 GHz) was approximately 72 minutes or about 3 s per frequency point.

## 4.2 Typical Periodic Solutions

In this section, periodic responses of the mass are presented at specific excitation frequencies. Unless otherwise indicated, the anisotropic stiffness case is considered. Figure 6 shows the response at the first resonant frequency of the linearized system, $\omega = 1$, for $\mu N = 8$. In this case, the mass is seen to undergo a harmonic-like resonant motion, with larger amplitudes of the displacement and velocity in the $i_1$ than the $i_2$ direction due to the smaller stiffness along $i_1$. There is no sticking and the mass slides at all times, which can be inferred from the velocity plot since there is no

<figure><img src="image_9.png" /><figcaption>**Fig. 6:** Response at the first resonance of the linearized ($\mu N = 0$) system, $\omega = 1$, for $\mu N = 8$ and $N_\phi = 160$. Displacement and velocity plots indicate the $i_1$-direction DOF [solid line] and $i_2$-direction DOF [dashed line].</figcaption></figure>

instant of time when both components of the velocity are zero. The trajectory of the mass on the plane of contact resembles an ellipse. In contrast to some existing methods which make assumptions about the trajectory or the slip angle [88], the friction point trajectory is obtained here directly from the E-WR procedure without any *a priori* assumption, which is an advantage of this approach. The components of the friction force are also depicted over one period of motion in Figure 6. The locus of the friction force vector is observed to remain on the perimeter of the Coulomb cone section (i.e., $\mu N = 8$) at all times, indicating that the oscillator undergoes pure sliding motion. This demonstrates that Coulomb's law is accurately satisfied by the E-WR solution.

Figure 7 shows the system response at $\omega = 0.17$ rad/s, corresponding to the largest super-harmonic resonance for $\mu N = 8$. It is observed that the mass is stuck for approximately the first 9.1 seconds of the period, as indicated by its zero velocity. Another sticking phase is seen from $t = 18.4$ s to 27.5 s. These two sticking phases per cycle are depicted in the plot of the trajectory of the mass as red circles. The particle starts from the top right end of the trajectory, then begins its motion on the lower branch until it arrives at the bottom left red circle for the second sticking phase. After the second rest, the oscillator follows the path on the upper branch to get back to the starting point. Also note the occurrence of two instants of near-sticking at $t = 14.4$ s and $t = 32.9$ s, when both velocity components vanish, but only instantaneously so. The plot of the locus of the friction force vector shows that it remains on the Coulomb cone during the sliding phases, as indicated by the circular sectors at the upper right and bottom left of the plot, while the lines connecting these sectors lie inside the friction cone and correspond to the sticking phases, where the friction force takes multiple values, all smaller than 8. Thus in this case again, Coulomb's law is captured with high accuracy in the solution procedure.

Decreasing the excitation frequency further, to the second largest super-harmonic resonance for $\omega = 0.1$, leads to additional sticking phases, as illustrated in Figure 8. There are four sticking phases per cycle: two long lasting ones
<figure><img src="image_10.png" /><figcaption>**Fig. 7:** Response at the largest subresonance, $\omega = 0.17$, for $\mu N = 8$ and $N_{\phi} = 160$. Red circles indicate sticking phases. Displacement and velocity plots indicate the **i**₁-direction DOF [solid line] and **i**₂-direction DOF [dashed line].</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8:** Response at the second largest subresonance, $\omega = 0.1$, for $\mu N = 8$ and $N_{\phi} = 160$. Displacement and velocity plots indicate the **i**₁-direction DOF [solid line] and **i**₂-direction DOF [dashed line].</figcaption></figure>

from $t = 0$ to 16 s and $t = 31$ to 47.7 s, and two brief ones from $t = 26.1$ to 27 s and from $t = 57.5$ to 58.5 s. These four sticking states are indicated with the red circles in the plot of the point mass trajectory, with those at the top right and bottom left corresponding to the long phases and the other two to the short phases. In the plot of the locus of the friction force vector, the two long diagonal lines inside the Coulomb cone designate the long sticking phases, while the small sharp dents away from the friction cone indicate the brief sticking phases. The circular sectors of the plot lie on the Coulomb friction cone and correspond to the sliding phases. Again, in this case of intricate response and trajectory of the mass, Coulomb's law is clearly satisfied by the E-WR method.

The next case, shown in Figure 9, is for excitation forces with a $\pi/2$ rad phase difference between the **i**₁ and **i**₂ directions, such that $f_1(t) = 10 \cos \omega t$ and $f_2(t) = 10 \sin \omega t$. The other parameters are the same as in Figure 8.

<figure><img src="image_12.png" /><figcaption>**Fig. 9:** Response to a forced excitation with components with a $\pi/2$ rad phase difference at $\omega = 0.1$, for $\mu N = 8$ and $N_{\phi} = 160$. Displacement and velocity plots indicate the **i**₁-direction DOF [solid line] and **i**₂-direction DOF [dashed line].</figcaption></figure>

Interestingly, in this scenario stick-slip behavior not observed, and the mass slides during the entire period. The mass trajectory on the contact plane resembles an ellipse, due to the phase difference between the excitation components. Also note the jagged nature of the time histories of the velocity components.

Another case, with different excitation frequencies in the two friction directions, is shown in Figure 10 for a system with isotropic stiffness, $k_1 = k_2 = 1$. Here the excitation forces are taken as $f_1(t) = 10 \cos \omega t$ and $f_2(t) = 10 \cos 3\omega t$. A complex behavior with multiple stick-slip transitions is observed. There are six sticking states of different durations, which are indicated with red circles on the plot of the trajectory of the mass. The oscillator starts the period at
<figure><img src="image_13.png" /><figcaption>**Fig. 10:** Response of a system with isotropic stiffness ($k_1 = k_2 = 1$) to a harmonic 1 excitation in the **i**₁ direction and a harmonic 3 excitation along **i**₂, for $\omega = 0.1$, $\mu N = 8$, and $N_\phi = 400$. Displacement and velocity plots indicate the **i**₁-direction DOF [solid line] and **i**₂-direction DOF [dashed line].</figcaption></figure>

the point with $u_1 = 8.8$ and $u_2 = 5.5$ and motion begins on the rightmost part of the path to reach the next stop at $u_1 = 6.7$ and $u_2 = -6.4$. It then goes through the trajectory, experiencing four more stops until it gets back to the initial point. The 3:1 frequency difference between the two components of the mass displacement is readily seen in the displacement time histories. The locus of the friction force vector shows more intricate behavior than previous cases, as indicated by the Gibbs oscillations near discontinuities just following an abrupt transition from sliding to sticking. These oscillations were found to be significant for $N_\phi = 160$, and in order to confine them to a smaller region a larger number of harmonics was used, namely $N_\phi = 400$, with a number of samples of $2^{12}$ in the FFT calculations. Interestingly, no convergence issues were observed for this very large number of harmonics, and as well, the E-WR procedure does not seem to be affected by the Gibbs phenomenon and is thus very robust. Here again, for this very intricate response behavior, Coulomb's law is accurately satisfied by the equality-based formulation.

The next case considered is for a system with isotropic stiffness in the two directions of friction, $k_1 = k_2 = 1$s, and with equal excitation components, $f_1(t) = f_2(t) = 10 \cos \omega t$. Figure 11 shows that, due to the symmetry in the system parameters and the excitations, the mass displacements and velocities are identical in the two directions. The

<figure><img src="image_14.png" /><figcaption>**Fig. 11:** Response of a system with isotropic stiffness ($k_1 = k_2 = 1$) for $\omega = 0.1$, $\mu N = 8$, and $N_\phi = 160$. Displacement and velocity plots indicate the **i**₁-direction DOF.</figcaption></figure>

mass trajectory is one-dimensional, reducing to a straight line, with four red circles denoting the positions where the oscillator sticks during its motion. Interestingly, the friction force components are also identical, meaning that the locus of the friction force vector reduces to a line: when sliding, the friction force lies solely on two points of the Coulomb cone, located at the extremities of a straight line on which the force lies during the four sticking phases. It is evident that in this case as well, Coulomb's law is fully captured by this method.

Finally, it is interesting that the 2D friction formulation developed herein can be readily be used to handle 1D friction, simply by equating to zero the excitation force in one direction. Figure 12 shows the response of a reduced one-dimensional system with $k_1 = 1$ and $f_1(t) = 10 \cos \omega t$ and $f_2(t) = 0$. As expected, the mass displacement, velocity, and friction force component in the **i**₂ direction are all equal to zero over the entire period, with motion taking place in the **i**₁ direction. A rich response with four sticking phases per cycle (two long ones and two very short ones) is shown, as indicated by four red circles for the mass trajectory. Here again the friction force locus reduces to a line, the extremities of which indicate sliding phases and lie on the Coulomb cone. This capability of the formulation to handle one-dimensional systems is a marked advantage which results from the selected nonsmooth equality, Equation (2). Other nonsmooth equalities, such as Equations (9) and (11), do not readily permit the handling of 1D friction, requiring the formulation of equalities restricted to that special case [81].
<figure><img src="image_15.png" /><figcaption>**Fig. 12:** Response of system for a 1D excitation with $k_1 = 1$ and $f_2(t) = 0$, for $\omega = 0.1$, $\mu N = 8$, and $N_\phi = 160$. Displacement and velocity plots indicate the $i_1$-direction DOF [solid line] and $i_2$-direction DOF [dashed line].</figcaption></figure>

## 4.3 Validation with Time Integration

In order to validate the E-WR method for 2D friction, the equations of motion are solved via numerical time integration using the implicit Euler method. The case shown in Figure 8 is selected due to the presence of multiple stick-slip transitions and rich behavior. A time step of $h = 10^{-3}$ is chosen for the Euler time integration, along with approximately $6.3 \times 10^4$ points per period.

Figure 13 compares the results obtained by time integration and those computed by the E-WR method. No

<figure><img src="image_16.png" /><figcaption>**Fig. 13:** Comparison of steady-state periodic response obtained by the E-WR method with $N_\phi = 160$ and by implicit Euler time integration ($h = 10^{-3}$), in the anisotropic stiffness case and for $\omega = 0.1$ and $\mu N = 8$. Displacements and velocities: E-WR [red] and Euler [black] in $i_1$ direction; E-WR [green] and Euler [blue] in $i_2$ direction. Friction force: E-WR [orange] and Euler [black].</figcaption></figure>

noticeable difference is observed between the two sets of results. Discontinuities in the friction force are captured with high accuracy by both methods, as well as the various sliding and sticking phases. Displacement and velocity plots also show an excellent match between time integration and the E-WR method.

However, the E-WR method is much more efficient in terms of computational cost. While it takes about 30 s for the Euler integration to go past the transient and reach the steady-state response, the E-WR method solves the same problem in 3.3 s, thus almost 10 times faster even for such a large number of harmonics as 160. Reasons for low computational cost are: the HB procedure directly yields the steady-state response, the E-WR formulation does not require switching between time and frequency domains as mixed methods do, and the Jacobian is calculated efficiently in a piecewise linear fashion as opposed to numerically using finite differences.

## 4.4 Convergence Study

In the E-WR method, the only approximation arises from the finite number of terms used in the Fourier series expansions, such that the accuracy of the solution increases with $N_\phi$. Figure 14 provides a typical example of how periodic responses change as $N_\phi$ increases. It is observed that for practical purposes, the time histories of the displacements, velocities and norm of the friction force are converged for $N_\phi = 50$. However, capturing Coulomb's law with high accuracy requires additional harmonics, say $N_\phi = 200$ or more. This is shown in the plots of the norm of the friction force versus that of the velocity and of the locus of the friction force, where the Gibbs oscillations corresponding to discontinuities are effectively confined to smaller regions by increasing the number of harmonics. One of the great capabilities of the E-WR method is the ability to handle very large numbers of harmonics (say more than 500) to
<figure><img src="image_17.png" /><figcaption>**Fig. 14:** Convergence analysis in terms of number of harmonics at $\omega = 0.1$ ($T = 62.83$) for $\mu N = 8$, $f_1 = 10 \cos \omega t$, $f_2 = 10 \cos 3\omega t$, and anisotropic stiffnesses $k_1 = 1$ and $k_2 = 2$. Displacement and velocity plots indicate the **i**₁-direction DOF [solid line] and **i**₂-direction DOF [dashed line].</figcaption></figure>

achieve high-accuracy solution in cases of complex 2D frictional motion and associated discontinuities, without facing any convergence issues or other numerical hardship.

A parametric study was conducted to evaluate the influence of the scaling parameter $\rho$ on the robustness of the E-WR method. Results show that the solution procedure remains stable and consistently yields accurate solutions over a wide range of $\rho$. The case in Figure 14 is considered, for $N_{\phi} = 200$, a solver tolerance of $10^{-5}$ on the norm of the residual (ie, the function value), and 4096 points per period. Table 1 presents a comparison of key performance metrics for various values of $\rho$: the final residual (i.e., the value of $\sqrt{\|g_{r_1}\|^2 + \|g_{r_2}\|^2}$ where $g_{r_1}$ and $g_{r_2}$ are calculated in Equation (29)), computational time, and the total number of iterations required for convergence. Note that each periodic solution was obtained by incrementally increasing the number of harmonics $N_{\phi}$ in the sequence
{20, 40, 60, 80, 160, 200}, as it leads to faster convergence than solving directly for the full harmonic content $N_φ = 200$. The number of iterations reported in Table 1 is the sum of the numbers of iterations for all values of $N_φ$ listed above.

The E-WR method is found to produce converged solutions over a wide range of $\rho$ spanning seven orders of magnitude, for $0.001 \le \rho \le 1000$. However, the number of iterations and computational time required by the nonlinear solve vary significantly depending on the magnitude of $\rho$. For the case considered, the best range was found to be $\rho \in [1, 100]$. This indicates that the approach is not highly sensitive to the choice of $\rho$, suggesting that it can be reliable for various applications. Results confirm that while variations in $\rho$ can affect computational cost, overall performance remains robust throughout the tested range.

<table><caption>Table 1: Effect of the scaling parameter $\rho$ on solution accuracy, computational time, and number of iterations.</caption><thead><tr><td>ρ</td><td>Residual</td><td>Iterations</td><td>Elapsed Time [s]</td><td>Status</td></tr></thead><tbody><tr><td>0.001</td><td>4.94×10<sup>-9</sup></td><td>790</td><td>310.8</td><td>Converged</td></tr><tr><td>0.01</td><td>1.25×10<sup>-11</sup></td><td>323</td><td>111.6</td><td>Converged</td></tr><tr><td>0.1</td><td>7.91×10<sup>-8</sup></td><td>139</td><td>54.9</td><td>Converged</td></tr><tr><td>1</td><td>3.22×10<sup>-13</sup></td><td>40</td><td>17.8</td><td>Converged</td></tr><tr><td>10</td><td>2.19×10<sup>-12</sup></td><td>32</td><td>15.3</td><td>Converged</td></tr><tr><td>100</td><td>3.33×10<sup>-11</sup></td><td>38</td><td>17.6</td><td>Converged</td></tr><tr><td>1000</td><td>8.99×10<sup>-11</sup></td><td>58</td><td>38.3</td><td>Inaccurate solution</td></tr></tbody></table>

## 4.5 Comparison with Piecewise Linear Approximation

The E-WR method can be readily applied to a smoothened Coulomb's law, such as the piecewise linear approximation described in Section 2.3. While the system equations of motion remain as Equations (18a) and (18b), the two nonsmooth friction conditions in Equations (18c) and (18d) are replaced with the two components of the approximated vector equality in Equation (16). Using a piecewise linear approximation means that the friction conditions are single valued rather than set valued, which generally makes them easier to solve when using existing methods such as the AFT (unless the equations become stiff due to a large slope of the piecewise linear function in the 'sticking' region). The E-WR formulation, however, remains completely transparent to whether exact or approximate Coulomb equalities are used. For all methods of course, any approximation or regularization of Coulomb's law results in discrepancies with the exact solution.

Figure 15 shows the system response (displacements, velocities, and friction force) for various values of the inverse slope parameter $\delta$ of the piecewise linear approximate friction law. The last column shows the E-WR result for the exact Coulomb's law, i.e. as $\delta \to 0$. Observe that as $\delta$ decreases, the approximate solution becomes closer to the exact solution: the quality of the piecewise linear approximation goes from poor for $\delta = 2$ and 1, to fair for $\delta = 0.5$, and to good for $\delta = 0.3$. However, note that jumps in the time history of the friction force are not captured using the piecewise linear approximation. As well, the variations of the norm of the friction force versus that of the sliding velocity show that true sticking is never experienced and thus that Coulomb's law is not satisfied. Furthermore, for $\delta < 0.3$, the equations become stiff due to large fluctuations in the gradient, resulting in convergence issues from which the exact E-WR method does not suffer. Regarding computational cost, for $N_φ = 160$ the exact formulation takes 26.3 s to converge, while the piecewise linear approximation requires 8.7 s when $\delta = 2$. Computational cost increases to 15.3 s for $\delta = 0.3$, which is the smallest value without convergence issues. Hence, the approximate method is slightly faster, but its downside is that it is an approximation of Coulomb's law. Overall, the piecewise linear approximation can be useful to obtain initial guesses for the nonlinear solver, but the exact formulation should be used for high-accuracy results.

# 5 Two-mass Oscillator with Two-Dimensional Friction

## 5.1 System of Interest and Governing Equations

To demonstrate the applicability of the E-WR formulation to larger systems, the planar motion of the two-mass model with four DOFs shown in Figure 16 is considered, where both masses are subject to periodic excitation. Mass $m_1$ is free to move frictionless on the surface, while $m_2$ is subject to a normal force $N$ and thus to friction in both $\mathbf{i}_1$ and $\mathbf{i}_2$ directions. The first mass can thus be regarded as a linear component and the second mass a nonlinear component,
<figure><img src="image_18.png" /><figcaption>Fig. 15: Periodic responses for various values of the inverse slope parameter δ (column-wise) in the piecewise linear law in Equation (16), compared to the exact E-WR formulation (far right column), for an anisotropic stiffness case with ω = 0.6 (or T = 10.47) and N<sub>φ</sub> = 160, subject to the forcing f<sub>1</sub> = f<sub>2</sub> = 10 cos ωt. Displacement and velocity plots indicate the i<sub>1</sub>-direction DOF [solid line] and i<sub>2</sub>-direction DOF [dashed line].</figcaption></figure>

each with two DOFs, a model which can be conceptually generalized to linear and nonlinear components of arbitrary
sizes, for example, a bladed disk sector (the first mass) with an attached dry friction damper (the second mass). For
simplicity, the stiffness matrices **K**<sub>1</sub> and **K**<sub>2</sub> are assumed to be diagonal, as well as the damping matrices **C**<sub>1</sub> and **C**<sub>2</sub>,
although this can be readily generalized. Thus, to linear order, motions in the two directions are uncoupled. Coupling
between the two directions of motion arises due to nonlinear friction terms only, resulting in 2D motions of the two
masses, which this study aims to investigate. The stiffness and damping matrices are

$$
\mathbf{K}_1 = \begin{bmatrix} k_{11} & 0 \\ 0 & k_{12} \end{bmatrix}, \quad \mathbf{K}_2 = \begin{bmatrix} k_{21} & 0 \\ 0 & k_{22} \end{bmatrix}, \quad \mathbf{C}_1 = \begin{bmatrix} c_{11} & 0 \\ 0 & c_{12} \end{bmatrix}, \quad \mathbf{C}_2 = \begin{bmatrix} c_{21} & 0 \\ 0 & c_{22} \end{bmatrix}. \qquad (34)
$$

For this four-DOF system, there are four equations of motion, alongside two nonsmooth equalities governing 2D
friction for mass 2, all of which can be written similarly to Equation (18) as the level-set functions

$$
\psi_{u11} = m_1 \ddot{u}_{11} + (c_{11} + c_{21}) \dot{u}_{11} - c_{21} \dot{u}_{21} + (k_{11} + k_{21}) u_{11} - k_{21} u_{21} - f_{11} = 0 \quad (35a)
$$

$$
\psi_{u12} = m_1 \ddot{u}_{12} + (c_{12} + c_{22}) \dot{u}_{12} - c_{22} \dot{u}_{22} + (k_{12} + k_{22}) u_{12} - k_{22} u_{22} - f_{12} = 0 \quad (35b)
$$

$$
\psi_{u21} = m_2 \ddot{u}_{21} + c_{21}(\dot{u}_{21} - \dot{u}_{11}) + k_{21}(u_{21} - u_{11}) - r_{21} - f_{21} = 0 \quad (35c)
$$

$$
\psi_{u22} = m_2 \ddot{u}_{22} + c_{22}(\dot{u}_{22} - \dot{u}_{12}) + k_{22}(u_{22} - u_{12}) - r_{22} - f_{22} = 0 \quad (35d)
$$
<figure><img src="image_19.png" /><figcaption>**Fig. 16:** Two-mass oscillator with mass 2 (also called the damper mass) subject to two-dimensional dry friction. Stiffness and viscous damping elements between the masses and between mass $m_1$ and the ground are shown conceptually and are assumed to reflect mathematically as in Equation (34). Directional coupling in the dynamics is only induced by friction. In the remainder, solid and dashed lines are used to indicate displacements and velocities in the two directions.</figcaption></figure>

$$ \psi_{r21} = r_{21} \max(\mu N, \|\mathbf{r} - \rho \dot{\mathbf{u}}_2\|) - \mu N(r_{21} - \rho \dot{u}_{21}) = 0 \qquad (35e) $$

$$ \psi_{r22} = r_{22} \max(\mu N, \|\mathbf{r} - \rho \dot{\mathbf{u}}_2\|) - \mu N(r_{22} - \rho \dot{u}_{22}) = 0 \qquad (35f) $$

where $\mathbf{r} = [r_{21}, r_{22}]^T$ and $\dot{\mathbf{u}}_2 = [\dot{u}_{21}, \dot{u}_{22}]^T$ are the friction force and interface (tangential) velocity vectors for mass 2, respectively. For periodic excitation of frequency $\omega$, periodic responses are sought by expanding the four displacement components and two friction force components into Fourier series in time. The Ritz-Galerkin projection procedure and the ensuing condensation to the nonlinear DOFs can be performed similarly to Section 3, resulting in the following two coupled nonlinear equations

$$ \mathbf{g}_{r_{21}}(\mathbf{u}_{21}(\mathbf{r}_{21}), \mathbf{u}_{22}(\mathbf{r}_{22}), \mathbf{r}_{21}, \mathbf{r}_{22}) = \mathbf{0}, \qquad (36a) $$

$$ \mathbf{g}_{r_{22}}(\mathbf{u}_{21}(\mathbf{r}_{21}), \mathbf{u}_{22}(\mathbf{r}_{22}), \mathbf{r}_{21}, \mathbf{r}_{22}) = \mathbf{0}. \qquad (36b) $$

where the Fourier coefficients for $u_{21}$, $u_{22}$, $r_{21}$, and $r_{22}$ are stored in the vectors $\mathbf{u}_{21}$, $\mathbf{u}_{22}$, $\mathbf{r}_{21}$, and $\mathbf{r}_{22}$, respectively, as in Section 3.3.

In the remainder, the system parameters are taken as $m_1 = m_2 = 1$ and $c_{11} = c_{12} = c_{21} = c_{22} = 0.02$. A case of anisotropic stiffnesses in the two directions of motion is considered, with $k_{11} = k_{21} = 1$ and $k_{12} = k_{22} = 2$.

## 5.2 Frequency Response and Typical Periodic Solutions

Figure 17 shows the frequency response curves for the second (nonlinear) mass for various $\mu N$ in the case of anisotropic stiffnesses. Only the first mass is subject to external excitation, such that $f_{11}(t) = f_{12}(t) = 20 \cos \omega t$ and $f_{21} = f_{22} = 0$. The total response amplitude of mass 2, $\|\mathbf{u}_2\|$, is plotted. The response of the linear system with $\mu N = 0$ is also shown. It consists of two uncoupled viscously damped 2-DOF systems in the two directions of motion, with resonant frequencies $\omega = 0.618$ and $\omega = 1.618$ in the $\dot{\mathbf{i}}_1$ direction and $\omega = 0.874$ and $\omega = 2.288$ in the $\dot{\mathbf{i}}_2$ direction. In Figure 17 only the response in the frequency range of the first two modes is shown. A relatively large number of harmonics, $N_\phi = 60$, is used in order to capture responses accurately. The scaling parameter is $\rho = 1$ in all cases, and no convergence issues were encountered. As the normal load $N$ increases, the expected decrease in the response amplitudes is observed, particularly at the primary resonances, until the resonant peaks vanish for $\mu N = 12$. Also
<figure><img src="image_20.png" /><figcaption>**Fig. 17:** Frequency response curves for the damper mass $m_2$, in a case of anisotropic stiffness ($k_{11} = k_{21} = 1$ and $k_{12} = k_{22} = 2$) for various values of the maximum friction force, for $N_{\phi} = 60$.</figcaption></figure>

note the apparition of a large number of smaller super-harmonic resonances and associated anti-resonances in the lower frequency range, indicating the occurrence of complex stick-slip responses.

For the same system parameters, Figure 18 depicts the response for $\mu N = 8$ at the primary resonance of the linearized system in the $\mathbf{i}_1$ direction, $\omega = 0.618$. Here the damper mass undergoes sliding during the entire period,

<figure><img src="image_21.png" /><figcaption>**Fig. 18:** Response at the primary resonance frequency $\omega = 0.618$, for $\mu N = 8$ and $N_{\phi} = 160$. Displacement and velocity plots indicate mass $m_1$ [blue] and mass $m_2$ [black]. Time histories of displacement and velocity indicate the $\mathbf{i}_1$-direction DOF [solid line] and $\mathbf{i}_2$-direction DOF [dashed line].</figcaption></figure>

and the trajectories of both masses are elliptical. The displacement amplitude of the first mass is smaller than that of the second mass. Also, displacements are much larger in the $\mathbf{i}_1$ direction, consistent with the resonance of the linearized system in that direction of motion. In the figure, the displacements of the two masses in the $\mathbf{i}_2$ direction are very close. The locus of the friction force vector consists of a circle on the perimeter of the Coulomb cone of friction, confirming the occurrence of sliding motion.

For the same system parameters, Figure 19 shows the response at a subresonance, for $\omega = 0.303$. The response of mass 2 features two sticking phases, shown in the figure with red circles. Sticking occurs in the time intervals [0.8, 3.1] and [11.1, 13.5], during which the velocity vanishes. While the friction mass moves back and forth between the two sticking positions during one period, the trajectory of the first mass is a more complicated and jagged, due to the absence of friction constraint. The locus of the friction force vector clearly shows sliding phases which lie on the perimeter of the Coulomb cone, and sticking phases inside it. Observe the small Gibbs oscillations for the friction force inside the Coulomb cone during the sticking phases, which could be confined to a smaller region by taking a larger number of harmonics.

At the lower subresonant frequency $\omega = 0.115$, additional phase transitions are observed. Figure 20 shows that the friction mass motion features four sticking phases per cycle: two long standing phases which are shown at the two ends of the ellipse-like trajectory, each with a duration of nearly 12 s, and two short stick phases of duration 1 s each. During the long sticking phases, the velocity of the first mass oscillates at a higher frequency. The locus of the friction force vector shows the sliding phases as circular arcs on the Coulomb cone, connected by two curves which lie inside the cone and denote the long sticking phases. The short sticking states are observable as two small kinks on the
<figure><img src="image_22.png" /><figcaption>**Fig. 19:** Response at the subresonance frequency $\omega = 0.303$, for $\mu N = 8$ and $N_{\phi} = 160$. Displacement and velocity plots indicate mass $m_1$ [blue] and mass $m_2$ [black]. Time histories of displacement and velocity indicate the $\mathbf{i}_1$-direction DOF [solid line] and $\mathbf{i}_2$-direction DOF [dashed line].</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 20:** Response at the subresonance frequency $\omega = 0.115$, for $\mu N = 8$ and $N_{\phi} = 160$. Displacement and velocity plots indicate mass $m_1$ [blue] and mass $m_2$ [black]. Time histories of displacement and velocity indicate the $\mathbf{i}_1$-direction DOF [solid line] and $\mathbf{i}_2$-direction DOF [dashed line].</figcaption></figure>

circular arcs of the force vector locus.

Figure 21 depicts a case of multi-harmonic excitations applied to the linear mass, such that $f_{11}(t) = 20 \cos \omega t$ and $f_{12}(t) = 20 \cos 3\omega t$. As was the case for the single-mass oscillator in Section 4.2, the multi-harmonic forcing results

<figure><img src="image_24.png" /><figcaption>**Fig. 21:** Response at $\omega = 0.115$, for $\mu N = 8$ and $N_{\phi} = 1800$, to multi-frequency excitation $f_{11}(t) = 20 \cos \omega t$ and $f_{12}(t) = 20 \cos 3\omega t$. Displacement and velocity plots indicate mass $m_1$ [blue] and mass $m_2$ [black]. Time histories of displacement and velocity indicate the $\mathbf{i}_1$-direction DOF [solid line] and $\mathbf{i}_2$-direction DOF [dashed line].</figcaption></figure>

in a complex response of the friction mass, with six sticking phases per cycle. Also, pronounced Gibbs oscillations were observed for the friction force near discontinuities. To reduce their extent to an acceptable level, at least 1200 harmonics in time had to be considered. In Figure 21, a very large number of harmonics, $N_{\phi} = 1800$, is taken and the Gibbs oscillations are barely visible. Also note the intricate pattern of sliding phases on the perimeter of the Coulomb cone and the sticking phases inside it. In this complex case as in all other cases considered, Coulomb's law is accurately satisfied by using the E-WR formulation, and no convergence issues were identified.

## 5.3 Validation with Time Integration

To validate the E-WR method for the two-component system, the case in Figure 19 is compared with results obtained using time integration, as depicted in Figure 22. Similarly to the single-mass system, the two sets of results show
<figure><img src="image_25.png" /><figcaption>**Fig. 22:** Comparison of the response of damper mass $m_2$ using the E-WR method with $N_\phi = 160$ and implicit Euler time integration, at $\omega = 0.303$, for $\mu N = 8$ and $N_\phi = 160$. Displacements and velocities: E-WR [red] and Euler [black] in $i_1$ direction; E-WR [green] and Euler [blue] in $i_2$ direction. Friction force: E-WR [orange] and Euler [black].</figcaption></figure>

excellent agreement, except for the small Gibbs oscillations observed near the discontinuities of the friction force in the E-WR method. These oscillations could be confined to a smaller region by considering a larger number of harmonics. The required computation time for the implicit Euler time integration to reach the periodic steady-state was approximately 2100 s, hence considerably larger in comparison to the 3 s it took to obtain the E-WR result. The efficiency of the E-WR method relative to time integration in terms of computational cost is expected to increase with the number of both linear and nonlinear DOFs. This is due to very small time steps required to capture state transitions and the unavailability of dynamic condensation when performing time integrations.

## 5.4 Sensitivity to Damper Mass

As mentioned in Section 1, the assumption of a massless frictional interface is commonly made to account for the flexible nature of contact surfaces, which also allows for the circumvention of technical difficulties associated with set-valued friction laws. Therefore, an analysis of the effect of damper mass $m_2$ on the system response is conducted here.

The E-WR formulation can handle interfaces that either have mass or are massless, in the latter case by simply setting $m_2 = 0$ in Equation (36). Figure 23 shows the responses of the friction mass as $m_2$ decreases from 0.5 to 0.1 to 0.01 to 0. Note that displacement amplitudes decrease as $m_2$ is reduced, due to the smaller inertia. The velocity time histories show an interesting trend: as the damper mass decreases, oscillations in velocity arise during the sliding phases, the frequency of which increases as $m_2$ decreases, culminating in a discontinuous jump in velocity and a corresponding impulse in acceleration for the massless damper ($m_2 = 0$). The jump in velocity occurs because the damper has zero inertia, but its physical validity may be questioned because in reality dampers do possess some amount of mass. In contrast, the friction force is seen to be continuous in time for a massless damper $m_2 = 0$, due to the absence of a discontinuous acceleration term in the damper mass equations of motion. In fact, the jump in the friction force, which is very significant for $m_2 = 0.5$, becomes barely visible for $m_2 = 0.1$ until it disappears for $m_2 = 0$. The plots of the locus of the friction force vector and the trajectory of the damper mass are also provided.

## 6 Conclusion

An equality-based weighted residuals formulation for periodic solutions of systems undergoing two-dimensional frictional occurrences has been presented. The two friction conditions on the two-dimensional contact surface are expressed as two implicit equalities, leading to a primal-dual mixed displacement-force formulation which can be solved using a standard Ritz-Galerkin projection in which the integrals in the nonsmooth friction terms arising from the weak formulation are computed numerically at each iteration of the nonlinear solver. The solution procedure is compact and entirely carried out in the frequency domain, alleviating the need for the time integration of the contact forces at each iteration of the nonlinear solver. The formulation allows for an exact representation of Coulomb's set-valued law for 2D frictional interfaces with mass, without resorting to common hypotheses such as those regarding the direction of motion on the contact surface, regularization, penalization, or massless flexible interfaces. The rate of convergence and accuracy of the method are solely dictated by the family of basis functions $\phi_k$ in the Ritz-Galerkin expansion and the number of these functions.

The new formulation has been shown to capture complex 2D behaviors such as multiple stick-slip transitions and intricate trajectories of the friction point, accommodating very large numbers of harmonics (up to 1800) when required
<figure><img src="image_26.png" /><figcaption>**Fig. 23:** Response of the two-mass system for values of the damper mass $m_2$ ranging from 0 to 0.5. Other system parameters are $f_{11}(t) = f_{12}(t) = 10 \cos \omega t$, $\omega = 0.303$ (or $T = 20.74$), $\mu N = 8$, and $N_{\phi} = 160$.</figcaption></figure>

to capture the friction force accurately, and thus demonstrating its effectiveness. Computational times were found to be reasonable and to decrease significantly by providing the system's quasi-analytical Jacobian to the nonlinear solver. Convergence and validation studies confirm the accuracy and robustness of the methodology and its ability to handle
large numbers of harmonics without numerical difficulties. The formulation has also been applied successfully to a
simple two-component system with a linear mass and a friction-damped mass, which features very rich sub-resonant
periodic responses. It is noteworthy that the E-WR formulation does not appear to be sensitive to the unavoidable
Gibbs phenomenon which arises in the solution when the friction force is a discontinuous function of time.

The equality-based weighted residuals formulation is versatile and, as such, offers a promising methodology for
analyzing complex or large-scale systems with friction and unilateral contact, including those with advanced friction
laws, three-dimensional coupled friction-contact for variable normal loads, complex interfaces with large numbers of
frictional contact points, multiple frictional interfaces, and large-scale structures such as bladed disk assemblies with
friction dampers and various frictional interfaces. The mixed primal-dual formalism of friction also holds promise for
the stability analysis of periodic solutions by computing the eigenvalues of the weak Jacobian of the system.

**Supplementary material**

The numerical code used in this work is shared publicly on a repository [89]. 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.

**References**

[1] V. Sessa, L. Iannelli, F. Vasca, and V. Acary. A complementarity approach for the computation of periodic oscillations in piecewise linear systems. *Nonlinear Dynamics*, 85:1255–1273, 2016. doi:10.1007/s11071-016-2758-5.

[2] M. Miyasaka, M. Haghighipanah, Y. Li, J. Matheson, A. Lewis, et al. Modeling cable-driven robot with hysteresis and cable-pulley network friction. *IEEE/ASME Transactions on Mechatronics*, 25(2):1095–1104, 2020. doi:10.1109/TMECH.2020.2973428.

[3] Y. Q. Sun and C. Cole. Finite element modeling and analysis of friction-wedge dampers during suspension pitch modes. *Journal of Vibration and Acoustics*, 130(2):021003, 2008. doi:10.1115/1.2827369.

[4] F. Bertails-Descoubes, F. Cadoux, G. Daviet, and V. Acary. A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies. *ACM Transactions on Graphics*, 30(1):1–14, 2011. doi:10.1145/1899404.1899410.

[5] J. Li, G. Daviet, R. Narain, F. Bertails-Descoubes, M. Overby, et al. An implicit frictional contact solver for adaptive cloth simulation. *ACM Transactions on Graphics*, 37(4), 2018. doi:10.1145/3197517.3201308.

[6] Y. L. Chen, M. Ly, and C. Wojtan. Primal-dual non-smooth friction for rigid body animation. In *ACM SIGGRAPH 2024 Conference Papers*, SIGGRAPH '24, New York, NY, USA, 2024. Association for Computing Machinery. doi:10.1145/3641519.3657485.

[7] B. Brogliato, F. Miranda-Villatoro, and A. Younes. Backstepping passivity-based trajectory tracking control of frictional oscillators: continuous and discrete time analyses. Preprint, 2024.

[8] J. Das and A. Mallik. Control of friction driven oscillation by time-delayed state feedback. *Journal of Sound and Vibration*, 297 (3):578–594, 2006. doi:10.1016/j.jsp.2006.04.013.

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

[10] R. Rathee. Numerical modeling and simulation of friction models for mechanical systems: A brief review. *Materials Today: Proceedings*, 2023. doi:10.1016/j.matpr.2022.12.185.

[11] Q. Wu, C. Cole, M. Spiryagin, and Y. Q. Sun. A review of dynamics modelling of friction wedge suspensions. *Vehicle System Dynamics*, 52(11):1389–1415, 2014. doi:10.1080/00423114.2014.943249.

[12] 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 *Multibody Dynamics: Computational Methods and Applications*, pages 107–130. Springer, 2011. doi:10.1007/978-90-481-9971-6_6.

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

[14] V. Acary and B. Brogliato. *Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics*. Springer, 2008. doi:10.1007/978-3-540-75392-6.
[15] L. Marino and A. Cicirello. Dynamic response of multi-degree-of-freedom systems with a coulomb friction contact under harmonic excitation. *Nonlinear Dynamics*, 106:1675-1709, 2021. doi:10.1007/s11071-021-06966-3.

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

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

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

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

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

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

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

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

[24] 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*, volume 5, pages 457-467. ASME, 2008. doi:10.1115/gt2008-50860.

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

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

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

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

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

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

[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] L. Woiwode, N. N. Balaji, J. Kappauf, F. Tubita, L. Guillot, et al. 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.

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

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

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

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

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

[43] K. Y. Sanliturk and D. J. Ewins. Modelling two-dimensional friction contact and its application using harmonic balance method. *Journal of Sound and Vibration*, 193(2):511–523, 1996. doi:10.1006/jsvi.1996.0299.

[44] C.-H. Menq, P. Chidamparam, and J. H. Griffin. Friction damping of two-dimensional motion and its application in vibration control. *Journal of Sound and Vibration*, 144(3):427–447, 1991. doi:10.1016/0022-460X(91)90562-X.

[45] B.-D. Yang. *Contact kinematics of friction interfaces and applications to the prediction of resonant response of frictionally constrained turbine blades*. PhD thesis, The Ohio State University, USA, 1996.

[46] P. E. Klauser. Modeling friction wedges: Part i—the state-of-the-art. In *ASME International Mechanical Engineering Congress and Exposition*, number IMECE2004-61094, pages 67–75, 2004. doi:10.1115/IMECE2004-61094.

[47] J. Argüelles, E. Casanova, and M. Asuaje. Harmonic response of a piping system considering pipe-support friction via hft method. In *Pressure Vessels and Piping Conference*, volume 44588, pages 339–347, 2011. doi:10.1115/PVP2011-57172.

[48] C. M. Firrone and S. Zucca. Modelling friction contacts in structural dynamics and its application to turbine bladed disks. In *Numerical Analysis - Theory and Application*, pages 301–334. IntechOpen, 2011. doi:10.5772/25128.

[49] J. H. Griffin and C.-H. Menq. Friction damping of circular motion and its implications to vibration control. *Journal of Vibration and Acoustics*, 113(2):225–229, 1991. doi:10.1115/1.2930173.

[50] C.-H. Menq and B.-D. Yang. Non-linear spring resistance and friction damping of frictional constraint having two-dimensional motion. *Journal of Sound and Vibration*, 217(1):127–143, 1998. doi:10.1006/jsvi.1998.1739.

[51] M. Afzal, I. L. Arteaga, and L. Kari. An analytical calculation of the jacobian matrix for 3D friction contact model applied to turbine blade shroud contact. *Computers & Structures*, 177:204–217, 2016. doi:10.1016/j.compstruc.2016.08.014.

[52] S. Quaeghebeur, 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:111277, 2022. doi:10.1016/j.ijsolstr.2021.111277.

[53] D. Laxalde, F. Thouverez, and J. P. Lombard. Forced response analysis of integrally bladed disks with friction ring dampers. *Journal of Vibrations and Acoustics*, 132(1):011013, 2010. doi:10.1115/1.4000763.

[54] E. Lemoine, D. Nélias, F. Thouverez, and C. Vincent. Influence of fretting wear on bladed disks dynamic analysis. *Tribology International*, 145:106148, 2020. doi:10.1016/j.triboint.2019.106148.

[55] T. Vadcard, F. Thouverez, and A. Batailly. Nonlinear normal modes-related isolated branches of subharmonic solutions for forced response blade-tip/casing contact problems. *Journal of Engineering for Gas Turbines and Power*, 147(5):051022-1-11, 2025. doi:10.1115/1.4066713.

[56] O. Poudou. *Modeling and analysis of the dynamics of dry-friction-damped structural systems*. PhD thesis, University of Michigan, USA, 2007.

[57] D. Charleux. *Etude des effets de la friction en pied d'aube sur la dynamique des roues aubagées*. PhD thesis, Ecole Centrale de Lyon, France, 2006.
[58] B.-D. Yang and C.-H. Menq. Characterization of 3D contact kinematics and prediction of resonant response of structures having 3D frictional constraint. *Journal of Sound and Vibration*, 217(5):909-925, 1998. doi:10.1006/jsvi.1998.1802.

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

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

[61] C. Li, Z. Chen, Z. Shen, and H. She. Analysis on vibration reduction of a rotating dovetailed blade with under-platform damper and installation preload. *Mechanics Based Design of Structures and Machines*, 50(3):877-896, 2022. doi:10.1080/15397734.2020.1732222.

[62] D. Charleux, C. Gibert, F. Thouverez, and J. Dupeux. Numerical and experimental study of friction damping blade attachments of rotating bladed disks. *International Journal of Rotating Machinery*, 2006(1):071302, 2006. doi:10.1155/IJRM/2006/71302.

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

[64] F. Xia. Modelling of a two-dimensional coulomb friction oscillator. *Journal of Sound and Vibration*, 265(5):1063-1074, 2003. doi:10.1016/S0022-460X(02)01444-X.

[65] F. Xia. Modelling of wedge dampers in the presence of two-dimensional dry friction. *Vehicle System Dynamics*, 37(sup1): 565-578, 2002. doi:10.1080/00423114.2002.11666263.

[66] Q. Wu, C. Cole, and M. Spiryagin. Train braking simulation with wheel-rail adhesion model. *Vehicle System Dynamics*, 58(8): 1226-1241, 2020. doi:10.1080/00423114.2019.1645342.

[67] J. R. Barber and X. Wang. Numerical algorithms for two-dimensional dynamic frictional problems. *Tribology International*, 80: 141-146, 2014. doi:10.1016/j.triboint.2014.07.004.

[68] 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_1.

[69] V. Acary, M. Brémond, and O. Huber. On solving contact problems with coulomb friction: formulations and numerical comparisons. In *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.

[70] W. Heemels, V. Sessa, F. Vasca, and M. Camlibel. Computation of periodic solutions in maximal monotone dynamical systems with guaranteed consistency. *Nonlinear Analysis: Hybrid Systems*, 24:100-114, 2017. doi:10.1016/j.nahs.2016.10.006.

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

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

[73] H. B. Dhia and C. Zammali. Level-sets fields, placement and velocity based formulations of contact-impact problems. *International Journal for Numerical Methods in Engineering*, 69(13):2711-2735, 2007. doi:10.1002/nme.1860.

[74] Y. Renard. Generalized newton's methods for the approximation and resolution of frictional contact problems in elasticity. *Computer Methods in Applied Mechanics and Engineering*, 256:38-55, 2013. doi:10.1016/j.cma.2012.12.008.

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

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

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

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

[79] 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%3C1001 AID-NME319%3E3.0.CO 2-A.
[80] 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.

[81] M. Legrand and C. Pierre. A compact, equality-based weighted residual formulation for periodic solutions of systems undergoing frictional occurrences. *Journal of Structural Dynamics*, 2:144-170, 2024. doi:10.25518/2684-6500.190.

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

[83] M. Colin, O. Thomas, S. Grondel, and E. Cattan. Very large amplitude vibrations of flexible structures: experimental identification and validation of a quadratic drag damping model. *Journal of Fluids and Structures*, 97:103056, 2020. doi:10.1016/j.jfluidstructs.2020.103056.

[84] B. Cochelin, D. Noureddine, and M. Potier-Ferry. *Méthode asymptotique numérique*. Hermès-Lavoisier, 2007. ISBN 9782746215672.

[85] M. Assidi, H. Zahrouni, N. Damil, and M. Potier-Ferry. Regularization and perturbation technique to solve plasticity problems. *International Journal of Material Forming*, 2(1):1-14, 2008. doi:10.1007/s12289-008-0389-1.

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

[87] L. Salles. *Etude de l'usure par fretting sous chargements dynamiques dans les interfaces frottantes : application aux pieds d'aubes de turbomachines*. PhD thesis, Ecole Centrale de Lyon; Bauman Moscow State Technical University, France/Russia, 2010.

[88] B.-D. Yang, J. J. Chen, and C.-H. Menq. Prediction of resonant response of shrouded blades with three-dimensional shroud constraint. *Journal of Engineering for Gas Turbines and Power*, 1999. doi:10.1115/1.2818504.

[89] Equal[iser]: A frequency-domain solver for the vibration analysis of structural mechanical systems subject to unilateral and frictional contact occurrences. 2025, URL https://archive.softwareheritage.org/swh:1:dir:376178718cf4bb2fc8d9ce0ccaa10534182679ab.