Tooth Position and Deformation of Flexspline Assembled with Cam in Harmonic Drive Based on Force Analysis

Deformation of the flexspline is the basis of analyzing tooth trajectory and designing tooth profile. Considering the tooth influence on the position of equivalent neutral layer, a piecewise method for calculating the deformation of flexspline assembled with a cam wave generator is presented in this paper. Firstly, a mechanic model of a ring of uniform thickness in contact with a rigid cam is established. The displacements of the ring inside and outside an unknown wrapping angle are determined by the geometric constraints of the cam profile and the equilibrium relationship, respectively. Meanwhile, the wrapping angle is solved according to the boundary conditions. The assembly forces are derived to investigate the circumferential elongation and strain. Then, considering the tooth effects on the neutral layer of flexspline, the tooth is positioned on the equivalent neutral layer, which is the non-elongation layer within one gear pitch but offset from the geometric mid-layer. The equivalent neutral layer is positioned by the empirical formula of the offset ratio, which is summarized by the orthogonal simulation on finite element models of racks. Finally, finite element models of a ring-shaped and a cup-shaped flexspline assembled with elliptical cam are established to verify the effectiveness and accuracy of the piecewise method. The results show that, compared with the geometric method, the tooth positioning deviation calculated by the piecewise method can be reduced by about 70% with a more accurate deformation description from the geometric condition and mechanic condition inside and outside the wrapping angle.


Introduction
Harmonic drive is a transmission mechanism invented by Musser [1]. It consists of three main components: wave generator (WG), flexspline (FS), and circular spline. The rotating WG drives the FS producing oval-shaped deformation periodically to realize the meshing between the FS tooth and the circular spline tooth and to transfer torque and motion. Harmonic drives possess a number of advantages: compact size, light in weight, large gear ratio in one stage, high load-carrying capacity, small moment which lays a theoretical foundation for the analysis of deformation and force of the FS. Based on this, Ivanov [7] and Shen [10] conducted a more detailed study on the FS deformation under the action of various WGs such as four-rollers WG, double-discs WG, and cam WG. Dong et al. [12] calculated the assembly deformation of a cupshaped FS and fitted the loaded deformation results with cubic spline curves. Ma et al. [13] developed an experimental apparatus that integrates a special-fabricated micro-displacement platform and a pair of laser sensors to measure the radial displacement of the FS and investigated the effect of the driving speed on the deformation characteristics of the FS.
At present, the theoretical analysis of FS deformation is still based on the assumption of small deformation and non-elongation of the neutral layer. Chen et al. [14] found that due to the assumption of small deformation, the arc length of the deformed neutral curve expressed by the displacements is a little longer than the original arc. Yang et al. [15] conducted a more rigorous discussion on the deformed neutral curve, which described by displacements, and pointed out that the correct curve length can be obtained by the integration on the curve's geometric shape of the deformed neutral layer, which can be expressed by the equidistant line of cam WG profile. The cam WG not only makes the FS maintain a more stable deformation and ideal meshing state under transmission torque but also produces the FS deformation to any shape through a proper design of the cam profile. Therefore, the cam WG is the most widely used form in the harmonic drive.
With the assumption of FS tightly fitted with the cam WG, the FS deformation under the action of cam is usually determined directly by the geometry of the cam, which dramatically simplifies the design process. Maiti [16] provided a cam formed by splicing arcs at the major axis and elliptical arcs at the minor axis to make the standard involute tooth work normally. Gravagno [17] used a cam follower system to study the influence of the cam WG profile on the meshing performance of the FS. In addition, Routh et al. [18] studied the effect of tooth pitch variation on the contact between the teeth pairs on the elliptic pitch curve of FS that formed by the action of elliptical cam. In fact, there is a certain gap between the cam and the FS. Routh et al. [19] and Li et al. [20] considered the clearance caused by the taper deformation of cup-shaped FS and by the roughness between the FS-cam contact surface, respectively, when studying the hydrodynamic lubrication.
In addition, in order to reduce the FS stresses, a clearance fit is usually designed between the FS and the WG. Therefore, the assumption that the WG and the FS contact over the entire circumferential area is inappropriate.
At the same time, slight circumferential elongation of the FS neutral layer under the action of the assembly force may lead to separation on a small scale. Therefore, it must be more realistic to consider that they are only in contact in some areas around the major axis of WG. Chen et al. [14,21] conducted theoretical studies on the neutral layer elongation of FSs under the action of doubledisks WG and four-roller WG and verified the elongation by finite element simulation. The research confirms that the circumferential elongation of the FS neutral layer exists, but it is far less than the extra length obtained by the integration the displacement under the assumption of small deformation. Considering a ring model in contact with WG in partial along the circumferential direction, mechanic analysis in contact region will help us to find a more realistic details of the deformation and the internal forces inside and outside of the more reasonable wrapping angle.
At present, the mid-layer radius of the ring model used to calculate the FS deformation is considering the same with the radius of the geometric mid-layer (GML) of FS tooth rim, which ignored the tooth. However, the tooth not only affects the stress [22,23] and meshing [24,25] of the FS but also causes the change of the actual neutral layer (ANL) on position and shape. Kondo [9] performed a bending test on an enlarged simple trapezoidal rack to obtain the ANL distribution. His research shows that the ANL of rack is not a straight line, and its position deviates from the GML. The shape and position change of the ANL makes the previous ring model have a certain model deviation. Considering the influence of the tooth to improve the ring model can reduce the deviation of the FS deformation results and improve the positioning accuracy of the tooth. This paper is aiming to calculate the displacement of FS accurately. The forces and displacements on the equivalent neutral layer (ENL) of the FS are studied under the action of a cam WG. In Section 2, the existing methods for tooth positioning based on the deformation of FS are described, and the causes of tooth position deviation are analyzed. In Section 3, the actual contact state and forces between the cam and the ring are analyzed, and a piecewise method of calculating the ring's displacement is proposed, which is verified by finite element analysis of a ring-cam assembly model with shell elements. In Section 4, a bending rack model is established with 2D solid elements to analyze the influence of the tooth geometric parameters on the shape and position of ANL. The concept of ENL is introduced, and an empirical formula of ENL offset ratio is fitted through orthogonal simulations. The mid-layer of the ring model is migrated to the position of FS's ENL. In Section 5, validation models of a ringshaped FS and a cup-shaped FS, which all assembled with an elliptical cam, are introduced. In Section 6, The simulations of the validation models are carried out, and the ring-shaped FS is used to verify the theoretical method of this paper. The deformation and force characteristics of cup-shaped FS are also discussed. Finally, some conclusions of this paper are summarized in Section 7. Ivanov [7] constructed the tooth profile coordinate system to position the tooth on the mid-layer of the FS and calculated the radial and circumferential displacements of the tooth based on the deformation and puts forward a graphical method for calculating the tooth trajectory and the backlash between the teeth pairs. Shen [10] proposed an envelope method for conjugate tooth profile design based on the deformation curve ρ(ϕ ) constructed by the displacement on the mid-layer of the FS to position the tooth. Figure 1 shows the correspondence relationship in the envelope method. Figure 1 illustrates that the mid-layer circle (dash line) of the FS with a radius of r m is deformed into a curve shape (heavy line). Coordinate system {O-XY} is fixed on the WG, and the Y and X axes coincide with the major and minor axes, respectively. Coordinate system {O 1 -X 1 Y 1 }, fixed on the FS tooth, is located on the deformed curve of the initial circle, and the Y 1 -axis coincides with the tooth symmetry line. Coordinate system {O-X 2 Y 2 }, located at the rotating center, is fixed on the circular spline. As the FS deformed, its tooth positioning point O 1 (at polar angle ϕ to the Y-axis) moves to O 1 ' (at polar angle ϕ 1 to the Y-axis) with the following displacements and orientation: radial displacement w(ϕ ), circumferential displacement v(ϕ ), and the normal rotation angle θ(ϕ ). By using the theories of material mechanics, Ivanov [7] and Shen [10] performed rigorous force analysis on various deformed models, such as four-rollers model and double-discs model, under the assumptions of small deformation and non-elongation on mid-layer, and derived the expressions and relationships of w, v, and θ.

Tooth Positioning Methods and Deviations
Generally, the deformed shape of the mid-layer is expressed in polar coordinates: The deformed curve of the mid-layer expressed by Eq.
(1) in the envelope method [10] is used to position the tooth. In general, Eq. (1) is also used to construct equidistant cam profiles to maintain the designed curve deformation. However, the expression of the curve by Eq. (1) is not accurate enough. It can be seen in Figure 1 that the polar angle corresponding to the polar diameter ρ should be ϕ 1 , not ϕ . Therefore, ρ(ϕ ) in Eq. (1) only describes the polar diameter accurately, but do not describe the polar angle clearly. Chen et al. [14] obtained the deformed curve length of the mid-layer by integrating ρ(ϕ ) with respect to ϕ and found that the deformed curve length is larger than the arc length calculated on the initial circle. Moreover, the geometric elongation is much more significant than the elongation caused by the circumferential assembly force. Considering this deviation, the current design manual [11] has made the following approximate correction to the polar angle: Chen et al. [26] derived the relationship between ϕ 1 and ϕ according to the assumption of non-elongation on the mid-layer, and substitute ϕ 1 for ϕ to position the tooth more accurately in the conjugate equation [10]. However, the object of integration in Eq. (3) is still ρ(ϕ ), so the deviation caused by geometric elongation still exists. Yang et al. [15] proposed that the deformed curve length of the mid-layer under the action of an elliptical cam should be integrated with polar diameter ρ(ϕ 1 ), which is directly expressed by its real polar angle, ϕ 1 , to determine the relationship between ϕ 1 and ϕ: The relationship between ϕ 1 and ϕ expressed by Eq. (4) is more accurate than that in Eq. (3). Based on the assumption of non-elongation of the mid-layer, the FS is assumed to be closely attached to the cam after deformation. Then the tooth can be positioned directly by ρ(ϕ 1 ). When calculating the tooth trajectory and backlash between tooth pairs, the tooth displacements are calculated as follows: The calculation in Eq. (5) can be performed according to the Geometric Equation [7] that satisfies the non-elongation assumption. It does not involve complicated force analysis and calculation, which can significantly simplify the design process. However, if the FS is not closely contacting with the cam, the deviations of displacements calculation will occur.
Besides, it is noted that the independent variable ϕ 1 of the displacements w 1 , v 1 , and θ 1 in Eq. (5) is the polar angle of ρ(ϕ 1 ). However, in the existing geometric relationship of the displacements, the independent variable of w, v, and θ is the polar angle ϕ on the mid-layer circle without deformation [7,10]. The descriptions of the displacements in terms of the deformed geometry will also cause displacement calculation deviations.
As mentioned above, factors such as the differences in the independent variables of deformation and displacement, the elongation of the FS under the assembly force, and the change in the contact state between the cam and the FS will cause deviations in the calculation of the FS deformation, and then affect the tooth positioning accuracy. Proposing a displacement calculation method based on force analysis of the FS and cam assembly might not only reveal the actual contact state of them but also reduce the deviations of displacement calculation and improve the accuracy of tooth positioning.

Model and Methodology Based on Force Analysis
At first, the model and methodology of deformation and force analysis of the FS is presented using the uniform thickness ring proposed by Ivanov [7], the tooth effect will be discussed in Section 4 later. As the calculation of the ring model is carried on the mid-layer, the cam profile is defined as contacting the ring on the mid-layer.

Ring Deformation Calculation Piecewise with Wrapping Angle
The analytical model of the ring deformed by the cam is shown in Figure 2 (due to the symmetry of the model, only the first quadrant is illustrated here). In Figure 2, the radius of the mid-layer of the ring (arc AC) is r m , the expression of the cam profile (curve A'C') is ρ(ϕ 1 ), the y-axis coincides with the cam's major axis, and the x-axis coincides with its minor axis. Under the assumption that the ring is in close contact with the cam, the deformed arc AC will coincide with curve A'C', so its deformed shape can also be expressed by ρ(ϕ 1 ). In this case, the tooth positioning can be achieved by ρ(ϕ 1 ) and Eq. (5). For the convenience of description, this method for calculating the displacement of the ring is named the Geometric Method. According to the content of Section 2, the ring and the cam may only be in contact partially, so the concept of wrapping angle is introduced to analyze the ring model.
Suppose that the segment AB (in Figure 2) on the ring continuously contacts the cam to form a curve A'B' after deformation, and the segment BC's deformed shape B'C' does not contact the cam. Then only the deformed shape A'B' can be determined with the expression of ρ(ϕ 1 ). Based on the assumption of non-elongation in the circumferential direction, the polar angle γ (corresponding to undeformed arc AB) and the polar angle γ 1 (corresponding to deformed curve A'B') have the relationship in Eq. (4). Define γ and γ 1 as the wrapping angles of ring AB before and after deformation, respectively. Therefore, the ring inside the wrapping angle (0 ≤ ϕ 1 ≤ γ 1 ) is in close contact to the cam, so the curve A'B' is constrained on the geometric profile of the cam, and its displacement can be calculated by Eq. (5); the arc BC, outside the wrapping angle (γ ≤ ϕ ≤ π/2), needs to solve the bending differential equation [7] to obtain the deformation calculation formula: where w 2 is the radial displacement, v 2 is the circumferential displacement, θ 2 is the normal rotation angle; E is the elastic modulus, and I is the moment of inertia of cross area. The unknown coefficients B 1 , B 2 , X 1 , and X 2, can be determined by boundary conditions (in Section 3.2.2). For the convenience of description, the method of calculating the ring displacement segmented with wrapping angle by the Eq. (5) and Eq. (6) is named the Piecewise Method.

Wrapping Angle and Ring Force Calculation
To determine the wrapping angle of the ring to the cam and the force condition of the ring, a force analysis of the ring deformed under the action of the cam is required. Inside the wrapping angle, the ring bears the radial external load caused by the contact of the cam, while there is no external load outside the wrapping angle, so the force analysis of the ring is segmented by the wrapping angle.

Force Analysis of Ring
According to the deformation analysis in Section 3.1, the shape change of the ring inside the wrapping angle (0 ≤ ϕ 1 ≤ γ 1 ) from AB to A'B' is known, then its curvature change can be obtained for calculating the bending moment M 1 (Figure 3a) on any cross-section by using the elastic equation [7]. In addition, there are also a circumferential force F N1 , a shear force F S1 , and a radial external load q r on any section within the A'B' segment. M γ , F Nγ , and F Sγ are the corresponding forces in the cross-section, which is on the wrapping boundary.
According to the equilibrium equation [7], the assembly forces on the ring inside the wrapping angle can be calculated as follows: where F N0 is the integral constant to be determined. Although the deformed shape of the ring outside the wrapping angle is unknown, it should meet the condition of force balance. Before deformed, this segment is an initial arc within γ ≤ ϕ ≤ π/2. When deformed, it has a circumferential force F N2 , a shear force F S2 , and a bending moment M 2 in any cross-section ( Figure 3b). The crosssection on the minor axis has bending moment X 1 and circumferential force X 2 . The forces of the ring outside the wrapping angle can be calculated as following [14]:

Determine Wrapping Angle with Boundary Conditions
Since the wrapping angle γ (or γ 1 ) is unknown, the abovementioned piecewise calculation conditions of deformation and force are currently incomplete. Based on the deformations and forces of the whole ring, the wrapping angle can be determined according to reasonable boundary conditions. Besides the wrapping angle, there are the other four unknown quantities: B 1 , B 2 , X 1 , and X 2 (Eq. (6)). Refer to the boundary conditions of analyzing of the wrapping angle of the ring to the double-discs, there are Substituting Eqs. (5)-(8) into the boundary conditions 1)-4), the following formulae can be obtained after simplification: where As γ and γ 1 has the relationship in Eq. (4), B 1 , B 2 , X 1 , and X 2 can be converted into variables about γ 1 . Then the wrapping angle can be solved with the following formula according to boundary condition 5): When γ 1 obtained from Eq. (10), the corresponding γ can be determined by Eq. (4). It can be noted from Eq. (10) that the wrapping angle only depends on the shape of the cam profile. As the wrapping angle is obtained, then B 1 , B 2 , X 1 , and X 2 can be determined by Eq. (9). Furthermore, the ring displacement outside the wrapping angle can be determined by Eq. (6), and the assembly forces expressed by Eq. (7) and Eq. (8) can also be determined.

Circumferential Strain and Elongation
In the theoretical calculation of the ring displacement, the circumferential non-elongation assumption is adopted, which is considered that the deformation is mainly formed by the shape change under bending moment rather than the elongation caused by the circumferential force. In fact, under the circumferential forces F N1 and F N2 shown in Figure 3, the ring will be extended with a specific elongation. This elongation is often ignored because it is too small, but it is enough to change the actual contact state of the FS to the WG, thus causing a specific tooth positioning deviation.
Since the circumferential force is continuous, the integral constant F N0 can be determined by the circumferential force in the cross-section at wrapping angle: According to the circumferential forces in Eqs. (7) and (8), the circumferential strain of the ring where χ 1 = −(ẅ 1 + w 1 ) r 2 m , (denotes the curvature change) The curvature change is the derivative of the normal angle to the arc length, i.e., χ 1 =θ 1 r m . Therefore, by integrating the circumferential strain ε N along the arc length, the circumferential elongation of the ring can be obtained

Forces and Displacements of Flexspline Assembled with Elliptical Cam
Segment analysis of ring deformation and force based on wrapping angle has previously always been the concept related to the double-disk WG. Because the radius of the eccentric disk is smaller than the radius of the ring, so the ring cannot be wrapped entirely on the disk. The formulae obtained from the present analysis is also applicable to solving the FS cases with a double-disk (Appendix 1).
To the best authors' knowledge, the deformation of FS assembled with cam has not been investigated based on force analysis. Take elliptical cam for calculation. The model parameters are tabulated in Table 1. The elastic modulus of the ring E = 196 GPa. The polar diameters on the major axis and minor axis of the elliptical cam in Table 1 are calculated by Cмиpнoв formula [10] and the expression of the cam profile is To save computer time, the contact model between elliptical cam and ring is established by shell elements in ANSYS environment, and the cam is simplified as a rigid ring with its elliptical shape. The contact elements are established between the inner surface of the flexible ring and the outer surface of the rigid elliptical ring (represents the cam). As the contact pressure of the model is illustrated in Figure 4, it can be seen that the contact between the elliptical cam and the ring exists only in a part of the area near the major axis of the cam, which indicates that the wrapping angle of the ring to the elliptical cam does exist.
According to Eq. (10), the wrapping angle γ 1 = 35.29° (γ = 35.56°) of the ring to the elliptical cam is obtained. Substituting it into Eq. (9), then B 1 = 2.553, B 2 = 3.486, X 1 = − 9.127 N·mm, X 2 = 0.58 N is determined. In the (14)    ρ a = r m + w 0 , finite element model (FEM), the bending moment of the ring at the minor axis is − 9.196 N·mm, and the circumferential force is 0.594 N. Compared with the finite element results, the deviations of the numerical solutions of X 1 and X 2 are 0.75% and 2.4%, respectively, which indicates that the numerical solutions are very close to the finite element results.

Displacements of Ring under Elliptical Cam
Calculate the displacements of the ring in Table 1 by the geometric method and the piecewise method and extract the displacement results from the FS for reference. Figure 5 illustrates the comparisons between two numerical and simulation results about the radial displacement, the circumferential displacement, and the normal rotation angle. In Figure 5, w t1 , v t1 , and θ t1 are numerical results of the radial displacement, the circumferential displacement, and the normal rotation angle calculated by the geometric method, respectively. w t2 , v t2 , and θ t2 are the numerical results of the corresponding displacements calculated by the piecewise method. UX, UY, and ROTZ are the simulation results of the FS. Figure 5 illustrates that the displacements calculated by the two theoretical methods are very close to the results of the FS. Also, the two numerical results of each displacement are entirely the same inside the wrapping angle, but slightly different outside the wrapping angle.
Based on the finite element results, two numerical results are used to minus the corresponding result of the FS to obtain the deviations of the displacements. As shown in Figure 6, Δw t1 , Δv t1 , and Δθ t1 illustrate the radial displacement deviation, the circumferential displacement deviation, and the normal rotation angle deviation of the geometric method. Δw t2 , Δv t2 , and Δθ t2 illustrate the corresponding deviations of the piecewise method.   Figure 6a illustrates that Δw t2 almost monotonously increases and Δw t1 (absolute value) increases rapidly outside the wrapping angle. Figure 6b shows that Δv t2 first increases and then decreases, reaching the maximum at the wrapping angle, while Δv t1 increases rapidly outside the wrapping angle. This is because v t2 is integrated from the minor axis to the major axis (Eq. (6)), so Δv t2 accumulates to zero at the minor axis. However, v t1 is still an integral starting from the major axis to the minor axis, even outside the wrapping angle (Eq. (5)), so its deviation Δv t1 is accumulated from the major axis to the minor axis. Figure 6c shows that Δθ t2 fluctuates near zero and finally reduced to zero, while Δθ t1 fluctuates greatly and does not decrease to zero at the minor axis. Figure 6 shows that, in general, the deviations of the two numerical results of each displacement are the same inside the wrapping angle, which are mainly caused by the independent variable deviation and circumferential elongation. Moreover, the calculation deviation of the piecewise method outside the wrapping angle is smaller than that of the geometric method. The calculation deviation of the piecewise method is mainly caused by the circumferential elongation, while the geometric deviation   is formed under the combined action of circumferential elongation, independent variable deviation and non-contact state outside the wrapping angle. It can be seen from Figure 6 that, overall, the piecewise method eliminates the displacement calculation deviation caused by these factors to a certain extent, so it has a higher calculation accuracy.

Assembly Force of Ring under Elliptical Cam
Calculate the numerical results of the assembly force of the ring, according to Eqs. (7), (8), and (11), and compared with the bending moment M11, tensile force N11, and shear force Q13 of the FS on unit width. There is no direct corresponding simulation result for the radial external load q r , so the load on the unit width of the ring is compared with the contact pressure of the FS. Figure 7 shows that the numerical results of the assembly forces of the ring are in good agreement with the finite element results. Among them, the bending moment M changes continuously and monotonously from the major axis ( ϕ = 0°) to the minor axis ( ϕ = 90°), and changes the direction at ϕ = 45° (Figure 7a); the shear force F S increases monotonously at first, then decreases monotonously after a sudden increase within a small range of about 2°-5° (Figure 7b), the area where the shear force suddenly changes is the wrapping angle position; the circumferential force F N increases monotonously from the major axis to the minor axis, and the growth rate increases at first and then slows down (Figure 7c). In Figure 7d, the external load q r of the ring decreases monotonously inside the wrapping angle, and most of the contact pressure of the FS of the ring inside the wrapping angle is basically consistent with that of q r . Different from the phenomenon that q r drops to zero directly outside the wrapping angle, the contact pressure increases sharply near the wrapping angle and then rapidly decreases to zero, which is due to the additional pressure caused by the shear jump at the wrapping angle of the FS.  The results in Figure 7 show that the piecewise method is basically accurate on the force analysis and calculation of the ring. Figure 6 shows that the elongation of the ring under the assembly force has a specific effect on the displacement deviation, especially on the circumferential displacement deviation. Calculate the circumferential strain and elongation of the ring, according to Eq. (12) and Eq. (13), and extract the corresponding results in the FS for comparison. The results are illustrated in Figure 8. The circumferential strain in Figure 8a has the same trend as the circumferential force in Figure 7c and its numerical result is basically consistent with the finite element result. The elongation of the ring obtained by integrating the circumferential strain along the arc length monotonically increases, and its numerical result is basically the same as the simulation result ( Figure 8b). In Figure 8b, when ϕ = 90°, the total elongation of the ring in the first quadrant under the assembly force is about 0.16 μm.

Circumferential Strain and Elongation of the Ring under the Elliptical Cam
Apply a 50 N tension force on the ENL (if the force is applied in other positions, the rack will bend) of the rack model (parameters in Table 2, Section 3.3), the elongation of the rack is 7.7 μm. Under the same conditions, the elongation of the smooth beam, which thickness is the same as the rack's wall thickness, is 8 μm. The elongation of rack is only about 3.75% smaller than that of the smooth beam. Therefore, the forces of the ring can be regarded as that of the FS.

Offset Correction of the Flexspline Equivalent Neutral Layer
In Section 3, the force and deformation analysis of the ring is based on the assumption of circumferential nonelongation of the mid-layer, in which the radius rm is an important parameter. As the FS is simplified to a uniform thickness ring, the ENL is exactly its GML. However, when considering the tooth structure, there is an offset of the ENL, or even it is no longer in a cylindrical shape. The offset of the GML to the ENL will cause the offset deviation in both the geometric method and the piecewise method.

Equivalent Neutral Layer of Flexspline
Kondo [9] conducted bending tests on simple trapezoidal racks with gear modulus of m = 32 mm and with different rack thickness and found that the ANL of the rack fluctuates with the tooth structure (Appendix 2). As for the FS, its ANL will also fluctuate along the cylinder, so it is impossible to effectively calculate the deformation or stress based on such an ANL shape. Therefore, find a cylindrical surface on the FS: when the FS is bent, although the circumferential strain of the cylindrical surface is not always constant to zero, the strain integral within a gear pitch equals to zero, i.e., the arc length of the layer within a tooth pitch is maintained constant. This cylindrical surface is defined as the ENL of the ring gear, and its radius is r en . Take a short section of the FS and reduce it to a straight rack model for bending analysis. Figure 9 shows the three layers when the rack model is bent: the ANL, the GML, and the ENL. In Figure 9, h n is the full tooth depth, h f (a) Results of circumferential strain (b) Results of circumferential elongation  is the dedendum height, s f is the tooth root thickness (defined as the thickness of the straight tooth profile extending to the tooth root), α is the tooth profile angle, r c is addendum arc radius, and r i is dedendum arc radius. d mn is the relative offset distance between the ENL and the GML. In order to express the degree of relative offset between the GML and the ENL, the ENL offset ratio is defined as the ratio of the offset distance to half the wall thickness As illustrated in Figure 10, Establish the FS of the rack by 2D solid elements in ANSYS environment, solve the model with a specific bending moment, for instance, 50 N·m. Then perform the following operations in the post-processing: 1) Establish multiple paths parallel to the Y-axis (Figure 9), and map the X-strain results to the path, connect the points with zero strain result on each path to obtain the ANL distribution of the rack; 2) Establish a path parallel to the X-axis according to the initial Δ mn , and map the X-strain results to the path; 3) integrate the strain in the X-direction to obtain the Elongation Δ p within a gear path; 4) iterative the Y-coordinate of the path to make the elongation within a gear pitch approach to zero.

Circumferential Strain and Elongation on Mid-layer due to Offset from Equivalent Neutral Layer
When the FS is deformed by inserting a cam, the assembly force, especially the circumferential force, will make both GML and the ENL producing a specified elongation. Besides, an amount of expansion or contraction produced by circumferential strain caused by bending will be superimposed on the GML due to the offset from the ENL. The circumferential strain on the GML caused by bending can be calculated according to the ENL offset ratio where the radial displacement w is applicable to be provided by any of the geometric method or the piecewise method. The addition elongation of the GML caused by ENL offset can be obtained by integrating ε mn along the arc length

Mid-layer Offset Characteristics of Rack
By using the finite element method shown in Figure 10, the simple trapezoidal rack tested by Kondo [9] is modeled and calculated and then compared with the (17) ε mn = � mn δ 0 2r 2 Figure 9 Schematic diagram of rack structure and equivalent neutral layer experimental results. The simulation results show good agreement with the experimental results (Appendix 2), which verified the accuracy of the present simulation method. In order to analyze the mid-plane characteristics of rack considering the ENL offset, a parametric rack model with detailed tooth profile features (Figure 9, Table 2) is established and the simulated according to the flowchart illustrated in Figure 10. Elastic modulus of rack E = 196 GPa. The position of the ANL of the rack obtained by the process shown in Figure 10 is illustrated in Figure 11. because of the symmetry of the structure, the horizontal axis is valued within half the tooth pitch (πm/2), from the tooth symmetry line to the tooth space symmetry line.
It can be seen from Figure 11 that: 1) the ANL of the rack has a wavy shape that deviates outward from the GML (> δ 0 /2); 2) the offset distance at different positions on the ANL varies continuously, and the maximum offset is achieved at the stress concentration position (near the connection between the dedendum arc and the root line); 3) the offset of the ANL at the position of the tooth slot symmetry line is more significant than that at the tooth symmetry line. These results show that the offset position of the ANL is not only related to the tooth structure but also affected by the stress concentration at the tooth space. Thus, the tooth dedendum arc, which can significantly influence the stress concentration, should also be considered when analyzing the offset of the neutral layer.

Induction of Equivalent Neutral Layer Offset Ratio Formula
According to Section 4.3, the geometric parameters of the tooth may affect the position of ANL. Analyze the sensitive parameters of the offset ratio of ENL by using the dimensional analysis, and then summarize the empirical formula of the offset ratio can provide a theoretical basis of the correction of the ring model for calculating the FS deformation and force.

Sensitive Parameter Analysis
According to the rack shown in Figure 9, the structural parameters that may affect the ANL offset include tooth depth, tooth thickness, dedendum arc radius, tooth profile angle, and rack thickness. Among them, in the univariate analysis, if the tooth thickness is constant, the change of tooth profile angle can be reflected by the change of tooth root thickness, so the change of tooth profile angle and tooth thickness are considered together as the change of tooth root thickness. Rack experiments [9] and simulations (Appendix 2) have shown that the rack thickness has a significant effect on the offset of the rack's ANL. Therefore, three structural parameters: tooth depth, tooth root thickness, and dedendum arc radius, are taken to be analyzed here. Taking the parameters in Table 2 as the basis and perform univariate simulations. The change of the rack's ANL with the three parameters is shown in Figure 12. Figure 12a shows that when the tooth depth coefficient h n * valued with 1.65, 2 and 3 (h n * = 3 exceeds the standard value of tooth height coefficient 2.35, here only valued for the purpose of analyzing the influence of tooth depth), the ANL distribution of the rack is basically the same, indicating that the change in tooth depth has almost no effect on the ANL when the tooth depth is enough. Figure 12b and Figure 12c show that when the tooth root thickness coefficient K s valued with 0.45, 0.55 and 0.65, or when the dedendum arc radius coefficient r i * valued with 0.02, 0.4, and 0.8, the ANL of the rack is significantly raised and the shape changes. Besides the rack thickness, it is illustrated in Figure 12 that the sensitive parameters of ANL offset also include tooth root thickness and dedendum arc radius.

Dimensional Analysis of Equivalent Neutral Layer Offset Ratio
According to the Π theorem in dimensional analysis [27], Δ mn as the dependent variable is a dimensionless quantity, and the three independent variables: δ 0 , s f , and r i , have the same dimension of L and with a power of 1.
Since the independent variables have the same dimension, one of them can be taken as the basic quantity. In this problem, it is appropriate to take δ 0 as the basic quantity. Then the functional relationship between the new dimensionless variables is where Π 1 (Δ mn ) represents the ENL offset ratio, Π 2 (s f /δ 0 ) reflects the relative thickness of the tooth (i.e., the width of the tooth root relative to the rack thickness), and Π 3 (r i /δ 0 ) reflects the size of the dedendum arc relative to the rack thickness. In order to obtain the function f of Eq. (19), 25 finite element rack models are designed and established for orthogonal simulation according to the parameter value recommendations [7,10]. The ENL offset ratio Δ mn , tabulated in Table 3, is obtained by the method illustrated in Figure 10. Due to the coupling constraint between the maximum dedendum arc radius and the maximum tooth root thickness, the items '-' in Table 3 indicate that the group of parameters exceed the geometric constraint. Table 3 shows that the ENL offset ratio of the rack increases monotonously with the relative tooth root thickness and increases monotonically with the relative dedendum arc radius. The results in Table 3 expressed as contour lines illustrated in Figure 13. The offset ratio has a linear relationship with the relative root thickness and the relative dedendum arc radius. When expressed in the Cartesian coordinate system, the offset ratio will appear as a spatial plane.
Fit the ENL offset ratio Δ mn in Table 3 with the general equation of the space plane as follows: After fitting, the parameters in Eq. (20) are: Π 0 = − 7.3191, A = 11.8458, B = 13.6256. The R-Square and Reduced Chi-Sqr of the fitting are 0.9891 and 3.455 × 10 −2 , respectively. The contour map of the fitted results of the ENL offset ratio is illustrated in Figure 14.
Then the ENL offset ratio Δ mn has the following relationship with tooth root thickness, dedendum arc radius, and rack thickness (or wall thickness of FS):  The empirical formula of the ENL offset ratio Δ mn makes it possible to calculate the position of the ENL of the FS numerically and to improve the ring model for calculating the displacement of FS. Alternatively, it can be used to evaluate the tooth circumferential positioning deviation when calculating the FS displacement on GML.

Solid Finite Element Models of Flexspline-cam Assembly for Verification
In Section 2, The assembly forces and deformation of the FS under the action of cam WG are analyzed, and then the offset position of the ENL of the FS considering the tooth effect is obtained in Section 3. To verify the accuracy of the piecewise method in calculating the displacement on the ENL of the FS, the solid FS models with detailed tooth information are established to assemble with ring-shaped cams for contact analysis.

Model of Ring-shaped Flexspline Assembling Elliptical Cam
The structural and tooth parameters of the ring-shaped FS are tabulated in Table 4. The FS is modeled by 3D-solid elements, and the cam, modeled by shell elements, is designed to have the same width of the FS. The profile of cam is calculated with the inner surface radius of FS, according to Eq. (14) and Eq. (15). Contact elements are arranged between the inner surface of the FS and the outer surface of the cam. The element size near the mid-layer of the FS is 0.035 mm. The assembly model is shown in Figure 15. the elastic modulus of the model material is E = 196 GPa, the Poisson's ratio and friction coefficient are set to zero.
(21) � mn = (−7.3191δ 0 + 11.8458s f + 13.6256r i ) 100δ 0 . Figure 16 illustrates the structure of a cup-shaped FS, which is formed by adding a cup to the ring-shaped FS in Section 5.1. In Figure 16, S f , S m , and S b are the front, middle, and rear cross-section positions on the tooth   The material settings are the same as Section 5.1 but considering the Poisson's ratio (υ = 0.3). As the cupshaped FS deforms, the outer ring of the flexible bearing will achieve a taper deformation accordingly. To avoid bearing deformation affecting the results, set the WG to an elliptical ring with a width of 0.4 mm. The parameters of the tooth rim of the cup-shaped FS are the same as those in Table 4, and the remaining parameters are tabulated in Table 5.

Results of Ring-shaped Flexspline Assembling Elliptical Cam
According to Eq. (21), the ENL offset ratio of the model is Δ mn = 8.86%. Referring to the parameters in Table 4, the radii of the GML and the ENL of the FS are: r m = 29.119 mm and r en = 29.152 mm, respectively. Set the maximum radial displacement w 0 = 0.375 mm and perform numerical and simulation calculations on the deformation and force of the ring-shaped FS.

Results of Circumferential Strain and Elongation
Extract the circumferential strain elongation of the FS on GML and ENL in the FS to compared with the numerical results. According to Eq. (12) and Eq. (13), there are only circumferential strain ε en = ε N and elongation Δl en = Δl N on the ENL of the ring-shaped FS. Referring to Eq. (17) and Eq. (18), besides ε N and Δl N , on the GML of the FS, there will superimpose a circumferential strain ε mn and elongation Δl mn caused by offset from the ENL. Therefore, on the GML, the circumferential strain ε m = ε N + ε mn , and the circumferential elongation Δl m = Δl N + Δl mn . Figure 17 illustrates the results of circumferential strain and elongation. Because the ANL of the FS is fluctuating, then the finite element results of the circumferential strain extracted along the circular path of both the ENL and the GML are also fluctuating ( Figure 17a). It is shown that the numerical solutions of the circumferential strain on the ENL and the GML are basically consistent with the corresponding finite element results. The maximum absolute value of ε m is about 10 times the maximum value of ε en .
The numerical results of circumferential elongation of GML and ENL shown in Figure 17b are in good agreement with the finite element results. Among them, the elongation result Δl en on the ENL shows a monotone increasing trend. Affected by the variation of circumferential strain ε m , the GML is compressed at first and then elongated along the circumferential direction. The maximum absolute value of Δl m is about 3 times of the maximum value of Δl en .
The above results show that the piecewise method for analyzing the ring model can accurately solve the forces of the ring-shaped FS. meanwhile, the empirical formula of the ENL offset ratio (Eq. (21)) is proofed to be accurate enough by the circumferential strain and elongation on the GML.

Results of Tooth Positioning
Use the intersection of the tooth symmetry line with the geometrical mid-layer or with the ENL of the FS as the tooth positioning point. On the FS model in Section 6.1.1, the radial displacement and circumferential displacement of the tooth positioning point are extracted on each tooth positioning point. However, the rotate angle of the tooth symmetry line cannot be extracted from the model directly. Then the angle can be calculated according to the radial displacement and circumferential displacement as follows: where R m , w m , and v m are the radial and circumferential displacements of the tooth positioning point, while w a and v a are the radial and circumferential displacements on the middle point of tooth top. R a is the outside radius of FS. Position the tooth on the GML and calculate its displacement by the geometric method. Then minus the displacement results by that of the FS on GML. The deviations are recorded as Δw g , Δv g , and Δθ g . Then position the tooth on the ENL and calculate its displacement by piecewise method. Then minus the displacement results by that of the FS on ENL. The deviations are recorded as Δw p , Δv p , and Δθ p . Figure 18 illustrates the comparison of these displacement deviations.
The results in Figure 18 are similar to those in Figure 6. The change of Δw g and Δw p in Figure 18a are consistent with those of Δw t1 and Δw t2 in Figure 6a, respectively, indicating that the offset of the neutral layer has little effect on the radial positioning of the tooth. In terms of maximum deviation, Δw p is about 75% smaller than Δw g .
The change of Δv p in Figure 18b is consistent with that of Δv t2 in Figure 6b, while circumferential deviation, which is caused by neutral layer offset, is superimposed on Δv g compared to Δv t1 , indicating that the neutral layer offset has a significant influence on the circumferential positioning of the tooth. In terms of maximum deviation, Δv p is about 77% smaller than Δv g .
The change of Δθ g and Δθ p in Figure 18c are consistent with those of Δθ t1 and Δθ t2 in Figure 6c, respectively. However, due to the change of the circumferential displacement deviation Δv g , Δθ g also increases to a certain extent, indicating that the neutral layer offset also has a specific effect on the positioning accuracy of tooth orientation. In terms of maximum deviation, Δθ p is about 70% smaller than Δθ g .
According to the above results, the method proposed in the present investigation to calculate the FS deformation on the ENL can significantly reduce the existing (22) calculation deviation and obtain a more accurate tooth positioning. Figure 19 illustrates the radial displacement of the cupshaped FS. It is noted that under the constraint of the cup bottom, the radial deformation of the cup-shaped FS has a taper characteristic, and the amount of radial deformation decreases linearly from the cup edge to the cup bottom. Establish paths at the intersection of the tooth rim's ENL (determined by Eq. (21)) with the three crosssections shown in Figure 16. Extract the circumferential strain and integrating the strain to obtain the circumferential elongation. The results are illustrated in Figure 20.

Results of Cup-shaped Flexspline Assembling Elliptical Cam
It is illustrated in Figure 20a that, from the major axis ( ϕ = 0°) to the minor axis ( ϕ = 90°), the average circumferential strain of ENL on the front section is positive and increases monotonously, the average strain on the middle section is also positive but decreases gradually, while the average strain on the rear section changes significantly from positive to negative and is close to zero around ϕ = 50°. On the major axis ( ϕ = 0°) of view, the circumferential strain of ENL is the largest in the rear section, the second in the middle section, and the smallest in the front section. These strains are all positive. On the minor axis ( ϕ = 90°) of view, the strain of ENL in the front section is the largest, the strain in the middle section is the second, and the strain in the rear section is the smallest and negative.
These phenomena indicate that ENL will gradually elongate on the front section and middle section, while ENL will elongate rapidly at the rear section at first, and then begin to shorten rapidly around ϕ = 50°. The elongation result in Figure 20b confirms this conclusion. The strain distribution and circumferential elongation illustrated in Figure 20 indicate that the equivalent neutral layer analyzed based on the ring-shaped FS in the present investigation has not accurately reflected the actual neutral layer characteristics of the cup-shaped FS.
Under the comprehensive influence of the cup-shaped structure, cylinder stiffness, and Poisson's ratio, the ENL of the cup-shaped FS tooth rim gradually changes from the front section to the rear section, forming a conical equivalent neutral layer. This tapered ENL makes the stress and deformation of the FS tooth rim more complicated. Therefore, besides the tooth effects, the influence of the FS cylinder on the stress and deformation of tooth rim deserves further study.

Conclusions
In this paper, the actual contact state and assembly forces between the cam and the flexspline are revealed based on mechanic analysis, and the empirical formula of the equivalent neutral layer offset ratio of the flexspline is fitted out. A piecewise method for calculating the flexspline deformation in the equivalent neutral layer is proposed. The accuracy of the theoretical calculation of tooth displacement is verified by establishing a solid finite element model of flexspline contacted elliptical cam. The main conclusions are as follows: (1) The flexspline has a wrapping angle far less than 90 to the cam. After deformation, the flexspline segment inside the wrapping angle is wrapped on the cam, while the outside segment deformed by force balance.