Dynamic Accuracy Analysis of a 5PSS/UPU Parallel Mechanism Based on Rigid-Flexible Coupled Modeling

In order to improve the low output accuracy caused by the elastic deformations of the branch chains, a finite element-based dynamic accuracy analysis method for parallel mechanisms is proposed in this paper. First, taking a 5-prismatic-spherical-spherical (PSS)/universal-prismatic-universal (UPU) parallel mechanism as an example, the error model is established by a closed vector chain method, while its influence on the dynamic accuracy of the parallel mechanism is analyzed through numerical simulation. According to the structural and error characteristics of the parallel mechanism, a vector calibration algorithm is proposed to reduce the position and pose errors along the whole motion trajectory. Then, considering the elastic deformation of the rod, the rigid-flexible coupling dynamic equations of each component are established by combining the finite element method with the Lagrange method. The elastodynamic model of the whole machine is obtained based on the constraint condition of each moving part, and the correctness of the model is verified by simulation. Moreover, the effect of component flexibility on the dimensionless root mean square error of the displacement, velocity and acceleration of the moving platform is investigated by using a Newmark method, and the mapping relationship of these dimensionless root mean square errors to dynamic accuracy is further studied. The research work provides a theoretical basis for the design of the parameter size of the prototype.


Introduction
Compared with the serial mechanism, the parallel mechanism has advantages of greater stiffness, higher precision, reduced inertia, and higher payload to weight ratio [1][2][3][4]. It has been widely used in aerospace, food processing, vehicle and ship industry, and medical equipment fields [5][6][7][8].
The geometric errors of mechanical parts are inevitably produced in the process of processing and assembly. The error calibration method can effectively reduce the influence of geometric errors on the accuracy of the parallel mechanism [9,10]. The calibration method is very important for geometric error compensation. The common calibration methods include external calibration, constraint calibration and self-calibration. The pose of the parallel mechanism is directly measured through external precision equipment in the external calibration method [11][12][13], which can identify the geometric parameters of the moving platform. However, it has several disadvantages such as low calibration efficiency, high cost, and poor external anti-interference ability, etc. In the constraint calibration method [14][15][16], the kinematic calibration is performed by imposing constraints on the system to limit the local motion capability of the mechanism components. Compared with the external calibration method, the constraint calibration method is relatively simple to operate and does not require expensive external measuring equipment, but it restricts the mechanical part motion characteristics, which makes it difficult to identify all kinematic parameters. The self-calibration method [11,14,17] only needs the extra redundant information of internal sensors of the mechanism to form full closed-loop control. In view of the economy of hardware and simplification of control system, the self-calibration method is adopted in our study. Meanwhile, the error calibration can improve the motion accuracy of the mechanism, which is of great significance to the dynamic accuracy analysis of the parallel mechanism [18,19].
Dynamic accuracy is an important performance index for a parallel mechanism [20][21][22], which requires that the actual motion state of the mechanism is as consistent as possible with the desired motion state. Accordingly, it is necessary to study the influence law of error factors on the dynamic accuracy of the parallel mechanism in the design stage of the mechanism. The performance of the system is often significantly reduced, while the failure rate will increase or even cause major losses if these error factors are not fully taken into account. Therefore, the dynamic accuracy of a precision mechanism is worthy of in-depth study. Yun et al. [23] established a kinematics model of an 8-PSS/SPS parallel mechanism system via the stiffness model and Newton-Raphson method, where Kane's method was used to build up the dynamics model for analyzing the workspace, motion precision and dynamic characteristics. The size ranges of leg of a 3-DOF UPU parallel mechanism were obtained from sensitivity analysis in Ref. [24]. Based on the design and sensitivity analysis, the parallel mechanism was developed and various experiments have shown that the manipulator exhibited high accuracy and precision. In Ref. [25], the elastic dynamic equations of a flat-shaped parallel robot were derived, and the effects of two commonly used static balancing techniques on the dynamic performance of closed-chain mechanisms were deeply analyzed. The finite element method [26] was used to derive the elastodynamic equation of a plane high-precision linkage, and its dynamic characteristics and also natural frequency were analyzed. Taking the plane 3-PRR parallel mechanism as the object [27], a mechanism dynamics model with the clearance of the rotating pair was established. The root mean squared error (RMSE) was proposed as the index to quantify the effect of joint clearance on the dynamics of the system. As discussed above, the rigid body models were usually adopted as the research object, while the finite element method was used for analyzing the dynamic accuracy of the planar parallel mechanism. Few related studies on dynamic accuracy, caused by the elastic deformation of the spatial parallel mechanism, were analyzed due to the complexity of elastic dynamic model.
The elastic deformation of the branch chain of the mechanism will cause elastic vibration, which leads to the pose error of the moving platform [28,29]. Therefore, the study of flexible deformation is very important to improve the dynamic accuracy of the mechanism. Taking a 5PSS/UPU coupling mechanism as an example [30,31], a closed vector chain method is proposed to reduce the influence of geometric errors on dynamic accuracy. Combing the finite element method with Lagrange method, the rigid-flexible coupling elastic dynamic equation of the system was established, and the effect of dynamic accuracy of component flexibility on the parallel mechanism was further analyzed. Finally, the design principle for the mechanism parameter of the prototype was obtained.
This paper is organized as follows. In Section 2, the influence of geometric error on the dynamic accuracy of the 5-PSS/UPU parallel mechanism is studied, and the calibration method of geometric errors for the parallel mechanism is proposed. Section 3 establishes the rigidflexible coupled model and analyzes the influence of the flexibility of the links on the dynamic accuracy of the parallel mechanism. Subsequently, numerical simulations are carried out. The conclusions are outlined in Section 4.

Error Modeling of 5-PSS/UPS Parallel Mechanism
The basic structure of the 5-PSS/UPU parallel mechanism is shown in Figure 1. It consists of a fixed base, a mobile platform and six branches connecting the two. The six branches include five PSS joint branches and one UPU joint branch. The symmetrically distributed PSS branched chains are the power input of the parallel mechanism, while the UPU branched chain provides constraints on the mechanism. Each drive branch is composed of 1 prismatic pair and 2 spherical pairs. The structure diagram of the 5-PSS/UPU parallel mechanism is shown in Figure 2, which depicts the fixed coordinate system O-XYZ and moving coordinate system D-X D Y-D Z D . D and O are the centers of the moving and fixed platforms respectively. The Z axis is perpendicular to the static platform, the X axis is from the coordinate origin to point A 1 , and the Y axis is determined by the right-hand rule. The Z D axis is perpendicular to the moving platform upwards, the X D axis is from the coordinate origin to point C 1 , and the Y D axis is determined by the right-hand rule. E i -X Ei Y Ei Z Ei is the local coordinate system where the origin coincides with A i , the Y Ei axis coincides with A i B i , the X Ei axis is perpendicular to OE i , and the Z Ei -axis is determined by the right-hand rule, where i = 1, 2, 3, 4, 5, 6.
In the fixed coordinate system O-XYZ, the error model of the 5-PSS/UPU mechanism is established by vector method. The vector closed loop equation of a branch of the parallel mechanism can be expressed as where D is the position vector of the coordinate origin of the dynamic coordinate system; O D R is the rotation matrix from the dynamic coordinate system to the fixed coordinate system; s i denotes the direction vector of the i-th driving branch; n i represents the direction vector of the link of the i-th branch. L i and l i denote the length of the i-th flexible link and displacement of the i-th actuator, respectively.
Assuming that all error sources are small variables, differentiating Eq. (1) can be rewritten as T represents the differentiation of the rotation angle of the moving platform around the X, Y and Z axes.
Eq. (2) can be rewritten by taking the dot products of both sides with vector n T i as The error can be approximated as a kind of differentiation since it is a small variable. The error mapping model of the parallel mechanism can be expressed as where δW denotes the position and pose error of the moving platform. δΛ represents the input error of driving displacement. δL is the length error of a fixed link. δ D C is the position error of the spherical hinge at the moving platform. δA is the error of the apical position of the linear actuator. δS is the direction error of the linear actuator. J W is the Jacobian matrix of error transmission. G L is the coefficient matrix of link length error. G C is the coefficient matrix of the position error of spherical hinge at the moving platform. G A is the coefficient matrix of apical position error of linear actuator. G S is the error coefficient matrix of the linear actuator direction.
The Jacobian matrix of error transmission is invertible if the parallel mechanism is in a nonsingular configuration. Eq. (4) can be abbreviated as ,where K represents the mapping matrix of geometric error. δp denotes the geometric error source of the parallel mechanism. In δ D C, δA and δS, there are errors in the x, y and z axis directions. δp of Eq. (5) includes 11 error sources of the parallel mechanism. Among them, the z-direction error of linear actuator apex δA z is 0. Since the direction vector of linear actuator along the x, y and z axes of the three errors are not independent of each other, two independent error parameters (declinations δϕ and δθ ) are used to facilitate the analysis of error sensitivity, as shown in Figure 3.
In Figure 3, The relationship between δϕ, δθ and δs x , δs y , δs z can be expressed as where both ϕ and θ are known quantity. The expressions of δϕ and δθ can be derived from Eq. (6), so 11 error sources can be reduced to 9 independent error sources. (3)

Calibration of Geometric Errors
According to the structural characteristics and driving mode of the 5-PSS/UPU parallel mechanism, a vector calibration algorithm based on the inverse position solution is used to calibrate the geometric error. The positioning accuracy of actuator can be ignored since the high-precision screw module has high repetitive position accuracy.
When considering the existence of errors, Eq. (1) can be written as where D, O D R and n i in Eqs. (1) and (7) are identical in the same pose of the moving platform. By subtracting the two equations, the new equation is taken the dot products of both sides with vector n T i as Error sources δ D C i and δs i contain two independent error parameters. There are eight unknown error parameters in Eq. (8), and the five branches include 40 unknown error parameters. Five equations can be obtained by measuring posture of the moving platform, thereby at least 8 groups of different mechanism position parameters, which are required to be different as large as possible between them to reduce the coupling effect of each error parameter. According to actual condition, several sets of position parameters are measured and took the average of that to obtain higher calibration accuracy.
The vector calibration algorithm is used to identify the parameters of the 5-PSS/UPU parallel mechanism. Adopting Eq. (5) and the error data of Figure 4 as the actual size error, different position parameters can be obtained in terms of the actual structure size of the parallel mechanism, which is used to replace the actual measured position parameters. The actual structure size parameters and mechanism size parameters obtained by identification are regarded as the input and output of the calculation of vector calibration. Finally, the input and output parameters are compared to evaluate the accuracy of the calibration algorithm. The actuator branched chain 1 is selected for simulation calculation, and the pose parameters given by the actual structure size of the parallel mechanism ( Table 1). As shown in Table 2, the error values were obtained by the vector calibration algorithm according to Eq. (8). One can find that most of the calibration values obtained by the vector calibration method have an accuracy of more than 90%, which indicates the correctness and effectiveness of the calibration method. The reason for affecting the accuracy is that the driving input error is ignored. In addition, it is noteworthy that sensors also have some accuracy problems in real situation, so the accuracy of the actual calibration results will be lower than that of the simulation results. On the whole, the vector calibration method is generally effective for the calibration of 5-PSS/ UPU parallel mechanism.
To further verify the effectiveness of calibration for improving dynamic accuracy, the position and pose errors of the moving platform before and after the calibration were compared. The specific process is that the joint displacement is obtained by using the inverse kinematics solution in a given the trajectory of the moving platform, and then the motion trajectory error before and after calibration is compared after substituting it into the actual model. The error calibration values of all error sources obtained by the vector calibration algorithm are shown in Figure 4.
The pose errors of the moving platform after calibration can be calculated along the trajectory in Eq. (9). As shown in Figure 5, the calibrated dynamic accuracy of  the parallel mechanism has been greatly increased. In the fixed coordinate system, the position errors along the X-axis, Y-axis and Z-axis directions are reduced by 93.5%, 92.5%, and 91.7% respectively, while the angle errors around the three directions are reduced by 93.6%, 91.2%, and 86.8% respectively. It illustrates that the necessity of calibration of the parallel mechanism before it is put into work. On the other hand, the reliability of the proposed calibration algorithm has also been illustrated.

Dynamic Modeling of Rigid-flexible Coupling
Before establishing elastodynamic modeling of parallel mechanisms, two assumptions are proposed: one is that the displacement of the mechanism obtained by elastodynamic analysis is much smaller than that of the rigid body dynamics. Second, whilst there is a coupling effect between the rigid and elastic motion of the mechanism, the coupling relationship terms between them can be ignored in terms of a mechanism with less flexibility [32].

Spatial Beam Element Model
The five-branch links of the parallel mechanism are identically regarded as flexible parts, while the moving platform, fixed base, modules and each motion pair are all regarded as rigid bodies. Meanwhile, the influence of clearance of the spherical pair is ignored. The flexible link is regarded as the spatial beam element, as shown in Figure 6. The coordinate system o ij -x ij y ij z ij of beam element is introduced where the subscript i and j represent the i-th branch link and the j-th element. Each element has 2 nodes where each node has 9 elastic displacement degrees of freedom. The generalized coordinate δ ij of the spatial beam element can be expressed as    where δ ij1 − δ ij3 and δ ij10 − δ ij12 represent element node displacements. δ ij4 − δ ij6 and δ ij13 − δ ij15 represent element node rotation angles, δ ij7 − δ ij9 and δ ij16 − δ ij18 represent element node curvatures. The angular displacement around the axis is interpolated by a cubic difference function. The displacement and rotation angle in the other directions are interpolated by a quintic function. If W x x ij , t , W y x ij , t , W z x ij , t , ϕ x x ij , t , ϕ y x ij , t and ϕ z x ij , t are used to denote the elastic displacement and elastic angular displacement of any point in the element along x-axis, y-axis and z-axis directions, other parameters can be represented by generalized coordinate δ ij .
According to the boundary conditions of the beam element and the corresponding interpolation function, the displacement functions of the beam element can be expressed as where N ij1 , N ij2 , N ij3 , and N ij4 represent the vector matrix obtained by interpolation.
where x A(B) , y A(B) and z A(B) represent the rigid body displacements of the element nodes A(B) along the X-axis, Y-axis and Z-axis directions. θ x , θ y and θ z represent the rotation angle of the two element nodes around the X-axis, Y-axis and Z-axis.
Assuming that the mass of each element is concentrated on the axis, the kinetic energy of the element can be expressed as where ρ and S are the material density and cross-sectional area of the element. I x represents the polar inertia moment of the element cross-section along the X-axis direction.
If the shear deformation of beam element and the coupling between the axial and lateral displacement are ignored, the deformation potential energy of the spatial beam element can be expressed as where I y and I z represent the polar inertia moments of the element cross-section of the element along the Y-axis and Z-axis. K ij denotes the stiffness matrix of the element.
According to the Lagrange dynamics equation [33], by applying the kinetic energy and deformation potential energy of the space beam element, the elastodynamic equation of the element in the local coordinate system can be expressed as (13)  where M ij denotes the mass matrix of the element. F ij denotes the generalized external force term including force and torque. P ij represents the interaction force term caused by the connection of the beam element. Q ij is the rigid body inertial force term.
To facilitate elastic modelling, the elastodynamic equations in the local coordinate system are converted to the base coordinate system. As shown in Figure 7, the generalized coordinates δ ij in the base coordinate system can be expressed as

Elastodynamic Modeling of Fixed-length Link
The fixed-length link is regarded as a flexible rod while the rest are assumed to be rigid parts. As shown in Figure 8, the link is divided into n elements, which are numbered 1, 2, ..., n + 1 in sequence. Each element is connected adjacently in turn. The elastic displacement of the right end of the j-th element is consistent with that of the left end of the (j + 1)-th unit, 1 ≤ j ≤ n − 1.
Connecting the left and right ends of the fixed-length link are spherical pairs, thereby the curvature of the end beam element is 0, that is, δ (i1)7 = δ (i1)8 = δ (i1)9 = 0, δ (in)16 = δ (in)17 = δ (in)18 = 0. Synthetically, the generalized coordinate q i of the link can be obtained from the displacement relationship between the elements.   where

Kinematics and Dynamics Constraint Equation
The kinematic constraint of the system is the deformation coordination condition. Compared with the five fixed-length links, the moving platform can be regarded as a rigid body due to relatively larger rigidity. In this parallel mechanism, the displacement of the hinge point of the fixed-length link connected to the moving platform (20) is consistent with that of the corresponding point of the moving platform, as shown in Figure 9.  The transformation matrix from coordinate system Eq. (24) is the kinematic constraint equation of the 5-PSS/UPU parallel mechanism. Meanwhile, the dynamic constraint equation needs to be meet, i.e., the external force and inertial force of the moving platform must be balanced with the force exerted by the fixed-length link on the moving platform. Ignoring the coupling relationship between the rigid body motion and elastic motion of the moving platform, the dynamic constraint equation of the moving platform can be expressed as where M D is mass matrix of the moving platform. q aD denotes the acceleration array of the moving platform. F i is the combined force array of the branch-chain on the (23) sin δα D ≈ δα D , sin δβ D ≈ δβ D , sin δγ D ≈ δγ D , cos δα D ≈ 1, cos δβ D ≈ 1, cos δγ D ≈ 1. Figure 9 Constraints of system kinematics moving platform, F w is the external force array on the moving platform. The dynamic constraint equation in term of q aD =q D +q rD can be rewritten as:

Elastic Dynamics Equation of Parallel Mechanism System
On the basis of the kinematics and dynamics constraint relations of each component of the system, the elastodynamic equations of each link are assembled from the elastodynamic equations of each branch chain. To facilitate the analysis, a generalized coordinate of branch chain q i (i = 1, 2, 3, 4, 5) is defined as The relation equation q i =R iqi between q i and q i can be obtained from the kinematics constraint relation equation. Substituting the equation into Eq. (22), the new equation is taken the dot products of both sides with vector R T i as Meanwhile, to facilitate the assembly of the system equations, the generalized system coordinate q i and q i can be written as The mapping relation between generalized system coordinates q and q i can be written as: By assembling and superposing all the elastodynamic equations of the branch links as well as considering the effect of the system damping, the elastodynamic equations of the whole system in the base coordinate system can be expressed as where M is the total mass matrix of the system. C is the coefficient matrix of the damping effect. K is the total stiffness matrix of the system, Q is the generalized force matrix of the system.

Analysis of Dynamic Accuracy
According to the force simulation results analyzed the Newton-Euler method of each rigid body component and the structure size of each component, one can know that the dynamic accuracy of the parallel mechanism is affected by the deformation of the five flexible links. Thus, the elastodynamic model is actually a rigid-flexible coupling dynamics model [34]. To verify the correctness of the elastodynamic model, the generalized forces are obtained via the co-simulation with HYPERMESH, ANSYS and ADAMS software. First, the flexible links are meshed by hexahedron elements in HYPERMESH, and the mnf files of the meshed links are obtained via ANSYS. Then, loading the link files, rigid links are replaced with flexible bodies in ADAMS. The simulation model is shown as Figure 10.
The structure and material parameters of the 5-PSS/ UPU parallel mechanism are shown in Table 3. The external force F w = [0 0 500 0 0 0] is imposed on the center of the moving platform. The motion trajectory of the moving platform is given in Eq. (32). The solution time and solution integration step are set to t = 2 s and t = 0.001 s.  The actuator displacements are obtained by substituting the motion and dynamic constraints into the deformation equation, which is used as the input of the simulation model to obtain the actuator forces. The comparison between theoretical actuator forces calculated via Eq. (31) and that of simulation is shown in Figure 11. One can find that the errors of the actuator force of the theoretical value and the simulation value are less than 4%, which proves the correctness of the rigid-flexible coupling dynamic model. The main reason for this error is that the number and form of the unit division of the simulation are different from the theory.
The Newmark method [35] is used to calculate the position and pose errors of the moving platform caused by the elastic deformation of the branch rod, as shown in Figure 12.
Obviously, it can be seen from Figures 12 and 13 that violent oscillations of the pose of the moving platform are emerged due to the elastic deformation of the branch link. The dynamic accuracy of the moving platform is greatly affected because of these oscillations. When the time t = 1.772 s and t = 1.264 s, the maximum displacement error occurs along the X-axis and Y-axis directions and the Z-axis direction. At the time t = 1.791 s and t = 1.248 s, the maximum angle error arises in the β and γ direction.
To quantitatively analyze the influence of factors, such as the cross-sectional area of the branch link, the elastic modulus and the mass of the moving platform on the dynamic accuracy of the flexible parallel mechanism, the evaluation of above factors is obtained by calculating the (32)    x = 0, β = π/18 sin (2πt), y = 0, γ = π/18 sin (2πt), z = 0.554 + 0.08 sin (πt).
dimensionless root mean square error [36]. The influence index of dimensionless root mean square error can be written as where x ai represents the output of the moving platform in real condition, x i denotes the output of the moving platform in ideal condition. RMS(x ai − x i ) represents the root mean square error of the output of the moving platform in real condition. RMS(x i ) represents the root mean square error of the output of the moving platform in ideal condition. N is the simulated sample size. To make the results more intuitive and representative, the dimensionless root mean square error influence index of the displacement, velocity and acceleration in the Z-axis direction of the centroid of the moving platform will be calculated. It can be seen from Eq. (33) that the dynamic accuracy is proportional to the dimensionless root mean square error. In terms of the different conditions of the cross-sectional area, elastic modulus and moving platform mass of the branch link, the quantitative calculation of the non-dimensional root mean square error of the output of the parallel mechanism is shown in Figure 14. Figure 14(a) shows that the non-dimensional rootmean-square error of the position, velocity and acceleration of the moving platform of the parallel mechanism are proportional to the cross-sectional area. When the cross-sectional area increases from 400π mm 2 to 900π mm 2 , the dimensionless root-mean-square error influence indexes of displacement, velocity and acceleration decrease from 0.0284%, 5.6235% and 59.6772% to 0.0019%, 0.8277%, and 8.6590%, respectively. It shows that the larger cross-sectional area of the branch rod corresponds to the smaller output error of the mechanism and the higher dynamic accuracy.
Similarly, the non-dimensional root-mean-square error of the position, velocity and acceleration of the moving platform of the parallel mechanism are proportional to the elastic modulus, as shown in Figure 14(b). When the elastic modulus increases from 70 GPa to 206 GPa, the dimensionless root-mean-square error influence indexes of displacement, velocity, and acceleration decrease from 0.0284%, 5.6235% and 59.6772% to 0.0029%, 0.7950% and 20.0719%, respectively. It illustrates that the larger elastic modulus of the branch rod corresponds to the smaller output error of the mechanism and the higher dynamic accuracy.
In Figure 14(c), the non-dimensional root-mean-square error of the position, velocity and acceleration of the moving platform of the parallel mechanism are proportional to the mass of moving platform. When the mass of moving platform increases from 20 kg to 60 kg, the dimensionless root-mean-square error influence indexes of displacement, velocity, and acceleration increase from 0.0284%, 5.6235% and 59.6772% to 0.0708%, 13.6946% and 124.250%, respectively. It shows that the greater mass of the moving platform corresponds to the greater the output error of the mechanism, the more unstable the vibration amplitude of the mechanism and the lower dynamic accuracy.

Conclusions
Taking a 5-PSS/UPU parallel mechanism as an example, on the basis of an error model and the mapping law of geometric error to the dynamic accuracy, a vector calibration algorithm was proposed to reduce the position and pose error along the whole motion trajectory. Then, the elastic dynamic model was established via analyzing the elastic deformation of the components. Furthermore, the effect of the flexibility of the components on the dynamic accuracy of the parallel mechanism was analyzed.
(1) The influence of each error source on the dynamic accuracy of the parallel mechanism is analyzed based on the geometric error model. After the main error sources of calibration with the closed-loop vector method, the position error along the x-axis, y-axis, and z-axis directions are reduced by 93.5%, 92.5%, 91.7% respectively, while the angle errors around the x-axis, y-axis, and z-axis directions are reduced by 93.6%, 91.2%, and 86.8% respectively. It shows that the dynamic accuracy of the parallel mechanism can be improved via reasonable compensation of geometric errors. (2) The elastic dynamics equation of the parallel mechanism is established by combining the finite element method with Lagrange method, which is solved by the Newmark direct integration method. Furthermore, the correctness of the rigid-flexible coupling dynamic model is verified by the co-simulation with HYPERMESH, ANSYS and ADAMS software. (3) Numerical simulation results show that the crosssectional area of the branch connecting rod and the elastic modulus are proportional to the dimensionless root mean square error of the output of the moving platform, while the mass of the moving platform is inverse proportional to the dimensionless root mean square error of the output of the moving platform. With the requirements of rigidity and assembly, a reasonable parameter combination can effectively improve the dynamic accuracy of the mechanism.