## 1 Introduction  
Frictional interfaces can be found in numerous aspects of everyday life: brakes in the automotive sector [1], joints in most engineering structures [2], friction damping devices in turbines (for energy production and aircraft engines alike) [3], and so on. Friction damping devices on turbine disks, either as integral parts of the disk (e.g., shrouds) [4; 5; 6] or as separate components [7; 8; 9; 10; 11; 12; 13], are considered strategic in the field, as, if well designed, they can play a key role in avoiding high cycle fatigue failures caused by large resonant stresses [2; 3; 14; 15]. A trustworthy design of structures with friction interfaces entails finding:  
1. a reliable model for the force-displacement relation at the contact interface;
2. an efficient yet accurate technique to solve the nonlinear equilibrium equations that obtain the nonlinear response in a timely manner, compatible with design purposes. Point 1 has been addressed by several authors, leading to different contact models and techniques. Some authors adopt a dynamic Lagrangian method to solve the contact patch [16; 17], i.e., the contact constraints are taken into account in their non-regularized form without additional compliance. Other authors, e.g., [12; 13; 18; 19; 20], apply a \"spring-slider\" contact element [3; 21; 22; 23] to each meshed node belonging to the contact area, introducing normal and tangential stiffnesses and a Coulomb friction law. In detail, the contact interactions can be geometrically divided in normal and tangential directions. A unilateral contact law is often considered in the normal direction (with or without normal contact stiffness \(k_{n}\)), and frictional law for the tangential contact. This last method is here preferred, as its calibration parameters (\(k_{n}\), \(k_{t}\) and \(\mu\)), however difficult to determine, represent a physical measurable property [24; 25; 26; 27; 28; 29]. These contact elements belong to the larger family of heuristic models, as opposed to microscale \"realistic\" models, where asperities and surface roughness are modeled using stochastic distributions [30]. Heuristic models are here preferred because of their simplicity, because of the ease of calibration, and, most importantly, because they have proven their adequacy in several experimental/numerical comparisons [31; 32]. In detail, this paper makes use of a 3D contact model [18; 22], most often applied because of its versatility. It should be noted that all the techniques here developed can be applied to any contact model belonging to the \"spring-slider\" family with minor adjustments. As for Point 2, the presence of friction-induced nonlinearities makes standard solution techniques (typically applied in commercial FE codes) inadequate. Whenever possible, the response is assumed to be periodic so that harmonic balance methods (HBMs) can be used in the frequency domain to compute the steady-state response [33; 34; 35; 36; 37]. It should be noted that the HBM is an approximate method [38]. In linear systems, the harmonic support is easily determined by the harmonic content of the forcing function (e.g., mono-harmonic). However, if nonlinearities are present, higher-order harmonics may need to be included to ensure a proper representation of friction forces and displacements, with a detrimental effect on the the computational time. Most importantly, the HBM cannot be applied effectively in all those cases where the response is non-periodic or the transient is of interest. This paper focuses on time marching methods, an alternative to HBM. Time marching methods start with an initial solution (defined at time zero) and march the equations in time while applying boundary conditions. Therefore, the solution is dependent on the solutions at previous times. An example of a time marching method is any direct time integration (DTI) algorithm, but additional alternatives will be explored within this work.