Published: 30 September 2014

# Mathematical models and dynamic contact analysis of involute/noninvolute beveloid gears

He Zeyin1
Lin Tengjiao2
Chenjia Chen3
Feiyu Geng4
1, 2, 3, 4State Key Laboratory of Mechanical Transmission, Chongqing University, Chongqing, 400030, P. R. China
Corresponding Author:
Lin Tengjiao
Views 219

#### Abstract

This study investigates an approach for parametric modeling and dynamic contact analysis of involute/noninvolute beveloid gears. Firstly, the mathematical models of involute/noninvolute beveloid gear pairs are derived based on the theory of gearing and the generation mechanism. Then the parametric modeling programs of involute/noninvolute beveloid gears are developed to automatically generate exact model via a Matlab code. Subsequently, a numerical example of intersecting axes beveloid gears is presented to evaluate the dynamic stress distribution and dynamic transmission error. Finally, the dynamic contact characteristics of involute and noninvolute beveloid gears are calculated by three-dimensional dynamic contact finite element method, respectively. The results show that the noninvolute beveloid gear pairs can relieve the high dynamic stress and contact shock problem of intersecting axes beveloid gear pairs.

## 1. Introduction

In modern gearing application, spiral bevel gears, hypoid gears, worm gears and face gears can easily realize the non-parallel axes power transmission, and they represent good meshing characteristic at large shaft angles, close to 90°. However, it is not practical to use them when the shaft angle is less than 45°, due to design and manufacturing difficulties. Beveloid gears, by contrast, also known as conical gears with perfect meshing characteristics, are suitable for power transmission not only between parallel axes but between non-parallel axes [1, 2].

Recently, a number of studies have been performed on the beveloid gears. Brauer derived the parametric equations for a straight conical involute gear tooth surface [3], and these formulas were used to create a finite element model [4], then reported on a theoretical study of transmission errors in involute conical gear transmissions [5]. Liu, Ye and Chen established the mathematical models of variable thickness involute gears and the contact characteristic was also studied [6-8]. Chen solved the tooth face equations of helical involute beveloid gear, and the precise tooth surface was generated utilizing Matlab [9]. Wu proposed an approach for geometrical design and contact stress analysis of skew conical involute gear drives in approximate line contact [10]. Li investigated the tooth profile equation of noninvolute beveloid gears and also calculated the tooth profile errors and axial errors [11]. In spite of the above-mentioned studies focusing on meshing theory and simulation of involute or noninvolute beveloid gears, there are very few published works aimed to solve the engagement equation and tooth profile equation of noninvolute beveloid gears between intersecting axes by means of tooth face equation of the imaginary helical rack cutter. Furthermore, none also made a comparison between dynamic contact characteristics of involute beveloid gear pair and noninvolute beveloid gear pair. These gaps will be the emphasis of this paper.

## 2. Mathematical model of the involute beveloid gear pair

### 2.1. Transverse tooth profile equation of helical rack cutter

Considering the meshing relationship between imaginary rack cutter and gear blank, we can obtain the tooth profile equation of involute beveloid gear basing on the tooth profile equation of the rack cutter in transverse cross section. Therefore, the parameters of rack cutter in the transverse cross section should be calculated first.

The relative position relationship between rack cutter and gear blank is shown in Fig. 1. Here, ${r}^{\text{'}}$ means gear pitch radius, ${r}^{\text{'}}\omega$ is linear velocity of gear and the pitch plane of imaginary rack cutter is set to form an inclination angle $\delta$ with respect to the pitch plane of gear.

Fig. 1Relative position relationship

Fig. 2 describes the geometry of normal tooth profile of rack cutter. ${S}_{n}\left({X}_{n},{Y}_{n},{Z}_{n}\right)$ means fixed coordinate system linked with normal cross section of rack cutter. The geometry of normal tooth profile of rack cutter, which mainly contains straight edge $BC$ and tool fillet curve $AB$, is represented. The $BC$ and $AB$ generate the working tooth surface and fillet surfaces of involute beveloid gear, respectively.

Fig. 2Normal tooth profile of rack cutter

The equation of $BC$ in coordinate system ${S}_{n}$ can be described as:

1
${\mathbf{R}}_{n}^{z}\left(l\right)=\left[\begin{array}{c}lcos{\alpha }_{n}-{h}_{an}^{*}{m}_{n}\\ ±\left(-l\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{n}+{h}_{an}^{*}{m}_{n}\mathrm{t}\mathrm{a}\mathrm{n}{\alpha }_{n}+\frac{\pi {m}_{n}}{4}\right)\\ 0\\ 1\end{array}\right],$

where ${\mathbf{R}}_{n}^{z}$ is position vector of arbitrary point on $BC$ in coordinate system ${S}_{n}$; ${\alpha }_{n}$ is the normal pressure angle; ${m}_{n}$ is normal module; ${h}_{an}^{*}$ represents the addendum coefficient; $l$ is the distance between moving point on $BC$ and point $B$; the upper and lower signs of “±” are used to describe the left and right straight edge of rack cutter.

The equation of $AB$ in coordinate system ${S}_{n}$ can be described as:

2
${\mathbf{R}}_{n}^{j}\left(\theta \right)=\left[\begin{array}{c}-{h}_{an}^{*}{m}_{n}-\rho cos\theta +\rho sin{\alpha }_{n}\\ ±\left({h}_{an}^{*}{m}_{n}\mathrm{t}\mathrm{a}\mathrm{n}{\alpha }_{n}+\frac{\pi {m}_{n}}{4}+\rho \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{n}-\rho \mathrm{s}\mathrm{i}\mathrm{n}\theta \right)\\ 0\\ 1\end{array}\right],$

where ${\mathbf{R}}_{n}^{j}$ is position vector of arbitrary point on $AB$ in coordinate system ${S}_{n}$; $\theta$ denotes the central angle between moving point on $AB$ and point $A$; $\rho$ represents the radius of fillet.

Using the coordinate transformation from ${S}_{n}\left({X}_{n},{Y}_{n},{Z}_{n}\right)$ to ${S}_{t}\left({X}_{t},{Y}_{t},{Z}_{t}\right)$ which obtains from rotating ${S}_{n}\left({X}_{n},{Y}_{n},{Z}_{n}\right)$ around the ${X}_{n}$ axis with an angle $\beta$, the coordinate system ${S}_{d}\left({X}_{d},{Y}_{d},{Z}_{d}\right)$ will be got though rotating ${S}_{t}$ around the ${Y}_{t}$ axis with an angle $\delta$, and they are illustrated in Fig. 2. Combining Eqs. (1) and (2), we can gain the rack cutter equations of straight edge $BC$ and tool fillet curve $AB$ in the ${X}_{t}{O}_{t}{Z}_{t}$ plane of ${S}_{t}$ as follows:

3
$\left\{\begin{array}{l}{\mathbf{R}}_{t}^{z}\left(l\right)={\left[\begin{array}{cccc}{R}_{nx}^{z}& \frac{{R}_{ny}^{z}}{\mathrm{c}\mathrm{o}\mathrm{s}\beta }& 0& 1\end{array}\right]}^{T},\\ {\mathbf{R}}_{t}^{j}\left(\theta \right)={\left[\begin{array}{cccc}{R}_{nx}^{j}& \frac{{R}_{ny}^{j}}{\mathrm{c}\mathrm{o}\mathrm{s}\beta }& 0& 1\end{array}\right]}^{T}.\end{array}\right\$

Then the rack cutter equations of straight edge $BC$ and tool fillet curve $AB$ in the ${X}_{d}{O}_{d}{Z}_{d}$ plane of coordinate system ${S}_{d}$ are as follows:

4
$\left\{\begin{array}{l}{\mathbf{R}}_{d}^{z}\left(l\right)={\left[\begin{array}{cccc}\frac{{R}_{nx}^{z}}{\mathrm{c}\mathrm{o}\mathrm{s}\delta }& \frac{{R}_{ny}^{z}}{\mathrm{c}\mathrm{o}\mathrm{s}\beta }+{R}_{nx}^{z}\mathrm{t}\mathrm{a}\mathrm{n}\beta \mathrm{t}\mathrm{a}\mathrm{n}\delta & 0& 1\end{array}\right]}^{T},\\ {\mathbf{R}}_{d}^{j}\left(\theta \right)={\left[\begin{array}{cccc}\frac{{R}_{nx}^{j}}{\mathrm{c}\mathrm{o}\mathrm{s}\delta }& \frac{{R}_{ny}^{j}}{\mathrm{c}\mathrm{o}\mathrm{s}\beta }+{R}_{nx}^{j}\mathrm{t}\mathrm{a}\mathrm{n}\beta \mathrm{t}\mathrm{a}\mathrm{n}\delta & 0& 1\end{array}\right]}^{T}.\end{array}\right\$

### 2.2. Coordinate transformation of rack cutter

Fig. 3 displays the coordinate systems of rack cutter. Herein, the coordinate system ${S}_{d}$ moving a given translational distance $u$ along negative direction of ${Z}_{d}$ axis and rotating it around the ${X}_{d}$ axis with an angle $\beta$, the coordinate system ${S}_{c}\left({X}_{c},{Y}_{c},{Z}_{c}\right)$ will be got by rotating ${S}_{p}$ around the ${Y}_{p}$ axis with $\delta$.

Fig. 3Coordinate systems of rack cutter

In the coordinate system ${S}_{c}$, the tooth surface equation of rack cutter can be expressed as:

5
${\mathbf{R}}_{c}={\mathbf{M}}_{cd}{\mathbf{R}}_{d}={\mathbf{M}}_{cp}{\mathbf{M}}_{pd}{\mathbf{R}}_{d},$

where ${\mathbf{R}}_{c}$ and ${\mathbf{R}}_{d}$ are position vector of rack cutter in coordinate system ${S}_{c}$ and ${S}_{d}$, respectively; ${\mathbf{M}}_{pd}$ and ${\mathbf{M}}_{cp}$ are coordinate transformation matrix.

Combining Eqs. (4) and (5), we can obtain the position vector of $AB$ and $AB$ in coordinate system ${S}_{c}$ as follows:

6
${\mathbf{R}}_{c}^{z}=\left[\begin{array}{c}{R}_{dx}^{z}+ucos\beta sin\delta \\ {R}_{dy}^{z}+usin\beta \\ ucos\beta cos\delta \\ 1\end{array}\right],{\mathbf{R}}_{c}^{j}=\left[\begin{array}{c}{R}_{dx}^{j}+ucos\beta sin\delta \\ {R}_{dy}^{j}+usin\beta \\ ucos\beta cos\delta \\ 1\end{array}\right].$

### 2.3. Meshing equation

In coordinate system ${S}_{c}$, equations of arbitrary point on $AB$ and $AB$ are the function which depends on the variables $l/u$ and $\theta /u$, respectively. The unit normal vectors of the rack cutter can be obtained as:

7
$\frac{\frac{\partial {\mathbf{R}}_{c}^{z}}{\partial l}×\frac{\partial {\mathbf{R}}_{c}^{z}}{\partial u}}{\left|\frac{\partial {\mathbf{R}}_{c}^{z}}{\partial l}×\frac{\partial {\mathbf{R}}_{c}^{z}}{\partial u}\right|}=\left[\begin{array}{c}{n}_{cx}^{z}\\ {n}_{cy}^{z}\\ {n}_{cz}^{z}\end{array}\right],\frac{\frac{\partial {\mathbf{R}}_{c}^{j}}{\partial \theta }×\frac{\partial {\mathbf{R}}_{c}^{j}}{\partial u}}{\left|\frac{\partial {\mathbf{R}}_{c}^{j}}{\partial \theta }×\frac{\partial {\mathbf{R}}_{c}^{j}}{\partial u}\right|}=\left[\begin{array}{c}{n}_{cx}^{j}\\ {n}_{cy}^{j}\\ {n}_{cz}^{j}\end{array}\right].$

Coordinate relationship between the rack cutter and involute beveloid gear is illustrated in Fig. 4. ${S}_{b}\left({X}_{b},{Y}_{b},{Z}_{b}\right)$ is a spatial fixed coordinate system. ${S}_{j}\left({X}_{j},{Y}_{j},{Z}_{j}\right)$ is a moving coordinate system attached to gear blank (involute beveloid gear), which forms an angle ${\phi }_{1}$ with respect to the coordinate system ${S}_{b}\left({X}_{b},{Y}_{b},{Z}_{b}\right)$.

Fig. 4Coordinate relationship between rack cutter and involute beveloid gear

The angular velocity vector revolving around the ${Z}_{j}$ axis of involute beveloid gear can be represented as:

8
${\mathbf{\omega }}_{1}=-{\omega }_{1}{\mathbf{k}}_{j}.$

Assuming that the position vector of one point $K$ on rack cutter is described as:

9
$\stackrel{\to }{{O}_{c}K}={\mathbf{r}}_{c}={x}_{c}{\mathbf{i}}_{c}+{y}_{\text{c}}{\mathbf{j}}_{c}+{z}_{\text{c}}{\mathbf{k}}_{c}.$

Then we can obtain that:

10
$\stackrel{\to }{PK}=\stackrel{\to }{P{O}_{c}}+\stackrel{\to }{{O}_{c}K}={x}_{c}{\mathbf{i}}_{c}+\left({y}_{c}-{r}^{\text{'}}{\phi }_{1}\right){\mathbf{j}}_{c}+{z}_{c}{\mathbf{k}}_{c}.$

Combining Eqs. (8) and (10), the relative velocity vector of one point $K$ on rack cutter is:

11
${\mathbf{v}}_{12}^{K}={\mathbf{\omega }}_{\text{1}}×\stackrel{\to }{PK}={\omega }_{\text{1}}\left[\left({y}_{c}-{r}^{\text{'}}{\phi }_{1}\right){\mathbf{i}}_{c}+{x}_{c}{\mathbf{j}}_{c}\right].$

According to differential geometry and the kinematics of gear geometry [12], the continuous tangency is detected by the equation of meshing which is formulated as:

12
${\mathbf{n}}^{K}\cdot {\mathbf{v}}_{12}^{K}=0,$

where the unit normal vectors of rack cutter ${\mathbf{n}}^{K}$ is ${\mathbf{n}}^{K}={\left(\begin{array}{ccc}{n}_{x}^{K}& {n}_{y}^{K}& {n}_{z}^{K}\end{array}\right)}^{T}$.

Then the motion parameter (rotation angle) for involute beveloid gear generation is:

13
${\phi }_{1}=\frac{{y}_{c}}{{r}^{\text{'}}}+\frac{{n}_{y}^{K}{x}_{c}}{{n}_{x}^{K}{r}^{\text{'}}}.$

Substituting Eqs. (6) and (7) into Eq. (13) enables us to solve the rotation angle of straight edge $BC$ and tool fillet curve $AB$ for rack cutter:

14
${\phi }_{1}^{z}=\frac{{R}_{cy}^{z}}{{r}^{\text{'}}}+{n}_{cy}^{z}\cdot \frac{{R}_{cx}^{z}}{{n}_{cx}^{z}\cdot {r}^{\text{'}}},\mathrm{}{\phi }_{1}^{j}=\frac{{R}_{cy}^{j}}{{r}^{\text{'}}}+{n}_{cy}^{j}\cdot \frac{{R}_{cx}^{j}}{{n}_{cx}^{j}\cdot {r}^{\text{'}}}.$

### 2.4. Tooth profile equation of involute beveloid gear

In coordinate system ${S}_{j}$, the tooth surface equation of involute beveloid gear is expressed as:

15
$\left[\begin{array}{c}{R}_{jx}^{z}\\ {R}_{jy}^{z}\\ {R}_{jz}^{z}\\ 1\end{array}\right]={\mathbf{M}}_{jc}\left[\begin{array}{c}{R}_{cx}^{z}\\ {R}_{cy}^{z}\\ {R}_{cz}^{z}\\ 1\end{array}\right]=\left[\begin{array}{c}\left({R}_{cx}^{z}+{r}^{\text{'}}\right)cos{\phi }_{1}^{z}+\left({r}^{\text{'}}{\phi }_{1}^{z}-{R}_{cy}^{z}\right)sin{\phi }_{1}^{z}\\ \left({R}_{cx}^{z}+{r}^{\text{'}}\right)sin{\phi }_{1}^{z}-\left({r}^{\text{'}}{\phi }_{1}^{z}-{R}_{cy}^{z}\right)cos{\phi }_{1}^{z}\\ {R}_{cz}^{z}\\ 1\end{array}\right],$

where ${\mathbf{M}}_{jc}$ is coordinate transformation matrix.

Similarly, we can also obtain the fillet equation of involute beveloid gear:

16
$\left[\begin{array}{c}{R}_{jx}^{j}\\ {R}_{jy}^{j}\\ {R}_{jz}^{j}\\ 1\end{array}\right]=\left[\begin{array}{c}\left({R}_{cx}^{j}+{r}^{\text{'}}\right)cos{\phi }_{1}^{j}+\left({r}^{\text{'}}{\phi }_{1}^{j}-{R}_{cy}^{j}\right)sin{\phi }_{1}^{j}\\ \left({R}_{cx}^{j}+{r}^{\text{'}}\right)sin{\phi }_{1}^{j}-\left({r}^{\text{'}}{\phi }_{1}^{j}-{R}_{cy}^{j}\right)cos{\phi }_{1}^{j}\\ {R}_{cz}^{j}\\ 1\end{array}\right].$

Hence, the mathematical models of the involute beveloid gear have been derived.

## 3. Mathematical model of noninvolute beveloid gear pair

### 3.1. Coordinate transformation

According to the differential geometry and the kinematics of gear geometry, the tooth profile equation and meshing equation of noninvolute beveloid gear can be mathematically generated from a mutually conjugate involute beveloid gear. Fig. 5 describes the relationship between involute beveloid gear cutter and gear blank (noninvolute beveloid gear), including the involute beveloid gear Gear_1, the noninvolute beveloid gear Gear_2 and the intersection angle $\mathrm{\Sigma }$.

The coordinate system ${S}_{3}\left({X}_{3},{Y}_{3},{Z}_{3}\right)$ and ${S}_{4}\left({X}_{4},{Y}_{4},{Z}_{4}\right)$ are fixed coordinate systems linked with Gear_1 and Gear_2, respectively, and the coordinate system ${S}_{1}\left({X}_{1},{Y}_{1},{Z}_{1}\right)$ and ${S}_{2}\left({X}_{2},{Y}_{2},{Z}_{2}\right)$ are motion coordinate systems joined to Gear_1 and Gear_2, respectively as illustrated in Fig. 6.

Fig. 5Position relationship between the involute beveloid gear cutter and gear blank

Under the initial position, the coordinate system ${S}_{1}$ overlaps with ${S}_{3}$ and the coordinate system ${S}_{2}$ coincides with ${S}_{4}$. The Gear_2 rotates clockwise around the ${Z}_{4}$ axis with an angle ${\phi }_{2}$, while the Gear_1 is revolving clockwise around the ${Z}_{3}$ axis with an angle ${\phi }_{1}$.

Fig. 6Coordinate systems of Gear_1 and Gear_2

Consequently, coordinate transformation matrix ${M}_{21}\text{,}$ which is utilized to perform the kinematical relationship between the involute beveloid gear cutter and the generated gear (noninvolute beveloid gear), can be expressed as ${M}_{21}={M}_{24}{M}_{43}{M}_{31}$, where:

17
${M}_{24}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{2}& -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{2}& 0& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{2}& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{2}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],{M}_{43}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }& 0& -\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }& 0\\ 0& 1& 0& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }& 0& \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }& 0\\ 0& 0& 0& 1\end{array}\right],$
${M}_{31}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}& -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}& 0& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right].$

### 3.2. Relative velocity and unit normal vectors of involute beveloid gear cutter

Assuming that angular velocity vector clockwise revolving around the axis of Gear_1 is ${\omega }_{1}$ and the angular velocity vector anticlockwise revolving around the axis of Gear_2 is ${\omega }_{2}$. Consequently, the angular velocity vector can be represented based on three components of unit vector in coordinate system ${S}_{3}\left({X}_{3},{Y}_{3},{Z}_{3}\right)$:

18
${\mathbf{\omega }}_{1}={\omega }_{1}{\mathbf{k}}_{3},\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathbf{\omega }}_{2}=-{\omega }_{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }{\mathbf{k}}_{3}-{\omega }_{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }{\mathbf{i}}_{3}.$

The position vector of one contact point $M$ is described as:

19
${\mathbf{r}}_{3}=\left({X}_{3},{Y}_{3},{Z}_{3}\right).$

Then the relative velocity at the point $M$ can be expressed as:

20
${\mathbf{v}}_{12}=\left({\mathbf{\omega }}_{1}-{\mathbf{\omega }}_{2}\right)×{\mathbf{r}}_{3}.$

Substituting Eqs. (18) and (19) into Eq. (20) enables us to obtain the relative velocity:

21
${\mathbf{v}}_{12}^{3}=\left(\begin{array}{ccc}{\mathbf{i}}_{3}& {\mathbf{j}}_{3}& {\mathbf{k}}_{3}\\ {\omega }_{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }& 0& {\omega }_{1}+{\omega }_{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\\ {X}_{3}& {Y}_{3}& {Z}_{3}\end{array}\right)=\left(\begin{array}{c}\left[-{Y}_{3}\left({\omega }_{1}+{\omega }_{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)\right]{\mathbf{i}}_{3}\\ \left[-{Z}_{3}{\omega }_{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }+{X}_{3}\left({\omega }_{1}+{\omega }_{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)\right]{\mathbf{j}}_{3}\\ \left[{Y}_{3}{\omega }_{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\right]{\mathbf{k}}_{3}\end{array}\right).$

Substituting ${i}_{21}={\omega }_{2}/{\omega }_{1}$ into Eq. (21), ${\mathbf{v}}_{12}^{3}$ can be further written as:

22
${\mathbf{v}}_{12}^{3}={\omega }_{1}\left(\begin{array}{c}\left[-{Y}_{3}\left(1+{i}_{21}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)\right]{\mathbf{i}}_{3}\\ \left[-{Z}_{3}{i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }+{X}_{3}\left(1+{i}_{21}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)\right]{\mathbf{j}}_{3}\\ \left[{Y}_{3}{i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\right]{\mathbf{k}}_{3}\end{array}\right).$

As a consequence, the relative velocity in coordinate system ${S}_{1}$ can be obtained as follows:

23
${\mathbf{v}}_{12}^{1}={\mathbf{M}}_{13}{\mathbf{v}}_{12}^{3}={\omega }_{1}\left(\begin{array}{c}\left[-{Y}_{1}\left(1+{i}_{21}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)-{Z}_{1}{i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\right]{\mathbf{i}}_{1}\\ \left[{X}_{1}\left(1+{i}_{21}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)-{Z}_{1}{i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}\right]{\mathbf{j}}_{1}\\ \left[\left({Y}_{1}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}+{X}_{1}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\right){i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\right]{\mathbf{k}}_{1}\end{array}\right),$

where matrix ${\mathbf{M}}_{13}$ describes coordinate transformation from ${S}_{3}$ to ${S}_{1}\text{;}$$\left({X}_{1},{Y}_{1},{Z}_{1}\right)$ means coordinate value of contact point $M$ in coordinate system ${S}_{1}$, and $\left({\mathbf{i}}_{1},{\mathbf{j}}_{1},{\mathbf{k}}_{1}\right)$ represents three components of unit vector in coordinate system ${S}_{1}$.

According to the position relationship between the involute beveloid gear cutter and noninvolute beveloid gear as shown in Fig. 5, the tooth surface of involute beveloid gear cutter in coordinate system ${S}_{1}$ can be derived as follows:

24
$\left[\begin{array}{c}{R}_{1x}\\ {R}_{1y}\\ {R}_{1z}\\ 1\end{array}\right]=\left[\begin{array}{c}{R}_{jx}^{z}\\ {R}_{jy}^{z}\\ {R}_{jz}^{z}+L\\ 1\end{array}\right],$

where $L$ is coordinate value on the ${Z}_{1}$ axis.

Hence, the unit normal vector of Gear_1 (involute beveloid gear) is represented as:

25
$\frac{\frac{\partial {R}_{1}}{\partial l}×\frac{\partial {R}_{1}}{\partial u}}{\left|\frac{\partial {R}_{1}}{\partial l}×\frac{\partial {R}_{1}}{\partial u}\right|}=\left[\begin{array}{c}\frac{{n}_{1x}}{\sqrt{{n}_{1x}^{2}+{n}_{1y}^{2}+{n}_{1z}^{2}}}\\ \frac{{n}_{1y}}{\sqrt{{n}_{1x}^{2}+{n}_{1y}^{2}+{n}_{1z}^{2}}}\\ \frac{{n}_{1z}}{\sqrt{{n}_{1x}^{2}+{n}_{1y}^{2}+{n}_{1z}^{2}}}\end{array}\right]=\left[\begin{array}{c}{n}_{1x}^{\text{'}}\\ {n}_{1y}^{\text{'}}\\ {n}_{1z}^{\text{'}}\end{array}\right].$

### 3.3. Meshing equation and tooth surface equation of noninvolute beveloid gear

According to the kinematics of gear geometry, the equation of meshing is defined as:

26
$\mathbf{n}\cdot {\mathbf{v}}_{12}=\text{0.}$

Combining Eqs. (23), (25) and (26), we can obtain the meshing equation:

27
$\left[\begin{array}{c}{n}_{1x}^{\text{'}}\\ {n}_{1y}^{\text{'}}\\ {n}_{1z}^{\text{'}}\end{array}\right]\cdot {\omega }_{1}\left(\begin{array}{c}\left[-{Y}_{1}\left(1+{i}_{21}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)-{Z}_{1}{i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\right]{\mathbf{i}}_{1}\\ \left[{X}_{1}\left(1+{i}_{21}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)-{Z}_{1}{i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}\right]{\mathbf{j}}_{1}\\ \left[\left({Y}_{1}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}+{X}_{1}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\right){i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\right]{\mathbf{k}}_{1}\end{array}\right)=0,$

where ${\phi }_{1}$ is angle parameter, which can be derived as:

28
${\phi }_{1}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\frac{\left(ca-\sqrt{{a}^{2}{b}^{2}+{b}^{4}-{b}^{2}{c}^{2}}\right)a}{\left({a}^{2}+{b}^{2}\right)b}-\frac{c}{b}}{-\frac{ca-\sqrt{{a}^{2}{b}^{2}+{b}^{4}-{b}^{2}{c}^{2}}}{{a}^{2}+{b}^{2}}}\right),$

herein:

$a={i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\left({Y}_{1}{n}_{1z}-{Z}_{1}{n}_{1y}\right),b={i}_{21}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\left({X}_{1}{n}_{1z}-{Z}_{1}{n}_{1x}\right),$
$c=\left(1+{i}_{21}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }\right)\left({X}_{1}{n}_{1y}-{Y}_{1}{n}_{1x}\right).$

Then the tooth profile equation of noninvolute beveloid gear can be obtained by coordinate transformation from ${S}_{1}$ to ${S}_{2}$:

29
$\left[\begin{array}{c}{R}_{2x}\\ {R}_{2y}\\ {R}_{2z}\\ 1\end{array}\right]={M}_{21}\left[\begin{array}{c}{R}_{1x}\\ {R}_{1y}\\ {R}_{1z}\\ 1\end{array}\right].$

As a consequence:

30
$\left\{\begin{array}{l}{R}_{2x}=\left(\begin{array}{c}\left(\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }-\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{2}\right){R}_{1x}\\ +\left(-\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }-\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{2}\right){R}_{1y}+\left(-\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\right){R}_{1z}\end{array}\right),\\ {R}_{2y}=\left(\begin{array}{c}\left(\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }+\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{2}\right){R}_{1x}\\ +\left(-\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }+\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{2}\right){R}_{1y}+\left(-\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }\right){R}_{1z}\end{array}\right),\\ {R}_{2z}=\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{1}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }{R}_{1x}+-\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{1}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Sigma }{R}_{1y}+\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Sigma }{R}_{1z}.\end{array}\right\$

Hence, the mathematical model of involute beveloid gear have been derived.

## 4. Modeling of involute and noninvolute beveloid gear pair

The main design parameters of involute beveloid gear pairs (pinion and wheel) are: number of teeth ${Z}_{1}/{Z}_{2}=$26/51, normal pressure angle ${\alpha }_{n1}/{\alpha }_{n2}=$ 17.5°, normal module ${m}_{n1}/{m}_{n2}=$ 5 mm, tooth width ${b}_{1}/{b}_{2}=$ 94 mm/89 mm, helical angle $\beta =$8°, shaft intersection angle $\delta =$10°; the main design parameters of noninvolute beveloid gear pairs (pinion and wheel) are: number of teeth ${Z}_{1}^{\text{'}}/{Z}_{2}=$ 26/51, the normal module ${m}_{n1}/{m}_{n2}=$ 5 mm, the tooth width ${b}_{1}/{b}_{2}=$ 94 mm/89 mm and shaft intersection angle $\delta =$10°.

The parametric modeling program is developed via a Matlab code to generate exact three-dimensional (3-D) discrete points on the tooth profiles. The Fig. 7(a) shows the discrete points on tooth surfaces of involute beveloid gear (wheel). Then the tooth surfaces of involute beveloid gear pair are established after we imported point set into Imageware software. Subsequently, the solid model will be created by importing tooth surfaces into UG to sew up. In this way, we can obtain the accurate three-dimensional solid model as shown in Fig. 7(b).

Applying the assembly relationship together with the mathematical models of beveloid gear pairs, the assembly models of involute and noninvolute beveloid gear pairs are generated as illustrated in Fig. 8.

Fig. 7Models of involute beveloid gears

a) Discrete points on tooth surface

b) Solid models

Fig. 8Assembly model of involute and noninvolute beveloid gears

## 5. Dynamic contact analysis of involute and noninvolute beveloid gear pair

### 5.1. Finite element model

Dynamic contact characteristics of the involute and noninvolute beveloid gear pair will be analyzed by LS-DYNA software. The pinion is rotating at 1600 rpm with a driving torque of 1194 N·m in numerical simulation analysis. The finite element meshes of gear pair are generated using solid element Solid164. In order to conveniently load the speed and torque [13], internal surfaces of gear pair are considered as rigid regions utilizing the Shell163 element. Afterwards, degrees of freedom of all nodes in rigid region are coupled to the centroid of rigid body, and apart from axial rotation, all translational and rotational degrees of freedom of rigid shell elements in pinion and wheel are constrained. The total number of elements is 82160 with 97608 nodes for the finite element model of involute beveloid helical gear pair, shown as Fig. 9.

Fig. 9Finite element model of involute beveloid helical gear pairs

### 5.2. Simulation results and discussion

Fig. 10 and Fig. 11 describe the dynamic contact stress contour and single tooth contact force of involute and noninvolute beveloid gear during the mesh cycle.

The results show that the contact area on flank of noninvolute beveloid gear is evidently larger than that of involute beveloid gear; the dynamic contact stress distribution of noninvolute beveloid gear is more homogeneous and the contact stress reduces significantly; the contact ratio of noninvolute beveloid gear is 3.74, greater than that of involute beveloid gear, which is 3.30; compared with the single tooth contact force of involute beveloid gear, the noninvolute beveloid gear has smaller single tooth contact force and smoother meshing impact.

Fig. 10Dynamic contact stress contour (Unit: Pa)

a) Involute beveloid gears

b) Noninvolute beveloid gears

Fig. 11Contact force of each mating teeth pair

a) Involute beveloid gears

b) Noninvolute beveloid gears

c) Comparison of single tooth contact force

Fig. 12 and Fig. 13 illustrate the dynamic fillet stress of involute and noninvolute beveloid gear, respectively. The results indicate that the dynamic fillet stress of involute beveloid gear mainly concentrate on the mid tooth and toe, by contrast, the heel almost do not bear fillet stress; the dynamic fillet stress distribution of noninvolute beveloid gear is more homogeneous on mid tooth, toe and heel.

Forecasting the dynamic transmission error (DTE) which is an important factor to cause the vibration, is another objective of dynamic contact analysis. The dynamic transmission error is usually expressed as:

31
$\text{DTE}\left(t\right)={n}_{\text{2}}-{n}_{\text{1}}\text{/}i,$

where ${n}_{1}$ and ${n}_{2}$ is the driving and driven shaft’s rotational speeds with $i$ being the transmission ratio of the gear set, respectively.

Fig. 12Dynamic fillet stress curves of involute beveloid gear pairs

a) Pinion

b) Wheel

Fig. 13Dynamic fillet stress curves of noninvolute beveloid gear pairs

a) Pinion

b) Wheel

Fig. 14 illustrates the DTE of involute and noninvolute beveloid gear. As observed in Fig. 14, the predicted DTE is approximate sine shape and the trend of peak-peak value of DTE for involute beveloid gear becomes higher than that of noninvolute beveloid gear due to smaller contact area for involute beveloid gear.

Fig. 14DTE of involute and noninvolute beveloid gear

## 6. Conclusions

1) The mathematical models of involute and noninvolute beveloid gears are derived and the parametric modeling programs of involute and noninvolute beveloid gears are developed to automatically generate exact model via Matlab code; then a numerical example of intersecting axes beveloid gear pair is presented to analyze the dynamic contact characteristics.

2) The contact area on flank of noninvolute beveloid gear is evidently larger than that of involute beveloid gear, and the dynamic contact stress distribution of noninvolute beveloid gear is more homogeneous and the contact stress reduces significantly.

3) The dynamic fillet stress of involute beveloid gear mainly concentrate on the mid tooth and toe, by contrast, the heel almost do not bear fillet stress, the dynamic fillet stress distribution of noninvolute beveloid gear is more homogeneous on mid tooth, toe and heel.

4) The DTE is approximate sine shape and the trend of peak-peak value for involute beveloid gear becomes higher than that of noninvolute beveloid gear due to smaller contact area for involute beveloid gear and the noninvolute beveloid gear has gentler meshing impact.

5) Hopefully, noninvolute beveloid gear pairs would relieve the high dynamic contact shock problem of intersecting axes beveloid gear pairs and improve range of application.

#### References

• Zhang Y.,Fang Z. Analysis of tooth contact and load distribution of helical gears with crossed axes. Mechanism and Machine Theory, Vol. 34, Issue 1, 1999, p. 41-57.
• Komatsubara H., Mitome K., Ohmachi T. Development of concave conical gear used for marine transmissions. JSME International Journal, Series C: Mechanical Systems, Machine Elements and Manufacturing, Vol. 45, Issue 1, 2002, p. 371-377.
• Brauer J. Analytical geometry of straight conical involute gears. Mechanism and Machine Theory, Vol. 37, Issue 1, 2002, p. 127-141.
• Brauer J. A general finite element model of involute gears. Finite Elements in Analysis and Design, Vol. 40, Issue 13-14, 2004, p. 1857-1872.
• Brauer J. Transmission error in anti-backlash conical involute gear transmissions: a global – local FE approach. Finite Elements in Analysis and Design, Vol. 41, Issue 5, 2005, p. 431-457.
• Liu C. C., Tsay C. B. Contact characteristics of beveloid gears. Mechanism and Machine Theory, Vol. 37, Issue 4, 2002, p. 333-350.
• Liu C. C., Chung B. Mathematical models and contact simulations of concave beveloid gears. Journal of Mechanical Design, Transactions of the ASME, Vol. 124, Issue 4, 2002, p. 753-760.
• Chen Y. C.,Liu C. C. Contact stress analysis of concave conical involute gear pair with non-parallel axes. Finite Elements in Analysis and Design, Vol. 47, Issue 4, 2011, p. 443-452.
• Chen Chen-Jia, Liu Wen, Lin Teng-Jiao, et al. Modeling of involute beveloid helical gears with intersection axes based on transverse profile of counterpart rack. Chinese Journal of Mechanical Research and Application, Vol. 26, Issue 125, 2012, p. 1-7.
• Wu S. H., Tsai S. J. Contact stress analysis of skew conical involute gear drives in approximate line contact. Mechanism and Machine Theory, Vol. 44, Issue 9, 2009, p. 1658-1676.
• Li Gui-Xian, Wen Jian-Min, Li Xiao, et al. Axial errors and profile errors of noninvolute bevoloid gears. Chinese Journal of Harbin Engineering University, Vol. 24, Issue 3, 2003, p. 302-304, 316.
• Litvin F. L.,Fuentes A. Gear Geometry and Applied Theory. Cambridge University Press, New York, 2010.
• Lin Teng-Jiao, Ou Heng-An, Li Run-Fang A finite element method for 3D static and dynamic contact/impact analysis of gear drives. Computer Methods in Applied Mechanics and Engineering, Vol. 196, Issue 9-12, 2007, p. 1716-1728.