Crack Propagation Analysis of Compression Loaded Rolling Elements

The problem of crack propagation from internal defects in thermoplastic cylindrical bearing elements is addressed in this paper. The crack propagation in these elements takes place under mixed-mode conditions—i.e., all three possible loading modes (tensile opening mode I and shear opening modes II and III) of the crack are combined together. Moreover, their mutual relation changes during the rotation of the element. The dependency of the stress intensity factors on the crack length was described by general parametric equations. The model was then modified by adding a void to simulate the presence of a manufacturing defect. It was found that the influence of the void on the stress intensity factor values is quite high, but it fades with crack propagating further from the void. The effect of the friction between the crack faces was find negligible on stress intensity factor values. The results presented in this paper can be directly used for the calculation of bearing elements lifetime without complicated finite element simulations.


Introduction
It has become quite common nowadays to use thermoplastic materials for the production of cylindrical rollers in roller bearings or linear guidance parts. Thermoplastics are often preferred for these applications because of their light weight and effective mass production [1]. Typical materials for these applications are polyamide (PA), polyoxymethylene (POM) and polyetheretherketone (PEEK) [1][2][3][4][5].
The cylindrical bearing elements are subjected to rolling in their operation. Parts subjected to rolling often fail due to cracks initiating on the surface (or slightly under the surface) that are in contact with another part (typically a metal rail, linear guide or a bearing ring). The mechanism is called rolling contact fatigue and it is usually associated with metal wheels running on rails [6][7][8][9], but it appears also in bearing elements such as balls and cylinders [10][11][12][13]. It is important to note that in the rolling contact fatigue the cracks grow under mixed-mode loading conditions, which means the shear modes II (in-plane shear) and III (out-of-plane shear) play a significant role in the crack propagation.
There is also another mechanism of failure of the bearing elements-the crack propagates from an internal defect until it reaches a critical length and the bearing element breaks. This mechanism is important for thermoplastic bearing elements, especially cylinders.

1.
Basic model with a flat crack-model of a cylinder compressed between two steel plates, that contains one circular flat crack in the middle. Although the faces of the crack can come in contact, no friction is considered between the faces. Different combinations of dimensions were considered.

2.
Model with a void-the crack in this model is not flat from the start, but it starts from a three-dimensional void in the middle of the cylinder. Although it starts from the void, the crack is still considered flat and circular. 3.
Model with a void and friction-the same as the model with a void, but friction is considered between the faces of the crack.
The details of these models are described further in this section. The results are discussed in Section 3.

Basic Model with a Flat Crack
The most important part of the model is a cylinder with a flat circular crack in the middle (the crack is modelled directly as a part of the geometry, refer to ANSYS APDL Help for more details about crack analysis [34]). Apart from the cylinder, two steel plates, one at the top and one at the bottom, compressing the cylinder were added to model the case of a cylinder installed between rails or bearing rings. The main dimensions of the basic model cylinder are denoted D (diameter) and L (length). The crack length is denoted a. The model is parametric, which means it is possible to simulate almost any set of dimensions and also the change in orientation of the crack during rotation of the cylinder. See the modelled situation in Figure 1. parameters-the stress intensity factors. The software ANSYS was used for the simulation. The model was created using a parametric macro programmed in the APDL (ANSYS Parametric Design Language). Three different modifications of the model, that differ in some features influencing the crack growth, were created: 1. Basic model with a flat crack-model of a cylinder compressed between two steel plates, that contains one circular flat crack in the middle. Although the faces of the crack can come in contact, no friction is considered between the faces. Different combinations of dimensions were considered. 2. Model with a void-the crack in this model is not flat from the start, but it starts from a three-dimensional void in the middle of the cylinder. Although it starts from the void, the crack is still considered flat and circular. 3. Model with a void and friction-the same as the model with a void, but friction is considered between the faces of the crack.
The details of these models are described further in this section. The results are discussed in Section 3.

Basic Model with a Flat Crack
The most important part of the model is a cylinder with a flat circular crack in the middle (the crack is modelled directly as a part of the geometry, refer to ANSYS APDL Help for more details about crack analysis [34]). Apart from the cylinder, two steel plates, one at the top and one at the bottom, compressing the cylinder were added to model the case of a cylinder installed between rails or bearing rings. The main dimensions of the basic model cylinder are denoted D (diameter) and L (length). The crack length is denoted a. The model is parametric, which means it is possible to simulate almost any set of dimensions and also the change in orientation of the crack during rotation of the cylinder. See the modelled situation in Figure 1. Only one half of the cylinder was modelled to take advantage of the symmetry and reduce the number of elements and computing time. Further reduction of the model was not possible, because other planes of symmetry are disrupted by the rotation of the crack inside the cylinder. Quadratic tetrahedral elements were used to mesh the most of the model. The mesh was made very fine and a little more regular around the crack front- Only one half of the cylinder was modelled to take advantage of the symmetry and reduce the number of elements and computing time. Further reduction of the model was not possible, because other planes of symmetry are disrupted by the rotation of the crack inside the cylinder. Quadratic tetrahedral elements were used to mesh the most of the model. The mesh was made very fine and a little more regular around the crack front-quadratic brick mesh was used in the area surrounding the crack front with special quarter-point elements directly at the crack front (see the mesh of the model with detailed view of the crack front in Figure 2). The material model of the cylinder was a linear elastic, quadratic brick mesh was used in the area surrounding the crack front with special quarter-point elements directly at the crack front (see the mesh of the model with detailed view of the crack front in Figure 2). The material model of the cylinder was a linear elastic, isotropic solid defined by the Young's modulus of 3.6 GPa and Poisson's ratio of 0.45. These values are typical for the POM used to make bearing elements. The steel plates compressing the cylinder at the top and at the bottom were meshed by regular quadratic bricks. The material model of these parts was linear elastic, isotropic, with Young's modulus of 210 GPa and Poisson's ratio of 0.3. Mutual contact was defined between the steel plates and the cylinder.
Boundary conditions were defined for the whole model. Symmetrical boundary conditions were defined on the plane of symmetry. Fixed support was defined at the bottom of the lower steel plate. The upper steel plate was loaded by the force of 350 N. The boundary conditions and loads are schematically depicted in Figure 2.
The growth of the crack inside the cylinder was simulated by modelling the cylinder with the crack length a going from 0.25 mm up to 1.75 mm by the steps of 0.50 mm with the addition of two extra crack lengths of 0.85 mm and 1.00 mm (this was done for a better comparison of the models with and without a void inside, which is discussed further in this paper). Note that, in this case, the crack is circular and the crack length a is its radius (as shown in Figure 1).
For every crack length, the stress intensity factors KI, KII and KIII were determined along the whole crack front. The position on the crack front is defined by the angle γ that goes from +90° across 0° to −90° (see Figure 1 for illustration).
The cylinder was considered rolling between the two compressing steel plates to simulate the operation of a bearing element. The rolling was achieved by changing the orientation of the crack in the model. It is important to remark, that contact was defined also between the faces of the crack, because the crack faces would be pressed against each other in some of the positions. The orientation of the crack is defined by the angle ρ between the plane of the crack and the vertical axis of the model (or the direction of forces acting on the steel plates, see Figure 3). The angle ρ goes from 0° up to 180° by the steps of 15°. For The steel plates compressing the cylinder at the top and at the bottom were meshed by regular quadratic bricks. The material model of these parts was linear elastic, isotropic, with Young's modulus of 210 GPa and Poisson's ratio of 0.3. Mutual contact was defined between the steel plates and the cylinder.
Boundary conditions were defined for the whole model. Symmetrical boundary conditions were defined on the plane of symmetry. Fixed support was defined at the bottom of the lower steel plate. The upper steel plate was loaded by the force of 350 N. The boundary conditions and loads are schematically depicted in Figure 2.
The growth of the crack inside the cylinder was simulated by modelling the cylinder with the crack length a going from 0.25 mm up to 1.75 mm by the steps of 0.50 mm with the addition of two extra crack lengths of 0.85 mm and 1.00 mm (this was done for a better comparison of the models with and without a void inside, which is discussed further in this paper). Note that, in this case, the crack is circular and the crack length a is its radius (as shown in Figure 1).
For every crack length, the stress intensity factors K I , K II and K III were determined along the whole crack front. The position on the crack front is defined by the angle γ that goes from +90 • across 0 • to −90 • (see Figure 1 for illustration).
The cylinder was considered rolling between the two compressing steel plates to simulate the operation of a bearing element. The rolling was achieved by changing the orientation of the crack in the model. It is important to remark, that contact was defined also between the faces of the crack, because the crack faces would be pressed against each other in some of the positions. The orientation of the crack is defined by the angle ρ between the plane of the crack and the vertical axis of the model (or the direction of forces acting on the steel plates, see Figure 3). The angle ρ goes from 0 • up to 180 • by the steps of 15 • . For a better description of the change in the values during rolling, eight steps of the angle ρ were added-3. The range 0 • -180 • simulates only a half turn of the cylinder, but the remaining part of the turn would be symmetrical. The crack growth from 0.25 to 1.75 mm and the stress intensity factor determination was carried out in every step of the rolling. More than 400 simulations in total were carried out. a better description of the change in the values during rolling, eight steps of the angle ρ were added-3.8°, 7.5°, 11.2° and 22.5° as well as 157.5°, 168.8°, 172.5° and 176.2. The range 0°-180° simulates only a half turn of the cylinder, but the remaining part of the turn would be symmetrical. The crack growth from 0.25 to 1.75 mm and the stress intensity factor determination was carried out in every step of the rolling. More than 400 simulations in total were carried out. In every step of the crack propagation (and in every version of the model), the stress intensity factors were evaluated using two different methods. The first method was the domain integral [35], which is fully integrated in the software ANSYS. The second method was the determination of stress intensity factors directly from the deformation of the nodes of the special crack tip elements [36].

Models with a Void and Friction
The models with the circular flat crack give an idea about the crack propagation and the mutual ratios of the crack propagation modes. However, the flat crack does not exactly reflect the reality of the bearing cylinders. The shrinkage defects, from which the cracks usually start, are not just simple flat discontinuities in the material. This was the reason to modify the basic model with a flat crack and create a model that contains a void in the middle. In the case of the models with a void and friction, the basic dimensions D × L were kept constant at 6 × 6 mm 2 . The variable parameters were the crack length, the orientation of the crack and the dimensions of the void.
In practice, the void was observed to be a rather complex three-dimensional, spheroidal structure with a strong variability in shape and surface roughness within the same class of bearing elements. Moreover, size, shape and surface topology depend significantly on processing parameters and the processed polymers. Hence, to keep it on a useful complexity level and in order to obtain generalizable data, the void was represented by a regular spheroid. The radius of the void rv was 0.75 mm (diameter dv = 1.5 mm). The height of the void hv was set as a variable to study its effect. See Figure 4 for the illustration of the situation.
To see what different shapes of the voids do to the stress distribution around the crack, the height of the void was varied in the model. The radius of the void was kept constant. The void height is described here in terms of a ratio of the height to the radius (further referred to as the void ratio)-hv/rv. The void ratio values were considered going from 0.125 up to 1.75 (with the step of 0.25 between 0.25 and 1.75). The void is shaped such as a perfect sphere in case of hv/rv = 1. In case of hv/rv < 1, the void is an oblate spheroid; in case of hv/rv > 1, the void has a shape of a prolate spheroid. Apart from the void, the model contains a circular crack similar to the model with a flat crack. The crack lengths a considered in the model with the void were the same as in the model with a flat crack, except the 0.25 mm and 0.75 mm. It was not possible to model these lengths due to the void in the middle. The crack length started at 0.85 mm and went up to 1.75 mm. In every step of the crack propagation (and in every version of the model), the stress intensity factors were evaluated using two different methods. The first method was the domain integral [35], which is fully integrated in the software ANSYS. The second method was the determination of stress intensity factors directly from the deformation of the nodes of the special crack tip elements [36].

Models with a Void and Friction
The models with the circular flat crack give an idea about the crack propagation and the mutual ratios of the crack propagation modes. However, the flat crack does not exactly reflect the reality of the bearing cylinders. The shrinkage defects, from which the cracks usually start, are not just simple flat discontinuities in the material. This was the reason to modify the basic model with a flat crack and create a model that contains a void in the middle. In the case of the models with a void and friction, the basic dimensions D × L were kept constant at 6 × 6 mm 2 . The variable parameters were the crack length, the orientation of the crack and the dimensions of the void.
In practice, the void was observed to be a rather complex three-dimensional, spheroidal structure with a strong variability in shape and surface roughness within the same class of bearing elements. Moreover, size, shape and surface topology depend significantly on processing parameters and the processed polymers. Hence, to keep it on a useful complexity level and in order to obtain generalizable data, the void was represented by a regular spheroid. The radius of the void r v was 0.75 mm (diameter d v = 1.5 mm). The height of the void h v was set as a variable to study its effect. See Figure 4 for the illustration of the situation.
To see what different shapes of the voids do to the stress distribution around the crack, the height of the void was varied in the model. The radius of the void was kept constant. The void height is described here in terms of a ratio of the height to the radius (further referred to as the void ratio)-h v /r v . The void ratio values were considered going from 0.125 up to 1.75 (with the step of 0.25 between 0.25 and 1.75). The void is shaped such as a perfect sphere in case of h v /r v = 1. In case of h v /r v < 1, the void is an oblate spheroid; in case of h v /r v > 1, the void has a shape of a prolate spheroid. Apart from the void, the model contains a circular crack similar to the model with a flat crack. The crack lengths a considered in the model with the void were the same as in the model with a flat crack, except the 0.25 mm and 0.75 mm. It was not possible to model these lengths due to the void in the middle. The crack length started at 0.85 mm and went up to 1.75 mm.  The model with the void was further modified by adding friction between the contact faces of the crack. It was assumed that the friction could slightly change the shear stress distribution in some of the crack positions, which could have an effect on the stress intensity factors KII and KIII. The coefficient of friction used in the contact was 0.32. This value was chosen as a conservative estimate based on relevant data about friction coefficients of POM and PEEK materials, which are usually peeking slightly above 0.3 at room temperature [37][38][39][40]. However, the friction coefficient can significantly vary for different blends of the same polymer.

Results and Discussion
At first, the stress intensity factors determined by the model with a flat crack were compared to values published by Berer et al. in [3] to validate the functionality of the new model described in the previous section. The numerical model by Berer et al. represents one eighth of the bearing cylinder and it was not used to determine the values of stress intensity factors KII and KIII, only KI was evaluated. The values for the case of the cylinder with D × L = 6 × 6 mm 2 in the rolling position ρ = 0° were taken from the paper and compared to the results that were produced by the model described here. The comparison is plotted in Figure 5. The values of KI are plotted as a function of the position on the crack front described by the angle γ. The discrepancy between the results of the two models is negligible. Note that Berer et al.'s values are in the range of 0°-90°, because the one-eighthtype of symmetry was used in their model. It was possible, because the values of KI are symmetrical, as shown by the newer results. However, it would not be possible to use this type of symmetry for the evaluation of KII and KIII.
It was mentioned above, that the stress intensity factor values were evaluated by two different methods in the model described here. The results produced by both these methods were also compared and they are plotted in Figure 5. The solid line in Figure 5 represents the domain integral method; the separate points (in the color of the line) represent the estimation from the node deformations. There is a nice match between the values produced by the two different methods, which also validates the functionality of the model. All the values presented further in this paper were obtained by the domain integral method.
In the following sections of this paper, results from the 3 versions of the model are presented and discussed. The model with the void was further modified by adding friction between the contact faces of the crack. It was assumed that the friction could slightly change the shear stress distribution in some of the crack positions, which could have an effect on the stress intensity factors K II and K III . The coefficient of friction used in the contact was 0.32. This value was chosen as a conservative estimate based on relevant data about friction coefficients of POM and PEEK materials, which are usually peeking slightly above 0.3 at room temperature [37][38][39][40]. However, the friction coefficient can significantly vary for different blends of the same polymer.

Results and Discussion
At first, the stress intensity factors determined by the model with a flat crack were compared to values published by Berer et al. in [3] to validate the functionality of the new model described in the previous section. The numerical model by Berer et al. represents one eighth of the bearing cylinder and it was not used to determine the values of stress intensity factors K II and KIII , only K I was evaluated. The values for the case of the cylinder with D × L = 6 × 6 mm 2 in the rolling position ρ = 0 • were taken from the paper and compared to the results that were produced by the model described here. The comparison is plotted in Figure 5. The values of K I are plotted as a function of the position on the crack front described by the angle γ. The discrepancy between the results of the two models is negligible. Note that Berer et al.'s values are in the range of 0 • -90 • , because the one-eighth-type of symmetry was used in their model. It was possible, because the values of K I are symmetrical, as shown by the newer results. However, it would not be possible to use this type of symmetry for the evaluation of K II and K III .
It was mentioned above, that the stress intensity factor values were evaluated by two different methods in the model described here. The results produced by both these methods were also compared and they are plotted in Figure 5. The solid line in Figure 5 represents the domain integral method; the separate points (in the color of the line) represent the estimation from the node deformations. There is a nice match between the values produced by the two different methods, which also validates the functionality of the model. All the values presented further in this paper were obtained by the domain integral method.
In the following sections of this paper, results from the 3 versions of the model are presented and discussed.

Model with a Flat Crack-Results and Discussion
It was mentioned above that the orientation of the crack to the acting load changes due to the rolling of the cylinder. The pure mode I appears only in a few instances during the turn. There is always rather a combination of modes I, II and III. The modes and their ratio change depending on the position of the cylinder (angle ρ) and they are also different for different positions on the crack front (given by the angle γ).
The results of the basic model with a flat crack provide an idea on how the modes change during one turn of the cylinder. The following figures contain plotted results for the cylinder with the dimensions of D × L = 6 × 6 mm 2 and loading force F = 350 N, unless stated otherwise. The results are plotted in Figure 6 (KI), Figure 7 (KII) and Figure 8 (KIII) for one of the simulated crack lengths-1.25 mm. 3D plots were chosen to visualize the values of stress intensity factors depending on both, the position on the crack front γ and the overall rolling orientation ρ of the crack. The results were also plotted in the form of 2D plots, where the crack length a and position on the crack front γ were fixed and the stress intensity factors are plotted as a function of the rolling orientation (angle ρ). Many 2D plots had to be created to illustrate the whole situation, because of many possible combinations of parameters (crack length and position on the crack front). The 2D plots are not included in the text of this paper for the sake of clarity. The most important 2D plots are included in Appendix A.
The 3D plot of KI ( Figure 6) shows that the KI values reach their maximum at the beginning and at the end of the turn of the cylinder. During the turn from 0° towards 180°, the values decrease up to the point when the two crack faces come into contact. The crack stays closed until the cylinder comes into a position where the opening stress starts acting on the crack again. The KI values are zero when the crack faces are in contact in the model. A simulation without the contact of crack faces was also carried out to investigate the exact moment of the crack closing and opening. If no contact is defined between the crack faces, the KI values become negative in the part of the cycle where the crack is closed (the negative values are also plotted in Figure 6). Negative values of KI cannot occur in reality, but this kind of simulation helps to evaluate the cycle of KI and its asymmetry, which can be important for a later use in lifetime estimations and for experimental testing of such a situation.
The points of crack closing and opening are slightly different for different positions on the crack front-for the plotted crack length of 1.25 mm, the crack closes a little earlier for γ = +90° (and −90°) than for γ = 0° (this can be better observed in Figure A5 in Appendix A). Their position also depends on the current crack length (the shift in the position can  [3] and from the newer model described in this paper. The values from the newer model were estimated by domain integral (solid lines) and node deformations (separate points).

Model with a Flat Crack-Results and Discussion
It was mentioned above that the orientation of the crack to the acting load changes due to the rolling of the cylinder. The pure mode I appears only in a few instances during the turn. There is always rather a combination of modes I, II and III. The modes and their ratio change depending on the position of the cylinder (angle ρ) and they are also different for different positions on the crack front (given by the angle γ).
The results of the basic model with a flat crack provide an idea on how the modes change during one turn of the cylinder. The following figures contain plotted results for the cylinder with the dimensions of D × L = 6 × 6 mm 2 and loading force F = 350 N, unless stated otherwise. The results are plotted in Figure 6 (K I ), Figure 7 (K II ) and Figure 8 (K III ) for one of the simulated crack lengths-1.25 mm. 3D plots were chosen to visualize the values of stress intensity factors depending on both, the position on the crack front γ and the overall rolling orientation ρ of the crack. The results were also plotted in the form of 2D plots, where the crack length a and position on the crack front γ were fixed and the stress intensity factors are plotted as a function of the rolling orientation (angle ρ). Many 2D plots had to be created to illustrate the whole situation, because of many possible combinations of parameters (crack length and position on the crack front). The 2D plots are not included in the text of this paper for the sake of clarity. The most important 2D plots are included in Appendix A.
The 3D plot of K I ( Figure 6) shows that the K I values reach their maximum at the beginning and at the end of the turn of the cylinder. During the turn from 0 • towards 180 • , the values decrease up to the point when the two crack faces come into contact. The crack stays closed until the cylinder comes into a position where the opening stress starts acting on the crack again. The K I values are zero when the crack faces are in contact in the model. A simulation without the contact of crack faces was also carried out to investigate the exact moment of the crack closing and opening. If no contact is defined between the crack faces, the K I values become negative in the part of the cycle where the crack is closed (the negative values are also plotted in Figure 6). Negative values of K I cannot occur in reality, but this kind of simulation helps to evaluate the cycle of K I and its asymmetry, which can be important for a later use in lifetime estimations and for experimental testing of such a situation.
The points of crack closing and opening are slightly different for different positions on the crack front-for the plotted crack length of 1.25 mm, the crack closes a little earlier for γ = +90 • (and −90 • ) than for γ = 0 • (this can be better observed in Figure A5 in Appendix A). Their position also depends on the current crack length (the shift in the position can be observed in Figures A1-A6 in Appendix A), which suggests that the crack might not propagate exactly as a regular circle in reality. be observed in Figures A1-A6 in Appendix A), which suggests that the crack might not propagate exactly as a regular circle in reality. The values of stress intensity factors for the shear modes, the KII and KIII, are higher in magnitude compared to the KI in terms of maximum values. On the beginning of the turn, both KII and KIII are zero along the whole crack front, as the crack is not subjected to shear loading at all. However, the conditions change with the turning. The middle of the crack (γ = 0°) is subjected to mode III type of loading and the mode II does not appear here at all during the turn, whereas the crack front ends (γ = 90° and −90°) develop mutually opposite values of KII during the cycle, and KIII remains equal to zero. In between these positions, the KII and KIII values follow different sine patterns. As the crack becomes perpendicular to the direction of loading (ρ = 90°), both of the shear modes disappear again. Then, in the following part of the turn, the values of KII and KIII appear again in the same places on the crack front, but with opposite signs-see   The values of stress intensity factors for the shear modes, the K II and K III , are higher in magnitude compared to the K I in terms of maximum values. On the beginning of the turn, both K II and K III are zero along the whole crack front, as the crack is not subjected to shear loading at all. However, the conditions change with the turning. The middle of the crack (γ = 0 • ) is subjected to mode III type of loading and the mode II does not appear here at all during the turn, whereas the crack front ends (γ = 90 • and −90 • ) develop mutually opposite values of K II during the cycle, and K III remains equal to zero. In between these positions, the K II and K III values follow different sine patterns. As the crack becomes perpendicular to the direction of loading (ρ = 90 • ), both of the shear modes disappear again. Then, in the following part of the turn, the values of K II and K III appear again in the same places on the crack front, but with opposite signs-see be observed in Figures A1-A6 in Appendix A), which suggests that the crack might not propagate exactly as a regular circle in reality. The values of stress intensity factors for the shear modes, the KII and KIII, are higher in magnitude compared to the KI in terms of maximum values. On the beginning of the turn, both KII and KIII are zero along the whole crack front, as the crack is not subjected to shear loading at all. However, the conditions change with the turning. The middle of the crack (γ = 0°) is subjected to mode III type of loading and the mode II does not appear here at all during the turn, whereas the crack front ends (γ = 90° and −90°) develop mutually opposite values of KII during the cycle, and KIII remains equal to zero. In between these positions, the KII and KIII values follow different sine patterns. As the crack becomes perpendicular to the direction of loading (ρ = 90°), both of the shear modes disappear again. Then, in the following part of the turn, the values of KII and KIII appear again in the same places on the crack front, but with opposite signs-see   The crack tip mixed-mode loading has a completely out-of-phase character. Again, refer to Figures A1-A6 in Appendix A for a more detailed view on the changes of K-values during the turn. When the KI reaches the maximum values, both KII and KIII are zero. There are short intervals at the beginning and the end of the turn, when the crack is open and loaded by a combination of KI and KIII (for the position of γ = 0°) or KI and KII (for γ = 90° and −90°). Between these, there are intervals in which the KII and KIII are non-zero and even reach their maxima (or minima) and the crack is closed (KI is zero). This means that the crack faces are being forced against each other and into one of the shear modes at the same time. It is quite likely that heat is generated by the friction of the crack faces in these intervals, which can have an influence on the crack propagation rate [31,32,41].
The instant values of K are not very practical for the description of the crack growth kinetics in the investigated situation. It is more practical to use the maximum values of stress intensity factor Kmax to describe the whole cycle. In Figure 9, these values are plotted as a function of the normalized crack length a/W (where a is the crack length and W corresponds to the radius or half of the length of the whole cylinder depending on position γ = 90° (−90°) or γ = 0°, respectively (see Figure 1). It is important to note here, that the rolling position ρ, at which the maxima and minima of KI and KIII are reached, are constant with the growing crack length a. The maximum of KI can be always found at ρ = 0° and ρ = 180°, the (theoretical) minimum at ρ = 90°. The maximum of KIII stays at ρ = 135° and minimum at ρ = 45°. However, the position where the KII reaches its extreme (minimum or maximum depending on the position, if γ = 90° or −90°) gradually shifts from 45° towards lower values of ρ with the crack length increasing and similarly the position of the other extreme shifts from 135° towards higher values. The shift can be observed in Figures  A1-A6 in the Appendix A. The cause of this is most likely that the crack becomes more influenced by a complicated stress state in the vicinity of the contact with the loading plates, which manifests itself the most at the positions γ = 90° and −90°, where KII reaches its maxima and minima. This shift in the position does not influence the characterization of the stress intensity factor cycles using the Kmax values though.
However, the magnitude of the stress intensity factor itself does not provide information about the character of the cycle. For the correct description of the loading cycle, it is necessary to specify also the asymmetry of the cycle in terms of R-ratio. The R-ratio is the ratio of the minimum value Kmin to the maximum value Kmax of the cycle.
In the investigated case, the R-ratio is dependent on the position at the crack front and also on the current crack length. The dependency is plotted in Figure 10. In case of mode I, the R-ratio is negative (this means that the minimum value is negative and the maximum is positive); it is lower in the 0° position on the crack front, where the minimum The crack tip mixed-mode loading has a completely out-of-phase character. Again, refer to Figures A1-A6 in Appendix A for a more detailed view on the changes of Kvalues during the turn. When the K I reaches the maximum values, both K II and K III are zero. There are short intervals at the beginning and the end of the turn, when the crack is open and loaded by a combination of K I and K III (for the position of γ = 0 • ) or K I and K II (for γ = 90 • and −90 • ). Between these, there are intervals in which the K II and K III are non-zero and even reach their maxima (or minima) and the crack is closed (K I is zero). This means that the crack faces are being forced against each other and into one of the shear modes at the same time. It is quite likely that heat is generated by the friction of the crack faces in these intervals, which can have an influence on the crack propagation rate [31,32,41].
The instant values of K are not very practical for the description of the crack growth kinetics in the investigated situation. It is more practical to use the maximum values of stress intensity factor K max to describe the whole cycle. In Figure 9, these values are plotted as a function of the normalized crack length a/W (where a is the crack length and W corresponds to the radius or half of the length of the whole cylinder depending on position γ = 90 • (−90 • ) or γ = 0 • , respectively (see Figure 1). It is important to note here, that the rolling position ρ, at which the maxima and minima of K I and K III are reached, are constant with the growing crack length a. The maximum of K I can be always found at ρ = 0 • and ρ = 180 • , the (theoretical) minimum at ρ = 90 • . The maximum of K III stays at ρ = 135 • and minimum at ρ = 45 • . However, the position where the K II reaches its extreme (minimum or maximum depending on the position, if γ = 90 • or −90 • ) gradually shifts from 45 • towards lower values of ρ with the crack length increasing and similarly the position of the other extreme shifts from 135 • towards higher values. The shift can be observed in Figures A1-A6 in the Appendix A. The cause of this is most likely that the crack becomes more influenced by a complicated stress state in the vicinity of the contact with the loading plates, which manifests itself the most at the positions γ = 90 • and −90 • , where K II reaches its maxima and minima. This shift in the position does not influence the characterization of the stress intensity factor cycles using the K max values though.
However, the magnitude of the stress intensity factor itself does not provide information about the character of the cycle. For the correct description of the loading cycle, it is necessary to specify also the asymmetry of the cycle in terms of R-ratio. The R-ratio is the ratio of the minimum value K min to the maximum value K max of the cycle.
In the investigated case, the R-ratio is dependent on the position at the crack front and also on the current crack length. The dependency is plotted in Figure 10. In case of mode I, the R-ratio is negative (this means that the minimum value is negative and the maximum is positive); it is lower in the 0 • position on the crack front, where the minimum has a lower value, and it gets slightly higher towards the 90 • (and −90 • ) position. The ratios in all the positions are rising with the growing crack. Additionally, the difference between the positions becomes more pronounced with the crack growing larger. The reason probably is that the position of γ = 90 • (−90 • ) starts to be more influenced by the stress distribution caused by the contact zone of bearing element and steel plates, when the crack length a grows larger. In case of mode II, the R-ratio of −1 (this denotes a symmetrical cycle) is the same for every position on the crack front, apart from the 0 • , where the K II is 0. Analogically to the mode II, the R-ratio for the mode III is equal to −1, apart from the 90 • and −90 • positions. has a lower value, and it gets slightly higher towards the 90° (and −90°) position. The ratios in all the positions are rising with the growing crack. Additionally, the difference between the positions becomes more pronounced with the crack growing larger. The reason probably is that the position of γ = 90° (−90°) starts to be more influenced by the stress distribution caused by the contact zone of bearing element and steel plates, when the crack length a grows larger. In case of mode II, the R-ratio of −1 (this denotes a symmetrical cycle) is the same for every position on the crack front, apart from the 0°, where the KII is 0. Analogically to the mode II, the R-ratio for the mode III is equal to −1, apart from the 90°and −90° positions. The Kmax values were determined for more combinations of dimensions D × L. The considered diameters D were 3, 4, 5, 6 mm, the considered lengths L were 3, 4, 5, 6, 9 and 12 mm. The entire range of crack lengths, as it is specified in Section 2, was considered for the 6 × 6 mm 2 type of cylinder only. Only some crack lengths were considered for the other combinations of length and diameter. These crack lengths were chosen with respect to the dimensions of the particular combination, because some of the lengths did not fit in the particular combination.
The obtained values of Kmax were then fitted with parametric functions that define the dependency of the stress intensity factors on the crack length and are generalized with respect to the loading force and dimensions of the cylinder. The fits were carried out for every point on the crack face, where the stress intensity factor reaches its maximum during the turn-this means there is one fit for the KIImax in the γ = 90° position and one for the KIIImax in the γ = 0° position. Two fits for the KImax were made, one for the 90° position and  has a lower value, and it gets slightly higher towards the 90° (and −90°) position. The ratios in all the positions are rising with the growing crack. Additionally, the difference between the positions becomes more pronounced with the crack growing larger. The reason probably is that the position of γ = 90° (−90°) starts to be more influenced by the stress distribution caused by the contact zone of bearing element and steel plates, when the crack length a grows larger. In case of mode II, the R-ratio of −1 (this denotes a symmetrical cycle) is the same for every position on the crack front, apart from the 0°, where the KII is 0. Analogically to the mode II, the R-ratio for the mode III is equal to −1, apart from the 90°and −90° positions.  The Kmax values were determined for more combinations of dimensions D × L. The considered diameters D were 3, 4, 5, 6 mm, the considered lengths L were 3, 4, 5, 6, 9 and 12 mm. The entire range of crack lengths, as it is specified in Section 2, was considered for the 6 × 6 mm 2 type of cylinder only. Only some crack lengths were considered for the other combinations of length and diameter. These crack lengths were chosen with respect to the dimensions of the particular combination, because some of the lengths did not fit in the particular combination.
The obtained values of Kmax were then fitted with parametric functions that define the dependency of the stress intensity factors on the crack length and are generalized with respect to the loading force and dimensions of the cylinder. The fits were carried out for every point on the crack face, where the stress intensity factor reaches its maximum during the turn-this means there is one fit for the KIImax in the γ = 90° position and one for the KIIImax in the γ = 0° position. Two fits for the KImax were made, one for the 90° position and The K max values were determined for more combinations of dimensions D × L. The considered diameters D were 3, 4, 5, 6 mm, the considered lengths L were 3, 4, 5, 6, 9 and 12 mm. The entire range of crack lengths, as it is specified in Section 2, was considered for the 6 × 6 mm 2 type of cylinder only. Only some crack lengths were considered for the other combinations of length and diameter. These crack lengths were chosen with respect to the dimensions of the particular combination, because some of the lengths did not fit in the particular combination.
The obtained values of K max were then fitted with parametric functions that define the dependency of the stress intensity factors on the crack length and are generalized with respect to the loading force and dimensions of the cylinder. The fits were carried out for every point on the crack face, where the stress intensity factor reaches its maximum during the turn-this means there is one fit for the K IImax in the γ = 90 • position and one for the K IIImax in the γ = 0 • position. Two fits for the K Imax were made, one for the 90 • position and another for the 0 • position, because the difference between the maxima in these two positions is not very pronounced, although technically the global maximum is reached only in the 0 • position. The position γ on the crack front is indicated by respective indices.
The functions are slightly different for every type of the stress intensity factor, but they always feature the loading force F, the square route of the crack length a, one of the dimensions of the cylinder (either the diameter D or the length L) and a dimensionless shape function Y that was found by curve fitting.
The shape function curve fits for the K Imax 90 • and K Imax 0 • are plotted in Figure 11a,b, respectively. The equation describing the dependencies are the following: where F is the loading force in N, a is the crack length in mm, D is the cylinder diameter in mm, L is the cylinder length in mm and Y I 90 • and Y I 0 • are the dimensionless shape functions that have been found in the following form: another for the 0° position, because the difference between the maxima in these two positions is not very pronounced, although technically the global maximum is reached only in the 0° position. The position γ on the crack front is indicated by respective indices. The functions are slightly different for every type of the stress intensity factor, but they always feature the loading force F, the square route of the crack length a, one of the dimensions of the cylinder (either the diameter D or the length L) and a dimensionless shape function Y that was found by curve fitting.
The shape function curve fits for the KImax 90° and KImax 0° are plotted in Figure 11a,b, respectively. The equation describing the dependencies are the following: where F is the loading force in N, a is the crack length in mm, D is the cylinder diameter in mm, L is the cylinder length in mm and YI 90° and YI 0° are the dimensionless shape functions that have been found in the following form: Even though the dimensions are in mm, the resultant unit of the stress intensity factors calculated using Equations (1) and (2) is in MPa·m 1/2 , which is the typical unit used for the stress intensity factors. This is ensured by the factor of √10 in the denominator of both equations.
The scatter of the points around the line in Figure 11 determines the level of generalization, which is achieved by the fit, and the difference between the parametric functions and the actual determined values. In cases of YImax 0° and YImax 90°, the difference is usually not more than 6%, but for some combinations it goes up to 30% (especially when there is a larger difference between D and L of the cylinder).
To describe the cycle of KI properly, the KImin values are also needed, because the Rratio does not stay constant for the KI cycle during the crack propagation (as illustrated in Figure 10). Even though the KImin values are only theoretical, because in practice the crack Even though the dimensions are in mm, the resultant unit of the stress intensity factors calculated using Equations (1) and (2) is in MPa·m 1/2 , which is the typical unit used for the stress intensity factors. This is ensured by the factor of √ 10 3 in the denominator of both equations.
The scatter of the points around the line in Figure 11 determines the level of generalization, which is achieved by the fit, and the difference between the parametric functions and the actual determined values. In cases of Y Imax 0 • and Y Imax 90 • , the difference is usually not more than 6%, but for some combinations it goes up to 30% (especially when there is a larger difference between D and L of the cylinder).
To describe the cycle of K I properly, the K Imin values are also needed, because the R-ratio does not stay constant for the K I cycle during the crack propagation (as illustrated in Figure 10). Even though the K Imin values are only theoretical, because in practice the crack is closed and the stress intensity factor is equal to zero, knowing these values makes it possible to describe the whole cycle in detail and most importantly to determine the precise moments of the crack closing and opening. The K Imin functions were created in the same manner as the K Imax functions above. The shape function fits are plotted in Figure 12a,b. The equations follow: The Y Imin 90 • and Y Imin 0 • are the dimensionless shape functions that have been found in the following form: is closed and the stress intensity factor is equal to zero, knowing these values makes it possible to describe the whole cycle in detail and most importantly to determine the precise moments of the crack closing and opening. The KImin functions were created in the same manner as the KImax functions above. The shape function fits are plotted in Figure  12a,b. The equations follow: The YImin 90° and YImin 0° are the dimensionless shape functions that have been found in the following form: The parametric functions were also found for the KIImax 90° and KIIImax 0° values. The fits for these values are plotted in Figure 13a,b, respectively. The generalized KIImax and KIIImax are much less scattered, which means that the parametric function provides a very good estimation of the real stress intensity factor values. The parametric functions have a similar form to the previous functions of KImax. The equations are the following: The equations describing functions YII 90° and YIII 0° were found to be the following: The parametric functions were also found for the K IImax 90 • and K IIImax 0 • values. The fits for these values are plotted in Figure 13a,b, respectively. The generalized K IImax and K IIImax are much less scattered, which means that the parametric function provides a very good estimation of the real stress intensity factor values. The parametric functions have a similar form to the previous functions of K Imax . The equations are the following: The equations describing functions Y II 90 • and Y III 0 • were found to be the following: Again, all the input dimensions in Equations (5)  All of the above stated equations are valid only in the range of 0 < 2a/D < 0.6 and 0 < 2a/L < 0.6, respectively, because the stress intensity factors were not evaluated outside this range. Care must be taken also if using the equations for a similar situation with different material parameters. Although the material parameters do not influence the stress intensity factors directly, they can influence the contact zone. This means that stress intensity values for larger crack lengths a may be influenced by the material parameters. However, the rest of the cases will not be influenced by a change in the material parameters, because the stress intensity factors in these cases only depend on the loading force and dimensions of the cylinders.

Models with Voids and Friction-Results and Discussion
This section summarizes the results from those versions of the model that contained an idealized void. All of the models with the voids were cylinders with D × L = 6 × 6 mm 2 and with a spheroidal void in the middle. The radius of the void rv was 0.75 mm and the height of the void hv was variable. In the following, the void dimensions are described by the void ratio hv/rv.
The general influence of the dimensions of the void on the stress intensity factors is illustrated by the results of a model with constant crack length of 0.85 mm and different void heights hv. Figure 14a   All of the above stated equations are valid only in the range of 0 < 2a/D < 0.6 and 0 < 2a/L < 0.6, respectively, because the stress intensity factors were not evaluated outside this range. Care must be taken also if using the equations for a similar situation with different material parameters. Although the material parameters do not influence the stress intensity factors directly, they can influence the contact zone. This means that stress intensity values for larger crack lengths a may be influenced by the material parameters. However, the rest of the cases will not be influenced by a change in the material parameters, because the stress intensity factors in these cases only depend on the loading force and dimensions of the cylinders.

Models with Voids and Friction-Results and Discussion
This section summarizes the results from those versions of the model that contained an idealized void. All of the models with the voids were cylinders with D × L = 6 × 6 mm 2 and with a spheroidal void in the middle. The radius of the void r v was 0.75 mm and the height of the void h v was variable. In the following, the void dimensions are described by the void ratio h v /r v .
The general influence of the dimensions of the void on the stress intensity factors is illustrated by the results of a model with constant crack length of 0.85 mm and different void heights h v . Figure 14a shows the K Imax values determined for the crack front positions γ = 90 • (−90 • ) and γ = 0 • . For the void ratio (h v /r v ) of 0.125, the K Imax values in both γ positions are the same as for the flat crack with the same crack length (a = 0.85 mm). After that, K Imax 90 • increases with increasing void ratio until the ratio of 1.75, where it is approximately three times higher. K Imax 0 • on the other hand, decreases to lower values than those calculated for the flat crack until it becomes zero. This is caused by the decrease in stiffness of the cylinder, when the void gets more prolate (h v /r v > 1)-the crack then  Figure 14b. The values of KIImax 90° and KIIImax 0° for the case of a flat crack of the same length a = 0.85 mm are plotted as solid lines for comparison. The KIImax 90° values for the void ratio of 0.125 are slightly lower than the KIImax 90° values for the flat crack length with the same crack length-by about 11%. With increasing void ratio, the KIImax 90° values decrease quite rapidly. For the void ratio of 1.75, they are approximately seven times lower than the values without the void. This is exactly the opposite tendency compared to KImax 90° in the same γ position. KIIImax 0° shows only a weak dependency on the void ratio. The KIIImax 0° for the void ratio of 0.125 is the same as for the flat crack. For the highest simulated void ratio of 1.75, the KIIImax 0° values are only slightly higher than those for the same flat crack length (by about 10%).
The results above suggest a strong influence of the void on the stress distribution around the crack tip (described by stress intensity factors). To see how the change in the void dimensions influences the actual crack propagation, crack propagation with the void was simulated. The void ratio of 0.25 (oblate spheroid) and 1 (perfect sphere) were simulated-these were chosen arbitrarily so that it was possible to describe a trend in the behavior. As was mentioned above, the crack growth from 0.85 mm up to 1.75 mm was considered. Figure 15 shows the values of Kmax from the flat crack model compared to the values from the models with the voids. The results agree with previous observations on the dependencies on the varying void ratio. The presence of the void increases the KI in the 90° position and decreases the KI in the 0° position. The KII is decreased substantially by the void, whereas KIII is subject to a mild increase. One feature in the results is common for all the stress intensity factors (no matter the loading mode nor the positions on the crack front): with the crack propagating further from the void, the stress intensity factors are less influenced by the presence of the void. This is given mainly by the stress concentration  Figure 14b. The values of K IImax 90 • and K IIImax 0 • for the case of a flat crack of the same length a = 0.85 mm are plotted as solid lines for comparison. The K IImax 90 • values for the void ratio of 0.125 are slightly lower than the K IImax 90 • values for the flat crack length with the same crack length-by about 11%. With increasing void ratio, the K IImax 90 • values decrease quite rapidly. For the void ratio of 1.75, they are approximately seven times lower than the values without the void. This is exactly the opposite tendency compared to K Imax 90 • in the same γ position. K IIImax 0 • shows only a weak dependency on the void ratio. The K IIImax 0 • for the void ratio of 0.125 is the same as for the flat crack. For the highest simulated void ratio of 1.75, the K IIImax 0 • values are only slightly higher than those for the same flat crack length (by about 10%).
The results above suggest a strong influence of the void on the stress distribution around the crack tip (described by stress intensity factors). To see how the change in the void dimensions influences the actual crack propagation, crack propagation with the void was simulated. The void ratio of 0.25 (oblate spheroid) and 1 (perfect sphere) were simulated-these were chosen arbitrarily so that it was possible to describe a trend in the behavior. As was mentioned above, the crack growth from 0.85 mm up to 1.75 mm was considered. Figure 15 shows the values of K max from the flat crack model compared to the values from the models with the voids. The results agree with previous observations on the dependencies on the varying void ratio. The presence of the void increases the K I in the 90 • position and decreases the K I in the 0 • position. The K II is decreased substantially by the void, whereas K III is subject to a mild increase. One feature in the results is common for all the stress intensity factors (no matter the loading mode nor the positions on the crack front): with the crack propagating further from the void, the stress intensity factors are less influenced by the presence of the void. This is given mainly by the stress concentration effect of the void, which is most important for short cracks initiated from void. The influence on the stress intensity factors for the shorter cracks is more pronounced with higher void ratios h v /r v . effect of the void, which is most important for short cracks initiated from void. The influence on the stress intensity factors for the shorter cracks is more pronounced with higher void ratios hv/rv.  Table 0. mm, hv/rv = 0.25 or 1).
The influence of the void on the asymmetry of the cycle R is the same as in the case without the void-the asymmetry in modes II and III is not affected by the crack propagation, the asymmetry in mode I increases slightly with longer cracks.
The results of the third version of the model that contains the void and also includes friction between the crack faces are plotted in Figure 16  The cycle of stress intensity factor KI has remained unchanged compared to the model with a void without friction. It is obvious, since the friction between crack faces cannot influence the tensile opening mode. The overall nature of the cycles of KII and KIII with friction is also very similar to the previous model without friction. The friction affects The influence of the void on the asymmetry of the cycle R is the same as in the case without the void-the asymmetry in modes II and III is not affected by the crack propagation, the asymmetry in mode I increases slightly with longer cracks.
The results of the third version of the model that contains the void and also includes friction between the crack faces are plotted in Figure 16 effect of the void, which is most important for short cracks initiated from void. The influence on the stress intensity factors for the shorter cracks is more pronounced with higher void ratios hv/rv.  Table 0. mm, hv/rv = 0.25 or 1).
The influence of the void on the asymmetry of the cycle R is the same as in the case without the void-the asymmetry in modes II and III is not affected by the crack propagation, the asymmetry in mode I increases slightly with longer cracks.
The results of the third version of the model that contains the void and also includes friction between the crack faces are plotted in Figure 16  The cycle of stress intensity factor KI has remained unchanged compared to the model with a void without friction. It is obvious, since the friction between crack faces cannot influence the tensile opening mode. The overall nature of the cycles of KII and KIII with The cycle of stress intensity factor K I has remained unchanged compared to the model with a void without friction. It is obvious, since the friction between crack faces cannot influence the tensile opening mode. The overall nature of the cycles of K II and K III with friction is also very similar to the previous model without friction. The friction affects mostly those values close to the positions of the closed crack-around the angle ρ = 90 • . The difference between the values is approximately 5-10%. This means that friction does not significantly influence the overall K max functions, which can be seen also from the plot in Figure 17 that illustrates the overall influence of the friction on K IImax and K IIImax values during the crack propagation. There is again a slight shift in the positions where K II and K III cycles reach their maxima and minima with the growing crack (observable in Figures A7-A10). The difference between the values is approximately 5-10%. This means that friction does not significantly influence the overall Kmax functions, which can be seen also from the plot in Figure 17 that illustrates the overall influence of the friction on KIImax and KIIImax values during the crack propagation. There is again a slight shift in the positions where KII and KIII cycles reach their maxima and minima with the growing crack (observable in Figures  A7-A10). It is important to remark here that heat is generated during friction, as the surfaces rub against each other. This effect is not considered in the described simulations. However, it may be significant, especially for cracks of longer lengths. Thermoplastic parts are generally not performing well in higher temperatures and even a slight increase in temperature can very negatively affect the resistance of the material against crack propagation. This is illustrated in the experimental study of mixed-mode I/III crack propagation in POM [31]. The work shows that the lifetime of the specimens decreases significantly with higher temperatures achieved as a result of friction in the vicinity of the crack. The high sensitivity of polymers to friction-generated heat is also reported in [32,40].

Conclusions
In the presented study, a parametrical model of a cylindrical bearing element with central crack made of the POM was developed. It was used to simulate and analyze the process of crack growth initiated from an internal defect, which is often the cause of failure of thermoplastic bearing elements.
The first version of the model contained just a circular flat discontinuity propagating from the middle of the cylinder towards its outer surface. It was observed that the gradual change in the orientation of the crack against the load during the rotation caused a combination of all three possible crack opening modes, characterized by stress intensity factors KI, KII and KIII, on the crack front. Originally, it was assumed that the opening mode characterized by the KI was the dominant mode, but it was found that the crack was actually closed with the load pressing the faces together during most of the cycle. Additionally, the maximum values of the shear mode stress intensity factors, KII and KIII, were higher than the maximum reached by KI.
A parametric study of the stress intensity factor values during crack propagation was carried out that resulted in the formation of parametric equations characterizing the course of the maximum values of stress intensity factors KImax, KIImax and KIIImax during the rotation. These equations contain the dimensions of the cylinders and can be used to quickly describe the crack tip stress situation in the cylinder upon crack propagation without having to carry out a lengthy FEM simulation.
The second version of the model was made by adding a spheroidal void. The void simulated a presence of a manufacturing defect. It turned out that the presence of the void It is important to remark here that heat is generated during friction, as the surfaces rub against each other. This effect is not considered in the described simulations. However, it may be significant, especially for cracks of longer lengths. Thermoplastic parts are generally not performing well in higher temperatures and even a slight increase in temperature can very negatively affect the resistance of the material against crack propagation. This is illustrated in the experimental study of mixed-mode I/III crack propagation in POM [31]. The work shows that the lifetime of the specimens decreases significantly with higher temperatures achieved as a result of friction in the vicinity of the crack. The high sensitivity of polymers to friction-generated heat is also reported in [32,40].

Conclusions
In the presented study, a parametrical model of a cylindrical bearing element with central crack made of the POM was developed. It was used to simulate and analyze the process of crack growth initiated from an internal defect, which is often the cause of failure of thermoplastic bearing elements.
The first version of the model contained just a circular flat discontinuity propagating from the middle of the cylinder towards its outer surface. It was observed that the gradual change in the orientation of the crack against the load during the rotation caused a combination of all three possible crack opening modes, characterized by stress intensity factors K I , K II and K III , on the crack front. Originally, it was assumed that the opening mode characterized by the K I was the dominant mode, but it was found that the crack was actually closed with the load pressing the faces together during most of the cycle. Additionally, the maximum values of the shear mode stress intensity factors, K II and K III , were higher than the maximum reached by K I .
A parametric study of the stress intensity factor values during crack propagation was carried out that resulted in the formation of parametric equations characterizing the course of the maximum values of stress intensity factors K Imax , K IImax and K IIImax during the rotation. These equations contain the dimensions of the cylinders and can be used to quickly describe the crack tip stress situation in the cylinder upon crack propagation without having to carry out a lengthy FEM simulation.
The second version of the model was made by adding a spheroidal void. The void simulated a presence of a manufacturing defect. It turned out that the presence of the void decreased the overall stiffness of the element and caused greater opening of the crack in mode I and, thus, higher K I in the position γ = 90 • (and −90 • ). Contrary to that, the K I 0 • and K II substantially decreased in the presence of the void. The K III remains almost unchanged even in cases with a rather large void. The common observation was that with the crack propagating further from the void, the void influence decreased up to a point where the stress intensity factor values were the same as for the flat crack.
The third version of the model included friction between the crack faces. It was found that the presence of friction did not substantially influence the values of stress intensity factor.
The potential of the results presented in this paper is in the possibility to use them for the calculation of lifetime estimations without complex numerical modelling. However, this requires reliable material data that characterizes the crack propagation rate under mixed-mode conditions.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
This appendix contains 2D plots of the course of stress intensity factors in the cylinder with a flat crack. The 2D plots are important for a better understanding of the process of changing the conditions of the crack propagation during the rolling of the bearing cylinder. However, a lot of 2D plots would be necessary to describe the whole process of rolling sufficiently. Only the most important situations are illustrated in the form of 2D plots in this part.
The unfilled squares in all the following plots are theoretical negative values of K I in the crack opening mode I, if the faces of cracks are pressed against each other. In reality, the stress intensity factor in such a situation is 0 and it cannot be negative. The values are plotted to illustrate the instances, in which the crack closes and later opens again.

Appendix B
This appendix contains 2D plots of the course of stress intensity factors in the cylinder with the void and friction.

Appendix B
This appendix contains 2D plots of the course of stress intensity factors in the cylinder with the void and friction.