Competitive Time Marching Solution Methods for Systems with Friction-Induced Nonlinearities

: Finding efﬁcient and accurate solution methods for nonlinear equilibrium equations is a challenging task. This is the case of systems with friction-induced nonlinearities, e.g., friction-damped turbomachinery assemblies and automotive applications such as brakes. In order to tackle this strategic task, several methods have been developed, both in the time and in the frequency domains. Time marching methods are regarded as the most accurate option, but their computational cost becomes prohibitive when friction nonlinearities are present. This poses a problem in all those cases where alternative frequency domain methods cannot be applied effectively, e.g., if transients, non-periodic excitation/solution, or highly nonlinear systems are of interest. The purpose of this paper is to propose three independent methods to make time-marching more competitive. Two of these methods can be applied to any existing direct integration scheme with minimal adjustments, but the computational time cut they introduce is signiﬁcant. The last method is instead tailored for systems where the inertia force contribution is negligible. All methods are thoroughly validated numerically using a standard Newmark- β integration scheme as a reference.


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 µ), 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.The purpose of this paper is to provide a series of methods to make time marching methods competitive and practically applicable for design purposes, a result of strategic importance in all those cases where frequency-domain solution methods cannot be applied.This task will be carried out in a series of substeps outlined below.

•
Section 2 will recapitulate the dynamic governing equations and will present the test case used as a demonstrator throughout this article.• Section 3 will recapitulate the main features of the contact model and present an alternative way to represent contact forces, as a sequence of linear states.

•
Section 4 will use the contact force representation from Section 3 to propose two methods applicable to any DTI algorithm.• Section 5 will use the contact force representation from Section 3 to propose a quasi-static time marching algorithm, tailored for systems where inertial effects are minimal, capable of reducing standard computational times of one order of magnitude.
All proposed methods will be validated and time savings will be quantified using a reference test case as a benchmark.Their applicability to different fields and situations will be discussed.

Governing Equations
All multi-degree of freedom (multi-DOF) systems subjected to an external dynamic forcing and to friction share similar equilibrium equations.In the time domain, they can be written as the dependence of the contact forces vector f C from the vector of displacements q is defined by the contact model [18,22], further discussed below.
Equation ( 1) can be tailored to represent an underplatform damper inserted between two adjacent turbine blades.The dissipative function of underplatform dampers has been discussed in the introduction, while Figure 1 shows a graphical representation of a turbine disk with underplatform dampers.The numerical code concerned with the solution of Equation ( 1), already applied and validated in [38,39], can be considered as a "plug-in" software used to represent any given underplatform damper between blades modeled using finite elements (see also Figure 2a).The damper is here considered as a rigid body.The rigid body assumption is well justified by the bulkiness of "solid-bar" underplatform dampers and has already been validated against specific experimental evidence on underplatform dampers and blades [23,31,32].The same conclusion has been drawn for similar applications in [4].Furthermore, this choice has two advantages: (1) damper mass, shape, and contact points are easily editable by the user, and (2) the number of physical damper DOFs is limited, regardless of the selected number of contact points.
In this paper, the system is further simplified, to be easily used as a demonstrator.A 2D version of the damper between a set of blades is considered, as shown in Figure 2b.This version of the damper system has been effectively used and validated in a series of previous works [19,39,40], where the blade platforms undergo bending (i.e., pure in plane) modes of vibration.Furthermore, it should be noted that the techniques described throughout Sections 3-5 can be easily extended to more complex systems.Their formulation is general.With reference to Figure 2b, Equation (1) becomes where M D = diag(m D , m D , I z ) is the mass matrix, f ED = (0, −|CF|, 0) is the vector of external forces on the damper, q D = (u D , w D , fi D ) is the vector of unknowns at the damper center of mass, and q P = (q L1P , q L2P , q RP ) is the vector of input platform displacements, where q iP = (u iP , w iP ) and i = L1, L2, R. The vector q P is here seen as an input because it comes from the displacement of the nodes on the blade platforms, provided by "higher-level" software [38] concerned with the blades' dynamics.The indication of the time instant τ has been dropped for the sake of brevity.
In this particular set of equations, there are no harmonic "external forces."The input platform displacements act on the contact elements and produce the contact forces that serve as the harmonic external excitation to the system.Further details on the contact model implementation are given in Section 3. The vector of contact forces f CD refers to the contact forces transposed at the damper center of mass: f CD is expressed in the global reference system, while f * i = (T i , N i ) , q * iP = (t iP , n iP ) , and q * iD = (t iD , n iD ) are expressed in the local reference system (in red in Figure 2b), specific for each contact point i.
Transformation matrices are used to switch between reference systems: where Γ iP is a 2 × 2 rotation matrix, and Γ iD is a 2 × 3 transformation matrix.

The Contact Model and Its Piecewise Linearity
The contact forces at a particular instant in time and at a given contact node pair f * i are computed starting from the relative displacements at the contact q * iP − q * iD using a friction contact model with normal load variation that is widely used in the nonlinear forced response analysis of frictional constraint structures [18,22], see also Figure 3.The time-discrete implementation facilitates the contact model application in time-marching solution methods.The main drawback is that contact state transitions are captured with finite accuracy, i.e., the length of the chosen time step ∆τ.This has, however, negligible effects on the solution for the range of investigated ∆τ. Figure 3 shows the contact element in case of 2D motion.In case of 3D motion, e.g., see Figure 2a, it is assumed that the in-plane tangential components of the contact force are independent from each other (i.e., quasi-3D).In other words, an uncoupled in-plane motion/force coordinate, here termed z − Z, orthogonal to t − T, is added.This assumption is commonly used in the calculation of the nonlinear forced response of structures with frictional constraints, see [20,[41][42][43] Each contact element can undergo, alternatively, three different contact states.
• Lift-off: There is a gap between the platform and damper contact points (i.e., n iP − n iD ≤ 0), and contact forces are null T = N = 0. • Stick: There is interference between the platform and damper contact points (i.e., n iP − n iD > 0).Furthermore, tangential relative displacement is limited, and a proportional relation between tangential forces and displacements exists: , where s i is the relative displacement of the tip of the slider with respect to the contact point on the damper.

•
Slip: There is interference between the platform and damper contact points (i.e., n iP − n iD > 0).Furthermore, the tangential relative displacement exceeds a threshold value, and a proportional relation between tangential and normal forces exists: The full set of conditions for the three contact states described above can be found in [11].
The standard implementation of this contact model sees a "black-box-like" subroutine, which receives as input the displacement vectors q * iP , q * iD at the time step τ and the slider displacement s i at time τ − ∆τ and gives, as output, the contact forces f * i = (T i , N, i) .Inside the sub-routine a predictor-corrector scheme (see Figure 3) checks the inputs against a series of "conditions" in order to assign the correct contact state and the corresponding values of contact forces.This paper proposes a novel alternative to implementing the contact model constitutive equations.The contact state transition conditions and contact forces expressions are unchanged, so the ultimate result is the same.However, this novel formulation brings out an important feature of the chosen kind of contact model: nonlinearity is obtained through a time sequence of linear systems.The nonlinear system is actually a piecewise-linear system.
With reference to Equation (3), let us consider a specific instant in time.The local vector of contact forces f * i (τ) is substituted by the product of the contact stiffness matrix K i (τ) and the vector of local relative displacement at the contact q * iP (τ) − q * iD (τ), plus a contribution from the slider.
The local contact stiffness matrix K i and the slider contribution K Si change their configuration with the contact state, as described in Table 1.Therefore, if one considers a time interval ∆τ during which the contact state of all the contact points remains unchanged, the system is perfectly linear.If the reader compares this formulation with the standard constitutive equations of the contact model, he/she will find the same ultimate result.
The presence of the contact forces, vector f CD , in Equation ( 2) can thus be substituted with a set of modified "effective" stiffness matrices.Considering Equations ( 2) and (3), Then, considering Equations (4) and 5, Table 1.Expression of the contact stiffness matrix depending on the contact state.

Contact State Condition
This formulation will be used throughout Sections 4-5 to ensure an optimal determination of the integration time step and as a basis for the development of novel time-marching solution methods.

Direct Time Integration Tailored to Systems with Friction Nonlinearities
In numerical analysis, numerical integration constitutes a broad family of algorithms for calculating the numerical value of a definite integral.By extension, the term is also used to describe the numerical solution of differential equations.
These methods rely on discretizing the continuous problem, i.e., ordinary differential equations (ODEs) such as Equation ( 1), since only rarely can it be solved analytically.
This involves, after the solution is defined at time zero (i.e., initial conditions), the attempt to satisfy dynamic equilibrium at discrete points in time.Most methods use equal time intervals at ∆τ, 2∆τ, ..., N∆τ.However, this is not mandatory; in some cases, a variable time step might be employed.An example will be given in the following paragraphs.

Selection of the Direct Time Integration Method
An overwhelmingly large number of algorithms has been proposed over the years with varying degrees of complexity, accuracy, and computational cost.It is sometimes difficult to understand which algorithms are better suited to a given purpose.Appendix A.1 presents a classification of all families of DTI methods available in the literature.The classification highlights the aspects relevant to the problem at hand (the solution of 2nd order ODEs), thus allowing for an informed choice.
The algorithm used in this work is the Newmark-β method, i.e., the constant average acceleration method [44].It was chosen for the following reasons: • it is reasonably accurate, especially for small ∆τ, and does not introduce spurious damping and period elongation; • it is suited for 2nd order ODEs, and no further processing is required; • the method is implicit, which has a positive impact on the required time step size ∆τ and on stability.

Adaptive Time Step Selection for Direct Time Integration Methods
As stated in [44], stability, convergence, and accuracy all depend on the size of the time step ∆τ.Specifically, the chosen Newmark-β method has the ∆τ/T < 0.318 limitation to ensure convergence.The quantity T is defined as the shortest natural period of the structural system, which is computed through a standard eigenvalue analysis.

Remark 1.
When computing T, one should not consider the "structural" stiffness of the system but rather the complete effective stiffness, contact springs included (see Section 3).In other words, the value of T changes with the contact state.This aspect is not stressed in the few publications where the DTI of friction-damped turbo-machinery systems is carried out.One way to define a value of ∆τ valid for all instants of the period of vibration is to use the highest effective stiffness that might be encountered, i.e., a full-stick state.With reference to Equation (2), the eigenvalue problem from which T is defined is where all K i (∀i) matrices have the "stick" configuration shown in Table 1.In this case, the structural stiffness matrix K D is not present because the damper is modeled as a rigid body.If it were, the overall stiffness matrix would be Typical values of the contact springs, i.e., ≈ 10 − 10 2 N/µm, unfortunately produce ∆τ ≈ 1 µs.This makes DTI cumbersome, especially at low frequencies.
The reader will have noticed that the "effective stiffness" of the system changes with time (i.e., with the contact state).Therefore, the limit on ∆τ depends on the specific contact state as well.For this reason, it is here proposed, as shown in Figure 4, that the time step is adapted to the contact state.The contact state is initially assumed to be equal to the previous one, i.e., ∆τ = τ n − τ n−1 : this value is fed into the integration scheme and produces the tentative value of the vector of displacements q D .The contact model run at instant τ n with q D produces a given contact condition and, consequently, the time step ∆τ.If ∆τ = ∆τ, the algorithm goes on; otherwise, the time step is updated; if ∆τ < ∆τ, the Newmark algorithm is re-executed with the updated value of ∆τ.There are typically a few transitions per period of vibration (i.e., 4-7).The time savings are quite consistent: the signal shown in Figure 4b is in full slip condition for 54% of the total period of vibration.If it is run with a unique time step compatible with the full stick condition, it will take 187 time steps to complete.The adaptive algorithm reduces this number to 148.The "relative" time step savings in the investigated cases are in the 20-40% range compared with the standard condition.These values depend on the amount of gross slip/contact point detachment and on the values of contact parameters.These results have been obtained using parameters experimentally estimated in [29].The adaptive time step can be applied to any DTI algorithm (e.g., Runge-Kutta and Euler) since all methods have limitations on ∆τ/T.

A Novel "Piecewise-Linear" Implementation for Direct Time Integration Methods
Figure 4a clearly shows that, since Newmark-β is an implicit scheme and the system at hand is nonlinear, an initial assumption on one of the unknown vectors at time τ n+1 , e.g., n+1 q D is required.The first version of the damper model presented in [11] assumes n+1 q D = n q D .In other words, it assumes that the displacement keeps constant over a period ∆τ.This assumption is not relied upon, as shown in Figure 4a and Appendix A.2: a predictor-corrector scheme (different from the one previously mentioned for the contact model) iterates until the assumed value is verified by the equilibrium equations.This assumption results in an average of ≈5 iterations per time step for frequencies <100 Hz and reaches as high as ≈11 for frequencies ≥200 Hz (see Figure 5a).Therefore, the impact on the computational time is quite negative and non-negligible.One may argue that better assumptions should be tried out, e.g., using the result of an explicit method as a starting guess.Many options have been tried out and the novel method proposed below was found to be the most suitable for our purpose: the number of iterations per time step is minimized thanks to an "ad-hoc" predictor-corrector.
The novel initial assumption (i.e., predictor) is based on two strategic observations: • the contact state changes only a few times per period of vibration (4-7 in cases encountered in this work); • if the contact state does not change from time step τ n to τ n+1 , the system can be considered perfectly linear (see Section 3) This novel formulation manages to derive a tentative value of all unknown vectors using the piecewise linearity of the contact model (Section 3), i.e., setting K i (τ n+1 ) = K i (τ n ) ∀i.If the assumption is verified (i.e., the contact state is unchanged), as is the case in >90% of cases, only one iteration is required for each step.If it is not, then the procedure is repeated using, as an assumption, the updated K i matrices corresponding to the forecast contact state.Typically, only one additional iteration is required (see Figure 5b).The complete formulation is reported in Appendix A.2.In other words, the Newmark-β method is applied as if on a linear system for the entire duration of a given contact state (therefore, no iterations are required as described in [45]).Iterations are necessary only when a contact state transition occurs.
The overall computational time is reduced by a factor of almost 50%, e.g., for a 100 Hz 3-period DTI run, the "standard" algorithm's the total time is 9.46 s, compared to the 5.09 s of the "Piecewise linear" implementation.Similar time savings have been found using different DTI integration schemes.Furthermore, it should be noted that this piecewise-linear implementation can be applied together with the adaptive time step to further increase time savings.

A Time-Marching Method for Systems with Negligible Inertial Effects
The concept of piecewise linearity of Coulomb friction contacts has secured significant improvements in the field of DTI algorithms and their implementation.Despite the improvements detailed in Section 4, the limitations on the value of the time step ∆τ = O(10 −6 ) s may still be forbidding in many applications.
It will be shown below how it is possible to reduce the computational time by more than one order of magnitude in all those cases where inertial effects are negligible.A typical example is an underplatform damper between a set of blades or a flexible strip damper [46,47] in the same condition.Alternatively, brake systems where decelerations produce forces that are negligible with respect to friction forces can be studied using the same algorithm.This paper will focus on the first example to develop a novel quasi-static solution method.

The Influence of Damper Dynamics on the Overall Damper Equilibrium
The equations solved through DTI are dynamic equations; i.e., they take into account inertial and viscous damping effects.However, given the reduced mass of the damper and the lack of a damping matrix C D , it is worth investigating what role the dynamics inside Equation ( 2) play and what the consequences of neglecting it are.Furthermore, in many applications, the presence of a viscous damping matrix can be neglected as the damping provided by friction is far larger.To properly answer the questions posed above, two aspects have to be considered: 1.
the magnitude of the damper inertia forces; 2.
the damper "inner resonance" (i.e., the natural frequencies of the damper oscillating on the contact springs).
To tackle the first point, let us consider a damper of mass m D .Both the centrifugal force (and therefore the contact forces) and the inertial effects are proportional to the mass of the damper; specifically, the centrifugal force is CF = m D • R D • Ω 2 , and the inertial effects are where Ω is the rotation speed of the bladed disk, and ω is the frequency of vibration, i.e., close to the mode under investigation ω B .Let us define a ratio between the two quantities: It should be noted that this ratio does not depend on the mass of the damper.Let us consider a damper mounted on a turbine for energy production; it is therefore reasonable to set Ω = 3600 rpm and R D ≈ 0.5 m.Based on the authors' experience, the magnitude of the damper displacement is at most 10 −5 − 10 −4 m.In this specific case, the frequency of vibration ω ≈ ω B (i.e., the mode being excited) would have to be >30 kHz in order to have the inertial forces at least 5% of the centrifugal force (i.e., r I/CF = 0.05).Consequently, the damper vibrating with the first modes of a bladed system, typically in the range of 150-600 Hz, will not be affected by the inertia forces.This aspect allows for the safe reconstruction of the damper equilibrium starting from experimental data at low frequencies [11,19,29] and at higher realistic frequencies (i.e., f B = 2πω B = 200-500 Hz) [39,40].Regarding the damper inner dynamics, the presence of the contact springs causes high-frequency oscillations of the damper: they are visible as superposed waves in Figure 6.The frequency of oscillation ω D , given by the mass of the damper (12 g in this example) and the values of contact springs (determined in ), is, in this case, ≈10 kHz.These oscillations are perfectly explained and predictable from a purely numerical point of view.However, they have not been encountered experimentally at either low [29] or higher frequencies [38][39][40].For this reason, they must be removed "a posteriori" in order to make results more legible.This can be done only when the damper internal resonant frequency is far from the mode of vibration of the blades ω D >> ω B , i.e., for ω B ≤ 1 kHz.

Remark 2.
Both ω B and ω D refer to modes of vibration of the blade/damper coupled system, but at ω = ω B the blade is vibrating (e.g., bending mode) and "driving" the damper with it.At ω = ω D , the blade itself may be far from its resonance (low mobility), and the damper is vibrating on the contact springs (local mobility mode).
If ω D >> ω B , the oscillations mentioned above can be filtered out thorough a low-pass filter (see Figure 6).Remark 3.These numerical oscillations cannot be removed safely if ω ≈ ω D , or, more generally, if a mode of vibrating high mobility for the blade is very close (or ≈ coincident) to a high mobility for the damper.Hysteresis cycles appear strongly distorted and hard to interpret: further numerical and experimental investigation on damper-blade dynamic coupling is needed.Since the frequency range typically investigated is ω B < 600 Hz, the damper acts only as a "spring-coupling" to the bladed system, so its dynamics can be safely neglected.

Quasi-Static Solution Method
The negligibility of the dynamic portion of the damper equilibrium equations is confirmed by the results of a newly developed solution method, here termed "quasi-static."In essence, the method, whose implementation algorithm is reported in Appendix A.2, solves a modified version of an equilibrium equation given above (Equation ( 7)): The strict limitations on the value of ∆τ (≈ 10 −6 s) in DTI algorithms are caused by the damper internal resonance, i.e., the stability and convergence of the algorithm can be ensured only if these high-frequency oscillations are correctly captured.However, these oscillations, which are not observed experimentally and may be an artifact due to an imperfect method chosen to model the contact, are of no interest to the designer, especially in the early design stage.
If one choses to neglect the inertial terms from the equilibrium equation, then

•
the spurious numerical oscillations will disappear (see Figure 7), and • inertial forces will cease to participate in the equilibrium, with negligible effects up to The time step used in this time-marching method is limited only by the precision the designer wants to obtain; i.e., if the user decides to discretize a 2π ω B period of vibration into 300 time steps, then a given contact state transition may be shifted by, at most, ±∆τ = ± 2π 300ω B .Typical force and kinematic results of the quasi-static method are shown along with DTI results in Figure 7.

A Critical Comparison of the Proposed Methods
This section is concerned with offering the user practical guidelines and an indication of expected results when one of the three methods proposed in this paper is applied.In order to achieve this goal, a test case has been defined.The test case at hand sees an underplatform damper between a set of blade platforms (as shown in Figure 2b) vibrating at 200 Hz for one period of vibration.The blade platform kinematics is large enough to cause the damper to reach the bilateral gross slip condition (both the right and left interfaces are sliding).This results in six changes in contact states during one period of vibration.The resulting hysteresis cycle at the left (flat-on-flat) contact interface is shown in Figure 6.
The test case defined above has been solved using the following four methods: 1. a standard Newmark-β algorithm, as described in Appendix A.3, to be used as a reference case, which applies standard techniques and no particular adaptation to friction damped systems; 2.
a standard Newmark-β algorithm with an adaptive time step (see Section 4.2); 3.
a piecewise-linear implementation of the Newmark-β algorithm, as described in Section 4.3 and Appendix A.4; 4.
a quasi-static solution method, as described in Section 5 and Appendix B with a user-defined 300 time steps per period.
All methods yield the same set of results; forces and displacements are perfectly superposed.Accuracy is not impaired in the slightest by the proposed methods.Computational time, on the other hand, varies from case to case.The computational time needed for the case described in Method 1 is used as a reference to normalize other values.
As shown in Figure 8 all proposed methods lead to a non-negligible decrease in the computational time.The adaptive time step method (Method 2) introduces a 20% reduction in the computational time and can be applied to ALL existing DTI algorithms.It enables the user to automatically set the largest possible time step that still ensures convergence.Method 3 is instead applicable to all DTI methods and no restrictions on the test case are observed.It is especially useful for implicit algorithms (in this case, Newmark-β) since it enables the user to make the best possible assumption when proceeding to the following time step (i.e., unchanged contact state) without unnecessary iterations.This leads to halved computational times.Methods 2 and 3 can be applied together, thus leading to an overall reduction in the computational time of The reader will notice that the quasi-static method (Method 4) reduces the overall computational time by more than one order of magnitude.This solution method, which is an alternative to DTI algorithms (and therefore to Methods 2 and 3) is, however, applicable only in those cases where inertial effects are negligible.When applicable, it is a very competitive alternative, especially in all those cases where displacements are large (leading to highly nonlinear systems) and the HBM needs large harmonic support, resulting in convergence issues.

Conclusions
Time marching solution methods are often regarded, especially in the turbomachinery field, as too cumbersome and therefore inapplicable to real industrial problems.This is especially true whenever nonlinearities, such as those induced by friction, are present.This paper proves that ad-hoc techniques tailored for the specific nonlinearity can make time-marching methods competitive.
All of them are based on one common idea: seeing the nonlinearity induced by friction as a time sequence of linear systems.As a result, only time steps corresponding to contact state transitions are truly non-linear.This idea has been applied to: • the development of an adaptive time step selection for DTI algorithms resulting in a 20% decrease in the computational time; • the development of a novel "piecewise-linear" formulation for the Newmark-β scheme, adaptable to other DTI algorithms, and resulting in a 50% decrease in the computational time; • the development of a novel quasi-static method to solve equilibrium equations where inertial effects are negligible, capable of reducing computational times by more than one order of magnitude.
The appendices to this paper offer a critical review of existing DTI algorithms that guide the reader through the selection of the most appropriate one for the specific test case at hand.Furthermore, 3.
The system of equations composed of Equation (A13) and the constitutive Newmark relations Equation (A7) is solved.Below is the expression for n+1 qD : The expressions for n+1 q D and n+1 qD are computed using Newmark's equations (Equation (A7)).

4.
Vector n+1 q D is transformed into local displacements at the contact n+1 q * iD using Equation ( 4).The same is done with the displacements of the platform contact points, thus obtaining q * iP (τ n + ∆τ) ∀i.

5.
The vectors of local displacements at the contacts (of both the platforms and the dampers), together with the slider variable s i are fed into a function (similar to the contact element function) that checks whether their values are compatible with the postulated contact condition.In other words, for each i, whether the structure of n+1 K i and n+1 K Si is compatible with the values of local platforms and damper displacements is checked (see Table 1).6.
CORRECTOR: If the compatibility is verified, then the equilibrium postulated in Equation (A13) is indeed satisfied, and the integration goes to τ n+2 .In this case, all slider variables are updated according to the verified contact condition ( n+1 s i , to be used at time step τ n+2 ).
Otherwise, the structure of all effective contact stiffness matrices is updated according to the indications of the function mentioned at Step 5-see also Table 1: n+1 Ki = F n+1 q * iD , n+1 q * iP , n s i n+1 KSi = F S n+1 q * iD , n+1 q * iP , n s i .

(A16)
The algorithm is then repeated starting from Step 3.

Appendix B. Quasi-Static Method Algorithm
The following makes reference to an equilibrium equation given above (Equation ( 10)): 1.
The values of n q D , n qD , n qD , f ED (τ n+1 ), q iP (τ n+1 ) ∀i, and the contact states' effective stiffness matrices n K i and n K Si are known either from initial conditions or the previous integration step.

2.
PREDICTOR: It is assumed that the contact state at τ n+1 is the same as that of the previous time step, i.e., It is therefore possible to compute the expression of n+1 q D using the tentative version of the equilibrium equations at time τ n+1 (see Equation ( 10)): where Vector n+1 q D is transformed into local displacements at the contact n+1 q * iD using Equation ( 4).The same is done with the displacements of the platform contact points, thus obtaining q * iP (τ n + ∆τ) ∀i.

5.
The vectors of local displacements at the contacts (of both the platforms and the damper), together with the slider variable s i are fed into a function (similar to the contact element function) that checks whether their values are compatible with the postulated contact condition.In other words, for each i, whether the structure of n+1 K i and n+1 K Si is compatible with the values of local platforms and damper displacements is checked (see Table 1).6.
CORRECTOR: If the compatibility is verified, then the equilibrium postulated in Equation ( 10) is indeed satisfied, and the integration goes to τ n+2 .In this case, all slider variables are updated according to the verified contact condition ( n+1 s i , to be used at time step τ n+2 ).
Otherwise, the structure of all effective contact stiffness matrices is updated according to the indications of the function mentioned at Step 5-see also Table 1: n+1 Ki = F n+1 q * iD , n+1 q * iP , n s i n+1 KSi = F S n+1 q * iD , n+1 q * iP , n s i

Figure 1 .
Figure 1.(a) CAD model of a turbine bladed disk with underplatform dampers (not directly visible in the figure).(b) Sketch representing a close up of three blades with interposed underplatform dampers.

Figure 2 .
Figure 2. (a) Contact model with three contact springs.(b) Contact patches on the blade platforms and on the strip damper.

Figure 3 .
Figure 3. (Left) Contact model representation and list of its calibration parameters.(Right) Standard implementation of the contact model routine through a predictor-corrector scheme.

Figure 4 .
Figure 4. (a) Scheme representing the Newmark-β implicit DTI scheme (standard implementation) with an adaptive time step.(b) Example of result of Newmark-β DTI with adaptive time step.

Figure 5 .
Figure 5. Histogram representing the number of iterations per time step of the predictor-corrector scheme implemented in the DTI algorithm: (a) constant displacement assumption; (b) unchanged contact state assumption.

Figure 6 .
Figure 6.(a) Low frequency hysteresis cycle at the left contact obtained through DTI without additional post-processing.(b) Low frequency hysteresis cycle at the left contact obtained through DTI and an a posteriori low-pass filter.

Figure 7 .
Figure 7. DTI vs. quasi-static results comparison.(a) Tangential displacement of the right damper contact point.(b) Vertical component of the right contact force.

Figure 8 .
Figure 8.Time savings introduced by the novel techniques and solution methods proposed in the paper.A standard Newmark-β computation is taken as a reference.