Theoretical and Experimental Study on Contact Characteristics of Spiral Bevel Gears under Quasi-Static and Large Loading Conditions

Abstract

The precise mathematical model for the tooth surface and transition surface of spiral bevel gears is derived. Taking a pair of spiral bevel gears of a heavy vehicle as an example of calculation and analysis, a finite element model of spiral bevel gears transmission system is established. Through the finite element tooth contact analysis under quasi-static loading and high loading condition, the influences of torque on the root stress distribution, contact stress, and transmission error are discussed, and the results are compared with the empirical formula results. Finally, a contact performance test bench of spiral bevel gear pair is developed, then the root bending stress, contact pattern, and transmission error tests are carried out. These experiment results are compared with analyzed ones, which showed a good agreement.

1. Introduction

Tooth contact analysis (TCA) technology is widely applied for assessing meshing quality duet to the increasing demand for vibration, strength, and mechanical efficiency of gear in various industrial fields such as aviation, automotive, and civil engineering [1,2]. There are some main evaluation items to reflect the tooth contact performances. For instance, the transmission error has a significant impact on gear vibration and noise, and the qualified contact pattern can avoid edge contact and premature failure [3,4,5]. Using a meshing simulation, the transmission errors and contact pattern of the gear pair can be determined.

Currently, the new optimization method based on the loaded tooth contact analysis (LTCA) for improving strength and contact performance of bevel gears has been a great concern to scholars worldwide. Simon analyzed the influences of machine settings and misalignments on tooth contact pressure and loaded transmission error, and then presented a method to improve gear meshing performance by reducing the loaded transmission error and the maximum tooth contact pressure of hypoid gears [6,7,8,9,10]. Antoni et al. studied the ease-off optimizing method for loaded transmission error and contact surface of bevel gears, which reduced the risk of vibration and noise [11,12,13]. Kolivand et al. [14,15] developed a novel formulation for LTCA using the ease-off topography. The study defined the surface of action and roll angle surface to simplify the task of locating the instantaneous contact curves. Yanming et al. proposed an ease-off flank modification method to solve the problem of the design and fabrication of high contact ratio spiral bevel gears with a seventh-order transmission error [16]. Nie et al. studied a new method of tooth surface topography modification for improving the meshing performance and correct contact area of spiral bevel gears [17].

Furthermore, a lot of research efforts have already gone into the analysis and optimization of individual gear pairs by experimental testing. De Vaujany et al. discussed the theoretical analysis and experimental research on load transmission error and actual coincidence of spiral bevel gears [18]. Kazumasa et al. investigated the tooth contact pattern and transmission errors of large-sized spiral bevel gears, and these analyzed results were compared with experimental ones [19]. Zhuo et al. investigated the contact characteristics of hypoid gears under quasi-static condition, and proposed a method for the global optimization of the tooth contact pattern and transmission error of spiral bevel and hypoid gears [20,21]. Cao et al. discussed a new method for gear tooth contact analysis, which concentrates on solving two considerable disadvantages of the generalized algorithm, which were verified with the gear-rolling test [22]. Liang et al. studied the influence of different grinding parameters on the surface topography and compared the simulation results with experiment results [23].

It can be seen that the research of loaded contact analysis and experimental testing of spiral bevel gear has not been sufficient so far. There are few comprehensive studies on its high loading comprehensive contact performance. In this work, an accurate finite element model of spiral bevel gear transmission system was established, which was based on the generating principle of spiral bevel gear and meshing theory. This model considered the influence of gear and shaft deformation, installation position, and bearing load, and the contact characteristics under quasi-static loading and high loading condition were analyzed through this model. In addition, a contact performance test bench of spiral bevel gear pair was developed in this paper, and the bending stress, contact pattern and transmission error test results were obtained for validating FEM results and the reliability of proposed model.

2. Theoretical Tooth Surface Modeling

2.1. Tooth Surface Modeling of the Gear

The method described in this paper is the generating method, in which both sides of the gear teeth are processed simultaneously with a double-sided cutter disc containing an inner cutter and an outer cutter. When machining large wheels, the cutter head rotates around its own axis while revolving around the axis of the production wheel. The position relationship between the large wheel blank and the cutter head is shown in Figure 1.

${\beta }_2$ is the coordinate system of the gear machine tool; $P_2$ is the unit vector of the axis direction of the gear; ${\alpha }_{02}$ is the pressure angle of the gear cutter head; ${\beta }_{02}$ is the coordinate system of the gear cutter head; $S_2$ is the radial cutter position; $q_2$ is the angular cutter position; $E_{02}$ is the vertical wheel position; $X_{02}$ is the axial wheel position; $X_{B2}$ is the bed position;${\delta }_{M2}$ is the root angle of machine tool; $r_0$ is the nominal radius of cutter head; $W_2$ is the angular cutter head radius; $r_{e2}$ is cutter tip arc radius; $O_e$ is the tool tip arc center; $M_0$ is the apex of the tool tip; ${\theta }_2$ is the angle between the section of cutter disk axis passing $\overline{O_0{M}_0}$ and ${\mathit{i}}_2$; $r_{02}$ is the radius of the tool tip.

Let $\left|\overline{M{M}_0}\right|$=$s_2$; then, in ${\beta }_{02}$, the radius vector ${\mathit{r}}_{C2}$ and normal vector ${\mathit{n}}_{C2}$ of M points on the cutting surface can be represented as follows:

$${\mathit{r}}_{c2}\left(s_2,{\theta }_2\right)=\left[\begin{array}{c}\left(r_{02}-s_{2}sin{\alpha }_{02}\right)cos{\theta }_2\\ \left(r_{02}-s_{2}sin{\alpha }_{02}\right)sin{\theta }_2\\ -s_{2}cos{\alpha }_{02}\end{array}\right]$$

(1)

$${\mathit{n}}_{c2}\left({\theta }_2\right)=\left[\begin{array}{c}cos{\alpha }_{02}cos{\theta }_2\\ cos{\alpha }_{02}sin{\theta }_2\\ -sin{\alpha }_{02}\end{array}\right]$$

(2)

${\theta }_{02}$ is the angle between $\overline{O_{e}M}$ and the axis of the cutter head, and the radius vectors ${\mathit{r}}_{e2}$ and normal vectors ${\mathit{n}}_{e2}$ of M on the cutter tip arc can be represented as

$${\mathit{r}}_{e2}\left({\theta }_{02},{\theta }_2\right)=\left[\begin{array}{c}\left(r_{Oe2}\mp r_{e2}sin{\theta }_{02}\right)cos{\theta }_2\\ \left(r_{Oe2}\mp r_{e2}sin{\theta }_{02}\right)sin{\theta }_2\\ r_{e2}\left(cos{\theta }_{02}-1\right)\end{array}\right]$$

(3)

$${\mathit{n}}_{e2}\left({\theta }_{02},{\theta }_2\right)=\left[\begin{array}{c}sin{\theta }_{02}cos{\theta }_2\\ sin{\theta }_{02}sin{\theta }_2\\ \mp cos{\theta }_{02}\end{array}\right]$$

(4)

The cutter head rotates $∆q_2$ around k, and the relationship between ${\beta }_{02}$ and ${\beta }_2$ can be obtained:

$${\mathit{M}}_{{\beta }_{02}}(∆q_2)=\left[\begin{array}{cccc}1& 0& 0& S_2\cos \left(q_2+∆q_2\right)\\ 0& 1& 0& S_{2}sin\left(q_2+∆q_2\right)\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(5)

Gear coordinate system ${\beta }_g$ moves $E_{02}$ along $\mathit{j}$:

$${\mathit{M}}_{E_{02}}=\left[\begin{array}{cccc}1& 0& 0& 1\\ 0& 1& 0& E_{02}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(6)

Move $X_{B2}$ along k:

$${\mathit{M}}_{X_{B2}}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& X_{B2}\\ 0& 0& 0& 1\end{array}\right]$$

(7)

Then, rotate (90°-${\delta }_{M2}$) around $\mathit{j}$,

$${\mathit{M}}_{{\delta }_{M2}}=\left[\begin{array}{cc}\begin{array}{cc}sin{\delta }_{M2}& 0\\ 0& 1\end{array}& \begin{array}{cc}cos{\delta }_{M2}& 0\\ 0& 0\end{array}\\ \begin{array}{cc}-cos{\delta }_{M2}& 0\\ 0& 0\end{array}& \begin{array}{cc}sin{\delta }_{M2}& 0\\ 0& 1\end{array}\end{array}\right]$$

(8)

Move $X_2$ along ${\mathit{P}}_2$:

$${\mathit{M}}_{X_2}=\left[\begin{array}{cc}\begin{array}{cc}1& 0\\ 0& 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}1& X_2\\ 0& 1\end{array}\end{array}\right]$$

(9)

The cutter head rotates $∆q_2$ around $\mathit{k}$ and the gear rotates ${\psi }_2=i_{02}∆q_2$ around ${\mathit{P}}_2$:

$${\mathit{M}}_{{\psi }_2}\left(∆q_2\right)=\left[\begin{array}{cccc}\begin{array}{c}cos\left({\psi }_2\right)\\ sin\left({\psi }_2\right)\\ 0\\ 0\end{array}& \begin{array}{c}-sin\left({\psi }_2\right)\\ cos\left({\psi }_2\right)\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 1\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\end{array}\right]$$

(10)

In summary, the homogeneous transformation matrix of the gear coordinate system ${\beta }_g$ and the machine tool coordinate system ${\beta }_2$ obtain:

$${\mathit{M}}_{{\beta }_g}={\mathit{M}}_{E_{02}}{\mathit{M}}_{X_{B2}}{\mathit{M}}_{{\delta }_{\mathit{M}2}}{\mathit{M}}_{X_2}{\mathit{M}}_{{\psi }_2}$$

(11)

The homogeneous transformation matrix of the cutter head coordinate system ${\beta }_{02}$ and the gear coordinate system ${\beta }_g$ can be calculated as:

$${\mathit{M}}_{02}\left(∆q_2\right)={\left({\mathit{M}}_{{\beta }_g}\right)}^{-1}{\mathit{M}}_{{\beta }_{02}}$$

(12)

In the gear coordinate system ${\beta }_g$, the radius vector ${\mathit{r}}_2$ of the conjugate tooth surface of the gear can be represented as:

$${\mathit{r}}_2(s_2,{\theta }_2,∆q_2)={\mathit{M}}_{02}{\mathit{r}}_{c2}$$

(13)

The unit normal vector ${\mathit{n}}_2$ of the conjugate tooth surface of the gear can be represented as:

$${\mathit{n}}_2({\theta }_2,∆q_2)={\mathit{M}}_{02}{\mathit{n}}_{c2}$$

(14)

The relative velocity at any point of the conjugate tooth surface of the gear can be calculated as:

$${\mathit{v}}_g\left(s_2,{\theta }_2,∆q_2\right)=\frac{\partial {\mathit{r}}_2}{\partial {\psi }_2}$$

(15)

The expressions of ${\mathit{v}}_g$ and ${\mathit{n}}_2$ are substituted for the meshing equation, and the following can be obtained:

$${\mathit{v}}_g·{\mathit{n}}_2=0$$

(16)

From Equation (16), it can be concluded that $s_2$ is a function with ${\theta }_2$ and $∆q_2$ as parameters. Plugging $s_2=f\left({\theta }_2,∆q_2\right)$ into Equation (13), the equation ${\mathit{r}}_2\left({\theta }_2,∆q_2\right)$ of the gear surface and transition surface is obtained. When $∆q_2$ changes, a series of contact lines can be calculated, and constituting the tooth surface of the gear, as shown in Figure 2.

2.2. Tooth Surface Modeling of the Pinion

The pinion is machined by tool tilting mechanism and one-sided method, and modeling of the pinion is more complicated than that of the gear. The forming process of the pinion is shown in Figure 3.

${\beta }_1$ is the coordinate system of the pinion machine tool; ${\mathit{P}}_1$ is the unit vector of the axis direction of the gear; ${\alpha }_{01}$ is the pressure angle of the pinion cutter head; ${\beta }_{01}$ is the coordinate system of the gear cutter head; $\mathit{b}$ is the projection of ${\mathit{k}}_{01}$ on plane ${\mathit{i}}_1-{\mathit{j}}_1$; $S_1$ is the radial cutter position; $q_1$ is the angular cutter position; $E_{01}$ is the vertical wheel position; $X_{01}$ is the axial wheel position; $X_{B1}$ is the bed position;${\delta }_{M1}$ is the root angle of machine tool; $j$ is the turning angle of the cutter; $i$ is the inclination angle of the cutter; $r_{01}$ is the radius of the tool tip; $r_{e1}$ is the cutter tip arc radius;${\theta }_1$ is the angle between the section of cutter disk axis passing $\overline{O_0{M}_0}$ and ${\mathit{i}}_1$. Let $\left|\overline{M{M}_0}\right|$=$s_1$; then, in ${\sigma }_{01}$, the radius vector ${\mathit{r}}_{C1}$ and normal vector ${\mathit{n}}_{C1}$ of M points on the cutting surface can be represented as follows:

$${\mathit{r}}_{c1}\left(s_1,{\theta }_1\right)=\left[\begin{array}{c}\left(r_{01}-s_{1}sin{\alpha }_{01}\right)cos{\theta }_1\\ \left(r_{01}-s_{1}sin{\alpha }_{01}\right)sin{\theta }_1\\ -s_{1}cos{\alpha }_{01}\end{array}\right]$$

(17)

$${\mathit{n}}_{c1}\left({\theta }_1\right)=\left[\begin{array}{c}cos{\alpha }_{01}cos{\theta }_1\\ cos{\alpha }_{01}sin{\theta }_1\\ -sin{\alpha }_{01}\end{array}\right]$$

(18)

${\theta }_{01}$ is the angle between $\overline{O_{e}M}$ and the axis of the cutter head, and the radius vectors ${\mathit{r}}_{e1}$ and normal vectors ${\mathit{n}}_{e1}$ of M on the cutter tip arc can be represented as:

$${\mathit{r}}_{e1}\left({\theta }_{01},{\theta }_1\right)=\left[\begin{array}{c}\left(r_{Oe1}\mp r_{e1}sin{\theta }_{01}\right)cos{\theta }_1\\ \left(r_{Oe1}\mp r_{e1}sin{\theta }_{01}\right)sin{\theta }_1\\ r_{e1}\left(1-cos{\theta }_{01}\right)\end{array}\right]$$

(19)

$${\mathit{n}}_{e1}\left({\theta }_{01},{\theta }_1\right)=\left[\begin{array}{c}sin{\theta }_{01}cos{\theta }_1\\ sin{\theta }_{01}sin{\theta }_1\\ \mp cos{\theta }_{01}\end{array}\right]$$

(20)

The cutter head rotates $∆q_1$ around k:

$${\mathit{M}}_{∆q_1}(∆q_1)=\left[\begin{array}{cccc}1& 0& 0& S_2\cos \left(q_2+∆q_2\right)\\ 0& 1& 0& S_{2}sin\left(q_2+∆q_2\right)\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(21)

${\beta }_{01}$ rotates ($q_1-j$) around ${\mathit{k}}_{01}$:

$${\mathit{M}}_j=\left[\begin{array}{cccc}\begin{array}{c}cos\left(q_1-j\right)\\ sin\left(q_1-j\right)\\ 0\\ 0\end{array}& \begin{array}{c}-sin\left(q_1-j\right)\\ cos\left(q_1-j\right)\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 1\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\end{array}\right]$$

(22)

Then, $i\mathrm{rotates}$ around ${\mathit{i}}_{01}$:

$${\mathit{M}}_i=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& cosi& sini& 0\\ 0& -sini& -cosi& 0\\ 0& 0& 0& 1\end{array}\right]$$

(23)

In conclusion, the relationship between ${\beta }_{01}$ and ${\beta }_1$ can be represented as:

$${\mathit{M}}_{{\beta }_{01}}\left(∆q_1\right)={\mathit{M}}_{∆q_1}{\mathit{M}}_j{\mathit{M}}_i$$

(24)

${\beta }_P$ moves $E_{01}$ along $\mathit{j}$:

$${\mathit{M}}_{E_{01}}=[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& E_{01}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}]$$

(25)

Then, $-X_{B1}$ moves along $\mathit{k}$:

$${\mathit{M}}_{X_{B1}}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& -X_{B1}\\ 0& 0& 0& 1\end{array}\right]$$

(26)

Rotating (90°-${\delta }_{M1}$) around $\mathit{j}$:

$${\mathit{M}}_{{\delta }_{M1}}=\left[\begin{array}{cc}\begin{array}{cc}sin{\delta }_{M1}& 0\\ 0& 1\end{array}& \begin{array}{cc}-cos{\delta }_{M1}& 0\\ 0& 0\end{array}\\ \begin{array}{cc}cos{\delta }_{M1}& 0\\ 0& 0\end{array}& \begin{array}{cc}sin{\delta }_{M1}& 0\\ 0& 1\end{array}\end{array}\right]$$

(27)

Moving $-X_1$ along ${\mathit{P}}_1$:

$${\mathit{M}}_{X_1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& -X_1\\ 0& 0& 0& 1\end{array}\right]$$

(28)

The cutter head rotates $∆q_1$ around $\mathit{k}$ and the pinion rotates ${\psi }_1=i_{01}∆q_1$ around ${\mathit{P}}_1$:

$${\mathit{M}}_{{\psi }_1}\left(∆q_1\right)=\left[\begin{array}{cccc}\begin{array}{c}cos\left({\psi }_1\right)\\ sin\left({\psi }_1\right)\\ 0\\ 0\end{array}& \begin{array}{c}-sin\left({\psi }_1\right)\\ cos\left({\psi }_1\right)\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 1\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\end{array}\right]$$

(29)

The homogeneous transformation matrix of the pinion coordinate system ${\beta }_P$ and the machine tool coordinate system ${\beta }_1$ is obtained:

$${\mathit{M}}_{{\beta }_p}={\mathit{M}}_{E_{01}}{\mathit{M}}_{X_{B1}}{\mathit{M}}_{{\delta }_{\mathit{M}1}}{\mathit{M}}_{X_1}{\mathit{M}}_{{\psi }_1}$$

(30)

The homogeneous transformation matrix of the cutter head coordinate system ${\beta }_{01}$ and the pinion coordinate system ${\beta }_P$ can be calculated as:

$${\mathit{M}}_{01}\left(∆q_1\right)={\left({\mathit{M}}_{{\beta }_p}\right)}^{-1}{\mathit{M}}_{{\beta }_{01}}$$

(31)

In the pinion coordinate system ${\beta }_P$, the radius vector ${\mathit{r}}_1$ of the conjugate tooth surface can be represented as:

$${\mathit{r}}_1(s_1,{\theta }_1,∆q_1)={\mathit{M}}_{01}{\mathit{r}}_{c1}$$

(32)

The unit normal vector ${\mathit{n}}_1$ of the conjugate tooth surface of the pinion can be represented as:

$${\mathit{n}}_1({\theta }_1,∆q_1)={\mathit{M}}_{01}{\mathit{n}}_{c1}$$

(33)

The relative velocity at any point of the conjugate tooth surface of the pinion can be calculated:

$${\mathit{v}}_p\left(s_1,{\theta }_1,∆q_1\right)=\frac{\partial {\mathit{r}}_1}{\partial {\psi }_1}$$

(34)

The expressions of ${\mathit{v}}_p$ and ${\mathit{n}}_1$ are substituted for the meshing equation, and the following can be obtained:

$${\mathit{v}}_p·{\mathit{n}}_1=0$$

(35)

From Equation (35), it can be concluded that $s_1$ is a function with ${\theta }_1$ and $∆q_1$ as parameters. Plugging $s_1=f\left({\theta }_1,∆q_1\right)$ into Equation (32), the equation ${\mathit{r}}_1\left({\theta }_1,∆q_1\right)$ of the pinion surface and transition surface is obtained. When $∆q_1$ changes, a series of contact lines can be calculated, and constitute the tooth surface of the pinion, as shown in Figure 4.

2.3. Meshing Contact Modeling

The mathematic models of the gear and the pinion are established in the gear coordinate system ${\beta }_g$ and the pinion coordinate system ${\beta }_p,$ respectively. In order to make the position relationship of spiral bevel gear pair consistent with the theoretical design, transforming the pinion coordinate system ${\beta }_p$ is into the gear coordinate system ${\beta }_g$. The gear rotates around ${\mathit{p}}_2$ axis and the pinion rotates around ${\mathit{p}}_1$ axis, when ${\mathit{r}}_2,{\mathit{n}}_2$ rotates ${\mathit{\phi }}_2$ around ${\mathit{p}}_2$, and ${\mathit{r}}_1,{\mathit{n}}_1$ rotates ${\mathit{\phi }}_1$ around ${\mathit{p}}_1$, the conjugate contact between M point on the gear tooth surface and M point on the pinion tooth surface is established. The equation ${\mathit{R}}_2$ and the normal vector ${\mathit{N}}_2$ of tooth surface of the gear at the contact point can be represented as:

$${\mathit{R}}_2=({\mathit{r}}_2·{\mathit{p}}_2){\mathit{p}}_2+sin{\phi }_2\left({\mathit{p}}_2\times {\mathit{r}}_2\right)+cos{\phi }_2\left({\mathit{p}}_2\times {\mathit{r}}_2\right)\times {\mathit{p}}_2$$

(36)

$${\mathit{N}}_2=({\mathit{n}}_2·{\mathit{p}}_2){\mathit{p}}_2+sin{\phi }_2\left({\mathit{p}}_2\times {\mathit{n}}_2\right)+cos{\phi }_2\left({\mathit{p}}_2\times {\mathit{n}}_2\right)\times {\mathit{p}}_2$$

(37)

The equation ${\mathit{R}}_1$ and the normal vector ${\mathit{N}}_1$ of the gear at the contact point can be represented as:

$${\mathit{R}}_1=({\mathit{r}}_1·{\mathit{p}}_1){\mathit{p}}_1+sin{\phi }_1\left({\mathit{p}}_1\times {\mathit{r}}_1\right)+cos{\phi }_1\left({\mathit{p}}_1\times {\mathit{r}}_1\right)\times {\mathit{p}}_1$$

(38)

$${\mathit{N}}_1=({\mathit{n}}_1·{\mathit{p}}_1){\mathit{p}}_1+sin{\phi }_1\left({\mathit{p}}_1\times {\mathit{n}}_1\right)+cos{\phi }_1\left({\mathit{p}}_1\times {\mathit{n}}_1\right)\times {\mathit{p}}_1$$

(39)

According to the meshing principle of gears, ${\mathit{R}}_2,{\mathit{R}}_1,{\mathit{N}}_2,$ and ${\mathit{N}}_1$ should satisfy the equations as follows:

$$\{\begin{array}{c}{\mathit{R}}_2={\mathit{R}}_1-\overline{O_1{O}_2}\\ {\mathit{N}}_2={\mathit{N}}_1\end{array}$$

(40)

There are six unknowns ${\theta }_1,{\theta }_2,∆q_1,∆q_2,{\phi }_1,$ and ${\phi }_2$ in the above equations; however, there are only five independent equations. When the pinion rotation angle ${\phi }_1$ is given as input according to the step size, the equations can be solved. The pinion rotation angle at each moment corresponds to a meshing point, and all meshing points from entry to exit form a contact trace. The following Figure 5 shows the contact point position of the two teeth surfaces (after the gear rotating angle ${\phi }_2$ and the pinion rotating ${\phi }_1$).

3. Solid Modeling and Finite Element Model

3.1. Examples of Calculation

In order to verify the correctness of the previous mathematical model of the gear pair and to analyze the loaded contact characteristics of the tooth surface, a pair of Grison spiral bevel gears which is used in the axle of a domestic heavy-duty vehicle was selected as the test specimen. The parameters are listed in

Table 1. Basic tooth geometric parameters.

Geometrical Parameter	Description
	Pinion	Gear
Number of teeth	15	46
Hand of spiral	Left	Right
Modification coefficient	0.35	−0.35
Pitch diameter (mm)	90	276
Tip diameter (mm)	117.48	423.28
Module (mm)	6
Shaft angle (°)	90
Pressure angle(°)	20
Tooth width(mm)	44
Spiral angle (°)	35

.

3.2. Constructing Solid Geometry Model

Based on the mathematical model established above, the calculated lattice data of tooth surface, transition surface, and root cone surface are imported into the software Pro/Engineer, and the discrete points are fitted with spline curve, which are shown in Figure 6 below.

Based on the spline curve, the gear surface (including concave, convex tooth surface, fillet surface and root cone surface) is constructed by some software function, such as variable section sweep, extend, boundary blend, trim, and so on. Then, the solid model of the big gear and the small gear can be obtained by cutting the surface of the gear blank model with the gear surface. The solid model of the spiral bevel gear pair can be obtained by virtual assembly according to the theoretical installation position of the gear pair, which is shown in Figure 7.

The gear modeling method is designed according to the gear forming method in the actual gear cutting process, and the consistency between the theoretical model and the actual gear is ensured as much as possible in the mathematical modeling.

3.3. Transmission System Model

The assembly model of spiral bevel gear was imported into ANSYS software for finite element analysis. Under actual conditions, the gear pair will not be an isolated system. The whole transmission system includes supporting shafts, bearings, Gearbox, and so on. However, if all parts of the transmission system are imported into software, the solving time will be greatly increased, and the results may not converge. Therefore, it is necessary to simplify the model properly.

As shown in Figure 8, other structural components are omitted, such as bearings and gearboxes, which have little influence on the analysis results, and key parts are retained, such as gear teeth and rotating shafts. In order to simplify the analysis, the gear and its rotating shaft are consolidated together. The gear adopts two-end span support structure, while the pinion adopts the common cantilever support structure. The bearing support position of the pinion spindle is not only installation and positioning boundary, but also input end of torque ${\mathit{T}}_1$, and the bearing support near the large wheel is load ${\mathit{T}}_2$ input end.

3.4. Material Property and Mesh Generation

The material of the spiral bevel gear pair is 20CrMnTiA, which was heat treated by carburization, quench, and low temperature tempering. The corresponding material properties are shown in

Table 2. Material properties of 20CrMnTiA.

Material	$\mathbf{E}/\left(\mathbf{G}\mathbf{P}\mathbf{a}\right)$	v	$\mathsf{\rho }/(\mathbf{k}\mathbf{g}/{\mathbf{m}}^3)$	${\mathit{\sigma }}_{\mathit{H}\mathit{l}\mathit{i}\mathit{m}}/\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	${\mathit{\sigma }}_{\mathit{F}\mathit{l}\mathit{i}\mathit{m}}/\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	${\mathit{R}}_{0.2}/\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$
20CrMnTiA	207	0.25	7800	1475	415	850

. Due to heavy load, the material of the rotating shaft adopts high strength carburized steel. In order to get better contact analysis results, the whole model uses eight-node linear hexahedron element for contact analysis. Meanwhile, in order to improve the calculation accuracy of tooth surface contact analysis, the mesh refinement of contact surface of gear pair is carried out, listed in Figure 9. We selected the Von Mises material model, which is suitable for metallic materials finite element analysis.

4. Simulation Results

4.1. Tooth Surface Contact Pattern

Different torques (500 Nm, 1500 Nm and 4500 Nm) are set to simulate light, medium, and heavy loads, respectively, the variation of contact pattern on the convex face of the gear is shown in Figure 10.

We can find in Figure 10 above that, when load T is 500 Nm, the instantaneous contact pattern is relatively small. Compared Figure 10b,c with Figure 10a, the directional angle γ becomes smaller, the contact line becomes longer, and the instantaneous contact pattern becomes larger with the load T increasing. Therefore, we can confirm that the load is an important factor affecting the contact area and the contact trace of the bevel gear tooth surface. In addition, the contact zone of tooth surface extends to the large end with the increase of load T, but not to the small end. When the torque T increases to 4500 Nm, the edge contact occurs on the tooth surface of the gear which will induce vibration, shock, and be harmful to transmission stability.

4.2. Tooth Surface Contact Stress

When load T = 1500 Nm, T = 4500 Nm, and T = 7500 Nm, the variation of contact stress on the convex face of the gear is shown in Figure 11.

As shown in Figure 11, the above results indicates that, with the increase of load T, the maximum contact stress ${\sigma }_H$ of the gear convex face becomes larger, and the worst contact zone appears in the small end teeth top of the gear. When load increases from 1500 Nm to 7500 Nm, the comparisons between the simulation results and the calculated results of empirical formula are listed in

Table 3. Results comparison of FEM and empirical formula of maximum contact stresses.

Load T/Nm	$\mathbf{Calculated}\mathbf{Results}{\mathit{\sigma }}_{\mathit{H}}/\mathbf{MPa}$	$\mathbf{Simulation}\mathbf{Results}{\mathit{\sigma }}_{\mathit{H}}/\mathbf{MPa}$	Values of Differences Δ%
1500	632.79	573.83	9.32
3000	815.48	758.25	7.02
4500	1045.32	914.49	12.51
6000	1324.21	1138.1	14.05
7500	1569.74	1331	15.21

below. The empirical formula is according to ISO 10300 and ISO 6336, which is shown in Equations (41) and (42) [24,25]. The calculation is based on the contact (Hertzian) stress, in which the load is distributed along the lines of contact. Calculations are to be carried out for pinion and wheel together, and some parameters for contact stress calculation were determined by engineering experience, processing technology, and actual operation conditions, such as gear geometry, manufacturing accuracy, rigidity of flange, bearing and bearing seat, and torque transmitted. We can find that, with the increase of loading torque, the change trend of FEM and calculated results are similar. In addition, the results of empirical formulas are about 10–15% larger than analyzed ones. The authors believe that the FEM accuracy, constraint condition setting, and the parameters selection of theoretical formula cause the difference value:

$${\mathit{\sigma }}_H={\mathit{\sigma }}_{H0}\sqrt{{\mathit{K}}_A{\mathit{K}}_V{\mathit{K}}_{H\beta }{\mathit{K}}_{H\alpha }}$$

(41)

$${\mathit{\sigma }}_{H0}=\sqrt{\frac{{\mathit{F}}_n}{{\mathit{l}}_{bm}{\mathit{\rho }}_{rel}}}{\mathit{Z}}_{M-B}{\mathit{Z}}_{LS}{\mathit{Z}}_E{\mathit{Z}}_K$$

(42)

$K_A$ is the application factor, ${\mathit{K}}_V$ is the dynamic factor, ${\mathit{K}}_{H\beta }$ is the face load factor for contact stress, ${\mathit{K}}_{H\alpha }$ is the transverse load factor for contact stress, ${\mathit{\sigma }}_{H0}$ is the nominal value of the contact stress, ${\mathit{F}}_n$ is the nominal normal force of the virtual cylindrical gear at mean point, ${\mathit{l}}_{bm}$ is the length of contact line in the middle of the zone of action, ${\mathit{\rho }}_{rel}$ is the radius of relative curvature vertical to the contact line, ${\mathit{Z}}_{M-B}$ is the mid-zone factor, ${\mathit{Z}}_{LS}$ is the load sharing factor, ${\mathit{Z}}_E$ is the elasticity factor, and ${\mathit{Z}}_K$ is the bevel gear factor.

4.3. Tooth Root Bending Stress

When load T = 500 Nm, T = 1500 Nm and T = 4500 Nm, the maximum root bending stress variation of the gear and the pinion is shown in Figure 12.

As shown in Figure 12, the maximum root bending stress appears near the transition line between tooth surface and tooth root. When the load T increases, the root bending stresses also increase, and the maximum root bending stress ${\sigma }_{F2}$ of the gear is less than the maximum root bending stress ${\sigma }_{F1}$ of the pinion. When torque increases from 1500 Nm to 7500 Nm, the maximum root bending stress ${\sigma }_F$ comparisons between the FEM results and the calculated results of empirical formula are listed in

Table 4. Results comparison of FEM and empirical formula of bending stresses.

Load T/Nm	Gear		Pinion
	FEM Results ${\mathit{\sigma }}_{\mathit{F}2}/\mathbf{MPa}$	Calculated Results ${\mathit{\sigma }}_{\mathit{F}2}/\mathbf{MPa}$	Values of Differencs Δ%	FEM Results ${\mathit{\sigma }}_{\mathit{F}1}/\mathbf{MPa}$	Calculated Results ${\mathit{\sigma }}_{\mathit{F}1}/\mathbf{MPa}$	Values of Differences Δ%
1500	296	343	13.71	338	389	13.11
3000	456	529	13.79	567	643	11.82
4500	621	710	12.54	791	883	10.42
6000	805	890	9.55	925	1039	13.86
7500	912	1067	12.52	1183	1352	12.5

below. The empirical formula is according to ISO 10300 and ISO 6336, which is shown in Equations (43) and (44). The calculation is based on the maximum bending stress at the tooth root. It is determined separately for pinion and wheel. Some parameters for root bending stress calculation were determined by engineering experience and actual operation conditions, such as geometric similarity, operation, and manufacturing of actual gears:

$${\sigma }_F={\sigma }_{F0}{K}_{\mathrm{A}}{K}_{\mathrm{V}}{K}_{F\beta }{K}_{F\mathsf{\alpha }}$$

(43)

$${\sigma }_{F0}=\sqrt{\frac{F_{\mathrm{vmt}}}{{\mathit{b}}_v{\mathit{m}}_{mn}}}{\mathit{Y}}_{Fa}{\mathit{Y}}_{Sa}{\mathit{Y}}_{\epsilon }{\mathit{Y}}_{BS}{\mathit{Y}}_{LS}$$

(44)

${\mathit{K}}_{F\beta }$ is the face load factor for bending stress, $K_{F\mathsf{\alpha }}$ is the transverse load factor for bending stress, ${\sigma }_{F0}$ is the nominal tooth root stress, $F_{\mathrm{vmt}}$ is the nominal tangential force of the virtual cylindrical gear, ${\mathit{b}}_v$ is the face width of virtual cylindrical gear,${\mathit{m}}_{mn}$ is mean normal module, ${\mathit{Y}}_{Fa}$ is the tooth form factor, ${\mathit{Y}}_{Sa}$ is the stress correction factor, ${\mathit{Y}}_{\epsilon }$ is the contact ration factor, ${\mathit{Y}}_{BS}$ is the bevel spiral angle factor, and ${\mathit{Y}}_{LS}$ is the load sharing factor.

As shown in Table 4, the results show that the results of finite element analysis and empirical formula are similar, and the empirical value is about 10–14% larger than that of finite element analysis. The authors believe that one reason is that the load T increases the deformation of the bevel gears, which increases the actual contact ratio of the bevel gears, but the empirical value does not take this into account. Another reason is that the finite element model constraint condition setting and the parameters selection of theoretical formula will affect the difference value.

In order to observe the distribution of the maximum tooth root bending stress in a meshing period, the FEM analysis of the root bending stress is carried out by rotating the gear step by step. When load increases from 1500 Nm to 7500 Nm, the variation of the maximum value of root bending stresses of the pinion and the gear is shown in Figure 13.

As shown in Figure 13, the horizontal axis represents single tooth meshing period of the gear and pinion due to instantaneous contact point from the small end to the large end, or vice versa. Longitudinal axis is the maximum bending stress distribution along the root cone line. We can find that, when the load T is small, the mean bending stresses of the small end teeth root of gear and pinion are greater than that of the big end teeth root during the same meshing period because the modulus of the small end are smaller, which causes the bending strength of the small end to be relatively weak. When the load T increases gradually, the mean bending stresses of the small end teeth root of gear and the pinion are still greater than that of the big end teeth root during the same meshing period; however, the reasons are quite different. From previous research of the contact pattern, we can know that the moving direction of the contact pattern on the gear teeth surface is from the small end teeth top to the big end teeth root with the torque T increasing, making the contact pattern closer to the big end of the gear. Thus, the bending stresses of the small end teeth root of the gear are slightly larger than that of the big end teeth root. In contrast, the moving direction of the contact pattern on the pinion teeth surface is from the big end teeth root to the small end teeth top with the torque T increasing, making the contact pattern closer to the small end of the pinion. Thus, the bending stress of the small end tooth root of the pinion is much larger than that of the big end tooth root. Meanwhile, because of the distance of the tooth direction of the pinion is larger, the distribution of the root bending stress of the pinion is more concentrated in the tooth direction.

4.4. Transmission Error

We simulated transmission error by KISSsoft software, which is widely used and reliable in the gear field. A precise spiral bevel gear transmission system model was established before simulation calculation, and KISSsoft can calculate the difference value between theoretical and actual rotation angle during the meshing process, which considers the influence of bearings and the deformation of shafts and gears. According to the actual load condition, different moments (500 Nm, 1500 Nm and 4500 Nm) are set to simulate light, medium, and heavy loads, respectively.

The transmission error curves of this gear drives under different moment are shown in Figure 14. In addition, all transmission error curves are translated to zero position for better comparison. We can know that the transmission error change period is synchronized with the pitch angle of the pinion. Due to the increasing elastic deformation of the gear and the rotating shaft with the load T increases, the absolute value and peak to peak of the transmission error curves become larger. The transmission error curve under the heavy load (4500 Nm) fluctuates much more severely because of the increase of gear deformation and meshing impact.

5. Verification Test

5.1. Construction of Test Bench

As shown in Figure 15 below, a contact performance test bench of spiral bevel gear pair was developed. The test scheme adopts servomotor and high deceleration ratio gearbox combined drive mode to realize quasi-static loading. Meanwhile, different loading torques are achieved through hydraulic torque loading equipment.

5.2. Tooth Contact Pattern Test

When the loading torque T = 500 Nm, 1500 Nm, and 4500 Nm is applied separately, the real images of the contact pattern of the teeth surface are compared with the FEM simulation result as shown in Figure 16.

As shown in Figure 16 above, we can find that the test and simulation results are basically consistent, and directional angle and contact trace change trends are similar, which preliminarily verifies the reliability of previously generated finite element models. The theoretical contact ratio of this spiral bevel gear pair is 2.26, which is calculated by Equation (45) according to gear design and manufacture parameters [26]. The theoretical contact ratio is the maximum contact ratio of gear pair can be achieved with the increase of load under the condition of no edge contact. Compared the real images of Figure 16a with Figure 16b above, we can find that, when the load T increases from 500 Nm to 1500 Nm, the contact pattern becomes larger but still within the tooth surface boundary, so the actual contact ratio of condition (a) and (b) are smaller than theoretical contact ratio, and the value are 1.53 and 1.74, which are calculated by Equations (46) and (47) [27,28]. When load T is 4500 Nm, compare with 500 Nm and 1500 Nm, the contact pattern extend to the large end of the boundary and occur edge contact, so the actual contact ratio of condition (c) exceeds theoretical contact ratio, the result is 2.31.

$${\epsilon }_{\gamma }=\sqrt{{\epsilon }_{\alpha }{}^2+{\epsilon }_{\beta }{}^2}$$

(45)

${\epsilon }_{\gamma }$ is theoretical contact ratio, ${\epsilon }_{\alpha }$ is transverse contact ratio, and ${\epsilon }_{\beta }$ is overlap contact ratio.

$${\epsilon }_{\gamma }{}^{\prime }=\frac{{\epsilon }_{\alpha }{{}^{\prime }}^2}{\sqrt{{\epsilon }_{\alpha }{{}^{\prime }}^2-I_1{}^2}}$$

(46)

$$I_1=\sqrt[3]{\frac{C_P{T}_1}{C_{PD}{T}_{1D}}}\frac{{\epsilon }_{\alpha }{\epsilon }_{\beta }}{{\epsilon }_{\gamma }}$$

(47)

${\epsilon }_{\gamma }{}^{\prime }$ is actual contact ratio, ${\epsilon }_{\alpha }{}^{\prime }$ is revised transverse contact ratio, $C_P$ is load distribution ratio, $C_{PD}$ is load distribution ratio under full load, $T_1$ is actual load and $T_{1D}$ is full load.

In summary, the actual contact ratio increases with the load T increases, which is beneficial to increase the meshing stiffness and improve the smoothness of transmission. When actual contact ratio is too large, the gear pair will occur severe edge contact, vibration, shock and noise. On the one hand, we need design theoretical contact ratio properly according to gear parameters and actual operating conditions, avoiding actual contact ratio greater than theoretical contact ratio. On the other hand, when we need loading large load, the actual installation position should make the contact zone of the tooth surface close to the small end when the light load is applied, so as to fully exert the loading capacity of the tooth surface, and avoiding edge contact under heavy load as much as possible.

5.3. Tooth Root Bending Stress Test

As established in Figure 17, the numbers of root bending stress signal channels are limited, and tooth pitch of pinion is too small to install the strain gauges. Eight strain gauges are installed on the adjacent teeth of the gear before the test, and the positions are arranged from the small end to the big end of the gear. Because of excessive torque (T ≥ 5000 Nm) may cause damage to strain gauges, load T = 1500 Nm, T = 3000 Nm and T = 4500 Nm are applied to the gear respectively during the test.

The results comparison of test and FEM of root bending stresses is shown in Figure 18 below. We can find that the test results and the FEM results showed a trend similar, which can verify the effectiveness of the finite element model. The test bending stress of the small end teeth root are slightly larger than that of the large end teeth root, and the maximum root bending stress spot moves towards the large end with torque increases, which indicates the reliability of the previous discussion. The difference between the test results and the FEM results is that the root bending stress solved by simulation is the geometric average value of the stress in three directions, while the strain gauge measures the strain in one direction of the tooth height. Due to the influences of installation position error, measurement accuracy, and temperature factors, there is a deviation between the test results and the FEM values.

5.4. Transmission Error Test

As established in Figure 19 below, this test bench adopts high precision circular gratings for measuring transmission error of gear pairs under quasi-static loading, and the circular gratings were installed on the both ends of spiral bevel gear pair, which can measure the relative angular displacement. We calculated transmission error according to Equation (48):

$$TE=∆{\epsilon }_1-i_T∆{\epsilon }_2$$

(48)

where TE is transmission error, $i_T$ is theoretical transmission ratio between input and output gear, and $∆{\epsilon }_1$, $∆{\epsilon }_2$ is rotation angle of input and output gear measured by circular gratings during the same period separately.

The test results of transmission errors after translation are shown in Figure 20 below. We can know that the absolute value and peak to peak of test results are larger than that of FEM results, but the variation trend and fluctuation period of the two are similar.

We can further analyze the difference between test results with simulation results. We define the rotation speed of gear as ${\mathit{w}}_1={\mathit{\tau }}_1$, and the rotation speed of gear ${\mathit{w}}_2=i\left(t\right){\mathit{\tau }}_2$. $i\left(t\right)$ is instantaneous transmission ratio. Relative speed of gear and pinion in meshing position can be represented as:

$${\mathit{v}}_{12}={\mathit{\tau }}_1\times {\mathit{R}}_2-i\left(t\right){\mathit{\tau }}_2\times {\mathit{R}}_1$$

(49)

For the meshing equation ${\mathit{v}}_{12}·{\mathit{N}}_2=0$, we can get $i\left(t\right)$:

$$i\left(t\right)=\frac{\left({\mathit{\tau }}_2,{\mathit{R}}_2,{\mathit{N}}_2\right)}{\left({\mathit{\tau }}_1,{\mathit{R}}_1,{\mathit{N}}_1\right)}=\frac{\left({\mathit{\tau }}_2,{\mathit{r}}_2,{\mathit{n}}_2\right)}{\left({\mathit{\tau }}_1,{\mathit{r}}_1,{\mathit{n}}_1\right)}$$

(50)

We can know from above Equation (50), instantaneous transmission ratio $i\left(t\right)$ changes with contact point changes, so it is not necessarily equal to theoretical transmission ratio $i_T$. When load T is small, the fluctuation amplitude of $i\left(t\right)$ is not obvious, and $i\left(t\right)$ fluctuates around the theoretical value. When load T becomes larger, the fluctuation amplitude of $i\left(t\right)$ increases with the deformation and the contact ratio of gear pair increases. One reason for the difference is that the instantaneous transmission ratio of spiral bevel gear pair is changing, while the test results used the unchanging transmission ratio for calculation. Another reason is that the measurement accuracy and signal interference of circular grating affect the correctness of test results.

6. Conclusions

This paper presents an accurate mathematical model of tooth surface, fillet surface, and root cone surface of spiral bevel gear pair. Based on the calculated lattice data, an FEM model of spiral bevel gears transmission system has been constructed which considers the deformation of gear and shaft, bearing support, and solving time. For further study, we develop a test bench for contact performance of spiral bevel gear pairs and compare the formula calculated results, FEM results, and experiment results together. The formula calculated value of contact stress is about 10–15% larger than that of FEM analysis value. The real images of contact pattern are basically with simulation results, and the actual contact ratio of 500 Nm and 1500 Nm load are smaller than theoretical values, and are 1.53 and 1.74, respectively. When the load is 4500 Nm, the test gear has edge contact and the actual contact ratio exceeds the theoretical value, reaching 2.31. The formula calculated value of tooth root bending stress is about 10–14% larger than that of the FEM analysis value, and the test results and the FEM results show variation trends that are similar and have some differences too. The absolute value and peak to peak of transmission error test results are larger than that of FEM results, but the variation trend and fluctuation period of the two are basically consistently similar. The reason for the difference is that the instantaneous transmission ratio $i\left(t\right)$ is changing, while the test results used the theoretical transmission ratio $i_T$ for calculation.

In the future, the target will focus on FEM model accuracy by comparing with real parts using 3D scanning or cross-sections, optimize the measurement scheme, and finally reduce the difference between FEM results with experiment results. In the meantime, we will thoroughly study the global optimization design method of tooth contact characteristics and the correlation between quasi-static and dynamic operation condition.