A quick stress calculation method for flexspline in harmonic actuators based on the finite element method

Abstract

The harmonic drive (HD) is increasingly used due to its outstanding advantages, especially for equipment requiring high precision and a large gear ratio. This research is based on the Finite Element Method (FEM) and mechanical analysis to investigate stress calculation methods for the flexspline (FS), a significant component of HD transmission, providing a basis for assessing the durability and failure of the FS. First, the FS deformation and stress analysis model based on the FEM is proposed. Since the two boundary elements of the model have the form of a rotation and translational joint, which have not been provided, based on the stiffness matrix of the fixed beam to build the stiffness matrix for these elements, which completes the proposed analytical model. Then, the structures of an HD reduced are designed and applying the proposed force analysis model, the deformation and stress distributions of the FS are obtained. Finally, the theoretical safety factor is calculated and analyzed. The results of the simulation completely coincide with the theoretical calculation and show that the most concentrated stress is at the root of the teeth, which is the leading cause of the FS failure. So, the study as reference material for rapidly evaluating the stress state and calculating the safety when designing and optimizing the FS.

Keywords:

• Harmonic drive
• Flexspline
• Mechanics analysis
• Finite Element Method
• Stress calculation

PUBLIC INTEREST STATEMENT

This research is based on the Finite Element Method (FEM) and mechanical analysis to investigate stress calculation methods for the flexspline (FS), a significant component of HD transmission, providing a valuable method for assessing the durability and failure of the FS. First, the FS deformation and stress analysis model based on the FEM is proposed. Then, the structures of an HD reduced are designed, and applying the proposed force analysis model, the deformation and stress distributions of the FS are obtained. Finally, the theoretical safety factor is calculated and analyzed. The simulation results coincide with the theoretical calculation and show that the most concentrated stress is at the root of the teeth, which is the leading cause of FS failure. So, the study, such as reference material for rapidly evaluating the stress state and calculating the safety when designing and optimizing the FS.

1. Introduction

Harmonic drives (HD) have the advantages of being lightweight and having a high transmission ratio and are widely used in precision transmission fields such as medical equipment, aerospace, and robotics (Brassitos & Jalili, 2017; Shen & Ye, 1985; Shuyan Wang et al., 2021; Ueura & Slatter, 1999). HD transmission uses elastic deformation formed in a flexible part called the flexspline (FS) to perform motion or power transmission (Yu et al., 2021). As an essential component of HD transmission, the FS is an elastic thin-walled structure and operates under alternating loads. As a result, fatigue failure of FS is one of the most common types of failure (Ma et al., 2017; Shi et al., 2010; Wang & Chen, 2005; Y.F. Wang et al., 2005). Therefore, studying the deformation and stress distribution of the FS is the basis for evaluating the FS’s safety, durability, and fatigue failure during the HD gear’s dynamic working process (Sahoo et al., 2020).

The experimental method is commonly used for stress testing of HD (Li, 2016; Pacana et al., 2017; Tjahjowidodo et al., 2013). Sticking strain gauges on the FS’s surface can measure stress (Guan et al., 2022; İrsel, 2021). However, strain gauges cannot obtain stress on teeth because of their narrow space with this approach. Furthermore, the experimental method is hard because of the extended test cycle and high cost (J. Wang et al., 2021). In addition, the deformation and stress distribution of FS is also determined by the empirical formula. M.H. Ivanov (Ivanov, 1987) and Shen et al. (Shen & Ye, 1985) provide a computational model using the column shell theory based on elastic theory and some additional assumptions for the model based on experimental results. In recent years, many works have been based on FEM to obtain a more accurate deformation and stress distribution of the FS (Cheng & Chen, 2022). Kayabasiet and F. Erzincan conducted a stress-strain analysis of the FS based on FEM to calculate the stress on the teeth of the FS and found the optimal tooth geometry to maximize the duration of the FS’s life (Kayabasi & Erzincan, 2005). Through finite element contact analysis, S.M. Zhang determines the boundary conditions and the initial load by surface contact between FS and the wave generator (WG; Zhang, 2009). The deformation and stress can be collected under different conditions, providing valuable data for optimizing the structures of the FS. Chen et al. based on the shell element and distribution of deformation in the transmission state obtained by nonlinear contact analysis set up an FS finite element model (Chen et al., 2011). Zou et al. established a surface contact finite element model of an HD through which they found that strain and stress in the cross-section of the FS gear changed geometry with high loads, even though the distribution of the increments in deformation and stress did not change (Zou et al., 2013). However, the above analysis cannot meet the structure’s integrity requirements because of the high reliability, safety, and long life.

Furthermore, technically FEM depends on computers and commercial software; modeling and preprocessing FEM is time-consuming to process and obtain results (S. Wang et al., 2019). Engineers needed a fast and reliable stress calculation method to evaluate the stress and safety of the FS. Therefore, a quick procedure with precise accuracy and reliability is practical for engineering.

This paper investigates the mechanical analysis methods (MAM) and rapid stress calculation for FS. The disadvantage of the original form is that it does not consider the FS ring length, so propose a new rapid method based on MAM and FEM, which is very helpful in computing the FS’s stress. After that, the finite element model for each factor is established; building each type of matrix in the proposed model is guided in detail. Next, a complete FS design is designed and manufactured, and through the proposed mathematical model, the deformation and stress of the FS have been obtained. Finally, the factor of safety of FS is also calculated, analyzed, and compared with the simulation results.

2. FEM model calculates flexspline stress

2.1. The basis of the method applied to the FS

The FEM divides the continuous geometric model into discrete elements and replaces them with displacement functions. The properties of the mechanical system, including the node displacement, the node force, and the stiffness, are described in a system of algebraic equations as follows:

(1) $K=1pt.\left[\delta \right]=\left[Q\right]$(1)

where:

$K=1pt$ is the stiffness matrix (SM) of the element system;

$\left[\delta \right]$ is the displacement matrix of all nodes in the model;

$\left[Q\right]$ is the node force matrix of all nodes in the model.

In the loaded state, the stress falls mainly at the root of the FS, since that is where the material concentration is the lowest relative to other points. Assume the law of the external force acting on the FS, as shown in Figure 1(a). It is possible to consider an element distributed over one tooth, ignoring the difference of the force distributed at both ends, to see that each element is a uniformly distributed force. There is a jump of the distributed force at the node between two elements (Figure 1(b)).

Figure 1. FEM model calculates FS stress.

Experiments with isotropic materials have shown that the dark bands (representing deformation) are also concentrated mainly on the gear’s root. The FS has a number of teeth over 90 in the assembled state with a WG whose deformation does not exceed 1/5 of the tooth height (Ginz Bun, 1966; Seit Lin & Sukermam, 1972). Furthermore, these deformations are only within the root clearance, not affecting the engagement process. In the worst case, the deformation on the teeth does not exceed 1/3 of the teeth height, and this strain height decreases as the number of FS’s teeth increases. The FS’s profile change when working is described as shown in Figure 2. The FS’s deformation in the direction of deformation caused by the WG is up to 2.m in magnitude, where m is the modulus of the FS that does not affect the change in the tooth’s original shape in the working area.

Figure 2. The change of the FS’s profile when working.

With the parameters of the material’s mechanical properties and the geometric relationship of the elements on the FS after deformation, there is enough basis for building the stiffness matrix of the element system. After determining the geometry of the FS’s neutral curve, based on the assumption of the above distribution of elements, the displacement matrix of the element system can be determined compared to before being affected by force.

Therefore, the “initial preliminary” solving process is completed after calculating the node force from Equation (1)(1) $K=1pt.\left[\delta \right]=\left[Q\right]$(1) . After that, there is enough data to solve the forward force problem for the displacement and stress of the mechanical system. The FEM software will solve the problem with the input data of the node force calculated from the above problem. Then the safety coefficient will confirm the functional ability of the FS. The following introduces the basis for determining the quantities present in Equation (1)(1) $K=1pt.\left[\delta \right]=\left[Q\right]$(1) .

2.1.1. Establish the stiffness matrix of the element system

$\bullet$ Recommend model

The FS shape often has two symmetry axes (over two symmetry axes are rarely used, so it is not considered). How to model general and show the relationship between the rest is a problem that needs to be solved. The following model is suggested for consideration:

The FS has only vertical and horizontal displacements when mounted on a WG. Therefore, the model’s first and last two elements must represent the connection with the rest of the FS. They are represented as a rotation and translational joint, with the translation direction coincident with the direction of the symmetry axis (Figure 3). This generates a link backlash at the node with the removed part. All remaining elements are represented as rotation joints. Because these rotation joints are in space, so each node has three displacements: vertical, horizontal, and angular.

Figure 3. FS node force model.

The elements in the kinematic series, apart from the angle in the standard frame of reference, are different. However, the mechanical and physical parameters, such as EF, EJ, and length L, are the same. Furthermore, the node-numbered elements continuously increase in one direction. Therefore, the calculated stiffness matrix is a band matrix (Table 1). These observations make the data entry process less time-consuming.

Table 1. The stiffness matrix is obtained by adding the extended matrices

	Column 1	Column 2	Column 3		Column n-1	Column n
	of matrix	of matrix	of matrix		of matrix	of matrix
Row 1 of matrix	$K_{11}^1$	$K_{12}^1$
Row 2 of matrix	$K_{21}^1$	$K_{22}^1+K_{22}^2$	$K_{23}^2$
Row 3 of matrix	$K_{32}^2$	$K_{33}^2+K_{33}^3$
				$\ddots$
Row n- 1 of matrix					$K_{n,n}^{n-1}+K_{n,n}^n$	$K_{n,n+1}^n$
Row n of matrix					$K_{n+1,n}^n$	$K_{n+1,n+1}^n$

$\bullet$ The local reference and global reference systems

With the parameters of the material’s mechanical properties and the geometric relationship of the elements on the FS after deformation, there is enough basis for building the stiffness matrix of the element system. After determining the geometry of the FS’s neutral curve, based on the assumption of the above distribution of elements, the displacement matrix of the element system can be determined compared to before being affected by force.

In Figure 1(b), the elements are denoted $e_1,e_2,\dots e_z$ respectively, in the anti-clockwise direction. Since the pitch on the neutral curve is constant, naturally, the elements formed based on this curve have the same length, the common symbol is $l_e$.

Choose z local reference system according to the right-hand rule, the origin of the reference system is from $O_1{x}_1{y}_1$ to $O_z{x}_z{y}_z$ as specified, which is located at the node with the smaller index of the element under consideration. The positive direction of the Ox axis is always from the small node to the large node of the element. The Oy axis is self-determining positive in the direction of positive rotation around Oz of the system. The global reference system of the chosen system coincides with the reference system with the smallest index of the local reference systems mentioned above.

There are two ways to convert the reference systems to each other: the direction cosine of the elements or the rotation matrix. We can choose the second way for convenience because the transpose components in the displacement matrix will be present more fully.

According to Figure 1(b), the angle between two consecutive elements is equal to that of two consecutive teeth in the distribution curve.

(2) $\theta =\frac{{360}^0}{z_{td}}$(2)

where $z_{td}$ is not the FS’s number of teeth, but the gear’s number of teeth that has a module equivalent to the FS, and the radius equal to the radius of the oval equivalent to the arc in consideration. (In the case of using a driving disc, the contact angle is as large as ${60}^0$, so this arc is naturally controlled by the driving disc and has the form of a circular arc, whose radius is derived from the radius of the driving disc).

The major arc radius of the oval must be calculated from smooth transition conditions and conditions such as the eccentricity of the two driving discs (details of the calculation are presented in Section 3.2).

The relationship of the stiffness matrices in the global and local reference system is built (Tinh, 1999):

(3) $\left[\bar{k}\right]=\left[R_e^T\right].\left[k\right].\left[R_e\right]$(3)

$\bullet$ Determine the mechanical properties of the elements such as: l, EF, EJ

$\bullet$ Determine the SM of the elements in the local reference system

The form of the SM of the bar element with two rotating ends in the local reference system can be found in the reference (Tinh, 1999). Unlike the elements on the inside, two elements on the model’s boundary have the form of a link as a rotation and translation joint. However, the stiffness matrix form of the bar with a connection like this is not provided. Hence, it is necessary to base on the fixed beam’s stiffness matrix form to build the bar’s stiffness matrix with a rotation and translational joint.

Consider a bar with a minor node (symbol i) connecting a rotation joint and a major node (symbol j) connecting a translational joint, as shown in Figure 4. The relationship between the internal forces of the two nodes and the displacements of the fixed beam is given by the following system of equations:

Figure 4. Bar with a rotation joint and a translation joint.

(4) $\left\{\begin{array}{c}{X}_i=-\frac{EF}{l}(u_j-u_i)\\ Y_i=-\frac{12EJ}{l^3}(v_j-v_i)+\frac{6EJ}{l^2}({\varphi }_i+{\varphi }_j)\\ M_i=-\frac{6EJ}{l^2}(v_j-v_i)+\frac{2EJ}{l}(2{\varphi }_i+{\varphi }_j)\\ X_j=-\frac{EF}{l}(u_i-u_j)\\ Y_j=-\frac{12EJ}{l^3}(v_i-v_j)-\frac{6EJ}{l^2}({\varphi }_i+{\varphi }_j)\\ M_j=-\frac{6EJ}{l^2}(v_j-v_i)+\frac{2EJ}{l}(2{\varphi }_j+{\varphi }_i)\end{array}\right.$(4)

At the minor node i, bending moment $M_i=0;M_j=0$, shear force $Y_j=0$, assign these values to the equations (4)(4) $\left\{\begin{array}{c}{X}_i=-\frac{EF}{l}(u_j-u_i)\\ Y_i=-\frac{12EJ}{l^3}(v_j-v_i)+\frac{6EJ}{l^2}({\varphi }_i+{\varphi }_j)\\ M_i=-\frac{6EJ}{l^2}(v_j-v_i)+\frac{2EJ}{l}(2{\varphi }_i+{\varphi }_j)\\ X_j=-\frac{EF}{l}(u_i-u_j)\\ Y_j=-\frac{12EJ}{l^3}(v_i-v_j)-\frac{6EJ}{l^2}({\varphi }_i+{\varphi }_j)\\ M_j=-\frac{6EJ}{l^2}(v_j-v_i)+\frac{2EJ}{l}(2{\varphi }_j+{\varphi }_i)\end{array}\right.$(4) and derive the corresponding quantities ${\varphi }_i;{\varphi }_j;v_j$ replacing all the remaining positions in the system. We have the relationship between the internal force and the displacement of the minor and major nodes. Then, converting the relationship to the formation matrix, the stiffness matrix obtained is as follows:

(5) $Stiffnessmatrix=\left|\begin{array}{cccccc}\frac{EF}{l}& 0& 0& -\frac{EF}{l}& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ -\frac{EF}{l}& 0& 0& \frac{EF}{l}& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(5)

Do the same for the bar with a major node connection, the result is the same. So, it can be seen that this is also the SM form of the bar with two ends of the rotation and the translational joint. The stiffness matrix of all elements present in the model is the same, so it is very convenient for calculations because these are sparse matrices and most coefficients are zero.

2.1.2. Establish the displacement matrix of the element system

Building displacement matrix will be the basis for calculating the intensity of the external force applied to each segment of the FS. The displacement determination is based on the measured data directly on the FS after putting it into a working state.

Figure 5 shows how to determine the displacement matrix of an FS with Z = 98, m = 1 (design details are presented in Section 3.1)

Figure 5. The process of building the displacement matrix of the element system.

Choose a reference system that does not coincide with the first element of the system, following the right-hand rule with an axis, as shown in Figure 5. Since the FS has two axes of symmetry, only 1/4 of the FS ring is considered, which reduces the size of the stiffness matrix. As shown in Figure 5(b), after determining the model with only a quarter of FS teeth, the elements are defined as connecting nodes as shown in Figure 5(c) by identical circles. The pre-deformation state divide command allows evenly dividing of the gear ring, replacing curved arcs with straight segments so that angular displacement can be measured, as shown in Figure 5(d). The displacement matrix shown here includes the signs of the quantities according to the standard frame of reference. The length displacement dimension is mm, and the dimension of the angular displacement is decimal degrees, measured with four zeros after the decimal point.

2.1.3. The reasonableness of the proposed model

With external gears, it is possible to experience tooth break due to load. However, with wave drive, the number of teeth in simultaneous engagement is much larger, and the total bending area of all teeth engaged at a time is much larger than the bearing area in a tooth groove.

Statistical practice shows that the major form of FS failure occurs at the root of the tooth, more precisely, the metal cup in which the teeth are distributed.

Because the modulus of transmission determines the thin wall thickness, a strong enough tooth cup is necessary and sufficient to ensure that the modulus has been appropriately chosen.

The computational model can be expanded if it is desired to add different force components to the model. For example, the modeling of the transmission of force on the teeth reduced to the root, which is represented by the Equation (6)(6) $\sum _{i=1}^n{q}_{i1}=\left[\begin{array}{c}0\\ \overline{Q_1}\\ 0\\ 0\\ \overline{Q_2}\\ 0\\ 0\\ \overline{Q_3}\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ \overline{Q_n}\\ 0\\ \hspace{0.17em}\end{array}\right]=\left[\begin{array}{c}0\\ \overline{{Q^1}_1}\\ 0\\ 0\\ \begin{array}{cc}\overline{{Q^1}_2}& +\overline{{Q^2}_2}\\ \hspace{0.17em}& \hspace{0.17em}\end{array}\\ 0\\ 0\\ \begin{array}{cc}\overline{{Q^2}_3}& +\overline{{Q^3}_3}\\ \hspace{0.17em}& \hspace{0.17em}\end{array}\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ \overline{Q_n}\\ 0\\ \hspace{0.17em}\end{array}\right]=\left[\begin{array}{c}0\\ \overline{{Q^1}_1}\\ 0\\ 0\\ \begin{array}{cc}\overline{{Q^1}_1}& +\overline{{Q^2}_2}\\ \hspace{0.17em}& \hspace{0.17em}\end{array}\\ 0\\ 0\\ \begin{array}{cc}\overline{{Q^2}_2}& +\overline{{Q^3}_3}\\ \hspace{0.17em}& \hspace{0.17em}\end{array}\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ \overline{Q_n}\\ 0\\ \hspace{0.17em}\end{array}\right]=\left[\begin{array}{c}0\\ \overline{{Q^1}_1}\\ 0\\ 0\\ \overline{{Q^1}_1}\\ 0\\ ⋮\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \hspace{0.17em}\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ \overline{{Q^2}_2}\\ 0\\ 0\\ \overline{{Q^2}_2}\\ 0\\ ⋮\\ 0\\ 0\\ 0\\ 0\\ 0\\ \hspace{0.17em}\end{array}\right]+\dots \dots +\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ ⋮\\ 0\\ \overline{{Q^{n-1}}_{n-1}}\\ 0\\ 0\\ \overline{{Q^{n-1}}_{n-1}}\\ 0\\ 0\\ 0\\ 0\\ \hspace{0.17em}\end{array}\right]+{\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ ⋮\\ 0\\ \overline{Q_n}\\ 0\\ 0\\ \overline{Q_n}\\ 0\\ \hspace{0.17em}\end{array}\right]}_{\hspace{0.17em}}$(6) node force matrix:

(6) $\sum _{i=1}^n{q}_{i1}=\left[\begin{array}{c}0\\ \overline{Q_1}\\ 0\\ 0\\ \overline{Q_2}\\ 0\\ 0\\ \overline{Q_3}\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ \overline{Q_n}\\ 0\\ \hspace{0.17em}\end{array}\right]=\left[\begin{array}{c}0\\ \overline{{Q^1}_1}\\ 0\\ 0\\ \begin{array}{cc}\overline{{Q^1}_2}& +\overline{{Q^2}_2}\\ \hspace{0.17em}& \hspace{0.17em}\end{array}\\ 0\\ 0\\ \begin{array}{cc}\overline{{Q^2}_3}& +\overline{{Q^3}_3}\\ \hspace{0.17em}& \hspace{0.17em}\end{array}\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ \overline{Q_n}\\ 0\\ \hspace{0.17em}\end{array}\right]=\left[\begin{array}{c}0\\ \overline{{Q^1}_1}\\ 0\\ 0\\ \begin{array}{cc}\overline{{Q^1}_1}& +\overline{{Q^2}_2}\\ \hspace{0.17em}& \hspace{0.17em}\end{array}\\ 0\\ 0\\ \begin{array}{cc}\overline{{Q^2}_2}& +\overline{{Q^3}_3}\\ \hspace{0.17em}& \hspace{0.17em}\end{array}\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ \overline{Q_n}\\ 0\\ \hspace{0.17em}\end{array}\right]=\left[\begin{array}{c}0\\ \overline{{Q^1}_1}\\ 0\\ 0\\ \overline{{Q^1}_1}\\ 0\\ ⋮\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \hspace{0.17em}\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ \overline{{Q^2}_2}\\ 0\\ 0\\ \overline{{Q^2}_2}\\ 0\\ ⋮\\ 0\\ 0\\ 0\\ 0\\ 0\\ \hspace{0.17em}\end{array}\right]+\dots \dots +\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ ⋮\\ 0\\ \overline{{Q^{n-1}}_{n-1}}\\ 0\\ 0\\ \overline{{Q^{n-1}}_{n-1}}\\ 0\\ 0\\ 0\\ 0\\ \hspace{0.17em}\end{array}\right]+{\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ ⋮\\ 0\\ \overline{Q_n}\\ 0\\ 0\\ \overline{Q_n}\\ 0\\ \hspace{0.17em}\end{array}\right]}_{\hspace{0.17em}}$(6)

3. Cases study and analysis

3.1. FS and CS geometry design

A practical case is design to illustrate the method:

The HD reduction ratio is 50, the module m = 1 and having the involute tooth profile. The main geometric design parameters are as follows:

Dual wave drive: ${\alpha }_{\mathrm{\partial }}={20}^o;i=50;Z_F=98;Z_C=100;m=1$

Thin wall thickness coefficient: $\upsilon =1$

FS profile shift coefficient: $\xi =0,022.Z_F=0,022.98=2,156$

Transmitter conversion module:

(7) $m_y=\frac{m.Z_C}{2\sqrt{{k^2}_S+1,3946.k_S+2,7255}}=\frac{1.100}{2\sqrt{{k^2}_S+1,3946.k_S+2,7255}}=0.9615(mm)$(7)

The radius of the neutral circle of the FS

(8) $r_{dt}=0,5.Z_F.m_y=0,\mathrm{5.98.0},9615=47,13(mm)$(8)

The diameter of the pitch circle of the FS $D_{cF}$$:$

$D_{cF}$=Z_F.m = 98.1 = 98(mm)

The radius of the end curve of the FS neutral curve $r_0$:

(9) $r_0=\frac{m_y.k_{\mathrm{\Delta }}}{0,214602}=\frac{0,9615.0,8197}{0,214602}=3,675(mm)$(9)

The radius of the addendum of the FS $R_{eF}$:

(10) $\begin{array}{rl}& R_{eF}=R_{eC}+0,10423r_0\\ & =48,61+0,10423.3,675=48,99(mm)\end{array}$(10)

Correspondingly with above the FS, the geometric parameters of circular spline (CS) are as follows:

$Z_c=100;m=1;f_0=1;c_0=0.2$

The diameter of the pitch circle of the CS $D_{cC}$

(11) $D_{cC}=Z_c.m=100.1=100(mm)$(11)

The radius of the addendum of the CS $R_{eC}$

(12) $\begin{array}{rl}& R_{eC}=m_y\left(\frac{Z_{Me}+\upsilon }{2}+k_{\mathrm{\Delta }}+c_0\right)\\ & =0,9615\left(\frac{98+1}{2}+0,8197+0,2\right)=48.61\end{array}$(12)

The diameter of the dedendum of the CS $D_{iC}$

(13) $D_{iC}=100+(1,25.1).2=102,5(mm)$(13)

The above geometric parameters allow us to determine the basic geometrical parameters of the FS and CS. If machining using the envelope contour method, there is no need to determine more specific factors about the tooth profile. Suppose machining on an EDM requires a detailed drawing of the tooth profile as input for the tool path on the machine. Table 2 shows that Autolist programming allows for defining the detailed drawing of the above tooth profile:

Table 2. Programmatically determine the detailed drawing of the tooth profile

(defun tk ()

(setq p0 (getpoint “center coordinates:”)

r (getdist “radius:”)

alpha 0

m1 (polar p0 (*90(/pi 180))r))

(command “circle” p0 r”“)

(command “pdmode” 34)

(command “point” m1)

(white ($<=$ alpha 30)

(setq p1 (polar p0 (*(/pi 180)(+alpha 90)) r)

p2 (polar p1(*(/pi 180) alpha)(*r alpha(/pi180))))

(command “pline” m1 p2”“)

(set m1 p2)

(setq alpha (+alpha 0.1))

))

A complete tooth profile is obtained after running the program from the Auto-CAD environment. Then, based on the above-calculated parameters, complete the drawing of the FS and CS as shown in Figure 6.

Figure 6. Completed FS and CS design drawing.

3.1.1. Application of the proposed method to calculate stress for FS

This section shows how to calculate the stress acting on FS with the number of teeth Z = 98 and the module m = 1 as in the design above.

Using the EDM method to machine the FS, after mounting the FS with the WG, the shape of the neutral curve can be determined using a three-dimensional coordinate measuring machine (CMM), as shown in Figure 7(a). The measured sample result has two semi-major axes and two semi-minor axes of oval shape, Figure 7(b).

Figure 7. The design, model, and neutral curve measurement results of the FS.

When taking the semi-axes of the FS to construct the oval instead, the list command of AutoCAD gives us the circumference of the scanned image as 302,7204(mm), giving the circumference of the oval shape as 301,9934(mm), so the oval equivalent is acceptable (Figure 7(c)). The list command for the circumference of arc AB of the oval is 70,1180(mm), and for the circumference of the arc, BC is 80,8787(mm). Because they are all arcs, so the teeth of the FS on each arc will be evenly distributed around its center, based on the relation of corresponding sides perpendicular.

$\bullet$ The angle between two consecutive teeth (two elements) on the AB arc is:

(14) $\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}=\frac{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{s}}=\frac{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}{\left[\frac{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}}{\mathrm{t}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{c}\mathrm{h}}\right]}$(14)

Taking an element through three teeth because of the small gear modulus, the rotation angle of the three teeth is relatively significant. Therefore, the tooth pitch here is three times the actual tooth pitch, with the modulus m = 1 of the FS.

(15) ${\alpha }_1=\frac{{97}^0}{\left[\frac{70.1180}{\mathrm{3.1.3},14}\right]}=13,{031}^0$(15)

Similarly, the angle between two consecutive teeth (two elements) on the AD arc is:

(16) ${\alpha }_2=\frac{{83}^0}{\left[\frac{80.8787}{\mathrm{3.1.3},14}\right]}=9,{6671}^0$(16)

Particularly for the first element, having the x-axis created with the Ox axis of the global reference system $+7^o$ angle, the following rotation matrix should be used:

(17) $R_1=Rot(z;7^0)=\left|\begin{array}{cccc}1& 0& 0& 0\\ 0& 0,992& 0,121& 0\\ 0& -0,121& 0,992& 0\\ 0& 0& 0& 1\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(17)

For the second element, the rotation angle is $7^o+13,{031}^o=20,{031}^o$, the corresponding rotation matrix is:

(18) $R_2=Rot(z;20,{031}^0)=\left|\begin{array}{cccc}1& 0& 0& 0\\ 0& 0,939& 0,342& 0\\ 0& -0,342& 0,939& 0\\ 0& 0& 0& 1\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(18)

On that principle, the subsequent rotation matrices need to add the angle ${\alpha }_1$.

(19) $R_3=Rot(z;33,{062}^0)=\left|\begin{array}{cccc}1& 0& 0& 0\\ 0& 0,838& 0,545& 0\\ 0& -0,545& 0,838& 0\\ 0& 0& 0& 1\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(19)

(20) $R_4=Rot(z;44,{093}^0)=\left|\begin{array}{cccc}1& 0& 0& 0\\ 0& 0,693& 0,72& 0\\ 0& -0,72& 0,693& 0\\ 0& 0& 0& 1\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(20)

From the fifth element onward, the rotation angle needs to be added ${\alpha }_2$

(21) $R_5=Rot(z;55,{76}^0)=\left|\begin{array}{cccc}1& 0& 0& 0\\ 0& 0,562& 0,826& 0\\ 0& -0,826& 0,562& 0\\ 0& 0& 0& 1\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(21)

(22) $R_6=Rot(z;65,{42}^0)=\left|\begin{array}{cccc}1& 0& 0& 0\\ 0& 0,415& 0,909& 0\\ 0& -0,909& 0,415& 0\\ 0& 0& 0& 1\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(22)

(23) $R_7=Rot(z;75,{08}^0)=\left|\begin{array}{cccc}1& 0& 0& 0\\ 0& 0,257& 0,966& 0\\ 0& -0,966& 0,257& 0\\ 0& 0& 0& 1\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(23)

(24) $R_8=Rot(z;84,{74}^0)=\left|\begin{array}{cccc}1& 0& 0& 0\\ 0& 0,091& 0,995& 0\\ 0& -0,995& 0,091& 0\\ 0& 0& 0& 1\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(24)

$\bullet$ Determine the mechanical properties of the elements: The FS designed is made of 40XC material and has the following properties:

(25) $l=3.m.\pi =9,42E-3(m)$(25)

Steel FS with the module of elasticity:

$E=2,1E11(N/m^2)$

Cross-sectional area: $F=1,16E-4(m^2)$

Based on the drawing, the Massprop command gives us: $J=3,75E-6(m^4)$; $EF=2,436E3(t)$; $EJ=0,078E3(t.m^2)$

After substituting numbers according to the mechanical properties specified above, all elements of the system have the form of a stiffness matrix associated with the local reference system as follows:

(26) $\left|K_i\right|=\left|\begin{array}{cccccc}258,598& 0& 0& -258,598& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 258,598& 0& 0& 258,598& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right|$(26)

The SM of the element system is built based on this matrix, and the rotation matrices from R1 to R8 mentioned above is a matrix of size 24*24 obtained by adding the elements. The construction of the SM of the element system is completed.

$\bullet$ Determine the matrix of displacement as presented in Section 2.3. The result is shown in Figure 5(d).

$\bullet$ Substituting the displacement matrix and the stiffness matrix into the Equation (1)(1) $K=1pt.\left[\delta \right]=\left[Q\right]$(1) , we get the node stress to result as follows:

(16,0764; 1,9598; 0; 57,7732; 24,9346; 0; 4,6034; 23,6091; 0; 85,8433; 119,3588; 0; −123,8549; 34,4285; 0; −72,3440; −58,2948; 0; −43,3082; −58,2622; 0; −23,3412; −87,7338; 0; 0; 0; 0)^T

3.1.2. Calculation and simulation to check the safety factor of the FS

As analyzed in Section 2.1, the safety coefficient will confirm the functional ability of the FS. Therefore, this section will calculate and simulate the functional ability of FS based on the above stress analysis results to verify the proposed method results.

The FS ring is part of working under heavy conditions, so its safety factor must be greater than 1,5. The following are the two main test procedures corresponding to the results of the geometrical and mechanical parameters of the above FS design.

$\bullet$ Check according to the traditional model

Specifically, the data of the FS designed above is made of 40XC material, FS is the active link$M_F=120(kg.m);i=50;HB=300;L=d_{cC};q_F=0,15;[p]=600(kg/c{m}^2)$$\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e},{\alpha }_{\mathrm{\partial }}={20}^0$. Incorporation of the FS geometry data when calculating kinematics: $D_{cC}=100(mm);m=1(mm);m_y=0,9615(mm);\vartheta =1;\xi =2,185;r_{dt}=47,13(mm)$

Determine the minimum normal stress acting on the wall of the suture on the tooth distribution segment:

(27) ${\sigma }_z=\frac{\delta E}{m_y{z}_F(0,1366\frac{z_{Me}}{k_{\mathrm{\Delta }}}-1)}=\frac{{0,3.2,1.10}^6}{0,9615.98(0,1366\frac{98}{0,8197}-1)}=427(kg/c{m}^2)$(27)

The value $\lambda { }^{\mathrm{\prime }}$ determined according to the graph is 2,75.

(28) ${\sigma }_{max}={\sigma }_z\lambda { }^{\mathrm{\prime }}=427.2,75=1174,25(kg/c{m}^2)$(28)

(29) ${\sigma }_{min}=-\frac{\delta E{{\alpha }^{\prime }}_{\sigma }}{z_F{m}_y(0,1366\frac{z_{Me}}{k_{\mathrm{\Delta }}}-1)}=-\frac{{0,3.2,1.10}^6.1,2}{98.0,9615(0,1366\frac{98}{0,8197}-1)}=-512(kg/c{m}^2)$(29)

The value of torsional stress on the wall of the FS in the limited groove portion of the toothed section is calculated as follows:

(30) ${\tau }_{max}=\frac{2M_{Me}}{\pi {d^2}_{cMe}\delta }=\frac{2.12000}{{3,14.9,8}^2.0,3}=265,2(kg/c{m}^2)$(30)

The frequency and cyclic stress are calculated:

(31) ${\sigma }_a=\frac{1}{2}({\sigma }_{max}-{\sigma }_{min})=\frac{1}{2}(1174,25+512)=843,125(kg/c{m}^2)$(31)

(32) ${\sigma }_m=\frac{1}{2}({\sigma }_{max}+{\sigma }_{min})=\frac{1}{2}(1174,25-512)=331,125(kg/c{m}^2)$(32)

(33) ${\tau }_a={\tau }_m=\frac{1}{2}{\tau }_{max}=\frac{265,2}{2}=132,6(kg/c{m}^2)$(33)

With 40XC crafting materials there are ${\sigma }_{-1}=4600(kg/c{m}^2);{\tau }_{-1}=2700(kg/c{m}^2)$.

The following formula determines the factor of safety under normal stress.

(34) $n_{\sigma }=\frac{{\sigma }_{-1}}{\frac{{\sigma }_a}{{\epsilon }_{\sigma }{\beta }_{\sigma }}+{\sigma }_m{\psi }_{\sigma }}=\frac{4600}{\frac{843,125}{0,6.0,9}+331,125.0,1}=4,885$(34)

(35) $n_{\tau }=\frac{{\tau }_{-1}}{0,5{\tau }_{max}(\frac{{{\alpha }^{\prime }}_{\tau }}{\epsilon {\beta }_{\tau }}+{\psi }_{\tau })}=\frac{2700}{0,5.265,2(\frac{1,95}{0,68.0,9}+0,05)}=6,292$(35)

Therefore, the combined factor of safety will be:

(36) $n=\frac{n_{\sigma }{n}_{\tau }}{\sqrt{{n_{\sigma }}^2+0,7{n_{\tau }}^2}}=\frac{4,885.6,292}{\sqrt{{4,485}^2+{0,7.6,292}^2}}=4,279$(36)

This safety factor is more significant than 1,5 to satisfy the requirements of working durability; the calculation is finished.

$\bullet$ Test with Cosmos Design Star simulation software

To perform the test, build a 3D model of the FS on AutoCAD Mechanical Desktop. Export files as (*.iges) or (*.igs). Import the model into the Cosmos Design Star environment, and do the following:

- Determine the material’s mark through the physical and mechanical parameters of the material.

- Set force based on the result of the node stress

- Control the displacement while working.

- Discrete continuous model of FEM model.

- Run interactive solver program.

- Select the test standard.

- Evaluation of result.

The software allows checking the factor of safety according to three basic criteria:

Mohr—Coulomb stress: FOS = 4,4

Test formula:

(37) $\frac{{\sigma }_1}{{\sigma }_{TensiveLimit}}+\frac{{\sigma }_2}{{\sigma }_{CompressiveLimit}}<1$(37)

Maximum Von Mises stress: FOS = 4,4

Test formula:

(38) $\frac{{\sigma }_{Vonmises}}{{\sigma }_{Limit}}<1$(38)

Maximum shear stress (Tresca): FOS = 4,4

Test formula:

(39) $\frac{{\tau }_{max}}{0,5{\sigma }_{Limit}}<1$(39)

Figure 8 shows the FS simulation results, specific deformation and stress contour plots.

Figure 8. Computational simulation results on cosmos design star

In Figure 8, the stress distribution of FS perfectly agrees with the proposed mathematical model. The deformation causing the internal bending moment has been discussed in Section 2, which is used to compute the stress via the MAM. According to the simulation results, the most significant stress occurs at the tooth root on the central axis of WG, which is also wholly consistent with many current research results and utterly correct according to the initial assessment of the research team. Soon, to further demonstrate the feasibility of the novel method, the research team will increase the number of fabricated teeth of FS to study the relationship between stress reduction and deformation at the tooth root.

The test results for the stress safety factor are all 4,4 quite close to the results calculated by theory 4,279. It shows that the proposed method helps quickly calculate the FS’s stresses. Therefore, engineers have one more reference method for stress assessment of FS during harmonic actuators’ design and structural optimization.

4. Conclusions

In this research, the stress computation methods for FS have been analyzed, a harmonic reducer was designed and manufactured, and the stress in FS was analyzed according to the proposed method. Some primary conclusions can be drawn as follows:

- The usual MAM can’t fully express the influence of the material’s stiffness and the FS’s displacement on stress. Therefore, a quick method based on the FEM model was proposed. The method obtained the analytical relationship between FS’s stiffness, deformation, and stress.

- Building a computational model of the stiffness matrix for the bar element at the border has connecting rotation and translational joints. Establish a complete model and a stress calculation process for the entire FS.

- The stress and safety factors of FS have been solved in the proposed model, and the simulation results are entirely consistent with the theoretical analysis results, which is approximately 4,4 almost three times greater than the standard safety factor of 1,5. The simulation results also show that the stress is most concentrated in the tooth root, which is the leading cause of fatigue failure of the FS. Furthermore, the experiment’s results on the actual FS set verified that the quick stress computation method can reach accurate results.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The research is supported by the China Postdoctoral Science Foundation, Grant No. 2021M703780; Postdoctoral Science Foundation of China.

Notes on contributors

ThanhTrung Trang

ThanhTrung Trang (1983, Male), PhD. He received a B.S and M.S in Mechanical Engineering from the Faculty of Mechanical Engineering, Thai Nguyen University of Technology, Vietnam in 2006 and 2011, respectively. In 2018, He received a Ph.D. degree in mechanical engineering from the Mechanical and Automotive Engineering, South China University of Technology, China. From 2018 to 2022, he was a Postdoctoral Researcher at the School of Automation Science and Engineering, South China University of Technology. His research interest covers mechatronics, robotics, automation in manufacturing, and optimization methods. He has published more than 18 international academic papers and is the main responsible for three research topics at the provincial and national levels in China.