# A Critical Analysis of Finite-Element Modeling Procedures for Radial Bearing Stiffness Estimation

Authors: Alberto Gabrielli; Mattia Battarra; Emiliano Mucchi
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Different strategies are commonly employed by researchers in order to decrease the computational effort associated with the finite-element analysis of rolling-element bearings. The purpose of this work is to review and analyze the procedures and hypotheses that may be exploited to evaluate the nonlinear radial stiffness of these components. Techniques are utilized to develop a meshing procedure aimed at balancing the computational effort and the accuracy of the results, to define a robust approach to the problem. The geometry is reduced by taking advantage of the available symmetry planes, by removing unloaded rollers, and by substituting the shaft with an equivalent sinusoidal load. In addition, the element dimensions are adapted to the applied load depending on the size of the contact area as computed by means of the Hertz theory. The proposed methodology may be applied to all bearing types provided that symmetry conditions and contact area dimensions are properly assessed. The estimated stiffness is compared against analytical formulae retrieved from the literature. Influence of different element types, roller position, cage, and clearance on accuracy and computational time is discussed.


---

## Body

## 1. Introduction

Rolling-element bearings (REBs) are one of the most frequently employed components in rotating machinery, where they cover a major role in transmitting vibrations [1]. Due to their relevance, REBs have been extensively studied, and vast literature regarding their characteristics and behavior is currently available. Nonetheless, the inherent complexity of these components still provides challenges to overcome for researchers and designers, such as the evaluation of their radial stiffness. This is usually a crucial input parameter for further analyses, e.g., in the dynamic modeling of systems containing bearings [2].Radial stiffness is not a constant value, as the relative displacement of the inner ring with respect to the outer ring depends on the applied load and the position of the rolling elements with respect to its direction. Over the years, the estimation of this nonlinear stiffness-load relationship has been faced by means of experimental, analytical, and numerical approaches. Experimental techniques are commonly divided into direct and indirect methods, depending on the employed procedure. For the former, the direct measurement of the displacement is required [3]. The latter requires other techniques, e.g., modal analysis [4]. Analytical formulae may be derived empirically [5, 6] or by providing a rigorous mathematical procedure typically based on the Hertz contact theory [7–9]. Among the numerical techniques, the most common approach involves the employment of in-house or commercial finite-element (FE) software to model the bearing under exam.FE simulations are usually exploited to determine rings displacement and contact stresses. A common issue is the complexity of the contact phenomenon, which leads to challenging modelization and large computational effort. Thus, researchers strived to simplify the problem using a variety of approaches. The most straightforward approach consists in reducing the original 3D model to a 2D problem. Within this framework, Zhao [10] described a contact algorithm to model roller bearing contacts in 2D simulations. Influence of various parameters on load distribution was evaluated. Demirhan and Kanber [11] investigated stress and displacement distribution on a roller bearing. A 2D mesh with plane strain option was used to simplify the model. Obtained results were validated against theoretical and experimental data. Hao et al. [12] proposed a 2D FE model that considered temperature effects and clearance change to examine their influence on displacement, stress, and bearing stiffness. The relative displacement between inner and outer ring obtained on a test rig was employed for validation. Good agreement was found between FE and experimental results, while analytical formulae departed from them. For 3D simulations, other methods were established. Kania [13] employed truss elements to replace rollers in simulating slewing bearings. This technique allowed great reduction of simulation time at the expense of additional effort in the preprocessing phase. A similar approach has been utilized by Daidié et al. [14] to analyze load distribution by employing nonlinear traction springs instead of truss elements. The efficiency of this method encouraged other researchers [15–18] to exploit it for slewing bearings modeling. Techniques devoted to removing contacts in lower size bearings were also studied by Molnar et al. [19]. In particular, two methods to replace contact algorithms between rollers and rings in needle bearings have been proposed. The first one is similar to the one employed by Daidié et al. [14], as rollers are replaced by springs. The second one consisted in substituting the entire volume between rings with a fictitious equivalent material replicating the same behavior of the row of rollers. Obtained radial displacements were compared with the full model involving contacts, demonstrating their capability to considerably reduce the computational time. Although these techniques were successful in speeding up the solution process, other researchers followed different strategies to model the phenomena occurring within the contact area. Guo and Parker [20] developed a procedure involving a combined surface integral and finite-element method to solve the contact problem in rolling and ball bearings. They computed radial, axial, and tilting stiffness to obtain a fully populated 6 × 6 matrix including cross-coupling terms. Results were compared against data available from the literature. Massi et al. [21] set up 2D and 3D simulations of ball bearings to compute the contact stress due to specific boundary conditions and relate them to bearing degradation. They reduced the size of the problem by modeling a portion of the bearing, with only one roller in contact with the two races, noticing that the relative error between 3D and 2D simulations was relevant due to conforming contact between ball and races. Lostado et al. [22] studied the contact stress in tapered roller bearings. They developed a procedure to adjust the original mesh by generating subsequent nonlinear submodels with increasingly smaller mesh densities. Relative displacement between raceways was also analyzed and compared to experimental data, showing good agreement. The procedure was utilized by Fernandez Martinez et al. [23] in combination with machine learning techniques to determine the optimal working conditions of the device. Murer et al. [24] presented a FE model of their experimental setup to assess the relevance of using capacitive probes for in situ measurements of bearing deflection. Li [25] developed software to compute the contact stress in ball and roller bearings by exploiting a novel contact algorithm. Stress distribution on contact areas was found to be different from results reported in previous studies [20] and analytical formulae.As it may be noticed from the literature review, there is a tendency to reduce the size of the computational domain of the problem whenever REBs modeling is involved. Typical strategies involve taking advantage of symmetries [10, 13–15, 18, 19, 22, 23], removing unloaded rollers [10], and replacing contacts with equivalent elements [13–18]. As far as the authors are concerned, only one symmetry plane is employed in 3D analysis, although for most bearings two symmetry planes are available under appropriate hypotheses. Unloaded rollers are kept by most researchers but removed by others, e.g., Zhao [10] and Murer et al. [24]. Cage is commonly neglected, but, e.g., Murer et al. [24] accounted for its effect by employing rigid connectors between rollers. Load is applied on the center of the shaft [10–12, 19], on rings [14, 18, 21], or on a central node connected with rigid elements to the inner ring [15–17]. Problems are solved by employing 2D and 3D approaches, usually exploiting quadrilateral or hexahedral elements, respectively. Meshing with tetrahedral elements is rare, and differences between linear and parabolic elements are not addressed. Convergence check is regularly performed, but some researchers such as Demirhan and Kanber [11] only tested it for one load value, although convergence rate depends on applied load, especially at low force values. Concerning the post-process, radial bearing stiffness may be computed considering the approach of bearing rings [11, 12, 22] or the displacement of the shaft axis [19].While every reviewed method has proved to be successful in different aspects of bearing analysis, there is a lack of uniformity in the employed approaches. Based on this observation, the purpose of this work is to review and analyze some modeling choices that may be utilized by bearing analysts to speed up and improve simulations involving REBs, in particular in scenarios where a reasonable estimation of the radial stiffness is needed. It is the authors’ opinion, in fact, that element choice, mesh size, exploited symmetries, load definition, and displacement evaluation location are not always accurately defined in research papers. While this work should not be considered by any means a comprehensive resource for bearing analysis, the aim is to provide guidelines and define a robust approach to this problem. In the light of this, only key points that have been considered of primary importance will be discussed, while more advanced topics are not covered here for the sake of brevity and clarity. By optimizing the considered aspects, a fast and efficient method for generating REBs mesh via FE simulations is detailed. The procedure involves the generation of a dedicated mesh for each load condition, in which the element size is determined based on the estimated contact area dimensions obtained by means of the Hertz contact theory. This methodology allows the estimation of reasonable dimensions for the element in the contact area, so that the contact phenomenon is captured while balancing the computational time and the accuracy of the results. The attained mesh should be utilized to determine radial bearing stiffness only, as the convergence check is performed only in terms of displacement of the inner ring.A roller bearing is taken as reference for the analyses. Despite the fact that the computations are carried out for this particular bearing only, the proposed methodology is still applicable for other bearings types. However, attention must be paid when analysing different geometries, as there could be only one symmetry plane available and not two as in the reference bearing, e.g., in tapered roller bearings and self-aligning double-row ball bearings. Concerning element type, mesh size, unloaded roller removal, and load application method, the proposed procedure may be adapted to all bearing types. Nonetheless, it is worth noting that for ball bearings different formulae must be employed when determining the size of the contact area, since its shape is an ellipse and not a rectangle as in the roller-raceway contact.The following section outlines general aspects regarding the simulations described throughout this paper. Section3 is focused on convergence analysis for different element types. Section 4 describes the procedures exploited to reduce the size of the computational domain, along with stiffness results determined by means of analytical and numerical methods. In addition, results of simulations considering different roller positions, inclusion of the cage and presence of diametral clearance are reported. Finally, Section 5 is devoted to concluding remarks.


## 2. General Considerations

### 2.1. Reference Geometry

A roller bearing, modelSKF NU 202 ECP, is taken as reference for the analyses. Figure 1 shows the 3D geometry of the mechanical component, whose geometrical data are reported in Table 1. The reported values refer to nominal dimensions, as clearance is not included in the reported CAD model. Rollers have a straight profile and are 6 mm wide. Since edges are rounded with a 0.2 mm radius, their effective length reduces to 5.6 mm, which is the length of the ideal contact line between roller and races. A 0.1 mm axial clearance between rollers and flanges is considered. For all components, material is steel with Young’s modulus E=210 GPa, density ρ=7760 kg/m3, and Poisson’s ratio ν=0.28. These properties are considered as constant values throughout the analyses.Figure 1
Geometry of the roller bearing taken as reference case study.Table 1
Roller bearing design parameters.
QuantityValueDescriptionZ11Number of rollersDr5.5 mmRoller diameterL6 mmTotal roller lengthLeff5.6 mmEffective roller lengthdi15 mmInner ring internal diameterdo19.3 mmInner ring external diameterDi30.3 mmOuter ring internal diameterDo35 mmOuter ring external diameterW11 mmRing width


### 2.2. Modeling Hypotheses and Considered Effects

Simulations are carried out by using Simcenter 3D as pre/post-processor and Simcenter Nastran as solver. In this context, it is worth clarifying that Nastran static solution SOL 101 is used to solve the model since it allows considering the nonlinear behavior given by the contact algorithm and concurrently reducing the computational burden by assuming a linear elastic material and small displacement.Concerning the contact algorithm, two methods are available, namely, segment-to-segment method and node-to-segment method. Previous work by El-Abbasi and Bathe [26] indicated that, while both methods provided stable results, the latter did not pass the patch test [27], leading to discretization errors that did not decrease with mesh refinement. In addition, proof of the successful employment of the segment-to-segment method in bearing simulation may be found in [22, 28]. These results support the choice of the first approach for REBs simulation. Bilinear Coulomb friction model is chosen to take into account friction-related contributions. Friction coefficient is set to 0.05, which indicatively corresponds to a greased contact condition [29]. This value is chosen for the sake of simplicity to represent a common value encountered when grease is employed for bearing lubrication. However, in case of particular lubricating conditions, more accurate methods may be utilized to determine the coefficient of friction and the resultant friction force, e.g., exploiting the methods reported in [30].The proposed investigation exclusively takes into account the major phenomena concurring with the determination of the bearing radial stiffness. Based on this approach, the hypotheses at the basis of the present study are hereinafter detailed.The first aspect concerns possible plasticization of the material in the neighborhood of the contact areas. Bearing races, in fact, permanently deform at sufficiently high loads, affecting the wear and degradation process of the races [21, 29]. In this work, only deformations in the elastic range are considered by employing a linear elastic material with undefined yield stress.Moreover, as this paper focuses on radial stiffness, estimation methods for other components of the bearing stiffness matrix will not be described, although axial and tilting stiffness may be important when out-of-plane vibrations are dominant [20]. In this scenario, it is suggested that the full 6 × 6 stiffness matrix including cross-coupling terms is determined.Particular attention should be paid to boundary conditions. In this work, the external outer race of the bearing is assumed to be connected to a rigid frame. As a matter of fact, the adopted hypothesis represents a common scenario where the frame is sufficiently rigid not to interfere with the bearing properties, and the applied loads stay within the elastic deformation range. For the sake of completeness, it has to be clarified that actual applications might involve the installation on compliant frames or excessive loads leading to large deformation of the outer race [8, 22, 31]. In such cases, shaft misalignment may appear [32, 33], causing a significant moment load on the bearing and the loss of one symmetry plane from the system. These peculiar conditions require dedicated analyses involving different solution schemes, and therefore they do not fall within the purpose of the present study.Furthermore, thermal effects and preloading are not considered. As a consequence, simulations are run considering constant properties at room temperature and applying a simple radial load only. It is worth noting, however, that temperature may affect stress distribution and stiffness [12], while preload is commonly employed to modify bearing stiffness characteristics [34].Finally, it should be noted that static methods are employable at low and moderate speed values only [20, 35] as inertia effects cannot be neglected at higher velocities [36].The described hypotheses underline that several additional aspects not included in the present work may influence the actual REBs stiffness. Nevertheless, they may be implemented starting from the concepts and insights given in this paper. It is also worth noting that the proposed analyses may be further extended, with appropriate modifications, to other bearing types. In fact, the concepts and methods hereinafter explained apply for all bearing types, although some modifications have to be considered. For instance, the main difference between ball and roller bearings resides in the different contact area generated by ball-raceway contact. As a result, only data obtained for the reference roller bearing are compared in this work.


## 2.1. Reference Geometry

A roller bearing, modelSKF NU 202 ECP, is taken as reference for the analyses. Figure 1 shows the 3D geometry of the mechanical component, whose geometrical data are reported in Table 1. The reported values refer to nominal dimensions, as clearance is not included in the reported CAD model. Rollers have a straight profile and are 6 mm wide. Since edges are rounded with a 0.2 mm radius, their effective length reduces to 5.6 mm, which is the length of the ideal contact line between roller and races. A 0.1 mm axial clearance between rollers and flanges is considered. For all components, material is steel with Young’s modulus E=210 GPa, density ρ=7760 kg/m3, and Poisson’s ratio ν=0.28. These properties are considered as constant values throughout the analyses.Figure 1
Geometry of the roller bearing taken as reference case study.Table 1
Roller bearing design parameters.
QuantityValueDescriptionZ11Number of rollersDr5.5 mmRoller diameterL6 mmTotal roller lengthLeff5.6 mmEffective roller lengthdi15 mmInner ring internal diameterdo19.3 mmInner ring external diameterDi30.3 mmOuter ring internal diameterDo35 mmOuter ring external diameterW11 mmRing width


## 2.2. Modeling Hypotheses and Considered Effects

Simulations are carried out by using Simcenter 3D as pre/post-processor and Simcenter Nastran as solver. In this context, it is worth clarifying that Nastran static solution SOL 101 is used to solve the model since it allows considering the nonlinear behavior given by the contact algorithm and concurrently reducing the computational burden by assuming a linear elastic material and small displacement.Concerning the contact algorithm, two methods are available, namely, segment-to-segment method and node-to-segment method. Previous work by El-Abbasi and Bathe [26] indicated that, while both methods provided stable results, the latter did not pass the patch test [27], leading to discretization errors that did not decrease with mesh refinement. In addition, proof of the successful employment of the segment-to-segment method in bearing simulation may be found in [22, 28]. These results support the choice of the first approach for REBs simulation. Bilinear Coulomb friction model is chosen to take into account friction-related contributions. Friction coefficient is set to 0.05, which indicatively corresponds to a greased contact condition [29]. This value is chosen for the sake of simplicity to represent a common value encountered when grease is employed for bearing lubrication. However, in case of particular lubricating conditions, more accurate methods may be utilized to determine the coefficient of friction and the resultant friction force, e.g., exploiting the methods reported in [30].The proposed investigation exclusively takes into account the major phenomena concurring with the determination of the bearing radial stiffness. Based on this approach, the hypotheses at the basis of the present study are hereinafter detailed.The first aspect concerns possible plasticization of the material in the neighborhood of the contact areas. Bearing races, in fact, permanently deform at sufficiently high loads, affecting the wear and degradation process of the races [21, 29]. In this work, only deformations in the elastic range are considered by employing a linear elastic material with undefined yield stress.Moreover, as this paper focuses on radial stiffness, estimation methods for other components of the bearing stiffness matrix will not be described, although axial and tilting stiffness may be important when out-of-plane vibrations are dominant [20]. In this scenario, it is suggested that the full 6 × 6 stiffness matrix including cross-coupling terms is determined.Particular attention should be paid to boundary conditions. In this work, the external outer race of the bearing is assumed to be connected to a rigid frame. As a matter of fact, the adopted hypothesis represents a common scenario where the frame is sufficiently rigid not to interfere with the bearing properties, and the applied loads stay within the elastic deformation range. For the sake of completeness, it has to be clarified that actual applications might involve the installation on compliant frames or excessive loads leading to large deformation of the outer race [8, 22, 31]. In such cases, shaft misalignment may appear [32, 33], causing a significant moment load on the bearing and the loss of one symmetry plane from the system. These peculiar conditions require dedicated analyses involving different solution schemes, and therefore they do not fall within the purpose of the present study.Furthermore, thermal effects and preloading are not considered. As a consequence, simulations are run considering constant properties at room temperature and applying a simple radial load only. It is worth noting, however, that temperature may affect stress distribution and stiffness [12], while preload is commonly employed to modify bearing stiffness characteristics [34].Finally, it should be noted that static methods are employable at low and moderate speed values only [20, 35] as inertia effects cannot be neglected at higher velocities [36].The described hypotheses underline that several additional aspects not included in the present work may influence the actual REBs stiffness. Nevertheless, they may be implemented starting from the concepts and insights given in this paper. It is also worth noting that the proposed analyses may be further extended, with appropriate modifications, to other bearing types. In fact, the concepts and methods hereinafter explained apply for all bearing types, although some modifications have to be considered. For instance, the main difference between ball and roller bearings resides in the different contact area generated by ball-raceway contact. As a result, only data obtained for the reference roller bearing are compared in this work.


## 3. Element Choice and Mesh Size

Analyses are performed to quantify the influence of adopting different mesh grids, with a specific focus on hexahedral and tetrahedral elements. Comparisons are carried out in terms of stiffness and associated computational time. The goal of the investigation is to determine an appropriate mesh size that allows attaining a reasonable compromise between the accuracy of the results and the computational effort. Within this framework, particular attention is devoted to the development of a methodology to determine a reasonable element dimension based on an analytical formulation, in order to reduce the time needed to check mesh convergence in further analyses. The problem with mesh size, in fact, is that the element dimensions in the neighborhood of the contact area should be sufficiently small to properly model the contact phenomenon. According to Hertz theory [37], the contact area for a roller of finite length Leff in contact with a curved surface and loaded with a force F is a rectangle of length Leff and half-width a. Let the contacting bodies be labeled with 1 and 2. Req is their equivalent contact radius:(1)Req=1R1+1R2−1,where R1 and R2 are the radii of curvature of the contacting surfaces. They may take positive or negative values if the surfaces are convex or concave, respectively.An equivalent Young’s modulusE∗, depending on material properties of both components, may also be defined as(2)E∗=1−ν12E1+1−ν22E2−1.ParameterE∗ is needed to compute the indentation depth d due to an applied load F:(3)d=4FπLeffE∗,which may be eventually utilized to determine the contact area half-width a as(4)a=Reqd.According to (4), the contact area decreases as the load reduces, leading to the need for smaller elements at lower loads. Moreover, contact area half-width is wider for roller-outer ring contact than for the inner ring contact, as the diameter of the outer raceway is larger. Within this framework, it has to be noted that Hertz theory does not account for the edge effects caused by the finite length of the components [17]. However, this theory is employed as it provides straightforward formulae that may be used for a rough estimation of the contact area. It is also worth noting that a rectangular contact area is produced only if the roller profile is straight. For fully or partially crowned profiles, a different formulation should be employed to estimate the contact area.For the reference bearing with parameters reported in Table1, an operating load range within Fr=0 N and Fr=10 kN is assumed based on its basic static load rating C0=10.2 kN. Fr is the radial load applied on the shaft and transmitted to the inner ring. The corresponding force on the most loaded roller may be computed, in absence of diametral clearance, by using the formula provided by Harris [8]:(5)Qmax=Fr0.2453Z,leading to a value Qmax=3706 N when Fr=10 kN. Figure 2 shows contact areas half-width for inner and outer ring contact with the most loaded roller.Figure 2
Contact area half-width for roller-raceway contacts.


### 3.1. Procedure for Mesh Size Performance Evaluation

Convergence analyses are carried out to determine a reasonable mesh size that allows obtaining a good estimation of the radial bearing stiffness while maintaining a relatively low computational time. An angular sector of the bearing is taken as reference in order to evaluate elements with progressively reduced dimensions. Geometry is cut to consider only half of the contact zone, whose length isLeff/2. Thus, only a quarter of the roller in contact with a portion of the two races is modeled, as the system is symmetrical with respect to two planes. The angle span of the reduced model, then, decreases to Δϕ=360/2Z=16.36∘. Figure 3 shows the dimensions of the adopted angular sector. The load is applied by providing a sinusoidal load distribution on the inner ring, which resembles the force exchanged with the shaft. Bodies are meshed with either tetrahedral or hexahedral elements, namely, CTETRA and CHEXA, in Nastran environment. For both elements, linear and parabolic formulations are available and identified by a number after their name, defining their number of nodes. As a consequence, element CTETRA4 indicates a linear tetrahedral element with four nodes, while CTETRA10 defines a parabolic element with ten nodes. Similarly, element CHEXA8 indicates a linear hexahedral element with eight nodes, while CHEXA20 defines a parabolic element with twenty nodes.Figure 3
Grids of the reduced geometry employed for convergence testing: (a) tetrahedral elements; (b) hexahedral elements and dimensions of the adopted angular sector.
(a)(b)Based on the necessity to determine the most convenient element, a procedure to determine element performance is established. The proposed method is applicable to all tested element types, with only few modifications needed to switch from tetrahedral to hexahedral elements. Convergence analysis is performed by decreasing the mesh size, which is quantified as the average length of element edges, so that the total number of nodes of the model progressively increases. Mesh dimension in the contact area,l, is smaller than the mesh size of the rest of the model, r. As the local dimension is reduced, r is kept proportional to l by means of a constant coefficient α so that l=αr. The size of the contact area is kept as a constant term for all load values. Its half-width is taken to be w=0.25 mm, which is a conservative value obtained through (4). This area must be geometrically defined because, for the contact algorithm to work, the zone where contact takes place must be specified. Parameter α is set to 0.14 for tetrahedral elements and 0.5 for hexahedral elements in order to provide a sufficiently smooth transition in passing from smaller to larger elements. Convergence is evaluated on radial stiffness corresponding to two different loads, specifically Q1=37 N and Q2=3706 N. According to (5), they correspond to two radial loads Fr equal to 0.1 kN and 10 kN. System stiffness ki corresponding to the applied load Qi is computed as(6)ki=Qiδavg,i,where δavg,i is the average displacement on the inner ring surface on the plane of maximum load, i.e., the plane sectioning the roller along the ideal line of contact.In addition, an improved mesh with load-dependent contact area and element dimension is proposed. This process is carried out based on the consideration that the critical area for these simulations is the contact zone, whose half-widtha decreases as the applied load Q reduces. Therefore, the involved elements should be sufficiently small to properly address the contact phenomenon for a wide load range. To tackle this problem, local mesh size l is changed at each load step based on the value of a computed by using the Hertz theory for that particular load. In this way, the contact phenomenon may be captured at each load, and grids become finer or coarser depending on the actual loading condition. Parameters w and α are also changed at each Q in order to maintain a low number of nodes. Results are compared with the initial meshes generated with constant w and α to assess advantages in terms of computational time. The improved meshing method will be employed to generate the grid for the reference bearing, as its mesh should be sufficiently light to attain a low computational time while taking into account the multiple contacts occurring in the bearing itself.It is worth noting that this work is focused on displacement evaluation and not on trying to obtain accurate pressure and stress distribution in the contact zones. To accurately estimate them, a different mesh optimization shall be performed. For this reason, the proposed convergence checks are only performed on the radial displacement, avoiding considerations regarding pressure and stresses.


### 3.2. Tetrahedral Elements: CTETRA4 vs. CTETRA10

Figure3(a) reports an example of the tetrahedral mesh applied to the bearing slice. Five meshes are generated, each denoted by different values of l and r. Table 2 reports their values along with the corresponding number of nodes for both 4-node and 10-node elements. Figures 4(a) and 4(b) show the estimated radial stiffness at Q1=37 N and Q2=3706 N, respectively. The corresponding computational time values for each simulation are depicted in Figures 4(c) and 4(d). At load Q2, the resulting stiffness is consistently lower for parabolic elements. Moreover, a constant trend is reached at a lower number of nodes. At load Q1, stiffness evaluated with linear elements is unstable for low number of nodes as the estimated value increases at first and then starts to descend. On the contrary, parabolic elements show a more stable trend which decreases as the number of nodes increases. From this analysis, CTETRA10 elements appear to be more efficient and reliable than their linear counterparts, in particular at higher loads.Table 2
Mesh dimensions and number of nodes for models meshed with CTETRA4 linear elements and CTETRA10 parabolic elements.
CTETRA4CTETRA10l (mm)r (mm)NodesDOFsl (mm)r (mm)NodesDOFs0.050.3579024240750.151.07112826348570.0250.17932406891210.10.71422506616770.020.143494281370500.050.357587671651340.010.0712133096057850.0250.1792204876306460.00750.05439387611262640.020.143341098980308Figure 4
Convergence for CTETRA elements: (a) stiffness at 37 N; (b) stiffness at 3706 N; (c) computational time at 37 N; (d) computational time at 3706 N.
(a)(b)(c)(d)The contact area half-width, namely,w, is modified based on the applied load according to the value estimated by (4). Table 3 shows theoretical half-width computed with (4), where subscript i refers to inner ring-roller contact, while subscript o refers to outer ring-roller contact. Contact area width w and local mesh size l are defined as proportional to computed a values on both sides of the roller. Parameter α is assumed to be equal to 0.06 at Q1 and 0.22 at Q2 to efficiently adapt the dimension of the elements of the rest of the model to local mesh size. Three possible combinations of w and l are proposed. They are reported in Table 4, along with the simulation results. Reported data show that the differences in computational time are negligible for load Q2, but they are relevant for load Q1, as the third combination allows to great reduction of the computational time in comparison with the other two combinations. Stiffness change among these approaches is relatively small, as the maximum percentage difference is 0.5% for Q1 and 0.2% for Q2. As a result, combination 3 is the approach chosen to generate the mesh for the roller bearing.Table 3
Contact area half-widtha estimated by (4).
Q1=37 NQ2=3706 Nai (mm)0.01270.127ao (mm)0.01590.159Table 4
Tested combinations ofw and l for the model meshed with CTETRA10 elements.
Combinationw (mm)l (mm)Q1=37 NQ2=3706 Nk (N/mm)t (s)k (N/mm)t (s)11.5a0.75a80710454131673821.25a1.25a81057447131446431.5a1.5a811502161317654Stiffness and computational time values obtained with the two improved meshes can be seen in Figure4. At applied load Q1=37 N, the maximum percentage difference with respect to the finest generated CTETRA10 grid is 0.1%, while the computational time reduces by 82.6%. At applied load Q1=3706 N, instead, stiffness deviates by 1.8% and time decreases by two orders of magnitude.


### 3.3. Hexahedral Elements: CHEXA8 vs. CHEXA20

The process carried out for tetrahedral elements is repeated for hexahedral elements. In addition to dimensionsl and r, a further parameter is needed, which is represented by the number of subdivisions along the shaft axis direction, nsubs. This value is assumed to be equal to 4, as shown in Figure 3(b). Such a value is chosen as a compromise between the number of nodes along the contact line and the local mesh size. It should be noted, in fact, that if obtaining a certain number of grid points is required, as Nsubs increases, l decreases accordingly. In addition, for two grids with the same value of parameter l, computational time will depend on Nsubs as it affects the total number of nodes and contacting elements. Thus, it is crucial to find a compromise between accuracy and computational time.As for tetrahedral elements, five meshes are generated for both linear and parabolic elements. Table5 shows the number of nodes for each one, as well as detailing element dimensions for the contact area and the rest of the model. Stiffness results are shown in Figures 5(a) and 5(b). At load Q2, both element types provide a stable trend even for a low number of nodes. At load Q1, parabolic and linear elements have opposite trends, since the former decrease as the number of nodes increases, while the latter behave in the opposite way. They eventually stabilize for a similar number of nodes. Parabolic mesh, however, allows obtaining stiffness values closer to stabilized values at lower number of nodes compared to linear elements. The computational time, in this case, is comparable for both element types. In fact, as shown in Figures 5(c) and 5(d), this parameter greatly increases for parabolic elements at higher number of nodes, but for coarser grids the difference is relatively low. Hence, CHEXA20 elements may be employed to generate the improved mesh.Table 5
Mesh dimensions and number of nodes for models meshed with CHEXA8 linear elements and CHEXA20 parabolic elements.
CHEXA8CHEXA20l (mm)r (mm)NodesDOFsl (mm)r (mm)NodesDOFs0.050.116405446030.10.20016783457970.0250.05619751708990.050.100615091711510.0150.031674654645260.0350.071221203419030.01250.0252423106732370.0250.0502338256572870.010.0237323510385410.020.0403640851025653Figure 5
Convergence for CHEXA elements: (a) stiffness at 37 N; (b) stiffness at 3706 N; (c) computational time at 37 N; (d) computational time at 3706 N.
(a)(b)(c)(d)Such a mesh, which will be later used to model the bearing, is generated by using the same method described inSection 3.2. Further load-dependent meshes are generated based on contact area dimensions calculated with (4) as reported in Table 3. Parameter α is assumed to be equal to 0.06 at Q1 and 0.5 at Q2. Table 6 reports tested combinations of parameters w and l along with the computed stiffness and the associated computational time values. The difference between these values and the ones reported in Table 4, which refers to tetrahedral elements, is worth noting. Hexahedral elements, in fact, allow considerable reduction of computational time. This is especially true at low load values, since the number of elements may be controlled more efficiently compared to tetrahedral elements.Table 6
Tested combinations ofw and l for the model meshed with CHEXA20 elements when nsubs=4.
Combinationw (mm)l (mm)Q1=37 NQ2=3706 Nk (N/mm)t (s)k (N/mm)t (s)13a0.5a81420661297191023aa8092031130258733a1.5a80604221302524Based on the results shown in Table6, the third combination is chosen to generate the improved mesh. Figure 5 shows evaluated stiffness obtained with the improved mesh. With respect to the finest CHEXA20 grid, the stiffness deviates by 0.5% for both loads. Computational time appears to be reduced by one to two orders of magnitude. These results show the benefits obtained by employing the proposed meshing method. By considering a local element size proportional to the contact area estimated through the Hertz theory, meshes showing a good compromise between accuracy and computational time may be obtained without the need to carry out dedicated convergence checks. The grids obtained in this manner, in fact, allow attaining results relatively close to converged ones but for a lower number of nodes.


## 3.1. Procedure for Mesh Size Performance Evaluation

Convergence analyses are carried out to determine a reasonable mesh size that allows obtaining a good estimation of the radial bearing stiffness while maintaining a relatively low computational time. An angular sector of the bearing is taken as reference in order to evaluate elements with progressively reduced dimensions. Geometry is cut to consider only half of the contact zone, whose length isLeff/2. Thus, only a quarter of the roller in contact with a portion of the two races is modeled, as the system is symmetrical with respect to two planes. The angle span of the reduced model, then, decreases to Δϕ=360/2Z=16.36∘. Figure 3 shows the dimensions of the adopted angular sector. The load is applied by providing a sinusoidal load distribution on the inner ring, which resembles the force exchanged with the shaft. Bodies are meshed with either tetrahedral or hexahedral elements, namely, CTETRA and CHEXA, in Nastran environment. For both elements, linear and parabolic formulations are available and identified by a number after their name, defining their number of nodes. As a consequence, element CTETRA4 indicates a linear tetrahedral element with four nodes, while CTETRA10 defines a parabolic element with ten nodes. Similarly, element CHEXA8 indicates a linear hexahedral element with eight nodes, while CHEXA20 defines a parabolic element with twenty nodes.Figure 3
Grids of the reduced geometry employed for convergence testing: (a) tetrahedral elements; (b) hexahedral elements and dimensions of the adopted angular sector.
(a)(b)Based on the necessity to determine the most convenient element, a procedure to determine element performance is established. The proposed method is applicable to all tested element types, with only few modifications needed to switch from tetrahedral to hexahedral elements. Convergence analysis is performed by decreasing the mesh size, which is quantified as the average length of element edges, so that the total number of nodes of the model progressively increases. Mesh dimension in the contact area,l, is smaller than the mesh size of the rest of the model, r. As the local dimension is reduced, r is kept proportional to l by means of a constant coefficient α so that l=αr. The size of the contact area is kept as a constant term for all load values. Its half-width is taken to be w=0.25 mm, which is a conservative value obtained through (4). This area must be geometrically defined because, for the contact algorithm to work, the zone where contact takes place must be specified. Parameter α is set to 0.14 for tetrahedral elements and 0.5 for hexahedral elements in order to provide a sufficiently smooth transition in passing from smaller to larger elements. Convergence is evaluated on radial stiffness corresponding to two different loads, specifically Q1=37 N and Q2=3706 N. According to (5), they correspond to two radial loads Fr equal to 0.1 kN and 10 kN. System stiffness ki corresponding to the applied load Qi is computed as(6)ki=Qiδavg,i,where δavg,i is the average displacement on the inner ring surface on the plane of maximum load, i.e., the plane sectioning the roller along the ideal line of contact.In addition, an improved mesh with load-dependent contact area and element dimension is proposed. This process is carried out based on the consideration that the critical area for these simulations is the contact zone, whose half-widtha decreases as the applied load Q reduces. Therefore, the involved elements should be sufficiently small to properly address the contact phenomenon for a wide load range. To tackle this problem, local mesh size l is changed at each load step based on the value of a computed by using the Hertz theory for that particular load. In this way, the contact phenomenon may be captured at each load, and grids become finer or coarser depending on the actual loading condition. Parameters w and α are also changed at each Q in order to maintain a low number of nodes. Results are compared with the initial meshes generated with constant w and α to assess advantages in terms of computational time. The improved meshing method will be employed to generate the grid for the reference bearing, as its mesh should be sufficiently light to attain a low computational time while taking into account the multiple contacts occurring in the bearing itself.It is worth noting that this work is focused on displacement evaluation and not on trying to obtain accurate pressure and stress distribution in the contact zones. To accurately estimate them, a different mesh optimization shall be performed. For this reason, the proposed convergence checks are only performed on the radial displacement, avoiding considerations regarding pressure and stresses.


## 3.2. Tetrahedral Elements: CTETRA4 vs. CTETRA10

Figure3(a) reports an example of the tetrahedral mesh applied to the bearing slice. Five meshes are generated, each denoted by different values of l and r. Table 2 reports their values along with the corresponding number of nodes for both 4-node and 10-node elements. Figures 4(a) and 4(b) show the estimated radial stiffness at Q1=37 N and Q2=3706 N, respectively. The corresponding computational time values for each simulation are depicted in Figures 4(c) and 4(d). At load Q2, the resulting stiffness is consistently lower for parabolic elements. Moreover, a constant trend is reached at a lower number of nodes. At load Q1, stiffness evaluated with linear elements is unstable for low number of nodes as the estimated value increases at first and then starts to descend. On the contrary, parabolic elements show a more stable trend which decreases as the number of nodes increases. From this analysis, CTETRA10 elements appear to be more efficient and reliable than their linear counterparts, in particular at higher loads.Table 2
Mesh dimensions and number of nodes for models meshed with CTETRA4 linear elements and CTETRA10 parabolic elements.
CTETRA4CTETRA10l (mm)r (mm)NodesDOFsl (mm)r (mm)NodesDOFs0.050.3579024240750.151.07112826348570.0250.17932406891210.10.71422506616770.020.143494281370500.050.357587671651340.010.0712133096057850.0250.1792204876306460.00750.05439387611262640.020.143341098980308Figure 4
Convergence for CTETRA elements: (a) stiffness at 37 N; (b) stiffness at 3706 N; (c) computational time at 37 N; (d) computational time at 3706 N.
(a)(b)(c)(d)The contact area half-width, namely,w, is modified based on the applied load according to the value estimated by (4). Table 3 shows theoretical half-width computed with (4), where subscript i refers to inner ring-roller contact, while subscript o refers to outer ring-roller contact. Contact area width w and local mesh size l are defined as proportional to computed a values on both sides of the roller. Parameter α is assumed to be equal to 0.06 at Q1 and 0.22 at Q2 to efficiently adapt the dimension of the elements of the rest of the model to local mesh size. Three possible combinations of w and l are proposed. They are reported in Table 4, along with the simulation results. Reported data show that the differences in computational time are negligible for load Q2, but they are relevant for load Q1, as the third combination allows to great reduction of the computational time in comparison with the other two combinations. Stiffness change among these approaches is relatively small, as the maximum percentage difference is 0.5% for Q1 and 0.2% for Q2. As a result, combination 3 is the approach chosen to generate the mesh for the roller bearing.Table 3
Contact area half-widtha estimated by (4).
Q1=37 NQ2=3706 Nai (mm)0.01270.127ao (mm)0.01590.159Table 4
Tested combinations ofw and l for the model meshed with CTETRA10 elements.
Combinationw (mm)l (mm)Q1=37 NQ2=3706 Nk (N/mm)t (s)k (N/mm)t (s)11.5a0.75a80710454131673821.25a1.25a81057447131446431.5a1.5a811502161317654Stiffness and computational time values obtained with the two improved meshes can be seen in Figure4. At applied load Q1=37 N, the maximum percentage difference with respect to the finest generated CTETRA10 grid is 0.1%, while the computational time reduces by 82.6%. At applied load Q1=3706 N, instead, stiffness deviates by 1.8% and time decreases by two orders of magnitude.


## 3.3. Hexahedral Elements: CHEXA8 vs. CHEXA20

The process carried out for tetrahedral elements is repeated for hexahedral elements. In addition to dimensionsl and r, a further parameter is needed, which is represented by the number of subdivisions along the shaft axis direction, nsubs. This value is assumed to be equal to 4, as shown in Figure 3(b). Such a value is chosen as a compromise between the number of nodes along the contact line and the local mesh size. It should be noted, in fact, that if obtaining a certain number of grid points is required, as Nsubs increases, l decreases accordingly. In addition, for two grids with the same value of parameter l, computational time will depend on Nsubs as it affects the total number of nodes and contacting elements. Thus, it is crucial to find a compromise between accuracy and computational time.As for tetrahedral elements, five meshes are generated for both linear and parabolic elements. Table5 shows the number of nodes for each one, as well as detailing element dimensions for the contact area and the rest of the model. Stiffness results are shown in Figures 5(a) and 5(b). At load Q2, both element types provide a stable trend even for a low number of nodes. At load Q1, parabolic and linear elements have opposite trends, since the former decrease as the number of nodes increases, while the latter behave in the opposite way. They eventually stabilize for a similar number of nodes. Parabolic mesh, however, allows obtaining stiffness values closer to stabilized values at lower number of nodes compared to linear elements. The computational time, in this case, is comparable for both element types. In fact, as shown in Figures 5(c) and 5(d), this parameter greatly increases for parabolic elements at higher number of nodes, but for coarser grids the difference is relatively low. Hence, CHEXA20 elements may be employed to generate the improved mesh.Table 5
Mesh dimensions and number of nodes for models meshed with CHEXA8 linear elements and CHEXA20 parabolic elements.
CHEXA8CHEXA20l (mm)r (mm)NodesDOFsl (mm)r (mm)NodesDOFs0.050.116405446030.10.20016783457970.0250.05619751708990.050.100615091711510.0150.031674654645260.0350.071221203419030.01250.0252423106732370.0250.0502338256572870.010.0237323510385410.020.0403640851025653Figure 5
Convergence for CHEXA elements: (a) stiffness at 37 N; (b) stiffness at 3706 N; (c) computational time at 37 N; (d) computational time at 3706 N.
(a)(b)(c)(d)Such a mesh, which will be later used to model the bearing, is generated by using the same method described inSection 3.2. Further load-dependent meshes are generated based on contact area dimensions calculated with (4) as reported in Table 3. Parameter α is assumed to be equal to 0.06 at Q1 and 0.5 at Q2. Table 6 reports tested combinations of parameters w and l along with the computed stiffness and the associated computational time values. The difference between these values and the ones reported in Table 4, which refers to tetrahedral elements, is worth noting. Hexahedral elements, in fact, allow considerable reduction of computational time. This is especially true at low load values, since the number of elements may be controlled more efficiently compared to tetrahedral elements.Table 6
Tested combinations ofw and l for the model meshed with CHEXA20 elements when nsubs=4.
Combinationw (mm)l (mm)Q1=37 NQ2=3706 Nk (N/mm)t (s)k (N/mm)t (s)13a0.5a81420661297191023aa8092031130258733a1.5a80604221302524Based on the results shown in Table6, the third combination is chosen to generate the improved mesh. Figure 5 shows evaluated stiffness obtained with the improved mesh. With respect to the finest CHEXA20 grid, the stiffness deviates by 0.5% for both loads. Computational time appears to be reduced by one to two orders of magnitude. These results show the benefits obtained by employing the proposed meshing method. By considering a local element size proportional to the contact area estimated through the Hertz theory, meshes showing a good compromise between accuracy and computational time may be obtained without the need to carry out dedicated convergence checks. The grids obtained in this manner, in fact, allow attaining results relatively close to converged ones but for a lower number of nodes.


## 4. Evaluation of the Bearing Radial Stiffness

This section is dedicated to the description of the procedure adopted to model the reference bearing and to estimate its radial stiffness. Procedures are exploited to generate the bearing mesh along with the meshing method introduced inSection 3. Various approaches to estimate the radial bearing stiffness are proposed, and results obtained through numerical and analytical procedures are reported. Influence of rollers’ position is assessed by varying their location with respect to the direction of maximum load. Additional simulations in which the cage is inserted into the model are run in order to determine if it is reasonable to neglect this component in REBs modeling. Lastly, the model is further improved to take into account the diametral clearance.


### 4.1. Model Reduction Strategy

The size of the computational domain is reduced by introducing a series of modeling strategies. The first one consists in taking advantage of the available symmetry planes to reduce the bearing geometry. The first plane is commonly used to decrease the dimension of the computational domain in bearing problems, and it is the one passing through the shaft axis on the plane of maximum load. It is available because of the geometry of the mechanical component. The second plane is the one normal to the shaft axis by halving the bearing width. It may be exploited only by assuming no misalignment effects and pure radial load.Unloaded rollers may be removed to reduce the number of contacting bodies and consequently the number of grid points. As a matter of fact, the number of rollers carrying the load from one ring to the other depends on the extent of the load zone, namely,ϕl, defined as [8](7)ϕl=2cos−1Pd2δr,where Pd is the diametral clearance and δr is the radial displacement. When Pd=0, the angular extent of the load zone is 180° as long as δr≠0. This means that only half of the rolling elements are effectively contributing to transfer the load. Therefore, unloaded rollers may be removed to reduce the number of nodes and decrease the quantity of contact areas, leading to a substantial reduction in computational time. For the reference bearing in Figure 1, the number of modeled roller reduces from 11 to 3. It is worth noting that, when Pd>0, ϕl reduces, while it increases when Pd<0. In these cases, attention must be paid to how many rollers are effectively loaded to avoid erroneously removing load-carrying rollers.Concerning the boundary conditions, the reference bearing is considered to be mounted on an infinitely rigid frame. Thus, it is clamped on the outer surface of the outer ring, meaning that all DOFs for that node are restrained.By considering the fixed boundary condition, the angular extent of the load zone, and the absence of some rollers, a portion of the outer and inner ring may also be removed. This is done by assuming that they do not make a significant contribution to the radial displacement in the plane of maximum load, i.e., the monitored displacement for subsequent stiffness assessment. As a consequence, the size of the computational domain is effectively reduced to 1/8th of the original geometry. In addition, since a 0.1 mm axial clearance between rollers and flanges is considered, the lateral flange is removed from the remaining geometry as there is no contact between the rollers and the flange itself. All these practical assumptions allow great reduction of the number of nodes needed to generate the mesh of the mechanical component.Finally, attention is devoted to alternative methods to apply the radial load. As seen inSection 1, there are various methods to apply the radial load, e.g., on the center of the shaft, on rings, or on a central node connected with rigid links to the inner ring. In order to decrease the number of grid nodes, it suggested that the shaft contribution is modeled by means of an analytical sinusoidal distribution which replicates the load applied by a radially loaded shaft on the inner ring. Such load also depends on the diametral clearance Pd of the bearing. In absence of clearance, the load insists on a 180° angular sector as depicted in Figure 6. This permits the removal of one contacting body while still taking into account its effect. The efficiency of this assumption is verified in the next section, after the definition of the bearing mesh.Figure 6
Sinusoidal load distribution acting on a180° angular sector of the inner surface of the inner ring.It is worth noting that these model reduction strategies are employed in this work for a cylindrical roller bearing but they may also be exploited to generate the grid of other bearing types. However, although assumptions concerning load application, unloaded rollers removal, element type, and mesh size are still applicable for all rolling-element bearings subjected to a radial load, the symmetry plane which halves the bearing width might not be always utilized, e.g., for tapered roller bearing and self-aligning double-row ball bearings. Consequently, the symmetry conditions should be always checked before reducing the model geometry as they are not the same for every bearing.


### 4.2. Mesh Generation

Figure7 shows examples of tetrahedral and hexahedral meshes exploiting all the model reduction strategies defined in Section 4.1. Due to the nonlinear relationship between load and displacement, a certain number of load values must be applied to obtain the stiffness curve. In this case, stiffness is evaluated for six different load values, in particular 100 N, 500 N, 1000 N, 2500 N, 5000 N, and 10000 N. Meshes are generated according to methods reported in Sections 3.2 and 3.3. Consequently, a mesh is defined for each load and characterized by different values of w, l, and α. Table 7 reports the values referring to both types of mesh for each loading condition. It is worth noting that the contact area semiwidth a is based solely on the value acting on the maximum loaded roller, i.e., computed through (5). The other two rollers carry a smaller load, but contact area dimensions are kept the same.Figure 7
Reduced bearing geometry meshed with (a) tetrahedral elements and (b) hexahedral elements.
(a)(b)Table 7
Contact area and average elements dimensions for CTETRA10 and CHEXA20 elements for different loading conditions.
Q (N)a (mm)CTETRA10CHEXA20w (mm)l (mm)αw (mm)l (mm)α1000.0130.0190.0190.060.0390.0190.065000.0280.0430.0430.080.0840.0430.1210000.0400.0600.0600.120.1200.0600.1725000.0640.0950.0950.140.1920.0950.2650000.0900.1350.1350.170.2700.1350.36100000.1270.1910.1910.220.3810.1910.5Hexahedral grids are chosen to verify the assumption made on the application of the radial load. The efficiency of the method is demonstrated by evaluating the load distribution on rollers. Results obtained by employing the sinusoidal distribution are compared to distributions computed with the other two methods. The first one is the analytical formulation reported in [8], for which the load distribution for Pd=0 mm may be computed as(8)Fϕ=Qmaxcosϕ10/9,where Qmax is calculated through (5). The second method consists in modeling the shaft and applying the total load in its central node. Size of the elements in the shaft mesh is the same as in the inner ring. Since the diametral clearance is not considered, shaft and inner ring touch on the entirety of the available contact surface. Two radial loads, Fr1=0.1 kN and Fr2=10 kN, are applied. For the two FE methods, load distribution is calculated for each roller as the sum of the contact forces in the corresponding contact areas. Results are shown in Figure 8, where the roller at angle 0° is the one whose center is aligned with the direction of maximum load. The sinusoidal load allows obtaining a load distribution which deviates by up to 3% from analytical results. The loaded shaft, on the other hand, shows a similar difference for the rollers at 0° and 32.7° while it deviates by 8.3% at maximum for the roller at 62.4°. Besides, the addition of the shaft causes the computational time to increase by 28.7% for Fr1 and by 56.8% for Fr2. These results show that it is possible to remove the shaft and substitute it with a sinusoidal load in order to eliminate one contacting body and lead to a consistent reduction in computational time.Figure 8
Load distribution on rollers evaluated for two different radial load values: (a)Fr1=0.1 kN; (b) Fr1=10 kN.
(a)(b)


### 4.3. Radial Stiffness Estimation: Analytical and Numerical Results

Radial stiffness may be evaluated in two different ways. The work in [9] suggests computing it as(9)k1=Qxr,where xr is the radial displacement due to a radial load Q; Guo and Parker [20], on the contrary, recommend using partial derivatives(10)k2=∂Q∂xr.Both methods require the evaluation of the displacement-load relationship for different loads to obtain the stiffness over a certain load range. Parameterk1 requires one evaluation for each load, since only the displacement at each load Q is needed. Parameter k2, on the other hand, requires at least two displacement evaluations for each load to compute the partial derivatives about load Q. These displacements correspond to applied forces F1=Q+δQ and F2=Q−δQ, where parameter δQ is a small disturbance which is arbitrarily assumed to be 0.01% of the load Q for each loading condition. Therefore, stiffness k2Q is computed as(11)k2Q=F1−F2xrF1−xrF2.This process may be repeated for any number of loads. This means that, ifN stiffness values are required, N simulations must be run for parameter k1 and 2N simulations have to be performed for parameter k2 to obtain stiffness estimations over the whole load range. Radial displacement xr is evaluated as the average displacement of the nodes on the plane of maximum load on the inner surface of the inner ring.Results are compared against three different stiffness-load laws obtained through methods found in the literature. They are chosen since each one exploits a different approach to determine bearing stiffness. Within this framework, Harris’ formula [8] is based on an analytical approach, Gargiulo’s [6] equation is based on empirical observations, and Petersen et al.’s [38] method is derived from a lumped parameter model of the bearing. With the exception of Gargiulo’s equation, which directly provides the stiffness k1 due to load Q, these methods are adapted to allow computing both k1 and k2.Figure9 compares stiffness values estimated with the analytical methods and the finite-element method: Figure 9(a) reports estimated k1 values, whereas Figure 9(b) presents k2 values. For both charts, stiffness estimated with two different element types provides similar results, as the maximum difference between hexahedral and tetrahedral elements is 3.7% for parameter k1 and 3.3% for parameter k2. Furthermore, they assume values consistently between Petersen’s and Harris’s estimations, which represent the upper and lower boundary, respectively. Parameter k2 is always higher than parameter k1 at all loads, for all considered methods. It is interesting to note the variety of stiffness results provided by analytical and numerical methods. FEM estimates deviate by up to 9.5% compared to Harris’ formula and up to 13.2% with respect to Petersen’s lumped parameter model. Moreover, Gargiulo’s equation provides a consistently lower estimation which deviates by up to 16.6%.Figure 9
Bearing radial stiffness estimated with analytical and numerical methods, for two different evaluation methods: (a)k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)Computational time due to tetrahedral and hexahedral element is shown in Figure10(a) for tetrahedral elements and in Figure 10(b) for hexahedral elements. They refer to a single simulation for each load. By employing the proposed procedure, the computational time at low loads greatly increases as tetrahedral elements are utilized. The simulation, in fact, is considerably demanding as the dimension of the elements in the contact area causes the number of nodes to rapidly increase. Hexahedral elements, on the contrary, are more suitable for simulations at low loads since the number of nodes may be controlled more efficiently. At higher loads, element type choice still has an influence on computational time but the order of magnitude is comparable. According to these results, hexahedral elements are well-fitted to estimate radial stiffness at low loads, while for higher loads the choice has less impact on the computational time. It is worth noting that the employed approach requires a longer setup time compared to traditional models as a different grid has to be generated for each load. However, this additional time is compensated by the decreased computational time needed to solve the problem and the lower time required to perform convergence checks, as the element dimension is analytically determined. Moreover, the setup process may be further accelerated by automatizing the procedure in FE software. Then, despite the increased number of generated grids, the computational time may be lessened compared to traditional approaches.Figure 10
Computational time for FE simulations: (a) tetrahedral elements; (b) hexahedral elements.
(a)(b)The validity of the assumptions made is verified by comparing the full bearing model without any geometrical modifications to the results obtained with the reduced model. The full bearing is meshed by employing hexahedral elements, with the same load-dependent dimensions utilized for the reduced model. Figure11 shows parameters k1 and k2 for both models. It is worth noting that there is some difference between the estimated stiffness, especially at lower loads where the deviation is up to 3.5%. The computational time, however, increases by two orders of magnitude at each load, leading to a considerable rise of the total time needed to run all the simulations. These results suggest that reducing the computational domain by adopting the proposed model reductions allows consistently decreasing the computational effort, even if there may be a slight change in the estimated stiffness. The proposed method, in fact, is developed in order to obtain an acceptable compromise between accuracy and simulation time.Figure 11
Comparison between the radial bearing stiffness evaluated for the reduced and the complete model, for two different evaluation methods: (a)k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)In addition to these results, it is worth noting that for this kind of bearings the position of the rollers with respect to the direction of maximum load has a slight influence on the estimated bearing stiffness values. This may be verified by considering two different roller positions, namely, position 1 and position 2, and comparing the resulting radial stiffness obtained with the proposed method. In position 1, the plane of maximum loads passes through the middle of one roller, as described inSection 4.2. In position 2, the plane is centered between two consecutive rollers. Examples of models meshed with parabolic hexahedral elements are depicted in Figures 12(a) and 12(b). Figure 13 shows the stiffness-load relationship for both positions. For parameter k1, the estimated stiffness in position 2 is consistently lower compared to position 1, up to 3.8% at Q=10 kN. Parameter k2 shows a less stable trend at position 1, with a maximum percentage difference of 6.7% at Q=1 kN. Whether different load directions should be considered or neglected in REBs analysis, then, depends on the required degree of accuracy.Figure 12
Models meshed with CHEXA20 elements with different applied load directions: (a) through the middle of one roller; (b) between two consecutive rollers.
(a)(b)Figure 13
Bearing radial stiffness evaluated for two different roller positions: (a)k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)


### 4.4. Inclusion of the Cage

The cage is often included in REBs in order to maintain a constant distance between rolling elements. Its presence, however, is commonly neglected in FE modeling, as detailed inSection 1. To investigate the influence of this component on radial bearing stiffness, its effect is evaluated by performing simulations in which the cage is effectively meshed and inserted into the model, as shown in Figure 14. The aim of the inquiry is to determine whether it is acceptable not to include this component in a FE model. The cage is modeled so that each slot contacts a roller on both sides in the circumferential direction, whereas an axial clearance equal to 0.1 mm prevents contact between pockets and roller ends. Tetrahedral elements are used since they allow modeling the complex cage geometry. Element dimensions within the cage are kept proportional to l through a coefficient, 2α, which varies with the load as reported in Table 7. Results for both k1 and k2 are depicted in Figure 15. The maximum percentage difference between simulations with and without cage is up to 2.1% for k1 and 3.7% for k2. It is worth noting, however, that 3.7% variation for k2 is obtained for only one load value, while for the other loading conditions the difference is below 1.4%. The drawback introduced by cage modeling, however, is the increase in computational time, as the average time needed to run the simulations increases by 2.4 times compared to models neglecting this component. Therefore, these results suggest that to avoid modeling the cage is a reasonable expedient to reduce computational effort without compromising the accuracy of the results.Figure 14
Alternative bearing model considering cage presence. The bodies are meshed with CTETRA10 elements.Figure 15
Bearing radial stiffness, numerically evaluated, taking into account the presence of the cage: (a)k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)Figure 16
Bearing radial stiffness computed by considering a diametral clearancePd equal to 0.04 mm: (a) k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)


### 4.5. Clearance Effect

Diametral clearance, namely,Pd, is an additional factor that affects radial bearing stiffness. Equation (7), in fact, shows that for values of Pd greater than 0 mm the extent of the load zone reduces, leading to a redistribution of the load among the rollers. In particular, the maximum load on roller increases according to the following equation [8]:(12)Qmax=Knδr−Pd2n,where Kn is the load-deflection factor for a roller in contact with two raceways, and the exponent n, associated with this kind of contact, is equal to 10/9. Factor Kn may be computed as(13)Kn=121/Kl1/nn,where Kl is the load-deflection factor for a steel roller in contact with a single steel raceway [8]:(14)Kl=8.06·104leff8/9.For the reference roller bearing, the total load-deflection factor isKn=172549 N/mm. Furthermore, the number of loaded rollers may decrease, depending on the applied load and the angular span Δϕr between two consecutive rollers. Consequently, the radial displacement δr increases, thus diminishing the radial bearing stiffness. The value of Qmax is needed to compute the contact area at each load in order to define the element size for the associated mesh, but (12) cannot be directly solved since the radial displacement is not initially known as it is equal to(15)δr=Pd2+QZKnJrPd,δr1/n,where Z is the number of rolling elements and Jr is a radial integral depending on Pd and δr. Since Jr depends on the radial displacement, (15) may be solved by means of an iterative procedure. Once δr is known for each radial load, Qmax is calculated by employing (12). Finally, contact area semiwidth is computed through (4) while the loaded angular sector is calculated via (7).Clearance effect is assessed by means of numerical and analytical approaches. A diametral clearancePd=0.04 mm is considered and inserted between the components. Load-dependent meshes to be employed for FE simulation are generated according to the methods shown in Section 4. Parabolic hexahedral elements are chosen and mesh dimensions l and r are selected according to the size of the contact area, as detailed in Section 4.2. Contact area half-width and angular extent of the load zone for different loads are reported in Table 8. It is worth noting that, as the applied load reduces, ϕl decreases since the radial displacement diminishes accordingly. Hence, the number of loaded rollers may reduce, depending on the ratio between ϕl and the angular distance Δϕr, which is equal to 32.72° for the reference bearing. For a bearing in position 1 as in Figure 12(a), three rollers are loaded at 1000 N, but only two of them are loaded between 500 and 5000 N and a single one is loaded for a radial force equal to 100 N. At each load, unloaded rollers may be removed from the FE model in order to lighten the mesh. The analytical sinusoidal load is not applied on a 180° angular sector as in the simulations performed for Pd=0 mm, but the angular span is changed at each load according to ϕl values reported in Table 8.Table 8
Clearance effect on contact area dimension and load zone extent. The number of loaded rollers refers to a bearing in position 1 as in Figure12(a), for which the maximum number of loaded rollers for Pd=0 mm is three.
Q (N)ai (mm)ao (mm)αϕl°Loaded rollers1000.0200.0250.0638.015000.0410.0510.1267.6210000.0530.0670.1782.4225000.0760.0950.26104.2250000.1010.1270.36121.62100000.1370.1720.50138.23Parametersk1 and k2 are evaluated, and FE results are compared with Harris’ equation and Petersen et al.’s lumped parameter model, as they both allow considering clearance effect. Figure 16 shows the estimated parameters for various loads. As expected, stiffness values evaluated with either k1 or k2 are lower than values obtained in absence of clearance, as the ones depicted in Figure 9. For parameter k1, FEM results deviate from Harris equation by 6.1% at most, while the difference with Petersen et al.’s model is higher, up to 10.7% at load 2500 N. For parameter k2, FEM results deviate from Harris equation by 7.5% at most; hence, the difference is similar to parameter k1. On the other hand, there is a relatively high difference with Petersen et al.’s model, up to 31.1% at load 1000 N. Influence of Pd on k2, in fact, is very low on Petersen et al.’s model. These results suggest that the FE method may be employed to quantify bearing radial stiffness when Pd>0 mm. In this case, however, analytical computations should be carried out to determine Qmax and ϕl, leading to more calculations and a longer setup time compared to the case in which Pd=0 mm.


## 4.1. Model Reduction Strategy

The size of the computational domain is reduced by introducing a series of modeling strategies. The first one consists in taking advantage of the available symmetry planes to reduce the bearing geometry. The first plane is commonly used to decrease the dimension of the computational domain in bearing problems, and it is the one passing through the shaft axis on the plane of maximum load. It is available because of the geometry of the mechanical component. The second plane is the one normal to the shaft axis by halving the bearing width. It may be exploited only by assuming no misalignment effects and pure radial load.Unloaded rollers may be removed to reduce the number of contacting bodies and consequently the number of grid points. As a matter of fact, the number of rollers carrying the load from one ring to the other depends on the extent of the load zone, namely,ϕl, defined as [8](7)ϕl=2cos−1Pd2δr,where Pd is the diametral clearance and δr is the radial displacement. When Pd=0, the angular extent of the load zone is 180° as long as δr≠0. This means that only half of the rolling elements are effectively contributing to transfer the load. Therefore, unloaded rollers may be removed to reduce the number of nodes and decrease the quantity of contact areas, leading to a substantial reduction in computational time. For the reference bearing in Figure 1, the number of modeled roller reduces from 11 to 3. It is worth noting that, when Pd>0, ϕl reduces, while it increases when Pd<0. In these cases, attention must be paid to how many rollers are effectively loaded to avoid erroneously removing load-carrying rollers.Concerning the boundary conditions, the reference bearing is considered to be mounted on an infinitely rigid frame. Thus, it is clamped on the outer surface of the outer ring, meaning that all DOFs for that node are restrained.By considering the fixed boundary condition, the angular extent of the load zone, and the absence of some rollers, a portion of the outer and inner ring may also be removed. This is done by assuming that they do not make a significant contribution to the radial displacement in the plane of maximum load, i.e., the monitored displacement for subsequent stiffness assessment. As a consequence, the size of the computational domain is effectively reduced to 1/8th of the original geometry. In addition, since a 0.1 mm axial clearance between rollers and flanges is considered, the lateral flange is removed from the remaining geometry as there is no contact between the rollers and the flange itself. All these practical assumptions allow great reduction of the number of nodes needed to generate the mesh of the mechanical component.Finally, attention is devoted to alternative methods to apply the radial load. As seen inSection 1, there are various methods to apply the radial load, e.g., on the center of the shaft, on rings, or on a central node connected with rigid links to the inner ring. In order to decrease the number of grid nodes, it suggested that the shaft contribution is modeled by means of an analytical sinusoidal distribution which replicates the load applied by a radially loaded shaft on the inner ring. Such load also depends on the diametral clearance Pd of the bearing. In absence of clearance, the load insists on a 180° angular sector as depicted in Figure 6. This permits the removal of one contacting body while still taking into account its effect. The efficiency of this assumption is verified in the next section, after the definition of the bearing mesh.Figure 6
Sinusoidal load distribution acting on a180° angular sector of the inner surface of the inner ring.It is worth noting that these model reduction strategies are employed in this work for a cylindrical roller bearing but they may also be exploited to generate the grid of other bearing types. However, although assumptions concerning load application, unloaded rollers removal, element type, and mesh size are still applicable for all rolling-element bearings subjected to a radial load, the symmetry plane which halves the bearing width might not be always utilized, e.g., for tapered roller bearing and self-aligning double-row ball bearings. Consequently, the symmetry conditions should be always checked before reducing the model geometry as they are not the same for every bearing.


## 4.2. Mesh Generation

Figure7 shows examples of tetrahedral and hexahedral meshes exploiting all the model reduction strategies defined in Section 4.1. Due to the nonlinear relationship between load and displacement, a certain number of load values must be applied to obtain the stiffness curve. In this case, stiffness is evaluated for six different load values, in particular 100 N, 500 N, 1000 N, 2500 N, 5000 N, and 10000 N. Meshes are generated according to methods reported in Sections 3.2 and 3.3. Consequently, a mesh is defined for each load and characterized by different values of w, l, and α. Table 7 reports the values referring to both types of mesh for each loading condition. It is worth noting that the contact area semiwidth a is based solely on the value acting on the maximum loaded roller, i.e., computed through (5). The other two rollers carry a smaller load, but contact area dimensions are kept the same.Figure 7
Reduced bearing geometry meshed with (a) tetrahedral elements and (b) hexahedral elements.
(a)(b)Table 7
Contact area and average elements dimensions for CTETRA10 and CHEXA20 elements for different loading conditions.
Q (N)a (mm)CTETRA10CHEXA20w (mm)l (mm)αw (mm)l (mm)α1000.0130.0190.0190.060.0390.0190.065000.0280.0430.0430.080.0840.0430.1210000.0400.0600.0600.120.1200.0600.1725000.0640.0950.0950.140.1920.0950.2650000.0900.1350.1350.170.2700.1350.36100000.1270.1910.1910.220.3810.1910.5Hexahedral grids are chosen to verify the assumption made on the application of the radial load. The efficiency of the method is demonstrated by evaluating the load distribution on rollers. Results obtained by employing the sinusoidal distribution are compared to distributions computed with the other two methods. The first one is the analytical formulation reported in [8], for which the load distribution for Pd=0 mm may be computed as(8)Fϕ=Qmaxcosϕ10/9,where Qmax is calculated through (5). The second method consists in modeling the shaft and applying the total load in its central node. Size of the elements in the shaft mesh is the same as in the inner ring. Since the diametral clearance is not considered, shaft and inner ring touch on the entirety of the available contact surface. Two radial loads, Fr1=0.1 kN and Fr2=10 kN, are applied. For the two FE methods, load distribution is calculated for each roller as the sum of the contact forces in the corresponding contact areas. Results are shown in Figure 8, where the roller at angle 0° is the one whose center is aligned with the direction of maximum load. The sinusoidal load allows obtaining a load distribution which deviates by up to 3% from analytical results. The loaded shaft, on the other hand, shows a similar difference for the rollers at 0° and 32.7° while it deviates by 8.3% at maximum for the roller at 62.4°. Besides, the addition of the shaft causes the computational time to increase by 28.7% for Fr1 and by 56.8% for Fr2. These results show that it is possible to remove the shaft and substitute it with a sinusoidal load in order to eliminate one contacting body and lead to a consistent reduction in computational time.Figure 8
Load distribution on rollers evaluated for two different radial load values: (a)Fr1=0.1 kN; (b) Fr1=10 kN.
(a)(b)


## 4.3. Radial Stiffness Estimation: Analytical and Numerical Results

Radial stiffness may be evaluated in two different ways. The work in [9] suggests computing it as(9)k1=Qxr,where xr is the radial displacement due to a radial load Q; Guo and Parker [20], on the contrary, recommend using partial derivatives(10)k2=∂Q∂xr.Both methods require the evaluation of the displacement-load relationship for different loads to obtain the stiffness over a certain load range. Parameterk1 requires one evaluation for each load, since only the displacement at each load Q is needed. Parameter k2, on the other hand, requires at least two displacement evaluations for each load to compute the partial derivatives about load Q. These displacements correspond to applied forces F1=Q+δQ and F2=Q−δQ, where parameter δQ is a small disturbance which is arbitrarily assumed to be 0.01% of the load Q for each loading condition. Therefore, stiffness k2Q is computed as(11)k2Q=F1−F2xrF1−xrF2.This process may be repeated for any number of loads. This means that, ifN stiffness values are required, N simulations must be run for parameter k1 and 2N simulations have to be performed for parameter k2 to obtain stiffness estimations over the whole load range. Radial displacement xr is evaluated as the average displacement of the nodes on the plane of maximum load on the inner surface of the inner ring.Results are compared against three different stiffness-load laws obtained through methods found in the literature. They are chosen since each one exploits a different approach to determine bearing stiffness. Within this framework, Harris’ formula [8] is based on an analytical approach, Gargiulo’s [6] equation is based on empirical observations, and Petersen et al.’s [38] method is derived from a lumped parameter model of the bearing. With the exception of Gargiulo’s equation, which directly provides the stiffness k1 due to load Q, these methods are adapted to allow computing both k1 and k2.Figure9 compares stiffness values estimated with the analytical methods and the finite-element method: Figure 9(a) reports estimated k1 values, whereas Figure 9(b) presents k2 values. For both charts, stiffness estimated with two different element types provides similar results, as the maximum difference between hexahedral and tetrahedral elements is 3.7% for parameter k1 and 3.3% for parameter k2. Furthermore, they assume values consistently between Petersen’s and Harris’s estimations, which represent the upper and lower boundary, respectively. Parameter k2 is always higher than parameter k1 at all loads, for all considered methods. It is interesting to note the variety of stiffness results provided by analytical and numerical methods. FEM estimates deviate by up to 9.5% compared to Harris’ formula and up to 13.2% with respect to Petersen’s lumped parameter model. Moreover, Gargiulo’s equation provides a consistently lower estimation which deviates by up to 16.6%.Figure 9
Bearing radial stiffness estimated with analytical and numerical methods, for two different evaluation methods: (a)k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)Computational time due to tetrahedral and hexahedral element is shown in Figure10(a) for tetrahedral elements and in Figure 10(b) for hexahedral elements. They refer to a single simulation for each load. By employing the proposed procedure, the computational time at low loads greatly increases as tetrahedral elements are utilized. The simulation, in fact, is considerably demanding as the dimension of the elements in the contact area causes the number of nodes to rapidly increase. Hexahedral elements, on the contrary, are more suitable for simulations at low loads since the number of nodes may be controlled more efficiently. At higher loads, element type choice still has an influence on computational time but the order of magnitude is comparable. According to these results, hexahedral elements are well-fitted to estimate radial stiffness at low loads, while for higher loads the choice has less impact on the computational time. It is worth noting that the employed approach requires a longer setup time compared to traditional models as a different grid has to be generated for each load. However, this additional time is compensated by the decreased computational time needed to solve the problem and the lower time required to perform convergence checks, as the element dimension is analytically determined. Moreover, the setup process may be further accelerated by automatizing the procedure in FE software. Then, despite the increased number of generated grids, the computational time may be lessened compared to traditional approaches.Figure 10
Computational time for FE simulations: (a) tetrahedral elements; (b) hexahedral elements.
(a)(b)The validity of the assumptions made is verified by comparing the full bearing model without any geometrical modifications to the results obtained with the reduced model. The full bearing is meshed by employing hexahedral elements, with the same load-dependent dimensions utilized for the reduced model. Figure11 shows parameters k1 and k2 for both models. It is worth noting that there is some difference between the estimated stiffness, especially at lower loads where the deviation is up to 3.5%. The computational time, however, increases by two orders of magnitude at each load, leading to a considerable rise of the total time needed to run all the simulations. These results suggest that reducing the computational domain by adopting the proposed model reductions allows consistently decreasing the computational effort, even if there may be a slight change in the estimated stiffness. The proposed method, in fact, is developed in order to obtain an acceptable compromise between accuracy and simulation time.Figure 11
Comparison between the radial bearing stiffness evaluated for the reduced and the complete model, for two different evaluation methods: (a)k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)In addition to these results, it is worth noting that for this kind of bearings the position of the rollers with respect to the direction of maximum load has a slight influence on the estimated bearing stiffness values. This may be verified by considering two different roller positions, namely, position 1 and position 2, and comparing the resulting radial stiffness obtained with the proposed method. In position 1, the plane of maximum loads passes through the middle of one roller, as described inSection 4.2. In position 2, the plane is centered between two consecutive rollers. Examples of models meshed with parabolic hexahedral elements are depicted in Figures 12(a) and 12(b). Figure 13 shows the stiffness-load relationship for both positions. For parameter k1, the estimated stiffness in position 2 is consistently lower compared to position 1, up to 3.8% at Q=10 kN. Parameter k2 shows a less stable trend at position 1, with a maximum percentage difference of 6.7% at Q=1 kN. Whether different load directions should be considered or neglected in REBs analysis, then, depends on the required degree of accuracy.Figure 12
Models meshed with CHEXA20 elements with different applied load directions: (a) through the middle of one roller; (b) between two consecutive rollers.
(a)(b)Figure 13
Bearing radial stiffness evaluated for two different roller positions: (a)k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)


## 4.4. Inclusion of the Cage

The cage is often included in REBs in order to maintain a constant distance between rolling elements. Its presence, however, is commonly neglected in FE modeling, as detailed inSection 1. To investigate the influence of this component on radial bearing stiffness, its effect is evaluated by performing simulations in which the cage is effectively meshed and inserted into the model, as shown in Figure 14. The aim of the inquiry is to determine whether it is acceptable not to include this component in a FE model. The cage is modeled so that each slot contacts a roller on both sides in the circumferential direction, whereas an axial clearance equal to 0.1 mm prevents contact between pockets and roller ends. Tetrahedral elements are used since they allow modeling the complex cage geometry. Element dimensions within the cage are kept proportional to l through a coefficient, 2α, which varies with the load as reported in Table 7. Results for both k1 and k2 are depicted in Figure 15. The maximum percentage difference between simulations with and without cage is up to 2.1% for k1 and 3.7% for k2. It is worth noting, however, that 3.7% variation for k2 is obtained for only one load value, while for the other loading conditions the difference is below 1.4%. The drawback introduced by cage modeling, however, is the increase in computational time, as the average time needed to run the simulations increases by 2.4 times compared to models neglecting this component. Therefore, these results suggest that to avoid modeling the cage is a reasonable expedient to reduce computational effort without compromising the accuracy of the results.Figure 14
Alternative bearing model considering cage presence. The bodies are meshed with CTETRA10 elements.Figure 15
Bearing radial stiffness, numerically evaluated, taking into account the presence of the cage: (a)k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)Figure 16
Bearing radial stiffness computed by considering a diametral clearancePd equal to 0.04 mm: (a) k1=Q/xr; (b) k2=∂Q/∂xr.
(a)(b)


## 4.5. Clearance Effect

Diametral clearance, namely,Pd, is an additional factor that affects radial bearing stiffness. Equation (7), in fact, shows that for values of Pd greater than 0 mm the extent of the load zone reduces, leading to a redistribution of the load among the rollers. In particular, the maximum load on roller increases according to the following equation [8]:(12)Qmax=Knδr−Pd2n,where Kn is the load-deflection factor for a roller in contact with two raceways, and the exponent n, associated with this kind of contact, is equal to 10/9. Factor Kn may be computed as(13)Kn=121/Kl1/nn,where Kl is the load-deflection factor for a steel roller in contact with a single steel raceway [8]:(14)Kl=8.06·104leff8/9.For the reference roller bearing, the total load-deflection factor isKn=172549 N/mm. Furthermore, the number of loaded rollers may decrease, depending on the applied load and the angular span Δϕr between two consecutive rollers. Consequently, the radial displacement δr increases, thus diminishing the radial bearing stiffness. The value of Qmax is needed to compute the contact area at each load in order to define the element size for the associated mesh, but (12) cannot be directly solved since the radial displacement is not initially known as it is equal to(15)δr=Pd2+QZKnJrPd,δr1/n,where Z is the number of rolling elements and Jr is a radial integral depending on Pd and δr. Since Jr depends on the radial displacement, (15) may be solved by means of an iterative procedure. Once δr is known for each radial load, Qmax is calculated by employing (12). Finally, contact area semiwidth is computed through (4) while the loaded angular sector is calculated via (7).Clearance effect is assessed by means of numerical and analytical approaches. A diametral clearancePd=0.04 mm is considered and inserted between the components. Load-dependent meshes to be employed for FE simulation are generated according to the methods shown in Section 4. Parabolic hexahedral elements are chosen and mesh dimensions l and r are selected according to the size of the contact area, as detailed in Section 4.2. Contact area half-width and angular extent of the load zone for different loads are reported in Table 8. It is worth noting that, as the applied load reduces, ϕl decreases since the radial displacement diminishes accordingly. Hence, the number of loaded rollers may reduce, depending on the ratio between ϕl and the angular distance Δϕr, which is equal to 32.72° for the reference bearing. For a bearing in position 1 as in Figure 12(a), three rollers are loaded at 1000 N, but only two of them are loaded between 500 and 5000 N and a single one is loaded for a radial force equal to 100 N. At each load, unloaded rollers may be removed from the FE model in order to lighten the mesh. The analytical sinusoidal load is not applied on a 180° angular sector as in the simulations performed for Pd=0 mm, but the angular span is changed at each load according to ϕl values reported in Table 8.Table 8
Clearance effect on contact area dimension and load zone extent. The number of loaded rollers refers to a bearing in position 1 as in Figure12(a), for which the maximum number of loaded rollers for Pd=0 mm is three.
Q (N)ai (mm)ao (mm)αϕl°Loaded rollers1000.0200.0250.0638.015000.0410.0510.1267.6210000.0530.0670.1782.4225000.0760.0950.26104.2250000.1010.1270.36121.62100000.1370.1720.50138.23Parametersk1 and k2 are evaluated, and FE results are compared with Harris’ equation and Petersen et al.’s lumped parameter model, as they both allow considering clearance effect. Figure 16 shows the estimated parameters for various loads. As expected, stiffness values evaluated with either k1 or k2 are lower than values obtained in absence of clearance, as the ones depicted in Figure 9. For parameter k1, FEM results deviate from Harris equation by 6.1% at most, while the difference with Petersen et al.’s model is higher, up to 10.7% at load 2500 N. For parameter k2, FEM results deviate from Harris equation by 7.5% at most; hence, the difference is similar to parameter k1. On the other hand, there is a relatively high difference with Petersen et al.’s model, up to 31.1% at load 1000 N. Influence of Pd on k2, in fact, is very low on Petersen et al.’s model. These results suggest that the FE method may be employed to quantify bearing radial stiffness when Pd>0 mm. In this case, however, analytical computations should be carried out to determine Qmax and ϕl, leading to more calculations and a longer setup time compared to the case in which Pd=0 mm.


## 5. Conclusions

In this paper, procedures and hypotheses commonly adopted in the FE modeling of rolling-element bearings were analyzed and discussed with the purpose of efficiently reducing the size of the computational domain and concurrently providing a robust approach to evaluate their radial stiffness. The proposed methodology consists in generating load-dependent grids, so that the mesh size in the contact area is properly tuned according to the dimensions of the analytically estimated contact zone between rolling elements and raceways. This strategy allows obtaining proper element dimensions to capture the contact phenomenon for a wide range of loads, while limiting the number of grid points and curbing the computational time. Bearing meshes were generated by employing either parabolic tetrahedral or parabolic hexahedral elements. The size of the computational domain was further reduced by taking advantage of the available symmetry planes, by removing unloaded rollers and replacing the shaft with an equivalent sinusoidal load distribution. Radial stiffness values were then computed by exploiting two different formulations, namely,k1 and k2, based on the chosen stiffness expression. The results showed that both considered element types provided similar radial stiffness estimates, although tetrahedral elements exhibited considerably higher computational time at low loads. Good agreement was also found by comparing these results with different analytical approaches retrieved from the literature, thus denoting the quality of the proposed methodology.The described bearing model was also employed to assess the variation in stiffness due to different roller positions with respect to the direction of maximum load. The analysis showed a deviation on the computed stiffness by up to6.7%, which indicates that the inclusion of this effect within the analysis might become mandatory depending on the degree of accuracy needed. It must be noted, however, that the total computational time increases if multiple roller positions are simulated. The FE model was further utilized to evaluate the stiffness variation caused by adding the cage as a meshed body. Inserting this component led to a slight change in stiffness, whereas the computational time consistently increased. Therefore, modeling the cage is not suggested if the goal is to attain a low computational time. Finally, clearance effect was assessed by physically inserting it between the components and by changing the sinusoidal load distribution in order to apply the load in a reduced angular sector. Unloaded rollers were also removed when the load angular span was sufficiently small. Numerical results were in good agreement with analytical formulations.



---

*Source: 9955398-2021-08-30.xml*