Study on Dynamic Coupling Behavior of End-Meshing Harmonic Reducer

1. Introduction

The utilization of industrial robots and aerospace exploration activities, among other sectors, necessitates the critical joint systems of their institutions to comply with the significant reliability and bearing capacity requirements of harmonic gear transmission devices [1]. To mitigate the inherent shortcomings of flexspline within traditional harmonic reducers, which are susceptible to fatigue damage, and to guarantee the precise operational state of the transmission system, the end-face engagement is employed to effectuate predetermined rigid thrust transmission. This solution effectively reconciles the conflict between deformation and load.

The dynamic characteristics of harmonic gears are inherently complex due to the periodic changes in meshing parameters such as tooth deformation, meshing stiffness, and tooth backlash during the meshing process. Such changes cause obvious nonlinear coupling effects and mechanical jitter, causing continuous vibration in the transmission, which can have a detrimental effect on transmission performance. Therefore, the study of the nonlinear dynamics of harmonic gear transmission systems can provide valuable engineering insights and a theoretical foundation for enhancing its performance.

Zhang Youlin et al. first postulated the concept of the end-face harmonic gear drive of oscillating teeth, elucidated its transmission principle, and subjected its kinematics law to comprehensive analysis. Subsequently, utilizing a virtual prototype model, they engineered and optimized the tooth profile to augment the transmission efficiency. The structural parameters pertinent to this drive system were calculated and optimized, however, the dynamic meshing behavior was not thoroughly explored.

Recognizing the influence of time-varying meshing stiffness, gear backlash, and other variables, various researchers have constructed a nonlinear dynamic model of planetary gear sets to depict the impact of internal and external excitations on vibration and shock responses. Furthermore, they offered strategies to mitigate vibration and shock. Zhu[9] et al. further extended the application of the harmonic balance method to the nonlinear dynamic modeling of compound planetary gear sets. Through their research, they were able to study the nonlinear dynamic characteristics of gear sets and identify the influence of multiple factors, including dimensionless backlash, meshing stiffness, and error excitation amplitude, on frequency response characteristics. Cui et al. [10-12] developed a coupled nonlinear dynamic model for compound planetary gears. This model was utilized to extract the natural frequency of the system and to evaluate the influence of moment of inertia, meshing stiffness, and other pertinent factors on the vibration response. Zheng[13] et al. proposed a translational-torsional coupled dynamics model of the RV reducer. This model was used to analyze the system's nonlinear time-varying behavior, yielding displacement responses of the various components under various conditions. Furthermore, the model was used to examine the sensitivity of the dynamic characteristics to both internal and external excitation factors. Tung[14] et al. developed a reduced-order time-variant numerical model for the compound reducer that enabled the prediction of the robot's dynamic stiffness with enhanced accuracy. Hu et al. [15] developed a dynamic analysis model of the three-stage planetary transmission within the wind turbine reducer. This model involved the utilization of the stiffness factor method to examine the mechanical properties of the interacting components. Subsequently, experimental validation was conducted to verify the reliability of the dynamic model. Yang [16] and colleagues established a dynamic model of planetary gear trains utilizing the lumped parameter approach and conducted in-depth research on the vibration responses of planetary gear trains when subjected to both deterministic and random loads. In consideration of the effects of backlash and time-varying stiffness, Saeed [17] et al. have refined the conventional spur gear dynamics model by incorporating Gaussian white noise into the loading terms. Liu Hui, Yang Wenguang[18, 19] et al. developed a non-linear dynamic model of planetary gear transmission systems, accounting for the effects of dynamic mesh parameters and stiffness. This model was subsequently subjected to experimental verification and analysis.

As per the findings from the literature review, current research on end-face harmonic gear drives predominantly involves kinematic analysis, with a limited focus on dynamic behavior. While numerous studies have been conducted on the system dynamics modeling of traditional harmonic and planetary gear reducers, it is evident that both internal and external excitations play a significant role in influencing the system's performance. Given the unique features of the end-face harmonic gear drive, the internal axial excitation has a significant effect on the dynamic behavior of the system, leading to significant buffeting phenomena that should not be overlooked. In contrast, traditional gear system dynamics modeling typically fails to analyze the axial meshing parameters based on the radial meshing form.

In light of this, this research aims to build a prototype of the specific end-face harmonic gear drive of oscillating teeth reducer. After a thorough examination of the dynamic meshing behavior law of the oscillating teeth, the influence of various dynamic meshing parameters, oscillating teeth deformation, time-varying stiffness, time-varying backlash, and other coupling factors on the nonlinear dynamic performance will be explored. Utilizing the lumped-mass method, a nonlinear dynamic model of the multi-tooth meshing system will be developed. Afterward, the theoretical model will be confirmed through experimental analysis to offer a theoretical basis for vibration suppression analysis of the end-face harmonic gear transmission mechanism.

2. System dynamic incentive analysis

In the process of gear transmission, the end-face harmonic gear drive of the oscillating teeth reducer is subjected to a combined effect of internal and external excitations. Specifically, internal excitation primarily involves the elastic deformation of gear teeth under load, deformation of the support system derived from the assembly relationship, and the combined influence of time-varying meshing stiffness as a result of the intermittent engagement of teeth. External excitation, on the other hand, is primarily affected by the combined effect of power sources at the input end and the fluctuation of load and torque at the output end. As a result, the gear transmission system experiences vibration and noise due to the combined action of internal and external excitations.

2.1 Multi-tooth Parameter Analysis of End-face Gear Drive

The principle of end-face harmonic gear drive of oscillating teeth [20] is shown in Figure 1-1. The specific configuration of the end-face harmonic reducer is shown in Figure 1-2.

The end-face cam rotates to push the oscillating teeth to move axially. The oscillating teeth always keep contact with the end face of the cam under the action of the spring, rise along the cam profile of the lift, and gradually engage with the gear teeth of the circular spline, which are forced to complete the rotation by the reaction force of the circular spline, and the power is output through the grooved pulley. This structure has both wheelbase and backlash adjustment functions.

To optimally ensure effective meshing between the oscillating teeth and the circular spline teeth, the tooth tops are appropriately shaped and trimmed to prevent potential interference during the meshing process. Due to the large number of oscillating teeth, the tooth extraction method is chosen, with each oscillating tooth being assigned three teeth to ensure optimal meshing strength and minimal motion vibration. The theoretical total number of oscillating teeth can be calculated by Eq.1:

$$Z_O=N\left(Z_A+Z_V\right)$$

(1)

In the equation, $N$ is the number of oscillating tooth blocks; $Z_A$ is the actual number of oscillating teeth; $Z_V$ is the number of teeth removed from the oscillating gear.

The difference between the number of circular spline teeth $Z_E$ and the theoretical total number of oscillating teeth $Z_O$ represents the number of dislocated teeth required for each rotation of the cam, while the number of cam waves represents the number of reciprocating movements of the oscillating teeth in a cycle. To ensure the correct meshing of the gear teeth, the two values must be guaranteed to be equal. In this design, it is a single-stage harmonic drive, and the transmission ratio can be expressed as Eq.2.

$$I=\frac{{\omega }_W}{{\omega }_g}=\frac{Z_E}{Z_E-Z_O}$$

(2)

The cam profile surfaces, which consist of two symmetrical helical surfaces, are integral to the determination of axial displacement of the oscillating teeth. In each cycle, the displacement of the oscillating teeth situated on the right-hand and left-hand helical surfaces of the cam can be expressed as Eq.3.

$$z=\left\{\begin{array}{ll}\frac{hnt}{30}\hfill & \left(0\le t\le \frac{30}{n}\right)\hfill \\ 2hn-\frac{hnt}{30}\hfill & \left(\frac{30}{n}\le t\le \frac{60}{n}\right)\hfill \end{array}\right.$$

(3)

Where $h$ is the cam lift, $n$ is the cam revolution speed, $t$ is the cam running time.

2.2 Analysis of Meshing Stiffness between Circular Spline and Oscillating Teeth

The oscillating tooth's force condition following its entry into the meshing phase is given. As depicted in Figure 2.

The meshing area and load point of the oscillating teeth are directly influenced by the displacement of the oscillating teeth. Given that the reaction force of the circular spline on the oscillating teeth is assumed to be an even load, the force generated by multiple teeth is equivalent to the concentrated load F of the center gear tooth. Based on the principles of force and torque equilibrium and input torque, the force equation of the oscillating teeth can be determined as Eq.4.

$$\left\{\begin{array}{l}{F}_{Ax}+F_{Bx}+F\cos \beta -N\sin {\theta }_W=0\hfill \\ F_{Ay}+F_{By}-F\sin \beta +N\cos {\theta }_W=0\hfill \\ F_{Bx}{L}_S-F_{By}{L}_g-F{L}_A-0.5N{L}_g\cos {\theta }_W=0\hfill \\ F_{Ay}{L}_g-F_{Ax}{L}_s-F{L}_B+0.5N{L}_g\cos {\theta }_W=0\hfill \\ 0.5F_{Ay}{L}_g\cos {\theta }_W-F_{Ax}\left(L_s+L_W\right)-F_{Bx}{L}_W-0.5F_{By}{L}_g\cos {\theta }_W-F{L}_C=0\hfill \end{array}\right.$$

(4)

In the formula, $L_S$ is the length of the groove; $L_W$ is the distance between the force point at the bottom of the oscillating tooth and the bottom of the groove; $L_M$ is the distance between the bottom of the oscillating tooth and the top of the tooth.

At the same time, the number of meshing teeth changes alternately with the cam angle. When the input torque is constant, the load on the oscillating teeth changes accordingly. The displacement of each oscillating tooth is denoted as $z_1$~ $z_7$, respectively. The meshing relationship between it and the gear teeth of the circular spline is shown in Figure 3.

The horizontal line in the figure represents the displacement required for the oscillating teeth to enter the meshing state. When the displacement of the oscillating teeth is higher than the horizontal line, the oscillating teeth enter the meshing. On the contrary, oscillating teeth are disengaged. That is the number of meshing teeth at the intersection of the horizontal line and the fold line changes.

The oscillating gear teeth can be simplified as a tapered beam on the flexible body of the oscillating teeth, and the load F is equivalent to the right end of the microelement, which is decomposed into force $F_x$ along the x-axis direction, force $F_y$ along the y-axis direction, and the equivalent bending M generated by F. As shown in Figure 4.

The effective length of the tapered beam is h, i.e. the distance between the base point N of the oscillating tooth and the top of the tooth; The gear tooth is divided into a series of rectangle micro-elements along the x-axis direction from the bottom of the tooth to the load point. Each microelement is represented by the symbol i, and its width is represented by the symbol $T_i$; F is the normal load on the oscillating tooth, $\beta$ is the angle between the load and the y-axis direction, $L_i$ is the distance between the microelement along the x-axis direction and the load point, $L_j$ is the distance between the load point along the x-axis direction and the tooth top, Y is the distance between the load point and the x-axis, $Y_N$ is the half-tooth width of the tooth root.

The calculation of deformation can be divided into three parts: bending deformation, shear deformation, and axial compression deformation of the oscillating tooth body; additional deformation caused by the elasticity of the oscillating tooth; and local contact deformation of the tooth at the meshing point of the oscillating tooth.

(1) Calculation of Deformation of Oscillating Tooth Body

Under load, the oscillating gear teeth produce axial compression deformation, shear deformation, and bending deformation along the equivalent deformation of the load direction. The deformation of a single microelement can be calculated and superimposed. It is assumed that the left end of each microelement i is fixed, and the portion connected to the right end of the microelement is regarded as a rigid body. The amount of compression deformation, shear deformation, and bending deformation of the microelement can be obtained as Eq.5.

$$\left\{\begin{array}{l}{\delta }_a={\displaystyle \sum _{i=1}^n{\delta }_{ai}}={\displaystyle \sum _{i=1}^n\frac{F_x{T}_i}{E_e{A}_i}}={\displaystyle \sum _{i=1}^n\frac{F{T}_i\sin \beta }{E_e{A}_i}}\\ {\delta }_s={\displaystyle \sum _{i=1}^n{\delta }_{si}}={\displaystyle \sum _{i=1}^n\frac{k{F}_y}{G{A}_i}}{T}_i={\displaystyle \sum _{i=1}^n\frac{12F(1+\nu )}{5E_e{A}_i}}{T}_i\cos \beta \\ {\delta }_b={\displaystyle \sum _{i=1}^n{\delta }_{bi}}={\displaystyle \sum _{i=1}^n\left({\omega }_1+{\theta }_1{l}_i+{\omega }_2+{\theta }_2{l}_i\right)}\end{array}\right.$$

(5)

In Eq.5, ${\delta }_a$, ${\delta }_s$, ${\delta }_b$ are the meshing point deformation caused by the compression, shear, and bending of the oscillating tooth body, and ${\delta }_{ai}$, ${\delta }_{si}$ and ${\delta }_{bi}$ are the deformation caused by the compression, shear, and bending of the microelement, respectively; $\nu$ is the Poisson's ratio of the material; $E_e$ is the equivalent elastic modulus of the tooth. According to Cornell's analysis, the ratio of the tooth width $b_i$ to the tooth thickness $s_p$, $b_i/s_p<5$, so it is a narrow tooth. At this time, the value of $E_e$ is the elastic modulus of the material; $A_i$ is the cross-sectional area of the gear teeth; ${\omega }_1$ is the deflection under the action of $F_y$, ${\theta }_1$ is the angle under the action of $F_y$, ${\omega }_2$ is the deflection under the action of M, ${\theta }_2$ is the angle under the action of M, respectively. As shown in Eq.6.

$$\left\{\begin{array}{l}{\omega }_1=\frac{F_y{T}_i^3}{3E_e{I}_i}=\frac{F{T}_i^3\cos \beta }{3E_e{I}_i}\\ {\theta }_1=\frac{F_y{T}_i^2}{2E_e{I}_i}=\frac{F{T}_i^2\cos \beta }{2E_e{I}_i}\\ {\omega }_2=\frac{M{T}_i^2}{2E_e{I}_i}=\frac{F{T}_i^2}{2E_e{I}_i}\left(L_i\cos \beta -Y\sin \beta \right)\\ {\theta }_2=\frac{M{T}_i^2}{E_e{I}_i}=\frac{F{T}_i}{E_e{I}_i}\left(L_i\cos \beta -Y\sin \beta \right)\end{array}\right.$$

(6)

In Eq.6, $I_i$ is the moment of inertia of the micro-element.

Substituting Eq.6 into Eq.5, results can be obtained as Eq.7.

$${\delta }_{bi}=\frac{F}{E_e{A}_i}{T}_i{\sin }^2\beta +\frac{12F(1+\nu )}{5E_e{A}_i}{T}_i{\cos }^2\beta +\frac{F{\cos }^2\beta }{E_e}\left(\frac{T_i^3+3T_i^3{L}_i+3T_i{L}_i^2}{3I_i}\right)-\frac{F\sin \beta \cos \beta }{E_e}\left(\frac{T_i^{2}Y+2T_{i}Y{L}_i}{2I_i}\right)$$

(7)

(2) Calculation of Oscillating Tooth Body Deformation

For narrow teeth, deal with the problem of plane stress as Eq.8.

$${\delta }_f=\frac{F{\cos }^2\beta }{Eb}\left[\frac{16.67}{\pi }{\left(\frac{L_f}{H_f}\right)}^2+2(1-v)\left(\frac{L_f}{H_f}\right)+1.534\left(1+\frac{{\tan }^2\beta }{2.4(1+v)}\right)\right]$$

(8)

In Eq.8, E is the elastic modulus of the material; $H_f$ is the tooth thickness at the tooth root N; $L_f$ is the equivalent arm of force.

(3) Calculation of Local Contact Deformation of Oscillating Teeth

The contact deformation of the meshing point of the oscillating tooth surface is caused by the contact and compression deformation of the gear meshing line. According to research by H. H. Lin and G. Lundberg, it can be expressed by Eq.9

$${\delta }_h=\frac{1.275F^{0.9}}{E_e^{0.9}{b}^{0.8}}$$

(9)

The total deformation of the meshing points of the movable teeth can be obtained by adding up each deformation as Eq.10.

$$\delta ={\delta }_a+{\delta }_s+{\delta }_b+{\delta }_f+{\delta }_h={\displaystyle \sum _{i=1}^n({\delta }_{ai}+{\delta }_{si}+{\delta }_{bi})}+{\delta }_f+{\delta }_h$$

(10)

When the normal load F is constant, when the oscillating tooth displacement $z>h_t-h$, the oscillating tooth and the circular spline tooth are meshing. As Eq.11.

$$\left\{\begin{array}{l}{L}_j=\left(z-h_t+h\right)/2\\ L_i=h-L_j-T_{i}i\\ Y=[L_j+h_t-h]\tan \beta \end{array}\right.$$

(11)

(4) Comprehensive Meshing Stiffness Calculation of Oscillating teeth

For a pair of intermeshing gear teeth, it can be regarded as a pair of springs in series. The compressive stiffness, shear stiffness, bending stiffness, deformation stiffness, and Hertz contact stiffness of oscillating gear teeth are expressed as Eq.12- Eq.16.

$$K_a=\frac{F}{{\delta }_a}={\displaystyle \sum _{i=1}^n\frac{E_e{A}_i}{T_i\sin \beta }}$$

(12)

$$K_s=\frac{F}{{\delta }_s}={\displaystyle \sum _{i=1}^n\frac{5E_e{A}_i}{12(1+v)T_i\cos \beta }}$$

(13)

$$K_b=\frac{F}{{\delta }_b}={\displaystyle \sum _{i=1}^n\left[\frac{T_i{\sin }^2\beta }{E_e{A}_i}+\frac{12(1+v)T_i{\cos }^2\beta }{5E_e{A}_i}+\frac{\left(T_i^3+3T_i^3{L}_i+3T_i{L}_i^2\right){\cos }^2\beta }{3E_e{I}_i}\right.}{\left.-\frac{\left(T_i^{2}Y+2T_{i}Y{L}_i\right)\sin \beta \cos \beta }{2E_e{I}_i}\right]}^{-1}$$

(14)

$$K_f=\frac{F}{{\delta }_f}=\frac{Eb}{{\cos }^2\beta }{\left[\frac{16.67}{\pi }{\left(\frac{L_f}{H_f}\right)}^2+2(1-v)\left(\frac{L_f}{H_f}\right)+1.534\left(1+\frac{{\tan }^2\beta }{2.4(1+v)}\right)\right]}^{-1}$$

(15)

$$K_h=\frac{F}{{\delta }_h}=\frac{E_e^{0.9}{b}^{0.8}}{1.275F^{0.1}}$$

(16)

In summary, the meshing stiffness of a single pair of gears at the load point can be obtained as Eq.17.

$$K=\frac{K_1{K}_2}{K_1+K_2}=1/\left(\frac{1}{K_{a1}}+\frac{1}{K_{s1}}+\frac{1}{K_{b1}}+\frac{1}{K_{f1}}+\frac{1}{K_h}+\frac{1}{K_{a2}}+\frac{1}{K_{s2}}+\frac{1}{K_{b2}}+\frac{1}{K_{f2}}\right)$$

(17)

In the formula, the meshing stiffness of the driving gear and the driven gear, namely the oscillating gear teeth and the circular spline teeth, at the load point, is a function of the position of the meshing point.

In the equation, $K_1$ and $K_2$ represent the mesh stiffness of the driving and driven gears, respectively. They are functions of the meshing point position, specifically the contact point between the oscillating teeth and the circular spline teeth.

Through the self-compiled software for example analysis and drawing curves, the stiffness of each part of the first oscillating tooth can be obtained as shown in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9.

The comprehensive stiffness is shown in Figure 10.

The remaining oscillating teeth have the same laws and only differ in phase differences. It can be obtained from the curve that the stiffness increases gradually with the oscillating tooth gradually entering the meshing, and decreases gradually with the oscillating tooth gradually withdrawing from the meshing; The bending stiffness, shear stiffness, and base stiffness are positively correlated with load angle, while compression stiffness is vice versa. The comprehensive stiffness of the oscillating tooth mesh can be improved by increasing $\beta$.

2.3 Contact Stiffness Analysis

The universal ball bearing is used between the oscillating teeth and the cam, and its contact point expands into a contact surface under the action of load N. The contact surface is projected onto the vertical surface of the contact normal, as shown in Figure 11-1. In the contact area between the universal ball and the circular spline, the contact stress is distributed in the semi-ellipsoid, as shown in Figure 11-2.

The contact between the cam and the universal ball is deformed as Eq.18.

$${\delta }_{13}=\frac{2K(e)}{\pi m_a}\sqrt[3]{\frac{1}{8}{\left(\frac{3}{E^{\prime }}\right)}^2{N}^2{\displaystyle \sum \rho }}$$

(18)

In Eq.18, N is the normal pressure at the contact point between the cam and the universal ball, and its direction is perpendicular to the cam profile; $E^{\text{'}}$ is the sum of the equivalent elastic modulus; $\Sigma \rho$ is the main curvature; ma is the long half-axis coefficient of the contact ellipse; $K\left(e\right)$ is the first type of complete elliptic integral related to the eccentricity of the ellipse e. The contact stiffness between the cam and the oscillating tooth can be obtained as Eq.19:

$$k_{13}=N/{\delta }_{13}$$

(19)

For the contact problem between the grooved pulley and the oscillating teeth, since the mass of the grooved pulley is significantly larger than the mass of the oscillating teeth, the grooved pulley can be regarded as a rigid body. The oscillating tooth is squeezed by the groove pulley during the meshing process, and the expression of the compressive stiffness between the groove pulley and the oscillating tooth can be obtained as Eq.20:

$$k_{23}=E_3{A}_3/b_3$$

(20)

In Eq.20, $b_3$ is the width of the oscillating tooth rod; $E_3$ is the elastic modulus of the oscillating tooth material; $A_3$ is the contact area between the oscillating tooth rod and the groove pulley.

2.4 Stiffness Analysis of Support System

Assuming the radial stiffness of the bearing is isotropic since the input shaft is supported by a single bearing, the radial support stiffness of the cam is the radial support stiffness of the bearing at the input end, and its value can be calculated by Eq.21.

$$k_{1r}=F_{b1}/({\delta }_{b11}+{\delta }_{b12}+{\delta }_{b13})$$

(21)

In Eq.21, $F_{b1}$ is the radial load on the bearing at the input end; ${\delta }_{b11}$ is the radial elastic displacement of the bearing at the input end; ${\delta }_{b12}$ is the contact deformation between the outer ring of the bearing at the input end and the box hole; ${\delta }_{b13}$ is the contact deformation between the inner ring of the bearing at the input end and the shaft diameter. The three expressions are respectively expressed by Eq.22, Eq.23, and Eq.24.

$${\delta }_{b11}={\beta }_{b1}{\delta }_{01}$$

(22)

$${\delta }_{b12}=H_{b11}{\Delta }_1$$

(23)

$${\delta }_{b13}=\frac{0.204F_{b1}{H}_{b12}}{\pi b_{b1}{d}_{b1}}$$

(24)

Where, ${\beta }_{b1}$ is the elastic displacement coefficient of the bearing at the input end, which is found from the standard according to the relative clearance $g_{b1}/{\delta }_{01}$; ${\delta }_{01}$ is the radial elastic displacement when the clearance in the bearing at the input end is zero; $g_{b1}$ is the clearance or preload in the bearing at the input end; ${\Delta }_1$ is the fit clearance in the diameter direction between the outer ring of the bearing at the input end and the inner hole of the frame at the input end; $H_{b11}$ is the elastic coefficient of the bearing at the input end; $F_{b1}$ is the radial load; $H_{b12}$ is the deformation coefficient; $b_{b1}$ is the width of the bearing ring at the input end; $d_{b1}$ is the inner diameter of the bearing at the input end.

Similarly, the radial stiffness of the bearing at the output end can be obtained by Eq.25.

$$K_{2r}=F_{b2}/({\delta }_{b21}+{\delta }_{b22}+{\delta }_{b23})$$

(25)

Since a pair of bearings are installed on the output shaft, the comprehensive radial support stiffness of the output terminal can be obtained by the method of spring parallel connection as Eq.26.

$$k_{2r}=K_{2r}/2$$

(26)

The two shaft structures are both stepped shafts, and the torque provided by the motor or the load is received during the movement. The segmentation and torque are shown in Figure 12-1 and Figure 12-2.

For the input shaft, the convex is simplified to a cylinder of equal width, ignoring the thin-walled support part between its inner diameter and the hub, and being processed in parallel with the matching shaft section, the total torsion angle of the input shaft can be obtained by Eq.27.

$${\phi }_1={\displaystyle \sum _{i=1}^n{\phi }_{1i}}=\frac{32T_1{l}_{11}}{G\pi d_{11}^4}+\frac{32T_1{l}_{12}}{G\pi d_{12}^4}+\frac{32T_1{l}_{13}}{\pi G(d_{13}^4+d_{15}^4(1-{(d_{14}/d_{15})}^4))}$$

(27)

In Eq.27, ${\phi }_{1i}$ is the torsion angle of the i-th shaft section of the input shaft; $T_1$ is the input torque; G is the shear modulus of the input shaft material, $l_{11}$, $l_{12}$, $l_{13}$ are the lengths of each shaft section of the input shaft; $d_{11}$, $d_{12}$, $d_{13}$ are the shaft diameters of each shaft section of the input shaft; $d_{14}$ and $d_{15}$ are the inner and outer diameters of the cam, respectively.

For the output shaft, the sheave is simplified into a cylinder of equal width, and its matching shaft section is regarded as a solid shaft as a whole, and the total torsion angle of the output shaft can be obtained by Eq.28.

$${\phi }_2={\displaystyle \sum _{i=1}^n{\phi }_{2i}}={\displaystyle \sum _{i=1}^n\frac{32T_2{l}_{2i}}{G\pi d_{2i}^4}}$$

(28)

In Eq.28, $T_2$ is the load torque; G is the shear modulus of the output shaft material; $l_{2i}$ is the length of the i-th shaft section of the output shaft; $d_{2i}$ is the shaft diameter of the i-th shaft section of the output shaft.

2.5 Other Excitation Parameters

The mass and moment of inertia of each component are shown in Table 1.

3. System Dynamics Analysis

3.1 System Dynamics Model

The following assumptions are made in the analysis:

(1) During the motion process of the transmission system, all transmission components remain on the same axis.

(2) Components such as the box are considered rigid bodies, and the contact between oscillating teeth and end-face cams, as well as that between the rigid circular spline and grooved pulley, is modeled as a spring-damping system.

(3) Additional vibration phenomena arising from factors such as assembly and transmission errors among components are not considered.

(4) The mass properties of each oscillating tooth are the same.

The relevant dynamic parameters are set in Table 2.

Due to the dynamic behavior of each oscillating tooth being entirely identical, there is a phase difference in kinematics. Therefore, only the case of a single oscillating tooth needs to be considered.

The relevant space coordinate systems are established as follows: $O-xyz$ is the fixed space coordinate system built on the rigid circular spline. $O-x_1{y}_1{z}_1$ is the following space coordinate system rotating with the cam. $O-x_2{y}_2{z}_2$ is the following space coordinate system rotating with the grooved pulley. $O-x_3{y}_3{z}_3$ the following space coordinate system rotating with the (i)-th oscillating tooth.

The dynamic model of the transmission system is illustrated in Figure 13.

Where the relevant parameters are listed in Table 3.

Taking the direction of the meshing point between the cam and the oscillating tooth toward the oscillating tooth as the positive direction, the projection of the cam displacement relative to the oscillating tooth along the direction of the mesh line can be obtained as Eq.29.

$${\delta }_{13i}=\left\{\begin{array}{l}-\frac{h\cos {\theta }_W}{\pi r_1}\left[x_1\sin {\tau }_i+y_1\cos {\tau }_i-x_{3i}+{\phi }_1{r}_1\right]-(z_{3i}+z_1)\cos {\theta }_W\left(0\le {\tau }_i<\pi \right)\\ \frac{h\cos {\theta }_W}{\pi r_1}\left[x_1\sin {\tau }_i+y_1\cos {\tau }_i-x_{3i}+{\phi }_1{r}_1\right]-(z_{3i}+z_1)\cos {\theta }_W\left(\pi \le {\tau }_i<2\pi \right)\end{array}\right.$$

(29)

Taking the direction of the meshing point between the grooved pulley and the oscillating tooth toward the oscillating tooth as the positive direction, the projection of the grooved pulley displacement relative to the oscillating tooth along the direction of the mesh line can be obtained as Eq.30.

$${\delta }_{23i}=\left\{\begin{array}{l}-x_2\sin \left({\tau }_i+{\phi }_W-{\phi }_G\right)-y_2\cos \left({\tau }_i+{\phi }_W-{\phi }_G\right)+x_{3i}+{\phi }_2{r}_2\left(0\le {\tau }_i<\pi \right)\\ x_2\sin \left({\tau }_i+{\phi }_W-{\phi }_G\right)-y_2\cos \left({\tau }_i+{\phi }_W-{\phi }_G\right)+x_{3i}+{\phi }_2{r}_2\left(\pi \le {\tau }_i<2\pi \right)\end{array}\right.$$

(30)

Taking the direction of the meshing point between the rigid circular spline and the oscillating tooth toward the oscillating tooth as the positive direction, the projection of the rigid circular spline displacement relative to the oscillating tooth along the direction of the mesh line can be obtained as Eq.31.

$${\delta }_{34i}=\left\{\begin{array}{l}-x_{3i}\cos {\theta }_E\cos \beta +y_{3i}\sin {\theta }_E\cos \beta -z_{3i}\sin \beta （0\le {\tau }_i<\pi ）\\ x_{3i}\cos {\theta }_E\cos \beta +y_{3i}\sin {\theta }_E\cos \beta -z_{3i}\sin \beta （\pi \le {\tau }_i<2\pi ）\end{array}\right.$$

(31)

3.2 System Dynamics Differential Equations

According to the system dynamics model, dynamics differential equations of transmission components can be obtained as Eq.32, Eq.34, and Eq.36.

For the cam:

$$\left\{\begin{array}{l}{m}_1{\ddot{x}}_1+c_{1r}{\dot{x}}_1+{\displaystyle \sum _{i=1}^z{c}_{13}}{\dot{x}}_{13i}+k_{1r}{x}_1+{\displaystyle \sum _{i=1}^z{k}_{13}}{x}_{13i}=0\hfill \\ m_1{\ddot{y}}_1+c_{1r}{\dot{y}}_1+{\displaystyle \sum _{i=1}^z{c}_{13}}{\dot{y}}_{13i}+k_{1r}{y}_1+{\displaystyle \sum _{i=1}^z{k}_{13}}{y}_{13i}=0\hfill \\ m_1{\ddot{z}}_1+c_{1a}{\dot{z}}_1+{\displaystyle \sum _{i=1}^z{c}_{13}}{\dot{z}}_{13i}+k_{1a}{z}_1+{\displaystyle \sum _{i=1}^z{k}_{13}}{z}_{13i}=0\hfill \\ J_1{\ddot{\phi }}_1+c_{1t}{\dot{\phi }}_1/r_1+k_{1t}{\phi }_1/r_1=\frac{T_1}{r_1}\hfill \end{array}\right.$$

(32)

Where $T_1$ is the input torque; $x_{13i}$, $y_{13i}$, $z_{13i}$ are the projections of ${\delta }_{13i}$ in the cam following coordinate system $O-x_1{y}_1{z}_1$ along the $x_1$, $y_1$, $z_1$ directions, respectively, and which can be expressed as Eq.33.

$$\left\{\begin{array}{l}{x}_{13i}=\pm {\delta }_{13i}\sin {\theta }_W\sin {\tau }_i\hfill \\ y_{13i}=\pm {\delta }_{13i}\sin {\theta }_W\cos {\tau }_i\hfill \\ z_{13i}={\delta }_{13i}\cos {\theta }_W\hfill \end{array}\right.$$

(33)

For the grooved pulley:

$$\left\{\begin{array}{l}{m}_2{\ddot{x}}_2+c_{2r}{\dot{x}}_2+{\displaystyle \sum _{i=1}^z{c}_{23}}{\dot{x}}_{23i}+k_{2r}{x}_2+{\displaystyle \sum _{i=1}^z{k}_{23}}{x}_{23i}=0\hfill \\ m_2{\ddot{y}}_2+c_{2r}{\dot{y}}_2+{\displaystyle \sum _{i=1}^z{c}_{23}}{\dot{y}}_{23i}+k_{2r}{y}_2+{\displaystyle \sum _{i=1}^z{k}_{23}}{y}_{23i}=0\hfill \\ m_2{\ddot{z}}_2+c_{2a}{\dot{z}}_2+{\displaystyle \sum _{i=1}^z{c}_{23}}{\dot{z}}_{23i}+k_{2a}{z}_2+{\displaystyle \sum _{i=1}^z{k}_{23}}{z}_{23i}=0\hfill \\ J_2{\ddot{\phi }}_2+c_{2t}{\dot{\phi }}_2/r_2+k_{2t}{\phi }_2/r_2=-\frac{T_2}{r_2}\hfill \end{array}\right.$$

(34)

Where $T_2$ is the output torque; $x_{23i}$, $y_{23i}$, $z_{23i}$ are the projections of ${\delta }_{23i}$ in the grooved pulley following coordinate system $O-x_2{y}_2{z}_2$ along the $x_2$, $y_2$, $z_2$ directions, respectively, and which can be expressed as Eq.35.

$$\left\{\begin{array}{l}{x}_{23i}={\delta }_{23i}\sin \left({\tau }_i+{\phi }_W-{\phi }_G\right)\hfill \\ \begin{array}{l}{y}_{23i}={\delta }_{23i}\cos \left({\tau }_i+{\phi }_W-{\phi }_G\right)\\ z_{23i}=0\end{array}\hfill \end{array}\right.$$

(35)

For the (i)-th oscillating tooth:

$$\left\{\begin{array}{l}{m}_3{\ddot{x}}_{3i}+c_{13}{\dot{x}}_{31i}+c_{23}{\dot{x}}_{32i}+c_{34}{\dot{x}}_{34i}+k_{13}{x}_{31i}+k_{23}{x}_{32i}+k_{34i}f(x_{34i})=0\hfill \\ m_3{\ddot{y}}_{3i}+c_{13}{\dot{y}}_{31i}+c_{23}{\dot{y}}_{32i}+c_{34}{\dot{y}}_{34i}+k_{13}{y}_{31i}+k_{23}{y}_{32i}+k_{34i}f(y_{34i})=0\hfill \\ m_3{\ddot{z}}_{3i}+c_{13}{\dot{z}}_{31i}+c_{34}{\dot{z}}_{34i}+k_{13}{z}_{31i}+k_{34i}f(z_{34i})=0\hfill \end{array}\right.$$

(36)

Where $x_{31i}$, $y_{31i}$, $z_{31i}$ are the projections of ${\delta }_{31i}$ in the (i)-th oscillating tooth following coordinate system $O-x_3{y}_3{z}_3$ along the coordinate axis, $x_{32i}$, $y_{32i}$, $z_{32i}$ are the projections of ${\delta }_{32i}$ in the (i)-th oscillating tooth following coordinate system $O-x_3{y}_3{z}_3$ along the coordinate axis, $x_{34i}$, $y_{34i}$, $z_{34i}$ are the projections of ${\delta }_{34i}$ in the (i)-th oscillating tooth following coordinate system $O-x_3{y}_3{z}_3$ along the coordinate axis, respectively. These above can be expressed as Eq.37.

$$\left\{\begin{array}{l}{x}_{31i}={\delta }_{13i}\sin {\theta }_W,y_{31i}=0,z_{31i}={\delta }_{13i}\cos {\theta }_W\hfill \\ x_{32i}={\delta }_{23i},y_{32i}=0,z_{32i}=0\hfill \\ x_{34i}=\pm {\delta }_{34i}\cos {\theta }_E\cos \beta ,y_{34i}={\delta }_{34i}\sin {\theta }_E\cos \beta ,z_{34i}=-{\delta }_{34i}\cos {\theta }_E\sin \beta \hfill \end{array}\right.$$

(37)

$f\left(x_{34i}\right)$, $f\left(y_{34i}\right)$ and $f\left(z_{34i}\right)$ are the nonlinear functions of tooth side clearance along the x-axis, y-axis, and z-axis directions, respectively. Assuming that the tooth side clearance is 2b, the displacement of the oscillating tooth along the tooth surface in the normal direction is $x_n$. These above can be expressed as (Eq.38, Eq.39, Eq40.), respectively.

$$f(x_{34i},b)=\left\{\begin{array}{cc}{x}_{34i}-b\cos {\theta }_E\cos \beta & (x_n>b)\\ 0& (-b\le x_n\le b)\\ x_{34i}+b\cos {\theta }_E\cos \beta & (x_n<-b)\end{array}\right.$$

(38)

$$f(y_{34i},b)=\left\{\begin{array}{cc}{y}_{34i}-b\sin {\theta }_E\cos \beta & (x_n>b)\\ 0& (-b\le x_n\le b)\\ y_{34i}+b\sin {\theta }_E\cos \beta & (x_n<-b)\end{array}\right.$$

(39)

$$f(z_{34i},b)=\left\{\begin{array}{cc}{z}_{34i}-b\cos {\theta }_E\sin \beta & (x_n>b)\\ 0& (-b\le x_n\le b)\\ z_{34i}+b\cos {\theta }_E\sin \beta & (x_n<-b)\end{array}\right.$$

(40)

From the above, the dynamics differential equations of the end cam input mechanism, oscillating teeth, and grooved pulley output mechanism have been obtained. Organizing and sequencing the equations of each component and the overall dynamic differential equations of the system in matrix form is obtained as Eq.41.

$$M\ddot{X}+C\dot{X}+KX=T$$

(41)

Where $M$ is the system generalized quality matrix; $X$ is the system generalized displacement matrix; $C$ is the system generalized damping matrix; $K$ is the system generalized rigidity matrix; $T$ is the system externally excited array matrix.

3.3 System Vibration Response Analysis

The nonlinear dynamics differential equations can be solved by using the Runger-Kutta method. Setting speed and load of the cam in three different states, and setting damping as Rayleigh damping. The transient vibration displacement response of the cam along each independent coordinate is obtained as shown in Figure 14, Figure 15, Figure 16 and Figure 17, respectively.

.From Figure 14, Figure 15, Figure 16 and Figure 17, the results show that the cam transverse displacement response is slightly larger than the radial displacement response. The reason is that the cam profile surface lift angle and return angle are along the transverse direction, there is a fixed angle of contact between the cam and the oscillating teeth, and most of the vibration is transmitted to the transverse direction of the cam.

The transient vibration displacement response of the grooved pulley along each independent coordinate is shown in Figure 18, Figure 19, Figure 20 and Figure 21, respectively.

.From Figure 18, Figure 19, Figure 20 and Figure 21, the results show that the grooved pulley transverse displacement response is slightly less than the radial displacement response and axial displacement response. The reason is that during the oscillating teeth meshing with the teeth of the rigid circular spline, there exists a central angle on the tooth side pointing to the center of rotation and a tooth angle on the tooth surface along the direction of the axis, which makes the vibration of the oscillating gear more generated in the radial direction of the grooved pulley and transmitted to the grooved pulley.

The vibration patterns between pairs of dissimilar oscillating teeth are the same and there is only a phase difference. Taking the first pair of oscillating teeth as an example, their transient vibration displacement responses along each independent coordinate are shown in Figure 22, Figure 23, and Figure 24, respectively.

.From Figure 22, Figure 23 and Figure 24, the results show that the transverse displacement of the oscillating teeth is the largest. The reason is that the setting of the angle between the oscillating teeth and the teeth of the rigid circular spline results in the line of action of the load being more directed in the transverse direction of the oscillating teeth. Therefore, the oscillating gear teeth are subjected to greater deformation and vibration response along the transverse direction on entering meshing.

The cam, the grooved pulley, and the oscillating teeth all generate regular periodic vibrations in the transverse, radial, and axial directions at their equilibrium positions. When the oscillating teeth enter into and disengage from engagement, the vibration response changes abruptly. With the increase of the input speed of the single-wave cam and the load borne by the output end, the vibration amplitude of the oscillating teeth in all directions increases gradually.

4. Dynamic behavioral validation

To validate the results of the nonlinear dynamics study of the end-face harmonic transmission system, building the test bench as shown in Figure 25. 220V AC power is stepped down by a transformer to supply power to the stepper motor, torque sensor, and magnetic particle brake. The inputs control the stepper motors for power input through the driver, and the development board used is powered by 5V DC. The torque sensor is mounted on the input side of the end-meshing harmonic reducer to measure speed and torque, which is then read by the dynamic torsion triple display meter. The encoder is mounted to the output side to measure the real-time rotational speed of the output shaft, and the signal conversion is accomplished by the debugger. Magnetic particle brakes are mounted to the end to provide load and are adjusted in size by the load controller. The resulting sensing signals are transmitted to the PC for data processing.

Measured at standard center distance under follows conditions: 88r/min, 0.65N·m; 88r/min, 0.98N·m and 220r/min, 0.65N·m. Reducing the center distance by 0.5mm and 1mm respectively, the date under 88r/min, 0.65N·m was measured. These above were recorded as Working Conditions (1-5). After data processing, the comparison of results between the output shaft circumferential torsional vibration response and the simulation results is shown in Figure 26, Figure 27, Figure 28, Figure 29 and Figure 30.

From Figure 26, Figure 27, Figure 28, Figure 29 and Figure 30, the results show that the theoretical and experimental curves are in general agreement in terms of trend and are numerically closer, which verifies the correctness of the theoretical model. In addition to the influence of the load at the output end and the rotational speed at the input end, the increase of the center distance makes the support position of the rod part of the oscillating teeth gear system change, which leads to the reduction of the gear teeth meshing rigidity and the vibration increases of the transmission system. According to the analysis of the meshing stiffness above, it can be seen that the increase in pressure angle makes the meshing stiffness increase, which can reduce the vibration of the system.

The main reasons for the error between the experimental and simulation data for the vibration response are as follows:

(1) Errors in the machining of components;

(2) Clearance between the rod part of the oscillating teeth gear system and the moving part of the grooved pulley.

(3) The box and bench support are not absolutely rigid.

(4) The sensors are subject to a certain amount of deviation during the measurement process with the vibration of the test bench.

5. Conclusions

To study the coupling mechanism and dynamic response of the end-face movable gear transmission system under complex excitation, a specific configuration of end-meshing movable gear reduction mechanism was used to achieve predetermined rigid thrust transmission and mismatched gear meshing functions, which solved the inherent defects of traditional harmonic gear mechanism thin-wall flexible wheel easy to be damaged by fatigue. The following are the main conclusions:

1. We have completed the development of a physical prototype of end-meshing harmonic reducer. Then, considering the effects of internal and external excitation parameters such as time-varying meshing stiffness, Hertzian contact stiffness, and gear tooth meshing side gap, we established the dynamic model of the end-meshing harmonic reducer, derived the vibration displacement response of each transmission member, and analyzed its influencing factors.

2. We built a comprehensive test bench of adjustable characteristics to realize the real-time transmission and acquisition of output shaft circumferential torsional vibration response data. Compared with the theoretical data, the correctness of the dynamics model of the end-meshing harmonic reducer was verified and the reasons for the errors were analyzed.

3. The results show that the magnitude of input speed and load torque has a large influence on the system vibration response. For the internal excitation characteristics, measures such as increasing the support stiffness of transmission components, shortening the center distance, increasing the meshing stiffness of the gear teeth, decreasing the length of the rod part of the oscillating teeth gear system, and increasing the tooth angle can produce vibration suppression effects.

6. Patents

This section is not mandatory but may be added if there are patents resulting from the work reported in this manuscript.