An enhanced formulation to model spatial revolute joints with radial and axial clearances

Mechanism and Machine Theory 116 (2017) 123–144

Contents lists available at ScienceDirect

Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmachtheory

Research paper

An enhanced formulation to model spatial revolute joints with radial and axial clearances Filipe Marques∗, Fernando Isaac, Nuno Dourado, Paulo Flores CMEMS-UMinho, Departamento de Engenharia Mecânica, Universidade do Minho, Campus de Azurém, Guimarães 4804-533, Portugal

a r t i c l e

i n f o

Article history: Received 20 February 2017 Revised 28 April 2017 Accepted 21 May 2017

Keywords: Revolute clearance joints Contact forces Friction forces Multibody dynamics

a b s t r a c t The primary objective of this work is to present a new formulation to model spatial revolute joints with radial and axial clearances. The methodology developed is based on the relative motion between the journal and bearing elements, which accounts for all the possible contact scenarios within a revolute clearance joint. In this process, normal and tangential contact force models are reexamined to evaluate their suitability to represent the contact-impact interactions between the mechanical parts connected by clearance joints. Also the equations of motion that govern the dynamic response of multibody systems and incorporate the intra-joint contact forces are presented. The Newton-Euler approach is applied here for these procedures. Finally, several numerical examples of the application are presented to discuss the main assumptions and procedures adopted in this work. The obtained results show that the clearance in mechanical joints can be quite significant in terms of the dynamic response of the system, and, consequently, plays a crucial role in the analysis, design and control of multibody systems. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Many authors agree that clearances in mechanical joints may affect performance of the mechanism significantly [1–6]. Extremely tight clearances in mechanical joints can make the performance close to the desired one; however the cost of such a mechanism may be prohibitive. Thus, a designer may wish to assign a clearance such that the output response is within the specified tolerances, and at the same time the clearance is as tight as possible [7]. Thus estimating the effects of clearances on the output characteristics are a prerequisite for proper analysis, design, control and optimization of multibody systems [8]. Over the last few decades, various researchers have investigated problems of modeling and analyzing mechanical systems with clearance joints. However, most of the published works focus on planar revolute joints [9–15], planar translational joints with clearance [16–20], spherical clearance joints [21–26], and cylindrical joints with clearance [27–30]. In sharp contrast, little research has been done on modeling spatial revolute joints with clearance. Both radial and axial clearances in spatial revolute joints play a crucial role in various applications, such as vehicle systems and components [31], robotic manipulators [32], space deployable systems [33], and industrial machinery [34], just to mention a few. In fact, even planar mechanical systems can exhibit out-of-plane motion due to misalignments, and thus justify the development of enhanced mathematical models to assess the influence of the clearance in spatial revolute joints that include both radial and axial gaps. ∗

Corresponding author. E-mail address: [email protected] (F. Marques).

http://dx.doi.org/10.1016/j.mechmachtheory.2017.05.020 0094-114X/© 2017 Elsevier Ltd. All rights reserved.

124

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

In fact, the problem of modeling and implementation issues associated with spatial revolute joints with clearance is a complex and demanding task, since it involves a large number of different contact scenarios between the journal and bearing elements [35]. Furthermore, particular attention must be given to the contact detection, as well as to the transition between the different contact scenarios. These aspects are of paramount importance in terms of accuracy and efficiency when dealing with the computational resolution of the equations of motion for multibody mechanical systems [36]. Dhande and Chakraborty [37] developed an analytical procedure to analyze the output error of function-generating spatial mechanisms considering random clearances in different mechanical joints. These authors were pioneers in investigating spatial joints with clearance, and the methodology proposed was illustrated with RSSR and RRSS spatial mechanisms; however, the wobbling of the pin was not taken into account in those investigations. In turn, Rhyu and Kwak [38] proposed a formulation for an optimal stochastic design of a four-bar mechanism considering tolerances on the lengths of the links and clearances at the joints in order to analyze the mechanical output error for different input angles. These authors implemented a multi-objective optimization procedure in which the objective functions were the mechanical error and the manufacturing cost. A four-bar function and a four-bar path generator were used to demonstrate the application. Dubowsky et al. [39] were also pioneers in the study of the spatial revolute joints with clearance. In their investigation, these authors presented a formulation for three-dimensional dynamics of spatial machine systems taking into account also the vibrations in the links, the supporting structures and the impacts at the clearance joints using a finite element approach. The results were used to establish design guidelines for high-performance machine systems and minimization of vibrations and noise. Later on, Kakizaki et al. [40] presented a model for spatial dynamics of robotic manipulators with flexible links and clearance joints, where the effect of the clearance was considered to control the robotic system. The results showed that the combined actions of clearance, flexibility of links and control system characteristics considerably affected the performance of the system. Deck and Dubowsky [41] theoretically and experimentally investigated the dynamic behavior of mechanical systems with flexible links and clearance joints. Again, the outcomes clearly showed the very high sensitivity of the joint operating parameters and large impact forces of industrial machines with clearance joints. Also, the flexibility of the links reduced the vibrations and impacts. Brutti et al. [42] are among the few investigators who modeled and analyzed spatial clearance joints with or without lubrication [43]. This type of investigation improved the accuracy of the dynamic performance evaluation of a spatial slider-crank mechanism. More recently, Yan et al. [44] proposed a comprehensive model for spatial revolute joints with clearance in multibody systems. These authors investigated a single journal-bearing system to reveal the main characteristics of the journal motion inside the bearing element. Liu et al. [45] presented a compliant force model for cylindrical joints with clearances, where Hertz’s law is only valid for large clearance sizes and small loads. This approach was compared and validated with results from FEM analyses [46]. Problems of modeling and analyzing mechanical joints with clearance have been also investigated by other researchers, namely Bauchau et al. [47,48], Parenti-Castelli et al. [49], Venanzi et al. [50], Erkaya et al. [51], Tian et al. [52] and Isaac and his co-authors [53]. The main goal of this work is to address the issues associated with the operating conditions of spatial revolute joints with radial and axial clearances that complement the ones proposed in the thematic literature. Basically, a journal inside a bearing is a typical composition of a revolute joint. Due to the necessary and unavoidable clearances, the journal can freely move in both radial and axial directions. In the present work, the journal and bearing elements are considered to be rigid and collide with each other, being the behavior of the joint is governed by contact-impact forces. In this process, the development of appropriate methodologies to handle the contact detection between the journal and bearing surfaces is of paramount importance and play a key role in the accuracy and efficiency of the response of the system. In this study, the intra-joint contact forces are based on the geometric and kinematic properties of the clearance joint, as well as on the characteristics of the contacting materials. These aspects represent the main core of this investigation. This paper is organized as follows. Section 2 presents a full characterization of spatial revolute clearance joints in the context of multibody systems formulation. The potential contact scenarios that can occur in this type of mechanical joint are also described in this section. Section 3 presents a general and comprehensive methodology to deal with the contact detection between the elements that constitute a spatial revolute joint with clearance. Section 4 describes the most suitable normal and tangential force models used to evaluate the intra-joint contact forces. In Section 5, the main ingredients to generate the equations of motion for constrained spatial multibody systems are revised. In Section 6, several numerical results from computational simulations of spatial multibody systems including revolute clearance joints are used to demonstrate the effectiveness of the proposed approach. Finally, the Summary and Conclusions drawn from this research work are presented in Section 7.

2. Spatial revolute clearance joint Fig. 1 shows a schematic representation of a typical dry revolute joint with clearance. The mechanical elements that compose this type of joint are the bearing or bushing and the journal or pin. In Fig. 1, the bearing is part of body i, while the journal is part of body j. In the present work, it is assumed that radial and axial clearances are of the same order of magnitude. Furthermore, in actual revolute joints the clearances are much smaller than the radii and lengths of the bearing and journal elements. In Fig. 1, Ri and Rj represent the radius of the bearing and of the journal, respectively, while Li and

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

125

Fig. 1. Schematic representation of a spatial revolute joint with radial and axial clearances.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Fig. 2. Different contact scenarios between the journal and bearing elements in a spatial revolute clearance joint.

Lj denote their corresponding lengths. The subscripts i and j refer to the bearing and journal, respectively, throughout this text. Both radial and axial displacements are possible in a spatial revolute joint with clearance. This is of paramount importance in practical mechanisms and machines due to misalignments, local elastic and plastic deformations, wear and thermal effects that can occur in the mechanical joints. Fig. 1 demonstrates that the radial and axial clearances can be expressed as follows

cr = Ri − R j

(1)

Li − L j 2

(2)

ca =

In the nominal and free flight motion, there is no contact between the journal and bearing elements. However, depending on the size of the radial and axial clearances, several different contact scenarios between the bearing and journal surfaces can take place, some of which are illustrated in Fig. 2 [54]. The contact configurations used in this study are as follows: (i) The journal contacts the bearing wall at a single point, as illustrated in Fig. 2(a) and (b);

126

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

Fig. 3. Nominal position of a spatial revolute joint with clearance.

(ii) (iii) (iv) (v) (vi) (vii)

The journal and bearing contact with each other at two different points, as shown in Fig. 2(c) up to (f); The journal contacts the bearing wall at three points, as shown in Fig. 2(g) and (h); The journal contacts the bearing wall at four different points, as represented in Fig. 2(i); The journal and the bearing elements contact each other along a line, as Fig. 2(j) illustrates; The bottom of the journal contacts with the bottom of the bearing, as shown in Fig. 2(k); A combination of the two previous scenarios can be seen in Fig. 2(l).

The 12 different contact possibilities between the journal and bearing elements illustrated in Fig. 2 affect the behavior of the joint and depend on the configuration of the dynamic system as well as on the radial and axial clearances. However, some of the contact cases presented in Fig. 2 may not occur in practice, since they are directly related to the dimensions of the journal and the bearing. Therefore, the size of the bodies may introduce some restrictions which do not allow the joint to have some particular contact configurations. Consequently, the methodology developed to model spatial revolute joints with clearance must be able to handle the potential contact configurations and the transition between them [36]. In the free flight or non-contact scenario, no forces are developed in the joint, because the journal moves freely inside the bearing boundaries. When the journal and bearing are in contact with each other, local deformations take place at the contact zone and, consequently, contact forces characterize the interaction between the two parts. The response of the system during the period of contact is obtained by adding the contact forces to the equations of motion of the multibody system as external generalized forces. This approach can provide accurate results, in the sense that the equations of motion are integrated over the period of contact, and take into account the changes in the configuration and kinematics of the system during that contact [35].

3. Contact detection methodology This section is devoted to a detailed description of the proposed methodology for contact detection between the journal and bearing elements in a spatial revolute joint with clearance. The procedure to determine the contact points is presented for all the previous contact scenarios mentioned. Fig. 3 shows a generic configuration of a revolute joint with clearance in which the bearing and the journal are in the nominal position, that is, fully aligned. In order to evaluate the location and orientation of these elements in each body, a point P located in the center of the cylinder axis and a unit vector a with its orientation are established, as depicted in Fig. 3. The vector ri defines the position of the mass center of body i in the global coordinate system, and si P is the position of P with respect to the body i in the same system [55]. A similar notation is used for the remaining bodies and points. In order to simplify the contact detection procedure, auxiliary points located in the center of each extremity of both elements are used, as shown in Fig. 3. Moreover, the top and bottom bases are nominated A and B, and the coordinates of these points can be determined according to the following equations

rAk = rPk +

Lk a , (k=i, j ) 2 k

(3)

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

127

Fig. 4. Configuration of the aligned spatial revolute joint with axial and radial contact scenario.

rBk = rPk −

Lk a , ( k =i , j ) 2 k

(4)

For the sake of simplicity, the contact scenarios described above might be divided into aligned (Fig. 2(j)–(l)) and misaligned cases (Fig. 2(a)–(i)). The alignment of the journal and bearing is evaluated by the parallelism of their axes, that is, their unit vectors must remain parallel, which can be expressed as

a˜ i a j = 0

(5)

where ã denotes the skew-symmetric matrix associated with the vector a [55]. Taking the aligned case into consideration, an arbitrary configuration where axial and radial contacts occur is shown in Fig. 4, and results in superficial and linear contact zones. The contact points must represent these zones and must be located at their geometrical center. These points are denoted by C, and the subscripts i and j stand for bearing and journal, respectively, while the superscripts a and r denote axial and radial contacts. The superscript also includes information relative to the extremity in which the contact occurs. For instance, Ci B,a represents the contact point of the bearing for axial contact at base B. Fig. 4 shows that the mid-point of the radial contact line is located at half way along the length of the journal. Moreover, the normal unit vector can be defined as follows

nv =

P   T r j − rPi − rPj − rPi ai ai    T   Pj   r − rPi − rPj − rPi ai ai 

(6)

Thus, the contact points can be determined using the following expressions



r

rCi = rPi + rPj − rPi

T

ai ai + Ri nv

r

rC j = rPj + R j nv

(7) (8)

The contact zone of the axial contact is characterized by a circle where the center point coincides with the center of the base of the journal. Therefore, considering the contact at base B, the coordinates of the contact points are given by

rCi



B,a

= rB j + rBi − rB j

B,a

= rB j

rC j

T

ai ai

(9) (10)

When the journal and bearing are not aligned, the contact detection is a more complex task as can be seen in Fig. 2(a)– (i), where the misaligned contact scenarios cause different points of contact. Thus, each extremity of the journal can collide with the lateral and top walls of the bearing. A generic misaligned revolute joint with axial and radial contact is represented in Fig. 5 to describe the proposed methodology.

128

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

Fig. 5. Configuration of the misaligned spatial revolute joint with axial and radial contact scenario.

The axial contact case here can be defined as the intersection of a circumference (base of the journal) with a plane (base of the bearing). These two geometric entities can be expressed as

  x − rB j 2 = R j 2

(11)

a j T x = a j T rB j ai T x = ai T rBi

(12)

where x represents the coordinates of a point belonging to each geometric element. Eqs. (11) and (12) can be solved together in order to obtain the intersection points. The solution of the intersection between a circumference and a plane can have zero, one or two points which represents, respectively, no contact, tangent contact or penetration. If the solution is given by two points, I1 and I2 , then the middle point and the center of the base of the journal form an auxiliary vector that will be used to compute the position of the contact point as

nd =

I  1 r 1 + rI2 − rB j  21     rI1 + rI2 − rB j 

(13)

2

Thus, the coordinates of the contact point in the journal can be given by

rC j

B,a

= rB j ± R j nd

(14)

where the sign must be chosen in such a way that it returns a point outside the bearing. Hence, the contact point in the bearing can be calculated as follows

rCi

B,a

= rC j

B,a



+ rBi − rC j

B,a

T

ai ai

(15)

The radial contact can be considered as the intersection of a circumference (base of the journal) with a cylindrical surface (lateral wall of the bearing), which are defined, respectively, as

  x − rA j 2 = R j 2 a j T x = a j T rA j

   a˜ i x − rPi 2 = Ri 2

(16) (17)

In contrast to the previous case, this set of equations form a nonlinear system, which cannot be solved analytically. Therefore, the intersection points must be obtained using an iterative procedure, such as Newton-Raphson algorithm. From the numerical point of view, the intersection between an arbitrary circumference and a cylindrical surface can have from zero to four points. Considering the physical restrictions of this problem, that is, the angle between the bearing and the journal is small, so the maximum number of solutions is reduced to two. Similarly to the axial contact detection, when penetration exists, an auxiliary vector is calculated using the following expression

nd =

I  1 r 1 + rI2 − rA j  21     rI1 + rI2 − rA j  2

(18)

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

129

Thus, the contact point in the journal can be obtained, in an analogous manner, yielding

rC j

A,r

= rA j ± R j nd

(19)

Again, the sign to be adopted must lead to a point outside the bearing. The normal unit vector can be expressed as:



nv =

rC j

A,r





− rPi − rC j

A,r

− rPi

T

ai ai

   T   A,r   rC j − rPi − rC j A,r − rPi ai ai   

(20)

Finally, the coordinates of the contact point in the bearing can be given as:

rCi

A,r



= rPi + rC j

A,r

− rPi

T

ai ai + Ri nv

(21)

The contact points have been evaluated according to the particular case of Fig. 5, that is, axial contact at base B and radial contact at base A. Similar analyses can be performed for the remaining contact cases. At every time step, the contact is evaluated for all the possible interactions, radial and axial at extremities A and B, and the contact case is defined through the effectively existing contacts. The contact kinematics for each collision can be characterized through the contact points in the bearing and the journal, respectively, Ci and Cj . Thus, the normal unit vector of the contact can be evaluated as follows:

nv =

C  r j − rCi    rC j − rCi 

(22)

The penetration depth, δ , can be expressed as:

  δ =  rC j − rCi 

(23)

In turn, the penetration or contact velocity, δ˙ , can be determined as

 T δ˙ = r˙ C j − r˙ Ci nv

where as [55]

r˙ Ci

Cj

and r˙

(24)

denote the translational velocities of the bearing and journal at the contact points which can be calculated

˜ k skCk , (k=i, j ) r˙ Ck = r˙ k + ω

(25)

˜ denotes the skew-symmetric matrix obtained from the angular velocity vector. in which ω The contact velocity determines whether the colliding points are approaching the contact surface or separating from each other. Finally, the relative tangential velocity necessary to evaluate the friction forces can be expressed as







vT = r˙ C j − r˙ Ci − r˙ C j − r˙ Ci

T

nv nv

(26)

4. Evaluation of the contact forces The contact-impact forces developed during the collisions between the journal and bearing elements play a fundamental role in the dynamic behavior of the system. Thus, in this section the normal and tangential contact force models used in this work are briefly described. The models selected take into account the geometric and kinematic characteristics of the joint, as well as the material and physical properties of the contacting parts. Furthermore, the contact force models used contribute to the stable resolution of the equations of motion [56]. In what concerns to the normal contact force evaluation, a model based on Hertzian contact theory together with a dissipative term is considered [57], and can be expressed as [58]

FN = K δ n + χδ n δ˙

(27)

where K is the generalized contact stiffness which depends on the materials and geometries in contact, χ represents the hysteresis damping factor, and n is an exponent that defines the degree of nonlinearity of the contact. The exponent n is equal to 1.5 for metallic contacts. In the context of this work, the contact force model proposed by Lankarani and Nikravesh is employed [59], where the hysteresis damping factor is defined as



χ=

3 1 − ce 2 4



K

(28)

δ˙ (−)

where ce is the coefficient of restitution and δ˙ (− ) denotes the initial impact velocity. Introducing Eq. (28) into Eq. (27), the generic expression to calculate the normal contact force yields



FN = K δ n 1 +



3 1 − ce 2 4



δ˙ δ˙ (−)



(29)

130

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

Alternative contact force models are available in the thematic literature, and can be found in the works by Machado et al. [60] and Banerjee et al. [61]. A wide range of friction force models has been developed due to the significant complexity of frictional phenomena. These models can be divided into two main groups, namely the static and dynamic friction models [62]. The most wellknown friction force model was presented by Coulomb [63], who stated that friction opposes the relative motion between two contacting surfaces, and the magnitude of the friction force is proportional to the normal contact force (FN ). This friction force model is dependent on the direction of the relative velocity (vT ), and can be described as

fF = FC sgn(vT )

(30)

FC = μk FN

(31)

where

in which μk denotes the kinetic coefficient of friction, sgn is a function which returns a unit vector with the direction of the relative velocity, and fF represents the friction force vector. Although this friction model is quite simple and straightforward, it exhibits some difficulties since it does not specify a friction force at zero velocity. This particular issue causes a discontinuity to determine friction forces leading to perturbations in the dynamic response of the system. In order to simplify and ensure computational efficiency, several researchers have proposed alternative methods, which replace the discontinuity at zero velocity by a finite slope [62]. An alternative approach is to use a hyperbolic tangent approximation, which is numerically stable because it has a continuous derivative. Thus, the friction force model can be expressed as

fF = FC tanh (kvT  )sgn(vT )

(32)

where k is a coefficient that defines the derivative for null relative velocity. A significant number of friction force models can be found in the literature, namely in the works by Marques et al. [62], Pennestrì et al. [64] and Pichler et al. [65]. 5. Equations of motion formulation A multibody system can be defined as a collection of bodies that are acted upon by generalized forces. The bodies are interconnected to each other by different types of kinematic pairs that constrain their relative motion in different forms. For a constrained multibody system, the kinematic joints can be described by a set of algebraic equations in the following form [55]

(q, t ) = 0

(33)

where q represents the generalized coordinates vector and t is the time variable. Each ideal or perfect spatial revolute joint contributes with five equations to the system of Eq. (33). Differentiating Eq. (33) with respect to time yields the velocity constraint equations. After a second differentiation with respect to time, the acceleration constraint equations are obtained as

q q¨ = γ

(34)

where q is the Jacobian matrix of the constraint equations, q¨ is the acceleration vector and γ is the right hand side of the acceleration equations, which contains the terms that are exclusively function of velocity, position and time [55]. In terms of the Cartesian coordinates, the translational and rotational equations of motion of an unconstrained multibody system are written as

Mq¨ = g

(35)

where M denotes the mass matrix of the global system, containing the mass and moments of inertia of all bodies and g is the generalized force vector that contains all the external forces and moments acting on the system including the intra-joint contact force developed at clearance joints. Using the Lagrange multipliers technique, the constraint Eq. (33) are added to the equations of motion (35). The equations of motion are written together with the second time derivative of constraint Eq. (34) yielding a system of equations written as



M

q

Tq



0

q¨

λ



=

g

γ

(36)

where λ is the vector of the Lagrange multipliers, which are physically related to the ideal joint reaction forces. The reaction forces, owing to the kinematic joints are expressed as [66].

g(c ) = −Tq λ

(37)

Eq. (36) represents a differential-algebraic equation that has to be solved and the resulting accelerations integrated in time. However, they do not use explicitly the position and velocity constraint equations allowing for a drift to develop

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

131

z

x

y

Fig. 6. Multibody model of the slider-crank mechanism with planar motion.

Table 1 Dimensional and inertia properties of each body of the slider-crank mechanism. Body

Length [m]

Mass [kg]

Principal Moments of Inertia [kg m2 ] Iξ ξ

Iηη

Iζ ζ

Crank Connecting Rod Slider

0.1524 0.3048 –

0.15 0.30 0.15

0.002 0.002 0.001

0.002 0.002 0.001

0.002 0.002 0.001

in the system constraints. In order to keep such constraint violations under control during the numerical integration, the Baumgarte stabilization technique is applied, and Eq. (36) is modified to



M

q

Tq 0

 q¨

λ



=

γ



g ˙ − β 2 − 2α 

(38)

where α and β are prescribed as positive constants that represent the feedback control parameters for the velocity and position constraint violations [67,68]. According to this formulation, the dynamic response of multibody systems involves the evaluation of the vectors g and γ , for each time step. Then, Eq. (38) is numerically solved for the system accelerations. These accelerations together with the velocities are integrated in order to obtain the new velocities and positions for the following time step. This process is repeated until the complete description of system motion is performed [55,69]. Ideal or perfect joints are introduced in Eq. (38) as kinematic constraints. However, joints with clearances are introduced in Eq. (38) as external forces on the right-hand-side. 6. Results and discussion The main goal of this section is to assess the validity of the proposed formulation for spatial revolute joints with clearance. Therefore, two different slider-crank mechanisms were used as examples to demonstrate the validity of the approach used here. 6.1. Slider-crank mechanism describing planar motion The first example of the application of this approach consists of a 3D slider-crank mechanism that describes a planar motion, and is represented in Fig. 6. Although the movement described by the mechanism is planar, the dynamic analysis of the multibody system is modeled through a spatial formulation. This multibody system has three moving bodies, namely the crank, the connecting rod and the slider. The bodies are connected through perfect kinematic joints, except for the case of the connection between the connecting rod and the slider, which is made by a revolute clearance joint. For the remaining ones, revolute, spherical and translational kinematic joints are considered for the couplings groundcrank, crank-rod, and slider-ground, respectively. The dimensional and inertia properties of each body are listed in Table 1. Moreover, the location of the center of mass is considered to be at the geometric center of each body as depicted in Fig. 6. In the initial configuration of the mechanism, the crank and connecting rod are aligned with each other, and the slider is located in the final stroke position. The initial angular velocity of the crank is 150 rad/s, and the remaining velocities are determined in order to fulfill the corresponding constraint equations. In the present analysis, two sets of initial conditions were

132

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144 Table 2 Parameters for the dynamic simulation with the revolute clearance joint. Parameter

Value

Parameter

Value

Radius of bearing, Ri Radius of journal, Rj Radial clearance, cr Length of bearing, Li Length of journal, Lj Axial clearance, ca Contact stiffness, K Coefficient of restitution, ce

10 mm 9.9 mm 0.1 mm 20.0 mm 19.8 mm 0.1 mm 1011 N/m3/2 0.9

Coefficient of friction, μk Tolerance parameter, k Baumgarte parameter, α Baumgarte parameter, β Integration algorithm Reported time step Integration tolerance Simulation Time

0.1 10 0 0 s/m 5 5 ode45 10−5 s 10−7 1s

tested in order to assess the influence of an initial misalignment of the revolute clearance joint in the dynamic response of the multibody system. In the first case, the local coordinate system of each body is aligned with the global coordinate system, i.e., the joint is completely aligned in the nominal position. In the second scenario, a rotation of the connecting rod of 0.2° over its own axis is considered. This rotation is possible due to the existence of an ideal spherical joint between the crank and the connecting rod. In terms of the revolute clearance joint, the journal is located at the connecting rod, while the bearing belongs to the slider. The dimensions of these elements are listed in Table 2, together with the contact properties to be used in the clearance joint formulation. For this example, the same clearance size is considered for both axial and radial directions, and the contact stiffness is identical for all contacts. Besides the contact forces, only the gravitational and inertial forces act on the system during the entire simulation. The numerical results obtained from the computational simulations of the two previously described scenarios are discussed and compared with the ideal joint case, in which the clearance joint is replaced by a kinematic spherical joint to avoid any redundant constraints. The total simulation time is 1 s. However, in order to facilitate the analysis, the outputs only report a portion of the total simulations performed, which represents a steady-state behavior following the initial impacts. The dynamic response of the mechanism is analyzed through the time plots of the position, velocity and acceleration of the slider, as depicted in Figs. 7 and 8 for each of the referred initial scenarios. Thus, the modeling of a clearance joint has a significant influence on the dynamic response of the system, particularly at the acceleration level. Regarding the position and velocity, the plots in Figs. 7 and 8 show the identical curve for the ideal and clearance joint with a slight lag due to the energy dissipation during the contact process. In terms of the different initial conditions, the macroscopic shapes of the position and velocity curves are similar. The acceleration plots show some differences, namely the aligned case has more peaks with smaller magnitudes (Fig. 7c), while the misaligned scenario exhibits fewer peaks with higher magnitudes (Fig. 8c). The phase portraits are used here to evaluate the dynamic response of the system. The slider velocity versus position and the slider acceleration versus velocity are the variables considered in the portraits to study the level of nonlinearity, as reported in Figs. 9 and 10. The modeling of a revolute clearance joint introduces nonlinear motion into the behavior of the mechanism. This phenomenon is undoubtedly related to the collisions associated with the clearance joint, which provoke sudden changes in the kinematics of the slider, which can be seen mainly in the velocity versus acceleration plots represented in Figs. 9(b) and 10(b). The effects of initial misalignment in a mechanism with planar motion are only felt during a preliminary phase when the type of intra-joint contacts is distinct, namely the contacts presented in Fig. 2(a), (c), (e), (f) and (g) have occurred. However, both scenarios present similar dynamic response for the rest of the simulation. Fig. 11 presents the evolution of the third Euler parameter of the connecting rod, which represents the rotation over its own axis. In the first case, when the joint is initially aligned, this Euler parameter is always null, since the mechanism has planar motion and no perturbations exist. This means that the mechanism performs a planar motion, and the joint does not have axial movement. In the second case (Fig. 11(a)), the dynamic analysis is initiated with a rotation of the connecting rod over η2 , and at the beginning of movement there are some oscillations in that direction. However, due to the initial impacts, the bearing accommodates the movement of the journal, which reduces the oscillations, and the mechanism tends to follow a planar motion. The same trend is shown by the evolution of the angle between the axes of the bearing and the journal, which is represented in Fig. 11(b). 6.2. Slider-crank mechanism describing spatial motion In this section, the study of the effectiveness of the revolute clearance joint formulation is extended to the dynamic analysis of the spatial slider-crank mechanism represented in Fig. 12. This mechanical system consists of three moving bodies, that is, the crank, the connecting rod and the slider, and their geometrical and inertial properties are listed in Table 3. The initial configuration of the mechanism exhibits the crank in a vertical position, and with an initial angular velocity of 2π rad/s in the positive direction of the y-axis. The remaining coordinates and velocities are calculated in such a way that

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

133

(a)

(b)

(c) Fig. 7. Kinematics of the slider for initial conditions with the aligned radial clearance joint: (a) position, (b) velocity, (c) acceleration.

Table 3 Geometrical and inertia properties of each body. Body

Crank Connecting Rod Slider

Length [m]

0.10 0.29 –

Mass [kg]

0.12 0.50 0.50

Principal Moments of Inertia [kg m2 ] Iξ ξ

Iηη

Iζ ζ

0.0 0 01 0.0040 0.0 0 01

0.0 0 01 0.0 0 04 0.0 0 01

0.0 0 0 01 0.0040 0.0 0 01

134

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

(a)

(b)

(c) Fig. 8. Kinematics of the slider for initial conditions with the misaligned radial clearance joint: (a) position, (b) velocity, (c) acceleration.

kinematic constraints of an ideal mechanism are fulfilled. Gravitational forces are considered in the negative direction of the z-axis. In this example, the revolute clearance joint is located at the connection between the crank and the ground. The bearing is placed in the crank, while the journal belongs to the ground. The remaining joints of this mechanical system are considered to be ideal, namely a spherical (crank-connecting rod), a universal (connecting rod-slider) and a translational joint

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

135

(a)

(b) Fig. 9. Phase portraits of (a) the slider velocity versus position and (b) the slider acceleration versus velocity for initial conditions with the aligned clearance joint.

(slider-ground). The properties that characterize the clearance joint of this multibody model are defined in Table 2 and are used here again. Furthermore, the multibody dynamics analysis was also performed for 1 s with the mechanism in operation. The effect of modeling a clearance joint is assessed through the analysis of the slider motion since it is the kinematic output of this mechanism. Thus, the slider position, velocity and acceleration are compared considering an ideal joint and a joint with clearance and friction, as depicted in Fig. 13(a)–(c). The existence of clearance involves the occurrence of several collisions between the journal and bearing, which highly affects the motion of the mechanism. These subsequent impacts result in peaks of the acceleration, as reported in Fig. 13(c). Fig. 13(b) depicts a staircase curve for the velocity of the slider, which is the consequence of high accelerations for short time intervals. The contact-impact phenomenon is also a cause for energy loss, as shown by the variation of mechanical energy (Fig. 13(d)) and by the reduction of amplitude in the position of the slider (Fig. 13(a)). The analysis and comparison of Figs. 13(c) and (d) show the larger losses of energy occur during the impacts with higher peaks of acceleration. Moreover, the phase portraits of velocity versus position and acceleration versus velocity demonstrate that, when a clearance joint is modeled, the system motion is highly nonlinear in contrast to the case of an ideal joint. This characteristic is more evident in the plots of Fig. 13(f) as shown by the dense overlapping lines. Fig. 14 shows the trajectories of the centers of the bases of the journal with respect to the bearing walls, only data on the first 0.20 s is displayed. Initially, the journal and the bearing are aligned. Then, several impacts and continuous contact between the journal and bearing take place. The analysis is now extended to the parameter study of the clearance size, either in radial or axial directions. First, the axial clearance size was set to 0.1 mm, and the dynamic simulation of the system was performed for four different values of radial clearance, namely 0.05, 0.1, 0.3 and 0.5 mm. An identical procedure was adopted to study the effect of the clearance size on the axial direction. Figs. 15 and 16 show the time evolution of the acceleration of the slider for different values of radial and axial clearance, respectively. The analysis of these figures demonstrates that the clearance size in either

136

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

(a)

(b) Fig. 10. Phase portraits of (a) the slider velocity versus position and (b) the slider acceleration versus velocity for initial conditions with the misaligned clearance joint.

direction has a similar influence, i.e., the increase of clearance is typically associated with the increase of magnitude of the acceleration peaks. Acceleration is directly related to the contact forces, therefore higher clearances produce higher impact forces in the joints. The influence of the clearance size can also be examined by the phase portraits of velocity versus acceleration presented in Figs. 17 and 18, respectively, for both radial and axial directions. These plots show that the motion of the mechanism tends to be more chaotic with the increase of clearance. Again, identical effects are observed for both radial and axial clearances. The performance of the spatial slider-crank mechanism was also analyzed by combing all of the radial and axial clearances cases considered above. Fig. 19 displays a 3D plot of the maximum acceleration of the slider as a function of the two clearance values. This analysis shows that the increase of clearance provokes higher accelerations, and this effect has similar significance for both types of clearance. These conclusions are valid for this particular mechanism, where the joint and radial and axial clearances have similar sizes. 7. Concluding remarks A comprehensive methodology for spatial dynamics modeling of three-dimensional revolute joints with clearance was presented in this paper. The formulation proposed here has been developed under the framework of multibody system methodologies and accounts for the combined effects associated with the radial and axial displacements, which is one of the innovative aspects of this work. The characterization of a three-dimensional revolute clearance joint is based on Cartesian coordinates, and the joint elements are modeled as two rigid bodies that can collide with each other. The methodology developed here takes into consideration all the possible contact scenarios between the journal and bearing parts that compose this type of joint. Moreover, the kinematics of the joint elements are fully described, being the basis for the evaluation of the intra-joint contact forces generated during the contact-impact process. Normal and tangential contact forces that incorporate the geometrical and material properties of the contacting bodies govern the dynamics of three-dimensional rev-

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

(a)

137

(b)

Fig. 11. (a) Evolution of the third Euler parameter (e2 ) of the connecting rod, and (b) angle between the journal and the bearing in the misaligned case.

z

x

y

Clearance joint

Fig. 12. Multibody model for the spatial slider-crank mechanism.

olute clearance joints. The normal contact force is determined as a function of the elastic local penetration together with a dissipative term that represents the energy dissipation during the contact-impact events. Furthermore, a modified Coulomb’s friction law is used to model the friction phenomena, which is characterized by good numerical stability in terms of numerical resolution of the equations of motion. Two different multibody systems were utilized to assess the effectiveness of the methodology proposed in this investigation. First, a slider-crank mechanism describing planar motion was used for two distinct scenarios. In the first scenario, the journal and bearing were positioned in the nominal position, which led to perfectly aligned contacts. This particular example of application showed that the modeling of clearance significantly affects the mechanism dynamic response. The dynamic analysis of this mechanism in the second scenario with a misalignment between the journal and bearing elements demonstrated different contact types, which represents a real event more clearly. A three-dimensional slider-crank mechanism was also used to study the applicability of the methodology. The respective numerical simulations again confirmed the paramount importance of taking into account the clearance in the dynamic modeling of a mechanical system. Furthermore, a parameter study on the clearance size was carried out for both radial and axial clearances. The numerical results showed that both types of clearance could have identical effects on the motion of the mechanism, that is, an increase in the clearances increases the nonlinearity of the system’s behavior.

138

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 13. Comparison of slider motion for ideal and revolute clearance joint (a) position, (b) velocity, (c) acceleration, (d) variation of mechanical energy, and phase portraits of (e) velocity vs position and (f) acceleration vs velocity.

ca

ca

orthogonal view

isometric view

orthogonal view

Fig. 14. Trajectory of the centers of the journal bases with respect to the bearing walls.

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

(a)

(b)

(c)

(d) Fig. 15. Effect of radial clearance size on the acceleration of the slider: (a) cr = 0.05 mm, (b) cr = 0.1 mm, (c) cr = 0.3 mm, (d) cr = 0.5 mm.

139

140

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

(a)

(b)

(c)

(d) Fig. 16. Effect of axial clearance size on the acceleration of the slider: (a) ca = 0.05 mm, (b) ca = 0.1 mm, (c) ca = 0.3 mm, (d) ca = 0.5 mm.

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

141

(a)

(b)

(c)

(d) Fig. 17. Effect of radial clearance size on the phase portraits of slider acceleration vs. velocity: (a) cr = 0.05 mm, (b) cr = 0.1 mm, (c) cr = 0.3 mm, (d) cr = 0.5 mm.

142

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

(a)

(b)

(c)

(d) Fig. 18. Effect of axial clearance size on the phase portraits of slider acceleration vs. velocity: (a) ca = 0.05 mm, (b) ca = 0.1 mm, (c) ca = 0.3 mm, (d) ca = 0.5 mm.

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

143

Fig. 19. Effect of radial and axial clearance on the peak acceleration of the slider.

Acknowledgments The first and second authors express their gratitude to the Portuguese Foundation for Science and Technology through the PhD grants (PD/BD/114154/2016 and PD/BD/128385/2017). This work has been supported by the Portuguese Foundation for Science and Technology with the reference project UID/EEA/04436/2013, by FEDER funds through the COMPETE 2020 - Programa Operacional Competitividade e Internacionalização (POCI) with the reference project POCI-01-0145-FEDER-006941. Finally, the authors are much indebted to the anonymous reviewers for useful comments, recommendations and suggestions. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

T.P. Goodman, Dynamic effects of backlash, Mach. Des. 35 (1963) 150–157. R.S. Haines, Survey: 2-dimensional motion and impact at revolute joints, Mech. Mach. Theory 15 (1980) 361–370. Y. Liu, Y. Yu, A survey of mechanism and robot with clearances, Mech. Sci. Technol. 23 (4) (2004) 454–457. O. Muvengei, J. Kihiu, B. Ikua, Dynamic analysis of multi-body mechanical systems with imperfect kinematic joints: a literature survey and review, Sustainable Res. Innov. Proc. 3 (2011) 16. A. Gummer, B. Sauer, Modeling planar slider-crank mechanisms with clearance joints in RecurDyn, Multibody System Dynamics 31 (2014) 127–145. G. Wang, H. Liu, Research progress of joint effects model in multibody system dynamics, Lixue Xuebao/Chinese J. Theor. Appl. Mech. 47 (1) (2015) 31–50 (2015). P. Flores, J. Ambrósio, J.C.P. Claro, H.M. Lankarani, C.S. Koshy, A study on dynamics of mechanical systems including joints with clearance and lubrication, Mech. Mach. Theory 41 (3) (2006) 247–261. Z. Wang, Q. Tian, H. Hu, P. Flores, Nonlinear dynamics and chaotic control of a flexible multibody system with uncertain joint clearance, Nonlin. Dyn. 86 (3) (2016) 1571–1597. P. Flores, H.M. Lankarani, Dynamic response of multibody systems with multiple clearance joints, J. Comput. Nonlin. Dyn. 7 (2012) 031003–031013. L.X. Xu, Y.G. Li, Modeling of a deep-groove ball bearing with waviness defects in planar multibody system, Multibody Syst. Dyn. 33 (3) (2015) 229–258. Z.F. Bai, Y.Q. Liu, Y. Sun, Investigation on dynamic responses of dual-axis positioning mechanism for satellite antenna considering joint clearance, J. Mech. Sci. Technol. 29 (2015) 453–460. J. Ma, L. Qian, G. Chen, M. Li, Dynamic analysis of mechanical systems with planar revolute joints with clearance, Mech. Mach. Theory 94 (2015) 148–164. M.A.B. Abdallah, I. Khemili, N. Aifaoui, Numerical investigation of a flexible slider-crank mechanism with multijoints with clearance, Multibody Syst. Dyn. 38 (2) (2016) 173–199. Y. Li, G. Chen, D. Sun, Y. Gao, K. Wang, Dynamic analysis and optimization design of a planar slider-crank mechanism with flexible components and two clearance joints, Mech. Mach. Theory 99 (2016) 37–57. J. Costa, J. Peixoto, P. Moreira, A.P. Souto, P. Flores, H.M. Lankarani, Influence of the hip joint modeling approaches on the kinematics of human gait, ASME J. Tribology 138 (3) (2016) 031201–031210. F. Farahanchi, S. Shaw, Chaotic and periodic dynamics of a slider-crank mechanism with slider clearance, J. Sound Vibration 177 (3) (1994) 307–324. P. Flores, J. Ambrósio, J.C.P. Claro, H.M. Lankarani, Translational joints with clearance in rigid multibody systems, J. Comput. Nonlin. Dyn. 3 (1) (2008) 0110071. P. Flores, R. Leine, C. Glocker, Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach, Multibody Syst. Dyn. 23 (2010) 165–190. Q.Z. Chen, V Acary, G. Virlez, O. Brüls, A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints, Int. J. Numerical Methods Eng. 96 (2013) 487–511. J. Zhang, Q. Wang, Modeling and simulation of a frictional translational joint with a flexible slider and clearance, Multibody Syst. Dyn. 38 (4) (2016) 367–389. P. Flores, J. Ambrósio, J.C.P. Claro, H.M. Lankarani, Dynamics of multibody systems with spherical clearance joints, J. Comput. Nonlin. Dyn. 1 (3) (2006) 240–247. C. Quental, J. Folgado, J. Ambrósio, J. Monteiro, Multibody system of the upper limb including a reverse shoulder prosthesis, J. Biomech. Eng. 135 (11) (2013) 111005.

144

F. Marques et al. / Mechanism and Machine Theory 116 (2017) 123–144

[23] E. Askari, P. Flores, D. Dabirrahmani, R. Appleyard, Nonlinear vibration and dynamics of ceramic on ceramic artificial hip joints: a spatial multibody modeling, Nonlin. Dyn. 76 (2014) 1365–1377. [24] Z. Jing, G. Hong-Wei, L. Rong-Qiang, D. Zong-Quan, Nonlinear characteristic of spherical joints with clearance, J. Aerosp. Technol. Manage. 7 (2) (2015) 179–184. [25] G. Wang, H. Liu, Dynamics analysis of 4-SPS/CU parallel mechanism with spherical joint clearance, J. Mech. Eng. 51 (1) (2015) 43–51. [26] Q. Tian, Y. Zhang, L. Chen, P. Flores, Dynamics of spatial flexible multibody systems with clearance and lubricated spherical joints, Comput. Struct. 87 (13-14) (2009) 913–929. [27] P. Flores, J. Ambrósio, J.C.P. Claro, H.M. Lankarani, Spatial revolute joints with clearance for dynamic analysis of multibody systems, Proc. Inst. Mech. Eng., Part-K J. Multi-body Dyn. 220 (4) (2006) 257–271. [28] Q. Tian, C. Liu, M. Machado, P. Flores, A new model for dry and lubricated cylindrical joints with clearance in spatial flexible multibody systems, Nonlin. Dyn. 64 (1-2) (2011) 25–47. [29] C. Brutti, G. Coglitore, P.P. Valentini, Modeling 3D revolute joint with clearance and contact stiffness, Nonlin. Dyn. 66 (2011) 531–548. [30] C. Liu, Q. Tian, H. Hu, Dynamics and control of a spatial rigid-flexible multibody system with multiple cylindrical clearance joints, Mech. Mach. Theory 52 (2012) 106–129. [31] J. Ambrósio, P. Verissimo, Improved bushing models for general multibody systems and vehicle dynamics, Multibody Syst. Dyn. 22 (4) (2009) 341–365. [32] A. Chaker, A. Mlika, M.A. Laribi, L. Romdhane, S. Zeghloul, Clearance and manufacturing errors’ effects on the accuracy of the 3-RCC spherical parallel manipulator, European J. Mech. - A/Solids 37 (2013) 86–95. [33] Z. Jing, G. Hong-Wei, L. Rong-Qiang, D. Zong-Quan, Nonlinear characteristic of spherical joints with clearance, J. Aerosp. Technol. Manage. 7 (2) (2015) 179–184. [34] C. Pereira, J. Ambrósio, A. Ramalho, Dynamics of chain drives using a generalized revolute clearance joint formulation, Mech. Mach. Theory 92 (2015) 64–85. [35] F. Marques, F. Isaac, N. Dourado, P. Flores, 3D formulation for revolute clearance joints, in: New Trends in Mechanism and Machine Science. Vol. 43 of the series Mechanisms and Machine Science, 2017, pp. 183–191. [36] P. Flores, J. Ambrósio, On the contact detection for contact-impact analysis in multibody systems, Multibody Syst. Dyn. 24 (1) (2010) 103–122. [37] S.G. Dhande, J. Chakraborty, Mechanical error analysis of spatial linkages, J. Mech. Des. 100 (4) (1978) 732–738. [38] J.H. Rhyu, B.M. Kwak, Optimal stochastic design of four-bar mechanisms for tolerance and clearance, J. Mech., Transm. Automation Des. 110 (3) (1988) 255–262. [39] S. Dubowsky, J.F. Deck, H. Costello, The dynamic modeling of flexible spatial machine systems with clearance connections, J. Mech., Transm. Automation Des. 109 (1) (1987) 87–94. [40] T. Kakizaki, J.F. Deck, S. Dubowsky, Modeling the spatial dynamics of robotic manipulators with flexible links and joint clearances, J. Mech. Des. 115 (4) (1993) 839–847. [41] J.F. Deck, S. Dubowsky, On the limitations of predictions of the dynamic response of machines with clearance connections, J. Mech. Des. 116 (3) (1994) 833–841. [42] C. Brutti, G. Coglitore, P.P. Valentini, Modeling 3D revolute joint with clearance and contact stiffness, Nonlinear Dynamics 66 (4) (2011) 531–548. [43] O. Pinkus, S.A. Sternlicht, Theory of Hydrodynamic Lubrication, McGraw Hill, New York, 1961. [44] S. Yan, W. Xiang, L. Zhang, A comprehensive model for 3D revolute joints with clearances in mechanical systems, Nonlin. Dyn. 80 (1) (2015) 309–328. [45] C. Liu, K. Zhang, L. Yang, The compliance contact model of cylindrical joints with clearances, Acta Mechanica Sinica 21 (2005) 451–458. [46] C.S. Liu, K. Zhang, R. Yang, The FEM analysis and approximate model for cylindrical joints with clearances, Mech. Mach. Theory 42 (2007) 183–197. [47] O.A. Bauchau, J. Rodriguez, Modeling of joints with clearance in flexible multi-body systems, Int. J. Solids Struct. 39 (1) (2002) 41–63. [48] O.A. Bauchau, C. Ju, Modeling friction phenomena in flexible multi-body dynamics, Comput. Methods Appl. Mech. Eng. 195 (50) (2006) 6909–6924. [49] V. Parenti-Castelli, S. Venanzi, Clearance influence analysis on mechanisms, Mech. Mach. Theory 40 (12) (2005) 1316–1329. [50] S. Venanzi, V. Parenti-Castelli, A new technique for clearance influence analysis in spatial mechanisms, J. Mech. Des. 127 (3) (2005) 446–455. [51] S. Erkaya, S. Dog˘ an, E. S¸ efkatlıog˘ lu, Analysis of the joint clearance effects on a compliant spatial mechanism, Mech. Mach. Theory 104 (2016) 255–273. [52] Q. Tian, J. Lou, A. Mikkola, A new elastohydrodynamic lubricated spherical joint model for rigid-flexible multibody dynamics, Mech. Mach. Theory 107 (2017) 210–228. [53] F. Isaac, F. Marques, N. Dourado, P. Flores, Recent developments on cylindrical contact force models with realistic properties, in: New Trends in Mechanism and Machine Science. Vol. 43 of the series Mechanisms and Machine Science, 2017, pp. 211–219. [54] F. Marques, Frictional contacts in multibody dynamics, MSc Dissertation, University of Minho, Guimarães, Portugal, 2015. [55] P.E. Nikravesh, Computer Aided Analysis of Mechanical Systems, Prentice Hall, Englewood Cliffs, New Jersey, 1988. [56] J. Alves, N. Peixinho, M.T. Silva, P. Flores, H. Lankarani, A comparative study on the viscoelastic constitutive laws for frictionless contact interfaces in multibody dynamics, Mech. Mach. Theory 85 (2015) 172–188. [57] H. Hertz, Über die Berührung fester elastischer Körper, J. reine und an-gewandte Math. 92 (1881) 156–171. [58] K.H. Hunt, F.R.E. Crossley, Coefficient of restitution interpreted as damping in vibroimpact, J. Appl. Mech. 42 (1975) 440–445. [59] H.M. Lankarani, P.E. Nikravesh, A contact force model with hysteresis damping for impact analysis of multibody systems, J. Mech. Des. 112 (3) (1990) 369–376. [60] M. Machado, P. Moreira, P. Flores, H.M. Lankarani, Compliant contact force models in multibody dynamics: evolution of the Hertz contact theory, Mech. Mach. Theory 53 (2012) 99–121. [61] A. Banerjee, A. Chanda, R. Das, Historical origin and recent development on normal directional impact models for rigid body contact simulation: a critical review, Arch. Comput. Methods Eng. 24 (2) (2017) 397–422. [62] F. Marques, P. Flores, J.C.P. Claro, H.M. Lankarani, A survey and comparison of several friction force models for dynamic analysis of multibody mechanical systems, Nonlin. Dyn. 86 (3) (2016) 1407–1443. [63] C.A. Coulomb, Théorie Des Machines simples, En Ayant égard au Frottement De Leurs parties, Et à La Roideur Des Cordages, Mémoire de Mathématique et de Physique, Paris, France, 1785. [64] E. Pennestrì, V. Rossi, P. Salvini, P.P. Valentini, Review and comparison of dry friction force models, Nonlin. Dyn. 83 (4) (2016) 1785–1801. [65] F. Pichler, W. Witteveen, P. Fischer, A complete strategy for efficient and accurate multibody dynamics of flexible structures with large lap joints considering contact and friction, Multibody Syst. Dyn. (2016), doi:10.1007/s11044- 016- 9555- 2. [66] J. García de Jálon, E. Bayo, Kinematic and Dynamic Simulations of Multibody Systems, Springer Verlag, New York, 1994. [67] J. Baumgarte, Stabilization of constraints and integrals of motion in dynamical systems, Comput Methods Appl. Mech. Eng. 1 (1972) 1–16. [68] P. Flores, M. Machado, E. Seabra, M.T. Silva, A parametric study on the Baumgarte stabilization method for forward dynamics of constrained multibody systems, J. Comput. Nonlin. Dyn. 6 (1) (2011) 0110191–0110199. [69] F. Marques, A.P. Souto, P. Flores, On the constraints violation in forward dynamics of multibody systems, Multibody Syst. Dyn. 39 (4) (2017) 385–419.