Coulomb friction effect on the forced vibration of damped mass-spring systems

This paper aims at assessing the effect of dry friction on the dynamic behaviour of a damped mechanical system subject to harmonic forcing. Previous work on friction damped systems highlighted that not including other forms of damping in the dynamic analysis can lead to unrealistic results such as the presence of infinite resonances. In this contribution, an exact solution is derived for the continuous steady-state response of multi-degree-of-freedom systems with a contact between one of the masses and an external wall, using Coulomb’s law to model the friction force and a modal damping model to account for the system’s damping. Closed-form expressions are also derived for the amplitude and phase of the continuous responses, while stick–slip responses are investigated by using an ad-hoc numerical approach. In addition, analytical and numerical results are used for exploring the features and the motion regimes of the dynamic response, leading to the following conclusions: (i) system’s damping has a limited effect on low-and high-frequency behaviours, on the presence of invariant points and inversions across the transmissibility curves and can therefore be neglected, in non-resonant conditions, in the analysis of structures where dry friction is the main source of dissipation; (ii) when the damping of the system is accounted for in the mechanical models along with Coulomb damping, finite resonant peaks are also obtained in continuous sliding regime and their amplitude decreases linearly with the friction force generated in the contact.


Introduction
The presence of dry friction in vibrating structures can be associated to either beneficial or detrimental effects.In fact, on the one hand, friction is known to cause problems such as noise, wear, loss of efficiency and even more serious consequences such as structural damage and component failures [1,2].However, on the other hand, friction can improve the performance of engineering structures by dissipating energy and reducing vibration amplitudes [3,4].For example, friction dampers are commonly used in turbomachinery [5,6], civil engineering [7,8], suspension systems [9,10], robotic devices [11,12], energy harvesters [13], airplane taxing systems [14] and satellites [15].Nonetheless, even in structures relying on friction dampers as the main source of dissipation, the overall damping usually results from multiple causes.These can include other localised sources of dissipation, such as viscous dampers or shock absorbers, or the general dissipation occurring, for instance, internally in the material or in loose joints.In particular, the latter form of damping is present, to some extent, in all real structures.Therefore, it is essential to understand how the damping levels of a structure can affect the dynamic performance of a friction damper.
Independently of the presence of other forms of damping, the dynamic analysis of friction damped systems is complicated by the nonlinear nature of the friction forces.Even when the simplest and most widely used model of dry friction, the Coulomb model, is considered, exact solutions for the forced vibration problem are difficult to obtain and are only available in a limited number of cases.Different techniques, including time integration, harmonic and multi-harmonic balance methods, have been proposed during the years to deal with complex systems (see review [16] for further details).However, the computational costs related to numerical approaches make them impractical for purposes such as parameter space investigations [17].Even high-order harmonic balance methods, which are extremely popular in the analysis of friction damped systems and offer a good approximation of their steady-state response [16], operate in the frequency domain and cannot easily be used to predict the motion regimes (continuous, stick-slip, permanent sticking) of the dynamic response across the parameter space.Therefore, the research of analytical solutions is of paramount importance for performing parameter space investigations and optimisations with a reduced computational cost, as well as for exploring the motion regimes and the response features of friction damped systems.
In some studies, closed-form solutions have been derived for the steady-state forced vibration of Coulomb friction oscillators assuming dry friction as the only source of damping.In particular, Den Hartog [18,19] and Hong and Liu [20,21] proposed exact solutions for the response of single-degree-of-freedom (SDOF) mass-spring systems with a contact between the mass and a fixed wall.More recently, analytical solutions have also been proposed by the authors for SDOF systems with a contact between the mass and an oscillating wall [22] and multi-degree-of-freedom (MDOF) systems with a friction contact [23,24].Further investigations on the features and stability of the response of Coulomb friction oscillator were also presented in [25,26].Experimental investigations carried out on single-storey [27] and two-storey [28] frame setups have shown that most of the behaviours described by these works can also be observed in laboratory conditions.Nonetheless, assuming Coulomb friction as the only source of dissipation in a vibrating system could lead to unrealistic results.In particular, it has been observed that Coulomb damped systems always present infinite resonant peaks, unless stick-slip or permanent sticking occur at resonance [19,24,25].Since one of the main goals of the use of friction dampers is to avoid large stresses at resonance [5], it is necessary to further investigate the resonant behaviour of these systems, also accounting for the damping levels of the structure in the mechanical model.
The dynamic behaviour of SDOF systems with mixed damping has been investigated by several authors throughout the years.In particular, since the exact nature of general dissipation in real structures is often unknown [29], the viscous damping has usually been considered in combination with Coulomb damping in the investigation of the SDOF behaviour.The seminal work in this area was performed by Den Hartog, who extended his studies on SDOF systems with Coulomb friction to also include viscous damping in [19].Den Hartog's results include the exact solution for the steady-state time response to forced harmonic excitation, as well as closed-form expressions for the response amplitude and phase and for the boundary between continuous and stick-slip motion regimes.Shaw [30] extended Den Hartog's solution to address different static and kinetic friction forces, also accounting for positive and negative viscous damping, where the latter can be viewed as a crude model of destabilising aerodynamic forces.Furthermore, he carried out a stability analysis of the periodic response of these systems, showing that instability can arise for negative values of the viscous damping coefficient.Further works concerned with the stability of systems with mixed viscous and Coulomb damping can also be found in [31,32].Different SDOF models have also been explored in the literature: Hundal [33] extended Den Hartog's solution to address base-excited systems, while the response of support-excited systems was investigated by Levitan [34].Various approaches have been proposed to identify the damping parameters from the vibration of systems with mixed damping: Tomlinson and Hibbert [35,36] determined the amplitude of the friction force and the hysteretic loss factor by observing the Nyquist plots of the frequency response function, while Liang and Feeny [37] were able to identify the Coulomb and the viscous damping parameters of lightly-damped SDOF systems by comparing experimental results with the analytical results obtained by Den Hartog [19] and Hundal [33].During the last years, the fundamental research on SDOF oscillators has mostly focused on the exploration of asymmetric and chaotic solutions [38,39], the introduction of more complex friction models [39,40] or different types of damping [41] in the dynamic analysis and the investigation of the energy dissipation [42][43][44].
Differently from the SDOF case, little information is available in the literature regarding the dynamic behaviour of MDOF systems with mixed damping.Yeh [45] extended Den Hartog's approach [19] to deal with a harmonically base excited 2DOF system with two viscous dashpots and a Coulomb friction contact between the first mass and the wall.Sachs [46] proposed an analytical approach based on a recurrence scheme to determine the transient and the steady-state response of MDOF systems with mixed damping subjected to periodic loading.While these approaches can in theory be applied to any systems with viscous and Coulomb damping, both methods become unpractical for larger numbers of DOFs and do not lead to closed-form solutions, thus offering a limited insight on the mixed damping effects on dynamic behaviour of these systems.
In Ref. [24], the authors derived a closed-form solution for the continuous steady-state response of MDOF systems with a Coulomb friction contact subjected to harmonic loading by using a modal superposition procedure in those time intervals where the governing equations of the problem are linear.Although the introduction of viscous damping does not alter the piecewise nature of the problem, it can lead, in general, to non-diagonal modal damping matrices and, therefore, to coupling between the normal modes of the system.This is typically the case when viscous damping is generated by localised sources, such as the shock absorbers in a car suspension, or when it is used for modelling the dissipation in structures with highly non-homogeneous damping levels, e.g., the building models accounting for the soil-structure interaction [29].However, the coupling between the normal modes is usually negligible when viscous damping is used to model the damping resulting in structures by internal dissipation in the materials, junctions or interfaces between parts of the structure and non-structural elements [29,47].These latter forms of damping are sometimes referred to as "boundary damping" and, in built-up structures, they are commonly found to be at least an order of magnitude higher than the intrinsic material damping [48].Material and boundary damping are usually defined at a system level rather than in terms of individual element properties [29] and can be modelled using specific damping models, such as modal damping, where the damping of each vibrating mode is expressed by a modal damping ratio, assigned by measurement or based on experience [47].
In this paper, a more general solution is derived for the continuous steady-state response of MDOF systems with a friction contact, also accounting for the damping of the structure in the mechanical model.The dry friction force generated by the sliding between one of the masses of the system and an external wall is modelled according to Coulomb's law, including different static and kinetic friction coefficients, while the system's damping is introduced using the modal damping model.This analytical solution, along with the numerical results obtained with an ad-hoc approach for stick-slip response, is used to investigate the overall damping effect on the features and the motion regimes of their dynamic response.In particular, this contributions aims at establishing: (1) the Coulomb friction effect on the response amplitudes and the motion regimes displayed by damped vibrating systems in resonant conditions; (2) how the damping of the system can alter features of the dynamic response such as the low-and the high-frequency behaviours, the presence of invariant points across the displacement transmissibilities and the stick-slip behaviour with respect to those described in [24] for MDOF systems with Coulomb damping only.The presented results are particularly relevant for engineering systems such as friction dampers in buildings [8], car suspensions [9], energy harvesters [13] and taxing of airplanes models [14], which can be modelled as MDOF mass-spring systems with a friction contact during the early design stages.
The analytical solution for the continuous steady-state response of damped MDOF systems with a friction contact between one of the masses and a grounded wall is derived Section 2, along with closed-form expressions for the displacement transmissibility and phase angle of each mass.Furthermore, an analytical formulation is provided for the boundary between continuous and stick-slip regimes.In Section 3, an ad-hoc numerical approach is used for validating the analytical results and providing an insight on the behaviour of these systems in stick-slip regime.Finally, the resonant, low-and high-frequency behaviours, as well as the dynamic response in permanent sticking regime, are explored in Section 4.

General formulation and assumptions
Let us consider a damped MDOF mass-spring system with a friction contact subjected to harmonic excitation, as shown in Fig. 1.This system is characterised by the  ×  mass and stiffness matrices: The harmonic load  cos() is applied to the mass   of the system and a friction contact occurs between the   mass and a fixed wall, generating the friction force: where: The so-defined sgn() function can assume any value between − and , with  > 1, when the relative velocity in the contact is zero.
The actual value will be such that the static friction force, whose amplitude is  , is equal and opposite to the sum of the other non-inertial forces acting on the mass in contact.The relation between friction force and velocity expressed by Eqs. ( 1) and ( 2) is depicted in Fig. 2. The governing equation of the system shown in Fig. 1 can be written as: where  is the linear damping matrix of the system,  is the -dimensional displacement vector and   is the unit vector of the th coordinate.It is worth mentioning that Eq. ( 4) could also be used to describe a MDOF mass-spring system of more general topology, provided that the matrices  and  of the system are specified and that a single contact is acting on the linear motion of a mass in a specified direction.Introducing the non-dimensional time and displacement vector as: and dividing all terms by  , it is possible to rewrite Eq. ( 4) in a non-dimensional form, as: where: In particular, the non-dimensional mass matrix can also be written as: where: is the th mass ratio of the system and: will be denoted as frequency ratio.In Eq. ( 6), the friction ratio  has also been introduced as: Due to the presence of the sgn() function, Eq. ( 6) is a piecewise linear equation.In Ref. [24], it was observed that, in the absence of the linear damping term, a continuous steady-state solution of this equation can be determined by using a generalised version of Den Hartog's approach for SDOF systems with Coulomb friction [19].Let us assume that the continuous steady-state is symmetric, in particular that the second half of response period matches the negative of the first half.Under this assumption, the solution of Eq. ( 6) can be derived by considering the half-period included between any pairs of subsequent stationary points of the displacement of the mass in contact.In fact, in each of these intervals, the sign of the velocity is constant and the governing equations become linear.For instance, if the non-dimensional time interval [0, ] included between a maximum and a minimum of x is considered, Eq. ( 6) reduces to: In the above equation, it has been assumed that a phase lag is generated between the excitation and the response x by Coulomb damping.Since the response is not assumed to be monoharmonic, it is worth underlining that the phase angle   is referred to the maxima of these two functions.To solve Eq. ( 12), the following coordinate transformation is introduced: where the modal matrix: is obtained from the eigenvalue problem [29]: Rewriting Eq. ( 12) in modal coordinates yields: Nonetheless, while the modal mass and stiffness matrices: are diagonal, the same is not generally true for the modal damping matrix Ĉ = Ψ  Ψ.When Ĉ is fully populated, as it is usually the case when viscous damping is generated by localised physical elements [29], Eq. ( 16) represents a system of coupled equations, which cannot be solved separately.However, when the linear damping matrix is only associated with the damping of the system, it is common in the engineering practice to use damping models leading to a diagonal modal damping matrix [29].This can be obtained, for instance, by expressing the physical damping matrix as a linear combination of the mass and the stiffness matrices or by simply neglecting the off-diagonal terms of the modal damping matrix [49].In general, a physical damping matrix  leads to a diagonal modal damping matrix if it satisfies the condition  −1  =  −1 , derived in Refs.[50,51].It is worth noting that the governing equations can also be decoupled when this condition is not verified, by introducing complex valued mode shapes (see, e.g., [52]); however, this case is not dealt within the present contribution.In this paper, the damping of the system will be modelled according to the modal damping model, i.e., it will be assumed that the matrix  is such that: where  1 , … ,   represent the modal damping ratios of the  vibrating modes of the system.When the modal damping ratios are known, the non-dimensional physical damping matrix  can be obtained from Ĉ using the following relationship [29]: It is worth mentioning that the modal damping assumption can also be justified mathematically by using perturbation calculus, provided that the system is lightly damped and its natural frequencies are well-separate [52].

Modal superposition procedure
Since the modal damping matrix Ĉ has been diagonalised, Eq. ( 16) now represents a set of  uncoupled equations and the generic th equation of the system can be written as: The general solution of Eq. ( 20) can be obtained as: where the modal frequency ratio   is defined as: (22) and: is the complex response function of the th vibrating mode of the system in the absence of Coulomb damping.In Eq. ( 21), the integration constants   and   can be evaluated from the initial conditions for the th modal displacement and velocity within the half-period [0, ], which can be introduced as: where the initial values  0 and  ′ 0 are still unknown at this stage.Due to the symmetry of the steady-state response, the final conditions can be written at  =  as: These further conditions can be used to determine a relation between the phase angle   and the initial values  0 and  ′ 0 .This procedure is similar to that introduced in [24] for MDOF systems with Coulomb damping only and is therefore reported in detail in Appendix.The resulting expressions for cos   and sin   can be written as: (26) sin (27) where the functions: and: are hereby introduced as the first and the second damping functions of the th mode of the system.In particular, in the limit case of   = 0, the first damping function reduces to the damping function formulated for the th mode of a MDOF system with Coulomb damping only in [23,24], while   = 0.When  = 1, the formulations of both damping functions reduce to those derived by Den Hartog for SDOF systems with combined viscous and Coulomb damping [19].

Response amplitude and phase of the mass in contact
Multiplying by   the numerators and the denominators of Eqs. ( 26) and ( 27), considering their sums from 1 to  and introducing the initial conditions for the displacement and the velocity of the mass in contact: it is possible to express cos   and sin   as: (31) and: In the above equations,   denotes the non-dimensional amplitude of the response of the mass   .Following the definition of nondimensional displacement given in Eq. ( 5)(b), this quantity also coincides with the th displacement transmissibility of the system.Furthermore, the generic th complex response function, the first and the second damping functions of the MDOF system have been introduced as the modal superpositions of the expressions provided for the th mode in Eqs. ( 23), ( 28) and ( 29) respectively: L. Marino and A. Cicirello It is necessary to observe that, due to the presence of the function   , the expressions derived for cos   and sin   are complex.However, expressing the complex response function as   = |  | ∠  and performing some algebraic manipulations, it is possible to rewrite Eqs. ( 31) and ( 32) as: and: Using the relation cos 2 (  + ∠  ) + sin 2 (  + ∠  ) = 1, it can be obtained that the non-dimensional response amplitude of the mass in contact   is given by: It can be observed that this expression is formally identical to that derived by Den Hartog for a SDOF system with combined viscous and Coulomb damping [19] and reduces to that formulation when  = 1.The phase angle   can finally be obtained from Eqs. ( 36) and ( 37) or, more synthetically, from: where the function atan2 is defined according to the expression used in [53].

Steady-state time response of all masses
In Appendix, the following expression has been derived for the th modal displacement: The time response of the generic mass   of the system in the non-dimensional time interval [0, ] can be obtained from Eq. ( 40), by multiplying both sides by   and introducing the th equation from Eq. ( 13), as: The response of the mass in contact   is completely known at this stage, since the initial values of its displacement and velocity are equal to   and to zero respectively, as expressed in Eq. ( 30).However, these values are still unknown for all the other masses of the system.In order to determine them, let us multiply by   the numerators and denominators of Eqs. ( 26) and ( 27) and consider their sums from 1 to .The following relationships are obtained: Comparing Eqs. ( 42) and (43) with Eqs. ( 31) and ( 32) respectively, it can be obtained that: and: Introducing these expressions into Eq.( 41) and considering the real part only, the response of the generic mass   can be finally written as:

Response amplitude and phase of the masses not in contact
The amplitude and the phase angle of the response of the generic mass   cannot be determined in a closed form from Eq. ( 46) for the masses not in contact with the external wall, i.e. for  ≠ .In general, it is possible to calculate them numerically, determining the amplitude as the maximum absolute value of x and the phase angle as: where  ,max is the time instant where such a maximum is reached.Nonetheless, in Ref. [24] it was shown that, when   = 0, approximated closed-form expressions for these quantities can be obtained by considering the monoharmonic approximation of the response: It has been verified that Eq. ( 48) also offers a very good approximation of the response x when the damping of the system is taken into account.In fact, it can be observed that, in continuous non-sticking regime, modal damping further reduces the nonmonoharmonic effects introduced by Coulomb friction.While these effects can still be significant for the response of the mass in contact, they become mostly negligible for the masses not in contact.Substituting Eqs. ( 44) and ( 45) into Eq.( 48), it can be obtained that: where: ] −   (50) and: Therefore, the response amplitude can be obtained as: while the phase angle between the maxima of the displacements of   and   within the time interval [0, ] is given by: The phase angle of the response of the generic th mass can be finally calculated as   =   +   , where   is obtained from Eq. ( 39).It can be demonstrated that Eq. ( 52) reduces to Eq. ( 38) when  = .

Domain of validity of the solution
As specified in Section 2.1, the solutions derived in the current section only hold if the steady-state response of the system is continuously non-sticking.In order to obtain a closed-form expression for the domain of validity of these solutions, it is therefore required to impose that the sticking conditions are never simultaneously verified in the half period [0, ].It is worth recalling that a stop occurs in the response of a Coulomb friction oscillator when the relative velocity in the contact is zero and the amplitude of the overall dynamic loading (including spring forces and external excitation) acting in such a contact does not exceed the value of the static friction force [24].In the current investigation, it has been assumed that x′  < 0 in all the internal points of the half-period [0, ], while the velocity is zero by definition at the ends of this time interval.Therefore, it will be necessary to explicitly impose that the second condition is not verified in these points.However, due to the symmetry of the steady-state response, it is sufficient to write this condition for  = 0. Thus, it can affirmed that stops will not occur if the following conditions are met: where   is the Kronecker delta.In order to express the domain of validity in terms of the friction ratio, let us substitute Eq. ( 46) into Eq.(54)a.The following inequality is obtained: where: and: From Eqs. ( 38) and ( 55), it can be obtained that the condition expressed in Eq. ( 54)a is verified when: In order to derive a similar expression from Eq. (54)b, it is possible to rewrite it in modal terms by introducing the coordinate transformation in Eq. ( 13) and considering the generic th modal coordinate: It is useful to express the term cos   as a function of  0 and  ′ 0 .From Eqs. ( 26) and ( 27), it is obtained that: Furthermore, it is easily obtained from the definition of   , provided in Eq. ( 23), that: Thus, it follows that: Introducing Eq. (62), it is possible to rewrite Eq. (59) as: Let us now consider the sums of both sides of Eq. ( 63), multiplied by   .Introducing the function: and considering that, from Eqs. ( 8) and ( 17) it is obtained that: it is possible to rewrite the above inequality as: The condition expressed in Eq. ( 66) is verified if: or: In order to obtain a continuous non-sticking response, the inequalities expressed in Eq. ( 58) and one between Eq. ( 67) and (68) must be simultaneously verified.However, it has been observed that Eqs. ( 58) and (68) never occur concurrently.Therefore, the domain of validity of the mathematical solution presented in this section can be expressed as: where the RHS represents the value  lim of the friction ratio at the boundary between continuous and stick-slip regimes.It is worth noting that, although the analytical solutions derived in this section for the continuous steady-state response are independent of the ratio  between the static and the kinetic friction forces, this parameter now appears in the boundary expressed in Eq. ( 69).This implies that, despite having no effect on the dynamic response in continuously sliding regime, larger values of the static friction force can affect (and reduce) the region of the parameter space where continuous motions are observed.Due to the presence of the function   , the evaluation of the boundary described by Eq. ( 69) requires the numerical calculation of the maximum of the time-dependant expression in Eq. ( 57) for each sets of parameters.In order to reduce the computational cost associated to this procedure, the following approximation can be considered.It has been observed numerically that, in most cases, the maximum value of the function   occurs at  = 0.This effect is due to the presence of the exponential term  −      in the expression of   and, therefore, exceptions are only observed when   is nearly zero.Thus, a first-order expansion of the function   around  = 0 can be considered: Neglecting the second-order terms in the above equation, it is obtained that: and, by substituting this expression into Eq.( 56), it is possible to write an approximated expression of the function   as: Introducing this approximated formulation into Eq.(69), it can be observed that the RHS of Eq. (72) coincides with   for  = 1, and becomes larger than   if  > 1.Therefore, Eq. (67) can directly be considered as an approximation of the boundary between continuous and stick-slip regimes.In particular, it has been observed that discrepancies between the exact and the approximated expressions can arise at low frequency ratios only when the modal damping ratios are nearly zero, consistently with the above observations.The case   = 0, discussed in [24], is therefore the case where the most significant disagreement is observed between the exact and the approximated boundary.

Numerical investigation of the steady-state response
This section introduces a numerical approach for the evaluation of the time response, displacement transmissibility and phase angle of damped MDOF systems with a Coulomb friction contact, with the scope of providing a numerical validation of the solutions presented in the previous section and achieving a complete overview of the dynamic behaviour of these systems, including stick-slip responses.

Description of the numerical approach
The main challenge in the numerical time integration of the set of governing equations reported in Eq. ( 6) is doubtlessly related to the occurrence of stick-slip in the response.In fact, while continuous responses can be calculated with standard numerical solvers, stick-slip motions are characterised by the presence of rapid variations due to the transitions between the sticking and the sliding phases and lead to numerical stiffness in the problem [3].The approach proposed in this section consists in using standard integration methods by setting explicit conditions to account for the transition between different regimes.A schematic representation of the implemented numerical algorithm is given in Fig. 3 and the major steps are detailed in what follows.• The minimum number of input parameters required to determine the response of the MDOF systems investigated in this paper is equal to 3 + 1, i.e. the frequency ratio , the friction ratio , the ratio , the  modal damping ratios, the  − 1 mass ratios ( 1 = 1 by definition) and the  − 1 stiffness ratios: where  1 = 1 by definition.The non-dimensional matrices ,  and  can be evaluated from the input parameters by implementing Eqs. ( 8), ( 19) and (7) or, alternatively, they can directly be specified by the user.The latter option is particularly practical when the dimensional mass, damping and stiffness matrices of the system are known; in this case, the mass, stiffness and frequency ratios will not be needed.Finally, the user can specify the number of DOFs of system, the masses in contact and under harmonic excitation respectively, the initial conditions for the displacement and the velocity and the number of excitation periods   .Since the length of each period is equal to 2 in the current formulation, the final time will be equal to 2  .• As a first step, the algorithm checks if the sticking conditions: are verified.In the above equations, the function: denotes the amplitude of the sum of the overall dynamic loading acting in the contact.A sticking or a sliding phase will take place depending on whether these conditions are met or not.• During the sliding phases, the motion of all masses is continuous and the solution is nonstiff.Therefore, the governing equations are integrated by using a variable-step Runge-Kutta (4, 5) method, implemented in the Matlab function ode45 [53].These phases are terminated by the event condition x′  = 0; when this happens, the algorithm checks if the second sticking condition from Eq. ( 74) is also verified to determine if the following phase will be sticking or sliding.
• During the sticking phases, the displacement of the mass in contact will remain constant and its velocity will be equal to zero.
However, all the other masses will keep oscillating and, therefore, numerical integration will also be needed at this stage.
Fig. 4. Steady-state time response of a 5DOF system with equal masses and springs where a friction contact and a harmonic load are applied to  1 for  = 1,  = 0.2,  = 0.9 and modal damping ratios 0.01 (a) and 0.1 (b): comparison between analytical (continuous lines) and numerical (round markers).
Fig. 5. Steady-state time response of a 5DOF system with equal masses and springs where a friction contact occurs on  3 and a harmonic load is applied to  1 for  = 1,  = 0.2,  = 0.9 and modal damping ratios 0.01 (a) and 0.1 (b): comparison between analytical (continuous lines) and numerical (round markers).
Although the motion of   during the stop is already known, it is convenient to integrate all the  equations of the system, disregarding the friction force, and imposing that x = x0 at each step, as shown in Fig. 3.In fact, this approach allows the use of the same mass and stiffness matrices used during the sliding stages.The stop will end when the dynamic force acting in the contact exceeds the static friction force, i.e. when  > ; this is again imposed as an event condition in the function ode45.
• The simulation will end when the final time 2  is reached.The response amplitudes can be determined as the maximum absolute values assumed by each mass displacement within the time interval [2(  − 1), 2  ], corresponding to the last excitation period.Similarly, given that the lower bound of this interval corresponds to a maximum of the excitation, the phase angle of each mass motion can be determined as the non-dimensional time instant corresponding to the maximum mass displacement.However, it is worth nothing that, when stick-slip occurs, the maximum displacement of the mass in contact usually coincides with a sticking phase rather than a single point.While this do not affect the evaluation of the response amplitude, in this case it is not possible to evaluate the phase angle between the maxima of the excitation and of x .Therefore, the numerical phase angles will only be evaluated for continuous responses.

Numerical validation of the analytical results
This subsection presents a numerical validation of the analytical solutions derived in Section 2. In particular, the validation of the continuous steady-state response is addressed in Section 3.2.1,considering as example a 5DOF system with equal masses and springs for two different locations of the friction contact.In Section 3.2.2, the closed-form expressions obtained for the response amplitude and phase are also validated for varying frequency and friction ratios, considering for simplicity the cases of a SDOF system and of Fig. 6.Displacement transmissibility and phase angle of a damped SDOF system with a Coulomb friction contact under harmonic excitation for  = 1 and varying frequency, friction and modal damping ratios: analytical vs numerical.Fig. 7. Motion regimes of a damped SDOF system with a Coulomb friction contact under harmonic excitation in the parameter space  −  for  = 1 and three different modal damping ratios.a 2DOF system with excitation and friction contact applied to different masses.For both systems, the analytical boundaries of the motion regimes are also shown in the parameter space  − .

Steady-state time response
The steady-state mass motions described by Eq. ( 46) have been compared to the numerical results obtained by using the algorithm introduced in Section 3.1 for different parameters sets.In particular, the numerical responses have been evaluated setting a number of cycles equal to   = 200 and zero initial conditions for the relative displacement and velocity.It has been observed that most responses converge to a steady-state solution for a significantly lower number of cycles.During the integration process, the absolute tolerance was set to 10 −6 .An excellent agreement has been observed for all the cases investigated.Figs. 4 and 5 show the comparison between analytical and numerical results for a 5DOF system with equal masses and spring for  = 0.9,  = 0.2,  = 1 and for two different modal damping ratios.In particular, in Fig. 4 the harmonic excitation and the friction contact are both applied to the mass  1 , while in Fig. 5 the contact occurs between  3 and a fixed wall.For simplicity, in both cases equal modal damping ratios have been considered for the all vibrating modes, but it is worth remembering that different values can be taken into account within the proposed analytical and numerical formulations.
In Ref. [24], it was observed that the dynamic behaviour of MDOF systems with a Coulomb friction contact is mostly affected by the location of the harmonic and of the friction forces and, more specifically, different behaviours are observed depending on whether these forces are applied or not to the same mass.The phenomenon can be briefly explained as follows: • when  = , all the masses oscillate in phase or in phase-opposition and, therefore, behave similarly to the mass in contact; • when, for instance,  > , the masses   , … ,   will still oscillate in phase or in phase-opposition, while the masses  1 , … ,  −1 will exhibit different phase angles and, in general, a richer dynamic behaviour.
In Figs.4a and 5a, it can be observed that this property is substantially maintained when   = 0.01.Differently, in Figs.4b and 5b, where all modal damping ratios have been set to 0.1, each mass oscillate with a different phase angle.Hence, it can be concluded that modal damping alters the phase angle among the mass motions, but this effect only becomes relevant for significant damping levels in the system.

Response amplitude, phase and motion regimes for varying frequency ratio
A numerical validation has also been carried out for the transmissibility and phase angle curves by using the algorithm presented in Section 3.1.As mentioned, different dynamic behaviours are observed in MDOF systems depending on whether the friction contact occurs on the excited mass or a different one.Since the MDOF behaviour has already been extensively investigated in [24], the purpose of this study is to explore how these curves are affected by Coulomb friction when the damping of the system is also taken into account.To this end, the simplest case-studies reflecting these two different configurations are considered in what follows: (i) a SDOF system and (ii) a 2DOF system with a contact on  2 and a harmonic loading on  1 .
Fig. 6 presents the transmissibility and phase angle curves for the SDOF case, obtained for the values [0.001, 0.01, 0.1] of the damping ratio  , while Fig. 7 shows the variation of the analytical boundaries of the motion regimes in the parameter space  − .In the latter figure, the boundary between continuous and stick-slip response has been evaluated from Eq. (69), while the boundary between sliding and permanent sticking regimes has been obtained using the approach introduced in Section 4.3.The response amplitude and phase of the 2DOF systems are plotted in Figs. 8 and 9 for the masses  1 and  2 respectively, while Fig. 10 shows the evolution of the boundaries of the motion regimes for increasing modal damping ratios.In this case, the values Fig. 8. Displacement transmissibility and phase angle of the mass  2 for a 2DOF system with equal masses and springs, a Coulomb friction contact on  2 and a harmonic excitation on  1 for  = 1 and varying frequency, friction and modal damping ratios: analytical vs numerical.  = [0.001∕, , 0.01∕ , , 0.1∕ , ], where  , is the th natural frequency ratio, have been considered in order to observe a similar damping effect for both vibrating modes in each case.
A very good agreement can be observed from the comparisons between the analytical and the numerical transmissibilities and phase angles presented in Figs. 6, 8 and 9.In the same figures, the numerical transmissibilities have also been represented for stickslip responses to provide a more complete insight of the dynamic behaviour of these systems across the different motion regimes.It has also been verified that continuous, stick-slip and permanent sticking regimes occur in the numerical response as predicted by Figs. 7 and 10.The results presented in Figs.6-10 will be further discussed in the following section, where the response features of these systems are explored.

Features of the dynamic response
In this section, the main features characterising the dynamic response of damped MDOF systems with Coulomb friction will be discussed.These features include the resonant, the low-and the high-frequency behaviours, the presence of invariant points and inversions across the transmissibility curves and the behaviour of these systems in permanent sticking regime.Although the results presented in the previous section and illustrated in Figs.6-10 will be referred to as example, this discussion will concern the more general case of the MDOF system of  masses, a contact on the mass   and a harmonic forcing applied to the mass   depicted in Fig. 1.In particular, this section aims at describing how the response features of a MDOF system with Coulomb damping only, dealt within [24], evolve when the damping of the structure is also taken into account.

Resonant behaviour
In the absence of other forms of damping, Coulomb friction cannot provide finite resonant peaks in discrete mechanical systems unless stick-slip or permanent sticking occur at resonance [24].However, the resonant behaviour of these systems is deeply affected by modal damping.In fact, in Figs. 6, 7 and 9, it can be observed that damped systems with Coulomb friction also present finite resonances in continuous motion regime.A procedure for determining an approximated expression for the amplitude of these finite resonant peaks is presented in what follows.
According to Craig [29], the values of the modal damping ratios typically lie in the range 0 ≤   ≤ 0.1; therefore, it can generally be assumed that   ≪ 1.Based on this consideration, it can also be assumed, at this stage, that resonant peaks occur in correspondence of the undamped natural frequencies of the system.Under these assumptions, an approximated expression can be derived for the amplitude of the th resonant peak of the response of the generic mass   by evaluating Eq. ( 52) for   → 1 and   ≪ 1.As a first step, let us determine how the complex response function   and the damping functions   and   are affected by these assumptions.
• Regarding the complex response function, in proximity of the th resonant peak, the th term of the summation in Eq. (33) becomes significantly larger than the other  − 1 terms.Therefore, it can be written that: where the expression of   is provided in Eq. ( 23).When the th resonance occurs, the amplitude of the complex response function is therefore given by: while the phase of this function will be equal to −∕2 or ∕2 depending on the sign of the product     .Thus, it can be demonstrated that, in Eq. ( 52): • Let us now evaluate the second damping function of the th vibrating mode for   → 1.If   ≪ 1, it is obtained from Eq. ( 29) that: It can be demonstrated that the above expression tends to infinity when   → 0. Therefore, since   ≪ 1, the th term of the summation in Eq. ( 35) will also grow much larger than the other terms when   → 1.Thus, in proximity of the th resonance, the second damping function will assume the value: • Finally, the first damping function   does not present any peculiar behaviours at resonance.In fact, although it can be shown from Eq. ( 28) that   = 0 for   → 1, the other  − 1 terms of the summation from Eq. (35) 52), it is obtained that, for  =  , : Since at resonance, as previously stated, the functions |  | and   become much larger than the damping functions   , the above equation can be approximated as: Finally, substituting Eqs. ( 77) and (80) into Eq.( 82), it is possible to write: Eq. ( 83) can be used for estimating the response amplitude of each mass of the system in correspondence of each resonant peak of the system.The most important implication of this expression is that, when   ≠ 0, Coulomb friction also reduces the resonant amplitudes without introducing stick-slip in the response.Furthermore, it can be observed that, for given modal damping ratios, the amplitude of these peaks decreases linearly with the friction ratio.To the best of the author's knowledge, this latter phenomenon has never been discussed in previous publications.However, it is worthwhile mentioning that Den Hartog investigated experimentally the response amplitude at resonance for a SDOF system with Coulomb friction and a small amount of viscous damping; in his results, the resonant amplitude shows an approximatively linear reduction with the friction ratio [19].
The dependence of the resonant peak amplitude on  and  is shown in Fig. 11 for the SDOF case; similar patterns have also been observed for the resonances of more complex systems.For a SDOF system, Eq. (83) reduces to: In Fig. 11, a boundary has also been represented to delimit the domain of validity of this formula, which is only valid in continuous motion regime.Based on the above observations regarding the response and damping functions, an approximated expression of this boundary can be determined by evaluating Eq. (69) for   → 1.After some algebraic manipulations, it is obtained that the response is continuous at the th resonance when the friction ratio is smaller than: consistently with the value  , = |  ∕  | (corresponding to ∕4 in the SDOF case) derived in [24] for   = 0. Fig. 11 shows how the value of  , remains approximatively constant for  ≤ 0.01 and only decreases significantly for larger values of the modal damping ratio.Therefore, while it is important to take the damping of the system into account when evaluating the amplitude of the resonant peaks, the aforementioned boundary expression from [24] can be considered a good estimation of the minimum value of the friction ratio for which stick-slip occurs at resonance in lightly-damped systems.In Fig. 12, the estimation of the resonant amplitude provided by Eq. ( 83) is compared to the amplitude of the response at  = 1 and to the actual resonant amplitude evaluated from Eq. (52), showing an excellent agreement with both quantities.It can also be observed that, in proximity of the boundary between continuous and stick-slip regimes, the actual value of the resonant amplitude becomes slightly larger than the estimated one.In Figs. 6, 8 and 9, it can be seen that, for large amounts of the friction ratio, the resonant peak is slightly shifted to the left and, therefore, it can be underestimated by Eq. (83).Nonetheless, the actual and the estimated values of the resonant peak have been compared for several MDOF systems and for different modal damping ratios; their agreement was very good in all the cases investigated.The linear dependency of   on the friction ratio is also shown in Fig. 12b for different values of the damping ratio.It can be observed that the slope of the curves decreases with  ; however, the rate of variation of this slope also becomes smaller as the damping ratio is increased.

Low-and high-frequency behaviours
In Figs.6-10, it is possible to observe that the dynamic behaviour of SDOF and MDOF systems at low and high frequency ratios is not significantly affected by the system's damping.Regarding the low-frequency behaviour, it can be easily verified that the boundary expressed in Eq. ( 69) tends to zero for  → 0; therefore, the quasi-static behaviour of damped systems with Coulomb friction will also be characterised by the occurrence of stick-slip in the response.However, as shown in Fig. 13, modal damping can have a smoothing effect on the quasi-static response of these systems; this is particularly evident from the curve corresponding to  = 0.1.The patterns displayed in the transmissibility curves in Figs. 6, 8 and 9 at low frequencies do not generally show significant differences for varying damping ratios.This implies that, in most cases, modal damping can reduce the number of stops but does not change significantly the amplitude of the stick-slip responses, in agreement with Fig. 13.
The behaviour of these systems at high frequency ratios is also similar to that observed in systems with Coulomb damping only in [23,24].In the above figures, it can be clearly observed that the transmissibilities always tend to zero when  → ∞, as expected.Therefore, also in this case, the discussion focuses on determining if stick-slip or permanent sticking can occur in these systems at high frequencies.By evaluating the limit of Eq. (69) for   → ∞, it can be verified that the boundary between continuous and stick-slip regimes tend to the same value determined in [24]: This means that only systems with the harmonic and the friction forces applied to the same mass can display a continuous response at high frequencies.This can also be observed from Fig. 7, where all the boundaries between continuous and stick-slip regimes displayed for varying  tend to the asymptotic value  ∞ = 0.537.Conversely, in Fig. 10 these boundaries tend to zero and, therefore, stick-slip is always expected for high frequency ratios.

Invariant points and stuck configurations
In Fig. 8, it is possible to observe how the inversion process of the transmissibility curves and the onset of the resonant peaks of the stuck configuration are not heavily affected by the system's damping.From a mathematical point of view, it can be observed that invariant points are not admitted by Eq. ( 52) unless   = 0.However, it can be observed that the inversions of the transmissibility curves still occur in a small region located in proximity of the invariant points of MDOF systems with Coulomb damping only, which can be determined from [24]: where   and   are the response and the first damping functions evaluated for   = 0. Therefore, the points determined from the above equation can also be considered, with good approximation, for damped systems.Permanent sticking between the mass and the wall in contact occurs in mechanical systems with a Coulomb friction contact when the maximum amplitude of the overall dynamic loading acting in the contact is not sufficient to overcome the static friction force.In MDOF systems, even when a contact is permanently stuck, some masses of the systems can still be excited by the dynamic loading; the subsystem which keeps oscillating in permanent sticking regime is usually defined as stuck configuration or mode [23,24].In the 2DOF case discussed in Section 3.2, the mass  1 will also respond to the harmonic excitation when  2 is stuck.The resonant peak associated to the stuck configuration of this system can be observed in Fig. 8, where it is represented by a black dotted line.Although the procedure for evaluating the response of system in stuck conditions and the boundary between the sliding and the permanent sticking regimes is similar to that described in [23], it must be considered that modal damping will also be present in the stuck configurations.Let us consider, for instance, the case of a MDOF system with  > .When permanent sticking occurs, the masses  1 , … ,  −1 will keep oscillating.Therefore, the response of the system can be determined from the linear equation: with a standard modal superposition procedure.In the above equation, the stuck mass, damping and stiffness matrices are determined considering the first  −1 rows and columns of the matrices ,  and .If needed, the modal damping ratios  * 1 , … ,  * −1 of the stuck configuration can be determined by using Eq. ( 18), after that the stuck natural frequencies and mode-shapes have been determined from Eq. ( 15), referring to the matrices  * and  * .Finally, the boundary  * lim can be evaluated from [24]: The same procedure can be used for the case  < , taking into account the masses  +1 , … ,   , while the system will be fully stuck if  = .From Figs. 8 and 10, it can be observed that modal damping only affects  *  and  * lim at resonance, while the starting values at  = 0 and the asymptotic behaviour at high frequencies remain unchanged.

Concluding remarks
In this paper, an analytical solution has been derived for the continuous response of damped MDOF systems with a friction contact subject to harmonic excitation, modelling the friction force according to Coulomb's law (with different static and kinetic friction coefficients) and the damping of the system according to the modal damping model.In particular, closed-form expressions have been obtained for the steady-state time response, the displacement transmissibility and the phase angle of all the masses of the system.Furthermore, exact and approximated formulations have been provided for the boundary between continuous and stick-slip regimes.
The results obtained from these analytical solutions have been validated with a numerical approach, taking into account: (i) a SDOF system; (ii) a 2DOF system excited on  1 and with a contact on  2 ; (iii) a 5DOF system with excitation and contact on  1 ; (iv) a 5DOF system excited on  1 and with a contact on  3 .For the systems (i) and (ii), a very good agreement has been observed for the transmissibilities and the phase angles, and the cases (iii) and (iv) have shown that an excellent agreement can be obtained between the analytical and the numerical time responses even when several DOFs are involved.
The features and the motion regimes exhibited by the dynamic response of these systems have been investigated, including the resonant, low-and high-frequency behaviours, the presence of invariant points and inversions of the transmissibility curves and the stuck configurations determined by the occurrence of permanent sticking.This investigation has revealed that, while most of the response features are not significantly altered by modal damping, the resonant behaviour is very different from that observed in SDOF and MDOF systems with Coulomb damping only.In fact, when the system's damping is accounted for in the mechanical models, the resonant peaks are also finite in continuous motion regime.An approximated expression was derived for the amplitude of the resonant peaks of a MDOF system with mixed damping.Besides enabling a quick estimation of the resonant amplitudes, this formula also highlights that, for given modal damping ratios, such amplitudes decrease linearly with the friction ratio and, therefore, with the intensity of the friction force.This implies that, in systems with non-negligible damping levels, dry friction can avoid large resonant peaks without introducing stick-slip or permanent sticking in the response.
Overall, the present investigation has shown that mechanical models with Coulomb damping only can provide an acceptable description of lightly-damped structures where a friction contact is the predominant source of damping.Even in the presence of higher levels of damping, certain features such as the high-frequency behaviour and inversions across the transmissibility curves can still be studied referring to the formulations derived in [24].However, in the latter case, the damping of the system should also be taken into account in the investigation of the resonant behaviour.In the above equations, the functions: which has also been written in Eq. ( 40) and represents the solution of the modal problem dealt within this appendix.It is worth observing that the initial values  0 and  ′ 0 are not known at this stage.Nonetheless, they can be obtained by imposing further conditions on the initial displacements and velocities of the masses, after considering the transformation from modal to physical coordinates written in Eq. ( 13).This is done in Sections 2.3 and 2.4.

Fig. 1 .
Fig.1.-dimensional MDOF system with a Coulomb friction contact between the mass   and fixed wall subjected a harmonic load acting on the mass   .

Fig. 2 .
Fig. 2. Friction force versus relative velocity in the contact according to Coulomb's law.

Fig. 3 .
Fig. 3. Flowchart of the numerical algorithm implemented for the calculation of the response of a damped MDOF system with a Coulomb friction contact under harmonic excitation.

Fig. 9 .
Fig.9.Displacement transmissibility and phase angle of the mass  2 for a 2DOF system with equal masses and springs, a Coulomb friction contact on  2 and a harmonic excitation on  1 for  = 1 and varying frequency, friction and modal damping ratios: analytical vs numerical.

Fig. 10 .
Fig.10.Motion regimes of a 2DOF system with equal masses and springs, a Coulomb friction contact on  2 and a harmonic excitation on  1 in the parameter space  −  for  = 1 and three different modal damping ratios.

Fig. 11 .
Fig. 11.Non-dimensional amplitude of the resonant peak of a SDOF system in continuous non-sticking regime for varying friction and modal damping ratio.

Fig. 12 .
Fig.12.Non-dimensional amplitude of the resonant peak of a SDOF system in continuous non-sticking regime for varying friction ratio: comparison between the estimated and the actual amplitudes (a) and estimated value for varying modal damping ratio (b).

Luca Marino :
Conceptualization, Methodology, Formal analysis, Investigation, Software, Writing -original draft, Writing -review & editing.Alice Cicirello: Conceptualization, Writing -review & editing, Supervision, Funding acquisition.Declaration of competing interest Luca Marino reports financial support was provided by Engineering and Physical Sciences Research Council.Alice Cicirello reports financial support was provided by Engineering and Physical Sciences Research Council.