A Stiffness Formulation for Spline Joints

Due to lack of knowledge in terms of their flexibility and deformation, spline joints are often assumed to be rigid in dynamic models of gearboxes, transmissions, and drivetrains. It was reported by several studies that compliance of a spline joint can substantially influence dynamic behavior, leading to rotor dynamic instability and nonsynchronous whirl [1–4]. Properly capturing the compliance of spline joints might be critical to the fidelity of dynamic analyses of drivetrains.

Only a few experimental studies have been reported in the literature to determine stiffness of spline joints. Among them, Ku et al. [5] determined the rotational stiffness (termed as tilting stiffness in this paper) of a spline joint by measuring the bending moments and angular deflections transmitted across an axial spline with a nonrotating shaft excited by an external shaker. They showed that the tilting stiffness decreases as the force and frequency of the excitation increase. Cura and Mura [6] measured torsional stiffness of a spline joint in both aligned and misaligned conditions to show that misalignment of a spline coupling would reduce torsional stiffness. These experimental studies provided experimental data of spline stiffness for the simplest torsional loading cases.

On the theoretical side, Marmol et al. [1] proposed an analytical model to calculate the torsional, radial, and tilting (rotational) stiffness values of spline joints based on the assumption that all spline teeth carry the same amount of load and the load is evenly distributed along the pitch line across face width. This assumption is often invalid since actual load distributions of splines usually vary from tooth to tooth and load on a specific tooth is usually unevenly distributed over a wide area as demonstrated clearly in Refs. [7–11]. Recently, Barrot et al. [12,13] proposed an analytical formula to calculate torsional stiffness of a spline joint. The model they employed for this purpose was limited to pure torsional loading since they also assumed that load distributions of all spline teeth are identical. Effects of tooth modifications, misalignments, combined loading conditions, and tooth indexing errors were also not accounted for in these theoretical studies.

This paper aims at proposing a general stiffness formulation for spline joints based on the semi-analytical spline load distribution model developed in Ref. [10]. The goal here is to define a fully populated stiffness matrix of a spline joint including radial, tilting, and torsional stiffness values as well as coupling between these motions defined by the off-diagonal terms of the stiffness matrix. With this stiffness formulation in place, a parametric study will be performed to investigate the influences of tooth surface modifications, misalignments, and tooth indexing errors on the spline stiffness matrix.

The stiffness of a spline joint relates reaction loads (forces and moments) acting on it to relative rigid body displacements between the external and internal components. In this study, relative displacement along the axial z direction (⁠$uz$⁠) between the external and internal spline is constrained and the axial load component $Fz$ is excluded. As such, there are five reaction load components (two radial forces and three moments) depicted by the reaction load vector $P=[FxFyMxMyMz]T$⁠, and five relative displacement components (two radial displacements and three rotational displacements) depicted by the relative rigid body displacement vector $Φ=[uxuyθxθyθz]T$ as shown in Fig. 1. Then, the stiffness matrix of a spline joint is defined as $K=∂P/∂Φ$⁠, which is written explicitly as

In the above equation, $C$ is the compliance matrix of a spline joint that includes tooth bending, tooth shear, tooth base flexibility, and torsional compliance of both the external and internal spline, as well as the contact compliance [10]. The n contact points are identified by discretizing each of the Z spline tooth contacts into P and Q contact cells along the profile and face width directions. $G$ is a matrix relating to tooth surface geometry at the contact points and $GT$ is the transpose of $G$⁠, $A$ is a vector of artificial variables needed by the simplex-type algorithm, I is the identity matrix, $n$ is the number of potential contact point pairs over the spline interface, and $Y$ is a vector of final gaps between all contact point pairs of which $Yi$ (⁠$i∈[1,n]$⁠) is the ith element. $F$ is a vector of contact forces between all contact point pairs of which $Fi$ (⁠$i∈[1,n]$⁠) is the ith element. Finally, $ε$ is a vector of the initial separations between all potential contact point pairs. Readers are referred to Ref. [10] for a detailed derivation of Eq. (2) including definitions of these submatrices and subvectors. With this semi-analytical model, two methods, namely, an analytical method and a numerical method, can be devised to compute the stiffness matrix $K$ of a spline joint.

In this equation, the matrices $Ĉ$ and $Ĝ$ depend on potential contact point pairs that are actually in contact, indicating that the stiffness of a spline joint may vary as spline load distribution changes. Factors that change the spline load distribution, such as torque level, tooth modifications, and indexing errors, will possibly influence the stiffness. Investigations on effects of these factors will be performed using a numerical procedure described subsequently.

In the above analytical stiffness formula, computing the inverse of the compliance matrix subset $Ĉ$ might require considerable computational time depending of the size of the compliance matrix, especially if the number of spline teeth is large and all or most spline teeth carry load. The flowchart of computation shown in Fig. 2 provides two avenues: one to handle splines having any arbitrary load distributions and another one to handle a special case where all spline teeth have identical contact conditions, i.e., all spline teeth have the same respective contact point pairs in contact as it is the case for pure torsional loading with no indexing errors. In this special case, the compliance submatrix $Ĉ$ has a unique structure

Accordingly, $Ĉ−1$ can be constructed from Q and W, requiring evaluation submatrices $[ZĈT+Ĉl]−1$ and $Ĉl−1$⁠. This can be expected to reduce computational time drastically compared to direct inversion of $Ĉ$⁠. Table 1 shows the CPU times taken for stiffness calculation using this blockwise inversion method and the direct inversion method for a few example cases. In this table, Z is the number of spline teeth, and P and Q denote the number of contact cells along the profile and face width directions, respectively. It is observed that as Z increases, the size of $Ĉ$ becomes larger and stiffness computation using the direct inversion approach takes much longer. However, if the blockwise inversion is used, the CPU time required becomes negligibly small. For instance, as the number of teeth increases from 25 to 60, CPU time using direct inversion increases by more than ten times from 3.3 s to 42.5 s, while CPU time using blockwise inversion only slightly increases from 0.03 s to 0.04 s.

The above blockwise inversion method applies only to cases where all spline teeth have identical contact conditions, such as pure torsion loading (with or without any nominal tooth modifications). In cases where misalignments or tooth indexing errors are present as well as cases when the spline joint is under a gear loading condition, load distributions vary from tooth to tooth such that blockwise inversion of $Ĉ$ in the form of Eq. (8) is no longer valid, and the entire $Ĉ$ matrix must be inverted directly. Fortunately, in these cases, many potential contact point pairs are not in contact so that the size of $Ĉ$ is reduced significantly, making computational demand for this still acceptable. Table 2 shows CPU time required to compute the K matrix of a helical gear loaded spline joint having $Z=60$ teeth. As helix angle increases, more contact point pairs are separated and the size of $Ĉ$ becomes much smaller. Consequently, the CPU time is reduced for larger gear helix angles. For instance, the CPU time required drops from 4.12 s to 0.81 s as the gear helix angle is increased from 15 deg to 25 deg. Sensitivity of computational time to misalignments of the spline is similar as shown in Table 3. As misalignment increases from 0.01 deg to 0.03 deg, required CPU time is reduced from 1.55 s to 0.73 s, solely because the size of $Ĉ$ is decreased. It should be noted that the computational times noted are for analyses of one particular spline load and one rotational position. As will be shown in Sec. 3, variance of load, along with analysis of many rotational positions, will be necessary to fully capture the variation of spline stiffness. This of course will further magnify the importance of the CPU times presented in Tables 1–3.

Elements of the K matrix in other columns are obtained using this procedure with perturbations applied to the other displacement terms. Using this first-order approximation, the load distribution problem must be solved six times to define K entirely. Higher order finite difference approximation formulas can also be used similarly, but they would require solving for the spline load distribution many more times for different perturbations. In this study, only the first-order approximation will be used since this numerical method is intended only for verification of the fidelity of the analytical solutions of Sec. 2.1.

A moderate contact grid density represented by $Q=24$ (number of contact grid cells in face width direction) and $P=10$ (number of contact grid cells along profile direction) is used here for both the numerical and analytical methods for a direct comparison. A perturbation step of 10⁻¹⁰ m is selected for $δux$ and $δuy$⁠, and $10−10 rad$ for $δθx,δθy$⁠, and $δθz$ in the numerical method. Splines under three different operating conditions of (i) pure torsional loading, (ii) combined loading, and (iii) torsional loading with misalignments are considered for this comparison. The parameters of an example spline joint used in this comparison are listed in Table 4.

Tables 5 and 6 compare the K matrices of the example spline under pure torsion, obtained using the analytical and numerical methods. Here, $Mz=2260 N⋅m$⁠, and all coefficients of the stiffness matrix are given in SI units, e.g., units for $∂Fx/∂ux$⁠, $∂Fx/∂θx$⁠, and $∂Mx/∂θx$ are N/m, N/rad, and N · m/rad, respectively. It is seen that K defined by the analytical method is symmetric, as expected, since the spline joint is conservative, while the stiffness matrix predicted by the numerical method is not. In addition, the off-diagonal stiffness terms predicted by these two methods do not match well. These deviations can be attributed to perturbation around zero relative rigid body displacement terms, i.e., $ux=uy=0$ and $θx=θy=0$⁠, in this pure torsion loading condition. Regardless of this, diagonal stiffness terms of the matrix, which are often of primary concern, predicted by these two methods agree well. The relative difference of the five diagonal stiffness terms ranges between 1.7% and 6.6%. Note that the CPU time is 0.16 s for the analytical method, while it takes 6.77 s when using the numerical method. Further investigation of Table 6 also identifies a weakness of using the numerical method showing that some off-diagonal terms of the matrix are calculated as not only nonsymmetric but also in some cases of opposite sign. This is an effect of attempting to apply a perturbation to a solution that is nominally zero (i.e., $ux=uy=θx=θy=0$⁠). This effect is not present when applying the analytical method.

Tables 7 and 8 provide a similar comparison for the combined loading defined by $Mz=2260 N⋅m,$$Mx=606 N⋅m,$$Fx=−29.6 kN$⁠, and $Fy=11.2 kN.$ Again, K defined by the analytical method is symmetric while the numerical stiffness matrix is not. Compared to previous pure torsion case, the off-diagonal stiffness terms predicted by these two methods have much better correlation, since the relative rigid body displacement terms are no longer zero in this combined loading condition. The diagonal stiffness terms predicted by these two methods also appear very close to each other. The relative errors of the five diagonal stiffness terms are all within 4.8%. Similar observations are extended to the case of torsional loaded splines with misalignments as shown in Tables 9 and 10, where the torque $Mz=2260 N⋅m$ and the misalignment terms include $ux=uy=10μm$ and $θx=θy=0.01 deg$⁠. Also note that in both cases, the CPU time required by analytical method is much smaller than that required by the numerical method.

The above comparisons clearly demonstrate the reliability and computational efficiency of the analytical method for stiffness calculation. Considering also the analytical method is independent of parameters such as order of finite difference approximation and perturbation steps, this analytical method will be employed in the following parametric studies for a thorough investigation on stiffness of spline joints.

As they influence load distribution of splines, tooth surface modifications should be expected to impact the corresponding stiffness matrix as well. Using the same example spline joint of Table 4, influences of the tooth surface modifications (profile and lead modifications) on K are investigated here.

Figure 3 shows the variation of the diagonal torsional (⁠$∂Mz/∂θz$⁠), radial (⁠$∂Fi/∂ui,$$i=x,y$⁠), and tilting (⁠$∂Mi/∂θi,$$i=x,y$⁠) stiffness elements of K with torque for four levels of external tooth profile crown values between 0 (no modification) and 9 μm. Here, the spline is under pure torsion and there is no lead modification. The torque range considered is $Mz=565$ N · m to 5650 N · m. In this pure torsion loading case, the radial stiffness values along the x and y axes are equal as shown in Fig. 3(b). The same is true in Fig. 3(c) for the tilting stiffness values about the x and y axes. It is seen that in the case of zero profile crown modification, all the diagonal stiffness terms remain constant regardless of the value of $Mz$ because the contact area (conformal contact of every spline tooth engagement) remains the same at all torque levels under this pure torsion loading condition. However, in presence of positive profile crown modification, these stiffness terms are seen to increase gradually as $Mz$ increases. For instance, for a 9 μm profile crown modification, torsional stiffness values of the spline joint are 10.0, 10.2, 10.3, and 10.4 $N⋅m/μrad$⁠, respectively, for $Mz=1130$⁠, 2260, 3390, and 4520 N · m. The main reason for this is that more contact point pairs come into contact as the torque increases. Aside from this, it is also observed that at a given $Mz$ level, these diagonal stiffness terms decrease with an increase in profile crown modification. For example, at $Mz=2260$ N · m, $∂Mz/∂θz=11.1$⁠, 10.5, 10.3, and 10.2 $N⋅m/μrad$ for profile crowns equal to 0, 3, 6, and 9 μm, respectively. This can be explained by the fact that a larger profile modification causes more contact point pairs to separate.

Figure 4 shows variation of $∂Mz/∂θz$⁠, $∂Fi/∂ui$⁠, and $∂Mi/∂θi$ (⁠$i=x,y$⁠) with torque $Mz$ for four levels of external tooth lead crown values between 0 (no modification) and 15 μm (no profile modifications). It is observed that the influence of lead crown modification is quite similar to that of profile crown modification: (i) stiffness terms are load-independent for the no lead crown case; (ii) stiffness terms of lead crowned splines become larger as torque increases; and (iii) at a given torque level, stiffness terms decrease with larger lead crown modifications.

Misalignments of spline joints were shown to alter load distribution significantly [8,9]. Influence of misalignments on stiffness of spline joints is investigated here using the same example spline. A nominal 10 μm lead crown modification is applied to the external spline teeth in this case since it is common in application to reduce load concentration of misaligned splines using lead crown modifications. Figure 5 shows variation of the diagonal terms of K (⁠$∂Mz/∂θz,$$∂Fx/∂ux,$$∂Fy/∂uy,$$∂Mx/∂θx$⁠, and $∂My/∂θy$⁠) as a function of $Mz$ applied to the spline. Here, spline misalignment is varied from 0 deg to 0.03 deg with an increment of 0.01 deg. Similar to previous cases, for a given misalignment value, these diagonal stiffness values increase as $Mz$ is increased for the same reason that the number of contact point pairs is increased. Meanwhile, at a given $Mz$ value, the change of the diagonal stiffness terms with misalignment becomes more complicated. At higher torque levels, say $Mz>2260$ N · m in Fig. 5, these stiffness values reduce as the misalignment magnitude is increased. For instance, at $Mz=3390$ N · m, $∂Mx/∂θx=2.23,$ 1.97, 1.62, and 1.44 N · m/μrad for misalignments of 0 deg, 0.01 deg, 0.02 deg, and 0.03 deg, respectively. However, at low torque levels, say $Mz≤2260$ N · m, some of the stiffness terms become larger as misalignment increases. For example, at $Mz=565$ N · m, $∂Mx/∂θx=0.30,$ 0.45, 0.81, and 1.06 $N⋅m/μrad$ for misalignments of 0 deg, 0.01 deg, 0.02 deg, and 0.03 deg, respectively. These can be explained by the effects of the lead crown modification. At low torque, lead crown modification results in a load distribution concentrated around the middle of the face width of the spline teeth. This load distribution pattern is poor at supporting tilting moment since the point pairs in contact are so close to the middle of the face width. As misalignment increases slightly, load distribution slightly moves away from the middle of the face width so that tilting stiffness also increases slightly in these low torque cases. Conversely, at high torque, the influence of lead crown modification is minor since point pairs in contact spread over almost the whole face width. In this case, as misalignment increases, more contact point pairs separate in the meantime lowering tilting stiffness as well as other diagonal stiffness terms.

In the presence of indexing errors, a spline joint was shown in to exhibit time-varying (rotation-dependent) load distribution characteristics [11]. This should influence its stiffness matrix as well. In order to demonstrate this influence, load distribution analysis of the example spline having a random indexing error sequence is performed for a complete rotation of the spline (from 0 deg to 360 deg with an increment of 6 deg). Here, a torque of $Mz=4525$ N · m is applied in addition to a misalignment of $θx=0.04 deg$ about the x axis. The spline teeth in this example have a 5 μm profile crown modification and a 10 μm lead crown modification. The random indexing error sequence applied to the spline is shown in Fig. 6 along with the corresponding spline contact stress distributions at four representative rotational positions of 0 deg, 90 deg, 180 deg, and 270 deg. In agreement with Ref. [11], a spline having indexing errors experiences significant changes in its load distribution as a function of its rotation.

Stiffness computations for this example spline joint are performed at these 60 rotational positions. The corresponding diagonal stiffness terms at different rotational positions are shown in Figs. 7(a)–7(e). It is evident that as the spline rotates, all the diagonal stiffness terms change significantly, resulting in a parametrically varying K at the rotational period of spline. For instance, in Fig. 7(a), the value of $∂Mz/∂θz$ term varies between 6.89 and 7.84 N · m/μrad, representing about a 15% peak-to-peak fluctuation over a full rotational cycle. This can be explained by the fact that the total contact area over the spline varies at different rotational positions. For the load distributions shown in Fig. 6, the number of grid cells that are in contact is 1437, 1441, 1353, and 1464 at rotational positions 0 deg, 90 deg, 180 deg, and 270 deg, respectively, where the total number of grid cells is n = Z × P × Q = 25 × 10 × 24 = 6000. Since the contact zone is uniformly discretized into individual grid cells of equal areas, a larger number of grid cells in contact represent a larger total contact area along the tooth surfaces, and hence a higher $∂Mz/∂θz$ value. In Fig. 7(a), $∂Mz/∂θz=7.41$ 7.48, 7.13, and 7.71 N · m/μrad at rotational positions 0 deg, 90 deg, 180 deg, and 270 deg, respectively, this correlates well with the corresponding contact areas. The same rotation-dependence is evident in Figs. 7(b)–7(e) as well for the other diagonal terms $∂Fx/∂ux$⁠, $∂Fy/∂uy$⁠, $∂Mx/∂θx$⁠, and $∂My/∂θy$⁠, which exhibit even larger peak-to-peak fluctuations of 35%, 36%, 18%, and 50%, respectively. The above analysis indicates that the indexing errors of spline joints should be accounted for in high-fidelity drivetrain dynamic models.

A general analytical stiffness formulation of spline joints was proposed in this paper based on a semi-analytical spline load distribution model. This analytical stiffness formulation was verified through comparisons to a numerical perturbation method. The analytical method was then used to perform a detailed parametric study on the impacts of tooth modifications, loading conditions, misalignments, and tooth indexing errors. Profile modifications and lead modifications were both found to reduce the torsional, radial, and tilting stiffness of spline joints under pure torsion loading. Misalignments of spline joints were also observed to significantly influence the stiffness values. In the presence of random tooth indexing errors, stiffness terms of spline joints under asymmetric loading were shown to be time varying with considerable peak-to-peak variations as the spline rotates.