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

Authors: Xian Zhang; Tao Tao; Gedong Jiang; Xuesong Mei; Chuang Zou
Journal: Shock and Vibration
 (2020)

---

## Abstract

To highlight the key factors which influence the dynamic performance of the harmonic drive, a refined harmonic drive model considering nonlinear stiffness, 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 identified. A series of experiments on the dynamic transmission error at different driving velocities are carried out to verify the proposed model. Based on the proposed model, the influence of different component stiffness on the velocity step response of the harmonic drive is analyzed. The results show that the influence of the component stiffness on the system dynamic response is more significant at high driving velocity, the increase of the stiffness of each component will decrease the dynamic transmission accuracy of harmonic drive, and the bearing radial stiffness is the most sensitive parameter to system’s dynamic response among all the stiffness factors.


---

## Body

## 1. 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. Therefore, 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 Figure1: flexspline, circular spline, and wave generator. The 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. The 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. The 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.Figure 1
Exploded view of harmonic drive and installation.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. Then, 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–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. Thus, 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. Therefore, 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. Through 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. Therefore, 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. These 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. The results of dynamic response analysis will be useful to the dynamic design and manufacturing of high-quality and high-performance harmonic drives.


## 2. Harmonic Drive Modeling

The structural relation and motion form of harmonic drive components are shown in Figure2. The 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. Then 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.Figure 2
Structural relation and motion form of harmonic drive components.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 Figure1, 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]. Thus, 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 Figure3). In Figure 3, the input axis rotates, and the elliptical cam moves. Then the elliptical cam motion forces the bearing to produce radial motion (x1). Through the bearing, x1 makes the flexspline teeth move x2 along the radial direction to engage into the circular spline teeth. The meshing effect produces the tangential motion of flexspline (y1) which, in turn, generates the tangential motion of output axis (y2) through the torsion of the flexspline cylinder. The term e is the equivalent synthetic eccentric error that leads to the pure kinematic error, Kb and Cb are radial stiffness and damping of the bearing, Ffb is friction in the bearing, Kmj is meshing stiffness of the jth gear pair, Ffm is the meshing friction, and Kt is torsional stiffness of the flexspline cylinder. θi, θo, rg, αt, Jg, and Jo 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(1)αn=arctan1Ntanαt,where N is gear ratio of the harmonic drive.Figure 3
The proposed harmonic drive model.(2)Km=∑j=1Kmj.As Figure2 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 Km is used to represent the nonlinear meshing stiffness asThen the relationship between input axis rotation and radial motion of bearing can be formulated as(3)x1=θirgtanαn.According to the Newtonian equations of motion, the ordinary differential equations for the radial motion and force transmission process are(4)Jgrg2x¨1x¨21+−CbCb−Cbe˙x˙1x˙21+−KbKb+Kmsin2αt−Kbex1x21=Kmsinαtcosαt−Ffmcosαty1y21.Meanwhile, the ordinary differential equations for the tangential motion and force transmission process are(5)Jgrg2y¨1y¨21+Ktrg2+Kmcos2αt-Ktrg2y1y21=Kmsinαtcosαt-Ffmsinαt-Ffbx1x21Jorg2y¨1y¨21+-Ktrg2Ktrg2y1y21=0The output rotations of the harmonic drive are(6)θo=y2rg.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.


## 3. Experimental Apparatus and Attributes Identification

### 3.1. Experimental Apparatus

A dedicated experimental apparatus based on double motor twisting (Figures4 and 5) was built. The 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. The input angle and input torque of harmonic drive are measured through the servo driver. All angle and torque signals were acquired synchronously. The loading system was composed of the torque sensor, twisted bar, attendant reducer, and loading servo motor. The 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.Figure 4
Schematic diagram of the experimental apparatus.Figure 5
The experimental apparatus.


### 3.2. Attribute Measurement and Parameters Identification

#### 3.2.1. 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. The stiffness curve of the harmonic drive was obtained, and it is shown in Figure6.Figure 6
Measured stiffness curve and simulation results of stiffness model.The relationship between the output torsional angle and the output torque could be described as follows:(7)Δθo=g1To1/3+g2To,where Δθo is output torsional angle, To is output torque, and g1 and g2 are nonlinear, unknown stiffness coefficients that could be identified through the least square method. The 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 stiffnessKm, flexspline cylinder torsional stiffness Kt, and radial stiffness of bearing Kb. Their corresponding output torsional compliance components of the harmonic drive are fm, ft, and fb, respectively. The output torsional compliance fHD of the whole harmonic drive could be obtained as follows:(8)fHD=fm+ft+fb.Kt could be calculated as a torsion problem of an ordinary cylinder and Kb could be estimated according to the method of stiffness calculation for an ordinary ball bearing, which will not be discussed in this paper. Their corresponding output torsional compliance components ft and fb could be calculated through the structure of the proposed model as follows:(9)ft=1Kt,(10)fb=1Kbrg2cot2αt.Then, according to (7)−(10), the compliance fm induced by gear meshing could be calculated as(11)fmTo=dΔθodTo−ft−fb.At last, the nonlinear meshing stiffnessKm could be determined according to the structure of the proposed model as follows:(12)KmTo=1fmrg2cos2αt.According to (7)−(12), the nonlinear relationship between the meshing stiffness and the output torque Km=fTo is obtained, and then the meshing stiffness Km could be adjusted according to To. In the model, the output torque To could be obtained as(13)To=Kty1−y2rg.


#### 3.2.2. Friction

Friction in the harmonic drive is a velocity-related attribute [29]; it has to be measured at different velocities. Thus, 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.Figure 7
Measured friction torque and simulation results of nonlinear friction model.The harmonic drive contains two main friction sources: the friction in bearingFfb and the meshing friction Ffm. Ffb could be modeled based on the classic bearing Palmgren empirical equation [32]; that is,(14)Ffb=10−7f0υ0nwg2/3dg3sgnnwgrg,where f0 is a coefficient determined by the bearing and lubrication style, υ0 is viscosity of the lubricant, nwg is input rotation velocity of harmonic drive, and dm is pitch diameter of the bearing.As the gear meshing is similar to the relative sliding of two flat surfaces,Ffm is evaluated through the Coulomb and linear viscous friction model as follows:(15)Ffm=Fc+cnwgsgnnwg,where Fc is Coulomb friction and c is viscous coefficient.Thus, the whole nonlinear friction torque model of harmonic drive is given by(16)Tfv=Ffbrg+Ffmrgsinαt.The parameters inFfb could be determined according to the corresponding references or manual [32]. The terms Fc and c have to be estimated through experimental data by a mathematical fitting method like the least square method. The simulation results of the nonlinear friction model are also shown in Figure 7.


#### 3.2.3. Kinematic Error

The 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:(17)θke=θo−θiN.This 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. The kinematic error could inevitably induce output velocity fluctuations and affect precision positioning and tracking operations of harmonic drive servo performance [15]. Therefore, 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:(18)θke_rec=a02+∑n=1kancosωnθi+bnsinωnθi,where θke_rec is reconstructed kinematic error and ωn is angular frequency which could be determined from the FFT of the harmonic. a0 is DC components; an and bn 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.Figure 8
Measured kinematic error and simulation results of the Fourier series method.The proposed harmonic drive model of Figure3 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:(19)e=θke_recrgtanαt.


## 3.1. Experimental Apparatus

A dedicated experimental apparatus based on double motor twisting (Figures4 and 5) was built. The 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. The input angle and input torque of harmonic drive are measured through the servo driver. All angle and torque signals were acquired synchronously. The loading system was composed of the torque sensor, twisted bar, attendant reducer, and loading servo motor. The 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.Figure 4
Schematic diagram of the experimental apparatus.Figure 5
The experimental apparatus.


## 3.2. Attribute Measurement and Parameters Identification

### 3.2.1. 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. The stiffness curve of the harmonic drive was obtained, and it is shown in Figure6.Figure 6
Measured stiffness curve and simulation results of stiffness model.The relationship between the output torsional angle and the output torque could be described as follows:(7)Δθo=g1To1/3+g2To,where Δθo is output torsional angle, To is output torque, and g1 and g2 are nonlinear, unknown stiffness coefficients that could be identified through the least square method. The 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 stiffnessKm, flexspline cylinder torsional stiffness Kt, and radial stiffness of bearing Kb. Their corresponding output torsional compliance components of the harmonic drive are fm, ft, and fb, respectively. The output torsional compliance fHD of the whole harmonic drive could be obtained as follows:(8)fHD=fm+ft+fb.Kt could be calculated as a torsion problem of an ordinary cylinder and Kb could be estimated according to the method of stiffness calculation for an ordinary ball bearing, which will not be discussed in this paper. Their corresponding output torsional compliance components ft and fb could be calculated through the structure of the proposed model as follows:(9)ft=1Kt,(10)fb=1Kbrg2cot2αt.Then, according to (7)−(10), the compliance fm induced by gear meshing could be calculated as(11)fmTo=dΔθodTo−ft−fb.At last, the nonlinear meshing stiffnessKm could be determined according to the structure of the proposed model as follows:(12)KmTo=1fmrg2cos2αt.According to (7)−(12), the nonlinear relationship between the meshing stiffness and the output torque Km=fTo is obtained, and then the meshing stiffness Km could be adjusted according to To. In the model, the output torque To could be obtained as(13)To=Kty1−y2rg.


### 3.2.2. Friction

Friction in the harmonic drive is a velocity-related attribute [29]; it has to be measured at different velocities. Thus, 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.Figure 7
Measured friction torque and simulation results of nonlinear friction model.The harmonic drive contains two main friction sources: the friction in bearingFfb and the meshing friction Ffm. Ffb could be modeled based on the classic bearing Palmgren empirical equation [32]; that is,(14)Ffb=10−7f0υ0nwg2/3dg3sgnnwgrg,where f0 is a coefficient determined by the bearing and lubrication style, υ0 is viscosity of the lubricant, nwg is input rotation velocity of harmonic drive, and dm is pitch diameter of the bearing.As the gear meshing is similar to the relative sliding of two flat surfaces,Ffm is evaluated through the Coulomb and linear viscous friction model as follows:(15)Ffm=Fc+cnwgsgnnwg,where Fc is Coulomb friction and c is viscous coefficient.Thus, the whole nonlinear friction torque model of harmonic drive is given by(16)Tfv=Ffbrg+Ffmrgsinαt.The parameters inFfb could be determined according to the corresponding references or manual [32]. The terms Fc and c have to be estimated through experimental data by a mathematical fitting method like the least square method. The simulation results of the nonlinear friction model are also shown in Figure 7.


### 3.2.3. Kinematic Error

The 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:(17)θke=θo−θiN.This 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. The kinematic error could inevitably induce output velocity fluctuations and affect precision positioning and tracking operations of harmonic drive servo performance [15]. Therefore, 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:(18)θke_rec=a02+∑n=1kancosωnθi+bnsinωnθi,where θke_rec is reconstructed kinematic error and ωn is angular frequency which could be determined from the FFT of the harmonic. a0 is DC components; an and bn 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.Figure 8
Measured kinematic error and simulation results of the Fourier series method.The proposed harmonic drive model of Figure3 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:(19)e=θke_recrgtanαt.


## 3.2.1. 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. The stiffness curve of the harmonic drive was obtained, and it is shown in Figure6.Figure 6
Measured stiffness curve and simulation results of stiffness model.The relationship between the output torsional angle and the output torque could be described as follows:(7)Δθo=g1To1/3+g2To,where Δθo is output torsional angle, To is output torque, and g1 and g2 are nonlinear, unknown stiffness coefficients that could be identified through the least square method. The 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 stiffnessKm, flexspline cylinder torsional stiffness Kt, and radial stiffness of bearing Kb. Their corresponding output torsional compliance components of the harmonic drive are fm, ft, and fb, respectively. The output torsional compliance fHD of the whole harmonic drive could be obtained as follows:(8)fHD=fm+ft+fb.Kt could be calculated as a torsion problem of an ordinary cylinder and Kb could be estimated according to the method of stiffness calculation for an ordinary ball bearing, which will not be discussed in this paper. Their corresponding output torsional compliance components ft and fb could be calculated through the structure of the proposed model as follows:(9)ft=1Kt,(10)fb=1Kbrg2cot2αt.Then, according to (7)−(10), the compliance fm induced by gear meshing could be calculated as(11)fmTo=dΔθodTo−ft−fb.At last, the nonlinear meshing stiffnessKm could be determined according to the structure of the proposed model as follows:(12)KmTo=1fmrg2cos2αt.According to (7)−(12), the nonlinear relationship between the meshing stiffness and the output torque Km=fTo is obtained, and then the meshing stiffness Km could be adjusted according to To. In the model, the output torque To could be obtained as(13)To=Kty1−y2rg.


## 3.2.2. Friction

Friction in the harmonic drive is a velocity-related attribute [29]; it has to be measured at different velocities. Thus, 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.Figure 7
Measured friction torque and simulation results of nonlinear friction model.The harmonic drive contains two main friction sources: the friction in bearingFfb and the meshing friction Ffm. Ffb could be modeled based on the classic bearing Palmgren empirical equation [32]; that is,(14)Ffb=10−7f0υ0nwg2/3dg3sgnnwgrg,where f0 is a coefficient determined by the bearing and lubrication style, υ0 is viscosity of the lubricant, nwg is input rotation velocity of harmonic drive, and dm is pitch diameter of the bearing.As the gear meshing is similar to the relative sliding of two flat surfaces,Ffm is evaluated through the Coulomb and linear viscous friction model as follows:(15)Ffm=Fc+cnwgsgnnwg,where Fc is Coulomb friction and c is viscous coefficient.Thus, the whole nonlinear friction torque model of harmonic drive is given by(16)Tfv=Ffbrg+Ffmrgsinαt.The parameters inFfb could be determined according to the corresponding references or manual [32]. The terms Fc and c have to be estimated through experimental data by a mathematical fitting method like the least square method. The simulation results of the nonlinear friction model are also shown in Figure 7.


## 3.2.3. Kinematic Error

The 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:(17)θke=θo−θiN.This 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. The kinematic error could inevitably induce output velocity fluctuations and affect precision positioning and tracking operations of harmonic drive servo performance [15]. Therefore, 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:(18)θke_rec=a02+∑n=1kancosωnθi+bnsinωnθi,where θke_rec is reconstructed kinematic error and ωn is angular frequency which could be determined from the FFT of the harmonic. a0 is DC components; an and bn 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.Figure 8
Measured kinematic error and simulation results of the Fourier series method.The proposed harmonic drive model of Figure3 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:(19)e=θke_recrgtanαt.


## 4. 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 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.Figure 9
Simulation block diagram of the harmonic drive model.Table 1
Parameters used in the model simulation.
AttributesValuesStiffnessKtNm/radKbN/mKmg1rad/Nm1/3g2rad/Nm3.25×1051×1081.756×10−51.95×10−5Frictionf0υ0mm2/sdmmmFc(N)cN⋅min/r1.4580586.436.72×10−4Kinematic errora0×10−3degω1ω20.461524a1×10−3degb1×10−3dega2×10−3degb2×10−3deg0.20574.50000.0110−1.0000Other parametersrgmmJg×10−7kgm2CbN⋅s/mαtdegJo×10−3kgm235.252.571.2×103202.55As 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. The measured input rotation angle from the motor encoder is used as the input rotation angle θi of the model simulation. Then the corresponding dynamic transmission errors are simulated through the proposed model, and the simulation results are compared with the experimental ones (Figure 10).Figure 10
Comparison between simulations and experiments of dynamic transmission error. (a) 60 rpm. (b) 600 rpm. (c) 1800 rpm. (d) 3000 rpm.
(a)(b)(c)(d)As shown in Figure10, the simulated dynamic transmission errors of the proposed model are generally in good agreement with the experimental results at each velocity. The 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.


## 5. 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 the analysis factor. Figures11 and 12 show the changes of output velocity with the component stiffness. When the stiffness coefficients of different components are changed to 1.2 times separately, the output velocity changes little under low velocity of 60 rpm, while the output velocity changes significantly under high velocity of 3000 rpm.Figure 11
Velocity step responses of different stiffness at 60 rpm.
(a)(b)(c)(d)Figure 12
Velocity step responses of different stiffness at 3000 rpm.
(a)(b)(c)(d)A stiffness ratio factorCK is defined for sensitivity analysis; it means that each component stiffness coefficient will change CK times. Because it is not easy to change the stiffness of the component significantly, the range of the stiffness ratio factor CK is set to 0.8∼1.2. The peak values of fluctuation of the velocity step responses are calculated, and their change rules with CK are shown in Figure 13.Figure 13
Peak values of velocity fluctuation with different stiffness ratio factors. (a) 60 rpm. (b) 3000 rpm.
(a)(b)According to the figure, whenCK 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. The velocity fluctuation of harmonic drive is more sensitive to the high driving velocity. As CK 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 Kb. Moreover, at high velocity, the influence of each stiffness on the velocity fluctuation is similar when CK is lower than 1, but the fluctuation becomes very sensitive to the bearing stiffness Kb when CK is higher than 1. On the whole, among the three stiffness Kb, Km, and Kt, the change of Kb has the most significant impact on system’s dynamic response.


## 6. 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. The 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. Through this model, the dynamic transmission error of harmonic drive under different driving velocities could be predicted accurately.(2)
The 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. Through 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. The 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.



---

*Source: 1841724-2020-06-29.xml*