Optimization of bias modification and dynamic behavior analysis of helical gear system

Abstract

A generalized finite element dynamic model of helical gear system considering time-varying mesh stiffness and loaded composite mesh error is developed. Bias modification is introduced into this model, and the corresponding time-varying mesh stiffness and loaded composite mesh error are determined based on improved loaded tooth contact analysis model. By introducing vibration exciting force, the time-variant differential equations of the system are transformed into time-invariant differential equations which can be solved faster by Fourier series method. Meanwhile, the internal excitation of the system can be captured clearly. Aiming to reduce vibration of the system, optimization method of bias modification is proposed with the objective of minimizing fluctuation of the vibration exciting force. The global optimal modification parameters are obtained using brute force optimization method. The sensitivity of dynamic response of the modified helical gear system to gear misalignment is then investigated. Simulation results suggest that bias modification has little effect on time-varying mesh stiffness but reduces system vibration significantly. Increasing gear misalignment will decrease time-varying mesh stiffness and increase the system vibration. Bias modification shows a well robustness in reducing vibration of helical gear system when gear misalignment is not too large.

Introduction

Due to higher contact ratio, larger load carrying capacity, lower vibration, and noise, helical gear transmission systems are widely used in many industry applications, such as automotive, wind turbine, mining, marine, and industrial power transmissions. With the increasingly tight restrictions on the vibration and noise of gear systems, the prediction and control of gear vibration and noise are a crucial issue in the design of the gear systems. It is well known that gear modification is an advanced technology in reducing vibration and noise of gear transmission. To obtain reasonable modification parameters for the gear design and processing is of great importance.

Since the profile modification method was proposed by H Walker ¹ for involute gear pair at first in 1938, numerous researches focused on the profile modification for spur gear system. Some of research works were related to the quasi-static contact analysis of spur gear system with profile modification,^2,3 aiming at minimizing the amplitude of static transmission error,^4,5 while others were devoted to investigate the effect of profile modification on dynamic behaviors of spur gear system,^6–11 aiming at minimizing the dynamic response of the system. ¹² Meanwhile, some of published literature focused on tooth surface modification for reduction of vibration and noise of helical gear system. Unlike spur gears, the calculation of mesh stiffness and load distribution are a complicated three-dimensional contact problem for a helical gear pair under load because the applied loads are not evenly distributed along the line of contact. Conry and Seireg^13,14 proposed a universal contact model of two elastic bodies based on compliance analysis method to the calculation of load distribution and elastic deflection of gear tooth and developed gear modification method for spur and helical gear. Simon ¹⁵ presented a method for the calculation of optimal tooth modifications for spur and helical gears aimed at minimizing load distribution factor. Umeyama et al. ¹⁶ developed a transmission error calculation model for cylindrical gear based on the finite element method and Hertzian contact analytical formula and studied the influences of gear-dimensional parameters and tooth surface modification on the transmission error. Rao and Yoon ¹⁷ presented a method of computing the three-dimensional variation in transmission error of helical gears on the basis of the analysis of the plane of action generated by a pair of mating teeth and used cubic splines to determine the optimal helical tooth profile. Wang and Zhang ¹⁸ proposed an analytical model for calculating time-varying mesh stiffness of helical gears with tooth profile error based on potential energy method and thin slice theory and analyzed the effects of tooth errors and gear modification on the loaded tooth contact characteristics of helical gears. Bruyère and colleagues^19,20 proposed a perturbation method to obtain approximate closed-form expressions for profile relief that minimize the amplitude of quasi-static transmission error under load and presented two original analytical formulations which give the amount and length of profile modifications minimizing the amplitude of quasi-static transmission error in narrow-faced spur and helical gears. Maatar and Velex ²¹ investigated the effects of profile modification and lead modification on static transmission error, load distribution, and dynamic transmission error of narrow-faced helical gears. Wu et al. ²² established a static contact finite element model of helical gear pair in ANSYS software to provide a profile modification scheme based on elastic deformation of gear teeth along the line of action and used dynamic contact finite element model to investigate the effect of profile modification on dynamic responses of the system.

In addition, some current researches are devoted to study robust optimization method for gear modification design of cylindrical gears. Park ²³ proposed a method to determine optimal lead curve and robust tooth surface design, using the response surface method and multi-objective optimization. Ghribi et al. ²⁴ combined analytical results and numerical simulations to design a robust profile modification to spur and helical gears. Artoni et al. ²⁵ proposed optimization method of tooth surface modification for cylindrical gear to achieve low sensitive of loaded static transmission error to applied torque and misalignment.

With published literature mentioned above, numerous researches have been carried out to investigate the effect of profile modification and lead modification on the quasi-static and dynamic characteristics for spur and helical gear system. However, as an effective method to reduce vibration and noise for helical gear system, bias modification has not been studied a lot at present, especially the loaded tooth contact characteristics of helical gear pair with bias modification and the effect of bias modification on reduction of vibration and noise of helical gear system have not been analyzed in detail. Meanwhile, the current researches lack the comprehensive and thorough investigation of the effects of gear misalignment on loaded tooth contact characteristics and dynamic behaviors of helical gear system with bias modification.

The main motivation of this study is to investigate the effect of bias modification on loaded tooth contact characteristics of helical gears, obtain the optimal bias modification parameters quickly, and analyze the sensitivity of dynamic responses of helical gear system with bias modification to gear misalignment based on a more precise generalized finite element dynamic model of the system. After this introduction, the rest of this article is organized as follows: section “Bias modification and dynamic model” gives the design of bias modification for helical gears and the calculation method of corresponding time-varying mesh stiffness and loaded composite mesh error based on the improved loaded tooth contact analysis (LTCA) model. Section “Vibration exciting force” develops a more precise generalized finite element dynamic model of helical gear system based on the finite element method and Timoshenko beam theory in order to consider shaft deflection, bearing location, and the direction of power flow. By introducing vibration exciting force, the time-variant motion differential equations of system are transformed into time-invariant equations and the dynamic responses can be obtained faster by Fourier series method. Section “Numerical results and discussion” gives a comparison of solution method in order to validate the proposed approximate transformation in section “Vibration exciting force.” Aiming to minimize the fluctuation of vibration exciting force, the optimal bias modification parameters are obtained by global brute force optimization method in potential ranges, and the effects of gear misalignment on dynamic responses of helical gear system with bias modification are investigated. Section “Conclusion” draws the main conclusions of the numerical simulation results of helical gear system with bias modification. The detailed simulation flowchart is shown in Figure 1.

Figure 1.

Flowchart of the simulation.

Bias modification and dynamic model

Description of bias modification for helical gear

As removing material from the bias of gear tooth along the contact line, bias modification consists of three basic parameters: the modification curve, the amount of modification, and the length of modification. Schematic diagram of bias modification of helical gears on tooth surface and plane of action are displayed in Figure 2, where the shaded region is the bias modification region, and B₁B₂B₃B₄ is the plane of action. The X₁OY₁ is the Cartesian coordinate system to define the plane of action, and O-X₂Y₂Z₂ is the Cartesian coordinate system for calculation. L is the normal length of bias modification, and La is the maximum potential length of bias modification. In this section, the symmetric bias modification on plane of action is adopted because tip modification of the pinion is equivalent to root modification of the wheel.

Figure 2.

Schematic diagram of bias modification for helical gears.

The equation of bias modification curve in the X₂OZ₂ coordinate system can be written as

$$e_i=\{\begin{array}{cc}{e}_{max}{\left(\frac{x_i+\overline{\mathrm{OG}}}{\overline{\mathrm{EG}}}\right)}^n,\hfill & x_i<-\overline{\mathrm{OG}}\hfill \\ \multicolumn{1}{c}{0,}& -\overline{\mathrm{OG}}\le x_P\le \overline{\mathrm{OH}}\hfill \\ \multicolumn{1}{c}{e_{max}{\left(\frac{x_i-\overline{\mathrm{OH}}}{\overline{\mathrm{FH}}}\right)}^n,}& x_i>\overline{\mathrm{OH}}\hfill \end{array}$$

(1)

where e_i is the modification amount of potential contact point i; x_i is the coordinate of contact point i along X₂-axis in the O-X₂Y₂Z₂ coordinate system; e_max is the maximum modification amount; n = 1 denotes that the modification curve is a straight line; and n = 2 denotes that the modification curve is a parabola curve, respectively. Parabola curve is adopted for bias modification in this article.

Mesh stiffness and loaded composite mesh error

The dynamic engagement process of a gear pair can be regarded as the quasi-static loaded contact process of two elastic bodies, as shown in Figure 3. When the external force P is applied, the two elastic bodies will be close to each other and contact. After the continuous contact lines are discretized into a series of contact points, the contact condition of possible contact point i for each engagement position will be described as

$${\delta }_i^{(1)}+{\delta }_i^{(2)}+{\epsilon }_i-\mathrm{LSTE}-d_i=0$$

(2)

where ${\delta }_i^{(1)}$ and ${\delta }_i^{(2)}$ denote deformation of pinion and wheel, respectively; ε_i is the initial clearance of contact point; and LSTE represents the rigid body approach of the two bodies. For a gear pair under load, LSTE is the loaded static transmission error, which consists of elastic deformation and gear errors. d_i refers to residual clearance after loading.

Figure 3.

Schematic diagram of loaded tooth contact of gear pair with errors.

The deformation of gears can be divided into two categories: global deformation which changes linearly with the applied force and local deformation which changes nonlinearly with the applied force. Thus, the equations can be rewritten as

$${\delta }_{\mathrm{Global}\_i}^{(1)}+{\delta }_{\mathrm{Global}\_i}^{(2)}+{\delta }_{\mathrm{Contact}\_i}+{\epsilon }_i-\mathrm{LSTE}-d_i=0$$

(3)

where ${\delta }_{\mathrm{Global}\_i}^{(1)}$ and ${\delta }_{\mathrm{Global}\_i}^{(2)}$ refer to the global deformation of pinion and wheel, respectively, and δ_{Contact_i} denotes the local contact deformation of mating gear teeth.

As global deformation of pinion and wheel is related to applied force P linearly, they can be expressed as

$$\{\begin{array}{c}{\delta }_{\mathrm{Global}\_i}^{(1)}=\sum _{j=1}^n{\lambda }_{\mathrm{Global}\_\mathrm{ij}}^{(1)}{F}_j\hfill \\ \multicolumn{1}{c}{{\delta }_{\mathrm{Global}\_i}^{(2)}=\sum _{j=1}^n{\lambda }_{\mathrm{Global}\_\mathrm{ij}}^{(2)}{F}_j}\end{array}$$

(4)

where ${\lambda }_{\mathrm{Global}\_\mathrm{ij}}^{(1)}$ and ${\lambda }_{\mathrm{Global}\_\mathrm{ij}}^{(2)}$ denote the global deformation of contact point i when a unitary force is applied on contact point j along the normal direction of tooth surface for pinion and wheel, respectively.

The global deformation compliance of gears can be obtained as

$${\lambda }_{\mathrm{Global}\_\mathrm{ij}}={\lambda }_{\mathrm{Global}\_\mathrm{ij}}^{(1)}+{\lambda }_{\mathrm{Global}\_\mathrm{ij}}^{(2)}$$

(5)

where the global compliance matrix of each contact position can be obtained by finite element method and substructure method. ²⁶ Considering the nonlinear relationship between local contact deformation and the applied force, the Hertzian contact deflection of interested contact point can be calculated by ²⁷

$${\delta }_{\mathrm{Contact}\_i}=\frac{4F_i}{\mathrm{dz}}\frac{1-{\nu }^2}{\pi E}[\ln \frac{2\sqrt{k_1{k}_2}}{a}-\frac{\nu }{2(1-\nu )}]$$

(6)

where F_i is the applied force and dz refers to the gear width. k₁ and k₂ are the distances on the pinion and the gear between the contact point and the tooth center line along the direction of applied force; E and ν are Young’s modulus and Poisson’s ratio, respectively; and a refers to the half width of contact zone along the profile direction, which can be determined by

$$a=\sqrt{\frac{8F/b{\rho }_1{\rho }_2}{({\rho }_1+{\rho }_2)(1-{\nu }^2)/\pi E}}$$

(7)

where ρ₁ and ρ₂ are the curvature radius of the pinion and gear, respectively.

The deformation compatibility condition of contact points can be written as

$$-\sum _{j=1}^n{\lambda }_{\mathrm{Global}\_\mathrm{ij}}{F}_j-{\delta }_{\mathrm{Contact}\_i}+\mathrm{LSTE}+d_i={\epsilon }_i$$

(8)

Each contact point can be coupled using LSTE, the contact equation of each engagement position will be written as

$$\begin{array}{c}-\left[\begin{array}{cccc}{\lambda }_{\mathrm{Global}\_11}\hfill & {\lambda }_{\mathrm{Global}\_12}\hfill & \cdots \hfill & {\lambda }_{\mathrm{Global}\_1n}\hfill \\ \multicolumn{1}{c}{{\lambda }_{\mathrm{Global}\_21}}& {\lambda }_{\mathrm{Global}\_22}\hfill & \cdots \hfill & {\lambda }_{\mathrm{Global}\_1n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\lambda }_{\mathrm{Global}\_n1}}& {\lambda }_{\mathrm{Global}\_n2}\hfill & \cdots \hfill & {\lambda }_{\mathrm{Global}\_\mathrm{nn}}\hfill \end{array}\right]\left\{\begin{array}{c}{F}_1\hfill \\ \multicolumn{1}{c}{F_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{F_n}\end{array}\right\}\hfill \\ -\left\{\begin{array}{c}{\delta }_{\mathrm{Contact}\_1}\hfill \\ \multicolumn{1}{c}{{\delta }_{\mathrm{Contact}\_2}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\delta }_{\mathrm{Contact}\_n}}\end{array}\right\}+\mathrm{LSTE}+\left\{\begin{array}{c}{d}_1\hfill \\ \multicolumn{1}{c}{d_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{d_n}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_1\hfill \\ \multicolumn{1}{c}{{\epsilon }_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\epsilon }_n}\end{array}\right\}\hfill \end{array}$$

(9)

In above, the matrix form of loaded tooth contact equation can be written as

$$\{\begin{array}{c}-{\left[\lambda \right]}_{\mathrm{Global}}\left\{F\right\}-{\left\{u\right\}}_{\mathrm{Local}}+\mathrm{LSTE}+\left\{d\right\}=\left\{\epsilon \right\}\hfill \\ \sum _{i=1}^n{F}_i=\left\{I\right\}\left\{F\right\}=P\hfill \\ \mathrm{if}{F}_i>0,d_i=0\hfill \\ \mathrm{if}{F}_i=0,d_i>0\hfill \end{array}$$

(10)

where [λ] _Global is the compliance matrix of global deflection of contact point, {F} is the external load vector, {u} _Local is the contact deflection of interested contact point, {d} is the odd clearance vector of contact point, and {ε} is the initial clearance vector of contact point.

The load distribution F and loaded static transmission error LSTE can be obtained using the iteration algorithm ²⁶ to solve the nonlinear matrix equations.

The time-varying mesh stiffness of helical gears with certain errors can be written as

$$k_m=\sum _{i=1}^N{k}_i=\sum _{i=1}^N\frac{F_i}{\mathrm{LSTE}-{\epsilon }_i}$$

(11)

The static load balance equation can be written as

$$\begin{array}{c}\sum _{i=1}^N{k}_i(\mathrm{LSTE}-{\epsilon }_i)=\sum _{i=1}^N{k}_i\mathrm{LSTE}-\sum _{i=1}^N{k}_i{\epsilon }_i\hfill \\ =k_m\mathrm{LSTE}-\sum _{i=1}^N{k}_i{\epsilon }_i=P\hfill \end{array}$$

(12)

Here, the loaded composite mesh error of helical gears can be defined as

$$e_m=\frac{1}{k_m}\sum _{i=1}^N{k}_i{\epsilon }_i$$

(13)

Therefore, the static load balance equation (12) can be rewritten as

$$k_m(\mathrm{LSTE}-e_m)=P$$

(14)

Obviously, the loaded composite mesh error is related to mesh stiffness, distribution form of gear errors, and the applied force.

Dynamic model of helical gear system

A simplified illustration of helical gear system and the corresponding generalized finite element model, which consists of three type nodes, namely, shaft node, bearing node, and power node, are shown in Figure 4. The parameters of each shaft are listed Table 1. In order to consider the effects of structure of shaft, the direction of power flow, and the supporting location of bearing on the dynamic responses of system, the helical gear system is divided into three elements: shaft element, mesh element, and bearing element. First, the dynamic model of shaft element is established by Timoshenko beam element. Second, the dynamic model of mesh element is developed by spring–damping–error element. Here, the stiffness of spring denotes time-varying mesh stiffness, and the error is loaded composite mesh error. In addition, the dynamic model of bearing element can be established by the stiffness and damping matrix theory of bearing. Finally, dynamic model of the system can be obtained by assembling dynamic model of each element based on the finite element method.

Figure 4.

Schematic diagram of helical gear system.

Table 1.

Basic parameters of each shaft.

Parameters of the shaft (mm)
Number of shaft segment	Shaft 1			Shaft 2
	Outer diameter	Inner diameter	Length	Outer diameter	Inner diameter	Length
1	70	0	100	90	0	80
2	90	0	110	110	0	110
3	100	0	260	120	0	55
4	90	0	110	160	0	150
5	70	0	60	120	0	55
6				110	0	110
7				90	0	120

Modeling of shaft element

The dynamic model of shaft element is established by Timoshenko beam element with 2 nodes and 12 degrees of freedom, as shown in Figure 5. The generalized coordinate of Timoshenko beam element is defined as

$$q_s={\left\{x_i,y_i,z_i,{\theta }_{xi},{\theta }_{yi},{\theta }_{zi},x_{i+1},y_{i+1},z_{i+1},{\theta }_{x(i+1)},{\theta }_{y(i+1)},{\theta }_{z(i+1)}\right\}}^T$$

(15)

where x_i, y_i, z_i, x_i+1, y_i+1, and z_i+1 denote the translational displacement of each node along coordinate axes, and θ_xi, θ_yi, θ_zi, θ_x(i+1), θ_y(i+1), and θ_z(i+1) are the rotational displacement of each node around each coordinate axis.

Figure 5.

Schematic diagram of shaft element.

Damping matrix of shaft element C_s can be calculated by Rayleigh damping

$$C_s={\alpha }_0{M}_s+{\alpha }_1{K}_s$$

(16)

where α₀ and α₁ denote the mass coefficient and stiffness coefficient, respectively, and M_s and K_s are the mass matrix and stiffness matrix of shaft element, respectively.

When the mass matrix, stiffness matrix, and damping matrix are obtained, the motion differential equation in matrix form for shaft element can be written as

$$M_s{\stackrel{··}{q}}_s\left(t\right)+C_s{\stackrel{·}{q}}_s\left(t\right)+K_s{q}_s\left(t\right)=0$$

(17)

where stiffness matrix K_s and mass matrix M_s can be obtained based on basic theory of elastic mechanics. ²⁸

Modeling of mesh element

The dynamic model of mesh element with 12 degrees of freedom is modeled by spring–damping–error element, as shown in Figure 6. The pinion and wheel are considered as two rigid disks that are connected by spring–damping at base circle, where β is the base helix angle; O_p and O_g denote the running center of pinion and wheel, respectively; r_p and r_g refer to the base circle of pinion and wheel; α is the mesh angle; ψ is the installation phase angle; x_i, y_i, z_i(i = p, g) denote the translational displacement of pinion and wheel along each coordinate axis; and θ_xi, θ_yi, θ_zi(i = p, g) denote the rotational displacement of pinion and wheel rotating around each coordinate axis (Table 2).

Figure 6.

Dynamic model of mesh element.

Table 2.

Basic parameters of the pinion and wheel.

Parameters	Pinion	Wheel
Tooth number	29	58
Module	5	5
Normal pressure angle (°)	20	20
Helix angle (°)	25	25
Addendum coefficient	1	1
Bottom clearance coefficient	0.25	0.25
Tooth width (mm)	50	50
Rotational speed (r/min)	2000	1000
Applied torque (N m)	750	1500

The generalized coordinate vector of mesh element is defined as

$$q_m={\{x_p,y_p,z_p,{\theta }_{\mathrm{xp}},{\theta }_{\mathrm{yp}},{\theta }_{\mathrm{zp}},x_g,y_g,z_g,{\theta }_{\mathrm{xg}},{\theta }_{\mathrm{yg}},{\theta }_{\mathrm{zg}}\}}^T$$

(18)

The relative displacement of pinion and wheel along the line of action can be expressed as

$$\tau =V{q}_m$$

(19)

where the vector V refers to the projective vector from the coordinate of gear node to the line of action, which can be given as

$$\begin{array}{c}\multicolumn{1}{c}{V=[\cos \beta \sin \phi ,\pm \cos \beta \cos \phi ,\mp \sin \beta ,}\\ r_p\sin \beta \sin \phi ,\pm r_p\sin \beta \cos \phi ,\pm r_p\cos \beta ,\hfill \\ -\cos \beta \sin \phi ,\mp \cos \beta \cos \phi ,\pm \sin \beta ,\hfill \\ r_g\sin \beta \sin \phi ,\pm r_g\sin \beta \cos \phi ,\pm r_g\cos \beta \hfill \end{array}]$$

(20)

where the plus sign in “±” and minus sign in “∓” are for anticlockwise rotation of pinion and the minus sign in “±” and plus sign in “∓” are for clockwise rotation of pinion.

Considering time-varying mesh stiffness and loaded composite mesh error, the motion differential equation of mesh element can be given as

$$\{\begin{array}{c}{m}_p{\stackrel{··}{x}}_p+[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]\cos \beta \sin \phi =0\hfill \\ \multicolumn{1}{c}{m_p{\stackrel{··}{y}}_p+[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]\cos \beta \cos \phi =0}\\ \multicolumn{1}{c}{m_p{\stackrel{··}{z}}_p+[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]\sin \beta =0}\\ \multicolumn{1}{c}{I_{\mathrm{xp}}{\stackrel{··}{\theta }}_{\mathrm{xp}}-[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]r_p\sin \beta \sin \phi =0}\\ \multicolumn{1}{c}{I_{\mathrm{yp}}{\stackrel{··}{\theta }}_{\mathrm{yp}}-[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]r_p\sin \beta \cos \phi =0}\\ \multicolumn{1}{c}{I_{\mathrm{zp}}{\stackrel{··}{\theta }}_{\mathrm{zp}}+[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]r_p\cos \beta =0}\\ \multicolumn{1}{c}{m_g{\stackrel{··}{x}}_g-[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]\cos \beta \sin \phi =0}\\ \multicolumn{1}{c}{m_g{\stackrel{··}{y}}_g-[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]\cos \beta \cos \phi =0}\\ \multicolumn{1}{c}{m_g{\stackrel{··}{z}}_g-[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]\sin \beta =0}\\ \multicolumn{1}{c}{I_{\mathrm{xg}}{\stackrel{··}{\theta }}_{\mathrm{xg}}-[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]r_p\sin \beta \sin \phi =0}\\ \multicolumn{1}{c}{I_{\mathrm{yg}}{\stackrel{··}{\theta }}_{\mathrm{yg}}-[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]r_p\sin \beta \cos \phi =0}\\ \multicolumn{1}{c}{I_{\mathrm{zg}}{\stackrel{··}{\theta }}_{\mathrm{zg}}+[c_m\stackrel{·}{\tau }+k_m(\tau -e_m)]r_p\cos \beta =0}\end{array}$$

(21)

where $\phi =\alpha \mp \psi$ , the minus sign in “∓” denotes anticlockwise rotation of pinion, and the plus sign in “∓” refers to clockwise rotation of pinion.

The time-varying stiffness matrix K_m and damping matrix C_m of mesh element can be expressed as

$$K_m=k_m{V}^{T}V$$

(22)

$$C_m=c_m{V}^{T}V$$

(23)

The motion differential equation in matrix form for mesh element can be written as

$$M_m{\stackrel{··}{q}}_m\left(t\right)+C_m{\stackrel{·}{q}}_m\left(t\right)+K_m\left(t\right)[q_m\left(t\right)-e\left(t\right)]=0$$

(24)

where M_m is the mass matrix, e(t) is the loaded composite mesh error vector of each degree of freedom, and q_m is the displacement vector of shaft node connected with gears.

Modeling of bearing element

Gear system is usually supported by gear housing and connected to foundation, ignoring the coupling effect between the flexibility of gear housing and gear system, and the bearing element can be considered as a spring–damping element which connects shaft node with foundation. The general stiffness matrix of bearing K_b can be given as

$$K_b=\left[\begin{array}{cccccc}{k}_{\mathrm{xx}}\hfill & k_{\mathrm{xy}}\hfill & k_{\mathrm{xz}}\hfill & k_{x{\theta }_x}\hfill & k_{x{\theta }_y}\hfill & 0\hfill \\ \multicolumn{1}{c}{k_{\mathrm{yx}}}& k_{\mathrm{yy}}\hfill & k_{\mathrm{yz}}\hfill & k_{y{\theta }_x}\hfill & k_{y{\theta }_y}\hfill & 0\hfill \\ \multicolumn{1}{c}{k_{\mathrm{zx}}}& k_{\mathrm{zy}}\hfill & k_{\mathrm{zz}}\hfill & k_{z{\theta }_x}\hfill & k_{z{\theta }_y}\hfill & 0\hfill \\ \multicolumn{1}{c}{k_{{\theta }_{x}x}}& k_{{\theta }_{x}y}\hfill & k_{{\theta }_{x}z}\hfill & k_{{\theta }_x{\theta }_x}\hfill & k_{{\theta }_x{\theta }_y}\hfill & 0\hfill \\ \multicolumn{1}{c}{k_{{\theta }_{y}x}}& k_{{\theta }_{y}y}\hfill & k_{{\theta }_{y}z}\hfill & k_{{\theta }_y{\theta }_x}\hfill & k_{{\theta }_y{\theta }_y}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(25)

Actually, the stiffness of each degree of freedom of bearing is time varying and coupled with each other. Here, ignoring the time-varying characteristic of bearing stiffness and the coupling relationship of different degrees of freedom, namely, the nonmain diagonal element of stiffness matrix, is zero.

The motion differential equation in matrix form for bearing element can be written as

$$M_b{\stackrel{··}{q}}_b\left(t\right)+C_b{\stackrel{·}{q}}_b\left(t\right)+K_b{q}_b\left(t\right)=0$$

(26)

where M_b denotes the mass matrix of bearing, C_b is the damping matrix of bearing which is identical with stiffness matrix of bearing, and q_b is the displacement vector of shaft node connected with bearings.

Modeling of gear system

Based on the basic theory of finite element method, the global dynamic model of the system can be obtained by assembling dynamic model of each element. The global motion differential equation in matrix form for the system can be written as

$$M\stackrel{··}{x}\left(t\right)+C\stackrel{·}{x}\left(t\right)+K\left(t\right)[x\left(t\right)-e\left(t\right)]=F$$

(27)

where M is the global mass matrix; C is the global damping matrix; K(t) is the global stiffness matrix; x(t) is the generalized coordinates of finite element nodes; e(t) refers to the vector of loaded composite mesh error, which is equal to zero at each degree of freedom except degrees of freedom of mesh element; and F is the external force vector.

The static balance equation of the system can be written as

$$K\left(t\right)[x_s\left(t\right)-e\left(t\right)]=F$$

(28)

where x_s(t) is the static displacement vector considering the deformation of shaft and bearing.

Vibration exciting force

Approximate transformation

The matrix form of motion differential equation (27) for the system can be rewritten as

$$M\stackrel{··}{x}\left(t\right)+C\stackrel{·}{x}\left(t\right)+K\left(t\right)x\left(t\right)=P_0+K\left(t\right)e\left(t\right)$$

(29)

Denote K(t) as the summation of mean value and fluctuation, which is

$$K\left(t\right)=K_0+\Delta K\left(t\right)$$

(30)

Substituting equation (30) into equation (29) and combining with equation (28) yields

$$M\stackrel{··}{x}\left(t\right)+C\stackrel{·}{x}\left(t\right)+K_{0}x\left(t\right)=[K_0+\Delta K\left(t\right)]x_s\left(t\right)-\Delta K\left(t\right)x\left(t\right)$$

(31)

As x(t) is not known yet, and the difference between x(t) and x_s(t) is small enough to be ignored for helical gears, replace x(t) on the right side of the equation by x_s(t). Then, the motion differential equation for the system can be written as

$$M\stackrel{··}{x}\left(t\right)+C\stackrel{·}{x}\left(t\right)+K_{0}x\left(t\right)=K_0{x}_s\left(t\right)$$

(32)

With the deduction above, the original time-variant differential equation (29) has been transformed into time-invariant differential equation (32). The right side of the equation K₀x_s(t) is time varying. Obviously, the K₀x_s(t) determines the vibration degree of gear transmission. Here, the K₀x_s(t) is defined as vibration exciting force.

The dynamic displacement x(t) can be expressed as Fourier series as

$$\begin{array}{c}x\left(t\right)=x_0+\sum _{i=1}^n[A_i\cos i{\omega }_{m}t+B_i\sin i{\omega }_{m}t]\hfill \\ \stackrel{·}{x}\left(t\right)=\sum _{i=1}^n[-i{\omega }_i{A}_i\sin i{\omega }_{m}t+i{\omega }_i{B}_i\cos i{\omega }_{m}t]\hfill \\ \stackrel{··}{x}\left(t\right)=\sum _{i=1}^n[-{\left(i{\omega }_i\right)}^2{A}_i\cos i{\omega }_{m}t-{\left(i{\omega }_i\right)}^2{B}_i\sin i{\omega }_{m}t]\hfill \end{array}$$

(33)

where ω_m is the mesh frequency and n denotes the maximum order of Fourier series.

As same as x(t), the vibration exciting force K₀x_s(t) can be denoted as Fourier series as

$$K_0{x}_s\left(t\right)=F_0+\sum _{i=1}^n[C_i\cos i{\omega }_{m}t+D_i\sin i{\omega }_{m}t]$$

(34)

Substituting equations (33) and (34) into equation (32) yields linear equations as

$$\left[\begin{array}{cc}-{\omega }_k^{2}M+K_0\hfill & {\omega }_{k}C\hfill \\ \multicolumn{1}{c}{-{\omega }_{k}C}& -{\omega }_k^{2}M+K_0\hfill \end{array}\right]\left[\begin{array}{c}{A}_i\hfill \\ \multicolumn{1}{c}{B_i}\end{array}\right]=\left[\begin{array}{c}{C}_i\hfill \\ \multicolumn{1}{c}{D_i}\end{array}\right]$$

(35)

where ω_k = iω_m (i = 1, 2, …, n). Substituting A_i and B_i into equation (33) after solving equations (35), the dynamic displacement x(t) of each node can be obtained.

The dynamic transmission error of the system can be determined as

$$\mathrm{DTE}=V{q}_m$$

(36)

The dynamic bearing force of the system can be calculated as

$$\{\begin{array}{c}\mathrm{DB}{F}_x=k_b{x}_b+c_b{\stackrel{·}{x}}_b\hfill \\ \multicolumn{1}{c}{\mathrm{DB}{F}_y=k_b{y}_b+c_b{\stackrel{·}{y}}_b}\\ \multicolumn{1}{c}{\mathrm{DB}{F}_z=k_r{z}_b+c_r{\stackrel{·}{z}}_b}\end{array}$$

(37)

where DBF_x, DBF_y, and DBF_z are dynamic bearing forces in the x, y, and z directions, respectively. k_b and k_r denote the radial stiffness and tangential stiffness of bearing. c_b and c_r are the radial damping and tangential damping of bearing. x_b, y_b, and z_b are dynamic displacements of bearing node along corresponding coordinate axis, respectively.

Comparison of solution methods

With the analysis and deduction above, the origin time-variant differential equation (29) is transformed into time-variant differential equation (32) by approximate transformation in order to improve calculation efficiency by Fourier series method. In order to validate the accuracy of approximate transformation, the dynamic responses of helical gear system with 10 µm of parabolic tip relief in three-tooth engagement zone are calculated by Newmark-β method and Fourier series method, respectively.

The root mean square (RMS) values of dynamic transmission error of the system in different speeds of pinion which are solved by the two methods are shown in Figure 7. It can be seen that the calculation result of the proposed method is very close to that of Newmark-β method, and the maximum relative error in all speeds of pinion is less than 4.2%. It indicates that the approximate transformation is feasible, and the precision of calculation is ensured perfectly.

Figure 7.

Comparison of Newmark-β method and proposed method.

Numerical results and discussion

Global brute force optimization

Traditional optimization algorithm has the problem of local convergence, and the effect of design variables on the trend of the objective function cannot be captured. However, the improved LTCA model is of high efficiency because the global deformation compliance of interested tooth surface needs to be calculated only once. Therefore, the optimal solution of bias modification can be obtained by the direct calculation of the global sweep in the potential range of the modification parameters.

It can be known that the degree of vibration of gear system is determined by vibration exciting force from section “Vibration exciting force.” Thus, the optimization method of bias modification parameters aiming to minimize the fluctuation of vibration exciting force is proposed in order to reduce vibration and noise of the system. The region of amount of modification is from 0 to Ca, and the region of length of modification is from 0 to La. Here, Ca denotes two times of the deformation of helical gears in mesh-in position. The schematic diagram of optimization of bias modification for minimizing the fluctuation of vibration exciting force is shown in Figure 8. It can be seen that there is a narrow zone in which each modification parameter can reduce the fluctuation of vibration exciting force significantly. It is easy to select the optimal bias modification parameters for the helical gear system from Figure 8. The optimal amount of modification is 0.85C_a, namely, 10.1 µm. The optimal length of modification is 0.75L_a, namely, 17.1 mm.

Figure 8.

Optimization of bias modification for minimizing the fluctuation of vibration exciting force.

The time-varying mesh stiffness and loaded composite mesh error of helical gears with bias modification are shown in Figure 9, and the corresponding load distributions are displayed in Figure 10. It can be observed that the load distributions of modified helical gears at mesh-in and mesh-out positions are smaller than that of unmodified helical gears but not equal to zero. That is to say, the tooth surface is still under fully contact condition. Therefore, the time-varying mesh stiffness curve of helical gears with bias modification is very close to that of unmodified gears. With the nonlinearity of contact deformation and change of load distribution, the maximum reduction of time-varying mesh stiffness in a mesh cycle is less than 0.8%. It is also observed that the shape of loaded composite mesh error is consistent with that of time-varying mesh stiffness, and the fluctuation of loaded composite mesh error 1.4 µm is far less than the maximum amount of bias modification, which is 10.1 µm. This is because loaded composite mesh error is not only related to amount of bias modification but also to the number of mating gear teeth and loaded tooth contact condition.

Figure 9.

Effect of bias modification on time-varying mesh stiffness and loaded composite mesh error.

Figure 10.

Effect of bias modification on load distribution of helical gears: (a) unmodified and (b) modified.

Dynamic responses analysis

The dynamic transmission error of the helical gear system and dynamic bearing force of bearing 2 in operating speed of pinion with 2000 r/min are shown in Figure 11. The results show that bias modification has significant effect on reduction of vibration of the system. As a result of intentional gear errors introduced by bias modification, the mean value of dynamic transmission error increases from 11.3 to 12.4 µm, about 9.8%, but the fluctuation of that decreases from 3.1 to 0.05 µm, about 98.5%. Meanwhile, the fluctuations of dynamic bearing force along the x and z directions reduce about 98.8% and 97.7%, respectively.

Figure 11.

Effect of bias modification on dynamic responses of helical gear system: (a) dynamic transmission error, (b) dynamic bearing force in the x direction, and (c) dynamic bearing force in the z direction.

Figure 12 shows the RMS values of dynamic transmission error in regions of speeds of pinion from 0 to 6000 r/min. It can be observed that the torsional resonance speed of the system is about 4000 r/min. Harmonic resonances of the system occur in many times when the speed of pinion is less than 4000 r/min, and the corresponding amplitudes of harmonic resonance are much less than that of torsional resonance. The RMS values of dynamic transmission error of helical gear system will reduce remarkably in all speeds of pinion owing to optimal bias modification.

Figure 12.

RMS of dynamic transmission error of helical gear system at different speeds of pinion.

Sensitivity to gear misalignment

It is well known that gear misalignment is a common and inevitable problem in gear transmissions that has an important influence on the contact condition of mating gear teeth and dynamic response of gear system. Some gear misalignments are caused directly by manufacturing and assembling errors, but others may be mainly related to bending and torsional deformation of shafts, deformation of gearbox during power transmission, or bearing clearances. It is obvious that the cause of gear misalignment is complicated, and it is difficult to predict and measure precisely.

Although the cause of gear misalignment is various, it can be divided into three categories: parallel misalignment due to center distance change, angular misalignment parallel to the plane of action (LOA misalignment), and angular misalignment perpendicular to the plane of action (OLOA misalignment).^25,29 LOA misalignment refers to relative angular displacement between the meshing tooth surfaces of pinion and wheel around an axis which is perpendicular to the plane of action, the effect of which on contact condition of mating gear teeth and dynamic response of gear system is most significant. Therefore, only LOA misalignment is considered in this section, and it can be modeled as helix deviation,^30,31 as shown in Figure 13.

Figure 13.

Definition of gear misalignment.

In total, four cases of gear misalignment which is from 0 to 36 µm with step of 12 µm are designed in order to investigate the effect of gear misalignment on quasi-static and dynamic behaviors of modified gear system. The effects of gear misalignment on time-varying mesh stiffness and loaded composite mesh error of modified helical gears are shown in Figure 14. The results show that both time-varying mesh stiffness and loaded composite mesh error curves are changed due to gear misalignment. With the increase in gear misalignment, the time-varying mesh stiffness will decrease, but loaded composite mesh error will increase notably. The mean values of time-varying mesh stiffness reduce about 3.8%, 17.9%, and 30.5% for 12, 24, and 36 µm gear misalignment, respectively, and the corresponding mean values of loaded composite mesh error increase about 478.1%, 788.1%, and 982.1%.

Figure 14.

Effect of gear misalignment on time-varying mesh stiffness and loaded composite mesh error (the corresponding straight dot line denotes mean value of time-varying mesh stiffness).

The effect of gear misalignment on the fluctuation of vibration exciting force is shown in Figure 15. It can be seen that the fluctuations of vibration exciting force increase nonlinearly with the increase in gear misalignment. The fluctuations of vibration exciting force of unmodified helical gear system increase about 4.1%, 7.9%, and 40.4% for 12, 24, and 36 µm gear misalignment. For modified helical gear system, the fluctuations of vibration exciting force increase 951.7%, 2824.4%, and 5122.7%, respectively. Compared with unmodified helical gear system, the fluctuations of vibration exciting force of modified helical gear system reduce 98.0%, 79.5%, and 45.0% for 0, 12, and 24 µm gear misalignment. When gear misalignment is 36 µm, the reduction is only 24.5%. This is because gear misalignment will change the load distribution of mating gear teeth, as shown in Figure 16. It can be observed that increasing gear misalignment will lead to appearance of a partial contact loss through the face width. When gear misalignment is larger than 24 µm, the actual length of contact line of mating gear teeth will reduce remarkably, and the contact loss will become very serious.

Figure 15.

Effect of gear misalignment on the fluctuation of vibration exciting force.

Figure 16.

Effect of gear misalignment on load distribution of modified helical gears: (a) without gear misalignment, (b) 12 µm gear misalignment, (c) 24 µm gear misalignment, and (d) 36 µm gear misalignment.

The effects of gear misalignment on dynamic transmission error of helical gear system with and without bias modification are displayed in Figure 17. It can be seen that increase in gear misalignment will bring larger vibration of both unmodified and modified helical gear systems. Due to decreased mesh stiffness, the lower resonance speed of pinion will become lower with the increase in gear misalignment. Bias modification has a well reduction of vibration and noise of system when gear misalignment is relatively small, namely, less than 12 µm. The corresponding RMS value of dynamic transmission error reduces more than 80% in operating speed of pinion. The bias modification has little effect on reduction of system vibration when gear misalignment is larger than 32 µm. In operating speed of pinion, the corresponding RMS value of dynamic transmission error decreases only less than 17%.

Figure 17.

RMS of dynamic transmission error of helical gear system.

Conclusion

The influence of optimal bias modification on quasi-static and dynamic behaviors of helical gear system and the sensitivity of vibration of the system with optimal bias modification to gear misalignment have been studied in this article. The dynamic model of the system is developed using the finite element method and Timoshenko beam theory, and the vibration exciting force is separated through approximate transformation which was validated by comparing with numerical integration method. Aiming to minimize the fluctuation of vibration exciting force, the optimal modification parameters are determined using global brute force optimization method based on improved LTCA model. The effects of optimal bias modification on time-varying mesh stiffness, loaded composite mesh error, and load distribution of helical gear system with and without gear misalignment are investigated. The influences of gear misalignment on the effectiveness of reduction of vibration of helical gear system with optimal bias modification are also analyzed. The main conclusions are summarized as follows:

1. The time-variant differential equations of dynamic model of helical gear system are transformed into time-invariant differential equations which can be solved faster by Fourier series method. The approximate transformation provides sufficient precision and greatly increases the computation efficiency. Meanwhile, the internal excitation of helical gear system, defined as vibration exciting force, can be observed clearly.

2. Optimal bias modification has little effect on time-varying mesh stiffness of helical gears. Owing to full contact of tooth surface of helical gears with optimal bias modification, the time-varying mesh stiffness curve of modified helical gears is very identical with that of unmodified ones. The corresponding RMS values of dynamic transmission error and fluctuation of dynamic bearing force will decrease remarkably.

3. Gear misalignment will increase the vibration of helical gear system and result in lower resonance speed because of decreased mesh stiffness. Bias modification shows a well reduction of vibration of helical gear system when gear misalignment is relatively small. When gear misalignment is too large, more serious contact loss will be observed, and robust optimization should be considered.