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

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.


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. is 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]. e 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][8][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.
us, researchers strived to simplify the problem using a variety of approaches. e 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. e 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. is 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. e efficiency of this method encouraged other researchers [15][16][17][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. e first one is similar to the one employed by Daidié et al. [14], as rollers are replaced by springs. e 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. ey 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. ey 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. ey 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. e 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][14][15][16][17][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][11][12]19], on rings [14,18,21], or on a central node connected with rigid elements to the inner ring [15][16][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. e 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. is 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. e 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 selfaligning 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. e following section outlines general aspects regarding the simulations described throughout this paper. Section 3 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.

Reference Geometry.
A roller bearing, model SKF 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. e 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/ m 3 , and Poisson's ratio ] � 0.28. ese properties are considered as constant values throughout the analyses.

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 nodeto-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]. ese 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].
is 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]. e 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. e 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. ese 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]. e 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 ballraceway contact. As a result, only data obtained for the reference roller bearing are compared in this work.

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. e 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. e 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 L eff in contact with a curved surface and loaded with a force F is a rectangle of length L eff and half-width a. Let the contacting bodies be labeled with 1 and 2. R eq is their equivalent contact radius: where R 1 and R 2 are the radii of curvature of the contacting surfaces. ey may take positive or negative values if the surfaces are convex or concave, respectively. An equivalent Young's modulus E * , depending on material properties of both components, may also be defined as Parameter E * is needed to compute the indentation depth d due to an applied load F: which may be eventually utilized to determine the contact area half-width a as 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 rollerouter 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 Table 1, an operating load range within F r � 0 N and F r � 10 kN is assumed based on its basic static load rating C 0 � 10.2 kN. F r is the radial load applied on the shaft and transmitted to the inner ring. e corresponding force on the most loaded roller may be computed, in absence of diametral clearance, by using the formula provided by Harris [8]: leading to a value Q max � 3706 N when F r � 10 kN. Figure 2 shows contact areas half-width for inner and outer ring contact with the most loaded roller.

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 is L eff /2. us, 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. e angle span of the reduced model, then, decreases to Δϕ � 360/(2Z) � 16.36 ∘ . Figure 3 shows the dimensions of the adopted angular sector. e 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. Based on the necessity to determine the most convenient element, a procedure to determine element performance is established. e 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. e 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). is 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 Q 1 � 37 N and Q 2 � 3706 N. According to (5), they correspond to two radial loads F r equal to 0.1 kN and 10 kN. System stiffness k i corresponding to the applied load Q i is computed as 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.
is process is carried out based on the consideration that the critical area for these simulations is the contact zone, whose half-width a decreases as the applied load Q reduces.
erefore, 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. e 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. Inner ring contact Outer ring contact

Mathematical Problems in Engineering
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. Figure 3(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 Q 1 � 37 N and Q 2 � 3706 N, respectively. e corresponding computational time values for each simulation are depicted in Figures 4(c) and 4(d). At load Q 2 , the resulting stiffness is consistently lower for parabolic elements. Moreover, a constant trend is reached at a lower number of nodes. At load Q 1 , 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. e 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 Q 1 and 0.22 at Q 2 to efficiently adapt the dimension of the elements of the rest of the model to local mesh size. ree possible combinations of w and l are proposed. ey are reported in Table 4, along with the simulation results. Reported data show that the differences in computational time are negligible for load Q 2 , but they are relevant for load Q 1 , 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 Q 1 and 0.2% for Q 2 . As a result, combination 3 is the approach chosen to generate the mesh for the roller bearing.

Tetrahedral Elements: CTETRA4 vs. CTETRA10.
Stiffness and computational time values obtained with the two improved meshes can be seen in Figure 4. At applied load Q 1 � 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 Q 1 � 3706 N, instead, stiffness deviates by 1.8% and time decreases by two orders of magnitude.

Hexahedral Elements: CHEXA8 vs. CHEXA20.
e process carried out for tetrahedral elements is repeated for hexahedral elements. In addition to dimensions l and r, a further parameter is needed, which is represented by the number of subdivisions along the shaft axis direction, n subs . is 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 N subs increases, l decreases accordingly. In addition, for two grids with the same value of parameter l, computational time will depend on N subs as it affects the total number of nodes and contacting elements.
us, 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. Table 5 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 Q 2 , both element types provide a stable trend even for a low number of nodes. At load Q 1 , 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. ey 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. e 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.
Such a mesh, which will be later used to model the bearing, is generated by using the same method described in Section 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 Q 1 and 0.5 at Q 2 . Table 6 reports tested combinations of parameters w and l along with the computed stiffness and the associated computational time values. e 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. is is especially true at low load values, since the number of elements may be controlled more efficiently compared to tetrahedral elements.
Based on the results shown in Table 6, 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. ese 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. e grids obtained in this manner, in fact, allow attaining results relatively close to converged ones but for a lower number of nodes.

Evaluation of the Bearing Radial Stiffness
is 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 in Section 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.

Model Reduction Strategy.
e size of the computational domain is reduced by introducing a series of modeling strategies. e first one consists in taking advantage of the available symmetry planes to reduce the bearing geometry. e 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. e 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] ϕ l � 2 cos −1 P d 2δ r , where P d is the diametral clearance and δ r is the radial displacement. When P d � 0, the angular extent of the load zone is 180°as long as δ r ≠ 0. is means that only half of the      rolling elements are effectively contributing to transfer the load. erefore, 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 P d > 0, ϕ l reduces, while it increases when P d < 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. us, 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. is 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 in Section 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 P d of the bearing. In absence of clearance, the load insists on a 180°angular sector as depicted in Figure 6. is permits the removal of one contacting body while still taking into account its effect. e efficiency of this assumption is verified in the next section, after the definition of the bearing mesh.
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. Figure 7 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). e other two rollers carry a smaller load, but contact area dimensions are kept the same.

Mesh Generation.
Hexahedral grids are chosen to verify the assumption made on the application of the radial load. e 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. e first one is the analytical formulation reported in [8], for which the load distribution for P d � 0 mm may be computed as where Q max is calculated through (5). e 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, F r1 � 0.1 kN and F r2 � 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. e sinusoidal load allows obtaining a load distribution which deviates by up to 3% from analytical results. e 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 F r1 and by 56.8% for F r2 . ese 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.

Radial Stiffness Estimation: Analytical and Numerical
Results. Radial stiffness may be evaluated in two different ways. e work in [9] suggests computing it as where x r is the radial displacement due to a radial load Q; Guo and Parker [20], on the contrary, recommend using partial derivatives Both methods require the evaluation of the displacement-load relationship for different loads to obtain the stiffness over a certain load range. Parameter k 1 requires one evaluation for each load, since only the displacement at each load Q is needed. Parameter k 2 , on the other hand, requires at least two displacement evaluations for each load to compute the partial derivatives about load Q. ese displacements correspond to applied forces F 1 � Q + δQ and F 2 � 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. erefore, stiffness k 2 (Q) is computed as Q max Q(ϕ)   is process may be repeated for any number of loads. is means that, if N stiffness values are required, N simulations must be run for parameter k 1 and 2N simulations have to be performed for parameter k 2 to obtain stiffness estimations over the whole load range. Radial displacement x r 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 stiffnessload laws obtained through methods found in the literature. ey 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 k 1 due to load Q, these methods are adapted to allow computing both k 1 and k 2 . Figure 9 compares stiffness values estimated with the analytical methods and the finite-element method: Figure 9(a) reports estimated k 1 values, whereas Figure 9(b) presents k 2 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 k 1 and 3.3% for parameter k 2 . Furthermore, they assume values consistently between Petersen's and Harris's estimations, which represent the upper and lower boundary, respectively. Parameter k 2 is always higher than parameter k 1 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%.
Computational time due to tetrahedral and hexahedral element is shown in Figure 10(a) for tetrahedral elements and in Figure 10(b) for hexahedral elements. ey 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. e 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. en, despite the increased number of generated grids, the computational time may be lessened compared to traditional approaches. e 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. e full bearing is meshed by employing hexahedral elements, with the same load-dependent dimensions utilized for the reduced model. Figure 11 shows parameters k 1 and k 2 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%. e 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. ese 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. e proposed method, in fact, is developed in order to obtain an acceptable compromise between accuracy and simulation time.
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. is 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 in Section 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 k 1 , the estimated stiffness in position 2 is consistently lower compared to position 1, up to 3.8% at Q � 10 kN. Parameter k 2 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.

Inclusion of the Cage.
e 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 in Section 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. e aim of the inquiry is to determine whether it is acceptable not to include this component in a FE model. e 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 k 1 and k 2 are depicted in Figure 15. e maximum percentage difference between simulations with and without cage is up to 2.1% for k 1 and 3.7% for k 2 . It is worth noting, however, that 3.7% variation for k 2 is obtained for only one load value, while for the other loading conditions the difference is below 1.4%. e 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. erefore, these results suggest that to avoid modeling the cage is a reasonable expedient to reduce computational effort without compromising the accuracy of the results.

Clearance
Effect. Diametral clearance, namely, P d , is an additional factor that affects radial bearing stiffness. Equation (7), in fact, shows that for values of P d 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]: where K n 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 K n may be computed as where K l is the load-deflection factor for a steel roller in contact with a single steel raceway [8]: K l � 8.06 · 10 4 l 8/9 eff .
For the reference roller bearing, the total load-deflection factor is K n � 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. e value of Q max 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 where Z is the number of rolling elements and J r is a radial integral depending on P d and δ r . Since J r depends on the radial displacement, (15) may be solved by means of an iterative procedure. Once δ r is known for each radial load, Q max 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 clearance P d � 0.04 mm is considered and inserted between the components. Loaddependent 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 halfwidth 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. e analytical sinusoidal load is not applied on a 180°a ngular sector as in the simulations performed for P d � 0 mm, but the angular span is changed at each load according to ϕ l values reported in Table 8. Parameters k 1 and k 2 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 k 1 or k 2 are lower than values obtained in absence of clearance, as the ones depicted in Figure 9. For parameter k 1 Figure 16: Bearing radial stiffness computed by considering a diametral clearance P d equal to 0.04 mm: (a) k 1 � Q/x r ; (b) k 2 � zQ/zx r . 16 Mathematical Problems in Engineering computations should be carried out to determine Q max and ϕ l , leading to more calculations and a longer setup time compared to the case in which P d � 0 mm.

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. e 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. is 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. e 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, k 1 and k 2 , based on the chosen stiffness expression. e 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. e 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. e analysis showed a deviation on the computed stiffness by up to 6.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. e 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. erefore, 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.

Data Availability
Data are available on request to the corresponding author. Table 8: Clearance effect on contact area dimension and load zone extent. e number of loaded rollers refers to a bearing in position 1 as in Figure 12(a), for which the maximum number of loaded rollers for P d � 0 mm is three.