Dynamic simulation analysis and experimental study of an industrial robot with novel joint reducers

Zero-backlash high-precision roller enveloping reducers (ZHPRER) possess a wide array of potential applications in the field of industrial robot because of their high precision and efficiency. This paper presents a pilot study to verify the applicability of this type of reducers in industrial robots. An industrial robot was designed such that its joint reducers were the roller enveloping reducers. A dynamic model for this robot which accounts for the dynamic responses of the reducers was also established using multibody dynamics. This dynamic model was then used for analyzing the dynamic behavior of the joint reducers and the body of the robot under different operating conditions. Simulation results yielded from the developed model were verified by comparing them with the data obtained from experiments. The present study for the first time confirms the evident advantages of ZHPRER for industrial robotic applications. Meanwhile, the proposed dynamic model for the industrial robot provides theoretical support for the subsequent design of error compensation control.


Introduction
Industrial robots as key equipment in manufacturing operations have been extensively applied in various industries such as mechanical engineering, healthcare, and aerospace. With recent advances in manufacturing processes, expectations to the industrial robots are becoming higher. Future robots are expected to have advantages such as high precision, high load bearing capacity, high efficiency, and low noise [1]. One component that critically influences the performance of an industrial robot is its joint reducer [2][3][4][5], whose precision determines the repeatability of the robot, whose output torque decides the robot's load bearing capacity, and whose efficiency regulates the power requirement of the robot's servo motor.
Y. Liu 1 The joint reducers used in industrial robots can be divided into two categories: rotary vector (RV) reducer and harmonic reducer [6,7]. Both reducers have high reduction ratio and high accuracy while the RV reducer has a higher load bearing capacity and is usually placed in the position of heavy load such as base, big arm, and shoulder. Such reducers are mainly used for the load-bearing body parts such as the trunk and arms of the robot. The harmonic drive is more compact with a simpler structure, and is usually placed in the small arm, wrist, or hand. These two reducers have attracted much attention in recent years and some works have been done to improve the performance of those reducers, therefore enhancing the overall quality of the industrial robots.
For the RV reducer, Yang et al. [8] proposed a modified advanced mean value method to optimize the design parameters of RV reducers to reduce their size while increasing their reliability. Wang et al. [9] presented a design method for slicing cutter aiming at the circular arc tooth in pinwheel housing of an RV reducer. The cutting experiment showed that this design method the requirement of the workpiece for the machining precision. Yang et al. [10] established a mathematical model for the analysis of the single secondary internal gear pairs with the tooth surface contact of an involute RV reducer based on a finite element linear programming method. Yang's model was verified by comparing the calculation results with the results obtained from finite element analysis. Other researchers, such as Yang [11], Li [12], and Wu [13], demonstrated how manufacturing errors, such as machining and assembly errors, influence the transmission accuracy of the RV reducer. The outcomes of those works provide theoretical support for improving the transmission accuracy of the RV reducer.
For the harmonic reducer, Hu et al. [14] developed a nonlinear torsional dynamic model for the harmonic gear reducer and used that model to analyze the effects of various factors on the torsional vibration of the reducer. Hu's results provide theoretical guidance for the optimal design and vibration reduction of the harmonic reducer. Li et al. [15] proposed a performance margin model that takes into account multisource uncertainties and the wear between meshed tooth surfaces. Calculation results showed that the developed method can facilitate the design and manufacturing of the harmonic reducer. Song et al. [16] put forward two methods for estimating the output torque of robotic joint with harmonic reducers via calibration of its existing flexibility. The estimated torque values using these two methods were in good agreement with the results measured using a commercial torque sensor.
Although the previous studies provided valuable approaches for optimizing the joint reducers of industrial robots, the actual operating conditions of the reducers as well as robots were not taken into account in those researches. Because of that, the computational or theoretical analysis results obtained in those studies differ noticeably from the real situation [11,17]. In addition, conventional reducers cannot completely overcome the friction and wear caused by sliding friction due to the contact and meshing of gear tooth surfaces [18][19][20], which have severely limited the efficiency and accuracy of the reducers. Those two issues are the main barriers that inhibit the enhancement of the industrial robots' performance.
To overcome those barriers, Deng et al. [21][22][23][24][25] developed a zero-backlash high precision roller enveloping reducer (ZHPRER). This new reducer coverts sliding friction to rolling friction during the gear engagement with the replacement of gear teeth with rollers, thereby reducing the tooth wear and increasing transmission efficiency and accuracy. This type of reducer has a wide array of potential applications including joint reducers for industrial robots, turning centers of CNC machines, and gearboxes of high-speed railway trains [26]. While manufacturing methods and lubrication mechanism of ZHPRER have been thoroughly investigated by Deng and coworkers [27][28][29][30][31][32], the applicability of such reducers in the industrial robots remains unstudied. To promote the application of ZHPRER in the field of industrial robot, a systematic understanding of the dynamic response of such reducers in an industrial robot and the robot's dynamic performance must be obtained.
In the present study, an industrial robot was first designed, in which the ZHPRERs were used as its joint reducers; a multibody dynamics-based model was also developed for this robot. Next, the dynamic performance of the industrial robot, including the effect of the joint reducers' dynamic response, was analyzed to obtain a set of performance parameters of the robot and its reducers under different operating conditions. Afterwards, the obtained parameters were employed to make a prototype for this robot. Experimental analysis results have once again confirmed the potential applicability of ZHPRER in the field of industrial robot and provided theoretical support for the futuristic design and optimization.
The rest of this paper is structured as follows: Sect. 2 introduces a kinematic model established for the robot based on DH method; Sect. 3 describes the dynamic model created using multibody dynamics method, which takes into account the dynamic response of ZH-PRERs; Sect. 4 calculates the transmission error of the ZHPRERs and the operation error of the output flange under different operation conditions using the established dynamic and kinematic models; Sect. 5 presents a prototype of the industrial robot and confirms the viability of using ZHPRER as the robot's joint reducers based on experimental results; finally, this paper is closed by the concluding Sect. 6. Figure 1 displays the configuration of the industrial robot designed by the authors, in which the first and second axes bear the maximum load. Thus, two ZHPRERs are placed in those locations as the joint reducers. Key design parameters of the ZHPRERs were previously determined by the authors as gear ratio i = 30, central distance A = 80 mm, and the diameter of rollers D = 10 mm [22,25]. An RV reducer and a harmonic reducer are placed on the third and four axes, which bear less load during the operation of the robot.

Kinematic model
The configuration of the robot (Fig. 1) can be simplified as a connecting rod system, for which a Denavit-Hartenberg (DH) coordinate system [33] is established, as shown in Fig. 2. The unlabeled coordinate axes in that figure can be determined using the right-hand rule. As illustrated in the figure, rotations of the wrist are realized through a linkage between the connecting rod and the triangular rocker, which ensures that the rotation axis at the output end faces downwards along the vertical direction. The correlations among the three rotational angles are indicated in the figure (the green arrows point in the positive directions): as the output end rotates from position 1 to position 2, the joint at the axis rotates by an angle β, and the wrist joint rotates by −β; as the output end rotates from position 2 to position 3, the joint at the axis and the wrist joint rotate by −γ and γ . Therefore, as the output end rotates from position 1 to position 3, the wrist joint would rotate by −(β + γ ). The relationship among θ 2 , θ 3 , and θ 4 can therefore be represented as θ 4 = −(θ 2 + θ 3 ). Table 1 lists main parameters for this connecting rod system. As indicated in Fig. 2, five local coordinate systems were established for the connecting rod system, where i = "1" denotes a ZHPRER joint reducer for the first axis; "2" the second joint reducer for the second axis; "3" an RV reducer for the third axis; "4" the wrist that connects the rod and the triangular rocker; and "5" a harmonic reducer. In Table 1, a i is the distance from z i to z i+1 along the x i axis; α i is the rotation angle from z i to z i+1 around the x i axis; d i is the distance from x i−1 to x i along the z i axis; and θ i is the rotation angle from x i−1 to x i around the z i axis.
Transformation matrices i−1 i T between a reference coordinate system "i −1" to a current system "i" can be obtained as:         With the θ i values for each reducer, the position of the output end in the global coordinate system can be determined using those matrices (Eqs. (1) to (5)) via [34] 0 Using this equation and the values in Table 1, the motion trajectory of the robot can be solved and compared with experimental data to determine the orientation error on the robot end [35].

Multibody dynamics -theory and modeling
During robot operation, the load on each joint reducer comes partly from the load at the output end and partly from the self-weight of the robot. The motion of robot body parts changes the location of the center of gravity of each body part, causing the load on each joint reducer to change dynamically [36]. The dynamic change of the load would affect the transmission error of the joint reducer [37], which in turn influences the precision of the robot's output end [38]. Kinetic analysis on the joint reducers or on the body of the robot is not enough to reveal the transmission performance of the joint reducers and the error of the output flange of the robot. In this study, we follow a multibody dynamics approach and use Recurdyn software [39] to create a dynamic model for the industrial robot which considers the dynamic responses of the joint reducers, and use this model to analyze key dynamic parameters such as the errors of the joint reducers and the output flange of the robot.
The entire modeling process includes the establishment of the coordinate system, the assignment of kinematic pairs, contact modeling, and the application of loading conditions, which have been minutely explained in our previous work [40]. Thus, in this paper, we skip those modeling details while only presenting core theoretical equations for the ZHPRER.

Theoretical analysis
During gear transmission, the contact pattern of tooth surfaces has a significant influence on the reducer's dynamic performance. According to the contact algorithm in Recurdyn [41], we can calculate the normal force f n and frictional force f f between two contacting objects based on the fine penetration depth when the models are in contact [42,43], as shown in Eqs. (7) and (8): where k denotes the elastic coefficient; c is the damping coefficient; m1, m2, and m3 are the stiffness, damping, and indentation exponent, respectively; δ and δ max denote the penetration depth and maximum penetration depth, respectively; δ max is used to limit abnormal f n calculated because of excessive penetration depth. The time derivative of δ,δ, represents the penetration rate. The coefficient of friction μ(ν) can be calculated based on the coefficient of static friction μ s , coefficient of dynamic friction μ d , as well as relative speed thresholds of static friction ν s and dynamic friction ν d as: hav sin (x, x 0 , y 0 , If the correlation coefficients and exponents are appropriately selected, the relationship among the contact force of tooth surface, penetration depth, and external load can be correctly determined. Here, the penetration depth between two contacting tooth surfaces can be approximated as the microscale contact deformation of the tooth surfaces under load to calculate the transmission error of the reducer, which is essentially a worm drive. Contact parameters used in this study were set according to Gummer and Sauer [44,45], as listed in Table 2. It is worth mentioning that the parameters listed in that table will be validated in Sect. 5.

Multibody dynamics modeling
A simplified model was created for the industrial robot to achieve desired computing efficiency and accuracy. For the robot body, the components that would not influence the analysis results such as threads, lead angles, and bolts were ignored. For the joint reducers, since the focus of this study is the ZHPRERs on axes 1 and 2, we only modeled the worm drives and ignored the bearings. For the reducers on axes 3 and 4, we only considered their mass and conducted a Boolean summation of the corresponding connecting rods. Bearing clearance and assembly error were also excluded when generating the simplified model (Fig. 3). Kinematic pairs were established for the connecting parts of the robot based on their relative motions following the method used by McPhee et al. [46] and Yang et al. [47]. For example, the kinematic pairs for the ZHPRERs for axes 1 and 2 considered the relationships among the worm, worm wheel, and other connecting components, as well as the relative rotation between the rollers and the worm wheel. The kinematic pairs for the reducers of axes 3 and 4 and other rotating components were defined as revolute joints. In particular, for a connecting rod, if both its ends were set as revolute joints that would lead to redundancy, which might affect the accuracy of subsequent calculation [48]. Therefore, the connection between a connecting rod and a triangular rocker was modeled as an inline joint not to overconstrain the entire system. The defined kinematic pairs for the simplified robot model are listed in Table 3.

Multibody dynamic simulation
The performance of the ZHPRERs used in the robot can be evaluated based on their load and rotation speed through rigid body dynamics simulation using the present multibody dynamics mode (Eqs. (1)-(10)). As displayed in Table 4, seven cases were simulated to reveal the operation conditions of the ZHPRERs at different rotation speeds and under different loads. The simulation results obtained from cases 1 to 4 demonstrate the effect of rotation speed on the performance of the ZHPRERs while the results acquired from case 1 and cases 5 to 7 reveal how the load of a ZHPRER influences its practical performance. Figure 4 illustrates the movement process of the robot. This mechanism is concurrently driven by axes 1 and 2 and the worm wheels of the ZHPRERs on those axes first rotate to their maximum rotation angle (85°) and then return to their original positions. How the ZHPRERs on axes 1 and 2 jointly affect the operation of the robot is elaborated through the simulation results.
In the simulation, if the rotation speed and load are applied suddenly, the contact status between tooth surfaces will subject to abrupt changes, which will cause errors in simulation results. To eliminate this influence, the speed and load are applied gradually following the curves displayed in Figs. 5(a) and 6, respectively. As shown in Fig. 5, no external load is applied within the first second of multibody dynamic analysis to let the entire system achieve an equilibrium status. After the first second, the speed of the ZHPRER increases linearly and achieves the maximum speed when the output shaft rotates by 15°, after which the ZHPRER undergoes a uniform rotation and starts to reduce linearly when its rotation angle reaches 70°until zero when the rotation angle is 85°. Next, the system stops for 1 second to maintain equilibrium and the ZHPRER starts to rotate in reverse direction following the same steps until the worm and worm wheel return to their original positions. Different rotation speeds are applied to reveal the effect of rotation speed on the two ZHPRERs' Fig. 11 Transmission errors of the ZHPRER of axis 2 during forward rotation at different rotation speeds performance, which leads to different operation times. Figure 5(b) plots the changes of the rotation angle of worm wheel corresponding to different rotation speeds. As shown in Fig. 6, the external load increases linearly to the maximum load within the first 0.5 seconds and remains constant for the remaining time. As the ZHPRER undergoes a forward rotation and a reverse rotation during the simulation, the simulation results are discussed based on these two phases. The dynamic processes were simulated on professional graphic workstations (graphics card: NVDIA Tesla K80, processor: Intel Core i5-12400, and the maximum RAM capacity was 16 GB). Each simulation took about six hours to complete.

Precision analysis of mechanical transmission of the ZHPRER on axis 1
Transmission errors of the ZHPRER of axis 1 in forward and reverse rotations at different rotation speeds are plotted in Figs. 7 and 8, respectively. Figure 7 shows that in the forward rotation phase, the transmission error increases as the speed increases. A relatively large error occurs during the start and stopping phase because of the abrupt acceleration/deceleration from a smooth operation (uniform rotation) in those phases, which would cause sudden changes in the contact status of the ZHPRER's worm drive. Thus, in practice, the acceleration and deceleration for the start and stopping phases should be as low as possible. Figure 8 indicates that in the reverse rotation phase, the transmission error also increases with the rotation speed. However, after the ZHPRER stops at the end of the forward rotation, the robot body would experience a small amplitude vibration because of inertia. The time interval of 1 second between the forward and reverse rotation could only allow the robot body to fully stop when the rotation speed is 200, 400, and 600 rpm, while the vibration when the rotation speed is 800 rpm or higher would not fully disappear within that 1 s time interval. In this case, the application of a speed along the reverse direction will result in Fig. 12 Transmission errors of the ZHPRER of axis 2 during reverse rotation at different rotation speeds robot body vibration with larger amplitudes, which will reduce the precision of the ZHPRER during the reverse rotation phase. Therefore, it is suggested that when the rotation speed is high (800 rpm or higher), the ZHPRER should stop for a longer time (>1 s) at the end of the forward rotation phase to eliminate the influence of the robot body vibration on the precision of the ZHPRER in the reverse rotation phase. Figures 9 and 10 reveal how the external load influences the transmission error of the same ZHPRER in both forward and reverse rotations, respectively. From these figures it can be seen that the transmission errors under different loads and in both rotation phases share similar features. All the errors experience large oscillations on the start and stopping stage and the errors occurred in the reverse rotation phase are larger than those in the forward rotation phase in general. The main reason for this phenomenon is that the variation of external load will not cause extra torque on the ZHPRER of axis 1, therefore the transmission errors under different loads appear to be consistent. Since the rotation speed was set as 800 rpm in those simulations, the errors occurred in the reverse rotation phase were higher than those in the forward rotation phase because of the robot body vibration at the beginning of that phase due to the high rotation speed, as explained in the last paragraph.
In summary, the transmission error of the ZHPRER on axis 1 is mainly influenced by rotation speed. The error increases with the speed and the error occurred in the reverse rotation phase is larger than that in the forward rotation phase. Under all the simulation cases, the error during an entire operation period varies within −1 to 1 , which confirms the high precision of the ZHPRER.  bears the weight of the robot parts beyond its upper arm, therefore the rotation speed only has a trivial effect on the ZHPRER's transmission and under all the simulation cases, the transmission error is within −0.5 to 0.5 . However, the transmission error undergoes highfrequency oscillations at low rotation speeds and the frequency of oscillation decreases as the speed increases. This phenomenon reveals that the influence of the change in the contact status between worm and worm gear tooth surfaces on the transmission error is more evident at low speeds than that at high speeds.

Precision analysis of mechanical transmission of the ZHPRER on axis 2
The features of transmission error of the ZHPRER on axis 2 are similar to those on axis 1, with the value of the error varies between −0.5 to 0.5 , as reflected in Fig. 12.
The transmission errors of the ZHPRER on axis 2 during forward and reverse rotations under different loads are plotted in Figs. 13 and 14, respectively. As shown in these figures, in both forward and reverse rotation phases, the transmission error does not change considerably as the external load increases. Unlike the ZHPRER on axis 1, the external load on the robot's output flange has direct influence on the magnitude of torque applied on the ZHPRER of axis 2. As listed in Table 4, in solving the transmission error under different external loads, the rotation speed was set as 800 rpm. The simulation results manifest that the effect of rotation speed on transmission error is stronger than that of external load, and the effect of external load on the transmission error cannot be identified from those figures.
From Figs. 11 to 14, we can summarize that the transmission error of the ZHPRER on axis 2 is considerably affected by rotation speed. The transmission error increases as the speed increases but does not evidently change as the external load increases. However, in all the simulation cases, the transmission error during the entire operation period is maintained between −0.5 and 0.5 , which again confirms the precision of ZHPRER.

Position precision analysis of robot body
The movement trajectory of the output flange of the industrial robot was calculated using the kinematic model developed in Sect. 2 and determined from the multibody dynamics simulation separately. The two trajectories were then compared with each other to decide the position error of the robot's output flange (Figs. 15 and 16). This position error is caused by the transmission errors of the ZHPRERs on axes 1 and 2. Figures 15 and 16 illustrates the variations of the position error in both forward and reverse rotations at different rotation speeds. In the forward rotation, the position error increases with the rotation angle of the worm gear. Meanwhile, as the rotation speed increases, the transmission error curves transition from high-frequency, narrow-range oscillations to low-frequency, large-range oscillations, leading to sudden changes in the position error on the stopping stage. This result is in alignment with the results of transmission error analysis (Sects. 4.1 and 4.2). Therefore, within the allowable conditions, low speeds should be chosen for those ZHPRERs to achieve high position precision of the output flange.
In the reverse rotation, as speed increases, the position error curves gradually change from high-frequency, narrow-range oscillations to low-frequency, large-range oscillations, which is similar to the forward rotation case. However, different from the forward rotation case, the position error decreases at large as the speed increases during the reverse rotation. The reason is that in the forward rotation phase, as the worm wheels of the ZHPRERs rotate, the angles between the external load and robot weight and the equivalent level arm increase, thereby leading to increased torque on the ZHPRER of axis 2. The maximum torque is achieved when the worm wheel rotates by 85°. On the contrary, those angles decrease in the reverse rotation phase and the maximum torque appears in the start phase. Therefore, the transmission error would increase during the forward rotation and decrease during the In the forward rotation, the maximum position error occurs in the stopping phase while in the reverse rotation, it occurs in the start phase. Likewise, this phenomenon is caused by the change of the torque on the ZHPRER of axis 2. Also, as mentioned before, the position error curves obtained for different external loads are quite similar to each other, this is because the influence of rotation speed on the transmission error of the ZHPRERs is much larger than that of external load.
In particular, it is noted from Figs. 15 to 18 that under all simulation conditions, the position error did not start from zero in the forward rotation. The initial error is generated because the actual position of the output flange deviates its theoretical position due to the influence of the robot body weight and external load. Moreover, the influence of servo motor control was not considered in those simulations. Therefore, in practice, an error compensation algorithm is needed to eliminate the effect of that initial error and improve the precision of the industrial robot. Moreover, as shown in the above figures, the final position errors of the output flange obtained from all the simulations are within 2.5 mm, which is desirable. However, during the process, the position error would undergo large fluctuations at high rotation speeds. Thus, low rotation speeds are recommended to smooth the position error during the operation of the robot.

Experimental analysis
The experimental study consists of two steps. At first, the accuracy of the contact parameters of the multibody dynamics model was verified to validate the reliability of the multibody dynamics simulation results. Secondly, a prototype of the industrial robot was manufactured and used for experimental analysis to decide its performance. The experimental results further confirmed the application potential of ZHPRERs in industrial robots.

Verification of contact parameters
The reliability of the multibody dynamics model mainly depends on the accuracy of the transmission error analysis for the ZHPRERs, whose transmission errors are determined by the contact parameters listed in Table 2. Here we designed and built a test platform for the performance assessment of the ZHPRERs on axes 1 and 2 based on their design parameters. As displayed in Fig. 19, circular gratings at the input and output end are used to measure the transmission error of the ZHPRERs during operation.
According to GB/T 35089-2018, a standard for performance test of precision gear transmission for robot [49], the rotation speed of the worm was chosen as 10 rpm and the load on the worm wheel was set as 100 N · m. The measured transmission error of the ZHPRER under that test condition is displayed in Fig. 20(a). Comparing the measured error with the simulation result ( Fig. 20(b)), the simulation results show a better regularity. As shown in  Fig. 20(b), at every instant when the worm wheel meshes in and meshes out, the transmission error would experience a sudden change while varying within ±0.2 during a mesh period. The experimental result is more complicated: in addition to the sudden changes occurred when the worm wheel meshes in and out, the transmission error also fluctuates more fiercely during normal meshing process. The difference between the measured and calculated error is caused by the error in machining the tooth surfaces and the error in assembling the ZHPRERs. However, the measured instantaneous error signals when the worm wheel meshes in and out are in consistent with the simulated error signals. Both the measured and calculated instantaneous error signals vary between −0.4 and 0.6 . The comparison result confirms the accuracy of the simulation results and the multibody dynamics model. Moreover, from both experimental and simulation results we can find that during a rotation period, the difference between the maximum and minimum error is only 1.24 , which proves that the ZHPRERs have decent precision.

Performance assessment of robot
A prototype of the presented industrial robot was built based on the design parameters (Fig. 21), whose joint reducers on axes 1 and 2 are ZHPRERs, joint reducer on axis 3 is an RV reducer, and that on axis 4 is a harmonic reducer.
To focus on the practical performance of the ZHPRERs in the robot, we locked axes 3 and 4 and only drove axes 1 and 2, on which the two ZHPRERs are located. The procedure includes the following steps: (1) Drive the two axes and let them rotate to specified locations   During the test process, it was found that if the rotation speeds of axes 1 and 2 are too high, the output flange will experience a large-amplitude vibration in the stopping phase due to inertia, which would cause inaccuracy in the micrometer data. Therefore, the rotation speed used for this test is only set as 80 and 240 rpm (10% and 30% of the maximum rotation speed 800 rpm, respectively). At each speed, the robot was repeatedly run 10 times and the micrometer data was recorded after each run. Those data reflect the position error of the output flange, as explained above. Figure 23 displays the position errors of the output flange when the axes 1 and 2 rotated at 80 and 240 rpm. As shown in Fig. 22(a), at 80 rpm, the position error is maintained in a small range, varying around 0.002 mm. Meanwhile, the variation of the errors measured after each run is quite small. This indicates that the ZHPRER applied in the industrial robot still maintains high precision. Figure 22(b) shows that the position error undergoes comparatively large fluctuations when the rotation speed is 240 rpm, which may be caused by the high inertia as the speed increases. Nevertheless, the maximum error is only 0.068 mm, which is still low. If we continue to increase the rotation speed, due to the influence of the meshing rigidity of the ZHPRERs and overall rigidity of the robot the amplitude of vibration of the output flange in the stopping phase will be too large for us to measure the position error using the micrometer. Therefore, in the future, measures need to be taken to enhance the meshing rigidity and the robot's structure stability so that the position error can be measured at high rotation speeds.
It can be found from the experimental results that the robot has high transmission precision (with the position error around 0.002 mm) when the rotation speed is 10% of the rated Fig. 21 A prototype of the industry robot Fig. 22 Measuring position error using a micrometer speed. As the speed increases to 30% of the rated speed, the position error renders irregular, relatively large fluctuations due to the influence of inertia but the maximum amplitude is only 0.068 mm. In summary, the ZHPRERs discussed in this paper are promising candidates as joint reducers for industrial robots.
Comparing the experimentally measured position error with the simulation results (Figs. 15 to 18), we found that the error measured from the experiment is far less than the position error calculated from the simulations. This is because the real robot is equipped  Due to their high precision, high efficiency, and long lifespan, the ZHPRERs have started to be used in CNC rotary tables for machine centers. Figure 24 shows a high-precision rotary table designed by the authors, in which the ZHPRERs are used as its reducers. In this study, for the first time we built a robot using the ZHPRERs as its joints and tested its performance. It is expected that in the future, the ZHPRERs will be applied in more industrial robots.

Conclusion
The applicability of ZHPRER to industrial robots is discussed in this paper. In this study, an industrial robot was manufactured, which used ZHPRERs as its joint reducers. A kinematic model and a dynamic model for the robot were established using the DH method and the multibody dynamics method, respectively. Simulations were then performed based on those models to evaluate the transmission precision of the joint reducers and the position error of the robot's output flange at different rotation speeds and under different external loads. Finally, an experimental analysis was conducted on a robot prototype to further confirm the applicability of the ZHPRER to industrial robots. Following conclusions have been drawn from this study: (1) The transmission error of the ZHPRER on axis 1 does increase as its speed increases due to the impact and impact-induced vibration caused by high speed. However, the transmission error does not change much as the external load goes up because the increase in the external load does not cause any extra torque on the ZHPRER. (2) For the ZHPRER on axis 2, at low speeds, its transmission error undergoes highfrequency, narrow-range oscillations; such oscillations then become low-frequency, large-range at high speeds. Like the ZHPRER on axis 1, the effect of external load on its transmission error is trivial and no obvious changes on the transmission error of that ZHPRER could be identified by raising the external load. (3) Because of the influence of the transmission error of the ZHPRERs on axes 1 and 2, the position error of the output flange would experience considerable fluctuations in the start and stop phase at high rotation speeds. Increasing the external load does not obviously affect this error. (4) The simulation results were validated through the experimental analysis. By testing the industrial robot prototype, the high precision of the ZHPRERs used in the robot has been verified. However, when the rotation speed is high, the output flange of the robot will experience a relatively large-amplitude vibration in the stopping phase due to inertia.
Future work will focus on three directions: (1) Optimizing the multibody dynamics model to take into account the influence of servo motor and error compensation. The upgraded model should be able to predict the performance of a servo-controlled robot. (2) Improving the meshing rigidity of the ZHPRERs and overall rigidity of the robot to reduce the influence of rotation speed on the transmission error. (3) Replacing the contact measurement (through a micrometer in this study) with the touchless measurement and develop corresponding devices to accurately measure the position error of robots at high speeds. (4) If a laser locator is available in the future, we will measure the transmission error of the ZH-PRERs and compare the measured results with those from the simulation to further verify the advantages of the present joint reducers and the developed multibody dynamics model.