Dynamic response analysis for multi-degrees-of-freedom parallel mechanisms with various types of three-dimensional clearance joints

For spatial multibody systems, the dynamic equations of multibody systems with compound clearance joints have a high level of nonlinearity. The coupling between different types of clearance joints may lead to abundant dynamic behavior. At present, the dynamic response analysis of the spatial parallel mechanism considering the three-dimensional (3D) compound clearance joint has not been reported. This work proposes a modeling method to investigate the influence of the 3D compound clearance joint on the dynamics characteristics of the spatial parallel mechanism. For this purpose, 3D kinematic models of spherical clearance joint and revolute joint with radial and axial clearances are derived. Contact force is described as normal contact and tangential friction and later introduced into the nonlinear dynamics model, which is established by the Lagrange multiplier technique and Jacobian of constraint matrix. The influences of compound clearance joint and initial misalignment of bearing axes on the system are analyzed. Furthermore, validation of dynamics model is evaluated by ADAMS and Newton–Euler method. This work provides an essential theoretical basis for studying the influences of 3D clearance joints on dynamic responses and nonlinear behavior of parallel mechanisms.


Introduction
Robots are becoming an indispensable part of modern industry and are playing an increasingly important role. 1,2 For some reasons, such as manufacturing tolerance, assembling errors, and wear phenomenon of the kinematic joints, clearance always exists in kinematics joints, and the performance of the robot is unsatisfactory. [3][4][5] Planarization of spatial clearance joints has limitations in studying the dynamics behavior for mechanical systems because the complex motion mechanism of spatial clearance joints is neglected. Moreover, frequent contact and collision events are the source of nonlinearity. It may exhibit chaotic responses under certain conditions and even affect the overall dynamic behaviors of the system.
Over the past decades, issues related to clearance joint in mechanism dynamics have been a hot topic. 6 It is noteworthy that most of the research were mainly associated with the issue of plane system. However, dynamics research of multi-degrees-of-freedom (DOF) spatial systems considering 3D compound clearance joint has not been reported as far as we know. Yu et al. 7 eliminated redundant constraints of an overconstrained mechanism and further investigated the dynamics model and performances of the system considering multiple revolute clearance joints. Zhan et al. 8 modeled joint clearance parameters as interval variables and provided fresh ideas for studying motion reliability of 3-RRR parallel manipulators considering joint clearances. Li et al. 9 built a rigid-flexible coupling dynamic model to analyze the dynamics behavior of solar array. A rigid frame and flexible parts are described by nodal coordinate formulation and absolute nodal coordinate formulation frames, respectively. Zheng et al. 10 developed a modeling technique based on nonsmooth dynamic method and have studied the effect of stiction on the dynamic response. Zhao et al. 11 presented a modeling approach of a piston skirt-liner system with lubricated pairs and evaluated lubrication and dynamic performances of the mechanism. It is worth mentioning that the above work takes into account different factors and provides ideas for establishing an accurate dynamic model with clearance joints. Different from the previous method, Ordiz et al. 12 used Simscape to model the industrial press. Then, ANSYS was used for stress analysis to explore the impact of clearance evolution on service life of the machines. Wu et al. 13 studied the nonlinear characteristics of a crank-slider mechanism with multiple clearance joints through correlation dimensions. Ambrósio and Pombo 14 and Marques et al. 15 improved pre-existing dynamic modeling methods of the mechanisms containing clearance joint, respectively. Both of these methods simplify the modeling process and improve modeling efficiency.
The above literature works only for planar systems. However, it is necessary to consider 3D models of the system in some cases. Erkaya 16 investigated the influence of spherical clearance joints and flexible joint connection on the vibration of linkage by experiment and numerical simulation. Chen and Sun 17 established a dynamics model of parallel mechanism containing spherical clearance pair. The clearance model in the above article has been extended to three dimensions, and the established dynamic model is more accurate than the plane model. However, the effect of revolute clearance pair is ignored. Using nonsmooth modeling method, Cavalieri and Cardona 18 developed a 3D model of revolute pair, in which axial clearance is neglected. Yan et al. 19 validated the existence of 3D revolute pair with clearance through experiments. They verified correctness of the proposed model, which can describe the axial and radial clearance. Tian et al. 20 investigated the influences of lubrication cylindrical joint and flexibility component on dynamics responses of crank-slider system. Then, Liu et al. 21 developed a four-point model to detect the contact between cylinder elements. Flores et al. 22 deducted a kinematic model for revolute clearance joints.
They carried out numerical simulations on the double pendulum mechanism. A new enhanced formulation was recently developed by Marques et al. 23 to simulate axial and radial contact of the revolute clearance joint. The simulation results show that it is necessary to consider radial and axial clearances in dynamic modeling process. Subsequently, Marques et al. 24 considered different friction models and explored the dynamic response of the RSSR mechanism. As can be seen from the above literature, many scholars have developed different 3D clearance joint frames for application in different fields and have achieved good results.
For multibody systems, to reveal dynamics behavior of the system more accurately and truthfully, various types of clearance joints should be considered in the dynamic model. In this article, a motion frame of 3D revolute clearance joint is proposed. A dynamics modeling method of parallel mechanism considering 3D compound clearance joint is given. The influence of the compound clearance joint on the dynamic response of the spatial parallel mechanism is discussed. Furthermore, evaluating influences of clearance joint on output responses are all prerequisites for analysis, design, optimization, and control of nonlinear systems. The purposes of the article are to analyze the influences of 3D compound clearance joint on dynamics responses of parallel mechanisms.
The structure of the remainder of this article is as follows: The second section provides a mathematical description of the 3D revolute clearance joint and 3D spherical clearance joint, including kinematics model and contact force model. Nonlinear dynamics model of the parallel mechanism with 3D revolute clearance joint and spherical clearance joint is presented in the third section. A computational procedure for the dynamic model is also presented in this section. In the fourth section, influences of compound clearance joint and axis misalignment on system response are explored, and dynamics model is verified by ADAMS and Newton-Euler method, respectively. Salient conclusions are presented in the fifth section.

Model for 3D clearances joints of parallel mechanism
Mathematical description for 3D spherical clearance joint Introduction of clearance into a spherical pair released three kinematic restrictions, which can be shown in Figure 1. It is generally considered that only one contact form, point contact, exists in dynamic simulation modeling because of its structural reasons. The geometric centers of the two elements (socket and ball) of the spherical joint are expressed as s i and b j . Vectors r si and r bj are positions of s i and b j in global coordinate frame. Two potential contact points are denoted as P si and P bj . Radii of ball and socket can be written as R bj and R si .
The clearance of the spherical joint can be described as The vector of eccentricity e is defined by Relative penetration depth is written as

Mathematical description for 3D revolute clearance joint
For idea revolute pair, the axis of journal and bearing always keep in coincidence during joint movement. However, when the robot is connected by the joint, clearance inevitably exists. In that case, axis misalignment may occur in the clearance joint. This results in different joint element configurations, namely, no contact, point contact, line contact, and plane contact. The above contact conditions can be considered as a combination of independent axial and radial contact conditions. The schematic diagram of four contact configurations can be found in other scholars' work. 19,22,23 A 3D revolute clearance joint is shown in Figure 2, and it consists of bearing and journal. The radius of the journal is R jo . R b1 and R b2 are the inner and outer radii of bearing, respectively. L j and L b denote the length of journal and bearing, respectively. The axial and radial clearances can be defined as For dealing with point contact, the contact impact detection of 3D revolute clearance joint can be regarded as the contact impact detection problem of two impact rings located at the top and bottom of bearing. Therefore, as shown in Figure 3, four potential contact points are located in the top impact ring, where P b and P j are used to detect contact between bearing top surface and journal flange, M b and M j are used to detect the impact between bearing inner wall and journal. Similarly, the four potential contact points corresponding to the bottom impact ring are Q b and Q j , N b and N j , respectively. Centers of journal and bearing are o j and o b , respectively. J T and B T are the centers of the top surface of journal and bearing, respectively. Unit vectors e j and e b in global coordinate frame are aligned with the axis of journal and bearing.
Position vectors of J T and B T can be given by where r oj and r ob are position vectors of o j and o b in global coordinate frame. The angle between journal flange and bearing is z, and the angle between vector e j and bearing is x. Construct a plane consisting of three points J T , B T , and J B , the unit vector of normal direction of the plane is Normal contact direction is derived as n M ¼ n 1 Âe j . The projection of vector e b on vector n 1 is Z 1 point, and its position vector can be given by The location vector of contact point M b of bearing can be written as Relative penetration depth d M caused by radial impact is evaluated as The position of potential contact point M j can be obtained by combining equations (7) and (8) For the opposite sides of contact state, the analysis process is similar. Figure 4 shows the schematic of axial contact state. It should be noted that potential contact points of axial contact and radial contact are not necessarily in the same plane. However, the position of potential axial contact is always related to e b . Moving vector e j to the starting point of vector e b , construct a plane that contains both vectors e j and e b . Normal vector of the plane is expressed by Normal contact direction of points P b and P j can be written as n p ¼ À e j .
The potential contact point at bearing top surface and the corresponding potential contact point on the journal flange are P b and P j , respectively. Their corresponding position vectors can be written as where d p , caused by contact between bearing top surface and journal flange, is the relative penetration depth and can be given by The analysis of contact states occurring in the bottom surface of bearing and journal flange is similar to the analysis above.
The above analysis is about point contact cases. When bearing and journal are parallel to each other, z ¼ x ¼ 0, the line contact state may occur in the radial direction, and the plane contact state may occur in the axial direction. This is a particular case of the above deduction.

Calculation of contact force for clearances joints
With the motion of clearance joint, contact-impact events frequently occur at clearance joints. An appropriate contact  force model is a vital part of restoring the dynamic behavior of kinematic pair and predicting global dynamics behaviors of the mechanism.
L-N continuous contact force model based on Hertz contact law is widely employed to describe normal contact force of contact-impact events. Dissipation term is introduced to reflect the energy dissipation in the contact process. 9,19,21 In this article, L-N model 25 is employed to describe normal contact force of clearance joints in point contact state where the first term Kd n k is estimated by Hertz contact model and used to describe elastic forces, the second term Hd n k _ d k is used to assess energy dissipation. The value of n is 3/2. d k represents the relative penetration depth at each contact point. _ d k is relative normal contact velocity of two contact bodies.
For the radial contact between two cylinders in a parallel state, based on Hertz pressure distribution, nonlinear models such as Johnson model, Radzimovsky model, and Goldsmith model have been developed. 26 These models are implicit functions that require numerical iteration to evaluate contact forces without explaining energy dissipation. It is noteworthy that point contact and line contact on the revolute joint does not have the same contact stiffness as in a spherical clearance joint in a real scenario, since they have different topologies. To simplify the calculation, this work adopts an approximate strategy. Hence, equation (12) is used to describe the relationship between force and penetration of line contact state. 23,27 The exponent n ranges from 1 to 3/2. A more general relation, equation (13), is used to describe normal contact force F N A under plane contact state where E jo and L b are Young's modulus and length of bear- Friction is a commonplace event in which the direction is relative to the direction of motion of two objects. Many friction models were developed to describe the friction phenomena reasonably. Coulomb friction model, one of the most classical models, is expressed as a function of normal force and relative tangential velocity. However, some drawbacks also limit the application of the model, such as discontinuity around zero velocity, neglecting Stribeck effect and viscous characteristics. After that, some bristlebased friction models have been reported. [28][29][30][31] These models can describe some complex phenomena such as microsliding but usually lead to inefficiency and introduce additional parameters that are difficult to identify. In this work, a continuous velocity-based friction model 32 is utilized to assess tangential friction force at the contact surface of clearance joint. This model can describe static, dynamic, and viscous friction phenomena without increasing computational time. 33,34 It has continuity and differentiability and is suitable for control and multibody system simulation, which can be written as where v t denotes relative tangential velocity between contact elements, v s is transition velocity, s and k are the static and dynamic friction coefficients. Contact force and moment of potential contact point acting on bearing are given by where t m is the tangential unit vector of potential contact point m and rd m is the radius vector of potential contact points m. Spherical joint only appears point contact state, contact force and moment can be calculated through equations (12), (14), and (15). Detection of the instant of contact is correlated with the penetration depth of two discrete integration steps t n and t nþ1 . It has also been reported in other literature. 9,17,21 Once the collision is detected, contact force should be added to motion equation of multibody system as an external force. Here, note that the normal force does not produce torque. The friction torque acting on the clearance joint should be calculated according to equation (15). Figure 5 shows a 4-UPS/RPS parallel mechanism sketch with 5DOF, including a fixed platform, a moving platform (end effector), and five variable-length kinematic chains. Five kinematic chains can be divided into two categories: UPS (U is universal joint, P represents prismatic joint, and S represents spherical joint) kinematic chains and RPS (R is the revolute joint) kinematic chain. In each kinematic chain, the prismatic joint situated at P i ¼ i ¼ 1; Á Á Á ; 5 ð Þ can be actuated. Five kinematic chains are linked at moving platform through five spherical joints located at

Dynamics model of parallel mechanism with two type clearance joints
The centers of revolute joint and universal joints are named R 1 and U i ¼ i ¼ 1; Á Á Á ; 5 ð Þ , respectively. In this model, the revolute joint R and spherical joint S 1 distributed on the RPS kinematic chain are regarded as imperfect joints. The geometrical and mass properties of the 4-UPS/RPS parallel mechanism are presented in Table 1.
As illustrated in Figure 5, global coordinate frame O A -X A Y A Z A is attached at mass center of the fixed platform. Furthermore, all the vectors mentioned later are represented in this framework. Similarly, the local coordinate frames o i -x i y i z i i ¼ 1; Á Á Á ; 11 ð Þare attached at the mass center of active components, respectively.
The position and orientation of active components are parameterized by variable where x i y i z i ½ T are position vectors of the origins of local coordinate frames in global coordinate frame. a i , b i , and g i are Euler angles.
Rotation matrix from local coordinate frames Universal joint allows relative rotation along perpendicular axes. The kinematic constraint can be written in the form where r j denotes position vectors of origins of local coordinate frames in global coordinate frame. R j is the rotation matrix from local coordinate frames to global coordinate frame. o u i ¼ r 2 0 cos i ð Þsin i ð Þ ½ T and o u j ¼ Àl s =200 ½ T are position vectors of U i in global coordinate frame and local coordinate frame, respectively. Ro i is a unit vector along z j direction in global coordinate frame. Ro j ¼ R j 0 1 0 ½ T is unit vector along y j direction. Prismatic joint, as the actuated joint, allows relative translational along the joint axis. In this mechanism, the constraint equation is defined as where the unit vectors p i and q j are aligned with joint axis of x i and x j . d n is a unit vector from the origin o i to o j , f i and f j are two unit vectors along axes z i and y j , respectively. The kinematic constraint of the perfect revolute joint can be written as Ro r ð Þ T Á Ro j2 2 6 4 3 7 5 ¼ 0 5Â1 (20) where o r denotes position vectors of the revolute joint in global coordinate frame. Ro r is a unit vector and can be written as Ro r ¼ 0 0 1 ½ T , Ro j2 ¼R j 0 1 0 ½ T , and Ro j1 ¼R j 1 0 0 ½ T . Kinematic constraints of the perfect spherical joint can be expressed as where r 11 is the position vector of o 11 in global coordinate frame. R 11 is the rotation matrix from o 11 -x 11 y 11 z 11 to global coordinate frame T are position vectors of S i in local coordinate frames o 11 -x 11 y 11 z 11 and o j -x j y j z j , respectively.
Driving constraint equation of five variable length kinematic chains can be given by where r i denotes the position vectors of origins of local coordinate frames in global coordinate frame and f i (t) is the preset law of motion for each actuated joint. In this work, the revolute joint (R 1 ) and spherical joint (S 1 )are assumed to be clearance joints, and then, they impose no kinematic constraints. Hence, they do not contribute to the overall constraint equation of the mechanism. Combining equations (18) to (22), one obtains an expression for constraint equations of the system as As the numerical iteration progresses, the accumulation of numerical errors will lead to violation of constraint and even nonconvergence of results. Based on the principle of feedback control, Baumgarte's approach 35 controls the numerical error by modifying acceleration constraint equations. Dynamics equation for 4-UPS/RPS parallel mechanism with both 3D revolute clearance joint and 3D spherical clearance joint utilizing Baumgarte stabilization method can be written in matrix form as where M is a 66 Â 66diagonal matrix, representing system mass matrix. Φ q is Jacobian matrix of constraint equations Φ. q includes generalized acceleration vectors of the mechanism. g represents right sides of acceleration constraint equations. l contains Lagrange multipliers associated with constraints. Q represents generalized external forces, such as gravity vector and contact forces. a 1 and b 1 are Baumgarte feedback parameters. The computational methodologies for solving the dynamic equations of the mechanism considering both the 3D revolute clearance joint and 3D spherical clearance joint are shown in Figure 6.

Dynamics response analysis of parallel mechanisms with two types of clearance joints
In this section, the numerical simulation results of the mechanism considering the 3D revolute clearance joint and 3D spherical clearance joint are analyzed. Simulation parameters of clearance joints are presented in Table 2. Initial positions and velocities of each component are deduced by kinematics analysis of this mechanism. The effect of initial axis misalignment of the 3D revolute clearance joint on dynamic response is assessed by comparing the two cases. In the first case, axis of journal and bearing coincides, and in the second case, the axis misalignment occurs at the initial moment, that is, the bearing rotated 0.035 clockwise around coordinate axis x 6 . The motion trajectory of end actuator is defined as (unit: s, m) As shown in Figures 7 and 8, the dynamic behavior at the 3D spherical and revolute clearance joints is given in preparation for the following analysis of global dynamic behavior of the parallel mechanism. Figure 7(a) and (b) shows center trajectory diagram of spherical joint and the corresponding time history diagram. The difference in the initial conditions (axis alignment and axis misalignment) caused the center trajectories not to coincide before 0.515 s. The system gradually stabilizes as it enters the continuous contact state (the trajectories are both concentrated in the rectangular area in Figure 7(a)). From the second cycle, the time history diagram of the two trajectories basically coincides. Besides, in the case of initial alignment, the spherical joint does not move in the Z direction, as shown in Figure 7(b). Figure 7(c) shows center trajectories of bearing relative to journal with initial axis alignment. Combining with Figure 7(d), it can be observed that center trajectory of top and bottom coincides basically. Figure 7(e) and (f) are center trajectory and time history diagrams in the case of initial axis misalignment. Unlike initial axis alignment, there are apparent differences between the top and bottom center trajectories, which continue until 0.515 s, especially in the Y direction. Contact forces of clearance joints are shown in Figure 8(a) and (b). Considering misalignment of axis, contact force of high amplitude appears in the initial stage. When the system enters the stable stage, gravity plays a dominant role. The contact forces of clearance joints are in a continuous contact state. In the stable stage, the collision time of the revolute is comparatively consistent with that of the spherical joint, the two clearance joints are at the upper and lower ends of the same kinematic branch, and the peak contact force of revolute joint is less than that of spherical joint. Driving force curves for prismatic actuators considering initial axis alignment and misalignment are shown in Figure 8(c) and (d). The maximum values of driving force are related to the transient response.
When considering axis misalignment, the peak values of the driving forces are higher than the case of axis alignment. Affected by the contact impact, oscillation positions of driving force curves have been aligned to the contact force curves. The enormous contact force at the initial stage is related to the initial conditions. For multibody dynamics systems with clearance joint, the initial stage response depends on the initial configuration and is transient. It can be observed that the contact force generated by the initial misalignment of the bearing axis is much greater than the initial alignment. The collision phenomenon described based on the penalty method takes extremely short. The forced displacement caused by the driving will cause a larger contact force. After the system enters a stable state, the contact force level returned to normal. Therefore, in dynamic analysis, the response after the initial transient stage is often more concerned. For example, in the literature, 16,[36][37][38] the initial transient response is filtered.
Next, ADAMS was used to verify the dynamic model. The virtual prototype model is imported into ADAMS and then materials are assigned to each component. Fixed platform is connected to the ground by a fixed joint. Four universal joints are added to the fixed platform according to the joint distribution angle. As shown in Figure 5, the ideal spherical joints and prismatic joints are added to the moving platform and kinematic branch chains,    respectively. The joint motions of prismatic joint are selected as translational. Create a cylindrical journal and a cylindrical bearing at the revolute clearance joint and add them to the fixed platform and the RP branch, respectively. The contact type is defined as cylinder to cylinder. Similarly, two balls with different radii are added to the PS branch chain and moving platform. The contact type is defined as sphere to sphere. Normal force and tangential force are defined as impact and Coulomb, and set the same simulation parameters as in the article. From the analysis of Figures 7 and 8, results indicated that initial axis misalignment only greatly impacts the initial stage. Therefore, axis alignment is considered in ADAMS model. From the observations in Figure 9(a) and (b), the compound clearance joint effect on the vibrational levels of displacement curves is inconspicuous. However, as shown in Figure 9(c) to (f), the velocity and acceleration curves of the moving platform fluctuate significantly. These fluctuations occur almost simultaneously with contact at clearance joints.
Compared with numerical simulation results, the curve of ADAMS has little difference in fluctuation degree and amplitude, and the vibration is more severe in Y direction, as displayed in Figure 9(c) and (d). From the acceleration curves (Figure 9(e) and (f)), it can be seen that there are several large fluctuations due to them is an alignment of the axis, which is also supported by contact force curve. The contact impact at the clearance joint is transmitted to the end effector, which can cause a high-frequency response of acceleration, and further to induce velocity fluctuations. It ultimately leads to a decrease in positioning accuracy. In short, clearance will bring a series of adverse effects, such as high-frequency vibration and low positioning accuracy. In this article, the dynamic model of 4-UPS/RPS parallel mechanism is established based on the Lagrange multiplier technology, and the influence of the compound clearance joint on the dynamic response of the system is further analyzed. The constraint equation library of corresponding types of joints can be established through the joint constraint equations given in the article, which can be easily applied to this type of parallel robot. For redundant constraints, LU factorization and SVD factorization techniques can be used to identify and delete redundant constraints in the system.
The dynamic model of the 4-UPS/RPS mechanism established by the constraint joint library using Lagrangian multiplier technology is compared with our previous work (Newton-Euler method). 39 Comparison curve of displacement response is shown in Figure 10. The maximum displacement errors of the two methods in the X and Y directions are 2.28 Â 10 À4 m and 9.5 Â 10 À4 m, respectively, as shown in Figure 10(a) and (b). Coefficient of determination (R 2 ) is used to evaluate agreement observed in the results and R 2 ¼ P y i À y ð Þ= Pŷ i À y ð Þ. y i andŷ i are the i'th data calculated by the Lagrange multiplier technique and Newton-Euler method, respectively, and y is the mean value of the data obtained by the Lagrange multiplier technique. As shown in Figure 10(a) and (b), R 2 of the displacement in the X and Y directions of the model established by the Newton-Euler method and Lagrange multiplier technique is 0.999988 and 0.999952, respectively. A similar method was used to obtain R 2 of the ADAMS axis alignment result and the Lagrange multiplier technique result. As shown in Figure 9(a) and (b), R 2 of the displacement in the X and Y directions are 0.999963 and 0.999927, respectively. The above results show that the agreement between the three methods is excellent. The results of Lagrange multiplier technique and Newton-Euler method are closer. However, the modeling and programming of Newton-Euler method is more complicated. According to the established constraint joint library, it can be directly integrated into the constraint equation of the corresponding mechanism. Although the Lagrangian multiplier introduction increases the computational cost, the modeling has better generalizability and easy to code in a computer program. Compound clearance joint effect on the displacement of end effector is inconspicuous, but it dramatically affects the motion state of the joints. At the same time, the existence of the clearance joint has caused a series of adverse effects such as noise, high-frequency vibration, and wear of the joint. This work has a guiding effect on the structure design and performance optimization of the mechanism.

Conclusions
As an essential nonlinear source, the contact phenomenon of clearance joint brings abundant dynamic behaviors to the system. The influences of the compound clearance joint on the dynamic responses of the system have received little attention in the previous studies.
(1) The kinematics models of 3D spherical joint and revolute joint with axial and radial clearances are established to this end. According to different contact conditions, the corresponding normal contact force model is provided. The friction phenomena between contact elements are described by a continuous velocity-based friction model. (2) The dynamics model of the parallel mechanism is derived using Lagrange multiplier technique and Baumgarte stabilization method. The compound clearance joint and the contact forces are embedded in the dynamics model within the framework of multibody dynamics. mechanisms. This work is of great significance to the structural design and performance optimization of the parallel mechanism. However, the work should further consider the influence of components' flexibility and temperature on the system to expand its application to more industrial robots.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: