A Semi-Analytical Load Distribution Model for Cycloid Drives with Tooth Proﬁle and Longitudinal Modiﬁcations

Featured Application: In this study, a new semi-analytical load distribution model is proposed which provides a powerful tool to predict the multi-tooth contact conditions and to determine the optimal tooth proﬁle and longitudinal modiﬁcations for cycloid drives. Abstract: The current load distribution model for cycloid drives based on the Hertz contact sti ﬀ ness typically assumes a two-dimensional planar problem without considering the tooth longitudinal modiﬁcation e ﬀ ects, which fails to comply with the practical situation. In this paper, this issue is clariﬁed by developing a semi-analytical load distribution model based on a three-dimensional and linear elastic solution. Unloaded tooth contact analysis is introduced to determine the instantaneous mesh information. The tooth compliance model considering tooth contact deformation is established by combining the Boussinesq force–displacement relationships in elastic half-space with an inﬂuence coe ﬃ cient method. With this, the loads, contact patterns, and loaded transmission error are calculated by enforcing the compatibility and equilibrium conditions. Comparisons to predictions made with the assumption of Hertz contact sti ﬀ ness are presented to demonstrate the e ﬀ ectiveness of the proposed model, which shows good agreement. At the end, the e ﬀ ect of tooth longitudinal modiﬁcations on load distributions is investigated along with various loading conditions. This study yields an in-depth understanding of the multi-tooth contact characteristics of cycloid drives and provides an e ﬀ ective tool for extensive parameter sensitivity analysis and design optimization studies.


Introduction
Cycloid drives have been widely applied in many industrial areas for speed, torque conversion purposes, and precision pointing, such as in robots, machine tools, and wind power, since they were invented by Mr. Lorenz Braren in the 1930s. Figure 1 illustrates the typical structure of a single-stage cycloid drive with its main elements including two cycloid gears, several rollers, crankshaft, and output mechanism. It is a complicated eccentric-rotation-type planetary transmission mechanism, which is unlike the more commonly understood involute gearing. The small tooth number difference and multi-tooth contact characteristics give cycloid drives the advantages of compact size, light weight, large gear ratio, high positioning precision, and high load carrying capability [1,2]. Despite these benefits, as of now, many design problems of cycloid drives have not been solved, and no known design standards for cycloid drives are available [3,4]. In particular, taking the effects of both the tooth profile and longitudinal modifications into consideration, the contact strength calculation and transmission To provide some guidelines for the optimal design of cycloid drives, many research works have focused particularly on the load distribution analysis of cycloid drives, as well as further influence analysis of tooth profile modification. Malhotra and Parameswaran [6] developed a method to compute the load distributed on various components of a speed reducer based on the assumption of perfect geometry and rigid bodies of gears. Blagojevic et al. [7] presented a new two-stage cycloid speed reducer and developed a method for force analysis based on Malhotra's method. Gorla et al. [8] proposed a calculation procedure to predict the loads distributed on the elements of a new cycloid speed reducer, and an experiment was performed to validate the prediction results. Wang et al. [9] proposed a tooth contact analysis model for an RV reducer by dividing the contact area into several differential elements to demonstrate the effects of the tooth profile modification and load on the transmission precision, mechanical properties, and transmission efficiency. Xu et al. [10,11] proposed a contact dynamics model to analyze the contact dynamics of the KHV-type cycloid speed reducer while considering the effects of a turning-arm cylindrical roller. Sun et al. [12] proposed a new tooth contact analysis (TCA) method based on the discrete point tooth profile for a new type of CBR (China Bearing Reducer) reducer. By using the new TCA model, the transmission error and contact force of the designed reducer can be predicted with different modified methods. Li et al. [13,14] proposed a meshing contact analysis method for RV cycloidal-pin gear transmission considering the effect of manufacturing errors. Li et al. [5,15] presented an analytical model for the KHV-type cycloid drive considering clearances of the cycloid gear pair and output mechanism at a system level, and they also proposed a mesh stiffness calculation method for cycloid gear pairs considering the tooth profile modifications and eccentric error. Sensinger et al. [16] considered that the cycloid profile may be generated by reducing the roller position to change the tips and roots of the tooth profile and qualitatively showed the difference between the topographic profiles modified by roller position and by roller radius. Lin et al. [17] investigated cycloid tooth profile modification quantitatively under individual and combined conditions. Lin et al. [18] developed a kinematic error analysis method and tolerance design procedure for cycloidal gear reducers considering geometry, manufacturing, and precision performance. To provide some guidelines for the optimal design of cycloid drives, many research works have focused particularly on the load distribution analysis of cycloid drives, as well as further influence analysis of tooth profile modification. Malhotra and Parameswaran [6] developed a method to compute the load distributed on various components of a speed reducer based on the assumption of perfect geometry and rigid bodies of gears. Blagojevic et al. [7] presented a new two-stage cycloid speed reducer and developed a method for force analysis based on Malhotra's method. Gorla et al. [8] proposed a calculation procedure to predict the loads distributed on the elements of a new cycloid speed reducer, and an experiment was performed to validate the prediction results. Wang et al. [9] proposed a tooth contact analysis model for an RV reducer by dividing the contact area into several differential elements to demonstrate the effects of the tooth profile modification and load on the transmission precision, mechanical properties, and transmission efficiency. Xu et al. [10,11] proposed a contact dynamics model to analyze the contact dynamics of the KHV-type cycloid speed reducer while considering the effects of a turning-arm cylindrical roller. Sun et al. [12] proposed a new tooth contact analysis (TCA) method based on the discrete point tooth profile for a new type of CBR (China Bearing Reducer) reducer. By using the new TCA model, the transmission error and contact force of the designed reducer can be predicted with different modified methods. Li et al. [13,14] proposed a meshing contact analysis method for RV cycloidal-pin gear transmission considering the effect of manufacturing errors. Li et al. [5,15] presented an analytical model for the KHV-type cycloid drive considering clearances of the cycloid gear pair and output mechanism at a system level, and they also proposed a mesh stiffness calculation method for cycloid gear pairs considering the tooth profile modifications and eccentric error. Sensinger et al. [16] considered that the cycloid profile may be generated by reducing the roller position to change the tips and roots of the tooth profile and qualitatively showed the difference between the topographic profiles modified by roller position and by roller radius. Lin et al. [17] investigated cycloid tooth profile modification quantitatively under individual and combined conditions. Lin et al. [18] developed a kinematic error analysis method and tolerance design procedure for cycloidal gear reducers considering geometry, manufacturing, and precision performance.
A review of the literature reveals that the current analytical methods for cycloid drives typically assume a two-dimensional planar problem without considering the tooth longitudinal modification effects, which fails to comply with the practical situation. In fact, the cycloid gear will be crowned to relieve the stress concentration occurring on the tooth edges, which will lead to serious contact failure of the cycloid drive. Therefore, the aim of this paper is to propose a computationally efficient, semi-analytical load distribution model for cycloid drives based on a three-dimensional and linear elastic solution while considering both the tooth profile and longitudinal modifications. It will have the ability to simulate quasi-statically the multi-tooth contact conditions and determine the contact stress and pattern, loaded transmission error, and gear ratio of cycloid drives. Another purpose of this paper is to propose several types of tooth longitudinal modifications of cycloid gear teeth to improve the contact stress distribution conditions by using the proposed model.
The rest of this paper is structured as follows: Section 2 presents the tooth surface geometry of the cycloid gear, its generation method is introduced, and several types of tooth profile and longitudinal modification are proposed. Section 3 presents the semi-analytical load distribution model based on the unloaded and loaded tooth contact analyses. In this section, the tooth compliance model is established by combining the Boussinesq force-displacement relationships in elastic half-space with an influence coefficient method. Section 4 gives a numerical example to verify the proposed model by comparison with the prediction results of the current model based on the Hertz contact stiffness method. The effects of tooth profile, longitudinal modifications, and applied torque levels on the contact stress and pattern, loaded transmission error, and gear ratio are investigated by using the proposed model. Finally, conclusions are drawn in Section 5.

Tooth Surface Generation
To the authors' knowledge, there are several approaches available in the open literature for the tooth surface generation of the cycloid gear, including the methods of circle-enveloping [19,20], inverted gear train coordinate systems [12], and instant velocity center [21]. In this study, the surface equation of cycloid gears is derived based on the circle-enveloping method and homogeneous coordinate transformation.
The coordinate systems corresponding to the cycloid drive are defined based on relative motion, as shown in Figure 2. There are two movable coordinate systems, S 1 (x 1 , y 1 , z 1 ) and S 2 (x 2 , y 2 , z 2 ), and a fixed coordinate system, S f (x f , y f , z f ), attached to the pinwheel, the cycloid gear, and the housing, respectively. The surface equation of the cycloid gear tooth and the dimensionless equation of meshing can be expressed as follows: where e is the eccentricity, u is the cycloid gear width, a is the roller position radius, and ρ is the roller radius. φ 2 and φ 1 are the rotational angles having the relationship φ 2 /φ 1 = n 1 /n 2 . n 1 and n 2 are the numbers of rollers and cycloid gear teeth with one tooth number difference.

Tooth Profile Modification
Tooth profile modification (TPM) of cycloid gears is necessary to compensate for errors in assembly and manufacture and to provide better lubrication conditions [8]. Profile modification is considered to have significant effects on the backlash, transmission error, contact strength, and torque ripple of cycloid drives [5,15,22], which should be investigated qualitatively and quantitatively via parametric sensitivity analysis. The equation of the four modified profiles of the cycloid gear by modifying two basic parameters, roller position a and roller radius ρ, in individual and combined conditions [17], according to Equations (1) and (2), can be rewritten as follows: where both of the two amounts of modification ∆a and ∆ρ can be negative or positive values with different relationships of absolute magnitude, resulting in four types of tooth profile modifications, namely, positive roller radius modification (PM), negative roller position modification (NM), negative roller position and negative roller radius modification (NNM), and positive roller position and positive roller radius modification (PPM), which are listed in Table 1.

Tooth Profile Modification
Tooth profile modification (TPM) of cycloid gears is necessary to compensate for errors in assembly and manufacture and to provide better lubrication conditions [8]. Profile modification is considered to have significant effects on the backlash, transmission error, contact strength, and torque ripple of cycloid drives [5,15,22], which should be investigated qualitatively and quantitatively via parametric sensitivity analysis.
The equation of the four modified profiles of the cycloid gear by modifying two basic parameters, roller position a and roller radius ρ , in individual and combined conditions [17], according to Equations (1) and (2), can be rewritten as follows: where both of the two amounts of modification ∆a and ρ ∆ can be negative or positive values with  Figure 3 shows the four modified profiles of a cycloid gear tooth with the same tip and root relief compared with the standard tooth profile. It seems that NNM and PPM obtain the smallest and largest backlashes, respectively, while the backlashes yielded by NM and PM in sequence lie somewhere in between. Since these two individual modifications can be incorporated into the combined ones, only NNM and PPM are compared and investigated in a later section. Figure 3 shows the four modified profiles of a cycloid gear tooth with the same tip and root relief compared with the standard tooth profile. It seems that NNM and PPM obtain the smallest and largest backlashes, respectively, while the backlashes yielded by NM and PM in sequence lie somewhere in between. Since these two individual modifications can be incorporated into the combined ones, only NNM and PPM are compared and investigated in a later section.

Tooth Longitudinal Modification
Research on cycloid drives concerning the geometrical design method [20,23], kinematic principle [18,24], and load distribution analysis [4,11,12,14,15,25] have all assumed a two-dimensional planar problem without considering the three-dimensional topological structure of the cycloid gear. However, in practice, like involute cylindrical gears [26][27][28], the cycloid gear is also crowned in the longitudinal direction to avoid edge contact between the teeth and rollers, considerably improving the contact stress distribution. Three types of tooth longitudinal modifications for cycloid gears are proposed in this paper: fully crowned modification (FCM), partially crowned modification (PCM), and logarithmically crowned modification (LCM), as shown in Figure 4a-c, correspondingly. The total amount of modification along the effective contact length L of a gear tooth in the longitudinal direction is defined by the following equations.
1. In the case of FCM, where f R is the radius of a circle with origin located on the centerline of a gear tooth.
2. In the case of PCM, where p R is the radius of a circle tangent to the line with length 0 L .
3. In the case of LCM, where f is the coefficient of modification.
Tooth longitudinal modification (TLM) of cycloid gears can be carried out by a computer numerical control (CNC) form grinding machine for finishing where the grinding wheel has an accurate modified profile representing the tooth space between two adjacent cycloid gear teeth. The grinding wheel performs a radial infeed motion according to the minute amount of longitudinal modification along the gear axis to remove material from the teeth.

Tooth Longitudinal Modification
Research on cycloid drives concerning the geometrical design method [20,23], kinematic principle [18,24], and load distribution analysis [4,11,12,14,15,25] have all assumed a two-dimensional planar problem without considering the three-dimensional topological structure of the cycloid gear. However, in practice, like involute cylindrical gears [26][27][28], the cycloid gear is also crowned in the longitudinal direction to avoid edge contact between the teeth and rollers, considerably improving the contact stress distribution. Three types of tooth longitudinal modifications for cycloid gears are proposed in this paper: fully crowned modification (FCM), partially crowned modification (PCM), and logarithmically crowned modification (LCM), as shown in Figure 4a-c, correspondingly. The total amount of modification along the effective contact length L of a gear tooth in the longitudinal direction is defined by the following equations.
1. In the case of FCM, where R f is the radius of a circle with origin located on the centerline of a gear tooth. 2. In the case of PCM, where R p is the radius of a circle tangent to the line with length L 0 . 3. In the case of LCM, where f is the coefficient of modification. Tooth longitudinal modification (TLM) of cycloid gears can be carried out by a computer numerical control (CNC) form grinding machine for finishing where the grinding wheel has an accurate modified profile representing the tooth space between two adjacent cycloid gear teeth. The grinding wheel

Load Distribution Model
In this section, a new load distribution model for cycloid drives is proposed, which is treated as a three-dimensional elastic body contact problem to demonstrate the effects of both the tooth profile and longitudinal modifications accurately. This semi-analytical model is developed and analyzed under the following assumptions: 1. The load distributed on one tooth pair does not affect elastic deformations on adjacent tooth pairs. 2. The rotational speed and inertia forces of moving elements are neglected in the quasi-static process. 3. The curvature of the gear tooth surface over the contact region is assumed to remain unchanged. 4. The friction between contacting tooth pairs is neglected.
With the above assumptions, the overall methodology for the proposed semi-analytical model is shown in Figure 5. The design dimensions, modification amounts, material properties, and applied torque are the input parameters for the proposed model. The cycloid gear tooth and roller profiles described by these parameters are transformed to a global coordinate system to perform the unloaded tooth contact analysis (TCA) and to determine the backlash, unloaded transmission error, contact points, and rotational angle in a mesh circle. Next the tooth compliance model considering tooth contact deformation is developed by combining the Boussinesq force-displacement relationships in elastic half-space with an influence coefficient method. Finally, a group of equilibrium and compatibility conditions are given and solved to calculate the loads, contact stress, contact pattern, loaded transmission error (TE), and gear ratio of the cycloid drive.

Load Distribution Model
In this section, a new load distribution model for cycloid drives is proposed, which is treated as a three-dimensional elastic body contact problem to demonstrate the effects of both the tooth profile and longitudinal modifications accurately. This semi-analytical model is developed and analyzed under the following assumptions: 1.
The load distributed on one tooth pair does not affect elastic deformations on adjacent tooth pairs. 2.
The rotational speed and inertia forces of moving elements are neglected in the quasi-static process.

3.
The curvature of the gear tooth surface over the contact region is assumed to remain unchanged. 4.
The friction between contacting tooth pairs is neglected.
With the above assumptions, the overall methodology for the proposed semi-analytical model is shown in Figure 5. The design dimensions, modification amounts, material properties, and applied torque are the input parameters for the proposed model. The cycloid gear tooth and roller profiles described by these parameters are transformed to a global coordinate system to perform the unloaded tooth contact analysis (TCA) and to determine the backlash, unloaded transmission error, contact points, and rotational angle in a mesh circle. Next the tooth compliance model considering tooth contact deformation is developed by combining the Boussinesq force-displacement relationships in elastic half-space with an influence coefficient method. Finally, a group of equilibrium and compatibility conditions are given and solved to calculate the loads, contact stress, contact pattern, loaded transmission error (TE), and gear ratio of the cycloid drive.

Unloaded Tooth Contact Analysis
The unloaded tooth contact analysis (TCA) is developed in the two-dimensional plane in order to instantly determine the instantaneous central points of potential contact regions. Figure 6 shows the coordinate systems for unloaded TCA. 1 where i is the number of rollers, = [0, 0,1] k is a unit vector, and θ ri is the angle parameter of the roller profile.
Similarly, the position vector

Unloaded Tooth Contact Analysis
The unloaded tooth contact analysis (TCA) is developed in the two-dimensional plane in order to instantly determine the instantaneous central points of potential contact regions. Figure 6 shows the coordinate systems for unloaded TCA. S 1 (x 1 , y 1 ) and S 2 (x 2 , y 2 ) are connected rigidly with the pinwheel and cycloid gear, and the fixed coordinate system S f (x f , y f ) coinciding with S 1 is fixed to the housing. Hence, the position vector r (1) f (θ ri ) and its unit normal vector n (1) f (θ ri ) of all the roller profiles can be expressed in S f : where i is the number of rollers, k = [0, 0, 1] is a unit vector, and θ ri is the angle parameter of the roller profile.
Similarly, the position vector r (2) f (θ ci , φ ci , φ in ) and its unit normal vector n where φ ci is the rotational angle of the cycloid gear tooth, φ in is the input angle of the crankshaft, and θ ci is the angle parameter of the cycloid tooth profile. To perform the unloaded TCA, two conditions must be satisfied, which are coincidence of the position vectors and collinearity of the normal vectors of the modified profile of the cycloid gear and roller profiles at the same contact points [15,25]. Thus, Since for the unit normal vectors, n (2) f = 1, three nonlinear equations with the four unknowns θ ri , θ ci , φ ci , and φ in can be developed. Given one, the input crankshaft angle φ in , the other three unknown parameters can be solved.
Then, the backlash φ bli and the rotational angle φ c of the cycloid gear can be expressed as follows [5,15]: where the backlash φ bli is defined as the required rotational angle for all the potential meshing tooth pairs, and the rotational angle φ c of the cycloid gear is the maximum value among the rotational angles of all the cycloid gear teeth. where φ ci is the rotational angle of the cycloid gear tooth, φ in is the input angle of the crankshaft, and θ ci is the angle parameter of the cycloid tooth profile.
To perform the unloaded TCA, two conditions must be satisfied, which are coincidence of the position vectors and collinearity of the normal vectors of the modified profile of the cycloid gear and roller profiles at the same contact points [15,25]. Thus, (1) Since for the unit normal vectors, = = (1) 1 f f n n , three nonlinear equations with the four unknowns θ ri , θ ci , φ ci , and φ in can be developed. Given one, the input crankshaft angle φ in , the other three unknown parameters can be solved.
Then, the backlash φ bli and the rotational angle φ c of the cycloid gear can be expressed as follows [5,15]: where the backlash φ bli is defined as the required rotational angle for all the potential meshing tooth pairs, and the rotational angle φ c of the cycloid gear is the maximum value among the rotational angles of all the cycloid gear teeth.

Tooth Compliance Model
With the central points of potential contact regions defined in the unloaded TCA, the next step is to determine the tooth compliances of both contacting members. Due to the larger tooth thickness of the dedendum and the shorter tooth height of the cycloid gear, contact deformation makes a much greater contribution to the tooth deformation than bending and shear deformations. Therefore, only the contact deformation is considered, which is sufficient for engineering applications [4]. The Boussinesq solution [29,30] for the contact deformation of tooth pairs in an elastic half-space is introduced in this study.

Tooth Compliance Model
With the central points of potential contact regions defined in the unloaded TCA, the next step is to determine the tooth compliances of both contacting members. Due to the larger tooth thickness of the dedendum and the shorter tooth height of the cycloid gear, contact deformation makes a much greater contribution to the tooth deformation than bending and shear deformations. Therefore, only the contact deformation is considered, which is sufficient for engineering applications [4]. The Boussinesq solution [29,30] for the contact deformation of tooth pairs in an elastic half-space is introduced in this study. Figure 7 shows the contact state of a cycloid gear pair. The two coordinate systems S 1 (x 1 , z 1 ) and S 2 (x 2 , z 2 ) are defined with an identical origin at the initial contact point (region) from the side view, as shown in Figure 7a. The tooth surface is defined in the S(x, y) coordinate system where x and y represent the gear profile and width directions, respectively. The potential contact region is discretized into a contact grid, with each grid cell represented by a contact point, as shown in Figure 7c. In the profile direction, the contact region of a tooth pair is divided into p grid cells of identical arc length. In the gear width direction, the same potential contact region is divided into q grid cells of equal length. The area of each contact grid cell is Ω = 2ξ × 2η. The contact deformation δ c at a point (x i , y j ) can be represented as where ν and E are the Poisson's ratio and Young's modulus, and subscripts 1 and 2 refer to the roller and the cycloid gear, respectively. P(x , y ) is the distributed pressure in the contact area Ω.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 20 Figure 7 shows the contact state of a cycloid gear pair. The two coordinate systems 1 1 1 ( , ) S x z and 2 2 2 ( , ) S x z are defined with an identical origin at the initial contact point (region) from the side view, as shown in Figure 7a. The tooth surface is defined in the ( , ) S x y coordinate system where x and y represent the gear profile and width directions, respectively. The potential contact region is discretized into a contact grid, with each grid cell represented by a contact point, as shown in Figure  7c. In the profile direction, the contact region of a tooth pair is divided into p grid cells of identical arc length. In the gear width direction, the same potential contact region is divided into q grid cells of equal length. The area of each contact grid cell is ξ η ′ Ω = × 2 2 .The contact deformation δ c at a point ( , ) i j x y can be represented as where ν and E are the Poisson's ratio and Young's modulus, and subscripts 1 and 2 refer to the roller and the cycloid gear, respectively. ′ ′ ( , ) P x y is the distributed pressure in the contact area Ω .
where ρ 1 is the curvature radius of the roller, and ρ 2 is the curvature radius of the cycloid gear teeth at the contact point [5,15]. For the case of FCM, according to Equation (5), the initial separations of the cycloid gear tooth surface can be rewritten as Also, for the case of PCM, according to Equation (6), the initial separations of the cycloid gear tooth surface can be rewritten as  The initial separations of the two surfaces are defined as ω 1ij and ω 1ij along the z axis, respectively, which can be expressed as follows: where ρ 1 is the curvature radius of the roller, and ρ 2 is the curvature radius of the cycloid gear teeth at the contact point [5,15]. For the case of FCM, according to Equation (5), the initial separations of the cycloid gear tooth surface can be rewritten as Also, for the case of PCM, according to Equation (6), the initial separations of the cycloid gear tooth surface can be rewritten as where H(t) = 1, t > 0 0, t ≤ 0 and t = (2j − 1)/(q + 1) − 1 − L 0 /L.
Similarly, for the case of LCM, according to Equation (7), the initial separations of the cycloid gear tooth surface can be rewritten as Then, the governing integral equation for the three-dimensional elastic body contact of the cycloid gear pairs with a modified flexibility method can be expressed as n r=1 P r f ij = δ t − ω 1ij − ω 2ij , r = 1, 2, . . . , n; i = 1, 2, . . . , p; j = 1, 2, . . . , q, where f ij is the influence coefficient [31], which can be represented as follows: In addition to displacement requirements, two restrictions are placed on force. The normal force F between the contacting tooth pair is calculated by n r=1 P r = F 4ξη , P r ≥ 0, r = 1, 2, . . . , n, where the pressure distributed in the contact region must be non-negative. The above linear Equations (20) and (22) with the boundary condition can be solved simultaneously to obtain the contact region, contact deformation, and contact stress. Thus, the tooth compliance of the cycloid gear pairs can be represented as Note that the tooth compliance is a nonlinear function of the forces, geometries, and material properties of the cycloid tooth pairs. The values of tooth compliance and the normal force can be calculated in an iterative process, which is described in the following section.

Loaded Tooth Contact Analysis
As shown in Figure 8, the compatibility conditions should be satisfied to calculate the number of contacting tooth pairs, which can be represented as follows: contact : no contact : where ∆φ c is the increment of the rotational angle of the cycloid disc, and α i is a small angular displacement of the cycloid gear tooth.
Then, the force and moment static equilibrium equation can be expressed as where the arm of force l i can be represented geometrically as l i = −e tan θ ri cos φ in − e sin φ in + a sin iψ r + a tan θ ri cos iψ r An iterative process is repeated to find the solutions, as shown in Figure 5. Firstly, two initial values of load and angular displacement are chosen to solve the compatibility and equilibrium Equations (24) and (25) and to calculate the tooth compliance according to Equations (15)-(23) in the meantime, after the backlash is calculated by unloaded TCA. All the values are updated each time the load and angular displacement enter into the system of equations. When the variations of the force and angular displacement converge within 1%, the iterative procedure is terminated to output the loads, contact stress, contact pattern, loaded transmission error (LTE), and gear ratio (GR).
An iterative process is repeated to find the solutions, as shown in Figure 5. Firstly, two initial values of load and angular displacement are chosen to solve the compatibility and equilibrium Equations (24) and (25) and to calculate the tooth compliance according to Equations (15)- (23) in the meantime, after the backlash is calculated by unloaded TCA. All the values are updated each time the load and angular displacement enter into the system of equations. When the variations of the force and angular displacement converge within 1%, the iterative procedure is terminated to output the loads, contact stress, contact pattern, loaded transmission error (LTE), and gear ratio (GR).

Numerical Example and Verification
An example cycloid gear pair was considered to demonstrate the capabilities of the proposed model, and its basic parameters are listed in Table 2. For convenience, it was assumed that the cycloid gear and rollers are made of the same material, which is steel with Young's modulus E = 206 GPa and Poisson ratio v = 0.3. Note that the numbers of discrete elements p and q were selected for each contact zone based on the convergence time and the calculation accuracy. As the number of discrete elements increases, the computational time consumption increases proportionally, while the calculation accuracy tends to be stable, since the variations of the calculated results converge within an acceptable error. With this, the numbers of discrete elements were all determined to be 40 in this numerical example.

Numerical Example and Verification
An example cycloid gear pair was considered to demonstrate the capabilities of the proposed model, and its basic parameters are listed in Table 2. For convenience, it was assumed that the cycloid gear and rollers are made of the same material, which is steel with Young's modulus E = 206 GPa and Poisson ratio v = 0.3. Note that the numbers of discrete elements p and q were selected for each contact zone based on the convergence time and the calculation accuracy. As the number of discrete elements increases, the computational time consumption increases proportionally, while the calculation accuracy tends to be stable, since the variations of the calculated results converge within an acceptable error. With this, the numbers of discrete elements were all determined to be 40 in this numerical example.

Verification of the Proposed Model
The calculated results including contact stress, LTE, and GR are presented and compared with those obtained by the current model based on the nonlinear Hertz contact stiffness method [5,15], which was developed specifically for loaded contact analysis of cycloid drives assumed as a two-dimensional plane problem. In the absence of experimental validation, this comparison serves to verify the proposed model. Without loss of generality, three different cases are considered here. One is the standard profile without tooth profile modification. The other two cases are NNM and PPM with the same radial clearances. The corresponding TPM parameters for the cases are listed in Table 3. The applied torque was set to 350 Nm. Table 3. Tooth profile modification (TPM) parameters of three cases for the example cycloid gear pair.  Figure 9. It can be seen that the proposed model agrees well with the current model in terms of the roller numbers which the load is distributed on and magnitudes of contact stress for the three cases, although edge contact leads to increasingly concentrated stresses, which could further accord with the practical situation. Aside from this, it was also observed that, in the presence of TPM, the contact stresses increase significantly as the number of contacting tooth pairs decreases. This can be explained by the fact that the tooth profile modification causes more contacting tooth pairs to separate. Figure 10 shows the LTE and the GR of the example cycloid gear pair calculated by the current model and the proposed model for the three cases. These time histories illustrate that the same shape and period of the LTE and the GR were obtained by both models, varying periodically as the crankshaft rotates, as expected, while the average absolute value of LTE and the peak-to-peak value of the GR obtained by the proposed model are a little smaller than those obtained by the current model.

Standard NNM PPM
The results shown in  Figure 10 shows the LTE and the GR of the example cycloid gear pair calculated by the current model and the proposed model for the three cases. These time histories illustrate that the same shape and period of the LTE and the GR were obtained by both models, varying periodically as the crankshaft rotates, as expected, while the average absolute value of LTE and the peak-to-peak value of the GR obtained by the proposed model are a little smaller than those obtained by the current model.  Figure 10 shows the LTE and the GR of the example cycloid gear pair calculated by the current model and the proposed model for the three cases. These time histories illustrate that the same shape and period of the LTE and the GR were obtained by both models, varying periodically as the crankshaft rotates, as expected, while the average absolute value of LTE and the peak-to-peak value of the GR obtained by the proposed model are a little smaller than those obtained by the current model. The results shown in Figures 9 and 10 by the proposed model are in good agreement with the predictions of the current model, which verifies the effectiveness of the proposed model. From  Figures 9 and 10, we can also see the effects of different types of TPM on the load distribution and contact stress, LTE, and GR.

Effect of Tooth Longitudinal Modifications
Using the same example cycloid gear pair with NNM, four cases were considered to demonstrate the effects of TLM on the contact stress, LTE, and GR. One is the NNM gear tooth without TLM. The other three cases are FCM, PCM, and LCM, the TLM parameters of which are listed in Table 4.

Effect of Tooth Longitudinal Modifications
Using the same example cycloid gear pair with NNM, four cases were considered to demonstrate the effects of TLM on the contact stress, LTE, and GR. One is the NNM gear tooth without TLM. The other three cases are FCM, PCM, and LCM, the TLM parameters of which are listed in Table 4. Figure 11 shows the contact stress distributions of the NNM cycloid gear pair with different types of TLM at a random crankshaft angle, illustrated in 3D representations. From the results, it can be clearly observed that without TLM, the maximum stress is 1229 MPa concentrated at the two edges of the face-width of the cycloid gear tooth at a 120 • crankshaft angle, while on other mesh phases, it could get even bigger, such as 1236 MPa at a 300 • crankshaft angle. Contact stresses decrease significantly with axial distance from this edge, as shown in Figure 9a. The minimum value is about 982 MPa, occurring at the middle of the contact region. Hence, the contact pattern looks like a saddle shape, as expected. Table 4. Tooth longitudinal modification (TLM) parameters of four cases for the example cycloid gear pair with NNM.

Without TLM FCM PCM LCM
-R f = 3000 mm R p = 4500 mm, L 0 = 13.6 mm f = 0.005 mm mm, mm Figure 11 shows the contact stress distributions of the NNM cycloid gear pair with different types of TLM at a random crankshaft angle, illustrated in 3D representations. From the results, it can be clearly observed that without TLM, the maximum stress is 1229 MPa concentrated at the two edges of the facewidth of the cycloid gear tooth at a 120° crankshaft angle, while on other mesh phases, it could get even bigger, such as 1236 MPa at a 300° crankshaft angle. Contact stresses decrease significantly with axial distance from this edge, as shown in Figure 9a. The minimum value is about 982 MPa, occurring at the middle of the contact region. Hence, the contact pattern looks like a saddle shape, as expected.
As shown in Figure 11, in the presence of TLM, a dramatic reduction in concentrated stress was observed, which demonstrates the significant effects of the TLM. The maximum stresses decreased to 1045 MPa, 1075 MPa, and 1030 MPa, corresponding to the three types of TLM. Among them, the case of LCM seems to show the best contact stress distribution wherein the contact pattern presents a drum shape where the contact stress is distributed more uniformly and the concentrated stress disappears; this demonstrates the effectiveness of TLM in improving the contact stress distributions of cycloid drives. Note that the TLM parameters were arbitrarily selected for the case studies in this paper; the optimization of parameters for the three types of TLM remains as a future task.  As shown in Figure 11, in the presence of TLM, a dramatic reduction in concentrated stress was observed, which demonstrates the significant effects of the TLM. The maximum stresses decreased to 1045 MPa, 1075 MPa, and 1030 MPa, corresponding to the three types of TLM. Among them, the case of LCM seems to show the best contact stress distribution wherein the contact pattern presents a drum shape where the contact stress is distributed more uniformly and the concentrated stress disappears; this demonstrates the effectiveness of TLM in improving the contact stress distributions of cycloid drives. Note that the TLM parameters were arbitrarily selected for the case studies in this paper; the optimization of parameters for the three types of TLM remains as a future task. Figure 12 shows the effects of TLM on the LTE and GR of the NNM cycloid gear pairs in a mesh circle. It was observed that the LTE for the four cases varied periodically as the crankshaft rotated with almost the same period, shape, and peak-to-peak value of 3.04 arcseconds, while their mean absolute values were slightly different due to the different types of TLM resulting in different deflection angles. The GR for the four cases showed good agreement in terms of the period, shape, and peak-to-peak value of 0.29, meaning that TLM has a rather negligible influence on the GR of cycloid drives. with almost the same period, shape, and peak-to-peak value of 3.04 arcseconds, while their mean absolute values were slightly different due to the different types of TLM resulting in different deflection angles. The GR for the four cases showed good agreement in terms of the period, shape, and peak-to-peak value of 0.29, meaning that TLM has a rather negligible influence on the GR of cycloid drives.

Effect of Applied Torque Levels
The applied torque should be limited within a certain range to prevent the contact stress exceeding the elastic limit of the material, which can be determined by using this verified model. Taking the case of an NNM and LCM cycloid gear pair in the previous section as an example, the effects of various applied torque levels from 50 to 550 Nm on the contact stress and contact pattern were investigated, as shown in Figure 13.

Effect of Applied Torque Levels
The applied torque should be limited within a certain range to prevent the contact stress exceeding the elastic limit of the material, which can be determined by using this verified model. Taking the case of an NNM and LCM cycloid gear pair in the previous section as an example, the effects of various applied torque levels from 50 to 550 Nm on the contact stress and contact pattern were investigated, as shown in Figure 13. It was observed that under lightly loaded conditions, the example cycloid gear pair had a lower number of tooth pairs in contact-only one or two. As the applied torque increased, this number increased due to the deflection of the contacting tooth pairs, to four under torque of 550 Nm in this case. Figure 13 also shows the contact pattern predicted by the proposed model. It can be seen that the contact was localized at the middle of the tooth width of the gear with no edge contact when the load distributed on the tooth was low, referring to the contact pattern on roller number 9 under torque It was observed that under lightly loaded conditions, the example cycloid gear pair had a lower number of tooth pairs in contact-only one or two. As the applied torque increased, this number increased due to the deflection of the contacting tooth pairs, to four under torque of 550 Nm in this case. Figure 13 also shows the contact pattern predicted by the proposed model. It can be seen that the contact was localized at the middle of the tooth width of the gear with no edge contact when the load distributed on the tooth was low, referring to the contact pattern on roller number 9 under torque of 150 Nm. Increasing torque caused the contact pattern to extend towards the edge of the tooth, accompanied by increasing contact stress. Figure 14 shows the effects of the various applied torque levels on the LTE and GR. It can be seen that as the torque level increased, the peak-to-peak and absolute mean values of LTE, which is represented by the error bars, as well as the peak-to peak values of GR, changed proportionally with torque, which is a result of the possibility of the contact deformation becoming significant with the increase of the load distributed on the tooth. For instance, at torque of 150 Nm, the peak-to-peak and absolute mean values of LTE were 1.70 arcseconds and 14.49 arcseconds, respectively, and the peak-to-peak value of GR was about 0.08. As the torque increased to 450 Nm, the values increased proportionally to 3.77 arcseconds, 24.19 arcseconds, and 0.26.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 20 Figure 14 shows the effects of the various applied torque levels on the LTE and GR. It can be seen that as the torque level increased, the peak-to-peak and absolute mean values of LTE, which is represented by the error bars, as well as the peak-to peak values of GR, changed proportionally with torque, which is a result of the possibility of the contact deformation becoming significant with the increase of the load distributed on the tooth. For instance, at torque of 150 Nm, the peak-to-peak and absolute mean values of LTE were 1.70 arcseconds and 14.49 arcseconds, respectively, and the peakto-peak value of GR was about 0.08. As the torque increased to 450 Nm, the values increased proportionally to 3.77 arcseconds, 24.19 arcseconds, and 0.26.

Conclusions
In this study, a semi-analytical load distribution model was proposed based on a threedimensional and linear elastic solution. Different types of tooth profile and longitudinal modifications of cycloid gears were considered in this model to simulate and predict the multi-tooth contact conditions of cycloid drives. Numerical examples of load distribution analysis of cycloid drives using the proposed model were presented and compared to the current Hertz contact stiffnessbased model to demonstrate the effectiveness of the predictions of this semi-analytical model. From the analysis results, the following specific observations and conclusions can be drawn: 1. Load distribution, the number of tooth pairs in contact, and the loaded transmission error are influenced significantly by different types of tooth profile modification having the same radial clearance. NNM tends to yield smaller LTE than PPM, while the contact stress is quite the contrary. Both of them need to be optimized with proper modification parameters for the best performance of the cycloid drive.
2. Tooth longitudinal modifications play a positive role in relieving the stress concentration at the two edges of cycloid gears and improving stress distribution along the axial direction, which should be considered seriously in the analysis and design of cycloid drives. Among the three proposed types of TLM, LCM exhibits better stress distribution characteristics, and further research is needed on the determination of modification parameters.
3. Applied torque levels were observed to have significant effects on the contact stress, LTE, and GR. The magnitudes of their maximum, peak-to-peak values change proportionally with torque. The contact pattern extends towards the edges of the tooth as the torque increases.
This proposed model can be applied to determine the limitation of the applied torque to prevent the contact stress exceeding the elastic limit of the gear material. Furthermore, due to the model's

Conclusions
In this study, a semi-analytical load distribution model was proposed based on a three-dimensional and linear elastic solution. Different types of tooth profile and longitudinal modifications of cycloid gears were considered in this model to simulate and predict the multi-tooth contact conditions of cycloid drives. Numerical examples of load distribution analysis of cycloid drives using the proposed model were presented and compared to the current Hertz contact stiffness-based model to demonstrate the effectiveness of the predictions of this semi-analytical model. From the analysis results, the following specific observations and conclusions can be drawn: 1. Load distribution, the number of tooth pairs in contact, and the loaded transmission error are influenced significantly by different types of tooth profile modification having the same radial clearance. NNM tends to yield smaller LTE than PPM, while the contact stress is quite the contrary. Both of them need to be optimized with proper modification parameters for the best performance of the cycloid drive.
2. Tooth longitudinal modifications play a positive role in relieving the stress concentration at the two edges of cycloid gears and improving stress distribution along the axial direction, which should be considered seriously in the analysis and design of cycloid drives. Among the three proposed types of TLM, LCM exhibits better stress distribution characteristics, and further research is needed on the determination of modification parameters.
3. Applied torque levels were observed to have significant effects on the contact stress, LTE, and GR. The magnitudes of their maximum, peak-to-peak values change proportionally with torque. The contact pattern extends towards the edges of the tooth as the torque increases.
This proposed model can be applied to determine the limitation of the applied torque to prevent the contact stress exceeding the elastic limit of the gear material. Furthermore, due to the model's better computational efficiency, it is very suitable for parameter sensitivity studies to optimize the design of cycloid drives. In the future, an experimental study will be carried out to validate the proposed model.