Rotation-joint stiffness modeling for industrial robots considering contacts

Abstract

Joint flexibility has a major impact on the motion accuracy of a robotic end effector, particularly at high speeds. This work proposes a technique of precisely modeling the torsional stiffness of the rotational joints for the industrial robots. This technique considers the contacts that exist in the joint system, which can have a significant effect on the overall joint stiffness. The torsional stiffness of the connections that commonly exist in the rotational joints, such as the belt connection, the connections using key, bolts, and pins, were modeled by combining the force analysis and the fractal theory. Through modeling the equivalent stiffness for the springs in serial and in parallel, the torsional stiffness of all joints for the ER3A-C60 robot were calculated and analyzed. The results show that the estimated stiffness based on the proposed technique is closer to the actual values than that based on the previous model without considering the contacts. The analysis is useful for controlling the dynamic characteristic of the industrial robots with the rotational joints while planning the trajectory for the end effector.

Introduction

The joint flexibility of a robot is dependent on its transmission system and can greatly influence motion accuracy of the end effector during high-speed motion, which is often used as an index in evaluating performance.^1,2 In the earlier studies, the flexibility of the robotic links has been considered in the robotic motion control by many scholars. Katic and Vukobratovic ³ presented two adaptive control techniques for the robots with flexible links. Chen ⁴ presented a linearized dynamic equation for the multi-link manipulators with considering the elastic deformations of all links. Filipovic et al. ⁵ established the kinematic and dynamic equations through assuming robotic links as elastic and stiff elements. However, since the industrial robotic joint is composed of numerous components, such as the motor used to output torque, and the reducer and synchronous belts as the transmission mechanism, the flexibility of the links is not the only or the most important factor that affects the whole joint stiffness. Therefore, in later studies, two methods were presented to obtain values for the whole joint stiffness: the identification method and theoretical modeling.

Due to its reliability, the identification method has been widely used in robotics research to predict joint stiffness. Abele et al. ⁶ identified the torsional stiffness of a particular joint experimentally by clamping the other joints. Using this method, the positioning error of the end effector with a certain load was estimated after identifying the stiffness of all joints. As such, the number of experiments is equal to the number of joints and therefore, as the load changes, more experiments are needed. By introducing the conservative congruence transformation, Dumas et al. ⁷ reduced the number of required experiments and presented a fast, robust procedure to obtain joint stiffness, which can be applied to any six-revolute serial robot. Based on the deflections measured by a laser tracker sensor and force measurements, Tyapin et al. ⁸ proposed an effective identification approach to obtain joint stiffness. Through identifying the stiffness of flexible joints, Tian et al. ⁹ modeled and analyzed the stiffness of a multiple coordinated robot system. While the identification method is valid for industrial robots, it is limited by the expensive experimental platform and inability to obtain stiffness values at the design stage.

To overcome these limitations, recent studies have proposed a number of approaches based solely on theoretical analysis. The joint stiffness for a six-revolute industrial robot was previously estimated by integrating the stiffness of the gear pairs, reducers, and rotational links in the transmission system and subsequently using a stiffness conversion formula related to the transformation ratios.^10–12 However, the results based on the analysis still produced high prediction errors in comparison to the experimental values for joint stiffness. Studies have shown that 50% of overall stiffness is caused by contacts formed in the connections; ¹³ therefore, contact stiffness in a joint system will have a major influence on the overall dynamic characteristics and is important to assess. ¹⁴ Vukobratovic ¹⁵ provided a review of the control approaches for the robots interacting with the dynamic environment by building a contact between the robotic system and the environment. To improve the accuracy of the predicted values for joint stiffness, the effects of contacts should also be considered in the estimation. Recent studies use the finite element method (FEM) to analyze the stiffness matrix of the revolute joint connection between two adjacent links, where the preload, external wrenches, and friction can be considered as the contact constraints.^16–18 However, FEM is suitable to model the joint stiffness for the redundant legs and parallel kinematic machines (PKMs), which has the simple joint structures limited. As mentioned previously, the industrial robotic joint system is a complex structure composed of several components including the motor, synchronous belt, and reducer, as shown in Figure 1. The use of FEM is limited due to the difficulties of the geometric modeling of the joint structure and the constraint settings. The Hertzian contact stiffness was linearized to build the stiffness for the roller screw mechanism. ¹⁹ This method solves this problem at the macro level, but does not consider the effects of the surface micro-profile parameters on the contact stiffness. Fractal theory is an effective approach to describe the topography of a machined surface and has been widely used in many engineering problems,^20–23 and herein it is introduced into the model of contact stiffness for a number of robotic connections including bolts and pins. In fractal theory, the contact between two rough surfaces is assumed to be a contact between an equivalent fractal rough surface composed of lots of asperities and a rigid surface, and is represented using two fractal parameters (the fractal dimension $\mathrm{D}$ and the fractal roughness parameter $\mathrm{G}$ ). The contact model of a single asperity is shown in Figure 2. In recent years, it has been modified to improve the prediction accuracy of the contact characteristics.^24–26 Moreover, by considering the elastic–plastic deformation at the contact surface, contact normal and tangential stiffness models were established based on three-dimensional (3D) fractal theory and validated by comparing the results to those based on previous methods and also experimentally. ²⁷

Figure 1.

The components of one classic joint system.

Figure 2.

The contact model of a single asperity.

In this article, a technique for modeling the joint stiffness is proposed by considering the influence of contacts within the joint structures, which have not yet been considered in the previous studies.^10–12 In this technique, the normal and tangential stiffness of the connections using keys, bolts, and pins, which commonly exist in the joint system of industrial robots, were modeled based on fractal theory and force analyses. Through building the stiffness transfer equation, the torsional stiffness model for each type of connection was obtained. Then, the total torsional stiffness of most rotation joints in industrial robotic system can be calculated by combining the stiffness of all connections and flexible structures (such as reducers) in parallel or in series. To clearly present this technique, the reminder of this work is organized as follows: section “Fractal contact stiffness” describes the modeling of the fractal contact stiffness, which establishes the mapping relation between the stiffness and the pressure distribution of the surface; section “Torsional stiffness for contacts” gives the torsional stiffness models of kinds of connections based on the fractal model proposed in section “Fractal contact stiffness”; section “Joint stiffness for ER3A-C60 robot” presents the technique of calculating the joint stiffness through combining the stiffness of all connections and the flexible components (such as reducers) in parallel or in series and using ER3A-C60 robot (a classic six-degree-of-freedom industrial robot) as an example; in section “Results and Discussion,” the results based on the presented model were compared with that based on the model without considering the contacts to show the significance of this work, and the effect of the output torque on the joint stiffness was analyzed, which are helpful in controlling the tracking precision of end effector.

Fractal contact stiffness

According to the 3D fractal model presented in the previous work, ²⁷ the normal force would lead to the elastic, elastic–plastic, and fully plastic deformations for the asperities on the rough surface. The largest truncated area ${\mathrm{a}}_{\mathrm{l}}^{\prime }$ of the asperities on the contact surface was used as an intermediate variable to represent the relationships between the total normal load $\mathrm{F}$ and the normal and tangential contact stiffness ${\mathrm{K}}_{\mathrm{n}}\mathrm{and}{\mathrm{K}}_{\mathrm{t}}$ , respectively, which can be expressed as

$$\begin{array}{c}F=\frac{2^{(11-2D)/2}}{3{\pi }^{(4-D)/2}}\cdot \frac{D-1}{(5-2D)}{(ln\gamma )}^{1/2}{G}^{(D-2)}E{\psi }^{(3-D)/2}{a}_{\prime (D-1)/2}^l\left({a^{\prime }}_1^{(5-2D)/2}-{a^{\prime }}_{1c}^{(5-2D)/2}\right)\hfill \\ +\frac{H_{G2}\left(D-1\right){\psi }^{(3-D)/2}{a^{\prime }}_1^{(D-1)/2}}{2(2.26-0.88D)}\cdot \left({a^{\prime }}_{1c}^{(2.26-0.88D)}-{a^{\prime }}_{2c}^{(2.26-0.88D)}\right)\hfill \\ +\frac{H(D-1){a^{\prime }}_1^{(D-1)/2}}{(3-D)}{\psi }^{(3-D)/2}{a^{\prime }}_{2c}^{(3-D)/2}\hfill \end{array}$$

(1)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{K}}_{\mathrm{n}}=\frac{2\sqrt{2}\mathrm{E}(4-\mathrm{D})(\mathrm{D}-1)}{3\sqrt{\pi }(3-\mathrm{D})(2-\mathrm{D})}{\psi }^{(3-\mathrm{D})/2}{{\mathrm{a}}^{\prime }}_1^{(\mathrm{D}-1)/2}}\\ \multicolumn{1}{c}{({{\mathrm{a}}^{\prime }}_1^{(2-\mathrm{D})/2}-{{\mathrm{a}}^{\prime }}_{1\mathrm{c}}^{(2-\mathrm{D})/2})}\\ \multicolumn{1}{c}{+\frac{{\mathrm{H}}_{\mathrm{G}1}{\pi }^{(3-\mathrm{D})/2}(\mathrm{D}-1){\psi }^{(3-\mathrm{D})/2}(1.76-0.38\mathrm{D})}{2^{4-\mathrm{D}}{\mathrm{G}}^{(\mathrm{D}-2)}{(\ln \gamma )}^{1/2}(3-\mathrm{D})(0.76-0.38\mathrm{D})}}\\ \multicolumn{1}{c}{{{{\mathrm{a}}^{\prime }}_{\mathrm{l}}}^{(\mathrm{D}-1)/2}({{\mathrm{a}}^{\prime }}_{1\mathrm{c}}^{(0.76-0.38\mathrm{D})}-{{\mathrm{a}}^{\prime }}_{2\mathrm{c}}^{(0.76-0.38\mathrm{D})})}\end{array}$$

(2)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{K}}_{\mathrm{t}}=\frac{4\mathrm{G}\prime (\mathrm{D}-1)}{\sqrt{2\pi }}{\psi }^{(3-\mathrm{D})/2}{{\mathrm{a}}^{\prime }}_1^{(\mathrm{D}-1)/2}\underset{{{\mathrm{a}}^{\prime }}_{1\mathrm{c}}}{\overset{{{\mathrm{a}}^{\prime }}_{\mathrm{l}}}{\int }}}\\ \multicolumn{1}{c}{{{\mathrm{a}}^{\prime }}^{(-\mathrm{D}/2)}{(1-{\mathrm{H}}_1{{\mathrm{a}}^{\prime }}^{(\mathrm{D}-2)/2})}^{1/3}{\mathrm{da}}^{\prime }}\\ \multicolumn{1}{c}{+\frac{8\mathrm{G}\prime (\mathrm{D}-1)}{\sqrt{\pi }{{\mathrm{H}}_{\mathrm{G}2}}^{1/2}}{\psi }^{(3-\mathrm{D})/2}{{\mathrm{a}}^{\prime }}_1^{(\mathrm{D}-1)/2}\underset{{{\mathrm{a}}^{\prime }}_{2\mathrm{c}}}{\overset{{{\mathrm{a}}^{\prime }}_{1\mathrm{c}}}{\int }}}\\ \multicolumn{1}{c}{{{\mathrm{a}}^{\prime }}^{(-\mathrm{D}/2)}{(1-{\mathrm{H}}_2{{\mathrm{a}}^{\prime }}^{(0.44-0.22\mathrm{D})})}^{1/3}{\mathrm{da}}^{\prime }}\end{array}$$

(3)

where $\mathrm{D}$ and $\mathrm{G}$ are the 3D fractal dimension and fractal roughness parameter, respectively, and are computed using the power spectral density method according to the data of surface profile, which is measured by the topography instrument; $\gamma$ is the dimension parameter of the spectral density which is $\gamma =1.5$ ; ${\mathrm{a}}^{\prime }$ is the truncated area of one asperity on the contact surface; ${{\mathrm{a}}^{\prime }}_{1\mathrm{c}}$ and ${{\mathrm{a}}^{\prime }}_{2\mathrm{c}}$ are the critical truncated areas that distinguish the elastic, elastic–plastic, and plastic regions, respectively; $\psi$ is the expansion coefficient; $\mathrm{E}$ and $G^{\prime }$ are the equivalent elastic and shear moduli of the connectors, respectively; $\mathrm{H}$ denotes the material hardness; and ${\mathrm{H}}_{\mathrm{G}1}$ , ${\mathrm{H}}_{\mathrm{G}2}$ , ${\mathrm{H}}_1$ , and ${\mathrm{H}}_2$ are constants and are determined by $\mathrm{D}$ , $\mathrm{G}$ , $\mathrm{E}$ , and $G^{\prime }$ , respectively. These calculations are not the focus of this work; however, a more detailed explanation of the computing processes can be found in Zhao et al. ²⁷

Torsional stiffness for contacts

Industrial robotic joints are always composed of several components including motors, synchronous belts, gear pairs, and reducers. In this section, the torsional stiffness of the contacts, which commonly exist in the rotational joints, were modeled as follows with the purpose of calculating the contact joint stiffness.

Torsional stiffness for a synchronous belt contact

Since belt materials (e.g. rubber) are less stiffer than pulley materials (e.g. aluminum alloys), the contact stiffness at the contact surface is mostly determined by the softer belt material, which has a normal stiffness ${\mathrm{K}}_{\mathrm{n}}$ and a tangential stiffness ${\mathrm{K}}_{\mathrm{t}}$ of $9\times 10^4$ and $6\times 10^3\mathrm{N}/\mathrm{m}$ , which are equal to the tensile strength and adhesive strength, respectively. ²⁸ The torsional stiffness is established based on both ${\mathrm{K}}_{\mathrm{n}}$ and ${\mathrm{K}}_{\mathrm{t}}$ , and according to the tooth type of the belt. To present the process of obtaining the torsional stiffness of a belt, the S3M-type circular-arc-tooth synchronous belt (MISUMI, Japan) was utilized and is also used in the ER3A-C60 robot. To simplify the modeling, we choose a point as the equivalent contact point, which is located on the circle with the radius as $({\mathrm{R}}_{\mathrm{belt}}-({\mathrm{h}}_{\mathrm{g}}/2))\mathrm{mm}$ . Assume that a normal force and a tangential force, which are equal to the values of the normal and tangential stiffness of the contact, respectively, are applied on the tooth.

As shown in Figure 3, the contact point B has the unit normal and tangential deformations as ${\epsilon }_{\mathrm{n}}={\epsilon }_{\mathrm{t}}=1\mathrm{m}$ . In $\Delta \mathrm{AOB}$ , the angle $\angle \mathrm{ABO}$ can be determined using $\angle \mathrm{ABO}=$ $\angle \mathrm{ABC}+\angle \mathrm{CBO}$ , where $\angle \mathrm{CBO}=\alpha -\beta$ . Since ${\epsilon }_n={\epsilon }_t=1\mathrm{m}$ , we can obtain $\angle \mathrm{ABC}={45}^{\circ }$ and $\mathrm{AB}=\sqrt{2}m$ . Finally, $\angle \mathrm{AB}O=45°+\alpha -\beta$ ${45}^{\circ }+\alpha -\beta$ . Then $\angle \mathrm{BAO}={180}^{\circ }-\angle \mathrm{AOB}-\angle$ $\mathrm{ABO}={135}^{\circ }-\delta -\alpha +\beta$ , where $\delta$ denotes the rotational deformation. According to Sine theorem, we obtain $\mathrm{OB}/\sin (\angle \mathrm{BAO})=\mathrm{AB}/$ $\sin (\angle \mathrm{AOB})$ , which can be rewritten as follows

$$\frac{({\mathrm{R}}_{\mathrm{belt}}-\left(\frac{{\mathrm{h}}_{\mathrm{g}}}{2}\right))}{\sin ({135}^{\circ }-\delta -\alpha +\beta )}=\frac{\sqrt{2}}{\sin \delta }$$

(4)

where ${\mathrm{R}}_{\mathrm{belt}}$ denotes the addendum-circle radius for the pulley, $\alpha$ and $\beta$ are the angles related to the equivalent contact point shown in Figure 3, and ${\mathrm{h}}_{\mathrm{g}}$ denotes the tooth height of the pulley.

Figure 3.

The force analysis of a single belt tooth and the deformation transition sketch.

The torque $\tau$ that causes unit deformations in the normal and tangential directions can be expressed in terms of ${\mathrm{K}}_{\mathrm{n}}$ and ${\mathrm{K}}_{\mathrm{t}}$ as

$$\tau =({\mathrm{R}}_{\mathrm{belt}}-\left(\frac{{\mathrm{h}}_{\mathrm{g}}}{2}\right))\times [{\mathrm{K}}_{\mathrm{n}}\cos (\alpha -\beta )+{\mathrm{K}}_{\mathrm{t}}\sin (\alpha -\beta )]$$

(5)

For a single belt, the torsional stiffness ${\mathrm{K}}_{\mathrm{belt}}$ can be deduced according to the definition of stiffness and is written as

$${\mathrm{K}}_{\mathrm{belt}}=\frac{{\mathrm{n}}_{\mathrm{belt}}\tau }{\delta }=\frac{{\mathrm{n}}_{\mathrm{belt}}({\mathrm{R}}_{\mathrm{belt}}-\left(\frac{{\mathrm{h}}_{\mathrm{g}}}{2}\right))\times [{\mathrm{K}}_{\mathrm{n}}\cos (\alpha -\beta )+{\mathrm{K}}_{\mathrm{t}}\sin (\alpha -\beta )]}{\arctan \left(\frac{\sqrt{2}\sin ({135}^{\circ }-\alpha -\beta )}{({\mathrm{R}}_{\mathrm{belt}}-\left(\frac{{\mathrm{h}}_{\mathrm{g}}}{2}\right))+\sqrt{2}\cos ({135}^{\circ }-\alpha -\beta )}\right)}$$

(6)

where ${\mathrm{n}}_{\mathrm{belt}}$ denotes the number of the engaged teeth.

For other types of synchronous belt contacts, the torsional stiffness can be obtained using the same modeling process whereby the rotational deformation $\delta$ and the corresponding torque are determined according to the geometrical parameters.

The process of calculating the torsional stiffness of the synchronous belt contact consists of four steps: first, applying the normal and tangential forces, which are, respectively, equal to the normal and tangential stiffness of the contact; second, according to the geometric relationships, the torsional deformations can be obtained; third, the torque can be expressed in terms of the normal and tangential stiffness; finally, the torsional stiffness, which is defined as the ratio of the torque and the torsional deformation, can be calculated. Although the gear pair contact and the synchronous belt contact are both the contacts between two teeth, for the former one the materials of the contact teeth have a similar strength. Unlike the latter one, its normal and tangential stiffness can be obtained according to the fractal model, which has been presented in section “Fractal contact stiffness.” Then following the four-step process, the torsional stiffness of the gear pair contact can be computed.

Torsional stiffness for a key connection

A key is an important connector used to transfer torque between two structures. The sketch of a key connection is shown in Figure 4, where $\mathrm{L}$ , $\mathrm{h}$ , and $\mathrm{w}$ are the length, height, and width of the key, respectively. ${\mathrm{R}}_{\mathrm{k}}$ is the shaft radius of the key connection structure and ${\mathrm{h}}_0$ denotes the height of the upper contact surface. While transferring the torque, the force-bearing contact has a significant impact on torsional stiffness, which thus determines the transmission accuracy of the assembly.

Figure 4.

Sketch of a key connection.

According to Figure 4, the equilibrium equations for the key connection system can be written as

$$\{\begin{array}{c}{\mathrm{F}}_2-{\mathrm{F}}_1={\mathrm{m}}_{\mathrm{k}}({\mathrm{R}}_{\mathrm{k}}^2-\frac{{\mathrm{w}}^2}{4}){\mathrm{a}}^2\hfill \\ \multicolumn{1}{c}{{\mathrm{F}}_2(\sqrt{{\mathrm{R}}_{\mathrm{k}}^2-\frac{{\mathrm{w}}^2}{4}}-\frac{\mathrm{h}-{\mathrm{h}}_0}{2})-{\mathrm{F}}_1(\sqrt{{\mathrm{R}}_{\mathrm{k}}^2-\frac{{\mathrm{w}}^2}{4}}+\frac{{\mathrm{h}}_0}{2})=\tau }\end{array}$$

(7)

where ${\mathrm{F}}_1$ and ${\mathrm{F}}_2$ denote the total normal loads on the contact surfaces; ${\mathrm{m}}_{\mathrm{k}}$ and $\mathrm{a}$ are the mass angular acceleration values of the key.

Substituting force ${\mathrm{F}}_1$ or ${\mathrm{F}}_2$ into the fractal stiffness model, the normal and tangential stiffness ${\mathrm{K}}_{\mathrm{n}}$ and ${\mathrm{K}}_{\mathrm{t}}$ of the contact surfaces can be calculated. To further obtain the torsional stiffness of the contact surface, the transfer equation needs to be modeled. Assuming that the points with the same radius at one cross-section of the structure have the same displacement along the circumferential direction, then the torsional stiffness can be defined as the ratio of the torque τ applied on this circumference and the angular displacement δ of the point, which is shown in Figure 5.

Figure 5.

Sketch of the torsional stiffness computation process.

This figure indicates that the force along the direction of the rotational axis has little effect on the angular displacement. When applying a circumferential force, which is equal to, in number, the stiffness ${\mathrm{K}}_{\mathrm{c}}$ in that direction, the torque and the angular displacement of the point on the contact surface can be expressed as follows

$$\tau ={\mathrm{K}}_{\mathrm{c}}\mathrm{R}$$

(8)

$$\delta =\frac{1}{\mathrm{R}}$$

(9)

where $R$ denotes the distance from the point to the rotational axis. Then the torsional stiffness is modeled as follows

$$\mathrm{K}=\frac{\tau }{\delta }={\mathrm{K}}_{\mathrm{c}}{\mathrm{R}}^2$$

(10)

As the normal direction of the key connection contact surface is perpendicular to the cross-section of the rotational axis, the stiffness ${\mathrm{K}}_{\mathrm{c}}$ in the torsional stiffness model is determined by ${\mathrm{K}}_{\mathrm{c}}={\mathrm{K}}_{\mathrm{n}}$ . Taking the center position of the contact surface as the equivalent point applied a normal force ${\mathrm{K}}_{\mathrm{n}}$ , the torsional stiffness of the two contact surfaces are ${\mathrm{K}}_{\mathrm{k}1}={\mathrm{K}}_{\mathrm{n}}{\mathrm{R}}_1^2$ and ${\mathrm{K}}_{\mathrm{k}2}={\mathrm{K}}_{\mathrm{n}}{\mathrm{R}}_2^2$ , where ${\mathrm{R}}_1=\sqrt{{\mathrm{R}}_{\mathrm{k}}^2-({\mathrm{w}}^2/4)}+({\mathrm{h}}_0/2)$ and ${\mathrm{R}}_2=\sqrt{{\mathrm{R}}_{\mathrm{k}}^2-({\mathrm{w}}^2/4)}-((\mathrm{h}-{\mathrm{h}}_0)/2)$ .

To obtain the total rotational stiffness of the contact surface for the key connection, the total rotational stiffness of two springs in series and in parallel were analyzed, as shown in Figure 6.

Figure 6.

Sketch of springs in series and in parallel.

Springs 1 and 2 are in series, while springs $1\prime$ and $2\prime$ are in parallel with the transmission ratios of ${{\lambda }^{\prime }}_1\mathrm{and}{{\lambda }^{\prime }}_2$ from the present spring to the output point $B^{\prime }$ , respectively. Given an input torque $\tau$ $({\tau }^{\prime })$ at point $\mathrm{A}$ $(\mathrm{A}\prime )$ , the rotational displacements ${\delta }_1,{\delta }_2$ $({{\delta }^{\prime }}_1,{{\delta }^{\prime }}_2)$ and torques ${\tau }_1,{\tau }_2$ $({{\tau }^{\prime }}_1,{{\tau }^{\prime }}_2)$ for the serial (or parallel) springs have the following relationships

$$\{\begin{array}{c}\begin{array}{c}{\delta }_1=\frac{{\tau }_1}{{\mathrm{K}}_1}\hfill \\ \multicolumn{1}{c}{{\delta }_2=\frac{{\tau }_2}{{\mathrm{K}}_2}}\\ \multicolumn{1}{c}{{\delta }_1+{\delta }_2=\delta }\\ \multicolumn{1}{c}{{\tau }_1={\tau }_2=\tau }\end{array}\},\mathrm{springs}\mathrm{in}\mathrm{serial}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}{{\delta }^{\prime }}_1=\frac{{{\tau }^{\prime }}_1}{{{\mathrm{K}}^{\prime }}_1}\hfill \\ \multicolumn{1}{c}{{{\delta }^{\prime }}_2=\frac{{{\tau }^{\prime }}_2}{{{\mathrm{K}}^{\prime }}_2}}\\ \multicolumn{1}{c}{{{\tau }^{\prime }}_1={\tau }^{\prime }}\\ \multicolumn{1}{c}{{{\tau }^{\prime }}_2=\frac{{{\lambda }^{\prime }}_1{\tau }^{\prime }}{{{\lambda }^{\prime }}_2}}\end{array}\},\mathrm{springs}\mathrm{in}\mathrm{parallel}}\end{array}$$

(11)

From this, the total rotational stiffness at point $\mathrm{B}$ ( $\mathrm{B}\prime$ ) can be defined as

$$\{\begin{array}{c}\mathrm{K}=\frac{\tau }{\delta }=\frac{\tau }{{\delta }_1+{\delta }_2}=\frac{1}{\left(\frac{1}{{\mathrm{K}}_1}\right)+\left(\frac{1}{{\mathrm{K}}_2}\right)},\mathrm{springs}\mathrm{in}\mathrm{serial}\hfill \\ \multicolumn{1}{c}{{\mathrm{K}}^{\prime }=\frac{{\tau }^{\prime }}{{\delta }^{\prime }}=\frac{{{\lambda }^{\prime }}_1{\tau }^{\prime }}{\left(\frac{{{\delta }^{\prime }}_1}{{{\lambda }^{\prime }}_1}\right)+\left(\frac{{{\delta }^{\prime }}_2}{{{\lambda }^{\prime }}_2}\right)}=\frac{1}{\left(\frac{1}{{{\mathrm{K}}^{\prime }}_1{{\lambda }^{\prime }}_1^2}\right)+\left(\frac{1}{{{\mathrm{K}}^{\prime }}_2{{\lambda }^{\prime }}_2^2}\right)},\mathrm{springs}\mathrm{in}\mathrm{parallel}}\end{array}$$

(12)

By extending this to n springs in series or parallel, which are general to other kinds of connections, the total stiffness can be defined as

$$\{\begin{array}{c}\mathrm{K}=\frac{1}{\frac{1}{{\mathrm{K}}_1}+\cdots +\frac{1}{{\mathrm{K}}_{\mathrm{n}}}},\mathrm{springs}\mathrm{in}\mathrm{serial}\hfill \\ \multicolumn{1}{c}{{\mathrm{K}}^{\prime }=\frac{1}{\frac{1}{{{\mathrm{K}}^{\prime }}_1{{\lambda }^{\prime }}_1^2}+\cdots +\frac{1}{{{\mathrm{K}}^{\prime }}_{\mathrm{n}}{{\lambda }^{\prime }}_{\mathrm{n}}^2}},\mathrm{springs}\mathrm{in}\mathrm{parallel}}\end{array}$$

(13)

Since the torsional stiffness of the two contact surfaces for the key connection (Figure 4) are in parallel, the total torsional stiffness can be modeled as

$${\mathrm{K}}_{\mathrm{k}}=\frac{1}{\left(\frac{1}{{\mathrm{K}}_{\mathrm{k}1}}\right)+\left(\frac{1}{{\mathrm{K}}_{\mathrm{k}2}}\right)}$$

(14)

Torsional stiffness for a bolt connection

In a bolt connection, the pre-tightening forces of the bolts are applied to the contact surface connecting the two parts, and their values greatly affect the contact quality. Modeling the relationship between the pre-tightening forces and rotational stiffness is important in guiding the bolt tightening of each joint.

With the transmission system of an industrial robotic joint, a bolt connection is always adopted to connect the reducer with the robotic body. For most reducers, the contact surface of the bolt connection is a looped region as shown in Figure 7, where ${\mathrm{n}}_{\mathrm{b}}$ is the number of bolts, ${\mathrm{R}}_{\mathrm{b}}$ denotes the nominal radius of each bolt, and ${\mathrm{R}}_{\mathrm{in}}\mathrm{and}{\mathrm{R}}_{\mathrm{ou}}$ represent the inner and outer radius of the circular contact surface, respectively.

Figure 7.

Sketch of a bolt connection.

For this contact surface, the total normal force is determined by the pre-tightening force ${\mathrm{F}}_{\mathrm{b}}$ of each bolt, which can be calculated as

$$\mathrm{F}={\mathrm{n}}_{\mathrm{b}}{\mathrm{F}}_{\mathrm{b}}$$

(15)

Substituting the normal force $F$ into the fractal model, the normal and tangential stiffness of the contact surface can be calculated. Since the normal direction of the bolt connection contact surface is parallel to the rotational axis, the stiffness $K_c$ in equation (10) is determined by ${\mathrm{K}}_{\mathrm{c}}={\mathrm{K}}_{\mathrm{t}}$ . Then the torsional stiffness of this contact surface is modeled as follows

$${\mathrm{K}}_{\mathrm{bo}}={\mathrm{K}}_{\mathrm{t}}{\mathrm{R}}_{\mathrm{bo}}^2$$

(16)

where ${\mathrm{R}}_{\mathrm{bo}}$ denotes the shaft radius of the bolt connection structure, which can be determined by ${\mathrm{R}}_{\mathrm{bo}}=({\mathrm{R}}_{\mathrm{in}}+{\mathrm{R}}_{\mathrm{ou}})/2$ .

Torsional stiffness for a pin connection

In some robotic joints, pins are used to connect two components. A sketch of this type of connection is shown in Figure 8, where ${\mathrm{n}}_{\mathrm{p}}$ denotes the number of pins and ${\mathrm{F}}_{\mathrm{p}}$ is the pre-tightening force of the pin.

Figure 8.

Sketch of a pin connection.

With this connection, several pins with external threads are used to apply pre-tightening forces on the cylindrical surface of the shaft leading to contacts. In the actual assembly work, the preloads of all pins in one pin connection are set to be the same to avoid the instability led by the unevenly distributed force on the contact surface. For each contact surface, the normal and tangential stiffness can be calculated by substituting the normal force ${\mathrm{F}}_{\mathrm{p}}$ into the fractal theory. As shown in Figure 8, as the normal direction at the center of each contact surface and the rotational axis intersect at one point, the stiffness ${\mathrm{K}}_{\mathrm{c}}$ in equation (10) is determined by ${\mathrm{K}}_{\mathrm{c}}={\mathrm{K}}_{\mathrm{t}}$ . Then the torsional stiffness of each contact surface is expressed as $\mathrm{K}={\mathrm{K}}_{\mathrm{t}}{\mathrm{R}}_{\mathrm{pi}}^2$ , where ${\mathrm{R}}_{\mathrm{pi}}$ is equal to the shaft radius of the pin connection radius. Since the stiffness of the pins are in series, the total stiffness of the pin connection can be written as

$${\mathrm{K}}_{\mathrm{pi}}=\frac{\mathrm{K}}{{\mathrm{n}}_{\mathrm{p}}}$$

(17)

Joint stiffness for ER3A-C60 robot

In the simulation analysis, the ER3A-C60 robot is used as the research object. For its first two joints, the harmonic reducer matches with the motor in structure, and they are assembled with a single key. In joints 3, 4, and 5, the synchronous pulley also matches with the motor; therefore, only one key connection exists between the connected components. However, in the sixth joint, a connecting shaft is used due to the mismatch between the reducer and the motor, such that the local assembly consists of two key connections. In addition, a shaft connects the synchronous pulley and the harmonic reducer so there are two connections between them. The harmonic reducers are connected to the robotic body using the bolts. A sketch of the robotic model, where the reducers and the contacts are represented as the torsional springs, is presented in Figure 9.

Figure 9.

Sketch of the model showing each joint and its connection types.

The problem of calculating torsional stiffness for each joint is broken down into several distinct problems of computing the torsional stiffness for each contact and reducer. The torsional stiffness of each contact can be calculated based on the model presented in section “.” The reducer is a complex structure and has been standardized by some manufacturing enterprises, such as Leaderdrive, an enterprise of China. The rotational stiffness would be calculated according to the stiffness characteristic map obtained in the experiment and given in the product manual, since the harmonic reducers of the ER3A-C60 robot are made in the Leaderdrive and meet the national standard. From the PRC national standard GB/T 30819-2014, the information and stiffness values of the harmonic drivers are listed in Table 1, where ${\mathrm{K}}_{\mathrm{r}1}$ , ${\mathrm{K}}_{\mathrm{r}2}$ , and ${\mathrm{K}}_{\mathrm{r}3}$ represent the torsional stiffness when the output torque ranges from 0 to ${\tau }_{\mathrm{r}1}$ , ${\tau }_{r1}$ to ${\tau }_{r2}$ , and greater than ${\tau }_{r2}$ , respectively.

Table 1.

The product models of harmonic reducers for all joints and their critical parameters.

Harmonic reducer	Reducer 12	Reducer 22	Reducer 36	Reducer 46	Reducer 56	Reducer 63
Model	LHS-25-120-C-I-S5	LHS-25-120C-I-S5	LHS-17-100C-I-S5	LHS-17-50C-III	LHS-14-100C-I-S5	LCD-14-50-C-I-S5
Transmissionratio	120:1	120:1	100:1	50:1	100:1	50:1
τr1 (N m)	14	14	3.9	3.9	2.0	2.0
τr2 (N m)	48	48	12	12	6.9	6.9
Kr1 (N m/rad)	3.1 × 104	3.1 × 104	0.81 × 104	0.81 × 104	0.47 × 104	0.34 × 104
Kr2 (N m/rad)	5.0 × 104	5.0 × 104	1.1 × 104	1.1 × 104	0.61 × 104	0.47 × 104
Kr3 (N m/rad)	5.7 × 104	5.7 × 104	1.3 × 104	1.3 × 104	0.71 × 104	0.57 × 104

Then, according to the output torque of each joint, a value of torsional stiffness for the corresponding reducer can be determined. The process of obtaining the joint stiffness is shown in Figure 10 and the parameters used are listed in Tables 1–6.

Figure 10.

Flowchart showing the process of obtaining joint stiffness.

Table 2.

The values of the main parameters for belt connections.

Belt connection		α $(\times \pi /180\mathrm{rad})$	Rbelt $(\times {10}^{-3}\mathrm{m})$	hg $(\times {10}^{-3}\mathrm{m})$	β $(\times \pi /180\mathrm{rad})$	nbelt	Transmission ratio
For joint 3	Belt con. 32	44	13.945	1.11	1.6768	15	1:1
	Belt con. 33		13.945		1.6768	15
For joint 4	Belt con. 42		12.035		1.9492	17	20:13
	Belt con. 43		18.72		1.2426	26
For joint 5	Belt con. 52		11.08		2.1216	12	1:1
	Belt con. 53		11.08		2.1216	12

Table 3.

The processing information of the contacts and the parameters used in fractal theory.

Parameters	Key connection (belt–key contact)	Key and pin connections (key(pins)–reducer contact)	Bolt connection
Processing mode	Grinding	Grinding	Grinding
Ra/μm	0.8	0.8	1.6
Oiled or not	No	No	No
$D$	2.3390	2.3390	2.3075
$G$ (m)	6.4960 × 10−11	6.4960 × 10−11	2.6927 × 10−11
$E$ (Pa)	5.5809 × 1010	1.1538 × 1011	1.1538 × 1011
$G\prime$ (Pa)	1.1442 × 1010	1.3053 × 1010	1.3053 × 1010
${\mathrm{a}\prime }_{1\mathrm{c}}$ (m2)	2.0554 × 10−8	3.1453 × 10−7	1.4607 × 10−6
${\mathrm{a}\prime }_{2\mathrm{c}}$ (m2)	5.7287 × 10−14	8.7661 × 10−13	1.0982 × 10−12
${\mathrm{H}}_{\mathrm{G}1}$	2.9610 × 106	5.4787 × 107	7.8422 × 107
${\mathrm{H}}_{\mathrm{G}2}$	73.2126	42.0346	23.8599
${\mathrm{H}}_1$	45.4996	9.9565	6.2109
${\mathrm{H}}_2$	0.6050	0.2576	0.3171

Table 4.

The values of parameters for key connections.

Key connection	h0 $(\times {10}^{-3}\mathrm{m})$	L $(\times {10}^{-3}\mathrm{m})$	h $(\times {10}^{-3}\mathrm{m})$	w $(\times {10}^{-3}\mathrm{m})$	Rk $(\times {10}^{-3}\mathrm{m})$
Key con. 11	2.3	16.6	5	5	7
Key con. 21	2.3	16.6	5	5	7
Key con. 31	1.2	14	3	3	4
Key con. 34	1.4	11	3	3	4
Key con. 35	1.4	12.8	3	3	4
Key con. 41	1.2	12	3	3	4
Key con. 44	2.8	8	6	6	9
Key con. 51	1.2	12	3	3	4
Key con. 54	1.4	11	3	3	4
Key con. 55	1	10.4	2	2	4
Key con. 61	1	6	2	2	4
Key con. 62	1	4.9	2	2	3

Table 5.

The values of parameters for bolt connections.

Bolt connection	nb	Rb $(\times {10}^{-3}\mathrm{m})$	Rin $(\times {10}^{-3}\mathrm{m})$	Rou $(\times {10}^{-3}\mathrm{m})$	Fb (N)
Bolts con. 13	16	4	33.5	42.5	4800
Bolts con. 23	16	4	33.5	42.5	4800
Bolts con. 37	16	3	30	40	2700
Bolts con. 47	16	3	22.5	32	2700
Bolts con. 57	8	3	25	35	2300
Bolts con. 64	4	5	10	20	6500

Table 6.

The values of parameters for pin connection.

Pin connection	np	Rp $(\times {10}^{-3}\mathrm{m})$	Rpi $(\times {10}^{-3}\mathrm{m})$	Fp (N)
Pins con. 45	3	3	20	2800

Based on that, the joint stiffness of the ER3A-C60 robot can be calculated when the output torque of each joint is known. To verify the effectiveness of the presented model, the results are compared with those based on the model without considering the contacts.^6–8 Besides that, to reveal the influence of the output torque on the joint stiffness, the ratio $\lambda$ of output torques of joints, which is defined as the ratio of the present torque to the allowable maximum torque as listed in Table 7, is taken as the variable to obtain the corresponding values for the torsional stiffness of each joint.

Table 7.

The maximum values of the output torque for each joint.

Joint	1	2	3	4	5	6
Maximum output torque (N m)	103	103	37	25	37	4.6

Results and discussion

The previous studies,^10–12 which proposed a technique of calculating joint stiffness without considering the contacts, show that the predicted joint stiffness is always larger than the experimental one due to the ignorance of the effects of some detailed structure, like contacts. To show the significance of the presented model, we define ${\eta }_1$ as the ratio of the experimental stiffness to the estimated joint stiffness using the previous model and ${\eta }_2$ as the ratio of the estimated joint stiffness using the presented model to that using the previous model. In the study by Qin, ¹¹ the ratios ${\eta }_1$ for the six joints of an industrial robot are 0.534, 0.313, 0.526, 0.212, 0.612, and 0.0724 and the ratios ${\eta }_2$ for the six joints of the ER3A-C60 robot are 0.845, 0.845, 0.579, 0.509, 0.685, and 0.222, respectively. The comparison of the ratios ${\eta }_1$ and ${\eta }_2$ is shown in Figure 11, where the ratio ${\eta }_3$ represents the level of the estimated values using the previous model and is always equal to 1. The values of the ratios ${\eta }_1$ and ${\eta }_2$ state that both the estimated stiffness for each joint using both the experiment and the presented model are much larger than those based on the previous model, that is, the stiffness of all joints based on the presented model are much closer to the actual value. Therefore, this comparison can prove that the presented model is more accurate than the previous one and is reliable in the calculation of joint stiffness.

Figure 11.

Comparison of the ratios ${\eta }_1$ and ${\eta }_2$ of the joint stiffness.

Figure 12 shows the variation of the joint stiffness for each joint with the increasing ratio of the output torques and the comparison of the results based on the previous and the presented models. The results indicate that the joint stiffness based on the presented model, which considers the contacts, is smaller than that based on the previous model, and as the joint output torque becomes closer to zero, the difference between the joint stiffness based on the two methods is larger. Moreover, the stiffness of all joints based on the presented model change nonlinearly due to the nonlinear relationship between contact stiffness and normal force, which is calculated using the torque value of each joint. The sudden changes in the joint stiffness curves are caused by variation in the reducer stiffness values as the output torque of the reducer is within different ranges. The figure also shows that the joint stiffness increases dramatically when the torque ratio is close to zero and tends to vary gently with increasing torque ratio, which implies that the too small output torques would lead to the large dynamic error of the end effector due to the high joint flexibility.

Figure 12.

Variation of joint stiffness with output torque for each joint.

Since the input torques are calculated based on the kinematic parameters (including the angular displacement, velocity, and acceleration) of all joints using the dynamic equations, ²⁹ the joint stiffness would be affected largely by these parameters. This effect would be used to improve the motion accuracy by maximizing the joint stiffness through the trajectory optimization. Furthermore, the surface pressures for the bolt and pin connections are determined by the pre-tightening forces of the bolts and pins. According to Zhao et al., ³⁰ the contact torsional stiffness and the surface pressure have the positive correlation, and the stiffness would be improved very slowly when the pressure increases to a certain level. This relation means that although the increase of the pre-tightening force can improve the joint stiffness, the pre-tightening forces should be optimized according to the actual requirements of the motion accuracy and the properties of the bolts or pins, which can guide the robotic assembly task.

Conclusion

In this work, a technique of modeling the joint stiffness for industrial robots with rotational joints is presented by considering the influence of the contacts. In this technique, the torsional stiffness of the contacts for the connections, which commonly exist in the rotational joint systems, were modeled. Through modeling the equivalent stiffness, the stiffness of the joints for an industrial robot, ER3A-C60, were calculated and analyzed. The results indicate that the estimated stiffness of all joints based on the previous model without considering the contacts are larger than those based on the presented model, which are much closer to the actual values. The comparison can prove that the presented model is more accurate than the previous one and is reliable in the calculation of joint stiffness. For the presented results, the joint stiffness of all joints change nonlinearly with the increasing torque ratio. Moreover, as the torque is closer to zero, the decrease of the joint stiffness is more dramatic. The analysis provides a foundation for improving the dynamic characteristics of the rotational joints while planning trajectory for the end effector.

In the future work, the stiffness of all joints in multiple robotic systems will be obtained using the experimental identification method and compared with that calculated using the presented model to further show the reliability of this model directly. Moreover, this model will be used to calculate the synthesis error of the robotic end effector that is caused by the joint flexibility and geometric errors by combining the dynamic analysis of the industrial robots, and then show its correctness by comparing the end effector motion errors obtained using laser tracker and the presented stiffness model.