#### Abstract

By analyzing the roller force, the nonlinear stiffness model of the double-row tapered roller hub bearing is derived, and the method of solving the hub bearing stiffness matrix is summarized: if the displacement between the inner and outer rings is known, the stiffness of the hub bearing can be directly calculated. If only the external load of the hub bearing is known, the displacement of the hub bearing needs to be solved by numerical method, and then the stiffness of the hub bearing can be calculated. The improved Newton-Raphson method is used to solve the stiffness matrix of the hub bearing. Three-dimensional FE model of DRTRBs is presented and validated the proposed the stiffness matrix of the hub bearing. It is found that the radial stiffness of the hub bearing is greater than the axial stiffness. The stiffness of the hub bearing is greatly affected by the vertical force of the ground and the wheel driving torque, showing obvious nonlinearity. The smaller the vertical ground load and wheel driving torque, the greater the influence of vertical ground load and wheel driving torque on the hub bearing stiffness.

#### 1. Introduction

Double-row Tapered Roller Bearings (DRTRBs) are designed to support the preload, axial load, radial load, and torque in static and dynamic conditions [1, 2]. DRTRBs have been particularly considered in high-load supporting applications in automobile wheel hub assembly. Stiffness is the essential factors that significantly affect the performance of DRTRBs; bearing stiffness is needed for analyzing the dynamic characteristics, for example, natural frequencies, mode shapes, and vibration amplitudes [3–10].

The first few major works on rolling element bearing stiffness were performed by Palmgren [11], Jones [12], Dowson [13], and Harris [14], which were adopted in numerous bearing analyses [8, 15–18]. De Mul et al. [17, 18] computed analytically the bearing stiffness matrix of roller bearings based on the theory of Jones [12] and used internally in the iterative bearing equilibrium calculation. Tong and Hong [4, 19] presented time-varying stiffness matrix of tapered roller bearings (TRBs) based on the bearing stiffness matrix based on the theory of De Mul et al. [17, 18]. Based on the Hertzian theory, Lim et al. [3, 9, 10] derivated a time invariant bearing stiffness matrix of dimension 5. Kabus et al. [20] presented a new multidegree of freedom frictionless quasi-static time-domain TRBs model which allows for an arbitrary stiffness of the outer ring supporting structure. In a CAD environment, Bourdon et al. [21] developed an alternative method to estimate the stiffness matrix and discussed the influence of the static loads on the stiffness matrices of TRBs in a rigid environment. Guo and Parker [22] developed a finite element contact mechanics model and obtained accurate bearing stiffness for a wide range of bearing types (cylindrical and ball bearings) and parameters, and the stiffness matrix is calculated numerically through finite differences to first, second, fourth, and sixth order, respectively. However, the research papers did not covered the stiffness derivation of DRTRBs.

Though DRTRBs are widely used in industrial, the study on stiffness derivation of DRTRBs is sparse. For instance, Gunduz et al. [23] computed stiffness matrix for double-row angular contact ball bearings. However, inertial loading due to the rotational speed effect was not considered in their bearing model. Bercea et al. [24, 25] presented the relationship of bearing stiffness given by Harris [14] in a DRTRBs example, and formulated the relative displacement between the bearing rings for various double-row bearing types such as DRTRBs [26]. Fernandez et al. [2] optimized presetting loads to improve working conditions for a double-row TRB. In many static and dynamic problems, a double-row bearing must be considered as an integrated unit. The stiffness matrices of two single-row bearings may not be simply superposed to obtain the stiffness matrix of DRTRBs [23]. So the stiffness of DRTRBs is derived.

To overcome the deficiency in the literature, an approach about stiffness for DRTRBs is proposed based on the Hertzian theory for back-to-back arrangements. By analyzing the roller force, the nonlinear stiffness model of the double-row tapered roller hub bearing is derived. Then the improved Newton-Raphson method is used to solve the stiffness matrix of the hub bearing. The proposed stiffness model was validated by a proposed finite element model of double-row tapered roller hub bearing.

#### 2. Structure of Hub Bearing

The hub bearing studied in this paper is DRTRBs that support the semifloating drive axle. The axle housing is a steel plate stamping and welding structure. There are six bearings in the drive axle; four single-row tapered roller bearings (SRTRBs) in the final reduction drive and two DRTRBs are hub bearing.

As shown in Figure 1, DRTRBs are often arranged in a back-to-back or face-to-face arrangement generally. The DRTRBs studied in this paper are installed in a back-to-back arrangement. Generally, DRTRBs are formed by a pair of inner raceways that are separated by a gap of constant distance, an outer raceway, and a group of tapered rollers that are located between the inner and outer raceways. The stiffness calculation method of DRTRBs is different from that of SRTRBs.

**(a) Back-to-back arrangement**

**(b) Face-to-face arrangement**

#### 3. Nonlinear Stiffness Model of Hub Bearing

##### 3.1. Contact Deformations

The difficulty is to calculate contact mechanical properties between roller and inner and outer ring of bearing. To solve this problem, the elastic deformation and stress distribution of the contact surface should be analyzed firstly. Hertz’s elastic contact theory [27] and Sjovall’s integrals [28] is the basis for predicting the stiffness of rolling bearings. Hertz’s theory makes the following assumptions:

The contact body is an isotropic linear elastic material, which obeys Hooke’s law and is in a small deformation state.

The length* a* of the contact area is much smaller than the radius R of curvature of the surface of the object, that is,* a* ≪* R*.

The contact surface is smooth and continuous, and there is no friction.

Since* a* ≪* R*, every object can be regarded as an elastic half space.

For the elastic deformation of the line contact problem, there is generally no accurate analytical solution, and an empirical formula or an approximate formula is usually used. Palmgren [11] proposed an empirical formula for the linear contact elastic deformation, expressed by the relation

where denotes the elastic deformation of bearing roller, denotes the normal load in the contact area of the roller, and describes the length of the roller. Equation (1) shows that the roller deformation is independent of the curvature radius of the roller and inner and outer ring and only depends on the roller length and the normal load, so the elastic deformation of bearing roller is only an approximation. When the diameter of the bearing roller and the bearing ring changes, the elastic deformation of bearing roller does not change, and sometimes it differs greatly from the true value.

Houpert [29] proposed an empirical formula for the elastic deformation, which can be described as

Ding Chang’an [30] proposed an analytical solution to elastic line contact of finite length, taking into account the curvature radius of the contact surface; the elastic deformation can be expressed by

To calculate the stiffness of the rolling bearing, Wu Hao [31] proposed that the elastic deformation formula is expressed as

Luo Jiwei [32] verified the Palmgren’s empirical formula (1) and calculated the elastic deformation of the cylindrical roller with =10mm when it was in contact with the half plane. The result shows that the elastic deformation is correct only when the radius of cylindrical roller is 20 mm, and the elastic deformation will change greatly when the radius of the roller changes. Considering the diameter of the bearing roller and the inner and outer rings of the bearing, Luo Jiwei proposed the modified Palmgren formula, which can be described as

where denotes the elastic deformation of rollers. and denote the diameter and length of the roller, respectively. denotes the normal load in the contact zone of roller. denotes the diameter of the raceway. Considering the unevenness of the curvature of the contact surface, “+” denotes that the roller is in contact with the inner ring, and “-” denotes that the roller is in contact with the outer ring.

Based on the modified Palmgren formula, the force and stiffness model of tapered roller bearings is derived. The stiffness between the roller and the inner and outer rings is and , respectively. The elastic contact deformation between roller and inner and outer rings is and , respectively. The elastic contact deformation between roller and inner and outer rings can be interpreted as two springs in series [20]. Because the taper angle of the tapered roller is small, for the sake of simplicity of calculation, it is treated as a cylindrical roller. is the average diameter of the tapered roller. According to (5), the contact deformation between roller and inner and outer rings can be expressed by

The contact deformation between the inner ring and the roller is projected to the direction of the contact force of the outer ring with the reference of the contact deformation of the outer ring of the bearing. The projection of in direction is , so the total relative displacement of the inner and outer rings of the bearing in the direction,

where denotes the cone angle of the roller. Substituting (6) and (7) into (8), the total relative displacement can be expressed by

##### 3.2. Bearing Force

As shown in Figure 2, the geometric center of the bearing is the coordinate origin . The distance from the center of the roller to the origin of the coordinate is . The projections of in the radial and axial directions are and , respectively.

Referring to the stiffness analysis of SRTRBs, assuming that the outer ring of the bearing is fixed, the inner ring and the roller are taken as a whole. The deformation and force of the left rollers and the right rollers are analyzed, respectively.

###### 3.2.1. Force Induced the Left-Hand Rollers

For the displacement components , the axial displacement , and the rotation angle of the inner ring relative to the outer ring, the radial displacement of the jth roller on the left side of the bearing (azimuth angle ) due to rotation angle is expressed by

The total radial displacement of the jth roller on the left-hand bearings can be expressed by

The axial displacement of the jth roller on the left-hand bearings (azimuth angle ) due to rotation angle is expressed by

The total axial displacement of the jth roller on the left-hand bearings can be expressed by

Since the radial displacement includes displacements and along* x* axis and* y* axis, the rotation angle includes rotation angles and around* x* axis and* y* axis, The total radial displacement and axial displacement of the jth roller on the left-hand bearings can be expressed by

The additional moment caused by the axial component of the roller contact force can be expressed by

The additional moment caused by the radial component of the roller contact force can be expressed by

The directions of force and displacement refer to the coordinate system in Figure 2, and the directions of rotation and torque obey the right-hand rule. The bearing forces and moments induced by the left-hand rollers can be expressed by

###### 3.2.2. Force Induced the Right-Hand Rollers

The above analysis is only for the left-hand rollers. The forces on the left-hand rollers and the right-hand rollers are different. The absolute value of the displacement and load of each roller on the right side is the same as the corresponding left roller, but the directions are not exactly the same. For example, when the inner ring is subjected to the axial displacement , the left-hand roller is close to the outer ring, and the right-hand roller is away from the outer ring. Then the contact force between the left-hand roller and the outer ring increases, while the contact force between the right-hand roller and the outer ring decreases. Referring to the force analysis of the left roller, the radial and axial displacements of the jth roller on the right can be expressed by

The bearing forces and moments induced by the right-hand rollers can be expressed by

###### 3.2.3. Bearing Force and Moments

According to (17) and (19), the force and moment in all directions of the double-row roller bearing can be expressed by

where and are the normal loads of the jth roller on the left-hand bearings and the right-hand bearings, respectively. According to (9), (14), and (18), the normal displacement of each roller can be expressed by

Substituting (21) into (9), the normal load between each roller and the outer ring can be expressed by

In (21) and (22), the directions of force and displacement refer to the coordinate system in Figure 2, and the directions of rotation and moments obey the right-hand rule. When the outer ring is fixed and the inner ring is subjected to force ( =* x, y, z*) and moments (*i* =* x, y, z*), the displacement and rotation angles of the inner ring are (*i* =* x, y, z*) and (*i* =* x, y, z*), respectively.

#### 4. Solution of Hub Bearing Stiffness Matrix

##### 4.1. Bearing Stiffness Matrix

Bearings are often simplified to spring elements with axial and radial stiffness in the finite element modeling of drive axle assembly. In fact, the bearing has five directions of stiffness except for the degree of freedom around the* z* axis, and the stiffness in each direction is coupled. Therefore, the bearing stiffness matrix can be more accurately expressed by

Diagonal terms in the stiffness matrix include radial stiffness , axial stiffness , and angular stiffness , where* i*=*x, y*. The off-diagonal cross-coupling terms in the stiffness matrix are the coupling stiffness values. The local coordinate system of the bearing is shown in Figure 2. The bearing is free to rotate about the* z* axis, so the sixth row and column in the stiffness matrix are zero. The stiffness matrix is symmetric because the bearing system is conservative. The load and displacement equations of bearings can be expressed by

where load vector is and displacement vector is .

##### 4.2. Force on Hub Bearing

Hub bearings studied in this paper are DRTRBs. The structure of the left and right hub bearings and the applied load are exactly the same. DRTRBs support the semifloating drive axle as shown in Figure 3. That is, the rear drive shaft is subjected to the moments transmitted from the vehicle transmission system and the vertical load of the wheel. The bearing is subjected to vertical and driving force from the road, as well as the bending moments induced by them.

To coincide with the coordinate system shown in Figure 2, the vertical force and the driving force of the road surface are represented by and , respectively. The force and moments of the hub bearing can be expressed by

where denotes the rear axle load. denotes the wheel driving torque.* r* denotes the tire rolling radius.* b* denotes the distance from the center of the wheel to the center of the hub bearing.

##### 4.3. Solution of Stiffness Matrix

The stress of hub bearing is relatively simple, so it is not necessary to adopt the method of finite element iteration. According to the external loads on the bearing, the displacement and the rotation angle between the inner and outer rings of the bearing are obtained. Then, the stiffness value of the bearing can be calculated by substituting and into (23).

According to (20) and (25), the relation of the displacement and rotation angle between inner and outer rings under external loads can be expressed by

where , external Loads Vector , displacement vector , bearing stiffness matrix is* K*(5×5), and then (26) can be expressed by

Equation (26) is a set of nonlinear equations for five variables ,,, , and . It is almost impossible to solve the exact solution by analytical method. In this paper, the improved Newton iteration method is used to solve this problem.

Newton-Raphson method is also called the tangent method. By continuously making the tangent of the curve and finding the coordinates of the intersection of the curve and the* x* axis, the solution of the equation can be obtained. The iterative formula of the one-dimensional equation can be expressed by

In (28), the selection of initial value will have a great impact on the convergence of the algorithm. This paper uses the improved Newton-Raphson method to calculate the displacement and stiffness of the hub bearing. The improved Newton-Raphson method of the one-dimensional equation can be expressed by

For multidimensional equations, the improved Newton-Raphson method can be expressed by

It can be seen from (30) that the improved Newton-Raphson method divides the iteration process into two steps. Firstly, the intermediate variable is obtained by evaluating . Secondly, is substituted into Newton-Raphson method, which can reduce the influence of initial value deviation on convergence. In order to improve the convergence, combined with Newton downhill method, the iterative formula can be expressed by

where ) is the downhill parameter or relaxation factor. The downhill condition can be expressed by

The downhill condition is satisfied by adjusting the magnitude of adaptively. When the downhill condition is satisfied, it is the classical Newton iteration method. When the downhill condition is not satisfied, is reduced by half until the downhill condition is satisfied. The calculation process of hub bearing displacement and rotation angle based on the improved Newton-Raphson method is shown in Figure 4.

#### 5. Finite Element Model Validation

A numeric model of the ball bearing is established using ABAQUS [33, 34]. In order to validate the method of bearing stiffness, three-dimensional FE model of DRTRBs is presented using ABAQUS, which is shown in Figure 5. Because of the nonlinearity of mechanical contact in the use of FEM, contact stresses are unrealistic when the FE models have limited contact surface and the size of mesh is very large. This problem is usually solved by making the mesh between contact pairs small, i.e., increasing the number of contact nodes. Demirhan and Kanber studied the effect of FE model element size on contact stress between the raceway and rollers [35]. Fernandez et al. used small mesh size to obtain fine FE models with contact areas. By using small mesh sizes of contact areas, the contact stress obtained by finite element method is not much different from that by theoretical analysis and experimental method [2]. Considering the central processing unit (CPU) running time increases exponentially with the element dimension decreases, mesh refinement on the contact areas of the rollers shown in Figure 6 was applied. To ensure high computation accuracy of contact stresses, the finite-element mesh is generated C3D8R elements, and the total number of elements is 1,757,080 using Hypermesh software. Raceways and tapered rollers are modeled using linear elastic and isotropic steel. The geometry and material parameters of the roller bearing are given in Tables 1 and 2.

**(a) Between the rollers and the inner raceway**

**(b) Between the rollers and the outer raceway**

The outer surface of the outer raceway is developed as fixed surface and static radial force was applied in the center of the bearing. A reference point is established in the center of rotation of the inner ring, the part of the inner surface of inner ring is coupled with the reference point. The interaction between surfaces is defined as surface-to-surface contact with no adjustment. Interaction properties are imposed using surface-to-surface hard normal contact between the master surface of roller and the slave surface of the outer raceway as well as between the master surface of roller and the slave surface of the inner raceway. The coefficient between the rollers and the inner raceway and between the rollers and the outer raceway is 0.001 [35].

The bearing stress distribution of FEM can be obtained when the radial force was applied. Figure 7 shows the stresses distribution on the outer raceway and inner raceway caused by the first and second row of rollers, which is in good agreement with the stresses distribution studied by Harris [36] and Lostado et al. [1] for DRTRBs. Figure 8 shows the variation of the radial force and the radial stiffness with respect to radial displacement of the bearing. Figure 9 shows variation of the moment and the tilting stiffness with respect to angular displacement of the bearing; Figures 8 and 9 are in good agreement with the load-deflection relationship that was studied by Harris [36] and Kania [37] for DRTRBs. The largest error shown on the Figures 8 and 9 is less than 10%. FEM verification of the proposed stiffness model is available.

**(a) The outer raceway**

**(b) The rollers**

#### 6. Results and Discussions

When the rear axle is fully loaded with 16kN and the drive torque is 200Nm. The relevant dimensions of the chosen hub bearing are listed in Table 2. The hub bearing displacement vector is obtained by the above method, . The five-dimensional stiffness matrix of hub bearing can be obtained by (20), (21), and (23). The stiffness matrix of hub bearing is expressed by

Equation (33) shows that the tilting stiffness is far greater than the stiffness in other directions. This is mainly because the relative rotation of* x* axis and* y* axis occurs between the inner and outer rings of the double-row roller, and the radial and axial components of the contact force between the roller and the outer ring cause a resistance moment around the center of the bearing, which is related to the structure of the DRTRBs.

Figure 10 shows the relationship between the hub bearing stiffness and the vertical load on the ground. When the bearing is only subjected to the vertical load, the stiffness of the bearing along the* x* axis and* y* axis is much greater than the stiffness along the* z* axis; that is, the radial stiffness of the bearing is much larger than the axial stiffness, and the relationship between the stiffness magnitudes is* Kx*>*Ky*>*Kz*. The stiffness around the y axis is slightly greater than that around the x axis because the vertical load of the drive wheel induces the moment around the* y* axis. When the vertical load on the ground is less than 2000N, with the increase of the vertical load on the ground, the bearing stiffness shows a significant nonlinear increase. When the vertical load is small, the bearing stiffness increases rapidly, and as the load continues to increase, the increased amplitude of bearing stiffness decreases gradually. When the vertical load exceeds 2000N, the bearing stiffness increases linearly with the increase of the vertical load. The vertical load variation on one wheel of the light passenger car studied in this paper ranges from 4000N to 8000N (no load to full load), and the bearing stiffness is basically within the linear variation range.

Figure 11 shows the relationship between the diagonal element of the bearing stiffness matrix and the wheel driving torque when the hub bearing is only subjected to the driving torque. The radial stiffness of the bearing is much larger than the axial stiffness, which is consistent with the bearing stiffness under the vertical load. The relationship between the magnitudes of radial stiffness and axial stiffness is* Ky* >* Kx* >*Kz*. The bearing stiffness around the* x* axis is significantly greater than the bearing stiffness around the* y* axis. This is because the ground driving torque acting on the wheel causes the bearing to generate a torque around the* x* axis, so the stiffness around the* x* axis is larger. When the driving torque of the wheel is less than 200Nm, the bearing stiffness increases nonlinearly with the increase of the driving torque, and the increasing amplitude of bearing stiffness gradually decreases. When the driving torque exceeds 200 Nm, the diagonal element of the bearing stiffness matrix increases linearly with the increase of the driving torque.

#### 7. Conclusions

In this paper, the stiffness matrix expression of double-row tapered roller hub bearings is derived, and the stiffness matrix is calculated by improved Newton-Raphson method. The conclusions are as follows:

The forces acting on the left rollers and the right rollers of DRTRBs are different. When calculating the overall stiffness of the hub bearing, the force acting on each roller should be considered separately.

If the displacement between the inner ring and outer ring has been determined, the stiffness of the hub bearing can be directly calculated. If only the external loads of the hub bearing have been determined, the displacement of the hub bearing should be calculated by numerical method, and then the stiffness of the hub bearing should be calculated.

The torsional stiffness of hub bearings around* x* axis and* y* axis is far greater than the stiffness in other directions. The radial stiffness of the hub bearing is much larger than the axial stiffness.

Affected by vertical load on the ground and wheel drive torque, the stiffness of the hub bearing exhibits significant nonlinearity. The smaller the vertical ground load and wheel driving torque, the greater the influence of vertical ground load and wheel driving torque on the hub bearing stiffness.

#### Data Availability

Climate and demand data used to support the findings of this study are available from the first author upon request.

#### Conflicts of Interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

#### Acknowledgments

The research work was supported by National Key R&D Program of China (Grant no. 2018YFB0106203).