Friction-induced vibrations in the framework of dynamic substructuring

In complex vibrating systems, contact and friction forces can produce a dynamic response of the system (friction-induced vibrations). They can arise when different parts of the system move one with respect to the other generating friction force at the contact interface. Component mode synthesis and more in general substructuring techniques represent a useful and widespread tool to investigate the dynamic behavior of complex systems, but classical techniques require that the component subsystems and the coupling conditions (compatibility of displacements and equilibrium of forces) are time invariant. In this paper, a substructuring method is proposed that, besides accounting for the macroscopic sliding between substructures, is able to consider also the local vibrations of the contact points and the geometric nonlinearity due to the elastic deformation, by updating the coupling conditions accordingly. This allows to obtain a more reliable model of the contact interaction and to analyze friction-induced vibrations. Therefore, the models of the component substructures are time invariant, while the coupling conditions become time dependent and a priori unknown. The method is applied to the study of a finite element model of two bodies in frictional contact, and the analysis is aimed to the validation of the proposed method for the study of dynamic instabilities due to mode coupling.


Introduction
In complex mechanical systems, the relative motion among components generates contact and friction forces at the contact interfaces. These forces produce a dynamic response of the system, known as frictioninduced vibrations [1,2], i.e., vibration and noise that in certain conditions can become relevant in terms of both structural integrity and comfort [3].
Friction-induced vibrations, and more generally contact problems, are naturally prone to be tackled using dynamic substructuring, because at least two different bodies in contact need to be considered. In fact, the aim of dynamic substructuring is to predict the dynamic behavior of a system made by component subsystems, starting from their dynamic behavior which is assumed to be known. A review on dynamic substructuring is provided in [4], where a general framework to set up substructuring problems is proposed. The main advantages of the substructuring approach are: to allow the use of a different kind of model for each substructure (finite element model or experimental model) in order to predict the dynamic behavior of the coupled structure [5]; to allow the use of reduced order models by condensing the system matrices with the only requirement of retaining the interface degrees of freedom [6,7]; to assess the effects of subsystem modifications on the dynamics of the entire model; and to identify the dynamic behavior of substructures from coupled system data [8][9][10][11]. With regard to friction-induced vibrations, a substructuring approach is applied in [12] to a time-invariant brake system showing its advantages in terms of computational burden without losing accuracy.
Although originally conceived for time-invariant systems, dynamic substructuring has been adapted with few modifications to a significant subset of time-variant systems built from time-invariant substructures subjected to time-variant coupling conditions [13][14][15] or to predict position-dependent dynamics of systems that change their configuration over time [16,17]. The numerical analysis of configuration-dependent problems in the framework of dynamic substructuring provides interesting results with relatively low computational effort.
However, in the approach presented in [15], the timevariant compatibility and equilibrium conditions arising from sliding contact are assumed to be known a priori. This is a limitation because it does not allow to account for the relative displacement at the contact interface due to the system deformation. In dynamic contact problems [18], the relative displacement due to the elastic deformation is the main cause of frictioninduced vibrations phenomena such as dynamic instabilities [19][20][21], stick-slip [22,23] or sprag-slip [24,25].
In this paper as in [15], sticking and detachment are not considered. However, the time-variant coupling conditions due to sliding are not anymore assumed to be known a priori, but they are estimated accounting for the deformation of the contacting bodies: This allows to account for the effects of such deformation on the relative displacement and, moreover, it introduces geometric nonlinearities. The sliding contact can be without or with friction. With friction, the set of degrees of freedom to which equilibrium conditions apply includes the tangential directions at the contact interface, which are not considered in the compatibility conditions because of sliding. This is in fact a non collocated interface as defined in decoupling problems [10,26].
The proposed approach is applied to a beamon-beam system [27] described by a finite element model and by considering a single contact point. This allows to evaluate the effectiveness of the approach in both stable and unstable conditions and to envisage possible applications to more complex systems.
The problem is tackled in the time domain using dual assembly: However, in this case a singularity problem arises that can be overcome by adapting the forward increment Lagrange multipliers method [28] to obtain a nonlinear forward increment dual assembly. The harmonic content of transient responses is compared with the results of time-dependent complex eigenvalue analysis and with the results of time-dependent frequency response function. Time-dependent complex eigenvalue analysis is performed using primal assembly. Time-dependent frequency response function (FRF) is defined by assuming that the rate of variation of system configuration is low with respect to the characteristic frequencies.
2 Contact problems using dynamic substructuring A coupled structural system composed of n connected subsystems is considered. The equation of motion of the r th linear time-invariant subsystem can be written as follows: where M (r ) , C (r ) and K (r ) are the mass, damping and stiffness matrices of subsystem r ; u(t) (r ) is the vector of displacements of subsystem r ; f (t) (r ) is the vector of external forces on subsystem r ; g(t) (r ) is the vector of connecting forces with other subsystems (internal constraint forces). The equation of motion of the n independent subsystems before coupling can be expressed in a block diagonal format as: with: (1) . . . (1) . . . (1) . . . (1) . . . (1) . . . (1) . . .
Compatibility requires that a given pair of matching DoFs at time t, e.g., DoF l on subsystem r and DoF m on subsystem s, must share the same displacement, that is u (r ) This condition can be generally written as: where each row of B C (t) defines compatibility between a pair of matching DoFs at time t. B C (t) can be written by splitting the contribution of the different subsystems: Equilibrium of internal constraint forces, associated with the compatibility conditions and possibly with sliding friction, requires that the sum of connecting forces at a pair of matching DoFs at time t is zero, e.g., g   k (t) = 0. When friction occurs, the equilibrium condition applies also to the tangential direction at the contact interface, which is not included among the compatibility DoFs whenever sliding occurs. This defines a so-called non-collocated interface.
The resulting set of equilibrium conditions can be expressed as: Note that the number of rows of L E (t) T is given by the sum of the number of non-interface DoFs and of the number of pairs of equilibrium interface DoFs.
Equations (2)(3)(4)(5)(6) can be gathered to obtain the socalled three-field formulation, i.e., the overall system of equations representing the coupled system: Note that the three-field formulation is very similar to that describing a coupled system with a timeinvariant interface, the difference being represented by B C and L E which now depend on time t. Since time-variant coupling conditions are involved, the timedomain approach is the most appropriate. Substructuring in time domain has been tackled by Rixen and van der Valk [29,30], by using impulse response functions, but the coupling conditions (compatibility and equilibrium) have been assumed as time invariant. If the solution is approached in the time domain, both primal assembly and dual assembly could be used [4,31].

Primal assembly in time domain
In the primal assembly, a unique set of interface DoFs is selected and the interface forces are automatically canceled by enforcing the equilibrium condition. The unique set of DoFs q includes also non-interface DoFs and satisfies: where L C (t) differs from the matrix L E (t) introduced previously when friction forces exist at the interface. Since Eq. (8) states that the DoFs of all subsystems are obtained from the unique set q, compatibility holds for any set q, i.e., Therefore, Since Eq. (9) is satisfied by the choice of the unique set q, the system of equations (7) becomes: By pre-multiplying the dynamic equilibrium equation by L E (t) T and by noting that L E (t) T g(t) = 0, Eq. (11) reduces to: i.e., wherẽ Equation (13) can be recast in state space form: i.e., Being A(t) a time-dependent matrix, Eq. (16) cannot be solved analytically but requires numerical time step integration that involves the computation and inversion of matrixM at each time step. The computational burden prevents from using primal assembly for time integration. Primal assembly, however, can be used to formulate a time-dependent complex eigenvalue problem starting from Eq. (13) with L E (t) and L C (t) computed without accounting for the geometric nonlinearity due to the elastic deformation of the system.

Dual assembly in time domain
In dual assembly, each interface DoF is considered as many times as there are substructures connected through that DoF. At time t, the equilibrium condition g (r ) l (t) + g (s) m (t) = 0 at a pair of interface DoFs can be ensured by setting g (r ) l (t) = −λ and g (s) m (t) = λ. Therefore, equilibrium at all the interface DoFs can be ensured by writing the connecting forces in the form: where λ are Lagrange multipliers corresponding to connecting force intensities and B E (t) is different from the matrix B C (t) previously defined to enforce the compatibility condition, because B E (t) must also account for the equilibrium of friction forces at the interface. The interface equilibrium condition (6) is thus written: Therefore, Since Eq. (18) is always satisfied by any set of connecting force intensities λ, the system of equations (7) becomes: (20) and in block matrix format: Equation (21) cannot be explicitly solved in the time domain due to its singularity. An approach similar to the forward increment Lagrange multipliers method [28] can be used to overcome the singularity problem.

Nonlinear forward increment dual assembly
This formulation relates constraints in the deformed position at time t n+1 = t n + h, where h is the integration time step, with Lagrange multipliers at time t n . In this case, the geometric nonlinearity due to the elastic deformation of the system is accounted for. The equation of motion becomes: In this notation, the subscripts n and n + 1 indicate the functions at time t n and t n+1 . The first equation of (22) can be rewritten as: For explicit integration of the equation, the Newmark β 2 scheme can be used, that for β 2 = 0.5 reduces to the central difference scheme: To ensure the stability of the Newmark central difference scheme, the time step h must satisfy the following inequality: where l is the minimum distance between two adjacent nodes, ρ is the mass density, E is the Young's modulus, and ζ Max is the maximum value of the modal damping factor.
By substituting the central difference expression of u n in the equation of motion (23), and by usingu n = (u n − u n−1 ) / h, one obtains: The displacement u at time n + 1 can be rewritten as the sum of a predictor u n+1 and a corrector u c n+1 , i.e., and The predictor u n+1 can be directly evaluated since it depends on the displacements of the system at times t n and t n−1 , but it does not account for contact forces λ. The corrector u c n+1 depends on the matrix B E evaluated at time t n+1 and represents an incremental displacement associated with the contact forces λ. It is unknown because contact forces λ n must be estimated and constraints B En+1 must be defined.
Contact forces can be estimated, by first substituting Eq. (27) in the second line of Eq. (22) that expresses the compatibility condition: then, by substituting Eq. (29) in Eq. (30): from which, λ n can be obtained as and it depends on matrices B E and B C defined at time t n+1 . In a previous work, the authors proposed an estimation of matrices B E and B C based on the assumption that the compatibility and equilibrium conditions, arising from sliding contact, are a priori known. In that case, the time dependency is obtained by assuming a rigid relative motion between component subsystems due to the time-dependent boundary conditions. In this paper, a better estimation is obtained by giving up the assumption of rigid relative motion and accounting also for system deformation. The system position U at time t n+1 is estimated using the predictor of displacements u n+1 , i.e., where X represents the undeformed position. Since the system deformation affects the coupling conditions, the coupled system is geometrically nonlinear. Moreover, by considering system deformation, the variation of the tangential displacement at the contact can be observed, thus allowing to account for friction-induced vibrations. Finally, by substituting λ n in Eq. (29), the incremental displacement u c n+1 due to the contact forces is evaluated and the displacement u n+1 is obtained.

Contact algorithm
To enforce compatibility and equilibrium, the position of the contact point must be defined. To this aim, a contact algorithm is proposed according to the well-known master element-slave node approach [32]. To identify the system configuration at time t n+1 , the system position U n+1 , in Eq. (33), is used. Figure 1 shows a target element defined by nodes T 1 and T 2 and a contact node. The position of contact node C is given only by the predictor; hence, it does not necessarily lie on the target element. It is assumed that the position of contact point C is given by the normal projection of the node C on the target element. The distance between node T 1 and point C can be expressed as: where l 12 is the length of the target element, i.e., the distance between nodes T 1 and T 2 . The target element in contact is the one that satisfies the following conditions: where l 1C is the distance between nodes T 1 and C . Finally, the adimensional relative position α of the slave node C along the target element is given by:

Compatibility and equilibrium for dual assembly
Using dual assembly, all the elements of u are retained. In order to evaluate the time response of the system, the explicit time integration method in Sect. 2.2.1 requires the time-variant matrices B C (t) and B E (t) to be defined. In fact, they are needed to evaluate the contact forces λ, in Eq. (32), and the incremental displacement u c associated with the contact forces, in Eq. (29). Figure 2a shows the connecting forces acting on the target element and on the contact node. The equilibrium of forces shown in Fig. 2a along the x and y directions and around node T 1 is: For a given friction coefficient μ > 0 and a relative velocity Fig. 2b shows the normal and tangential component, respectively, F N and F T , of the contact force between the contact node and the target element.
The directions of F N and F T in Fig. 2b are those arising when the two bodies are pressed against each other, and the relative velocity −→ v r is directed as shown. The connecting force components g Cx and g Cy on the contact nodes can be expressed by the projection of F N and F T on the x and y directions.
where the ± accounts for the inversion of relative velocity −→ v r that causes the inversion of the tangential force F T too.
If a constant Coulomb friction law is considered, i.e., |F T | = μ|F N |, the connecting forces on the contact node can be expressed as: where Γ 1 and Γ 2 collect all the friction and geometric parameters. By combining Eqs. (36), (37) and (39), the nonzero elements of the connecting force vector g can be expressed as: Hence, by considering λ = F N , the set of equations (40) can be rewritten in matrix form as: and according to Eq. (17), the vector on the right-hand side of Eq. (41) contains the nonzero elements of B E T .
As stated in Sect. 2.2, the matrix B E is different from the matrix B C because it accounts for friction forces at the interface, while the compatibility is enforced only along the direction normal to the contact (i.e., normal to the target element). However, for μ = 0 the two matrices are the same as in the collocated problem, The factors Γ for a nil friction coefficient are: By substituting the previous equation into the Eq. (42), the nonzero elements of the matrix B C T , used to enforce compatibility, are: According to the definition of λ before Eq. (41), the contact assumption is verified, if the normal contact force F N acting on node C is positive, that is if λ > 0.

Time-dependent frequency response function
The time-dependent frequency response function, introduced by the authors in [15], gives an idea of the influence of the time-dependent interface on the dynamics of the system at different frequencies. For systems composed of time-invariant substructures with time-variant contact interfaces, the matrices B C and B E are configuration dependent, i.e., they depend on the relative position of the contacting bodies. If the elastic deformation of the system is neglected, the matrices B C (t) and B E (t) can be estimated a priori as function only of the time-dependent boundary conditions [15]. The time-dependent frequency response function matrix is expressed as [15]: It provides, with a much smaller computational effort than that required to compute the time response, a result that is qualitatively equivalent to the time-frequency analysis of the linear system response.

Application
The approach proposed in the previous section is here applied to a beam-on-beam system, similar to that described in [27]. Beam-on-beam models are often used as benchmark systems for contact problems. The objective of the study is the verification of the ability of the substructuring-based procedure to cope with friction-induced vibrations. The system is relatively simple, but the procedure here proposed could be applied to systems composed of more complex substructures to deal with engineering applications. In fact, engineering systems with sliding contact interfaces are often affected by friction-induced vibrations. Therefore, numerical methods must necessarily be able to reproduce these phenomena in order to be able to reproduce the actual dynamic behavior of mechanical system with sliding contact interface.
3.1 Numerical model The mechanical system shown in Fig. 3a is made by two contacting beams forming an angle ϑ of 30 • . The horizontal beam 2 is made of polycarbonate and is fixed at the left end B. The oblique beam 1 is made of aluminum alloy. The contact between the beams at point C is ensured by a vertical load F y applied at the upper end A of beam 1 , whereas the relative motion between the two beams is provided by setting the upper end A in motion with a horizontal velocity v x . The amount of the horizontal velocity is selected to be much lower than 15.24 m/s, that is the critical speed of a moving load in the horizontal direction. A sliding friction coefficient μ is assumed between the two beams. Furthermore, the initial position of the contact point C is denoted as s 0 . Geometrical dimensions, mechanical properties and boundary conditions are shown in Tables 1 and 2. Different values of friction coefficient μ are used in order to induce both stable and unstable system responses.
Moreover, in order to simulate roughness at the contact interface, a vertical force-defined as a bandlimited white noise-is applied on point C. The two beams are modeled by 350 plane stress elements, 531 nodes and 1062 DoFs. The substructure 1 has 456 DoFs, while the substructure 2 has 606 DoFs. To reduce the computational effort, a modal reduction is performed on the two substructures. Width, a (mm) 10 20 Angle, ϑ (•) 3 0 0  (Fig. 3b). Furthermore, 20 fixed interface modes are selected for substructure 1 (including modes up to 150 kHz) and 20 free interface modes are selected for substructure 2 . Overall, after the reduction, substructure 1 has 28 DoFs (u 1 ) and substructure 2 has 220 DoFs (u 2 ). To solve the contact problem, master elements are the top edges of the elements of the horizontal beam 2 and the slave node is the bottom corner C of the inclined beam 1 .

Results
The substructuring method presented in Sect. 2 is used to investigate the dynamic behavior of the beam-onbeam system in different conditions. Once the geometrical characteristics and mechanical properties of the system have been selected, the stability of the system may basically depend on the friction coefficient and on the position of the contact point. Generally when dealing with time-or positiondependent dynamic behavior, a system is stable when the real parts of its eigenvalues are negative at any time instant or at any position. Conversely, a system is unstable when the real parts of its eigenvalues become positive during a time interval or within a range of positions. In the proposed case study, the system's configuration changes with time due to the displacement of the contact point from the fixed end to the free end of the horizontal beam. The continuous configuration change is sampled to obtain a finite set of configurations. For each configuration, the two substructures are coupled using primal assembly as shown in [15] and complex eigenvalue analysis is performed, providing the corresponding complex eigenvalues. Therefore, the definition of stable or unstable system is extended to the whole set of considered configurations. If the real parts of the system's eigenvalues are all negative, the system is stable; otherwise, if at least one of the system's eigenvalues has a positive real part for at least one configuration, the system is unstable.   Figure 4 shows the locus plots of the eigenvalues of the coupled system, corresponding to all the sampled positions of the contact node C during the sliding from the constrained end B to the free end T of the horizontal beam. The variation of each complex eigenvalue due to the displacement of the contact node C along the horizontal beam is represented with a colored line on the complex plane. Dotted lines represent the locus of the complex eigenvalues of each subsystem: They depend only on proportional viscous damping, i.e., on the values of coefficients α and β in Table 2, and the eigenvalues of the uncoupled subsystems would lie on that lines. It is expected that the eigenvalues of the coupled system lie between the dotted lines unless mode coupling is present.
Two typical results are shown in Fig. 4, one corresponding to a friction coefficient μ = 0.1 and the other corresponding to a friction coefficient μ = 0.3. The locus plot of the eigenvalues highlights that all eigenvalues have negative real part for μ = 0.1 (stable behavior), while the real part of some eigenvalues may reach positive values when μ = 0.3 (possible instabilities at some positions of the contact point). Specifically, instabilities may occur around 7 kHz and 10.5 kHz.
Hence, transient analyses have been performed for two different values of the friction coefficient μ in order to show stable or unstable behavior. Moreover, analyses have been performed both without and with roughness. When roughness is present on the upper edge of the beam, the contact node C follows the asperities of the profile during the sliding. Hence, vertical acceleration arises at the contact interface and inertia forces appear on both substructures that cause a variation of the normal contact force to ensure the dynamic equilibrium. Therefore, roughness is simulated as a force acting on the contact node C in vertical direction having a band-limited white spectrum with standard deviation σ as a fraction of the vertical preload F y [27]. In the following, four different cases are considered: 1. Stable smooth, μ = 0.1 and σ = 0 2. Stable rough, μ = 0.1 and σ = 0.5% of F y 3. Unstable smooth, μ = 0.3 and σ = 0 4. Unstable rough, μ = 0.3 and σ = 0.5% of F y . Figure 5 shows the vertical accelerations of the free end of the horizontal beam in the four considered cases. When the friction coefficient is low, the system is stable and the accelerations in Fig. 5a remain on limited values for both the smooth and the rough case. For a higher value of the friction coefficient (μ = 0.3), during sliding the system runs through instability zones, the acceleration in Fig. 5b reaches the limit values and then decreases when the contact point changes its position. The maximum acceleration depends both on the presence of roughness, modeled as a wide-band random force, and on the amount of the preload F y . Figure 6 shows the spectrogram of the vertical acceleration of the free end of the horizontal beam for the two cases with roughness. The roughness acts as a wideband excitation, and the spectrogram in Fig. 6a shows, for the stable case, the resonances of the system and how they change when the contact point travels along the horizontal beam. For the case in Fig. 6b, besides the response to the wide-band excitation, the spectrogram highlights the frequency and time intervals within which the response increases because the system runs through instabilities. The spectrogram of the response highlights two instabilities of the system, one at around   Fig. 4b, and another one at 3.5 kHz (at around 2.8 s) that was not identified by the complex eigenvalue analysis. Recalling that complex eigenvalue analysis is linear, i.e., it considers the coupled system in the undeformed position, the effect can be ascribed to the geometric nonlinearities due to the deformation of the system, because it modifies the direction of the normal and tangential contact forces. Furthermore, no instability can be observed around 7 kHz probably because the real part of the eigenvalue is too low to trigger it. Figure 7 shows the spectrum of the vertical acceleration of the free end of the horizontal beam in the unstable rough case for three time intervals corresponding to three bursts of the response in Fig. 5b. The results agree with those observed in the spectrogram. Furthermore, during the last burst it can be clearly observed that the frequency around 10.5 kHz appears in the spectrum besides the main frequency around 3.5 kHz.    relative velocity crosses the zero; (ii) if detachment tends to occur, i.e., if λ < 0. The horizontal velocity of the contact point C is a good estimation of the relative velocity. For the unstable case with roughness, the horizontal velocity of the contact point C is shown in Fig. 8a and the normal contact force is shown in Fig. 8b. The normal contact force is positive during the entire simulation time, i.e., the contact hypothesis is always verified and no detachment occurs. Moreover, the horizontal velocity of the contact point oscillates around the value of 30 mm/s, that is the horizontal velocity v x imposed on the upper beam, and never changes its sign. It means that the sliding occurs always in the same direction and that unstable vibrations do not produce an inversion of motion.
Finally, Eq. (45) is used to compute time-dependent FRFs of the coupled system when the upper beam is sliding along the horizontal beam with the same speed v x as in the transient simulation and a friction coeffi-  Figure 9 shows the time-dependent FRFs by considering an input force acting at the contact point C in vertical direction and the vertical acceleration of the tip T of the horizontal beam. The results show a pattern with reference to the high levels that is very similar to the pattern of the two spectrograms in Fig. 6a, associated with the response of the stable system to the wide-band excitation. Furthermore, the time-dependent FRF highlights a low-level pattern that accounts for time-dependent antiresonance locations. Nevertheless, the time-dependent frequency response function computed for a friction coefficient μ = 0.3 is very similar to that shown in Fig. 9 and it does not provide clear information about the instabilities of the system: For this reason, it is not shown.

Convergence check of the solution
In order to check the convergence of the obtained results, the response has been compared with that obtained using a finer discretization of the horizontal beam. The refinement is applied only on the finite element model of the horizontal beam. In fact, the Craig-Bampton reduction applied to the inclined beam model maintains only its end nodes and the first 20 fixed interface modes, that are derived from a finite element model that ensures an adequate spatial discretization. Therefore, a finer discretization of the inclined beam would not significantly modify its reduced model. Conversely, the increase in the number of elements used to     discretize the horizontal beam results in an increased number of target elements of the reduced model, thus affecting the proposed numerical method, whose convergence is hereafter verified. The results previously discussed, obtained using 100 target elements (i.e., using 100 plane elements along the longitudinal direction of the beam), are here compared with the results of the same simulation using 150 target elements. The time step adopted to solve the analysis with the finer model must be reduced in order to respect the stability condition in (25). The comparison is performed considering the "Unstable smooth" case, i.e., friction coefficient μ = 0.3 and absence of roughnessrelated contact noise. Figure 10 compares vertical acceleration of the free end of the horizontal beam obtained using the two mod-els. The results highlight a substantial similarity of the two responses: Except at t 1.5 s, the vibration bursts associated with the crossing of the instability regions appear at the same time and have very similar amplitude.
However, the small burst around 1.5 s is not very significant, as confirmed by the time-frequency analysis of the acceleration responses in Fig. 11, where the frequencies of the main bursts in the two spectrograms are the same.

Concluding remarks
In this paper, the general framework for dynamic substructuring is extended to time-variant interfaces among coupled substructures in order to observe the onset of friction-induced vibrations. Specifically, a time-variant interface due to a sliding contact is considered that is able to take into account the oscillations of the contact point in the formulation of the compatibility and equilibrium conditions used in the transient analysis. The proposed solution approach in time domain uses a nonlinear forward increment dual assembly, adapting the forward increment Lagrange multiplier method to the substructuring framework. The main feature of the proposed method is the definition of the compatibility and equilibrium matrices that are estimated according to the framework of dynamic substructuring but accounting for the deformation of the contacting bodies. This technique highlights some mode coupling instabilities that are not detected using the linear timedependent complex eigenvalue analysis because they are due to geometric nonlinearities. At this stage of development, the proposed method is valid only if the substructures are in relative sliding. The sliding condition is verified during the simulation, and if sticking or detachment occurs the simulation is stopped. Results on the beam-on-beam system described by finite element model highlight the effectiveness of the proposed approach in both stable and unstable scenarios, and the possible application to more complex systems.