Research on High Precision Stiffness Modeling Method of Redundant Over-Constrained Parallel Mechanism

Traditional stiffness modeling methods do not consider all factors comprehensively, and the modeling methods are not unified, lacking a global stiffness model. Based on screw theory, strain energy and the virtual work principle, a static stiffness modeling method for redundant over-constrained parallel mechanisms (PMs) with clearance was proposed that considers the driving stiffness, branch deformation, redundant driving, joint clearance and joint contact deformation. First, the driving stiffness and branch deformation were considered. According to the strain energy and Castiliano’s second theorem, the global stiffness matrix of the ideal joint mechanism was obtained. The offset of the branch was analyzed according to the restraint force of each branch. The mathematical relationship between the joint clearance and joint contact deformation and the end deformation was established. Based on the probability statistical model, the uncertainty of the offset value of the clearance joint and the contact area of the joint caused by the coupling of the branch constraint force was solved. Finally, taking a 2UPR-RR-2RPU redundant PM as an example, a stiffness simulation of the mechanism was carried out using the finite element method. The research results show that the high-precision stiffness modeling method proposed in this paper is correct, and provides an effective method for evaluating the stiffness performance of the PM.


Introduction
The stiffness reflects the displacement of the end of the mechanism relative to the expected position under external forces. The stiffness of the parallel mechanism (PM) often determines the kinematic positioning accuracy and kinematic characteristics of the end effector. In addition, insufficient stiffness may reduce the natural frequency and dynamic performance of the system. Therefore, the stiffness performance should be evaluated in the product design phase.
In terms of stiffness modeling methods of PM, the main methods include screw theory [1], the semi-analytical method [2], the virtual joint method [3], strain energy [4,5], the structural stiffness matrix method [6], etc. The factors considered in the modeling method are gradually developing toward moving platform flexibility [7], driving stiffness and branch flexibility [8,9] and gravity [10][11][12], making the accuracy of analytical model building more accurate and closer to perfection. Compared with non-redundant mechanisms, redundant mechanisms have many advantages, such as avoiding kinematic singularity, improving stiffness distribution, enhancing bearing capacity and improving dynamic characteristics. Therefore, redundant driving branches are usually added to the mechanism to improve the stiffness distribution in its workspace [13][14][15][16][17].
When the mechanism has clearance joints, the end position depends on the attitude, load, geometric parameters of the joint and the joint force, and the relative position inside (2) The driving stiffness is analyzed from the perspective of material mechanics. The strain energy of each branch under the constraint force is calculated and the functional relationship between the strain energy of the branch and the constraint force of the branch is derived. (3) Based on Castigliano's second theorem, the partial derivatives of the strain energy of each branch on the constraint force of the branch are calculated, respectively. The functional relationship between the deformation of the end of the branch and the constraint force of the branch is obtained, namely the flexibility matrix of the branch. (4) According to the deformation compatibility equation of the mechanism and virtual work principle, the stiffness matrix of the whole mechanism and the constraint force of each branch are solved and the end deformation caused by the driving stiffness and branch deformation is obtained. (5) Based on the constraints of each branch, the position distribution of the clearance joints and the contact area of the joints are predicted by the probability statistical model. It solves the problem of the modeling error caused by coupling constraints of over-constrained mechanisms. The deformation caused by the joint clearance and joint contact deformation along the constraint direction at the connection between each branch and the moving platform end is analyzed. Based on the principle of virtual work, the end deformation interval caused by joint clearance and joint contact deformation is obtained. (6) Assuming that there is no coupling between rigid body displacement and elastic deformation, the four factors are regarded as mutually independent and can be linearly superposed. A high-precision static stiffness model considering multiple factors is obtained.

General Equation of High-Precision Stiffness Model of Redundant Over-Constrained PMs Considering Multiple Factors
Based on the above assumptions and modeling process, the high-precision static stiffness model considering multiple factors can be expressed as ∆X = ∆X ς + ∆X cle + ∆X con = (k a,c ) −1 w + J * nac K cle nac + J * rac K cle rac + J * nac K con nac τ nac + J * rac K con rac τ rac (1) In the Equation, ∆X represents the total deformation interval in the direction of six DOFs of the end, and ∆X ς represents the end deformation caused by the driving stiffness, branch deformation and redundant driving. ∆X cle represents the end deformation interval caused by the joint clearance corrected based on the statistical model. ∆X con represents the end deformation interval caused by the joint contact deformation corrected based on the statistical model. w represents the six-dimensional load of the moving platform.
In the Equation, J * na represents the mapping matrix of the non-redundant branch drive force and external load. J * nc represents the mapping matrix of the non-redundant branch constraint force and external load. J * ra represents the mapping matrix of the redundant branch drive force and external load. J * rc represents the mapping matrix of the redundant branch constraint force and external load. K cle na represents the offset matrix of the end of the branch caused by the driving force of the non-redundant branch with clearance. K cle nc represents the offset matrix of the end of the branch caused by the constraint force of the non-redundant branch with clearance. K cle ra represents the offset matrix of the end of the branch caused by the driving force of the redundant branch with clearance. K cle rc represents the offset matrix of the end of the branch caused by the constraint force of the redundant branch with clearance. K con na represents the contact stiffness matrix of the joint at the end of the branch caused by the driving force of the non-redundant branch. K con nc represents the contact stiffness matrix of the joint at the end of the branch caused by the restraint force of the non-redundant branch chain. K con ra represents the contact stiffness matrix of the joint at the end of the branch caused by the driving force of the redundant branch. K con rc represents the contact stiffness matrix of the joint at the end of the branch caused by the restraint force of the redundant branch chain. τ na represents the non-redundant branch driving force matrix. τ nc represents the non-redundant branch constraint matrix. τ ra represents the redundant branch driving force matrix. τ rc represents the redundant branch constraint matrix.

Stiffness Matrix Derivation and Restraint Force Analysis of Redundant PMs with Ideal Joints
The force analysis of over-constrained PMs is statically indeterminate, and the structural stiffness of each branch needs to be considered. In this paper, the moving platform was assumed to be rigid, and the spatial composite elastic deformation of the branch bending, stretching and torsion was considered. The driving force/restraint stiffness matrix of the branch is defined as the mapping relationship between the driving force/restraint amplitude and the elastic deformation of the end along the driving force/restraint direction.
The moving platform of the mechanism is acted on by a six-dimensional external force w. According to the static equilibrium equation of the moving platform, we can obtain τ na , τ nc represent the column vectors composed of the driving force and constraint force amplitude of each non-redundant branch, respectively. τ na , τ nc represent the column vectors composed of the driving force and constraint force amplitude of each redundant branch, respectively. J a = J na J ra = $ T a,1 · · · $ T a,k · · · $ T a,i T (4) In Equations (4) and (5), $ a,i represents the driving force screw of branch i acting on the moving platform, and $ c,i represents the constraint force screw of branch i acting on the moving platform. n represents the number of all constrained screws of non-redundant branches in the mechanism, m, and i represents the number of all constraint forces of redundant branches in the mechanism.
Assuming that the branch driving force τ ai and the branch terminal deformation ∆q ai caused by the driving force and the branch restraint force τ ci and the deformation ∆q ci of the branch terminal caused by the constraint force satisfy the basic linear elasticity assumption, we can obtain τ na = k na ∆q na ; τ ra = k ra ∆q ra (6) τ nc = k nc ∆q nc ; τ rc = k rc ∆q rc where k na = diag k a1 k a2 · · · k ak , k nc = diag k c1 k c2 · · · k cn and k aj (1 ≤ j ≤ k) are defined as the stiffness of the non-redundant branch in the direction of the driving force. k cj (1 ≤ j ≤ n) is defined as the stiffness of the non-redundant branch in the direction of the constraint force. k ra = diag k ak k a(k+1) · · · k ai , k rc = diag k cn k c(n+1) · · · k cm and k aj (k ≤ j ≤ i) are defined as the stiffness of the redundant branch in the direction of the driving force. k cj (n ≤ j ≤ m) is defined as the stiffness of the redundant branch in the direction of the constraint force.
In Equation (8), ∆X ς = ∆x ∆y ∆z ∆α ∆β ∆γ T represents the micro-deformation of the moving platform under the external forces when only the driving stiffness and branch deformation are considered. Sorting Equation (3) and Equations (6)-(8) simultaneously, we can obtain In order for Equation (9) to be established, it is necessary to satisfy Combining Equations (3), (6), (7), (10) and (11), the relationship between the external force on the moving platform and the micro-deformation of the moving platform when considering the driving stiffness and the deformation of the branch can be obtained.
In the Equation, k a,c represents the stiffness matrix of the over-constrained PM considering the driving stiffness of the mechanism and the deformation of the branches.
Combining Equations (5), (6) and Equations (9)-(11), the expressions of the driving force, constraint force and external force of the moving platform can be obtained.
where subscript a represents the driving force, c represents the constraint force (i = a, c), n represents a basic branch and r represents a redundant branch. Rearranging the above equation can obtain

Stiffness Modeling Method of Over-Constrained Redundant PMs Considering Joint Clearance and Joint Contact Deformation
The offset/deflection of a branch with joint clearance under an external force is affected by multiple factors, such as the clearance size, attitude of mechanism, redundant branch and branch constraint force/couple. Over-constrained redundant PMs generally belong to statically indeterminate structures, and the constraints between branches are complex. Under a certain static attitude and load, the joint with clearance will produce tiny line Sensors 2023, 23, 5916 6 of 33 deflection/angular deflection in its constraint direction under the action of constraint force/torque. It is difficult to judge the joint contact area caused by the branch compound constraint, and it is difficult to establish an accurate mathematical model. In this section, based on certain assumptions, the approximate change interval of the mechanism end is obtained when the joint clearance and joint contact deformation are considered. The specific process is as follows: (1) Firstly, the drive force f i a , constraint force f i c and constraint couple f i m of redundant and non-redundant branches are solved by using the system stiffness matrix of the ideal joint.
(2) It is assumed that the additional motion of the branch with clearance is in the same direction as the restraint force of the branch. According to the size and direction of the restraint force of each branch, the additional motion direction of the branch with gaps can be determined. The static stiffness modeling process of redundant over-constrained PMs when considering joint clearance is shown in Figure 1. (3) First, the coupling effect of the constraint and constraint couple is ignored and the constraint is considered as a pure force or a pure couple. In most working conditions, the drive force of each branch is much greater than the constraint force, so it is assumed that the clearance joint offset/deflection and joint contact area in the drive force direction reach the maximum. It is also assumed that the ratio of the drive force f i a in the branch with clearance to the offset ∆Xi cle f a along the drive force direction and the ratio of the constraint force f i c to the offset ∆Xi cle f c along the constraint direction are equal and regarded as constant (which can be approximately equivalent to a spring system with equal stiffness)   (3) First, the coupling effect of the constraint and constraint couple is ignored and the constraint is considered as a pure force or a pure couple. In most working conditions, the drive force of each branch is much greater than the constraint force, so it is assumed that the clearance joint offset/deflection and joint contact area in the drive force direction reach the maximum. It is also assumed that the ratio of the drive force   Therefore, the offset of each branch along the constraint direction is The offset/deflection of the end of each branch caused by the joint clearance and joint contact deformation is mapped to the deformation of the moving platform by the principle of virtual work. When the deviation of the restraint direction of the clearance joint does not reach the maximum value, the direction also does not reach the contact state and there is no contact deformation between the joints. The static stiffness modeling process of the redundant PM when considering the contact deformation of the joint is shown in Figure 2 (∆Ki cle f c represents the contact deformation caused by the constraint of branch i, ∆Ki cle f a represents the contact deformation caused by the driving force of branch i and ∆Ki cle m c represents the angular contact deformation caused by the constraint couple of branch i.).
a a c c c c con con con con con (5) In order to solve the uncertainty of the gap joint offset value and joint contact area caused by the branch-coupling force, a random method can be used to describe the variation characteristics of the gap and contact area. i ς is the correction coefficient of random physical variables. Based on the maximum value of the clearance joint offset value and joint contact area, a specific probability distribution attribute is assigned. Using the Monte Carlo method, the probability distribution of the stiffness of the mechanism under a given position and attitude is obtained. The influence of joint clearance and joint contact deformation on the end stiffness is predicted through the probability statistical model, and the variation range of the end stiffness is approximately determined. (5) In order to solve the uncertainty of the gap joint offset value and joint contact area caused by the branch-coupling force, a random method can be used to describe the variation characteristics of the gap and contact area. ς i is the correction coefficient of random physical variables. Based on the maximum value of the clearance joint offset value and joint contact area, a specific probability distribution attribute is assigned. Using the Monte Carlo method, the probability distribution of the stiffness of the mechanism under a given position and attitude is obtained. The influence of joint clearance and joint contact deformation on the end stiffness is predicted through the probability statistical model, and the variation range of the end stiffness is approximately determined.

Example Analysis of 2UPR-RR-2RPU Redundant PM Stiffness Considering Multiple Factors
In this paper, a 2UPR-RR-2RPU redundant over-constrained PM was taken as an example to verify the accuracy of the proposed stiffness modeling method. The mechanism is a loop decoupled PM that adopts the method of adding two redundant branches to improve the overall bearing capacity and deformation resistance of the output end.

Configuration Analysis of 2UPR-RR-2RPU Redundant PM
As shown in Figure 3, the 2UPR-RR-2RPU PM consists of two UPR branches, one RR branch and two RPU branches for a total of five branches. The fixed coordinate system o − x 0 y 0 z 0 is established at the center point of the fixed platform. The x 0 axis along the oA 4 direction, the y 0 axis along the oA 1 direction, and the z 0 axis is determined according to the right-hand rule. The moving coordinate system o 1 − x 1 y 1 z 1 is established at the center point of the fixed platform. The x 1 axis along the o 1 B 4 direction, the y 0 axis along the o 1 B 1 direction, and the z 1 axis is determined according to the right-hand rule. The middle branch of RR is a just-constrained branch and its axis of rotation connected to the fixed platform is along the direction of x 0 in the fixed coordinate system o − x 0 y 0 z 0 . The revolute pair provides the x-direction DOF of rotation of the moving platform around the axis. Its axis of rotation connected to the moving platform is along the direction y 1 in the moving coordinate system o 1 − x 1 y 1 z 1 . The R pair provides the y-direction DOF of rotation of the moving platform around the axis. Among them, the R pairs of the two UPR branches whose respective U pairs are connected to the fixed platform and the R pairs of the RR just-constrained branch that are connected to the fixed platform are coaxial. The respective R pairs of the two UPR branches are parallel to the R pairs of the RR just-constrained branch, which is connected to the moving platform. Among them, the R pair of the two RPU branches whose U pairs are connected to the moving platform is coaxial with the R pair that is connected to the RR just-constrained branch and moving platform, and the respective R pairs of the two RPU branches are, respectively, parallel to the R pairs that are connected to the RR just-constrained branches and fixed platform. The R pair provides the y-direction DOF of rotation of the moving platform around the axis. Among them, the R pairs of the two UPR branches whose respective U pairs are connected to the fixed platform and the R pairs of the RR just-constrained branch that are connected to the fixed platform are coaxial. The respective R pairs of the two UPR branches are parallel to the R pairs of the RR just-constrained branch, which is connected to the moving platform. Among them, the R pair of the two RPU branches whose U pairs are connected to the moving platform is coaxial with the R pair that is connected to the RR just-constrained branch and moving platform, and the respective R pairs of the two RPU branches are, respectively, parallel to the R pairs that are connected to the RR just-constrained branches and fixed platform.

Full Jacobian Matrix Solution for 2UPR-RR-2RPU Redundant PM
In order to obtain the full Jacobian matrix of the mechanism, the constraint Jacobian matrix and the driving Jacobian matrix of the mechanism should be solved first.   In order to obtain the full Jacobian matrix of the mechanism, the constraint Jacobian matrix and the driving Jacobian matrix of the mechanism should be solved first.
A i B i (i = 1, 2, 3, 4) represents the direction vector of the kinematics pair of the four driving branches, i = 1, 2 is the non-redundant branch and i = 3, 4 is the redundant branch. A i indicates the coordinates of the joint center between each branch and the fixed Sensors 2023, 23, 5916 9 of 33 platform, and B i indicates the center coordinates of the joint that connects each branch and the moving platform.
It can be represented by a screw as In Equation (16), Let the rotation angle of the mechanism around the rotation axis x o be α and the rotation angle around the rotation axis y 1 be β, and let a, b be the radius of the fixed platform and the moving platform, respectively.
According to the screw theory, the screw that is reciprocal to the non-redundant UPR branched motion screw is its constraint screw.
Redundant UPR branches and non-redundant UPR branches are symmetrically arranged in the plane of the base coordinate system with respect to the RR just-constrained branch. In order to obtain the constraint screw f $ 3 c1 , m $ 3 c2 and driving force screw $ a3 of the redundant UPR branch, it is only necessary to invert the constraint screw f $ 1 c1 , m $ 1 c2 and a 1 , b 1 in the driving force screw $ a1 of the non-redundant UPR branch.
According to the screw theory, the screw that is reciprocal to the non-redundant RPU branched motion screw is its constraint screw The redundant RPU branches and the non-redundant RPU branches are symmetrically arranged on the base coordinate system y o z o plane with the RR just-constrained branches. In order to obtain the constraint screw f $ 4 c1 , m $ 4 c2 and driving force screw $ a4 of the redundant RPU branch, it is only necessary to invert the constraint screw f $ 2 c1 , m $ 2 c2 and a 1 , b 1 in the driving force screw $ a4 of the non-redundant UPR branch.
The RR intermediate branch-constrained screw is expressed in the fixed coordinate system as From Equation (3), it can be obtained that the non-redundant constrained Jacobian matrix J nc and the redundant constrained Jacobian matrix J rc of the 2UPR-RR-2RPU PM are The following assumptions are made for the stiffness model: the weight of all components is ignored; all joint models are frictionless; the moving platform is assumed to be a rigid body; and the spatial composite elastic deformation of the branch is considered, including tensile, shear, bending and torsional deformation.
The branch driving stiffness can be simplified as a series spring system, and the driving stiffness coefficient k ai can be expressed as In the Equation, it is expressed as a U pair along the rod axial stiffness k ai,1 , pendulum rod stiffness k ai,2 , telescopic rod stiffness k ai,3 , and R pair seat along the rod axial stiffness k ai,4 . Among them, k ai,2 , k ai,3 are constants. k ai,2 can be calculated by k ai,2 = EA/(q i − l 1 ), EA represents the tensile modulus of the pendulum rod, l 1 represents the length of the telescopic rod and q i is the distance between the joint center connected with the moving platform and the joint center connected with the fixed platform.

Solution of Branch Restraint Stiffness Matrix of 2UPR-RR-2RPU Redundant PM
The force diagram of the UPR branch is shown in Figure 4. The branch coordinate system o 1 − x 1 y 1 z 1 is established at the rotation center of the branch close to the moving platform. The y 1 axis is parallel to the axis direction of the R pair and the z 1 axis is along the axis of the rod. x 1 is determined by the right-hand screw rule. In order to analyze the deformation of the branch under the restraint force, the restraint force f 11 passing through the U pair center and parallel to the R pair direction is translated to o 1 − x 1 y 1 z 1 . According to the force translation theorem, there will be an additional force couple passing through the origin of the branch coordinate system o 1 − x 1 y 1 z 1 and the direction parallel to the axis x 1 with magnitude m 12 . There will also be an additional constraint force passing through the origin of the branch coordinate system o 1 − x 1 y 1 z 1 and the direction parallel to the axis y 1 with magnitude m 11 . The branch is also constrained by a couple whose axis is parallel to the U pair normal and whose magnitude is m 11 . f 1x represents the force in the x-axis direction at the cross-section of the RR branch, f 1y represents the force in the y-axis direction at the cross-section, f 1z represents the force in the z-axis direction at the cross-section, m 1x represents the force couple in the x-axis direction at the cross-section, and m 1z represents the force couple in the z-axis direction at the cross-section.
The internal force of the section where the l i is displaced from the negative direction of the coordinate system z 1 to the UPR branch can be expressed as In the Equation, i l 1 is the length of the linear actuator and L 1 is the total length of the UPR branch and joint center that connect the moving and fixed platforms. When i l = L 1 , In the Equation, m ς is the direction vector of the constraint couple m 11 , e 1 is the direction vector of the axis x 1 in the coordinate system o 1 − x 1 y 1 z 1 , e 2 is the direction vector of the axis y 1 in the coordinate system o 1 − x 1 y 1 z 1 , e 3 is the direction vector of the axis z 1 in the coordinate system o 1 − x 1 y 1 z 1 and m 12 = f 11 L 1 , where L 1 is the length of the chain rod.
The strain energy of UPR branches can be expressed as In the Equation, l i1 is the length of the linear actuator and L 1 is the total length of the UPR branch and joint center that connect the moving and fixed platforms. When l i = L 1 , In the Equation, are the elastic modulus and shear modulus of the telescopic rod and lead screw. A i1y , A i2y are the effective shear area along the axis y i of the telescopic rod and the lead screw. I i1x , I i2x are the inertia moment of the cross-section of the telescopic rod and the lead screw about the axis x i . J i1 , J i2 are the cross-sectional polar inertia moments of the telescopic rod and the lead screw, respectively.
According to Castigliano's theorem, we can obtain Arranging the above Equation into matrix form, we can obtain In the Equation the stiffness matrix K UPR of the UPR branch is the inverse of the compliance matrix C UPR .
It can be seen from the UPR branch compliance matrix that the compliance matrix is not a diagonal matrix. There is a coupling relationship between the constraint force f 11 and the constraint couple m 11 and their corresponding elastic deformations. The size of the compliance matrix is related to the direction of the normal vector m ς of the U pair in the UPR branch.
The force diagram of the RPU branch is shown in Figure 5. The branch coordinate system o 2 − x 2 y 2 z 2 is established at the U pair center of the branch close to the moving platform. The axis y 1 is parallel to the axis direction of the R pair of the U pair close to the moving platform and the axis z 1 is along the axis direction of the rod. x 1 is determined by the right-hand screw rule. According to the screw theory, the RPU branch is subject to a constraint force passing through the origin of the branch coordinate system o 2 − x 2 y 2 z 2 and the direction parallel to the axis y 1 with a magnitude of f 21 . It is also subject to a constraint couple whose axis is parallel to the U pair normal with a magnitude of m 21 . f 2x represents the force in the x-axis direction at the cross-section of the RR branch, f 2y represents the force in the y-axis direction at the cross-section, f 2z represents the force in the z-axis direction at the cross-section, m 2x represents the force couple in the x-axis direction at the cross-section, and m 2z represents the force couple in the z-axis direction at the cross-section.
It can be seen from the UPR branch compliance matrix that the compliance matrix is y is parallel to the axis direction of the R pair of the U pair close to the moving platform and the axis z 1 is along the axis direction of the rod. x 1 is determined by the right-hand screw rule. According to the screw theory, the RPU branch is subject to a constraint force passing through the origin of the branch coordinate system o -x y z 2 2 2 2 and the direction parallel to the axis y 1 with a magnitude of f 21 . It is also subject to a constraint couple whose axis is parallel to the U pair normal with a magnitude of m 21 .   The internal force of the section where the displacement from the RPU branch to the negative direction of the branch coordinate system z 2 is l i can be expressed as In the Equation, m ξ is the direction screw of the constraint couple m 21 . e 1 is the direction vector of the axis x 1 of the coordinate system o 2 − x 2 y 2 z 2 . e 2 is the direction vector of the axis y 1 of the coordinate system o 2 − x 2 y 2 z 2 . e 3 is the direction vector of the axis z 1 of the coordinate system o 2 − x 2 y 2 z 2 .
The strain energy of RPU branches can be expressed as In the Equation, l i1 is the driving displacement of the driving pair when l i = l 2 : According to Castigliano's theorem, we can obtain In the Equation, Arranging the above Equation into matrix form, we can obtain In the Equation, The stiffness matrix of the RPU branch is the inverse of the compliance matrix It can be seen from the RPU branch compliance matrix that the compliance matrix is not a diagonal matrix. There is a coupling relationship between the constraint force f 21 , the constraint couple m 21 and their corresponding elastic deformations. The size of the compliance matrix is related to the direction of the normal vector m ξ of the U pair in the RPU branch.
The force figure of the RR just-constrained branch is shown in Figure 6. The branch coordinate system o 3 − x 3 y 3 z 3 is established at the revolute pair center of the branch close to the moving platform. The axis y 1 is parallel to the direction of the R pair axis, and the axis z 1 is along the axis of the rod. x 1 is determined by the right-hand screw rule. In order to facilitate an analysis of the deformation of the rod under the restraint force, the restraint force f 31 that passes the R pair center near the fixed platform and is parallel to the R pair axis direction that is close to the moving platform is translated to o 3 − x 3 y 3 z 3 . According to the force translation theorem, there will be an additional force couple passing through the origin of the branch coordinate system o 3 − x 3 y 3 z 3 and the direction parallel to the axis x 3 with a magnitude of m 32 . There will also be an additional constraint force passing through the origin of the branch coordinate system o 3 − x 3 y 3 z 3 and the direction parallel to the axis y 3 with a magnitude of f 31 . In addition, the branch is also subjected to a constraint couple whose magnitude is m 31 , with the axis parallel to the direction. It is also subjected to a constraint force f 33 and a constraint force whose axis is parallel to the direction of y 3 with a magnitude of f 32 . f 2x represents the force in the x-axis direction at the cross-section of the RR branch, f 2y represents the force in the y-axis direction at the cross-section, f 2z represents the force in the z-axis direction at the cross-section, m 2x represents the force couple in the x-axis direction at the cross-section, and m 2z represents the force couple in the z-axis direction at the cross-section. The internal force of the section where the displacement from the RR just-constrained branch to the negative direction of the branch coordinate system z 1 is i l can be expressed as  The internal force of the section where the displacement from the RR just-constrained branch to the negative direction of the branch coordinate system z 1 is l i can be expressed as In the Equation, m 32 = f 31 L 3 , where L 3 is the rod length of the branch. The strain energy of the RR-constrained branch can be expressed as According to Castigliano's theorem, we can obtain Arranging the above Equation into matrix form, we can obtain In the Equation, The stiffness matrix K RR of the RR just-constrained branch is the inverse matrix of the flexibility matrix C RR : It can be known from the RR branched compliance matrix that the compliance matrix is a diagonal matrix and that there is no coupling relationship between the constraint force, constraint couple and their corresponding elastic deformation.

Stiffness Analysis of 2UPR-RR-2RPU Redundant PM Considering Joint Clearance and Joint Contact Deformation
The precise mathematical model of the over-constrained redundant PM with joint clearance is very complicated. The tiny deformation transmitted by each branch with joint clearance to the moving platform is affected by many factors, such as the driving force, posture and attitude, clearance size and external load. The position and attitude of the moving platform is not completely calculated by the inverse kinematic solution, so the moving platform has an active status in the tiny working space, resulting in an uncertain position of the moving platform. An excessive constraint force between joints under large-load conditions will lead to slight penetration deformation, and the penetration deformation of joints will also affect the stiffness of the moving terminal. This cannot be ignored in the field of high precision.
Because most parallel robots have no large deformation during operation, it is assumed that there is no coupling between rigid body displacement and elastic deformation in the calculation. This assumption also ensures that the end stiffness changes caused by joint clearance can be calculated independently. Usually, the R pair consists of two parts: the rotating shaft and the outer cylinder of the clearance. It is assumed that there is an assembly clearance between the rotating shaft and the rod in the axial and radial directions of the rotating shaft. A local coordinate system o − xyz is established. The axis y is defined as the axis direction of the R auxiliary shaft. The maximum radial clearance is c 1 and the maximum axial clearance is c 2 . The length of the outer cylinder of the rod is l R and the diameter of the inner wall of the outer cylinder of the rod is d R . When the joint is subjected to pure axial and radial constraints, as shown in Figure 7, an offset of the maximum clearance c 1 , c 2 , occurs in the direction of the constraint force.  When the constraint force is the pure constraint couple around the y-axis, as shown in Figure 8. This axis is in the direction of rotational degrees of freedom and will not undergo angular deformation; When the constraint force is a pure constraint moment around the x-axis or z-axis, it can be inferred from the geometric relationship of the clearance model Solving the above equation, we can obtain When the constraint force is the pure constraint couple around the y-axis, as shown in Figure 8. This axis is in the direction of rotational degrees of freedom and will not undergo angular deformation; When the constraint force is a pure constraint moment around the x-axis or z-axis, it can be inferred from the geometric relationship of the clearance model Without considering the coupling effect of the constraint force and constraint force couple, the simplified approximate clearance model of the R pair can be expressed as Solving the above equation, we can obtain Without considering the coupling effect of the constraint force and constraint force couple, the simplified approximate clearance model of the R pair can be expressed as The UPR branch is subject to the constraint couple m c in the U pair normal direction, the constraint force f c passing through the U pair center and the direction parallel to the R pair axis direction, and the driving force f a along the P pair direction. The deflection of the terminal of the UPR branches when they are individually constrained by forces/couples was analyzed separately. Figure 9 shows the movement of the U pair with clearance when it is subjected to the driving force along the rod. The clearance model of the U pair can be simplified as a superposition of two revolute pair clearance models with orthogonal axes. The U pair coordinate system is established in Figure 9. Taking the center of the U axis as the coordinate origin and establishing two local coordinate systems o − u v w and o − uvw through the right-hand rule, the axis u and the axis v are along the two orthogonal axis directions of the U axis, respectively, and the coordinate system o − uvw is obtained by rotating the angle α around the axis v of the coordinate system o − u v w . The axis w direction of the coordinate system o − uvw is the same as the P pair direction in the UPR branch. Without considering the coupling effect of the constraint force and constraint force couple, the simplified approximate clearance model of the R pair can be expressed as The UPR branch is subject to the constraint couple c m in the U pair normal direction, the constraint force c f passing through the U pair center and the direction parallel to the R pair axis direction, and the driving force a f along the P pair direction. The deflection of the terminal of the UPR branches when they are individually constrained by forces/couples was analyzed separately. Figure 9 shows the movement of the U pair with clearance when it is subjected to the driving force along the rod. The clearance model of the U pair can be simplified as a superposition of two revolute pair clearance models with orthogonal axes. The U pair coordinate system is established in Figure 9. Taking the center of the U axis as the coordinate origin and establishing two local coordinate systems Assuming that the radial and axial clearances of the two orthogonal R pairs in the U pair are c 1 , c 2 respectively, when the UPR branch is only subjected to the driving force along the direction of the P pair, the rotating shaft where the U pair and the P pair are connected will generate a clearance c 1 along the radial axis w direction of the rotating shaft and a contact force f u along this direction. The spatial motion of this axis is expressed in matrix form as ∆U The first lines of three represent the offsets along the three axes of the coordinate system o − u v w , and the last lines of three represent the deflections around the three axes.
Then, the offset/deflection of the U pair in the coordinate system o − uvw under the action of the driving force f a is: In the Equation, o o R is the rotation transformation matrix of the coordinate system o − u v w relative to the coordinate system o − uvw, and, at this time, The offset ∆U f a w of the UPR branch U pair along the direction of the driving force f a and under the action of the driving force f a can be obtained.
The R pair is acted on by the driving force f a , and the offset ∆R f a w in the direction of the driving force is radial clearance c 1 .
Therefore, the total offset of the UPR branch in the direction of the driving force produced by the driving force f a is The UPR branch is subjected to a constraint couple along the normal direction of the U axis. The deflection generated by the U pair can be equivalent to the superposition of the deflections of the two R pairs with mutually orthogonal axes around the direction of the constraint couple; that is, ∆U m c w = 2m θ . The deflection generated by the R pair is the deflection ∆R m c w , which is around the direction of the constraint couple, ∆R m c w = m θ . Then, the total deflection of the UPR branch along the direction of the constraint couple m c produced by the constraint couple m c is The RPU branch clearance model is the same and will not be introduced here.
The RR intermediate branch clearance model can be equivalent to a U pair, and the branch terminal offset/deflection is We can obtain The method of the virtual work principle is suitable for the prediction of the position and attitude accuracy of the mechanism under complex working conditions, and it can achieve a more reasonable description between the joint clearance and the end position and attitude error. The virtual work principle is used to determine the relationship between the position error of each branch and the position error of the platform caused by the clearance. The superposition method is used to quantify the position and attitude error of the branch caused by the clearance of each joint. It is assumed that the offset/deflection of the end random numbers. It includes 16 random variables related to joint clearance. Equation (60) is amended as follows: (63) In the Equation, C m×N = a ⊗ B m×N is defined as multiplying the real number a by each element in the matrix B to obtain the matrix C.
The contact force is approximately equal to the constraint force and is much larger than the friction.
Contact deformation is linear elasticity, regardless of the coefficient of restitution, contact deformation speed and initial impact speed.
The contact force of the joint is approximately equal to the constraint reaction force. The influence of the coupling of the constraint force and constraint couple on the contact deformation is ignored. The contact deformation under the independent action of the restraint force/moment is calculated separately, and then the calculated results are superposed.
When there is clearance, the restraint force of each branch is approximately equal to that obtained by solving the stiffness matrix of the mechanism when there is no clearance.
There are three cases of the R-pair contact model: (a) plane-plane contact along the axis of the shaft caused by axial restraint; (b) cylinder-cylinder contact along the radial direction of the shaft caused by radial restraint; (c) line-plane contact caused by the radial restraint couple.
The kinematics pair has two parallel contact surfaces. The contact areas of the contact surface are known. The contact force model is expressed as The corrected elastic modulus can be expressed as In the Equation, the elastic modulus and Poisson's ratio of one contact body are E 1 and µ 1 , and the elastic modulus and Poisson's ratio of the other contact body are E 2 and µ 2 .
It can be obtained that the contact deformation caused by the axial constraint force is Two cylinders whose axes are parallel to each other are in contact, and the contact force can be expressed as The contact deformation caused by the radial constraint force is For the contact deformation caused by the constraint couple in the radial direction, the contact force can be expressed as The functional relationship between the contact force and the deformation at the contact point is R is the equivalent radius of the two contact bodies, and is shown as follows: It can be seen from the above Equation that the functional relationship between the contact angle deformation of the joint and the constraint couple is not linear. In order to simplify the contact deformation model, the constraint couple within a certain range of the contact deformation model is universal. A correction factor κ =τ 1/3 c is introduced to linearize the function of the contact angle deformation as a function of the constraining couple.
The R pair approximate contact deformation model can be defined as In the Equation, f c R is the six-dimensional constraint matrix of the R pair. The U pair joint contact deformation model can be equivalent to the U pair joint clearance model, and the approximate contact deformation model in the U-axis coordinate system o − u v w can be defined as = diag k con a + k con b k con a + k con b 2k con b 0 0 2k con In the Equation, f c o is the six-dimensional constraint matrix in the U-axis coordinate system o − u v w .
Assuming that each clearance joint of the mechanism reaches a permanent contact state under a static large load, the contact force is approximately equal to the constraint force. Taking the UPR branch as an example, the contact deformation of the driving pair is ignored. Firstly, the contact deformation caused by the driving force f 11 along the direction of the driving pair is analyzed. The driving pair is regarded as a two-force rod. The R pair and U pair are connected with the rotating shaft of the driving pair to produce radial contact deformation, and the rotating shaft that the U pair and the fixed platform are connected to produce radial contact and axial contact.
In order to analyze the contact deformation of U pair axis 1 in the direction of the driving force f 11 , the constraint force f 11 is transformed from the branch coordinate system o − uvw to the U axis coordinate system o − u v w .
In the Equation, According to Equation (77), the six-dimensional contact deformation of the U pair shaft 1 in the coordinate system o − u v w under the action of the driving force f 11 is The contact deformation ∆R con 1 generated by the U pair axis 1 is transformed into the branched coordinate system o − uvw.
Therefore, the total contact deformation of the U pair in the direction of the driving force f 11 is the superposition of the deformation ∆U con 1 (3, 1) of the shaft 1 in the direction of the driving force f 11 and the deformation ∆U con 2 (3, 1) of the shaft 2 in the direction of the driving force f 11 .
When the UPR branch is only subjected to the driving force f 11 , the R pair produces radial contact deformation ∆R con f 11 in the direction of the driving force f 11 The contact offset of the UPR branch in the direction of the driving force caused by the driving force f 11 is ∆UPR con f 11 = ∆U con f 11 +∆R con f 11 = f 11 k con a sin 2 α+ f 11 k con b cos 2 α+2 f 11 k con b = ∆KUPR con The UPR branch is subject to a constrained couple m 11 along the direction of the U-axis normal. The contact deflection of the U pair produced by the constraint couple m 11 can be equivalent to the superposition of the angular contact deformation of the two R pairs whose axes are orthogonal to each other around the direction of the constraint couple; that is, ∆U con m 11 = 2θ con c . The contact deflection of the R pair produced by the constraint couple m 11 is ∆R con m 11 = θ con c . Then, the total angular contact deformation of the UPR branch along the direction of the constraint couple m 11 produced by the constraining couple m 11 is ∆UPR con m 11 = ∆U con m 11 +∆R con m 11 = 3θ con c = ∆KUPR con m 11 m 11 (83) The UPR branch is subjected to a constraint force f 12 passing through the U pair center and in a direction parallel to the R pair axis. Among them, the U pair produces the contact deformation d con a in the axial direction of the rotating shaft parallel to the R pair and the contact deformation d con b in the radial direction of the rotating shaft orthogonal to the R pair. The R pair produces contact deformation d con a in the axial direction. Then, the contact offset of the UPR branch along the direction of the constraint force f 12 produced by the constraint force f 12 is ∆UPR con f 12 = ∆U con f 12 +∆R con f 12 = 2d con a + d con b = ∆KUPR con The RPU branch clearance model is the same and will not be introduced here. mean value. Therefore, the modified random physical variable considering joint contact deformation can be expressed as the multiplication of initial physical parameters and random numbers. It includes 16 random variables related to joint contact. The correction of Equation (90) is as follows: K con na = ς N ⊗ K con na ; K con ra = ς N ⊗ K con ra ; K con nc = ς N ⊗ K con nc ; K con rc = ς N ⊗ K con rc ; (91)

Global Stiffness Analysis of Ideal Joints
In order to reflect the end deformation of a 2UPR-RR-2RPU redundant PM when driving stiffness and branch deformation are considered, the square root of linear deformation in three directions is taken as the stiffness evaluation index.
In the Equation, ∆X The structural parameters of the 2UPR-RR-2RPU redundant PM are as follows: the external diameter of the swing rod is 40 mm; the internal diameter of the swing rod and the external diameter of the push rod are 25 mm; the length of the swing rod is 205 mm; the diameter of the intermediate restraining chain is 45 mm and its height is 320 mm; the radius of the fixed platform is 260 mm; the radius of the moving platform is 210 mm; the elastic modulus is 200 GPa; and the Poisson's ratio is 0.3. In the analysis of joint clearance, it is assumed that the axial clearance of the joint is 0.1 mm and the radial clearance is 0.05 mm. In the analysis of joint contact deformation, the length of the joint outer cylinder is 105 mm, and the diameter of the joint outer cylinder inner wall is 25 mm given a 7500 N external load downward along the normal direction of the moving platform. This parallel mechanism is applied in the field of wheel coupling fatigue durability testing. The mechanism is connected in series with the z-direction actuator cylinder to simulate the force acting on tires on rough, potholed and inclined roads. It compensates for the deficiency of a four-channel road simulation platform that can only provide z-direction excitation. The quarter vehicle load is 7500. Through the above stiffness performance evaluation indicators, the global performance distribution of a 2UPR-RR-2RPU redundant PM under four factors was studied, respectively, as shown in Figure 10.
Through the defined stiffness index in this paper, the following conclusions can be drawn.
After the introduction of redundant branches, the stiffness characteristic distribution of the PM is improved. The structural stiffness has good symmetry in the workspace. The linear stiffness and angular stiffness are obvious. The stiffness of the mechanism is less affected by the rotation angle β of the moving platform and the platform is more stable when the attitude angle is output. Through the defined stiffness index in this paper, the following conclusions can be drawn.
After the introduction of redundant branches, the stiffness characteristic distribution of the PM is improved. The structural stiffness has good symmetry in the workspace. The linear stiffness and angular stiffness are obvious. The stiffness of the mechanism is less affected by the rotation angle β of the moving platform and the platform is more stable when the attitude angle is output.

Comparison between Theoretical Model and Finite Element Model of Static Stiffness for Ideal Joint
In order to verify the correctness of the theoretical model under the ideal joint, the Solidworks model was imported into the Ansys Workbench. The material performance parameters were defined as follows: density ρ = 3 7850kg / m , Poisson's ratio μ = 0 .3 , elastic modulus  E= 2 11 10 Pa . The theoretical model and finite element model of the mechanism under two positions and attitudes were established for comparison. The first attitude: The joint was set as an ideal kinematics pair, and fixed constraints were added to the fixed platform. The moving platform was set as rigid, and the following six-dimensional load was applied at the center point of the moving platform coordinate system: The deformation cloud diagram obtained by finite element simulation is shown in Figure 11, and the comparison between theoretical and simulation results is shown in Table 1.

Comparison between Theoretical Model and Finite Element Model of Static Stiffness for Ideal Joint
In order to verify the correctness of the theoretical model under the ideal joint, the Solidworks model was imported into the Ansys Workbench. The material performance parameters were defined as follows: density ρ = 7850 kg/m 3 , Poisson's ratio µ = 0.3, elastic modulus E = 2 × 10 11 Pa. The theoretical model and finite element model of the mechanism under two positions and attitudes were established for comparison. The first attitude: α = 0, β = 0; The second attitude: α = 15 • , β = 15 • . The joint was set as an ideal kinematics pair, and fixed constraints were added to the fixed platform. The moving platform was set as rigid, and the following six-dimensional load was applied at the center point of the moving platform coordinate system: The deformation cloud diagram obtained by finite element simulation is shown in Figure 11, and the comparison between theoretical and simulation results is shown in Table 1.

Model Correction Considering Hinge Clearance and Hinge Contact Deformation Stiffness
The given load of the moving platform is as follows: The revised stiffness model considering the clearance and joint contact deformation was calculated 5000 times. The 32 random variables involved in the calculation each time meet the normal distribution. A histogram was drawn and the distribution was fitted for the 5000 groups of final end deformations-that is, the distribution shape of the end stiffness-taking 50 sample spaces within the maximum and minimum range of linear deformation/angular deformation at the end of the 2UPR-RR-2RPU redundant PM. The stiffness distribution of the initial position and attitude considering joint clearance and joint contact deformation is shown in Figures 12 and 13.   According to the stiffness distribution diagram, the following conclusions are drawn. The fluctuation range and probability distribution of the output end stiffness caused by the random parameters related to 16 joint gaps and 16 joint contacts are clearly shown. According to the law of large numbers, although the linear stiffness/angular stiffness distribution in each direction is slightly different, it basically conforms to the normal distribution.

Comparison between Theoretical Model and Finite Element Model Considering Joint
The shape of the stiffness distribution is determined by the predetermined pose and given parameters (structure size, clearance value, contact parameters, load and variance).

Comparison between Theoretical Model and Finite Element Model Considering Joint Clearance and Joint Contact Deformation
In order to facilitate the convergence of the model, the calculation time was saved. In the simplified model, the connecting rod was set as a rigid body, and the clearance joint was set as frictionless contact. The mesh of the contact area was refined, the contact parameters and iteration steps were set reasonably and the fixed constraints were added to the fixed platform. The deformation of the moving platform with clearance joints under two different loads in the attitude The above load on the moving platform was applied and the contact parameters of the clearance joint (U pair, R pair) were set. The angular deformation around the x, y, z axis tends to zero under the given boundary conditions. Therefore, only the line deformation along the x, y, z direction was analyzed. When the standard deviation is 1/16, the finite element simulation of the stiffness prediction interval at the output end of the attitude and the center of the moving platform is obtained. The simulation results are shown

Comparison between Theoretical Model and Finite Element Model Considering Joint Clearance and Joint Contact Deformation
In order to facilitate the convergence of the model, the calculation time was saved. In the simplified model, the connecting rod was set as a rigid body, and the clearance joint was set as frictionless contact. The mesh of the contact area was refined, the contact parameters and iteration steps were set reasonably and the fixed constraints were added to the fixed platform. The deformation of the moving platform with clearance joints under two different loads in the attitude α = 0, β = 0 is discussed. The above load on the moving platform was applied and the contact parameters of the clearance joint (U pair, R pair) were set. The angular deformation around the x, y, z axis tends to zero under the given boundary conditions. Therefore, only the line deformation along the x, y, z direction was analyzed. When the standard deviation is 1/16, the finite element simulation of the stiffness prediction interval at the output end of the attitude and the center of the moving platform is obtained. The simulation results are shown in Table 2, and the total deformation cloud diagram is shown in Figure 14. The deformation in each direction under two different loads is shown in Figures 15 and 16.  Table 2, and the total deformation cloud diagram is shown in Figure 14. The deformation in each direction under two different loads is shown in Figures 15 and 16.         As shown in Table 2, when considering the joint clearance and joint contact deformation, the finite element simulation results under the given standard deviation are within the stiffness variation range of the theoretical model. The maximum error of theory and simulation is 5.79% when measuring the error by sample mean. The correctness of the theoretical model was verified through a comparative analysis. The model can be applied to the static stiffness performance evaluation of a 2UPR-RR-2RPU redundant PM at the design stage, laying the foundation for stiffness optimization.

Stiffness Model Analysis of 2UPR-RR-2RPU Redundant PM Considering Multiple Factors
In order to evaluate the accuracy of the joint clearance model and joint contact deformation model under different external loads, the deformation of the end under two different external loads is discussed. The force amplitude of load 1 along the z direction is far greater than that along the x, y directions, whereas the force amplitude of load 2 along the x, y, z directions is equal. Taking t as the step size, the simulation time is given as 2 s, and the output motion law and applied dynamic load of the PM are as follows: Plan attitude and load 1: Plan attitude and load 2: Based on the analysis method proposed in this paper, the driving stiffness, bar deformation, joint clearance and joint contact deformation were considered. Using Matlab software to calculate the 2UPR-RR-2RPU redundant PM, the deformation of the PM output end in the six-DOF direction is caused by various factors under the above load conditions and given motion laws. Statistical analysis was carried out on the modified clearance model and joint contact deformation model. The stiffness performance considering clearance and joint contact deformation was evaluated by means of sample mean and the quality of the statistical simulation results was evaluated. The maximum/minimum values of the samples were used to judge the stiffness fluctuation range considering the clearance and joint contact deformation. The approximate linear/angular stiffness interval of the mechanism in each direction can be determined. The maximum, minimum and average values of line deformation/angular deformation in each direction considering different factors under two different planning attitudes and loads are shown in Figures 17 and 18.

Conclusions
Firstly, a high-precision stiffness modeling method of an over-constrained redundant PM was proposed. Taking a 2UPR-RR-2RPU over-constrained redundant PM as an example, the ideal joint was derived by using strain energy and Castigliano's second theorem. The stiffness matrix of the mechanism with its driving stiffness and branch deformation and the restraint force of each branch should be considered.
Based on the principle of virtual work, a stiffness modeling method considering the joint clearance and joint contact deformation was proposed. The probabilistic model was used to predict the position distribution of the clearance joints, and the problem of the modeling error caused by coupling constraints of over-constrained mechanisms was solved. The stiffness variation interval considering the joint clearance and joint contact deformation under different external loads was obtained, and the influence degree of each factor (driving stiffness, bar deformation, joint clearance, joint contact deformation) on the end deformation was analyzed.
Finally, the correctness of the proposed stiffness model was verified using the finite element method. This method can help robot designers and manufacturing engineers to fully consider the influence of multiple factors (driving stiffness, bar deformation, joint clearance, joint contact deformation) on the static stiffness of the mechanism. Through this method, a high-precision static stiffness model can be established to achieve the best performance of the mechanism.