Dynamic characteristic analysis of spur gear system considering tooth contact state caused by shaft misalignment

The effect of gear contact state change due to shaft misalignment on meshing stiffness is usually neglected in the traditional stiffness calculation model with misalignment error, the further influence mechanism of shaft misalignment on gear dynamic characteristics is also unclear. To address these shortcomings, a new mesh stiffness calculation model with misaligned gear considering the effects of tooth contact state is proposed by combining the improved loaded tooth contact analysis (LTCA) model. Then the effects of tooth contact state changes aroused by shaft misalignment on the meshing stiffness excitation are investigated. Moreover, a dynamic model of the misaligned gear system with 8 degrees of freedom (DOF) is established, and based on which the dynamic characteristics of the gear system are investigated and verified by experiment. The study results indicate that the proposed model can be used to evaluate the stiffness excitation and dynamic characteristics of the misaligned gear system with the tooth contact state taken into consideration. This study can provide a theoretical method for evaluating and identifying shaft misalignment errors.

Abstract The effect of gear contact state change due to shaft misalignment on meshing stiffness is usually neglected in the traditional stiffness calculation model with misalignment error, the further influence mechanism of shaft misalignment on gear dynamic characteristics is also unclear. To address these shortcomings, a new mesh stiffness calculation model with misaligned gear considering the effects of tooth contact state is proposed by combining the improved loaded tooth contact analysis (LTCA) model. Then the effects of tooth contact state changes aroused by shaft misalignment on the meshing stiffness excitation are investigated. Moreover, a dynamic model of the misaligned gear system with 8 degrees of freedom (DOF) is established, and based on which the dynamic characteristics of the gear system are investigated and verified by experiment. The study results indicate that the proposed model can be used to evaluate the stiffness excitation and dynamic characteristics of the misaligned gear system with the tooth contact state taken into consideration. This study can provide a theoretical method for evaluating and identifying shaft misalignment errors.

List of symbols a
Gear central angle corresponding to the micro-section at x from the dedendum circle a 1 Gear central angle corresponding to the force point a r /a q Gear central angle corresponding to the force point on the tooth #r/#q (r = 1, 2; q = 1, 2, r = q) a 3 Half of the tooth angle on base circle Cross-sectional area of the aligned/ misaligned gear A Deformation influence coefficient matrix A c Deformation influence coefficient matrix of tooth pair #1 (c = 1) or #2 (c = 2) c kx /c ky /c kz Bearing damping of driving gear (k = 1) or driven gear (k = 2) in x/y/z direction d Distance between tooth force point and dedendum circle Clearance matrix corresponding to unit grid of common tangent plane of tooth pair #1 (c = 1) or #2 (c = 2) H fr Clearance matrix of tooth #r (r = 1, 2) caused by the fillet foundation deformation under the meshing force on the tooth #r (r = 1, 2) H frq Clearance matrix of tooth #r caused by fillet foundation deformation under the meshing force on the adjacent tooth #q (r = 1, 2; q = 1, 2, r = q) I Unit vector matrix I x /I x ' Section moment of inertia of the aligned/misaligned gear tooth I px Polar moment of inertia of the microsection J k Rotational inertia of driving gear (k = 1) or driven gear (k = 2) k m Mesh stiffness k h /k b /k s /k a / k f /k t Stiffness of contact/bending/shear/axial compression/fillet foundation/torsion k kx /k ky /k kz Bearing stiffness of driving gear (k = 1) or driven gear (k = 2) in x/y/z direction L 0 Effective contact line length L eq Distance between applied point of equivalent force and gear tooth center m k Mass of driving gear (k = 1) or driven gear (k = 2) n Number of unit grid in tooth common tangent plane p j Applied force at unit grid j P Pressure matrix of common tangent plane of contact tooth pair P c Pressure matrix of common tangent plane of tooth pair #1 (c = 1) or #2 (c = 2) r b Base circle radius of gear r bk Base circle radius of driving gear (k = 1) or driven gear (k = 2) s Area of unit grid of tooth common tangent plane T Torsional torque of misaligned gear tooth T k Torque of driving gear (k = 1) or driven gear (k = 2) U b /U s /U a / U t Deformation energy of bending/shear/ axial compression/torsion w k Deformation of unit grid node k W Tooth width x k = _ x k = € x k Lateral displacement/velocity/ acceleration of driving gear (k = 1) or driven gear (k = 2) along x direction y k = _ y k = € y k Lateral displacement/velocity/ acceleration of driving gear (k = 1) or driven gear (k = 2) along y direction z k = _ z k =€ z k Lateral displacement/velocity/ acceleration of driving gear (k = 1) or driven gear (k = 2) along z direction h k = _ h k = € h k Angular displacement/velocity/ acceleration of driving gear (k = 1) or driven gear (k = 2) d mis Comprehensive deformation of gear Dh x /Dh y Misalignment error in x/y direction t Poisson's ratio

Introduction
Due to the advantages of high transmission efficiency, compact structure and high reliability [1,2], gear transmission system is widely applied in the modern industry, e.g., wind turbine, helicopter, automobile [3,4]. Shaft misalignment error is one of the most common errors in the gear transmission system, which can lead to variation of tooth surface contact state, including the change of contact stress distribution and effective contact area. This variation is not conducive to tooth surface lubrication and will accelerate surface wear of gear teeth [5]. In addition, shaft misalignment will also cause changes in the internal excitation of gear pairs, thereby increasing vibration and noise level and reducing the service life of gear system [6]. Therefore, studying the influences of gear contact state change caused by shaft misalignment on the internal excitation and dynamic responses of gear system is of great significance for monitoring misalignment errors and improving gear transmission life. Tooth contact analysis model is a common method adopted by many scholars in the analyzation of gear contact state [7,8]. For instance, Hartnett [9,10] established the cylinder contact analysis calculation model using the Boussinesq half-space force-displacement relations and verified the effectiveness of the proposed model through experiment. Based on the model proposed by Hartnett, Ahmadi [11] developed an analysis and calculation model for the tooth surface contact state of spur gear. Then, Kolivand and Kahraman [12] improved Ahmadi's model by considering the gear tooth profile surface deviations, and then analyzed the contact characteristic of the cone tooth surface. Considering the influences of tooth contact deformation and tooth bending deflection, Wu et al. [13] presented an LTCA model of skew conical gear drives in approximate line contact. Chao et al. [14] put forward an LTCA model for turbine screw by combining the gear space meshing theory and the curvature change law of gear space surface. Taking into account both the deformation of loaded teeth and the twist deformation of gear shaft, Ye et al. [15] proposed an LTCA model to investigate the contact state of high-contact-ratio spur gears based on the influence coefficient method. Although LTCA model and its improved versions see their success in lots of studies, the fillet foundation deformation of the singletooth and the coupling deformation between the neighboring teeth are usually ignored in the traditional LTCA model, while the fillet foundation deformations have an effect on the tooth contact state of the meshing gear [15].
It is known that the shaft misalignment can change the gear internal excitation, and thus the dynamic characteristics of the gear system will be affected [16]. Some scholars have investigated the influence of shaft misalignment on the internal excitation and dynamic characteristics of gear system through establishing theoretical models. For example, Wang et al. [17] analyzed the influence of offset load caused by shaft misalignment on gear mesh stiffness and vibration responses based on the finite element models. Considering the impact of eccentricity due to the shaft misalignment, Cao et al. [18] established a dynamic model of planetary gear train, and investigated the effects of shaft misalignment on the gear vibration and meshing force. Saxena et al. [19] considered the influence of the direction change of meshing force caused by misalignment error, and improved the calculation model of mesh stiffness of misaligned gear based on the energy method. Xu et al. [20] established a three-dimensional dynamic model of misaligned spur gear, and discussed the effect of shaft misalignment on gear dynamic responses. Jiang et al. [21] investigated the dynamic characteristics of misaligned helical gears including sliding friction. Above all, the influence of tooth contact state caused by shaft misalignment error on the stiffness excitation is neglected in most of the existing stiffness calculation models, and the influence mechanism of the contact state variation aroused by shaft misalignment on the gear dynamic characteristics is still unclear.
As can be seen from the above review, there are two shortcomings of the existing simulation models with shaft misalignment: (1) The influence of fillet foundation deformation on gear contact state is ignored in the traditional LTCA model; (2) The effect of tooth contact state change on gear mesh stiffness is not considered in the traditional stiffness calculation model of misaligned gears. And the underlying influence mechanism of the gear contact state variation aroused by shaft misalignment on the dynamic characteristics of gear system needs to be further studied.
Therefore, the aim of this current paper is to realize the accurate calculation of the tooth contact state of misaligned gear, then propose a mesh stiffness calculation model that can consider the effects of tooth contact state affected by shaft misalignment, and finally reveal the influence mechanism of the tooth contact state variation caused by misalignment error on the dynamic characteristics of gear system. The detailed research contents of this paper are as follows: (1) An improved LTCA model is proposed by considering the influences of fillet foundation deformation on the gear tooth contact state, then the effects of shaft misalignment on gear contact characteristics are studied. The structure of this paper is as follows: The problem needs to be addressed in this paper is analyzed in Sect. 2. The improved LTCA model for calculating the gear contact state parameters is introduced in Sect. 3. The proposed new mesh stiffness calculation model is described in Sect. 4. The 8-DOF dynamic model of misaligned gear system and the setting of the verification experiment are presented in Sect. 5. Then the results obtained from the above three models and experiment are discussed and analyzed in Sect. 6. Finally, the main conclusions are summarized in Sect. 7.

Problem formulation
The deviation of one gear pair with misalignment error is depicted in Fig. 1, where, z 1 and z 2 represent the ideal axis direction of gear #1 and gear #2 severally. The plane I is perpendicular to the plane II, in which the x direction is parallel to plane II and the y direction is perpendicular to plane II. The ideal axis of gear #1 will be offset when the shaft is misaligned, that is, compared with the ideal axis z 1 , the axis of gear #1 has a misalignment angle Dh x and Dh y in the x and y directions, respectively.
Shaft misalignment will cause the change of gear tooth contact state (including contact stress distribution and contact region), which may further affect the internal excitation and dynamic characteristics of gear system. Therefore, the accurate modeling of stiffness excitation of the misaligned gear pair considering the tooth contact status is necessary, which is beneficial to investigate the influence of the gear contact state changes caused by shaft misalignment on the dynamic characteristics of gear system, and also helpful to the accurate identification and prediction of misalignment errors.
As displayed in Fig. 2a, in the traditional mesh stiffness calculation model of misaligned gear [19], it is assumed that the contact force distribution is a specific parabolic distribution, and the length of contact line is consistent with the change of misalignment error, and is equal to the tooth width. However, as presented in Fig. 2b, when the misalignment error becomes larger, not only the distribution amplitude of gear contact stress will change, but also the contact line length will decrease sharply in the direction of Fig. 1 Deviations of the misaligned shaft [15] (a) (b) Fig. 2 Contact state of misaligned gear: a fixed contact line length [19], b variable contact line length tooth width [7,15]. According to Ref. [21], the tooth contact changes aroused by misalignment error have a large impact on the mesh stiffness and dynamic responses of gear system, which is commonly ignored in the traditional stiffness calculation model. Therefore, combined with the improved LTCA model, a stiffness calculation model considering the contact state change caused by shaft misalignment is proposed in this paper. Then, the effects of the contact state alterations on the dynamic characteristics of gear system are investigated base on the proposed model. The specific research content of the paper is given in the following chapters.

Model for calculating tooth contact state parameters
In this section, in order to accurately calculate the contact state parameters of misaligned gear and provide the parameter basis for the later calculation modeling of gear mesh stiffness, the traditional LTCA model is improved by considering the influence of fillet foundation deformation on the gear contact characteristic.

Traditional LTCA model
Traditional LTCA model [13,15] is established by using the influence coefficient method and considering the relationship between the contact area deformation and load. As shown in Fig. 3a, the relationship of the contact deformation d, the deformation w M1 and w M2 at points M 1 and M 2 and the separation distance h M can be expressed as [13,15]: if pointsM 1 and M 2 are beyond contact: ð2Þ As shown in Fig. 3b, the common tangent plane of the gear tooth contact region is divided into n grid units with an area of 2A 9 2B. Since the size of the grid unit is much smaller than that of the actual tooth contact region, the stress distribution in the grid unit can be simplified to a constant stress distribution, and then the contact process of the grid unit can be approximately regarded as satisfying the ideal Hertz contact condition. Therefore, the deformation w k of any grid unit k can be decomposed into the sum of the deformation w kj (contact deformation w H , kj , bending deformation w B , kj and shear deformation w S , kj ) generated by any grid unit j in the contact region when there is a force p j acting on the grid unit j, it can be deduced as: (a) (b) Fig. 3 a Contact between a tooth pair; b discretization of the contact region Through the conversion of above equations, the contact stress calculation expression matrix of the single-tooth pair can be expressed as [13], Under the assumption that the influence coefficients of two neighboring tooth pairs are independent from each other, the expression matrix for calculating contact stress of the two tooth pairs is proposed as follows [13,15], where, the tooth clearance matrix H can be calculated according to the spatial geometric relationship in Ref. [15]. The calculation expression of the deformation influence coefficient A of bending, Hertzian contact and shear can be found in Refs. [9,10,22,23], respectively. The contact stress of tooth pair can then be calculated iteratively by combining the non-negative condition of pressure and the correction method of the contact elastic deformation d [13,15].

Improved LTCA model
As shown in Fig. 4, when the gear is in meshing, the gear tooth #2 will not only have the Hertzian contact, bending and shear deformation under the action of meshing force F 2 , but also a certain fillet foundation deformation (represented by d f2 in Fig. 4) due to the flexibility of gear body. In addition, the coupling effect of the fillet foundation deformation between the neighboring gear teeth appears during the doubletooth meshing stage, that is, the fillet foundation deformation of gear tooth #1 under the meshing force F 1 will cause certain deflection of the gear tooth #2 (denoted by d f21 in Fig. 4). It can be seen that the fillet foundation deformation and the coupling deformation will directly influence the parameters in contact clearance matrix H between the meshing tooth pairs. However, the effects of the above two kinds of deformations on the contact characteristics of gear pair are ignored in the traditional LTCA model [13,15]. Therefore, the traditional LTCA model is improved in our study considering the influences of the fillet foundation deformation on the contact characteristics of gear pairs. The uneven distribution of meshing forces along the tooth width of misaligned gear will lead to the inconsistent distribution of the fillet foundation deformation along the tooth width direction. In this study, it is still assumed that the fillet foundation deformations of all tooth cross sections in the direction of tooth width keep the same based on the equivalence principle for the convenience of modeling. Thus, the improved LTCA model of the double-tooth pairs can then be written as, where the clearance variable matrix caused by the fillet foundation deformation can be calculated as:  In this paper, the fillet foundation deformation d fr (r = 1, 2) of the single tooth is calculated by the following equation [24][25][26], where, the parameters u f , S f , L, M, P and Q can be found and obtained in Ref. [24]. In the LTCA model of double-tooth pairs, the coupling deformation d frq (r = 1, 2; q = 1, 2, r = q) caused by the fillet foundation of adjacent gear tooth can be calculated as [27,28], where, F q (q = 1, 2) refers to the mesh force on gear tooth pair #1 (q = 1) or tooth pair #2 (q = 2) calculated by the load distribution during the gear meshing. The parameters in Eq. (11) can be seen in [27,28]. Shaft misalignment may cause the reduction in contact line length of gear tooth, which will also affect the calculation of the fillet foundation deformation. In this paper, at the beginning of the stress iterative calculation, the contact state of the gear is set as fulltooth contact, that is, the contact line length equals to the tooth width W. After the completion of a new round of iterative calculation, the new contact line length L 0 obtained from the calculation will replace the contact line length W in Eqs. (10) and (11), and then it will be replaced into the improved LTCA model for next iterative calculation. The iteration will end until the change value of the contact line length is within the allowable error range. Mesh stiffness of misaligned gear mainly includes the contact stiffness, the tooth stiffness (bending stiffness, shear stiffness and axial compression stiffness), the fillet foundation stiffness and the torsional stiffness. Therefore, considering the changes of contact state parameters caused by misalignment error, the above four stiffness calculation formulas are modified in this section, and the mesh stiffness calculation model of misaligned gear is finally established.

Calculation of contact stiffness
The non-ideal Hertz contact stiffness algorithm is a contact stiffness calculation method based on the elastic body force-deformation relationship, which is applicable for the calculation of contact stiffness of any elastic body that does not meet the ideal Hertz contact conditions [13]. Therefore, based on the improved LTCA model, the acting load F and the corresponding contact elastic deformation d under different misalignment errors and different meshing positions can be calculated, and then the contact stiffness k h can be obtained as, 4.2 Calculation of tooth stiffness Figure 5 illustrates the physical cantilever beam model of normal tooth. The bending energy U b , shear energy U s and axial compression energy U a stored in the tooth cantilever beam under the action of meshing force F can be calculated as [29][30][31][32], in which, As displayed in Fig. 6, the existence of misalignment error will lead to the direction change of meshing force, which is different from the force conditions of the aligned gear. When the gear exists the misaligned angle deviation Dh x and Dh y shown in Fig. 1, the component force of the meshing force of misaligned gear tooth in Eq. (14) can be deduced as, where, tan Df ¼ tan Dh x Á cos Dh y [15].
The effective contact line length L 0 of the tooth pair can be determined according to the length of the stress distribution region calculated by the improved LTCA model proposed in Sect. 3 when the misalignment error parameters are determined. Then, for the misaligned gear tooth, the parameters I x and A x in Eq. (14) can be redefined as, Combined with Eqs. (13), (16) and (17), the bending stiffness k b , shear stiffness k s and axial compression stiffness k a of the misaligned gear can be finally deduced as follows, Þsin a 1 þ cos a 1 À cos a 3 À cos a þ a 3 À a ð Þsin a þ cos a 3 ½ È À sin a 1 cos Dn À cos a 1 sin Dh y sin Dn Þcos a sin a 1 cos Dn À cos a 1 sin Dh y sin Dn

Calculation of fillet foundation stiffness
In this paper, the fillet foundation deformation formula of the traditional aligned tooth is improved based on the tooth contact line length L 0 obtained from the improved LTCA model. In the improved calculation formula, the fillet foundation deformation of the single-tooth and the coupling fillet foundation deformation aroused by the nearby teeth are all considered.
In the calculation process, the tooth width W in Eqs. (10) and (11) should be replaced by the actual contact line length L 0 in the last iteration. Then the above two kinds of fillet foundation stiffness can be obtained as,

Calculation of torsional stiffness
As shown in Fig. 6, the shaft misalignment will lead to uneven distribution of gear contact stress along the tooth width direction, which will further generate torsional torque T of gear tooth. The torsional energy stored in the tooth cantilever beam can be represented mathematically as [19], According to the contact line length L 0 obtained from the improved LTCA model, I px can be calculated as, The torsional torque T can be denoted as, Thus, the final calculation expression of the torsional stiffness k t can be written as, 5 Dynamics modeling of gear system and verification experiment setup

Dynamic modeling of gear system with shaft misalignment
In this section, an 8-DOF dynamic model of misaligned spur gear system is established to investigate the effects of tooth contact state variation caused by shaft misalignment on dynamic characteristics of gear system, which is presented in Fig. 7. Different from the 6-DOF model of the spur gear system proposed in Refs. [33,34], the translational DOF in the axial direction due to the angle tilt generated by the shaft misalignment is considered in the proposed 8-DOF model. Combined with the spatial coordinate transformation matrix M P in [15], the comprehensive deformation d mis of the misaligned gear can be calculated as follows, d mis ¼Dx cos a sin Dh y þ Dy sin a cos Dn À cos a sin Dh y sin Dn À Á þ Dz sin a sin Dn þ cos a sin Dh y cos Dn À Á þ h 1 cos Dh y r 1 À Á þ h 2 cos Dnr 2 À sin Dn sin Dh y r 1 À Á where, Dx¼x 1 À x 2 , Dy ¼ y 1 À y 2 , Dz ¼ z 1 À z 2 . Then the dynamic meshing force can be expressed as, According to Eq. (16), the components of the dynamic meshing force along x, y and z directions can be deduced as, Then, the four motion equations of the driving gear along x, y, z and rotation directions can be obtained as follows, The motion equations of the driven gear can be expressed as,

Verification experiment setup
In order to verify the accuracy of the simulation results obtained from the 8-DOF dynamic model, a verification experiment setup was established in this paper. The structure of the verification experiment system is presented in Fig. 8a, which is mainly composed of the two-stage variable speed gear system (including gear, shaft, bearing, box, etc.), the signal acquisition system (including triaxial acceleration sensor, computer, signal acquisition and processing software, etc.), and the power control system (including control cabinet, motor, magnetic powder brake, etc.). The system signal acquisition and processing system is mainly responsible for collecting and processing the output signals of the operation gear system. The power control system is mainly used to control the input speed of the motor and the torque of the load end and other parameters. Figure 8b shows the set position of the misalignment error in the gear transmission system, in which the bearing outer ring is nested in the bearing end cover, and the bearing end cover is fixed with the gear box by bolts. In the experiment, the shaft misalignment of gear system is simulated by seeding the misalignment error with 0.5°on the inner hole surface of bearing end cover. The triaxial acceleration sensor is mounted on the bearing end cover with misalignment error.
The basic parameters of the spur gear pairs in the two-stage transmission gearbox are displayed in Table 1. In this experiment, the input speed and load torque are set as 700 rpm and 25 NÁm, separately. And the rotation frequency of each shaft and the meshing frequency of gear pairs under input speed with 700 rpm are listed in Table 2.

Results and discussions
Based on the theoretical models proposed in the previous sections, the characteristics of the tooth contact, stiffness excitation and dynamic responses of the misaligned gear system are simulated and analyzed in this section. The gear parameters used in the simulation calculation are presented in Table 3.

Tooth contact state analysis under different misalignment error amplitudes
At the meshing position a 1 = 16°, the contact stress distribution of the driving gear tooth surface under different misalignment errors is shown in Fig. 9. It is seen from Fig. 9a that the tooth contact stress of the aligned gear is basically evenly distributed in the tooth width direction, and only a small stress concentration appears at the tooth edge. Compared with the aligned gear, the misaligned gear is still in contact with fulltooth width when the misalignment error is relatively small, but the tooth edge stress concentration of the misaligned gear is greater than that of the normal gear, which is seen in Fig. 9b. As illustrated in Fig. 9c and d, with the misalignment error going up, the degree of stress concentration increases gradually, and the effective tooth contact line length gradually becomes shorter and is no longer equal to the full-tooth width. Compared with the normal case, the maximum distributed stress of gear tooth increases by about 60%, 150% and 250% severally when the misalignment error Dh x = Dh y = 0.05°, 0.1°and 0.15°, respectively. It can be concluded that the contact stress distribution and the effective contact line length of the gear tooth will change obviously when the gear tooth have larger misalignment errors. Mass moment of inertia (kg m 2 ) 5.764 9 10 -4 9.12 9 10 -4 Bearing stiffness(N/m) 6.4 9 10 6 6.4 9 10 6 Bearing damping (N /m) 5.8 9 10 3 5.8 9 10 3 Input speed (rpm) 1000 -Input torque (N m) 100 -

Tooth contact state analysis under different misalignment directions
At the meshing position a 1 = 16°, the tooth contact stress distributions of the gear with misalignment error Dh x = 0.1°, Dh y = 0°and Dh x = 0°, Dh y = 0.1°are shown in Fig. 10a and b, respectively. It can be calculated that compared with the normal case, the maximum distributed stress of gear tooth increases by about 100%, 35% severally when the misalignment error only appears in x and y directions, respectively. Compared with the error in y direction, the misalignment error in x direction has a more significant influence on the tooth contact state, and can lead to greater stress concentration and smaller effective contact line length under the same misalignment error. This phenomenon can be explained as: compared to the error in y direction, the misalignment error in x direction can lead to a great change in the curvature of tooth meshing position, then results in a rapid change in the tooth contact state. Figure 11 presents the contact stress distribution of misaligned gear tooth under different meshing positions when the misalignment error is Dh x = Dh y = 0.1°. Figure 11a displays the tooth contact stress distribution at the meshing position a 1 = 3.6947°. At this meshing position, the gear just enters the doubletooth meshing stage, and the corresponding meshing force is 946 N. It can be found that the stress concentration degree is obvious and the effective contact line length is shorter due to the smaller meshing force when the gear teeth just enter into engagement. When the meshing angle ups to a 1 = 13.6039° (Fig. 11b) corresponding meshing force is 1892 N. We can find that compared with the case in Fig. 11a, the tooth stress concentration decreases and the contact line becomes longer at this position. The results reveal that the tooth involute curvature becomes smaller and the contact meshing force goes up with the meshing angle increasing, which leads to the lengthening of the contact line and the decrease in the stress concentration.

Tooth contact state analysis under different meshing positions
The corresponding meshing position angle in Fig. 11c and d is a 1 = 13.6947°and 18.0947°, respectively. At the two positions, the gear has entered the single-tooth meshing state after double-tooth meshing, and the meshing force reaches 2875 N. In the single-tooth meshing state, the tooth bears the largest meshing force, which leads to longer contact line and relatively small stress concentration. The meshing position in Fig. 11e and f is a 1 = 18.1004°a nd 28.0039°severally, and the corresponding meshing forces are 1892 N and 946 N, respectively. In this stage, the rear adjacent tooth enters meshing, and the stress concentration increases and the contact line length becomes shorter with the increase in the meshing angle. It can be obtained that the change of the stress concentration degree and the contact line length of the meshing tooth are mainly affected by the meshing force and the tooth surface curvature at each meshing position.

Stiffness analysis of misaligned gear considering tooth contact state
It is seen in Sect. 6.1 that the contact stress distribution parameters of gear with different misalignment errors at different meshing positions can be calculated based on the improved LTCA model. Therefore, according to the stiffness calculation model proposed in Sect. 4, the effects of the tooth contact state changes caused by shaft misalignment on gear stiffness are simulated and analyzed in this section. Figure 12 illustrates the influences of different misalignment errors on the gear contact stiffness. It can be found that the contact stiffness raises with the increase in misalignment error, which is due to the fact that the higher misalignment error can cause higher contact stress concentration, and then leads to the increase in tooth contact stiffness. It can also be found that the misalignment error direction has different influences on the tooth contact stiffness. The misalignment error in the x direction leads to the maximum contact stiffness near the tooth tip position, while the misalignment error in the y direction causes the maximum contact stiffness at the tooth root position. Under the same misalignment error, the contact stiffness in the x direction is up to about 10 times that in the y direction, which is owing to the fact that the misalignment error in the x direction has an obvious effect on the tooth contact state.  Fig. 13. It is seen from Fig. 13 that the comprehensive mesh stiffness calculated by the traditional model is about 30% larger than the proposed model, which is attributed to the fact that the influence of tooth contact state change on gear stiffness is ignored in the traditional model. When the gear pair engages in the alternating position between the single and doubletooth meshing areas, a step phenomenon exists in the mesh stiffness of single-tooth pair calculated by the proposed model, which is not found in the results obtained by the traditional model. The main explanation for this step phenomenon is given in Fig. 14 that when the gear is in the alternating position between the single and double-tooth meshing areas, the effective tooth contact line length changes sharply due to the sudden change of meshing force, which eventually results the step phenomenon of mesh stiffness. It is also seen from Fig. 14 that the maximum change value of contact line length caused by shaft misalignment yields 26.9% compared to the full-tooth contact, and the maximum change occurs in the alternating meshing stage of single and doubletooth meshing areas, which is the main reason for the mesh stiffness differences between the proposed and traditional models. Figure 15 presents the effects of different misalignment errors on the mesh stiffness. It can be observed that the increase in misalignment error will lead to the decrease in gear mesh stiffness. Compared with the case with error 0.05°, the comprehensive mesh stiffness decreases by about 8.3%, 16.7% and 38.3% severally when the error equals to 0.08°, 0.1°and 0.15°, respectively, which is owing to the fact that the effective tooth contact line length becomes shorter with the increasing of misalignment error, and further leads to the decrease in gear mesh stiffness. Figure 16 shows the effects of misalignment error directions on the mesh stiffness. It can be found that the mesh stiffness decreases with the raising of misalignment error in the x direction, while yields an opposite regular in the y direction. This phenomenon can be explained as follows: the increase in misalignment error in the y direction has little effect on the effective tooth contact line length, and mainly leads to the increase in gear contact stiffness, which will cause  the increase in comprehensive mesh stiffness. Different from the error in y direction, the misalignment error in x direction will not only increase the contact stiffness, but also significantly shorten the contact line length. Due to the fact that the reduction influence of the shortened contact line length on mesh stiffness is greater than the augment effect of the increasing of contact stiffness on mesh stiffness, the misalignment error in x direction eventually leads to the raising of meshing stiffness. We can also observe that the step phenomenon of the mesh stiffness of single-tooth pair caused by the misalignment error in x direction is more obvious than that in the y direction, which is also attributed to that the misalignment error in x direction has a more significant impact on the tooth contact state.

Dynamic response analysis of misaligned gear system considering tooth contact state
Time domain and frequency domain signal analysis methods are often adopted to analyze the fault vibration characteristics, then providing a basis for evaluating the evolution degree of faults and putting forward effective fault diagnosis methods [35][36][37]. In this section, the time and frequency domain dynamic responses of gear system with different misalignment errors are simulated and analyzed based on the  Figure 17 shows the comparison of vibration displacement responses of the gear with misalignment error Dh x = Dh y = 0.1°and normal gear in x, y and z directions. The displacement responses of the misaligned gear in three directions are all larger than that of the normal gear, which is owing to the lower mesh stiffness of the misaligned gear. It is worth noting that the normal gear has no vibration response in the z direction, while the misaligned gear produces vibration responses in the z direction, which is due to that the axial component of the meshing force generates vibration excitation on the axial direction (z direction) of the misaligned gear.

Comparison of dynamic responses between aligned and misaligned gears
In order to analyze the frequency domain response differences between the misaligned and aligned gear, the vibration displacement frequency spectra of the two type gears in the y direction are given in Fig. 18. According to the gear parameters shown in Table 3, the gear meshing frequency can be calculated as f m = 416.67 Hz under input speed 1000 rpm. As can be seen from Fig. 18, the response frequencies of both aligned and misaligned gears are around the meshing frequency f m and its harmonic frequencies (2 fm, 3 fm, …). The frequency response amplitude of the misaligned gear is larger than that of the aligned gear, which indicates that the misalignment error can raise the vibration energy level in the frequency domain.

Dynamic response analysis under different misalignment error amplitudes
The gear displacement responses under different misalignment error amplitudes in x, y and z directions are presented in Fig. 19. It can be observed that the vibration responses in x and y directions increases at first and then goes down with the raising of misalignment error, while the vibration response in z direction keeps increasing. The reason of this phenomenon may due to that: when the misalignment error is relatively smaller, it has no effect on the contact line length, but mainly affects the magnitude of contact stress, which will increase the contact stiffness and eventually lead to smaller vibration response. In contrast, the contact line length will be significantly shorter when the misalignment error becomes larger, which will result in the reduction in the gear mesh stiffness, thus causing a larger vibration response. With the further raising of misalignment error, the stiffness increment of the single and double teeth alternation decreases, which will reduce the vibration responses. As the misalignment error increases, the gear axial deflection will keep increasing, which results in the raising of excitation force in the z direction, and finally leads to the continued growth of vibration responses in the z direction. The displacement frequency spectra of gear in the y direction under different error amplitudes are displayed in Fig. 20. It can be observed that the response frequencies of gears with different increase with the raising of misalignment error. It is necessary to point out that with the frequency increasing, the amplitude of harmonic frequencies will reach maximum at case Dh x = Dh y = 0.1°even  Figure 21 shows the time domain displacement responses of gear under different misalignment error directions. From Fig. 21 we can find that the vibration responses caused by the misalignment error in the x direction are greater than that in the y direction, which is attributed to the fact that the influence of tooth contact line length change generated by the misalignment error in the x direction is more significant than that in the y direction and has a greater impact on the mesh stiffness.

Analysis of vibration statistical indicators under different misalignment errors
In order to investigate the effects of different misalignment error on the vibration response statistical indicators, the misalignment errors with 0°, 0.05°, 0.08°, 0.1°and 0.15°are randomly combined in x and y directions for simulation analysis. And then the yz peak ratio, RMS value and kurtosis value of the vibration displacement responses under 25 working conditions are obtained, which are displayed in Fig. 22. It is observed from Fig. 22a that with the misalignment error in the y direction increasing, the yz peak ratio decreases significantly. As the misalignment error in the x direction increases, the y-z peak ratio changes slowly, indicating that the misalignment error in the y direction has a greater impact on the vibration response in the z direction. In Fig. 22b and c, the RMS value goes up with the raising of misalignment error in both the x and y directions, while the kurtosis value changes slowly. Therefore, according to the performances of the three indicators, the RMS is recommended to evaluate and identify the shaft misalignment error. Figure 23 displays the comparison between the experimental and simulation results of frequency domain signals of the aligned and misaligned gears in y direction when the input speed is 700 rpm. In Fig. 23, both the meshing frequency (f m2 = 172 Hz) and its harmonic frequencies (2f m2 , 3f m2 ) of the simulation results can be found in the experiment results. It should be noted that there is no side frequency in the simulation results, while some side frequencies with an interval of intermediate shaft rotation frequency f s2 (f s2 = 3.2 Hz) and its multiple frequencies (2f s2 , 3f s2 , …) appear in the experiment results, which may be due to the modulation phenomenon caused by other excitations in the periodic meshing of gear pairs [38][39][40]. By comparing the experimental frequency spectra of the aligned and misaligned gears, it can be seen that the appearance of misalignment error will lead to the aggravation of vibration responses and the increase in the response amplitude corresponding to  Figure 24 shows the comparison of time domain response signals of the aligned and misaligned gears at 700 rpm. It can be found that the vibration responses increase obviously when the gear is misaligned. And the axial vibration responses of the misaligned gear are obviously larger than that of the aligned gear. The experimental results are also in good agreement with the simulation results, and thus validate the correctness of the simulation.

Comparative analysis of simulation and experiment results
As shown in Table 4, the RMS values of gear vibration acceleration signals obtained from the experiment are calculated and analyzed. It can be obtained that the RMS values of vibration responses of the misaligned gear are greater than that of the normal gear in x, y and z directions. The RMS value of z direction changes obviously and the maximum increase reaches 884% compared to the aligned case, while the increase rate of x and y directions is close to each other. The above results are consistent with the simulation results, which further verifies the accuracy of the proposed simulation models.

Conclusions
In order to reveal the influence mechanism of the tooth contact state change due to shaft misalignment on the stiffness excitation and dynamic characteristics of gear system, the traditional LTCA model is improved considering the effects of fillet foundation deformation on the tooth contact state, and the influences of shaft misalignment on the tooth contact state are studied. In addition, based on the improved LTCA model, a new mesh stiffness calculation model considering tooth contact state alterations is proposed, and then the effects of misalignment error on gear stiffness are investigated. Finally, according to the established dynamic model of misaligned gear system, the impacts of shaft misalignment errors on the dynamic characteristic of gear system are simulated and analyzed, and finally the correctness of the simulation results is verified by experiment. Through the above research, some conclusions can be drawn as follows: (1) The increasing of misalignment error will reduce the tooth contact line length and increase the stress concentration degree. The misalignment error in the x direction has a larger influence on the gear contact state comparing to the error in the y direction (The x and y direction definition can be seen in Fig. 1). The curvature radius and load of the meshing position are the main factors affecting the tooth contact state. (2) Mesh stiffness obtained by the proposed stiffness calculation model is significantly lower than that calculated by the traditional model, which indicates that the influence of tooth contact state alteration caused by shaft misalignment on mesh stiffness cannot be ignored. It is worth noting that the tooth contact state caused by shaft misalignment affects the gear internal excitation and dynamic characteristics as an intermediate variable, which is difficult to be measured and presented in the experiment. The accurate measurement of the contact state parameters in the experiment is worthy of further study, which is also our follow-up research works. of Chongqing, China (Grant No. CYB21014) and the Chongqing Science and Technology Innovation and Application Development Special Project, China (cstc2020jscx-msxmX0194).
Data availability The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.