Introduction

Spur bevel gears are key components in gear transmission systems and are widely used in various fields such as automobiles, precision machine tools, aerospace, etc. The time-varying meshing stiffness is an important basis for evaluating the transmission performance of bevel gears, and accurate calculation of its time-varying meshing stiffness is an important prerequisite for ensuring the stability of the transmission performance of bevel gears. Therefore, it is of great significance to calculate the time-varying meshing stiffness of bevel gears efficiently and accurately.

Many scholars at home and abroad have conducted a lot of research on the calculation of time-varying meshing stiffness of gears. There are more studies on cylindrical gears, and the problem of calculating their meshing stiffness has basically been solved. However, there are fewer studies on the time-varying meshing stiffness of spatial bevel gears, and most scholars use the finite element method to calculate it. Although the finite element method has high calculation accuracy, it has a long calculation time and the calculation efficiency is limited by many factors. The traditional analytical calculation method has a short calculation time, but it cannot consider the coupling effect between tooth pairs, which will lead to a large calculation error of multi-tooth meshing stiffness.

In response to the above problems, this paper combines the advantages of the finite element method and the analytical method, and proposes a slice calculation method for the meshing stiffness of spur bevel gears considering the coupling effect of tooth pairs. The meshing stiffness, the angular displacement of the center point of each slice, and the load distribution of the spur bevel gears under the assembly error conditions are calculated and analyzed, and the accuracy of the method proposed in this paper is verified by comparing with the finite element calculation results.

**1. Meshing Stiffness Calculation Model of Spur Bevel Gears**

**1.1 Spur Bevel Gear Slice Model**

Considering the three-dimensional spatial structure and spatial meshing characteristics of the spur bevel gear, the gear tooth of the bevel gear is discretized along the tooth width direction by the back cone, and the gear body is discretized along the circular plane of the axis direction to form a series of thin slice bevel gears with the same tooth width. Based on the back cone equivalence principle, the large end back cone of each thin slice bevel gear is expanded. When the width of the thin slice is small enough, the bevel gear can be approximated as a series of spur gears with the same number of teeth (but different sizes). The slice model of the spur bevel gear is shown in Figure 1.

Symbol | Definition |
---|---|

d6 | Width of the thin slice gear tooth |

ri, r5, rij (j = 1, 2) | Radius of the pitch circle, radius of the addendum circle, and radius of the dedendum circle of the bevel gear (i and j represent the driving gear and the driven gear, respectively) |

δj, B, R, rij | Cone angle, tooth width, outer cone distance, and radius of the shaft hole of the bevel gear |

Based on the back cone equivalence principle, the equivalent number of teeth zij, the equivalent radius of the pitch circle rq’i, the equivalent radius of the addendum circle rqi, and the equivalent radius of the dedendum circle ryi of the ith thin slice can be expressed as:

where i is the thin slice number, i = 1, 2,…, n; z is the number of teeth; δa and δf are the tooth top angle and tooth root angle of the bevel gear, respectively.

In addition, the gear body part of the thin slice bevel gear is divided into a series of thin gear body slices with a small enough width along the circular plane perpendicular to the axis direction. The relationship between the width di of the thin slice gear body, the width ds of the gear tooth, the tooth width b, and the number of thin slices n can be expressed as:

**1.2 Slice Meshing Stiffness Calculation**

Each thin slice can be approximated as a spur gear, and the spur bevel gear can be approximated as a number of spur gears connected in parallel. The stiffness of the thin slices participating in the meshing at each instant is obtained, and finally the total meshing stiffness at each meshing position can be obtained by accumulation. Under the action of an external load, the equivalent meshing stiffness of the gear caused by the deformation of the gear tooth includes the Hertz contact stiffness, bending stiffness, shear stiffness, and radial compression stiffness. In addition, there is also the equivalent meshing stiffness caused by the deformation of the gear body, namely the gear body stiffness. Therefore, the single-tooth linear meshing stiffness kik^ij of the ith thin slice can be expressed as:

where kh, kb, ks, k., and kf are the contact stiffness, bending stiffness, shear stiffness, radial compression stiffness, and gear body stiffness, respectively. The Hertz contact stiffness can be calculated by Equation (10) in Reference [9], the gear body stiffness can be calculated by Equation (17) in Reference [10], and the bending stiffness, shear stiffness, and radial compression stiffness can be calculated by Equations (2.26), (2.30), and (2.32) in Reference [11], respectively. The superscripts 1 and 2 represent the driving gear and the driven gear, respectively.

The traditional analytical method calculates the multi-tooth meshing stiffness by adding the single-tooth meshing stiffness of the tooth pairs, which has a large calculation error. This is because the traditional analytical method only considers the single-tooth meshing situation when calculating the gear body stiffness, ignoring the coupling effect between the tooth pairs. To this end, this paper introduces a correction coefficient of the gear body stiffness in the calculation of the double-tooth meshing stiffness, and expresses the double-tooth linear meshing stiffness kli’ as:

where j represents the meshing tooth pair; λj’ and λj^2 are the correction coefficients of the gear body stiffness of the driving wheel and the driven wheel, respectively, and this coefficient will be determined in Section 1.3.

In addition, the time-varying meshing stiffness can also be evaluated by the torsional meshing stiffness. According to the relationship between the linear meshing stiffness and the torsional meshing stiffness, the torsional meshing stiffness kt of the bevel gear can be expressed as:

where rbn is the base circle radius corresponding to the ith thin slice.

**1.3 Calculation of the Gear Body Stiffness Correction Coefficient**

Considering the three-dimensional spatial structure of the bevel gear, a three-dimensional finite element model of each thin slice bevel gear is established in the finite element software Ansys. The finite element model of the thin slice i is shown in Figure 2.

Symbol | Definition |
---|---|

E – 10000E | Elastic modulus of the gear tooth set as 10,000 times that of the gear body |

O1 – x1y1z1, O2 – x2y2z2 | Local coordinate systems of the meshing points |

Fij | Force applied at the meshing position Oij (i = 1, 2; j = 1, 2) |

δij^1, δij^2 | Displacement along the action line at the meshing position Oij |

To calculate the displacement along the action line induced by the gear body at the meshing point, the elastic modulus of the gear tooth of the bevel gear thin slice is set to 10,000 times that of the gear body. At this time, compared with the gear body, the gear tooth can be regarded as a relatively rigid region, so the displacement caused by the local Hertz contact and the elastic deformation of the gear tooth can be ignored. A mass21 element is established at the midpoint of the tooth width of the axis of the thin slice, coupled with the inner ring nodes and subjected to full constraints. In the figure, the axes x1 and x2 are respectively along the direction of the action line. At the double-tooth meshing position O1j, a force Fij is applied along the axis x1, and at the meshing position O2j, a force F2j (F1j = F2j = 0.5N) is applied along the axis x2, and the displacement generated at Oij (i = 1, 2; j = 1, 2) along the action line direction is measured. According to the finite element results, we can get:

The load distribution coefficients Lsf1 and Lsf2 can be expressed as:

Based on the linear superposition assumption, Kn and Kr2 can be expressed as:

Combining Equations (7) and (8), the gear body correction coefficients λ1 and λ2 can be obtained as:

By calculating the gear body correction coefficient λ1 and λ2 corresponding to each thin slice and substituting them into Equation (5), the double-tooth meshing stiffness of the bevel gear can be obtained. Compared with the finite element model of the bevel gear pair, the thin slice finite element model is very simple and does not need to set the contact pair, the normal penalty stiffness coefficient, the solution method, etc., so most of the calculation time can be saved.

**2. Meshing Stiffness Calculation Model of the Error Tooth Profile**

**2.1 Deviation of Each Thin Slice Tooth Profile Caused by the Assembly Error of the Bevel Gear**

There are many factors that affect the meshing stiffness of the spur bevel gear, in addition to the design deviation and manufacturing error, the installation error has a greater impact on the meshing stiffness of the spur bevel gear. Therefore, in this paper, only the deviation of the tooth profile position (thin slice back cone error tooth profile) caused by the installation error is considered, and its meshing stiffness model is established. According to the national standard, the installation error can be divided into the shaft intersection angle error and the shaft intersection point error. The relationship diagram of the assembly error coordinate system is shown in Figure 3.

Symbol | Definition |
---|---|

θ, ε | Shaft intersection angle error and shaft intersection point error |

S1(O1 – x1y1z1), S2(O2 – x2y2z2) | Coordinate systems fixed with the driving and driven gears, respectively |

S1i(O1i – x1iy1iz1i), S2i(O2i – x2iy2iz2i) | Auxiliary coordinate systems at the ith thin slice |

δ1, δ2 | Angles between z1i, z2i and z1, z2, respectively |

ψ | Angle between the center line of the driving and driven gear tooth slices and the x1i direction of the driving gear slice |

From Figure 3(a), under the shaft intersection angle error, the coordinate system S1 rotates clockwise around the axis x1 by an angle θ to the illustrated position S10′. From Figure 3(b), under the shaft intersection point error, the coordinate system S2 translates along the axis z1 by ε to the illustrated position S2′. Therefore, the tooth profile deviation E11 and E12 caused by the shaft intersection angle and the shaft intersection point of the ith thin slice are:

where δj (j = 1, 2) is the cone angle of the driving and driven wheels; Ri is the outer cone distance of the thin slice i. When the driving gear rotates clockwise, ψ can be expressed as:

where α is the pressure angle; α’ is the angle between the center line of the driving and driven gear tooth slices and the x1i direction of the driving gear slice.

In addition, the rotation deviation e11 and e12 of the ith thin slice caused by the shaft intersection angle and the shaft intersection point are:

where rbn is the base circle radius corresponding to the ith thin slice.

**2.2 Meshing Stiffness Model of the Error Tooth Profile Bevel Gear**

Under the installation error condition, the axis of the bevel gear will shift and tilt, which makes the meshing situation of each thin slice tooth pair different, often leading to edge contact, which seriously affects the transmission performance of the gear. Therefore, the calculation of the meshing stiffness under the error condition needs to judge whether each thin slice participates in the meshing. The load-bearing contact process of the bevel gear thin slice under the error condition is shown in Figure 4.

Symbol | Definition |
---|---|

Ei (i = m, k, l,…) | Tooth profile deviation of the thin slice i |

Ei’ | Sum of the deformation amount and the tooth profile deviation of the thin slice |

ei | Rotation deviation of the thin slice |

e | Sum of the angular displacement and the rotation deviation of the thin slice under the error condition |

T | Applied external torque |

Under the action of the load, as the driving gear rotates, the rotation deviation of the thin slice with the smallest rotation deviation amongthe load, the angular displacement θm of the center point of the thin slice m starts to increase. When the sum of θm and ei is greater than the rotation deviation ek of the thin slice k, the thin slice k starts to contact and deform, bearing the load. If the bearing capacity at this time does not reach the load of the bevel gear, θm will continue to increase until the sum of θm and em is greater than the rotation deviation el of the tooth pair l, and the tooth pair l starts to participate in the meshing. And so on, until the bearing capacity of all the thin slice tooth pairs participating in the meshing is equal to the external torque T received by the bevel gear, the gear is in a state of force balance, and the deformation ends.

According to Figure 4, the angular displacement of the bevel gear after loading can be expressed as:

From Equation (14), when the center point angular displacement of the thin slice tooth pair m is known, the center point angular displacement of the thin slice tooth pair i can be obtained as:

According to the principle of force balance, the bearing torque generated by the whole bevel gear should be the sum of the bearing capacities of each thin slice tooth pair in the meshing state, that is:

where T is the applied external torque; Q is the number of thin slices participating in the meshing under ideal conditions; Hi is the judgment factor, which takes 1 when the thin slice participates in the meshing, and 0 otherwise; ki is the torsional meshing stiffness of the thin slice i; θi is the center point angular displacement of the thin slice i.

The calculation process of the angular displacement of the thin slice under the error condition is shown in Figure 5. Input i and θseo, and by judging whether θi is greater than 0, the thin slices participating in the meshing can be obtained. If the relative error between the bearing capacity of the thin slice tooth pairs participating in the meshing and the external torque T is greater than the acceptable value ξ, increase the step size of θm and continue the loop calculation; otherwise, output θm and θ1, and the loop ends. In this paper, the initial value of i is set to 0, θm0 = 0, Δ = 10^(-5) rad, and ξ = 0.01%.

When the bevel gear pair is regarded as a whole, under the action of an external load, the angular displacement of the center point of the bevel gear is the same as that of the thin slice tooth pair that first contacts, that is, the center angular displacement of the thin slice tooth pair m. Therefore, the comprehensive torsional meshing stiffness of the bevel gear pair under the error condition can be expressed as:

In addition, according to the deformation amount δi of each thin slice, the load borne by each thin slice can be obtained as:

**3. Example Analysis**

Take a pair of spur bevel gears as an example (parameters are shown in Table 2) for finite element analysis, and compare the finite element calculation results with the method in this paper and the traditional analytical method to verify the accuracy of the method in this paper.

Table 2 Spur Bevel Gear Parameters

Parameter | Value | Parameter | Value |
---|---|---|---|

Module/mm | 1.75 | Tip Height Coefficient h | 1 |

Number of Teeth of Driving Gear z1 | 35 | Clearance Coefficient c | 0.25 |

Number of Teeth of Driven Gear z2 | 35 | Hole Radius r/mm | 15 |

Tooth Width b/mm | 10 | Poisson’s Ratio | 0.3 |

Pressure Angle (°) | 20 | Elastic Modulus/CPa | 206 |

**3.1 Finite Element Model of the Spur Bevel Gear**

The finite element model of the spur bevel gear is shown in Figure 6.

Symbol | Definition |
---|---|

mass21 | Element established at the midpoint of the tooth width of the axis of the driving and driven gears |

T | Torque applied to the driving gear |

θ | Angular displacement of the coupling center point of the driving gear |

A mass21 element is established at the midpoint of the tooth width of the axis of the driving and driven gears, and coupled with the inner surface nodes of the gear wheel body; the axial rotation degree of freedom of the driving gear is retained and a torque T = 100 N·m is applied; the driven gear is subjected to full freedom constraints. The normal penalty stiffness coefficient is set to 1, the penetration coefficient is set to 0.1, and the augmented Lagrangian method is used for the solution.

In this model, by calculating 10 meshing positions within one meshing cycle, the torsional meshing stiffness of the spur bevel gear can be obtained, that is:

**3.2 Verification of the Meshing Stiffness of the Spur Bevel Gear**

During the meshing cycle, taking the meshing stiffness when the rolling angle is 0°, that is, when the spur bevel gear just enters the double-tooth meshing, as the benchmark for grid independence verification, as shown in Figure 7.

Symbol | Definition |
---|---|

k1 | Torsional meshing stiffness |

It can be seen from Figure 7 that the convergence of this method is good, and when the number of thin slices is greater than 60, the meshing stiffness value tends to be stable.

When the number of thin slices is 60, the comparison of the meshing stiffness calculation results of the spur bevel gear by the method in this paper, the traditional analytical method, and the finite element method is shown in Figure 8.

It can be seen from Figure 8 that the error of the double-tooth meshing stiffness of the spur bevel gear obtained by the analytical method is relatively large, about 40%; the single and double-tooth meshing regions of the bevel gear meshing calculated by the method in this paper are basically consistent with the finite element method, and the calculated results are also relatively close, with the error of the double-tooth meshing stiffness being about 2% and the error of the single-tooth meshing stiffness being about 3%. This shows that the method in this paper can accurately calculate the meshing stiffness of the bevel gear pair.

The comparison of the calculation time between the method in this paper, the analytical method, and the finite element method is shown in Table 3.

It can be seen from Table 3 that, under the condition of similar calculation accuracy, the calculation time of the method in this paper (including the time for calculating the correction coefficient of the gear body stiffness) is much less than that of the finite element method, about 0.3% of the calculation time of the finite element method; in addition, the accuracy of the double-tooth meshing stiffness calculated by the method in this paper is much higher than that of the analytical method, about 95% higher. In summary, it shows that the method in this paper can quickly and accurately calculate the meshing stiffness of the bevel gear.

Table 3 Comparison of Calculation Time

Calculation Method | Relative Error of Single-Tooth Meshing Stiffness/% | Relative Error of Double-Tooth Meshing Stiffness/% | Calculation Time/min |
---|---|---|---|

Method in This Paper | 3 | 2 | 5 |

Analytical Method | 3 | 40 | 0.4 |

Finite Element Method | – | – | 1440 |

**3.3 Analysis of the Meshing Stiffness and Load Distribution Considering the Assembly Error**

To verify the accuracy of the calculation method in this paper under the assembly error condition, take the shaft intersection point error ε as 80 μm, and the comparison of the meshing stiffness calculation results of the bevel gear by the method in this paper and the finite element method is shown in Figure 9.

It can be seen from Figure 9 that the calculation error of the double-tooth meshing stiffness is about 3%, and the calculation error of the single-tooth meshing stiffness is about 2%, indicating that this method also has accuracy in calculating the meshing stiffness when considering the assembly error.

According to the requirements of the national standard GB/T 11365 – 2019 for the installation accuracy of bevel gears, take the shaft intersection angle error θ as 0°, 0.01, 0.015, and 0.02, and the shaft intersection point error ε as 0, 30, 60, and 90 μm, respectively. The meshing stiffness calculation results of the corresponding spur bevel gear pair are shown in Figures 10 and 11, respectively.

It can be seen from Figures 10 and 11 that as the shaft intersection angle error and the shaft intersection point error of the gear increase, the meshing stiffness of the bevel gear decreases. With the decrease of the meshing stiffness of the bevel gear, the meshing deformation and transmission error of the bevel gear will increase, which will affect the transmission performance and service life of the bevel gear.

Figure 12 shows the angular displacement of the center point of the bevel gear and the load distribution on the tooth surface under the ideal state. It can be seen from Figure 12(a) that in the case of no error, the angular displacement of the center point of each thin slice of the bevel gear first increases and then decreases from the tooth root to the tooth tip. This is because only one pair of tooth pairs bears the load during single-tooth meshing, while two pairs of tooth pairs bear the load during double-tooth meshing, resulting in a larger angular displacement during single-tooth meshing. It can be seen from Figure 12(b) that the place that bears the maximum load is the large end during the single-tooth meshing of the bevel gear, because the meshing stiffness of the large end thin slice is greater than that of the small end, resulting in a larger torque at the large end.

Figure 13 shows the angular displacement of the center point of each thin slice of the bevel gear and the load distribution on the tooth surface when the shaft intersection angle error θ = 0.01°. It can be seen from Figure 13(a) that the angular displacement of the center point of the small end thin slice of the bevel gear under load is the largest, because the rotation deviation of the small end is the smallest under the shaft intersection angle error state, and it participates in the meshing first in the meshing process, resulting in the largest angular displacement among all the participating thin slices. It can be seen from Figure 13(b) that the maximum load borne by the bevel gear is the small end during the single-tooth meshing.

Figure 14 shows the angular displacement of the center point of each thin slice of the bevel gear and the load distribution on the tooth surface when the shaft intersection point error ε = 90 μm. It can be seen from Figure 14(a) that the angular displacement of the center point of the large end thin slice of the bevel gear under load is the largest, because the rotation deviation of the large end is the smallest under the shaft intersection point error state, and it participates in the meshing first in the meshing process, resulting in the largest angular displacement among all the participating thin slices. It can be seen from Figure 14(b) that the place that bears the maximum load is the large end during the single-tooth meshing of the bevel gear.

In summary, the assembly error of the bevel gear causes serious unbalanced loading during the transmission process of the gear, resulting in increased meshing deformation and transmission error, which will reduce the bearing capacity of the spur bevel gear pair and seriously affect its service life. Therefore, in actual engineering, the installation accuracy of the bevel gear pair should be improved as much as possible to ensure its service life.

**4. Conclusion**

In this paper, the following research is carried out for the efficient calculation of the time-varying meshing stiffness of the spur bevel gear:

- Based on the back cone equivalence principle, the bevel gear is equivalent to a series of spur gears, and the energy method is used to calculate the meshing stiffness of each thin slice. In the calculation of the double-tooth meshing stiffness, the coupling effect between the tooth pairs during multi-tooth meshing is considered, and the finite element method is used to introduce the correction coefficient of the gear body stiffness to reduce the calculation error of the double-tooth meshing stiffness.
- By comparing the calculation results of the method in this paper, the analytical method, and the finite element calculation, it is shown that the error of the single-tooth meshing stiffness calculated by the method in this paper is about 3%, the error of the double-tooth meshing stiffness is about 2%, and the calculation time is much less than the finite element calculation time, which verifies the accuracy and efficiency of the method in this paper.
- By solving the angular displacement at the center point of each thin slice tooth pair of the spur bevel gear after loading, the meshing stiffness, the angular displacement of the center point of each thin slice, and the load distribution of the bevel gear under the error condition are calculated; compared with the finite element, the accuracy of the meshing stiffness calculation under the error condition is verified.

The method in this paper combines the advantages of the finite element method and the analytical method, which not only greatly improves the calculation efficiency of the time-varying meshing stiffness of the bevel gear, but also has high calculation accuracy. In addition, it is also verified that the method in this paper is still accurate under the installation error condition. However, it should be noted that in addition to the installation error, the time-varying meshing stiffness of the bevel gear is also affected by factors such as the tooth profile error caused by the surface roughness and the friction and lubrication during the tooth surface meshing, which cannot be considered by the method in this paper and still needs to be further improved.