Analysis of Nonlinear Dynamic Characteristics of Involute Spline Drive Considering Misalignment and Instantaneous Flash Temperature

Tengfei Xue, Mingjie Shen^*, Xiangzhen Xue

College of Mechanical & Electrical Engineering, Shaanxi University of Science & Technology, Xi’an, China.

DOI: 10.4236/wjet.2025.134054                

Abstract

For ease of installation, the involute spline used in aerospace, automobile and other fields will have clearance. Due to the presence of backlash, involute splines often experience misalignment during torque transmission. Under the influence of misalignment, the load distribution is highly concentrated on individual spline teeth, and frictional heat generation is also severe, increasing the risk of adhesive wear and tooth surface bonding of spline teeth. At the same time, these risks will lead to permanent damage to the spline tooth profile and further increase of the tooth gap, thus introducing harmful vibration, generating abnormal noise, and reducing transmission accuracy and efficiency. Based on this, this study constructed a dynamic model of involute spline transmission considering misalignment. On the basis of Blok flash temperature theory, the contact temperature of the spline working tooth surface and the time-varying backlash affected by temperature were obtained under misaligned working conditions. On this basis, the dynamic load characteristics of involute spline teeth under the combined influence of misalignment and contact temperature, as well as the effect of misalignment on instantaneous flash temperature, were explored. The research shows that the dynamic load in the process of spline transmission increases with the increase of contact temperature. Under the condition of a misalignment of 1 × 10⁻⁴ m and a temperature of 100˚C in the spline body, the maximum meshing force of the spline teeth is 183.4 N, and the maximum dynamic meshing torque is 46.81 Nm. The spline tooth working surface’s instantaneous flash temperature peak and its frequency of occurrence are positively correlated with the degree of misalignment.

Keywords

Involute Spline, Misalignment, Instantaneous Flash Temperature, Backlash, Dynamic Load

1. Introduction

Owing to their superior torque transmission and inherent self-alignment, spline couplings are extensively utilized in mechanical systems. In consideration of manufacturing error and installation convenience, the spline tooth side usually has side clearance. The presence of backlash will cause misalignment of the shaft between the splined shaft and the splined hub during transmission [1]. However, the presence of misalignment will also result in uneven load distribution and not all teeth will engage. During the meshing transmission process, the gear with small backlash shall be engaged first, and the gear shall be engaged successively according to the increase of backlash [2]-[4]. Moreover, the bending moment caused by misalignment further affects the contact behavior of the teeth [5].

Liu et al. [6]-[8] investigated the dynamics of an aero-engine’s two-rotor system by developing a model of the unbalanced rotor. Through analysis and simulation, they revealed the resulting vibration characteristics and their transmission law. Xiao [9] studied the action mechanism of misalignment on spline pair wear through self-made spline pair test stand. The test results show that the existence of misalignment increases the strength of the vibration signal of the splined coupling and also increases the wear of the gear teeth. For the solution of mesh stiffness, the traditional stiffness calculation method has errors when misalignment exists. For this purpose, Al-Hussain et al. [10]-[12] analyzed the effect of misalignment on torsional and lateral rotor system response and derived the corresponding stiffness matrix and force vector coupling. Their results confirm that misalignment has a pronounced influence on system dynamics. Huang et al. [13] developed an enhanced calculation model for the lateral and angular stiffness of spline couplings that accounts for parallel misalignment. The accuracy of this stiffness calculation method is satisfactory.

On the other hand, the dynamic engaging force, key to characterizing the system’s response, is therefore paramount in splined rotor system dynamics research. Niu et al. [14] analyzed the dynamic mesh stiffness and engagement force between teeth in harmonic gear transmission, and found that the engagement force generated by the teeth in the middle part was larger, while the engagement force generated by the teeth in the front and rear parts was smaller. Zhang et al. [15] developed an improved involute spline pair dynamics model. This model accounts for two types of misalignment: static, due to imperfections in manufacturing and installation, and dynamic, generated by the oscillatory relative motion of the involute spline coupling components. The results indicate that the existence of misalignment makes the spline teeth bear uneven load, especially under the condition of static parallel misalignment. Nevertheless, the dynamic misalignment has a minimal effect on the spline tooth load distribution. Xue et al. [16]-[18] proposed a dynamic model considering spline backlash and misalignment and studied the impact of misalignment on spline dynamics. The results show that the system enters into chaotic state with the increase of misalignment. Both the meshing forces between the teeth of involute splines and the dynamic load coefficient during spline transmission exhibit an increasing trend with greater misalignment.

The above studies have investigated the influence of misalignment on spline transmission in various aspects. However, the influence of variation in spline contact temperature on the dynamic characteristics is neglected. Variations in contact temperature will inevitably cause corresponding variations in spline backlash. Consequently, the backlash exhibits time-varying characteristics and is no longer a constant value, as traditionally assumed. In this case, the dynamic load in the splined transmission process is influenced. These are not currently considered for studies of misalignment. Therefore, this work aims to indirectly study the effect of contact temperature on dynamic loads (such as the dynamic engagement force) in the spline transmission process by examining the role of contact temperature on clearance. The effect of misalignment on spline contact temperature is also investigated. It can provide more accurate numerical analysis results of dynamic load for the wear prediction of involute spline pair.

2. Time Varying Backlash Calculation

2.1. Contact Temperature Calculation

Contact temperature ∆_B of spline working face is sum of the bulk temperature ∆_M and the instantaneous flash temperature ∆_f of the tooth surface [19].

${\Delta }_B={\Delta }_M+{\Delta }_f$ (1)

Assuming that the external and internal splines have identical tooth surface temperatures during the transmission process, the tooth profile clearance for any pair of spline teeth is defined by the following formula at any instant.

$c_j\left(t\right)=c_0-2\lambda \cdot s\cdot {\Delta }_{Bj}\left(t\right)\cdot \cos \alpha$ (2)

where, c₀ is the initial spline clearance, m. λ is the linear expansion coefficient, 1/˚C. s is the spline tooth thickness, m. α is the pressure angle.

Instantaneous flash temperature of working tooth surface of spline is caused by relative sliding of internal and external spline caused by vibration. The energy consumed due to friction during relative sliding will be converted into heat. From the Blok flash temperature theory, the instantaneous flash temperature ∆_f of the key teeth can be obtained [19].

${\Delta }_f=\frac{u{f}_e{f}_m\left|v_1-v_2\right|}{\sqrt{B}\left(\sqrt{g_1{\rho }_1{c}_1{v}_1}+\sqrt{g_2{\rho }_2{c}_1{v}_1}\right)}$ (3)

In the formula, u is the temperature rise coefficient, and for a pair of spline pairs, u = 0.83 is taken. f_e is the normal load on each tooth width (N/m). f_m is the friction coefficient. v_i (i = 1, 2) is the tangential velocity on the outer and inner spline tooth surfaces (m/s). g_i (i = 1, 2) is the thermal conductivity, respectively (J/(m∙s∙˚C)). ρ_i (i = 1, 2) is the density (kg/m³). c_i (i = 1, 2) is the specific heat capacity.

B, the semi-width (m) of the Hertzian contact band, is determined from Hertz contact theory under the assumption of identical tooth widths for both gears [20].

$B\left(t\right)=\eta \sqrt{\frac{2\left(1-{\mu }^2\right)}{E}\cdot \frac{F}{b}\cdot \frac{R_1\left(t\right)R_2\left(t\right)}{R_1\left(t\right)+R_2\left(t\right)}}$ (4)

Herein, η is calculation coefficient and η = 1.127 for a pair of spline pairs. E is Young’s modulus (GPa). μ is Poisson ratio. F is the normal load (N). b is the tooth surface width of the external and internal splines (m). R_i(t) (i = 1, 2) is the curvature radius of the external and internal splines at the meshing point (m) and is represented by the equations below.

$R_i\left(t\right)=r_{ci}\left(t\right)\sin {\alpha }_{mi}\left(t\right)$ (5)

Here, r_ci(t) (i = 1, 2) is the distance from the engagement-point to the rotational center of the external and internal splines, respectively (m). α_mi(t) (i = 1, 2) is the angle between the involute and the generating line [19].

${\alpha }_{mi}\left(t\right)=\arccos \left(r_{bi}/r_{ci}\left(t\right)\right)$ (6)

Among them, r_bi (i = 1, 2) is the base circle radius of the involute external and internal splines, and for involute splines, r_b1 = r_b2. The tangential velocity is formulated as follows:

$v_i\left(t\right)=\omega r_{ci}\left(t\right)\cos \alpha$ (7)

Here, α is the pressure angle(˚). ω is the angular velocity of the involute splines, respectively (rad/s). r_ci(t) is solved in section 2.2.

2.2. The Distance between the Engagement-Point and the Spline’s Rotational Center

The position of the engagement-point on the spline tooth profile is constantly changing due to vibration. Therefore, to obtain the length from the point of engagement to the spline’s rotational center, r_ci(t), the vibration displacement in the tooth profile direction shall be obtained first. In solving spline dynamic engagement parameters, most researchers approximate the involute tooth profile with a trapezoidal shape, resulting in a clear and computationally efficient solution. As shown in Figure 1(a), an external spline tooth is an example.

The angle of the jth tooth of the outer spline at a certain moment is defined as θ_j.

${\theta }_j=\frac{2\pi j}{Z}+\omega t+{\theta }_0$ (8)

where Z represents the number of spline teeth, ω represents the spline angular velocity, t is time, and θ₀ is the half angle (rad) of the tooth thickness at the reference circle, given by θ₀ = π/2Z. Then, the angle between the working tooth flank line of the involute spline and the X-direction is defined as φ_j,

Figure 1. Engagement line displacement solution diagram. (a) Angle diagram (b) Position of tooth profile lines before and after vibration.

${\phi }_j={\theta }_j-\alpha =\frac{2\pi j}{z}+\omega t+{\theta }_0-\alpha$ (9)

As shown in Figure 1(b), the intersection point between the spline working tooth profile and the pitch circle is located at point A at the initial moment. During the rotation of spline, due to the existence of lateral vibration displacement, the point of engagement moves from point A to point A₃ (if there is no vibration, the meshing point shall move to the intersection point A₂ between the reference circle and the spline tooth profile). It can be seen from the translation law that the straight lines of line JA and line KA₁ are parallel to each other, and the corresponding working tooth profile side lines are also parallel to each other. Draw the straight line LJ perpendicular to AA₂ and the straight line LK perpendicular to A₁A₃ through point L. From the geometric relation, it can be obtained that: ∠LA₃A₂ = φ_j, ∠KLA₃ = ∠LA₃A₂ = φ_j, ∠A₂AL = ∠LA₃A₂ = ∠KLA₃.

Assume that the outer spline vibrates for a distance of x₁ in the X-axis direction at a certain time, and the vibration distance in the Y-axis direction is y₁. That is, the length of LA₃ is x₁, and the length of LA is y₁. Therefore, the displacement in the meshing direction is the length of line segment AA₂ or line segment A₁A₃. As can be seen from the figure, the length of segment AA₂ is equal to the length of JA minus the length of JA₂. That is, AA₂ =JA – JA₂. In ▭JA₂A₃K, ∆LKA₃, ∆LJA: JA₂ = KA₃, KA₃ = LA₃ sinφ_j = x₁sinφ_j, JA = LA cosφ_j = y₁cosφ_j.

The vibration displacement of Involute spline teeth along the meshing line direction can be obtained from the geometric relationship as follows:

$n_1=l_{{\text{AA}}_2}=l_{\text{JA}}-l_{{\text{JA}}_2}=y_1\cos {\phi }_j-x_1\sin {\phi }_j$ (10)

From this, we can obtain the total relative displacement of the spline along the meshing line as follows [21]:

$\Delta n_j\left(t\right)=\left(x_1-x_2\right)\sin {\phi }_j-\left(y_1-y_2\right)\cos {\phi }_j+r_b\left({\theta }_2-{\theta }_1\right)$ (11)

It is known that AA₂ = y₁sinφ_j – x₂cosφ_j. In triangle ∆A₂AH, tanφ_j = A₂H/AA₂, A₂H = (y₁sinφ_j – x₂cosφ_j)tanφ_j. In triangle ∆LHA₃, HA₃ cosφ_j = LA₃ = x₁, so HA₃ = x₁/cosφ_j. Therefore, the vibration displacement of spline in the tooth profile direction is:

${\Delta }_{rj}={\text{A}}_{\text{2}}{\text{A}}_3={\text{A}}_{\text{2}}\text{H}+{\text{HA}}_3=\left(y_1\cos {\phi }_j-x_1\sin {\phi }_j\right)\tan {\phi }_j+x_1/\cos {\phi }_j$ (12)

Due to the influence of vibration during the rotation process, the meshing point of the spline may move towards the root circle or the tip circle of the spline. In both cases, the length from the point of engagement to the center of rotation of the spline is different.

The vibration displacement in the tooth profile direction determines the position of the engagement-point. So the regulations are as follows:

When the first derivative of the vibration displacement in the tooth profile direction is less than 0, the meshing point moves towards the tooth root circle. When the first derivative of the vibration displacement in the tooth profile direction is greater than 0, the meshing point moves towards the tooth tip circle.

As shown in Figure 2, it is a schematic diagram taking the movement of the engagement-point towards the tooth root circle as an example. Caused by vibration, the meshing point shifts from A₂ to A₃.

Figure 2. Position of the engagement point after vibration impact.

The distance between the engagement-point and the spline center can be obtained from geometric relationships and cosine theorem,

$r_{ci}\left(i=1,2\right)=\left|{\text{OA}}_3\right|=\left|\sqrt{\begin{array}{l}{\left(\left(y\cos {\phi }_j-x\sin {\phi }_j\right)\tan {\phi }_j+\frac{x}{\cos {\phi }_j}\right)}^2+{\left(\frac{d}{2}\right)}^2\\ -d\left(\left(y\cos {\phi }_j-x\sin {\phi }_j\right)\tan {\phi }_j+\frac{x}{\cos {\phi }_j}\right)\cos \alpha \end{array}}\right|$ (13)

Among them, d/2 is the radius of the spline indexing circle, which is also the length of OA₂. Similarly, when the point of engagement moves towards the tooth tip circle, the distance between the engagement-point and the spline center can be obtained from geometric relationships and the sine theorem as follows:

$r_{ci}\left(i=1,2\right)=\frac{\frac{d}{2}\sin \left(\pi -\alpha -\arcsin \left(\frac{2\left(\left(y\cos {\phi }_j-x\sin {\phi }_j\right)\tan {\phi }_j+\frac{x}{\cos {\phi }_j}\right)\sin \alpha }{d}\right)\right)}{\sin \alpha }$ (14)

The parameters of the involute spline are shown in Table 1. At this point, the instantaneous flash temperature and clearance change of the spline meshing tooth surface can be obtained.

Table 1. Involute spline parameter table.

Parameter	Value
Number of teeth/Z	16
Module/m (mm)	2
Pressure angle/α (˚)	30
Tooth surface width/b (m)	0.03
linear expansion coefficient/λ (1/˚C)	1.2 × 10−5
Thermal conductivity/g (J/(m∙s∙˚C)	46.47
Density/ρ (kg/m3)	7800
tooth thickness/s (m)	3.14 × 10−3
specific heat capacity/c (J/kg∙˚C)	481.48
Modulus of elasticity/E (N/m2)	2.1 × 1011
Poisson ratio/μ	0.3
Backlash/c0 (m)	2 × 10−4

3. Dynamic Load

As a critical factor characterizing the dynamic behavior of involute spline subsystems, the meshing force plays an essential role in research. The engagement force F_nj of each spline tooth along the meshing line primarily depends on the deformation of the contact surfaces and the mesh stiffness of the teeth [22].

$F_{nj}\left(t\right)=k_m\Delta +c_m\dot{\Delta }$ (15)

$\Delta$ is spline single-pair tooth engaging deformation function, and $\dot{\Delta }$ is derivative function of deformation function of single pair teeth engagement of spline. The expression is as follows:

$\Delta =\{\begin{array}{l}\Delta n_j\left(t\right)-c\left(t\right)\Delta n_j\left(t\right)>c\left(t\right)\\ 0c\left(t\right)\le \Delta n_j\left(t\right)\le c\left(t\right)\\ \Delta n_j\left(t\right)+c\left(t\right)\Delta n_j\left(t\right)<-c\left(t\right)\end{array}$ (16)

$\dot{\Delta }=\{\begin{array}{l}\dot{\Delta }{n}_j\left(t\right)\Delta n_j\left(t\right)>c\left(t\right),\Delta n_j\left(t\right)<-c\left(t\right)\\ 0-c\left(t\right)\le \Delta n_j\left(t\right)\le c\left(t\right)\end{array}$(17)

k_m is the mesh stiffness, and c_m is the mesh damping. The relatively small vibration amplitude, combined with the low thermal expansion coefficient of the spline material (AISI 9310), results in a minimal impact on the meshing stiffness. Therefore, this study adopts static meshing stiffness for solution calculation. Using the potential energy method to calculate the single tooth mesh stiffness of splines [23].

Finally, complete content and organizational editing before formatting. Please take note of the following items when proofreading spelling and grammar:

$k_m=1/\left(\frac{1}{k_b}+\frac{1}{k_s}+\frac{1}{k_h}+\frac{1}{k_a}\right)$ (18)

In the formula, k_b is the bending stiffness, k_s is the shear stiffness, k_h is the Hertz contact stiffness, and k_a is the axial compression stiffness.

The instantaneous meshing force for each spline tooth along the coordinate axis direction is defined as below.

$\{\begin{array}{l}{F}_{njx}\left(t\right)=F_{nj}\left(t\right)\sin \alpha \\ F_{njy}\left(t\right)=F_{nj}\left(t\right)\cos \alpha \end{array}$ (19)

The meshing moment is expressed as:

$T\left(t\right)=r_b{F}_{nj}\left(t\right)$ (20)

Hence, the component force of the total spline meshing force along the coordinate axis is:

$\{\begin{array}{l}{F}_{mx}\left(t\right)={\displaystyle {\sum }_{j=1}^z{F}_{njx}\left(t\right)}={\displaystyle {\sum }_{j=1}^z{F}_{nj}\left(t\right)\sin \alpha }\\ F_{my}\left(t\right)={\displaystyle {\sum }_{j=1}^z{F}_{njy}\left(t\right)}={\displaystyle {\sum }_{j=1}^z{F}_{nj}\left(t\right)\cos \alpha }\end{array}$ (21)

The total meshing moment resulting from spline engagement is:

$T_m\left(t\right)=r_b{\displaystyle {\sum }_{j=1}^z{F}_{nj}\left(t\right)}$ (22)

4. Model Establishment and Solution

Based on the lumped mass method, a three-dimensional dynamic model for involute spline transmission under misaligned conditions is developed, as illustrated in Figure 3. The numerical labels 1 through 4 in the model correspond to the motor, outer spline, inner spline, and load, respectively.

Figure 3. Dynamic model of involute spline.

The motor is coupled to the external spline, with a torsional stiffness of k_T1 and a torsional damping coefficient of c_T1 at the connection interface. Similarly, the internal spline is connected to the load, exhibiting a torsional stiffness of k_T2 and a torsional damping coefficient of c_T2 at its junction. The mesh stiffness k_m, mesh damping C_m and the time-dependent backlash caused by temperature changes are considered. The support stiffness and support damping along the coordinate axis at the outer spline are respectively k_px1, k_py1, c_px1 and c_py1. Similarly, k_px2, k_py2, and c_px2, c_py2 are the support stiffness and support damping of the inner spline in the X and Y directions, respectively. Motor input torque is T_d and load torque is T_L.

For misalignment settings, an external spline is an example. The coordinates of each axis of the external spline are changed from x₁, y₁ to ${x^{\prime }}_1$ , ${y^{\prime }}_1$ after being affected by misalignment [18].

Where, ${x^{\prime }}_1=x_1+l_x,{y^{\prime }}_1=y_1+l_y$ . l_x, l_y are the misalignment values of X axis and Y axis respectively. The derivation in Chapter 2.2 is carried out with l_x = l_y = 0. For the derivation of other misalignments, it is clear that the contents of Chapter 2.2 remain true.

Based on this, the differential equation of bending-torsion coupling dynamics of involute spline pair is established, such as Equation (23). See Table 2 for parameter values in the equation.

$\{\begin{array}{l}{J}_M{\ddot{\theta }}_M+k_{\text{T1}}\left({\theta }_M-{\theta }_1\right)+c_{\text{T1}}\left({\dot{\theta }}_M-{\dot{\theta }}_1\right)=T_d\\ m_1{{\ddot{x}}^{\prime }}_1+c_{\text{p}1}{{\dot{x}}^{\prime }}_1+k_{\text{p}1}{x^{\prime }}_1=F_{\text{mx}}\\ m_1{{\ddot{y}}^{\prime }}_1+c_{\text{p}1}{{\dot{y}}^{\prime }}_1+k_{\text{p}1}{y^{\prime }}_1=F_{\text{my}}-m_{1}g\\ J_1{\ddot{\theta }}_1+k_{\text{T}1}\left({\theta }_1-{\theta }_M\right)+c_{T1}\left({\dot{\theta }}_1-{\dot{\theta }}_M\right)=T_m\\ m_2{\ddot{x}}_2+c_{\text{p}2}{\dot{x}}_2+k_{\text{p}2}{x}_2=-F_{\text{mx}}\\ m_2{\ddot{y}}_2+c_{\text{p}2}{\dot{y}}_2+k_{\text{p}2}{y}_2=-F_{\text{my}}-m_{2}g\\ J_2{\ddot{\theta }}_2+k_{\text{T2}}\left({\theta }_2-{\theta }_L\right)+c_{\text{T2}}\left({\dot{\theta }}_2-{\dot{\theta }}_L\right)=-T_m\\ J_L{\ddot{\theta }}_L+k_{\text{T2}}\left({\theta }_L-{\theta }_2\right)+c_{\text{T2}}\left({\dot{\theta }}_L-{\dot{\theta }}_2\right)=-T_L\end{array}$(23)

Table 2. Basic parameters table.

Parameter	Value
Motor moment of inertia/JM	7.54 × 10−3 kg∙m2
Load moment of inertia/JL	7.52 × 10−3 kg∙m2
External spline moment of inertia/J1	3.43 × 10−3 kg∙m2
Internal spline moment of inertia/J2	3.73 × 10−3 kg∙m2
X-axis and Y-axis support stiffness/kp1, kp2	5 × 106 N/m
X-axis and Y-axis support damping/cp1, cp2	5 N∙s/m
Torsional stiffness of external spline shafts/kT1	7.855 N∙m/rad
Torsional stiffness of internal spline shafts/kT2	7.856 N∙m/rad
Torsional damping of external spline/cT1	0.191 N∙m∙s/rad
Torsional damping of internal spline/cT2	0.196 N∙m∙s/rad
Mesh damping between internal and external splines/cm	0.053 N∙s/m

Of all the system’s degrees of freedom, the unconstrained torsional degree of freedom permits rigid-body displacement of the spline about the Z-axis. This results in a non-positive-definite differential equation, preventing numerical convergence. Thus, the removal of rigid-body displacement is a critical requirement for the iterative numerical solving of the system’s equations of motion. To address this issue, both sides of the torsional equations in the system are divided by the moments of inertia J_M, J₁, J₂ and J_L, respectively, and a new set of degrees of freedom—Δ₁, Δ₂ and Δ₃ is introduced, yielding the following form:

$\{\begin{array}{l}{\Delta }_1=r_b\left({\theta }_1-{\theta }_M\right)\\ {\Delta }_2=r_b\left({\theta }_2-{\theta }_1\right)\\ {\Delta }_3=r_b\left({\theta }_2-{\theta }_L\right)\end{array}$

r_b is the base radius of the inner and outer splines (m). Presuming the prime mover’s angular velocity ${\dot{\theta }}_M$ remains constant and matches the system angular velocity, eliminate the rigid body displacement from the above equations. Combine the prime mover and outer spline torsional vibration equations, the outer spline and inner spline torsional vibration equations, and the inner spline and load torsional vibration equations in pairs. The above system of equations can be simplified to:

$\{\begin{array}{l}{{\ddot{x}}^{\prime }}_1+\frac{c_{\text{p}1}}{m_1}{{\dot{x}}^{\prime }}_1+\frac{k_{\text{p}1}}{m_1}{x^{\prime }}_1=\frac{F_{\text{mx}}}{m_1}\\ {{\ddot{y}}^{\prime }}_1+\frac{c_{\text{p}1}}{m_1}{{\dot{y}}^{\prime }}_1+\frac{k_{\text{p}1}}{m_1}{y^{\prime }}_1=\frac{F_{\text{my}}}{m_1}-g\\ {\ddot{x}}_2+\frac{c_{\text{p}1}}{m_2}{\dot{x}}_2+\frac{k_{\text{p}1}}{m_2}{x}_2=-\frac{F_{\text{mx}}}{m_2}\\ {\ddot{y}}_2+\frac{c_{\text{p}1}}{m_2}{\dot{y}}_2+\frac{k_{\text{p}1}}{m_2}{y}_2=-\frac{F_{\text{my}}}{m_2}-g\\ {\ddot{\Delta }}_1+k_{\text{T}1}\left(\frac{1}{J_1}+\frac{1}{J_M}\right){\Delta }_1+c_{\text{T}1}\left(\frac{1}{J_1}+\frac{1}{J_M}\right){\dot{\Delta }}_1=-\frac{r_b{T}_d}{J_M}+\frac{r_b}{J_1}{T}_m\\ {\ddot{\Delta }}_2+\frac{k_{\text{T2}}}{J_2}{\Delta }_3-\frac{k_{\text{T}1}}{J_1}{\Delta }_1+\frac{c_{\text{T}2}}{J_2}{\dot{\Delta }}_3-\frac{c_{\text{T1}}}{J_1}{\dot{\Delta }}_1=-r_b{T}_m\left(\frac{1}{J_1}+\frac{1}{J_2}\right)\\ {\ddot{\Delta }}_3+k_{\text{T2}}\left(\frac{1}{J_2}+\frac{1}{J_L}\right){\Delta }_3+c_{\text{T2}}\left(\frac{1}{J_2}+\frac{1}{J_L}\right){\dot{\Delta }}_3=-\frac{r_b{T}_m}{J_2}-\frac{r_b{T}_L}{J_L}\end{array}$ (24)

In SI units, when the magnitudes of the parameters in a differential equation differ greatly, solving the equation becomes extremely difficult, and it is hard to control the step size and error. Therefore, in actual engineering, it is necessary to perform dimensionless processing on such differential equation systems. Introducing dimensionless reference parameters ω and l, the equation can be made dimensionless as follows:

$\{\begin{array}{l}{{\ddot{\overline{x}}}^{\prime }}_1+\overline{c_{\text{p}1}}{{\dot{\overline{x}}}^{\prime }}_1+\overline{k_{\text{p}1}}{{\overline{x}}^{\prime }}_1=\overline{F_{\text{mx1}}}\\ {{\ddot{\overline{y}}}^{\prime }}_1+\overline{c_{\text{p}1}}{{\dot{\overline{y}}}^{\prime }}_1+\overline{k_{\text{p1}}}{{\overline{y}}^{\prime }}_1=\overline{F_{\text{my1}}}-\overline{g}\\ {\ddot{\overline{x}}}_2+\overline{c_{\text{p}2}}{{\dot{\overline{x}}}^{\prime }}_2+\overline{k_{\text{p2}}}{\overline{x}}_2=\overline{F_{\text{mx2}}}\\ {\ddot{\overline{y}}}_2+\overline{c_{\text{p}2}}{\dot{\overline{y}}}_2+\overline{k_{\text{p2}}}{\overline{y}}_1=\overline{F_{\text{my2}}}-\overline{g}\\ {\ddot{\overline{\Delta }}}_1+\overline{k_{\text{T}1}}{\overline{\Delta }}_1+\overline{c_{\text{T1}}}{\dot{\overline{\Delta }}}_1={\overline{T}}_1\\ {\ddot{\overline{\Delta }}}_2+\overline{k_{\text{T}23}}{\overline{\Delta }}_3-\overline{k_{\text{T21}}}{\overline{\Delta }}_1+\overline{c_{\text{T}23}}{\dot{\overline{\Delta }}}_3-\overline{c_{\text{T}21}}{\dot{\overline{\Delta }}}_1=-{\overline{T}}_2\\ {\ddot{\overline{\Delta }}}_3+\overline{k_{\text{T}3}}{\overline{\Delta }}_3+\overline{c_{\text{T3}}}{\dot{\overline{\Delta }}}_3=-{\overline{T}}_3\end{array}$(25)

In the formula:

$\overline{k_{\text{p1}}}=\frac{k_{\text{p1}}}{m_1{\omega }^2},\overline{c_{\text{p1}}}=\frac{c_{\text{p1}}}{m_1\omega },\overline{k_{\text{p2}}}=\frac{k_{\text{p2}}}{m_2{\omega }^2},\overline{c_{\text{p2}}}=\frac{c_{\text{p2}}}{m_2\omega },\overline{k_{\text{T1}}}=\frac{k_{\text{T1}}}{{\omega }^2}\left(\frac{1}{J_1}+\frac{1}{J_M}\right),$

$\overline{k_{\text{T3}}}=\frac{k_{\text{T2}}}{{\omega }^2}\left(\frac{1}{J_2}+\frac{1}{J_L}\right);\overline{c_{\text{T1}}}=\frac{c_{\text{T1}}}{\omega }\left(\frac{1}{J_1}+\frac{1}{J_M}\right),\overline{c_{\text{T3}}}=\frac{c_{\text{T2}}}{\omega }\left(\frac{1}{J_2}+\frac{1}{J_L}\right),\overline{k_{\text{T21}}}=\frac{k_{\text{T1}}}{J_1{\omega }^2},$

$\overline{k_{\text{T23}}}=\frac{k_{\text{T2}}}{J_2{\omega }^2};\overline{c_{\text{T21}}}=\frac{c_{\text{T1}}}{J_1\omega },\overline{c_{\text{T23}}}=\frac{c_{\text{T2}}}{J_2\omega };L=\frac{l}{r_b},{{\ddot{\overline{x}}}^{\prime }}_1=\frac{{{\ddot{x}}^{\prime }}_1}{l{\omega }^2},{{\dot{\overline{x}}}^{\prime }}_1=\frac{{{\dot{x}}^{\prime }}_1}{l\omega },{{\overline{x}}^{\prime }}_1=\frac{{x^{\prime }}_1}{l};$

${\ddot{\overline{x}}}_2=\frac{{\ddot{x}}_2}{l{\omega }^2},{\dot{\overline{x}}}_2=\frac{{\dot{x}}_2}{l\omega },{\overline{x}}_2=\frac{x_2}{l};{{\ddot{\overline{y}}}^{\prime }}_1=\frac{{{\ddot{y}}^{\prime }}_1}{l{\omega }^2},{{\dot{\overline{y}}}^{\prime }}_1=\frac{{{\dot{y}}^{\prime }}_1}{l\omega },{{\overline{y}}^{\prime }}_1=\frac{{y^{\prime }}_1}{l};{\ddot{\overline{y}}}_2=\frac{{\ddot{y}}_2}{l{\omega }^2},{\dot{\overline{y}}}_2=\frac{{\dot{y}}_2}{l\omega },$

${\overline{y}}_2=\frac{y_2}{l};{\ddot{\overline{\Delta }}}_1=\frac{{\ddot{\Delta }}_1}{l{\omega }^2},{\dot{\overline{\Delta }}}_1=\frac{{\dot{\Delta }}_1}{l\omega },{\overline{\Delta }}_1=\frac{{\Delta }_1}{l};$

$\overline{F_{\text{mx1}}}=\frac{F_{\text{mx}}}{m_{1}l{\omega }^2};\overline{F_{\text{my1}}}=\frac{F_{\text{my}}}{m_{1}l{\omega }^2};\overline{F_{\text{mx2}}}=\frac{F_{\text{mx}}}{m_{2}l{\omega }^2};\overline{F_{\text{my2}}}=\frac{F_{\text{my}}}{m_{2}l{\omega }^2};$

${\overline{T}}_1=\frac{r_b{T}_m}{J_{1}l{\omega }^2}-\frac{r_b{T}_d}{J_{M}l{\omega }^2};{\overline{T}}_2=r_b{T}_m\left(\frac{1}{J_{1}l{\omega }^2}+\frac{1}{J_{2}l{\omega }^2}\right);{\overline{T}}_3=\frac{r_b{T}_m}{J_{2}l{\omega }^2}-\frac{r_b{T}_L}{J_{L}l{\omega }^2}.$

5. Dynamic Response Analysis under Various Operating Conditions

5.1. Abbreviations and Acronyms

Figure 4(a) shows the relative displacement of involute spline with time in the direction of the line of engagement.

The results show a smooth sinusoidal time-domain response. The displacement-velocity phase diagram (Figure 4(b)) presents a smooth oval shape with the corresponding Poincare mapping (Figure 4(c)) as a single point. The distribution of meshing forces across individual spline teeth is presented in Figure 4(d). The engaging force between teeth is similar, about −170 N. All values are negative because the applied torque is clockwise (the counterclockwise torque is positive) [18]. The meshing component force distribution of each tooth of involute spline in X-axis direction is presented in Figure 4(e). The stress distribution exhibits a cosine pattern, and the resultant meshing component forces on each tooth are in a state of equilibrium (sum to zero). As shown in Figure 4, when misalignment and temperature-induced clearance variations are not considered, the dynamic response of the spline remains stable, and the forces are uniformly distributed among the spline teeth.

Figure 4. Dynamic response under ideal operating conditions. (a) Relative displacement along the engaging line. (b) Phase plane diagram. (c) Poincare map. (d) Meshing force of each tooth. (e) The meshing force of each tooth on the X-axis.

5.2. Analysis of Spline Response When Considering Misalignment

This section examines the dynamic response of the spline using X-direction parallel misalignment as a case study.

1) Dynamic response under parallel misalignment l_x = 5 × 10⁻⁵ m

When the parallel misalignment l_x = 5 × 10⁻⁵ m, as presented in Figure 5(a), the maximum relative displacement value of involute spline is 9.37 × 10⁻⁴ m. It increases slightly compared to the ideal condition and shows a certain fluctuation. The displacement-velocity phase diagram (Figure 5(b)) presents a plurality of annular patterns with 10 discrete points with periodic motion in the corresponding Poincare map (Figure 5(c)).

Additionally, Figure 5(d) shows that the engagement force distribution across the spline teeth is no longer uniform. The maximum and minimum engagement forces are about 175 N and 163 N. The engagement component force F_nx distribution of the involute spline tooth is still cosine (Figure 5(e)), but the engagement component force F_nx of individual teeth is almost zero. The dynamic load coefficient k_v fluctuation range of spline during transmission is 0.689 - 1.223, as shown in Figure 5(f).

2) Dynamic response under parallel misalignment l_x = 1 × 10⁻⁴ m

As the parallel misalignment l_x increases further, it reaches l_x = 1 × 10⁻⁴ m. The maximum relative displacement value of the spline is further increased to 9.66 × 10⁻⁴ m, as shown in Figure 6(a). The change of relative displacement with time is more disordered and no longer a smooth sine time domain response. The

Figure 5. Dynamic response under parallel misalignment l_x = 5 × 10⁻⁵ m. (a) Relative displacement along the engaging line. (b) Phase plane diagram. (c) Poincare map. (d) Meshing force of spline teeth. (e) The meshing force acting along the X-axis for each tooth. (f) Dynamic load factor.

displacement-velocity phase diagram (Figure 6(b)) shows more annular patterns and inflection points. Five more discrete points with periodic motion are added to the corresponding Poincare map (Figure 6(c)) compared to the case with misalignment l_x = 5 × 10⁻⁵ m. Furthermore, Figure 6(d) reveals that the load imbalance among the spline teeth is more pronounced. The maximum and minimum engagement forces are about 182 N and 157 N respectively, while the meshing component force F_nx of spline teeth is still distributed in cosine form (Figure 6(e)). The dynamic load coefficient k_v fluctuation range of spline during transmission is 0.69 - 1.225, slightly increased, as shown in Figure 6(f).

This section analyzes the spline pair’s dynamic response by treating parallel misalignment as a variable. Compared with the ideal working condition, the misalignment not only increases the vibration displacement of spline pair, but also causes the unbalance load of spline pair teeth. Both the spline pair’s maximum engagement force and the meshing transmission’s dynamic load coefficient grow with increasing misalignment. Therefore, the existence of misalignment destroys the stability and symmetry of the load transfer during the rotation process of the spline pair, and introduces additional excitation, which greatly affects the stability of the spline transmission.

5.3. Dynamic Load Response of Involute Spline Pairs at Different Bulk Temperature

In this part, the dynamic response of spline pair is studied by changing the temperature of spline duplicate under the condition of parallel misalignment l_x = 1 × 10⁻⁴ m as the reference condition. However, only the influence of temperature on the tooth profile clearance is considered in this study, and the influence of temperature change on the deformation and rigidity of profile is ignored. Therefore, only the dynamic load change caused by temperature change is analyzed in this part.

Figure 6. Dynamic response under parallel misalignment l_x = 1 × 10⁻⁴. (a) Relative displacement along the engaging line. (b) Phase plane diagram. (c) Poincare map. (d) Meshing force of spline teeth. (e) The meshing force acting along the X-axis for each tooth. (f) Dynamic load factor.

Dynamic response analysis of the involute spline at bulk temperature ∆_M = 50˚C with Misalignment quantity l_x = 1 × 10⁻⁴ m.

When the effect of temperature on tooth profile clearance is incorporated, the meshing force on each spline tooth exhibits minor variations, and the engagement force of each tooth increases by about 0.71 N, as shown in Figure 7(a).

Figure 7. Dynamic load analysis of the spline at bulk temperature ∆_M = 50˚C. (a) Meshing force. (b) X-axis meshing force of spline teeth. (c) Meshing torque. (d) Dynamic load factor.

The distribution of engaging force is consistent with that of misalignment l_x = 1 × 10⁻⁴ m. The influence on the meshing component force F_nx is almost negligible (Figure 7(b)). However, according to Figure 7(c), the dynamic meshing torque fluctuation range of spline pair is greatly improved. When only considering the influence of misalignment, the fluctuation range of meshing torque is 21.15 Nm -37.54 Nm, and the direction of meshing torque is clockwise. With the influence of temperature taken into account, the mesh moment fluctuates from 21.51 Nm to 46.64 Nm and no longer exhibits a sinusoidal time-domain distribution. Correspondingly, the dynamic load coefficient of the meshing transmission increases accordingly, and the fluctuation range is 0.702 - 1.521, as presented in Figure 7(d). The underlying cause is that thermal-induced expansion reduces spline backlash. This is due to the increased tooth thickness of the inner and outer splines caused by thermal expansion. This also reduces the slot width of the inner spline, and increases the tooth thickness of the outer spline, so the backlash of the spline decreases. Given that relative displacement in the meshing direction is constrained, the meshing displacement must consequently rise. An increase in meshing displacement necessarily leads to a corresponding rise in both meshing force and meshing torque. Therefore, the temperature change in the process of spline transmission is also one of the sources of external excitation of the system.

Dynamic response analysis of the involute spline at bulk temperature ∆_M = 100˚C with Misalignment quantity l_x = 1 × 10⁻⁴ m.

As shown in Figure 8(a), when the spline body temperature ∆_M = 100˚C, the variation in meshing force among the spline teeth remains minimal.

Figure 8. Dynamic load analysis of the involute spline at bulk temperature ∆_M = 100˚C. (a) Meshing force. (b) X-axis meshing force of spline teeth. (c) Meshing torque. (d) Dynamic load factor.

Each tooth engagement force is increased by about 1.4 N compared to the case considering only misalignment. The numerical variation of meshing component force F_nx is similar to that of meshing force (Figure 8(b)). As presented in Figure 8(c), the dynamic meshing torque of the spline pair fluctuates within a range of 21.68 Nm to 46.81 Nm, while the direction of torque transmission remains unchanged, maintaining a clockwise orientation. Correspondingly, the variation range of the dynamic load coefficient in the meshing transmission exhibits a marginal increase of approximately 0.006, as depicted in Figure 8(d). When the spline body temperature rises, the dynamic load between the spline pair tooth surfaces changes little because the spline temperature rise and the thermal expansion coefficient of the material itself are small. Therefore, the variation of spline backlash is also very small, and the impact on dynamic load is correspondingly reduced.

In this part, the effect of contact temperature on the dynamic load of involute spline pair is studied. The results show that the change of temperature has little influence on the dynamic engagement force of single tooth of spline pair. The temperature variation significantly affects the meshing moment and dynamic load coefficient, but this result is based on the existence of misalignment. This also indicates that misalignment has a greater impact on the instantaneous flash temperature of the spline drive.

5.4. Instantaneous Flash Temperature Change of Spline Tooth Surface under Different Misalignment

This section investigates the effect of misalignment on the instantaneous flash temperature under an input torque of 30.7 Nm and a rotational speed of 6000 r/min.

Instantaneous flash temperature change of the spline when misalignment quantity l_x = 0 m.

Irrespective of misalignment, the instantaneous flash temperature variation of the spline is presented in Figure 9, with the maximum instantaneous flash temperature of about 0.011˚C.

Figure 9. Without considering the instantaneous flash temperature during misalignment.

Irrespective of misalignment, the instantaneous flash temperature variation of the spline is presented in Figure 9, with the maximum instantaneous flash temperature of about 0.011˚C. The peak instantaneous flash temperature is attained at the maximum relative displacement between the internal and external splines along the tooth profile. At this instant, the relative tangential velocity between the external and internal splines reaches its maximum at the meshing point. The moment when the instantaneous flash temperature is 0 is the moment when the relative displacement of the inner spline and the outer spline along the tooth profile direction reaches the minimum. At this time, the engaging point is on the reference circle of spline, and the tangential velocity of inner spline and outer spline is equal. The temperature curve is very smooth.

Instantaneous flash temperature changes of splines under different misalignment quantities.

As shown in Figure 10, the presence of misalignment significantly increases the instantaneous flash temperature between the spline tooth flanks, which also indicates that misalignment significantly increases the frictional heat effect between the tooth flanks.

Figure 10. Instantaneous flash temperature changes under different misalignment quantities. (a) l_x = 2.5 × 10⁻⁵ m. (b) l_x = 5 × 10⁻⁵ m. (c) l_x = 7.5 × 10⁻⁵ m. (d) l_x = 1 × 10⁻⁴ m.

The maximum peak values of the instantaneous flash temperatures shown in Figures 10(a)-(d) also show a significant upward trend with increasing parallel misalignment. The misalignment varies from 2.5 × 10⁻⁵ m to 1 × 10⁻⁴ m, and the maximum instantaneous flash temperature is 4.6˚C, 9.6˚C, 15.91˚C and 20.96˚C. This trend clearly demonstrates that greater misalignment leads to elevated peak instantaneous flash temperatures on the tooth surface and intensifies frictional heating. In addition, the flash temperature curve becomes more complex as misalignment causes the splined mesh drive to become unstable. The temperature profile is no longer smooth, with a sudden rise in temperature. As indicated in Figure 10(c) and Figure 10(d), under conditions of significant misalignment, the occurrence frequency of elevated temperature peaks rises noticeably. The rapid increase and persistent presence of this flash temperature significantly elevate the risk of thermal damage failures, such as surface gluing, on the spline teeth, particularly under severe misalignment conditions. Therefore, in the design, manufacture and installation of spline transmission system, strict control of misalignment is essential to ensure its thermal reliability and service life.

The results of the analysis in this part show that misalignment increases the frictional heat effect in the process of splined transmission. As misalignment grows, the peak instantaneous flash temperature of the spline exhibits a continuous rise and occurs more frequently, thereby compromising the operational safety of the spline drive.

6. Conclusions

In this study, the involute spline transmission dynamics model considering misalignment is established based on the existing research. Through the Blok flash temperature theory, the instantaneous flash temperature change of spline working surface is obtained under different parallel misalignment conditions. (In the future, we will focus on the impact of angle misalignment on it.) Based on this foundation, the effect of contact temperature on the dynamic load during spline pair transmission is investigated.

Taking the parallel misalignment as an example, the dynamic response of spline is studied when the parallel misalignment l_x = 5 × 10⁻⁵ m and l_x = 1 × 10⁻⁴ m. Under the working condition of parallel misalignment l_x = 10⁻⁴ m, the relative displacement of spline caused by vibration reaches 9.66 × 10⁻⁴ m. With the increase of misalignment, the relative displacement and maximum engagement force of spline increase, and the stress of spline teeth is more uneven. In addition, the motion law of spline is more complicated and disordered. Obviously, the existence of misalignment will introduce external excitation to the spline transmission system, and have a greater impact on the stability of the spline transmission.

Contact temperature exerts only a minor influence on both the meshing force and the force components of spline teeth. However, as the contact temperature increases, the engagement force demonstrates a consistent upward tendency. Under the working condition of misalignment l_x = 1 × 10⁻⁴ m and spline body temperature ∆_M = 100˚C, the maximum engagement force of spline teeth is about 183.4 N. Under the working condition of comprehensive consideration of misalignment and contact temperature, the dynamic engagement moment and dynamic load coefficient (fluctuation range reaches 0.708 - 1.527) in the process of spline transmission have obvious improvement, compared with only considering misalignment working condition. Meanwhile, the change of contact temperature makes the change rule of meshing torque more irregular.

With the increase of misalignment, the instantaneous flash temperature peak value between spline tooth surfaces increases continuously and the temperature rise curve becomes more complex. The frequency of instantaneous flash temperature peak is also increased. The existence of misalignment greatly increases the frictional heat effect between tooth surfaces, and poses a threat to the safety of spline transmission.

Therefore, to ensure reliability in practical applications, the strict control of misalignment is imperative. This can be achieved through enforcing tighter manufacturing tolerances on spline geometry, adopting precision alignment tools during assembly to guarantee coaxiality, and implementing real-time condition monitoring for critical systems.

Acknowledgements

This work was supported by the Youth Innovation Team of Shaanxi Universities (2024), Shaanxi Province Qin Chuangyuan “Scientist + Engineer” Team construction of No. 2024QCY-KXJ-112.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Curà, F. and Mura, A. (2018) Theoretical and Numerical Evaluation of Tilting Moment in Crowned Teeth Splined Couplings. Meccanica, 53, 413-424.[CrossRef]
[2]	Barrot, A., Paredes, M. and Sartor, M. (2006) Determining Both Radial Pressure Distribution and Torsional Stiffness of Involute Spline Couplings. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 220, 1727-1738.[CrossRef]
[3]	Cuffaro, V., Curà, F. and Mura, A. (2013) Test Rig for Spline Couplings Working in Misaligned Conditions. Journal of Tribology, 136, Article 011104.[CrossRef]
[4]	Curà, F., Mura, A. and Gravina, M. (2012) Load Distribution in Spline Coupling Teeth with Parallel Offset Misalignment. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 227, 2195-2205.[CrossRef]
[5]	Curà, F. and Mura, A. (2013) Experimental Procedure for the Evaluation of Tooth Stiffness in Spline Coupling Including Angular Misalignment. Mechanical Systems and Signal Processing, 40, 545-555.[CrossRef]
[6]	Liu, B.Y. (2021) Nonlinear Dynamics Analysis of Aero-Engine Misaligned Dual Rotor System. M.A. Thesis, Shenyang University of Technology.
[7]	Ma, P.P. (2021) Vibration Characteristics of Aero-Engine Twin Rotor system. M.A. Thesis, Dalian University of Technology.
[8]	Chen, Y.M. (2018) Study on the Vibration Response Characteristics of Aero-Engine Twin Rotor System without Alignment. M.A. Thesis, Hunan University of Science and Technology.
[9]	Xiao, L., Xu, Y., Sun, X., Xu, H. and Zhang, L. (2021) Experimental Investigation on the Effect of Misalignment on the Wear Failure for Spline Couplings. Engineering Failure Analysis, 131, Article 105755.[CrossRef]
[10]	AL-Hussain, K.M. and Redmond, I. (2002) Dynamic Response of Two Rotors Connected by Rigid Mechanical Coupling with Parallel Misalignment. Journal of Sound and Vibration, 249, 483-498.[CrossRef]
[11]	Patel, T.H. and Darpe, A.K. (2009) Vibration Response of Misaligned Rotors. Journal of Sound and Vibration, 325, 609-628.[CrossRef]
[12]	Patel, T.H. and Darpe, A.K. (2009) Experimental Investigations on Vibration Response of Misaligned Rotors. Mechanical Systems and Signal Processing, 23, 2236-2252.[CrossRef]
[13]	Huang, W., Tian, H., Ma, H., Wang, P., Yang, Y. and Han, Q. (2023) An Improved Method for Calculating the Lateral and Angular Stiffness of Spline Couplings Considering Parallel Misalignment. Mechanism and Machine Theory, 189, Article 105436.[CrossRef]
[14]	Niu, H.B. (2020) Study on the Mesh Stiffness Analysis and Mesh Force Distribution Law of Harmonic Gear Transmission. M.A. Thesis, Tianjin University of Technology.
[15]	Zhang, C., Zhu, R., Chen, W., Wang, D., Yin, X. and Song, D. (2023) An Improved Dynamic Model of the Spline Coupling with Misalignment and Its Load Distribution Analysis. International Journal of Mechanics and Materials in Design, 20, 393-408.[CrossRef]
[16]	Xue, X., Huo, Q., Liu, J. and Jia, J. (2021) Nonlinear Dynamic Load Analysis of Aviation Spline Coupling with Mass Eccentricity and Misalignment. Advances in Mechanical Engineering, 13, 1-19.[CrossRef]
[17]	Xue, X., Huo, Q. and Liu, J. (2021) Nonlinear Vibration Characteristics of the Involute Spline Coupling in Aeroengine with the Parallel Misalignment. International Journal of Aerospace Engineering, 2021, 1-19.[CrossRef]
[18]	Xue, X., Liu, J., Li, Y., Sui, L., Lin, K. and Yu, W. (2023) Nonlinear Dynamic Analysis of Floating Involute Splines Considering Comprehensive Misalignment. Journal of Sound and Vibration, 555, Article 117720.[CrossRef]
[19]	Gou, X., Zhu, L. and Qi, C. (2017) Nonlinear Dynamic Model of a Gear-Rotor-Bearing System Considering the Flash Temperature. Journal of Sound and Vibration, 410, 187-208.[CrossRef]
[20]	AGMA 925-A03 (2003) Effect of Lubrication on Gear Surface Distress. https://www.renrendoc.com/paper/86834250.html
[21]	Yue, K., Wang, L., Lan, H., Yu, W., Shao, Y. and Chen, Z. (2024) Nonlinear Characteristics of Floating Involute Spline Coupling in a Hydraulic Disc Brake with Inner Hub Rim Crack Propagation. International Journal of Non-Linear Mechanics, 160, Article 104636.[CrossRef]
[22]	Xue, X., Yu, W., Lin, K., Zhang, N., He, X. and Jiang, Y. (2024) Nonlinear Dynamic Meshing Force Analysis of Peek-Based Polymer Composite Involute Spline Coupling System. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-Body Dynamics, 238, 194-208.[CrossRef]
[23]	Fu, C.G. (2000) Aircraft Engine Design Manual Volume 19: Rotor Dynamics and Whole Machine Vibration. 5th Edition, Aviation Industry Press.