Local grid refinement in multigrid method for point contact problems of polycrystalline anisotropic material under dry and lubricated conditions

Predicting rolling bearing fatigue life requires knowledge of the three-dimensional (3D) stress fields in the roller and raceway near the lubricated contact. Owing to the increasingly severe operating conditions, the effect of localized features such as surface roughness, subsurface inclusions, and even the crystallographic structure of the material becomes important. Achieving such detail requires (locally) extremely dense gridding in simulations, which in 3D is a major challenge. Multigrid techniques have been demonstrated to be capable of solving such problems. In this study, multigrid techniques are shown to further increase the efficiency of the solution by exploiting local grid refinement while maintaining the simplicity of a uniform discretization. This is achieved by employing increasingly finer grids only locally, where the highest resolution is required. Results are presented for dry contact and elastohydrodynamically lubricated contact cases, circular as well as elliptic, with varying crystallographic structure, and with surface roughness. The results show that the developed algorithm is very well suited for detailed analysis, with also excellent prospects for computational diagnostics involving actual material crystallographic structure from electron backscatter diffraction measurements.


Introduction
Rolling bearing steels are usually exposed to millions of Hertzian-type stress cycles during their service. As a result, the bearing raceways and rolling elements can undergo rolling contact fatigue (RCF). Depending on the operating conditions, such as load and material microstructure, there are different manifestations of RCF, which are summarized in a recent review paper [1]. Under well-lubricated conditions, RCF usually appears in the form of material spalling, driven by material heterogeneity. It has been recognized that the steel chemical composition, heat treatment, and cleanliness significantly affect the bearing life [2]. Therefore, understanding the effect of the steel microstructure on bearing performance is of crucial importance to meet the constant demand for reduced carbon dioxide emissions, higher energy efficiency, reduced friction, and cost across industries.
While often approximated as homogeneous for easier stress calculations, it is important to point out that in general the microstructure of bearing steels is very heterogeneous. Different aspects contribute to the overall heterogeneity. First, the chemical composition and heat treatment of the steel affect its texture as manifested through different martensite grain sizes (from laths to packets), grain orientation, amount of retained austenite, and the size, shape, and chemistry of residual and tempered carbides [3]. In addition, the metallurgical processes used in steel www.Springer.com/journal/40544 | Friction Nomenclature a Half width of Hertzian contact in x direction (m) b Half width of Hertzian contact in y direction (m) c 11 , c 12 , c 44 Elastic constants for cubic material (Pa) C i, j, k, l The elastic stiffness matrix (Pa), i, j, k, l are the indexes of the fouth-rank tensor d The

II
High order interpolation operator production influence the overall steel cleanliness (e.g., the presence of nonmetallic inclusions, material segregations, etc.). All of these aspects affect the overall subsurface stress state, which can potentially lead to faster or slower RCF. Therefore, a better understanding of the effect of the steel microstructure is of crucial importance for predicting and improving bearing performance. This means that the development of fast and accurate numerical tools for stress calculations that can handle the types of heterogeneity mentioned above is of high practical importance. The type of contact in a bearing (point versus line) depends on the geometry of the rolling elements. For roller bearings, the contact geometry can be simplified to a finite or an infinite line contact, whereas for ball bearings, the contact geometry is typically a point or elliptical contact. Both line and circular contact may be treated as special cases of an elliptical contact. Line contact problems have been studied intensively over the last few decades, and closedform analytical solutions are available [4]. However, analytical solutions are usually limited to smooth contact between isotropic homogeneous materials, and therefore, they cannot easily be applied to more realistic problems. Hence, different numerical techniques have been developed to address various contact problems. The numerical implementation for line contact problems is computationally inexpensive, and various interesting problems have been addressed [5,6]. Warhadpande and Sadeghi [7] investigated the effect of surface defects on RCF in heavily loaded lubricated contacts. The results in Ref. [7] showed that the time required for crack nucleation is negligible compared to the time required for crack propagation during spalling, which is in line with experimental observations. In the work of Bomidi and Sadeghi [8], both stress and accumulated plastic strain-based damage evolutions were considered in order to predict fatigue failure initiation and propagation in heavily loaded contacts. The model was demonstrated to be a phenomenological tool for reliable analyses of RCF life, scatter, and subsurface spalling. Moghaddam et al. [9] investigated the butterfly wing formation around non-metallic inclusions, and studied the effects of inclusion size, stiffness, and depth on crack formation. The results in Ref. [9] showed that the crack initiation locations and final spall shapes were similar to the observations of failed bearings.
It is well known that the isotropic elastic properties of bearing steels emerge on a macroscale when assuming a random distribution of the orientation of anisotropic crystals (martensite blocks). Locally, however, the different material regions are anisotropic with individual orientations, leading to a complex stress distribution that requires a detailed microscopic analysis using the concept of anisotropic elasticity. Hence, a full understanding of experimental data and accurate fatigue life predictions requires more detailed stress analysis using the concept of anisotropic elasticity at the scale of single grains. The analytic solution for an anisotropic homogenous material in plane strain can be derived from using Green's function for a half space and the Stroh formalism [10,11]. However, because it is rather difficult to handle these functions in the most general case, and because they are limited to homogenous materials, the numerical approach is considered to be more efficient. Paulson et al. [12] studied the effect of polycrystalline anisotropy on RCF using the finite element method (FEM), and found that the variation in crystal elasticity affects the fatigue life, which may help to explain the life scatter observed in experiments. Later, this work was extended to consider the influence of the lubricating film and crystal elasticity [13], which makes the simulation results much closer to experimental results than when employing isotropic material.
For a point contact, the stress field analysis requires much more computational effort because of the extension of the dimension from two-dimensional (2D) to three dimensional (3D). Wang et al. [14][15][16] applied discrete convolution and fast Fourier transform (DC-FFT) [17] in a numerical method to speed up the computation of the surface deformation and subsurface stress fields for layered materials in a point elastohydrodynamic lubricated (EHL) contact.
Zhang et al. [18] used an equivalent inclusion method [19] to investigate the effect of inhomogeneous materials on the lubrication performance. Essentially, in these cases, a point contact problem is solved as all material aspects are captured in the influence functions that describe the surface deformation. The related stresses are computed a posteriori using an FFT-based approach. Their results showed the influence of the position, size, and properties of the inclusions on the resulting stress field, film thickness, and pressure distribution. Golmohammadi and Sadeghi [20] developed a 3D finite element model to investigate the influence of refurbishing on RCF, and conducted a parametric study to evaluate the effect of material removal depth, load cycles, and applied load on the fatigue life. In Ref. [20], some simplifications were used to reduce the time cost, such as assuming a Hertzian pressure distribution and performing the computation in half of the 3D domain. Boffy et al. [21][22][23] developed a multigrid solution for the 3D stress field analysis of heterogeneous materials (granular material and fiber-reinforced material). Later, Zhang and Venner applied the multigrid method for the stress field analysis of polycrystalline anisotropic material under dry contact [24] and lubricated conditions [25], and investigated the influence of material anisotropy on RCF [26]. Some recent studies have shown that the multigrid method is also suitable for the stress field analysis of composite materials [27]. Finally, parallel computing can be used to further enhance the efficiency of multigrid methods [28,29].
To summarize, the numerical solution of the contact stress field has been achieved using semi-analytical methods [30], such as the integral transform combined with multilevel multi integration (MLMI) [31], fast Fourier transform (FFT) [14,15,18], FEM [7,9,12,13,20], and the finite difference method combined with multigrid methods [21][22][23][24][25][26][27][28]. The former two methods are suitable for some relatively simple cases in which the influence coefficient (the function relating contact pressure and stresses) is available. The integral transform combined with an FFT is suitable for the displacement and stress field analysis of layered, inhomogeneous, viscoelastic materials [32], and its advantage is the fast calculation of the top boundary displacement. In previous studies, multigrid methods have been shown to have high efficiency in solving the stress field for heterogeneous [21][22][23], composite [27], and polycrystalline anisotropic [24][25][26] materials. The efficiency of the multigrid method can be further increased by incorporating local grid refinement (LGR), in which high-resolution solutions are computed only locally where necessary. It should be noted that in multigrid solvers based on Ref. [31], local grid refinement is possible; e.g., Ref. [33]. However, in Ref. [33], the standard discretization of the semiinfinite domain elastic deformation integrals yields influence coefficients that are mesh size dependent, and hence requires special measures to be able to make the MLMI algorithm work efficiently [34].
To facilitate a faster and more efficient analysis of the stress conditions in heterogeneous anisotropic materials, in this study the multigrid method developed in Refs. [24][25][26] was extended to incorporate the LGR technique [35,36]. The new implementation was tested on both homogeneous isotropic and polycrystalline anisotropic materials under 3D dry and lubricated contact conditions. The accuracy of the result was verified by comparing with solutions obtained without LGR. The remainder of this paper is organized as follows. In Section 2, we introduce the theoretical foundation for the multigrid elastic solver, and explain the implementation of grid refinement. In Section 3, we verify the current implementation for isotropic and anisotropic materials under both dry and lubricated contact, resembling typical contact cases in bearings. Section 4 discusses the computational efficiency of the new implementation. The new implementation is then tested on two different contact phenomena in Section 5. Finally, in Section 6, we revisit the main findings and further research directions.

Local grid refinement (LGR): Background and theory
The sides of the 3D body are traction-free: where σ ij is the stress tensor, and j n is the surface normal. Finally, zero-displacement boundary conditions are set on the bottom surface (see Fig. 1). The deformation behavior of the elastic body is described by the Navier-Cauchy equations in three dimensions [37]:  (4) in which C global is the stiffness matrix in the global coordinate system X global , Y global , Z global , and is defined by where  2 ,  ,  are three Euler rotation angles for each grain, R x , R y , and R z are three rotation matrices, and where c 11 , c 12 , and c 44 are three material constants. To solve these partial differential equations numerically, they are discretized using a discretization scheme such as the finite-difference method (FDM) or FEM. In this study, FDM was used. The system of equations was discretized with second-order accuracy on a uniform grid with the same mesh size in each spatial direction. The discrete equations are given in Refs. [24,25,27]. The calculation domain was further divided into small cells (volumes) by Voronoi tessellation: A grain (Voronoi cell) consists of all grid points of the calculation domain closer to that seed (crystal nucleus) than to any other. Each Voronoi cell was assigned three Euler rotation angles (  2 ,  ,  ) to mimic the grain topology of the bearing material [12]. For isotropic material, the elastic stiffness matrix is the same after rotation, which means that the global stiffness matrix of each grid point is the same. For anisotropic material, the elastic stiffness matrix becomes full after rotation, which introduces additional terms in the stress-strain relation compared to isotropic material. This means that the global stiffness matrix varies between grains because of the difference between the Euler rotation angles. The size of the grains varies depending on the Voronoi tessellation. The shape of the grains can be reflected by the variation of C global of each grid point. The discretization of the above equations leads to a system of algebraic equations for the displacements at each grid point. A simple way to solve this system is through the use of iterative methods such as Gauss-Seidel relaxation. Such methods generally reduce high-frequency (with respect to the grid) errors efficiently in the solution. However, low-frequency (smooth) error components are reduced at a much slower rate, and dominate the asymptotic convergence behavior [38].
Multigrid methods accelerate the convergence of such iterative methods by solving for these lowfrequency errors on coarser grids. This improves efficiency in two ways: 1) Low-frequency errors appear more oscillatory on coarse grids, allowing them to be reduced effectively by relaxations. 2) Performing relaxations on coarser grids requires much less computational work owing to the reduced number of grid points. Assuming standard coarsening with the coarse grid having twice the mesh size of the fine grid, the number of points is reduced by a factor of eight in 3D problems. The key to the successful development of a multigrid algorithm is a relaxation process that provides good error smoothing. For a system of equations, as considered here, this also depends on the coupling between the equations. For the Navier-Cauchy equations, if the Poisson ratio is away from the incompressible case of 0.5, a simple sequential relaxation per unknown works sufficiently well to yield grid-independent convergence of the multigrid cycle. Regarding the boundary condition, a collective relaxation is required for dry contact, and a line relaxation in the boundary plane is required for the elastohydrodynamically lubricated contact. This is discussed in detail in Refs. [21][22][23][24][25][26].
Suppose a solution is sought on the finest (target) grid with mesh size h f . A multigrid (MG) cycle consists of the following steps.
1) Starting with an initial approximation on the finest grid,  1 relaxations are performed on the target grid, yielding an improved approximation from which 2) The problem is transferred to the next-coarser grid using a suitable choice of intergrid transfers and coarse-grid operators.
3) The coarse grid problem is solved recursively by performing a number of multigrid cycles. The coarse-grid solution is then interpolated to provide a correction to the next-finer grid. 4) Finally,  2 relaxations are performed on the fine grid to reduce high-frequency errors that may have been introduced by the interpolation.
When appropriately designed, the multigrid cycles yield a convergence rate independent of the target grid mesh size, in an amount of computational work that is equivalent to just a few iterations carried out on the target grid. In addition, coarser grids can be used to generate a good first approximation, so that on the target grid, one already starts with a small error. This is the essence of a full multigrid (FMG) algorithm, which consists of performing MG cycles on each of the grids starting from the coarsest grid and interpolating the solution to supply a good approximation to the initial solutions of the nextfiner grid. Figure 2 illustrates the algorithm for three grid levels when one multigrid W cycle is used to solve the coarse-grid problem. In principle, the FMG algorithm can solve the finest-grid problem in O(N) operations, where N is the number of unknowns in the algebraic system of equations [38]. For a more detailed description of the equations being solved by the FMG algorithm, see the Appendix. B. Additionally, interested readers are referred to Refs. [36,38] for more information on multigrid methods.
In our previous work, a FMG algorithm based on a finite-difference discretization was developed for dry contact [24] and lubricated point contact problems [25]. These problems require a large computational domain to justify the assumption of zerostress or zero-displacement boundary conditions far from the contact. Additionally, small mesh sizes must be used to capture the intricacies of the stress field in polycrystalline anisotropic materials. This leads to many unknowns, especially as it is a 3D problem with three equations and unknowns (displacements). The total work of FMG algorithm for 3D using W cycles is ideally twice the total work involved in the relaxations on the finest grid in a cycle. As this number is usually small O (3)(4)(5), the total work of solving the problem is the equivalent of approximately (6-10) relaxations on the fine grid to solve the problem to the level of the discretization error, which is extremely effective. Further optimization can be obtained by parallelizing the algorithm using multiple processors [28,29], or, as is the topic in this study, reduce the dominant contribution to the amount of work, i.e., the work done on the finest grid. Figure 3 shows the solution time spent on each grid for the studied problems in the case of five grid levels, taken from solving an actual problem. The graph shows that the total work is indeed approximately twice the work actually done on the finest grid. To reduce the amount of work further, one should reduce the amount of work done on the finest grid, which can be achieved by exploiting this finest grid only where absolutely necessary for solution accuracy.
Discretization errors in finite-difference approximations depend on the chosen mesh size and the size of the gradients in the solution. In the contact problems considered, large gradients in the displacements are expected only in the proximity of the contact. This suggests that the computational cost can be reduced by using fine meshes only in a small region surrounding the contact and a coarser  | https://mc03.manuscriptcentral.com/friction mesh elsewhere. Such LGR can greatly reduce the computational time and memory usage while retaining accurate results for the stress field, and contact pressure close to the contact. Figure 4 shows a schematic demonstration of a 2D grid without and with LGR. The red (Fig. 4(f)) and blue (Fig. 4(g)) dashed lines indicate regions covered by the finer mesh of a higher grid level. These patches of increasingly finer grids are used only at locations where they are required for the accuracy of the solution. Consequently, the total number of grid points that must be relaxed on the fine grids is significantly reduced compared to the full-domain grids, as shown in Figs. 4(c) and 4(d). Significant savings in computational time and memory cost can therefore be obtained if the fine-grid patches are small relative to the size of the domain. This is the case for the contact problems studied in Refs. [24,25], and the implementation of LGR in the developed solvers can be expected to make them more efficient and hence even more suitable for engineering applications.
The multigrid formulation in the full approximation  scheme allows for the straightforward implementation of LGR. The advantage here is that, combined with the FMG algorithm, a solution is generated on an effectively non-uniform grid that is realized by increasingly fine uniform local grids. A more detailed explanation is given in the Appendix C. Further background on the use of multigrid with LGR can be found in Refs. [35,36,38]. In this study, the LGR technique is incorporated into the previously developed multigrid methods [24,25,39], and it is shown that this leads to a more efficient algorithm, while the accuracy of the solutions is almost unaffected. The results are verified against solutions obtained using grids covering the full domain for heterogeneous grain-structured material with varying crystallographic orientation and for a few other cases characteristic of (lubricated) contact mechanics simulations. It is emphasized that the purpose of this paper is to explain and demonstrate the natural way in which local grid refinement can be implemented in a multigrid solution algorithm, and to show LGR accuracy, capability, and robustness. Parameter studies regarding various aspects of the problem have already been presented in Refs. [24][25][26]. In this study, we consider only the case of a single refinement region with a rectangles shape. However, it is straightforward to have multiple refined regions (patches) at each grid level, which in turn can be refined in multiple regions or complex refinement shapes. Multiple refined regions with complex shapes can be used for applications such as cluster or inclusion cases [40]. For more details, please refer to the original papers of Brandt et al. [35,41], to Chapter 3 of Ref. [38], and to Ref. [42].

Verification
In this section, the results of MG with LGR are compared with those without LGR and other numerical methods under dry contact and lubricated conditions. The contact and operating parameters are listed in Table 1. The isotropic and cubic anisotropic material parameters are listed in Table 2. The isotropic stiffness tensor is chosen such that its Young's modulus matches the experimental value of 206 GPa measured for 100Cr6 (SAE-AISI 52100) martensitic bearing steels.
The anisotropic cubic elastic constants corresponding to the body-centered cubic (bcc) Fe were taken from Ref. [43]. Note that it is not uncommon to use some of the averaging methods, such as Voigt or Reuss, to compute the approximate isotropic elastic constants (see Appendix A) [44]. For the sake of completeness, the Voigt and Reuss averages for Young's modulus are estimated to be  v 227 GPa E and  R 194 GPa E , respectively. Because 100Cr6 is a low-alloyed steel, the two values are the upper and lower bounds of the experimental value, as expected thermodynamically [45,46]. The calculation domains and number of grid points on each level are listed in Table 3. In this study, LGR starts from grid level 3 to ensure a good approximation from the coarse grids, i.e., the MG approach is carried out on the whole domain for grid levels below 3, and from level 3 calculations are performed on locally defined regions. The borders of LGR are listed in Table 3. The effect of mesh size is given in Appendix D. Verification was performed for both isotropic and polycrystalline anisotropic materials  of the bearing body material. The dimensionless results are as follows.

Isotropic material results
In this section, the numerical results with LGR are compared with those without LGR and other numerical methods for dry and lubricated contact problems assuming isotropic material behavior of the bearing body. The results are summarized in LGR for both contact problems. Second, we compare the von Mises stress (VMS) as a function of X and Z, for Y = 0 (compare Figs. 5(a) and 5(e) with Figs. 5(b) and 5(f), respectively) and of Y and Z, for X = 0 (compare Figs. 5(c) and 5(g) with Figs. 5(d) and 5(h), respectively). We note that the reduced mesh density in the part of the domain outside of the LGR region locally led to step-like solutions displaying minor quantitative differences with the solution of a fine grid over the entire computational domain. However, this part of the domain is typically located far away from the zone of the maximum Hertzian contact, and thus is typically less important. In addition, RCF is mostly controlled by the stress in the zone, which is very well described by LGR, and it does not obviously affect the solution results at the surface. Furthermore, for the sake of quantitative clarity, Figure 6 presents the plots of surface pressure, surface displacement, film thickness, and VMS. Under dry contact, the maximum VMS locates along the central line (X = Y = 0). Under lubrication, the value of the maximum VMS is larger than that of the dry contact, and its position moves to the right owing to the second pressure spike. The results of using LGR agree well with those obtained without using LGR and also with the results of the numerical methods in Refs. [16,38]. The maximum percent errors of the surface pressure, surface displacement/film thickness, and VMS with and without LGR on level 5 of the central XZ (Y = 0)   www.Springer.com/journal/40544 | Friction and 6(h)). For more detailed comparisons of dry contact and lubrication conditions on pressure and subsurface stress field distribution, interested readers are referred to Ref. [26].

Polycrystalline anisotropic material results
To further test the accuracy and reliability of the LGR technique, solutions were also obtained considering a polycrystalline anisotropic material. In this example, we include the most important length scale of the high-C martensite 100Cr6. To approximate the topology of the bearing material, a granular approach was adopted. Millions of grains with various local stiffness tensors were used to fill in the calculation domain. The typical length scale incorporated in the simulation is the average grain size, which has a value of 10 μm specifically chosen to match the block size observed in high-C martensitic bearing steels [47]. While martensite has a very hierarchical microstructure with the smallest unit being the submicron laths of a single variant, the main strengthening contribution is related to the block size [48]. The randomly assigned rotation angles (    2 , , ) of grains were assumed to be uniformly distributed from 0 to π/2 . The rotation of a cubic stiffness tensor from local material coordinates to the global coordinate system leading to a situation in which every grain has a different stiffness along the X, Y, or Z direction. Figure 7 shows the elastic modulus along the Z coordinate, computed as where (l, m, n) is an arbitrary lattice plane, and s 11 , s 12 , and s 44 are the elastic compliance constants for the cubic material [49]. This very wide distribution in stiffness between the grains can lead to local stress concentrations between grains with large orientation differences and therefore faster RCF. However, the average stress cannot be greater or less than the values computed using isotropic elasticity and Voigt and Reuss average constants, which act as upper and lower bounds, respectively.  Figures 9(a), 9(c), 9(e), and 9(g) clearly demonstrate that the LGR approach is sufficiently accurate to model polycrystalline anisotropic material. For the stress distributions shown in Figs. 9(b), 9(d), 9(f), and 9(h), values calculated with LGR agree very well with those obtained on a uniform grid with the finest grid mesh size. Note that the final LGR consists of the solution values in all grid points in the local area of the domain, which are in the locally finest grid. As long as the grid point is from the top grid level where the finest mesh is used (Z < 4.375 and X < 2.1875), the solutions are virtually identical to the full domain solution of the finest grid. For regions where the grid is coarser, differences occur, as in these regions, the grid is effectively coarser. Yet, the existence of a finer grid in the higher gradient regions also improves the solution in these areas. Nevertheless, in the regions not covered by the finest mesh, the stress distribution with LGR has some differences compared to that on the finest grid without using LGR. There are two reasons for explaining the differences: 1) The accuracy resolution is different for the two solutions. For a grid point on a coarse grid, the calculation accuracy is lower than that on a fine grid. 2) The variation of the material topology (i.e. the orientation and   dimensions of the grains) is comparatively less well represented with coarser grids. As a result, the final displacement, strain, and stress distributions have some differences in regions where only a coarser grid is used. This also explains why the stresses of both solutions agreed well for isotropic materials in the previous section.
Although both solutions have some differences in the regions not covered by the finest mesh for the LGR approach, these differences are relatively small. Moreover, the solution in these regions is due to the existence of the finer grid, which is more accurate than that it would have been when using only the coarser grid. This phenomenon can also be seen from the variation of the maximum VMS, as shown in Fig. 10. Here, the domains of LGR on grid levels 3 and 4 are the same as those shown in Table 3. On the top grid level (level 5), the range of X, Y, and Z decreases gradually by the same distance. The horizontal axis of the graph is the starting point of LGR in the X direction. The calculation domain on level 5 is [X-(−X), X-(−X), 0-(−2X)], which means that the size of the LGR volume decreases with the increase of X start . If X start equals −5, LGR is not employed. From the graphs shown in Fig. 10, it can be observed that the effect of LGR on the maximum VMS becomes large when its border is close to the location of the maximum VMS, and causes the largest deviation of 0.79% (Fig. 10(a) black line) and 1.4 % (Fig. 10(b) blue line) for dry and lubricated contacts, respectively.
The minimum film thickness computed for the EHL case is sensitive to the local mesh size when the so-called sidelobes where they occur are small, as in that case, a small mesh size is required to capture it accurately. For a detailed evaluation of the influence of the LGR region on the numerical results, the variation of the minimum film thickness as a function of LGR borders is shown in Fig. 11. It can be seen that if the LGR border is far away from the location of the minimum film thickness, the results are very close to those without LGR. When the LGR border is close to the actual location of the minimum film thickness location, some deviations can be observed. However, with the largest deviation of 1.4% compared to the results without using LGR for the studied case, it is concluded that the proposed MG method with LGR  can meet the accuracy requirement when the domains of interest, such as the location of the maximum VMS or the minimum film thickness, are covered by a fine mesh, and the point of interest is not close to the LGR border.
Finally, to demonstrate the robustness of the developed LGR approach, 51 cases were considered for different microstructures. In each microstructure domain, there were more than two million grains with random orientations. The LGR domains are the same as those listed in Table 3. The random variations of the crystallographic orientation and grain topology result in a 10% scatter of the (dimensionless) maximum VMS and the (dimensionless) minimum film thickness, as shown in Fig. 12. Note that the horizontal axis here has no physical meaning other than a case number, which represents a different distribution of the crystallographic orientation over the grains. The results with and without using LGR agreed well. The maximum error percentages are 0.063%, 0.046%, and 0.024% for Figs. 12(a)-12(c), respectively, which demonstrates the robustness and accuracy of the proposed LGR. Figure 13 shows the solution time as a function of FMG W-cycles. The Voronoi tessellation, size of the LGR patches, and other working conditions are the same as those in Section 3.2. The average ratios between the cases with and without LGR are 2.47  LGR is approximately eight. The large decrease in memory cost is caused by the fact that the number of grid points on fine grid levels is greatly reduced. The reduction ratio of the calculation time and memory costs depends on the size of the patches and the number of grains (isotropic, homogenous isotropic, and polycrystalline anisotropic) used in the calculation. It should be pointed out that although the adoption of the LGR decreases the time cost, it is still slower compared to the point-lubricated contact problem solvers, which use influence functions for the surface deformation in the solution process with fast evaluation of this deformation by FFT, e.g., presented in Refs. [14][15][16] for the displacement calculation of the layered and inhomogeneous materials, or in the solvers described in Ref. [35]. This is because the present algorithm is a full 3D solver of the displacement equations combined with the contact problem. A detailed comparison of computing time is given in Ref. [24]. However, when full 3D topological and crystallographic orientation data are to be evaluated, this can only be achieved through the present type of solvers.

Application example
The application of LGR enables previously developed MG methods to exhibit similar accuracy resolutions at reduced computational resources, or yields higher accuracy solutions with reduced additional com-putational cost. The ability of LGR to exploit varying resolutions for heterogeneous anisotropic materials, as presented in this paper, can also be used in combination with other aspects in contact mechanics and EHL. In this section, this is illustrated using two typical cases. The first case shown in Fig. 14 is a contact problem that includes the effect of a single surface roughness (a half-sine wave). A locally refined grid mesh is used around the location of the single surface roughness, and the result is shown in Fig. 14(b).
Compared to the result shown in Fig. 14(c) without LGR, the grain boundaries can be observed more clearly because of the locally refined mesh (without LGR: h x, y, z = 0.0195; with LGR: local h x, y, z = 0.0049). The second case is an elliptical contact problem. Elliptical contacts require a larger calculation domain in the direction of the major axis compared to that of the minor axis. Here, an ellipticity ratio (b/a) of 1.667 was used. The dimensionless calculation domain in the Y direction ranged from −10 to 10. The remaining calculation domains and average grain size were the same as before. For the dry contact and lubrication conditions, surface roughness was used to test the robustness of the developed solvers. The waviness of the roughness was the same for both cases. Under dry contact, the surface roughness results in some isolated contact regions (Fig. 15(a)), which is very harmful for bearing fatigue life as it causes a larger local stress concentration (Fig. 15(c)) close to the contact surface, and plastic deformation may occur if the local stress is sufficiently high. When a lubricating film separates the two contacting bodies, pressure www.Springer.com/journal/40544 | Friction fluctuations ( Fig. 15(b)) and local stress concentrations ( Fig. 15(d)) are also observed. This illustrates the importance of the surface finish for extending the fatigue life of rolling bearings [50].

Conclusions and future research
In this study, the use of local grid refinements embedded in multigrid techniques is adopted to further enhance the efficiency of the previously developed multigrid methods for displacement field solution in 3D with stress field analysis of polycrystalline anisotropic material under dry and lubricated point contacts. It has been proven that with the local grid refinement incorporated into previously developed programs [24][25][26], the accuracy of the results remains almost unaffected, while the solution time and memory storage requirements are further reduced. As effective non-uniformity is achieved by exploiting only uniform grid discretization formulations, the application of the local grid refinement is relatively straightforward, and makes the multigrid-based solvers more efficient and suitable for calculations using normal computers and for engineering applications such as material optimizations and multi-scale simulations. Of course, a further reduction in computing time can be achieved by parallelization.
In this study, the orientation angles of the grains were assumed to be random, and their distribution was uniform. This is more or less a worst case. In reality, grain orientations are not necessarily random, and can be affected by manufacturing. Using electronbackscatter diffraction techniques, the rotation angles of grains and compositions (Fe 3 C, martensite, and austenite) [51] in actual bearing steel can be obtained, which can serve as an input to the proposed algorithm. Together with the efficiency of the algorithm, this enables detailed and accurate predictions of the fatigue life, and brings actual computational diagnostics within engineering reach.

Appendix A: Voigt and Reuss average elastic constants
Here we provide equations for computing the Voigt and Reuss average elastic constants [44]. Considering a stiffness tensor c of a material with cubic symmetry expressed using matrix notation, the Voigt average  The average Young's modulus is then computed as Next, Reuss elastic constants are computed using material compliance tensor   1 s c , which then yields 11 As mentioned above, the average elastic constant can be computed for polycrystalline material if the condition for random grain orientation is satisfied. The Voigt and Reuss average constants act as upper and lower bounds, respectively, over the grain stiffness expressed in the global (loading) coordinate system. In general, the Voigt average provides a better estimate for random grains under constant strain, whereas the Reuss average is more suitable for grains under constant stress.
www.Springer.com/journal/40544 | Friction The FMG algorithm additionally makes use of the coarser grids to obtain a good initial approximation on the finest (target) grid. Starting from the coarsest grid level 1, multigrid cycles are performed, and their solutions are interpolated to supply an initial approximation for calculations performed on the next-finer grid. The resulting algorithm is depicted in Fig. 2 for three levels and W-cycles.
www.Springer.com/journal/40544 | Friction values of the maximum pressure, deformation, VMS, and minimum film thickness on the central XZ and YZ planes are given. The maximum percentage errors between the two grid resolutions were 2.67857% and 2.84587% for dry contact and lubrication, respectively.
Owing to the small difference between the results using two different mesh sizes, the grid resolution of 10/512 is used in this study.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.