A Refined Dynamic Model of Harmonic Drive and Its Dynamic Response Analysis

To highlight the key factors which inﬂuence the dynamic performance of the harmonic drive, a reﬁned harmonic drive model considering nonlinear stiﬀness, kinematic error, and friction of the critical components is established. A dedicated experimental apparatus based on double motor twisting is constructed to measure the characteristics of harmonic drive, and the attribute parameters of the proposed model are identiﬁed. A series of experiments on the dynamic transmission error at diﬀerent driving velocities are carried out to verify the proposed model. Based on the proposed model, the inﬂuence of diﬀerent component stiﬀness on the velocity step response of the harmonic drive is analyzed. The results show that the inﬂuence of the component stiﬀness on the system dynamic response is more signiﬁcant at high driving velocity, the increase of the stiﬀness of each component will decrease the dynamic transmission accuracy of harmonic drive, and the bearing radial stiﬀness is the most sensitive parameter to system’s dynamic response among all the stiﬀness factors.


Introduction
Since harmonic drive was invented in the 1950s [1], it has been widely used in various types of important electromechanical equipment such as robots, machine tools, radars, satellites, space rovers, spatial manipulators, and other weapon systems because of their desirable characteristics of near-zero backlash, light weight, compactness, and high gear ratio [2]. In the meantime, as a critical transmission component, some unexpected dynamic behaviors of the harmonic drive are presented in the previous literature. For example, Taghirad and Belanger [3] demonstrated that harmonic drive systems suffer from high flexibility, resonance vibration, friction, and structural damping nonlinearities. Yamamoto et al. [4] showed that the nonlinear, elastic attribute of the harmonic drive might deteriorate the static control accuracy due to the transmission error with hysteresis. Zhu et al. [5] designated the high friction and the dynamics of the flexspline in harmonic drive as the main issues that significantly challenge the control systems. erefore, it is essential to have a detailed understanding of the dynamic behaviors in harmonic drive.
A harmonic drive consists of the three main components identified in Figure 1: flexspline, circular spline, and wave generator. e transmission form is mainly a reduction drive; that is, the wave generator is connected to the motor axis and the flexspline is connected to the output, while the circular spline is fixed. e elliptical wave generator which rotates within the flexspline deflects the flexspline slightly from its natural circular form into an elliptical shape continuously. Simultaneously, the teeth on the flexspline will mesh with that of the circular spline at regions near the major axis of the wave generator, which would create a counterrotation at the flexspline output with respect to the wave generator. e counterrotation of flexspline is a combined effect of meshing between flexspline and circular spline and pushing driven by the wave generator, as demonstrated in literature [6]. Since the teeth of the flexspline are usually two fewer than that of the circular spline, for a full revolution of the wave generator, the flexspline will only counterrotate a very small angle of two teeth, so the harmonic drive generates a high transmission ratio.
Initiated in 1985, Good et al. [7] noted the effect of nonlinear friction in harmonic drive on system transmission performance, while Aliev [8] studied the dynamic behavior of flexible gear in harmonic drive in 1986. en, a substantial research on dynamic behaviors of harmonic drive is followed. Seyfferth et al. [9,10] developed a mechanical model that considers compliance and friction in the harmonic drive, the Coulomb-type friction at the input shaft, and the viscous-type damping, and the model parameters were identified by a series of quasistatic and stationary dynamic experiments. Similarly, Taghirad and Belanger [3] proposed a model for compliance, hysteresis, and friction in harmonic drive; Kircanski and Goldenberg [11] modeled the nonlinear stiffness, hysteresis, and friction in robot joints with harmonic drives. All necessary parameters in the models above were identified by experimental results. Since the harmonic drive has many dynamic attributes, some researches focused on a specific attribute. For example, Gandhi et al. [12], Han et al. [13], and Liao et al. [14] studied the modeling, identification, and compensation of friction in harmonic drive. Ghorbel et al. [15] and Yamamoto et al. [4] investigated the kinematic error and transmission error in the harmonic drive, respectively. Dhaouadi et al. [16] proposed a new dynamic model of hysteresis in harmonic drive using the heredity concept. Moreover, Tjahjowidodo et al. [17] and Rhéaume et al. [18] studied the nonlinear modeling of torsional behavior and torsional stiffness in the harmonic drive, respectively. In 2011, Preissner et al. [19] presented a comprehensive harmonic drive model, which included four dominant aspects of harmonic drive behavior: nonlinear viscous friction, nonlinear stiffness, hysteresis, and kinematic error. Many researchers [5,[20][21][22][23][24][25][26][27] established the dynamic models of different transmission systems including harmonic drives. Based on the research above, most mentioned models described the harmonic drive as a gray or black model with certain dynamic attributes. High attention has been paid to the input and output of harmonic drive rather than the three interacting components, and the identified attributes of the black model are the dynamic attributes of the whole harmonic drive. us, the influence of the three components on the dynamic performance of harmonic drive is not fully investigated. Moreover, the adjustments of the three components in the design and manufacturing process cannot be implemented in the harmonic drive to achieve a better dynamic performance. erefore, for a detailed understanding of the dynamic behaviors, a refined dynamic model that considers the specific attributes of the three components in harmonic drive is needed.
Tuttle and Seering [28,29] proposed an intuitive model based on an equivalent gear pair which could be used to explain the transmission mechanism of harmonic drive, but the attributes of the model are still integrated expressions of the whole harmonic drive system. Zhang et al. [30] modeled the compliance of the flexspline and the wave generator instead of the overall harmonic drive transmission, while the friction of the flexspline and the wave generator was not included.
rough FEM, Dong et al. [31] established a harmonic drive dynamic model. However, the modeling is complicated, and it is hard to solve the FEM model because of the difficulty of convergence. Zou et al. [6] presented a more complex model of harmonic drive based on equivalent gear pair, but the evaluation of the meshing stiffness needed too many structural parameters, and the influence of the components' attributes on the system dynamic response has not been analyzed. erefore, up to now, the existing models are not refined enough to explain the dynamic behavior of harmonic drive, so that the influence of the components' parameters on the system dynamic response is not clear. ese models cannot be used to modify the design and manufacturing process for a better dynamic performance of the harmonic drive.
Given this, this paper establishes a more refined, new harmonic drive model that includes all of the nonlinear stiffness, kinematic error, and friction of the key components in the harmonic drive. A specialized experimental apparatus based on double motor twisting is constructed to measure the characteristics of harmonic drive, and the attribute parameters are identified. A series of experiments to evaluate the dynamic transmission error of the model at different velocities are carried out, and the influencing rules of the components' attributes on the harmonic drive dynamic response are researched. e results of dynamic response analysis will be useful to the dynamic design and manufacturing of high-quality and high-performance harmonic drives.

Harmonic Drive Modeling
e structural relation and motion form of harmonic drive components are shown in Figure 2. e wave generator of harmonic drive is composed of an elliptical cam and a thin-wall ball bearing. When harmonic drive works, the elliptical cam rotates causing the radial deflection of the bearing. en the radial motion of bearing forces the flexspline to engage into the circular spline in the radial direction. Due to the meshing effect, the flexspline teeth move along the tangential direction which causes the output rotation of harmonic drive. According to the mechanism analysis above, the compliance of the bearing and the meshing compliance are important factors influencing the transmission of harmonic drive. According to Figure 1, as the flexspline is a flexible thin-wall cylinder, the torsional compliance of the cylinder also should be considered. Also, when the harmonic drive works, the friction in the bearing and that in the meshing teeth pairs influence the smooth transmission of the harmonic drive, so the frictions should be evaluated in the model. Additionally, the synthetic eccentric error caused by manufacturing and assembly errors can lead to kinematic error, which is the main excitation source of the harmonic drive dynamic performance according to [6,15,29]. us, the compliance, friction, and error factors above are all important attributes which may cause system kinematic inaccuracy, energy dissipation, and velocity fluctuation. Since they can influence the dynamic behavior or reduce the dynamic transmission accuracy of the harmonic drive, they need to be measured, evaluated, and included in the proposed model.
A refined harmonic drive model considering the components' attributes is established (see Figure 3). In Figure 3, the input axis rotates, and the elliptical cam moves. en the elliptical cam motion forces the bearing to produce radial motion (x 1 ). rough the bearing, x 1 makes the flexspline teeth move x 2 along the radial direction to engage into the circular spline teeth. e meshing effect produces the tangential motion of flexspline (y 1 ) which, in turn, generates the tangential motion of output axis (y 2 ) through the torsion of the flexspline cylinder.
e term e is the equivalent synthetic eccentric error that leads to the pure kinematic error, K b and C b are radial stiffness and damping of the bearing, F fb is friction in the bearing, K mj is meshing stiffness of the jth gear pair, F fm is the meshing friction, and K t is torsional stiffness of the flexspline cylinder. θ i , θ o , r g , α t , J g , and J o are input rotation angle, output rotation angle, radius of flex spline cylinder, equivalent gear tooth angle (constant), the inertia of the flexspline, and load inertia of the output axis, respectively. α n is angle of the wedge standing for the elliptical cam; it ensures that the harmonic drive model is merely an ideal reducer model in the absence of any friction, compliance, or kinematic error. α n can be written as where N is gear ratio of the harmonic drive.
As Figure 2 shows, the engagement extents of the gear teeth are different at different positions, so the meshing stiffness of different gear pairs is also not the same. Moreover, with the increase of the load at the output of harmonic drive, more flexspline teeth will touch the circular spline teeth and more gear pairs will start to mesh, making the meshing stiffness change nonlinearly. As it is hard to evaluate the stiffness of each gear pair, an integrated meshing stiffness K m is used to represent the nonlinear meshing stiffness as en the relationship between input axis rotation and radial motion of bearing can be formulated as According to the Newtonian equations of motion, the ordinary differential equations for the radial motion and force transmission process are Meanwhile, the ordinary differential equations for the tangential motion and force transmission process are Input axis Output axis F Flexspline teeth (Tangential motion of output axis) ...
Elliptical cam motion Circular spline teeth If the attribute parameters like stiffness, friction, and kinematic error are determined and incorporated into the above-mentioned equations and the input rotation angle is given, the output rotation angle of the harmonic drive then can be predicted.

Experimental Apparatus.
A dedicated experimental apparatus based on double motor twisting (Figures 4 and 5) was built. e driving servo motor with an encoder (whose resolution is 8,192 lines) was responsible for the input motion of harmonic drive. A harmonic drive of type CSF-25-120 manufactured by HD System Inc. was chosen as the measured harmonic drive. A HEIDENHAIN high-precision RON786 encoder (whose resolution is 18,000 lines) and an HBM T22 torque sensor were used to measure the output angle and output torque of the measured harmonic drive. e input angle and input torque of harmonic drive are measured through the servo driver. All angle and torque signals were acquired synchronously.
e loading system was composed of the torque sensor, twisted bar, attendant reducer, and loading servo motor. e attendant reducer, which was another harmonic drive with a larger gear ratio, amplified the loading torque. When the output end of the measured harmonic drive needed to be loaded, the driving servo motor could be controlled to pause at home position, and the loading servo motor was controlled to rotate until the output of harmonic drive was loaded with target torque. In the case of no loading of the measured harmonic drive, the coupling between the high-precision encoder and the torque sensor could be dismantled.

Stiffness.
When the stiffness of harmonic drive was measured, the driving servo motor paused, and the loading servo motor rotated back and forth. Meanwhile, the output torsional angle and the output torque were acquired. e stiffness curve of the harmonic drive was obtained, and it is shown in Figure 6. e relationship between the output torsional angle and the output torque could be described as follows: where Δθ o is output torsional angle, T o is output torque, and g (1) and g (2) are nonlinear, unknown stiffness coefficients that could be identified through the least square method. e simulation results of the nonlinear stiffness model are also shown in Figure 6. Since nonlinear stiffness and friction dissipation both result in the stiffness hysteresis loop and the friction dissipation is not included in (7), the hysteresis is not expressed in the simulation data compared with the experimental data.
In the harmonic drive, there are three stiffness terms: meshing stiffness K m , flexspline cylinder torsional stiffness K t , and radial stiffness of bearing K b . eir corresponding output torsional compliance components of the harmonic drive are f m , f t , and f b , respectively. e output torsional compliance f HD of the whole harmonic drive could be obtained as follows: K t could be calculated as a torsion problem of an ordinary cylinder and K b could be estimated according to the method of stiffness calculation for an ordinary ball bearing, which will not be discussed in this paper. eir corresponding output torsional compliance components f t and f b could be calculated through the structure of the proposed model as follows: en, according to (7)−(10), the compliance f m induced by gear meshing could be calculated as At last, the nonlinear meshing stiffness K m could be determined according to the structure of the proposed model as follows: According to (7)−(12), the nonlinear relationship between the meshing stiffness and the output torque K m � f(T o ) is obtained, and then the meshing stiffness K m could be adjusted according to T o . In the model, the output torque T o could be obtained as

Friction.
Friction in the harmonic drive is a velocityrelated attribute [29]; it has to be measured at different velocities. us, the driving servo motor was controlled at different constant velocities from 5∼3000 rpm. In every constant velocity, the input driving torque of harmonic drive was regarded as the friction torque and acquired through the servo driver. Figure 7 shows the input friction torque of harmonic drive in two directions at different velocities. e harmonic drive contains two main friction sources: the friction in bearing F fb and the meshing friction F fm . F fb could be modeled based on the classic bearing Palmgren empirical equation [32]; that is,    Shock and Vibration where f 0 is a coefficient determined by the bearing and lubrication style, υ 0 is viscosity of the lubricant, n wg is input rotation velocity of harmonic drive, and d m is pitch diameter of the bearing. As the gear meshing is similar to the relative sliding of two flat surfaces, F fm is evaluated through the Coulomb and linear viscous friction model as follows: where F c is Coulomb friction and c is viscous coefficient. us, the whole nonlinear friction torque model of harmonic drive is given by e parameters in F fb could be determined according to the corresponding references or manual [32]. e terms F c and c have to be estimated through experimental data by a mathematical fitting method like the least square method. e simulation results of the nonlinear friction model are also shown in Figure 7.

Kinematic Error.
e manufacturing and assembly errors of components of the harmonic drive induce the kinematic error; it shows the deviation between the actual output rotation angle and the ideal output rotation angle of the harmonic drive as follows: is periodic kinematic error will act as an exciter and produce undesirable vibration effects. At high velocity, these vibrations will become dominant, particularly at the resonant frequency. e kinematic error could inevitably induce output velocity fluctuations and affect precision positioning and tracking operations of harmonic drive servo performance [15]. erefore, this equivalent synthetic eccentric error should be understood and identified for error compensation of harmonic drive servo control system. If the pure kinematic error needs to be measured, the measurement experiment should be carried out at a low driving velocity and without load to avoid inducing the stiffness-related error component [15]. So the driving servo motor was controlled at the constant velocity of 15 rpm, and the input and output angles of the harmonic drive were acquired through encoders. Figure 8 shows the measured kinematic errors for two cycles. Such a harmonic of integer frequencies with respect to the input revolutions could be expressed by Fourier series method as follows: θ ke rec � a 0 2 + k n�1 a n cos ω n θ i + b n sin ω n θ i , where θ ke_rec is reconstructed kinematic error and ω n is angular frequency which could be determined from the FFT of the harmonic. a 0 is DC components; a n and b n are coefficients of the harmonic components that could be estimated through the numerical integration of the experimental data. Figure 8 also shows the simulation results of the Fourier series method of the kinematic error. e proposed harmonic drive model of Figure 3 represents the kinematic error by an equivalent synthetic eccentric error e. According to the geometric relation of the model, the relationship between e and θ ke_rec is as follows:

Simulation and Experimental Verification
According to (3)−(6), when the attributes are identified and the input rotation angle θ i is given, the output rotation θ o could be predicted. Figure 9 shows the simulation block  Shock and Vibration diagram for predicting the harmonic drive dynamic performance. Runge-Kutta method with varying steps is used to solve the mathematical equations, and the simulations of system dynamic performance with different operating conditions or parameters are performed through Matlab/ Simulink commercial software [21]. Table 1 summarizes the parameters used in the model simulation.
As the existence of the dynamic factors such as compliance and friction, the deviation between the actual output rotation angle and the ideal output rotation angle of the harmonic drive under high velocity is not the same as the kinematic error. Here, the deviation is called dynamic transmission error and its calculation method is as (17). Four experiments on the dynamic transmission error were carried out under driving velocities of 60, 600, 1800, and 3000 rpm. e measured input rotation angle from the motor encoder is used as the input rotation angle θ i of the model simulation.
en the corresponding dynamic transmission errors are simulated through the proposed model, and the simulation results are compared with the experimental ones ( Figure 10). As shown in Figure 10, the simulated dynamic transmission errors of the proposed model are generally in good agreement with the experimental results at each velocity. e prediction accuracy of the proposed model is verified; then it could be used to analyze the influence of the components' attributes on the system response.

Dynamic Response Analysis
As the proposed model considers the components' attributes, their influence on the system response could be captured here. Because the position information is not sensitive to the change of the attribute parameters, the velocity step response of harmonic drive is used as the analysis object, and the stiffness of every component is regarded as  Kinematic error Other parameters 25 2.57 1.2 × 10 3 20 2.55    A stiffness ratio factor C K is defined for sensitivity analysis; it means that each component stiffness coefficient will change C K times. Because it is not easy to change the stiffness of the component significantly, the range of the stiffness ratio factor C K is set to 0.8∼1.2. e peak values of fluctuation of the velocity step responses are calculated, and their change rules with C K are shown in Figure 13.
According to the figure, when C K is changed from 0.8 to 1.2, the velocity fluctuation peak value changes little (<15%) at the low velocity of 60 rpm, while it changes significantly (>35%) at 3000 rpm. e velocity fluctuation of harmonic drive is more sensitive to the high driving velocity. As C K increases, the velocity fluctuation of harmonic drive enlarges, which means that the increase of the stiffness of each component will decrease the dynamic transmission accuracy of harmonic drive. At low velocity, the velocity fluctuation is more sensitive to the bearing stiffness K b . Moreover, at high velocity, the influence of each stiffness on the velocity fluctuation is similar when C K is lower than 1, but the fluctuation becomes very sensitive to the bearing stiffness K b when C K is higher than 1. On the whole, among the three stiffness K b , K m , and K t , the change of K b has the most significant impact on system's dynamic response.

Conclusions
To figure out the factors which influence the running performance of the harmonic drive and their influencing rules, a refined harmonic drive model is presented, and the sensitivity analysis of the components' stiffness on the system response is carried out. e following conclusions can be drawn: (1) A refined harmonic drive model considering the radial stiffness of bearing, the meshing stiffness, the flexspline cylinder torsional stiffness, the friction in bearing, the meshing friction, and the kinematic error is established. rough this model, the dynamic transmission error of harmonic drive under different driving velocities could be predicted accurately.
(2) e stiffness, kinematic error, and friction of harmonic drive are explained by power function, Fourier series method, and empirical and classic friction model. A dedicated experimental apparatus based on double motor twisting was constructed. rough the apparatus, the stiffness, kinematic error, and friction of the whole harmonic drive are measured, and the attribute parameters are identified by the corresponding mathematical fitting methods.
(3) Based on the proposed model, the influence of different component stiffness on the velocity step response of harmonic drive is analyzed. e analysis results show that the component stiffness greatly influences the system dynamic response at high driving velocity. Increasing the stiffness of each component will enlarge the dynamic transmission error of harmonic drive; moreover, among all the components' stiffness, the bearing stiffness is the most sensitive dynamic parameter.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.