Crack Propagation Analysis of Compression Loaded Rolling Elements

materials Article Crack Propagation Analysis of Compression Loaded Rolling Elements Pavol Dlhý1,2 , Jan Poduška 1,3,* , Michael Berer 4, Anja Gosch 5, Ondrej Slávik 1,3 , Luboš Náhlík1and Pavel Hutaˇr 1   Citation: Dlhý, P.; Poduška, J.; Berer, M.; Gosch, A.; Slávik, O.; Náhlík, L.; Hutaˇr, P. Crack Propagation Analysis of Compression Loaded Rolling Elements. Materials 2021,14, 2656. https://doi.org/10.3390/ma14102656 Academic Editor: Christian Motz Received: 15 April 2021 Accepted: 17 May 2021 Published: 19 May 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Institute of Physics of Materials, Czech Academy of Sciences, Žižkova 22, 616 00 Brno, Czech Republic; [email protected].cz (P.D.); [email protected] (O.S.); [email protected] (L.N.); [email protected] (P.H.) 2Central European Institute of Technology, Brno University of Technology, Purkyˇnova 656/123, 612 00 Brno, Czech Republic 3Faculty of Mechanical Engineering, Brno University of Technology, Technická2896/2, 616 69 Brno, Czech Republic 4 Polymer Competence Center Leoben GmbH, Roseggerstrasse 12, 8700 Leoben, Austria; [email protected] 5Material Science and Testing of Polymers, Montanuniversitaet Leoben, Otto Gloeckel-Straße 2, 8700 Leoben, Austria; [email protected] *Correspondence: [email protected] Abstract: 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. Keywords: fracture mechanics; crack propagation; stress intensity factor; thermoplastic material; bearing element; finite element method 1. 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–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 – 9 ], but it appears also in bearing elements such as balls and cylinders [ 10 – 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. Materials 2021,14, 2656. https://doi.org/10.3390/ma14102656 https://www.mdpi.com/journal/materials Materials 2021,14, 2656 2 of 24 Thermoplastic cylindrical bearing elements are typically manufactured by injection molding, which is a very common way well known for its high processing speed. However, internal defects occur in the cylinders as a result of significant shrinkage from the solidification of the molten polymer in the mold after the injection [ 3 ]. When a defected cylinder is installed in a bearing and starts rolling under load, a crack can initiate from the shrinkage defect and start propagating. This process eventually results in failure of the cylinder. Similar to the previously mentioned rolling contact fatigue process, the crack grows under mixed-mode conditions that are caused by the ever-changing orientation of the crack to the load acting during rolling. As the production defects are very difficult to eliminate [ 3 ], it is important to describe the process that leads to failure and thus enable reasonable lifetime estimation of these parts. The crack propagation in this case can be described by the linear elastic fracture mechanics (LEFM) parameters, e.g., the stress intensity factors. The relation between fatigue crack growth rate and stress intensity factor for polymer materials is commonly used [ 14 – 22 ]. Then, an estimation of the total lifetime can be performed if the crack growth rate is known from experiments and the considered crack is accurately described by fracture parameters—see details in [ 23 ]. This process is well established and used, e.g., for polyethylene (PE) pipes [ 23 – 26 ], or recently for polymer parts produced by additive manufacturing [ 27 , 28 ]. Approaches for lifetime estimation based on fracture mechanics are used also for the delamination in composites with polymer matrix [ 29 ] and plasticencapsulated microelectronics [ 30 ]. However, all the referenced studies considered mode I crack propagation only. Frequent presence of mixed-mode conditions in fracture mechanics problems concerning parts such as bearing elements or gears made of thermoplastics such as POM, PEEK and other is fueling recent interest in testing the crack growth conditions in these materials in mixed-mode conditions—see [ 31 ]. It was found that mixed-mode loading conditions lead to a clear lifetime reduction in the case of mode I/mode III loading. Additionally, during fatigue testing in mixed-mode conditions important temperature increase, caused by friction between crack flanks, was observed. This temperature increase depends on the loading frequency and can significantly deteriorate mechanical properties of the polymer. Therefore, any lifetime prediction that neglects the effect of the mixed-mode conditions leads to strongly non-conservative results. The detailed knowledge of the mode mixity and accurate estimation of the stress intensity factors during crack propagation is necessary for further lifetime prediction of the component. A thorough finite element simulation (FEM) of a real component where cracks propagate in mixed-mode conditions is also essential for the further practical application of the findings of previous research on the mixed-mode crack propagation in thermoplastics [31–33]. Detailed FEM simulation of the crack propagation in cylindrical bearing elements made of POM during rolling is described in this paper. This study is a direct follow-up to the work of Berer et al. [ 3 ], but it takes into account the presence of the mixed-mode conditions at the crack front. The stress intensity factors for the mode I (tensile crack opening), mode II (in-plane shear crack opening) and mode III (out-of-plane shear crack opening) were determined and their mutual relationship was assessed. General equations that allow the assessment of stress intensity factors for a wide range of dimensions of the cylinders were established. Different conditions of the crack propagation were considered in the simulation. Attention was paid to the shape of the initial defect and the friction between crack faces, as they come into contact in some positions of the rolling cylinder. 2. Numerical Models A parametrical numerical model was developed in order to simulate a growing crack inside a rolling cylinder under different conditions and to determine fracture mechanics 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). Materials 2021,14, 2656 3 of 24 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. 2.1. 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. Materials 2021, 14, x FOR PEER REVIEW 3 of 24 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. 2.1. 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. Figure 1. The situation of the bearing cylinder with an internal defect, compressed between two steel plates and rolling. 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— Figure 1. The situation of the bearing cylinder with an internal defect, compressed between two steel plates and rolling. 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, Materials 2021,14, 2656 4 of 24 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. Materials 2021, 14, x FOR PEER REVIEW 4 of 24 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. Figure 2. The mesh of the whole symmetrical FEM model of the bearing cylinder. The refined area in the vicinity of the crack tip is pictured in the detail view (note the special crack tip elements in the middle). 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 same basic set of 7 combinations of dimensions D × L as in Berer et al.’s previous work [3] was chosen to be modelled (dimensions are in mm): 3 × 3, 4 × 4, 5 × 5, 6 × 6, 6 × 3, 6 × 9, 6 × 12. Nine more combinations were added to make the set more robust for the intended parametric study: 3 × 4, 3 × 5, 3 × 6, 3 × 9, 3 × 12, 4 × 5, 4 × 6, 5 × 4 and 6 × 4. 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 Figure 2. The mesh of the whole symmetrical FEM model of the bearing cylinder. The refined area in the vicinity of the crack tip is pictured in the detail view (note the special crack tip elements in the middle). 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 same basic set of 7 combinations of dimensions D × Las in Berer et al.’s previous work [ 3 ] was chosen to be modelled (dimensions are in mm): 3 × 3, 4 × 4, 5 × 5, 6 × 6, 6×3 , 6 × 9, 6 × 12. Nine more combinations were added to make the set more robust for the intended parametric study: 3 × 4, 3 × 5, 3 × 6, 3 × 9, 3 × 12, 4 × 5, 4 × 6, 5 × 4 and 6×4. The growth of the crack inside the cylinder was simulated by modelling the cylinder with the crack length agoing 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 ais 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 1for 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.8 ◦ , 7.5 ◦ , 11.2 ◦ and 22.5 ◦ as well as 157.5 ◦ , 168.8 ◦ , 172.5 ◦ and 176.2. Materials 2021,14, 2656 5 of 24 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. Materials 2021, 14, x FOR PEER REVIEW 5 of 24 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. Figure 3. Schematic illustration of the changing orientation of the crack during rolling and the angle ρ that describes the rolling position. 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]. 2.2. 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. Figure 3. Schematic illustration of the changing orientation of the crack during rolling and the angle ρ that describes the rolling position. 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]. 2.2. 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 × Lwere 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 4for 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 aconsidered 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. Materials 2021,14, 2656 6 of 24 Materials 2021, 14, x FOR PEER REVIEW 6 of 24 Figure 4. Model of a bearing cylinder with an initial spheroidal void in the middle. 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–40]. However, the friction coefficient can significantly vary for different blends of the same polymer. 3. 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 K III , 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-eighthtype 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. Figure 4. Model of a bearing cylinder with an initial spheroidal void in the middle. 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–40] . However, the friction coefficient can significantly vary for different blends of the same polymer. 3. 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 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 5represents 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. Materials 2021,14, 2656 7 of 24 Materials 2021, 14, x FOR PEER REVIEW 7 of 24 Figure 5. Comparison of K I values from Berer et al. [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). 3.1. 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 Figure 5. Comparison of K I values from Berer et al. [ 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). 3.1. 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 aand 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 Materials 2021,14, 2656 8 of 24 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. Materials 2021, 14, x FOR PEER REVIEW 8 of 24 be observed in Figures A1–A6 in Appendix A), which suggests that the crack might not propagate exactly as a regular circle in reality. Figure 6. 3D plot of K I as a function of both, the circumferential position on the crack front γ and the rolling orientation of the crack ρ, crack length a = 1.25 mm. The upper surface reflects the situation with contact defined between crack faces (highlighted by solid red line). The lower surface are theoretical negative values of K I acquired from a model where no contact was defined between crack faces (highlighted by dashed red line). 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 Figures 7 and 8. Figure 7. 3D plot of K II as a function of both, the circumferential position on the crack front γ and the rolling orientation of the crack ρ, crack length a = 1.25 mm. Figure 6. 3D plot of K I as a function of both, the circumferential position on the crack front γ and the rolling orientation of the crack ρ , crack length a= 1.25 mm. The upper surface reflects the situation with contact defined between crack faces (highlighted by solid red line). The lower surface are theoretical negative values of K I acquired from a model where no contact was defined between crack faces (highlighted by dashed red line). 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 Figures 7and 8. Materials 2021, 14, x FOR PEER REVIEW 8 of 24 be observed in Figures A1–A6 in Appendix A), which suggests that the crack might not propagate exactly as a regular circle in reality. Figure 6. 3D plot of K I as a function of both, the circumferential position on the crack front γ and the rolling orientation of the crack ρ, crack length a = 1.25 mm. The upper surface reflects the situation with contact defined between crack faces (highlighted by solid red line). The lower surface are theoretical negative values of K I acquired from a model where no contact was defined between crack faces (highlighted by dashed red line). 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 Figures 7 and 8. Figure 7. 3D plot of K II as a function of both, the circumferential position on the crack front γ and the rolling orientation of the crack ρ, crack length a = 1.25 mm. Figure 7. 3D plot of K II as a function of both, the circumferential position on the crack front γ and the rolling orientation of the crack ρ, crack length a= 1.25 mm. Materials 2021,14, 2656 9 of 24 Materials 2021, 14, x FOR PEER REVIEW 9 of 24 Figure 8. 3D plot of K III as a function of both, the circumferential position on the crack front γ and the rolling orientation of the crack ρ, crack length a = 1.25 mm. 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 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 Figure 8. 3D plot of K III as a function of both, the circumferential position on the crack front γ and the rolling orientation of the crack ρ, crack length a= 1.25 mm. The crack tip mixed-mode loading has a completely out-of-phase character. Again, refer to Figures A1–A6 in Appendix Afor 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 Kare 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 ais 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 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 Materials 2021,14, 2656 16 of 24 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). Materials 2021, 14, x FOR PEER REVIEW 16 of 24 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). Figure 17. K IImax and K IIImax depending on the normalized crack length; comparison of model with void without friction and model with void with friction; void ratio h v /r v = 0.25, r v = 0.75 mm. 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]. 4. 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 Figure 17. K IImax and K IIImax depending on the normalized crack length; comparison of model with void without friction and model with void with friction; void ratio hv/rv= 0.25, rv= 0.75 mm. 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]. 4. 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 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 K Imax ,K IImax and K IIImax during the rotation. These equations contain the dimensions of the cylinders and can be used Materials 2021,14, 2656 17 of 24 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. Author Contributions: Conceptualization, P.D., J.P. and M.B.; methodology, P.D.; formal analysis, P.D. and O.S.; investigation, P.D.; resources, M.B. and P.H.; data curation, P.D. and J.P.; writing— original draft preparation, J.P. and P.D.; writing—review and editing, M.B., A.G. and P.H.; supervision, P.H. and L.N.; project administration, M.B.; funding acquisition, M.B. and L.N. All authors have read and agreed to the published version of the manuscript. Funding: The research has obtained support from the Polymer Competence Center Leoben GmbH (PCCL, Austria) within the framework of the COMET-program of the Federal Ministry for Climate Action, Environment, Energy, Mobility, Innovation and Technology and the Federal Ministry for Digital and Economic Affairs with contributions by the Department of Polymer Engineering and Science, Montanuniversitaet Leoben (Austria). The PCCL is funded by the Austrian Government and the State Governments of Styria, Lower Austria and Upper Austria. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. 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. Materials 2021,14, 2656 18 of 24 Materials 2021, 14, x FOR PEER REVIEW 18 of 24 Figure A1. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm, load F = 350 N, a = 0.25 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A2. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 0.75 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A1. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm, load F= 350 N, a= 0.25 mm, position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Materials 2021, 14, x FOR PEER REVIEW 18 of 24 Figure A1. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm, load F = 350 N, a = 0.25 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A2. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 0.75 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A2. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 0.75 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Materials 2021,14, 2656 19 of 24 Materials 2021, 14, x FOR PEER REVIEW 19 of 24 Figure A3. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 0.85 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A4. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.00 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A3. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 0.85 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Materials 2021, 14, x FOR PEER REVIEW 19 of 24 Figure A3. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 0.85 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A4. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.00 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A4. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 1.00 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Materials 2021,14, 2656 20 of 24 Materials 2021, 14, x FOR PEER REVIEW 20 of 24 Figure A5. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.25 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A6. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.75 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Appendix B This appendix contains 2D plots of the course of stress intensity factors in the cylinder with the void and friction. Figure A5. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 1.25 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Materials 2021, 14, x FOR PEER REVIEW 20 of 24 Figure A5. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.25 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Figure A6. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.75 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Appendix B This appendix contains 2D plots of the course of stress intensity factors in the cylinder with the void and friction. Figure A6. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 1.75 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Appendix B This appendix contains 2D plots of the course of stress intensity factors in the cylinder with the void and friction. Materials 2021,14, 2656 21 of 24 Materials 2021, 14, x FOR PEER REVIEW 21 of 24 Figure A7. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 0.85 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Model including friction. Figure A8. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.00 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Model including friction. Figure A7. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 0.85 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Model including friction. Materials 2021, 14, x FOR PEER REVIEW 21 of 24 Figure A7. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 0.85 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Model including friction. Figure A8. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.00 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Model including friction. Figure A8. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 1.00 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Model including friction. Materials 2021,14, 2656 22 of 24 Materials 2021, 14, x FOR PEER REVIEW 22 of 24 Figure A9. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.25 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Model including friction. Figure A10. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.75 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Model including friction. References 1. Avalle, M.; Romanello, E. Tribological characterization of modified polymeric blends. Procedia Struct. Integr. 2018, 8, 239–255, doi:10.1016/j.prostr.2017.12.026. 2. Harrass, M.; Friedrich, K.; Almajid, A. Tribological behavior of selected engineering polymers under rolling contact. Tribol. Int. 2010, 43, 635–646, doi:10.1016/j.triboint.2009.10.003. 3. Berer, M.; Mitev, I.; Pinter, G. Finite element study of mode I crack opening effects in compression-loaded cracked cylinders. Eng. Fract. Mech. 2017, 175, 1–14, doi:10.1016/j.engfracmech.2017.03.008. 4. Berer, M.; Major, Z. Characterization of the global deformation behaviour of engineering plastics rolls. Int. J. Mech. Mater. Des. 2010, 6, 1–9, doi:10.1007/s10999-010-9111-9. 5. Berer, M.; Major, Z. Characterisation of the Local Deformation Behaviour of Engineering Plastics Rolls. Strain 2011, 48, 225–234, doi:10.1111/j.1475-1305.2011.00816.x. 6. Bonniot, T.; Doquet, V.; Mai, S.H. Mixed mode II and III fatigue crack growth in a rail steel. Int. J. Fatigue 2018, 115, 42–52, doi:10.1016/j.ijfatigue.2018.01.010. Figure A9. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 1.25 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Model including friction. Materials 2021, 14, x FOR PEER REVIEW 22 of 24 Figure A9. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.25 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Model including friction. Figure A10. Course of stress intensity factors during rotation in the cylinder D × L = 6 × 6 mm 2 , load F = 350 N, a = 1.75 mm, position on the crack front; (a) γ = 0°; (b) γ = −90°. Model including friction. References 1. Avalle, M.; Romanello, E. Tribological characterization of modified polymeric blends. Procedia Struct. Integr. 2018, 8, 239–255, doi:10.1016/j.prostr.2017.12.026. 2. Harrass, M.; Friedrich, K.; Almajid, A. Tribological behavior of selected engineering polymers under rolling contact. Tribol. Int. 2010, 43, 635–646, doi:10.1016/j.triboint.2009.10.003. 3. Berer, M.; Mitev, I.; Pinter, G. Finite element study of mode I crack opening effects in compression-loaded cracked cylinders. Eng. Fract. Mech. 2017, 175, 1–14, doi:10.1016/j.engfracmech.2017.03.008. 4. Berer, M.; Major, Z. Characterization of the global deformation behaviour of engineering plastics rolls. Int. J. Mech. Mater. Des. 2010, 6, 1–9, doi:10.1007/s10999-010-9111-9. 5. Berer, M.; Major, Z. Characterisation of the Local Deformation Behaviour of Engineering Plastics Rolls. Strain 2011, 48, 225–234, doi:10.1111/j.1475-1305.2011.00816.x. 6. Bonniot, T.; Doquet, V.; Mai, S.H. Mixed mode II and III fatigue crack growth in a rail steel. Int. J. Fatigue 2018, 115, 42–52, doi:10.1016/j.ijfatigue.2018.01.010. Figure A10. Course of stress intensity factors during rotation in the cylinder D × L= 6 × 6 mm 2 , load F= 350 N, a= 1.75 mm , position on the crack front; (a)γ= 0◦; (b)γ=−90◦. Model including friction. References 1. Avalle, M.; Romanello, E. Tribological characterization of modified polymeric blends. Procedia Struct. Integr. 2018 ,8, 239–255. [CrossRef] 2. Harrass, M.; Friedrich, K.; Almajid, A. Tribological behavior of selected engineering polymers under rolling contact. Tribol. Int. 2010,43, 635–646. [CrossRef] 3. Berer, M.; Mitev, I.; Pinter, G. Finite element study of mode I crack opening effects in compression-loaded cracked cylinders. Eng. Fract. Mech. 2017,175, 1–14. [CrossRef] 4. Berer, M.; Major, Z. Characterization of the global deformation behaviour of engineering plastics rolls. Int. J. Mech. Mater. Des. 2010,6, 1–9. [CrossRef] 5. Berer, M.; Major, Z. Characterisation of the Local Deformation Behaviour of Engineering Plastics Rolls. Strain 2011 ,48, 225–234. [CrossRef] Materials 2021,14, 2656 23 of 24 6. Bonniot, T.; Doquet, V.; Mai, S.H. Mixed mode II and III fatigue crack growth in a rail steel. Int. J. Fatigue 2018 ,115, 42–52. [CrossRef] 7. Bold, P.; Brown, M.; Allen, R. Shear mode crack growth and rolling contact fatigue. Wear 1991,144, 307–317. [CrossRef] 8. Liu, H. Material Modelling for Structural Analysis of Polyethylene. Master’s Thesis, University of Waterloo, Waterloo, ON, Canada, 2007. 9. Liu, Y.; Stratman, B.; Mahadevan, S. Fatigue crack initiation life prediction of railroad wheels. Int. J. Fatigue 2006 ,28, 747–756. [CrossRef] 10. Avanzini, A.; Donzella, G.; Mazzù, A.; Petrogalli, C. Wear and rolling contact fatigue of PEEK and PEEK composites. Tribol. Int. 2013,57, 22–30. [CrossRef] 11. Kadin, Y.; Rychahivskyy, A. Modeling of surface cracks in rolling contact. Mater. Sci. Eng. A 2012,541, 143–151. [CrossRef] 12. Berer, M.; Major, Z.; Pinter, G. Elevated pitting wear of injection molded polyetheretherketone (PEEK) rolls. Wear 2013 ,297, 1052–1063. [CrossRef] 13. Koike, H.; Kida, K.; Honda, T.; Mizobe, K.; Oyama, S.; Rozwadowska, J.; Kashima, Y.; Kanemasu, K. Observation of Crack Propagation in PEEK Polymer Bearings under Water-Lubricated Conditions. Adv. Mater. Res. 2012,566, 109–114. [CrossRef] 14. Berer, M.; Pinter, G. Determination of crack growth kinetics in non-reinforced semi-crystalline thermoplastics using the linear elastic fracture mechanics (LEFM) approach. Polym. Test. 2013,32, 870–879. [CrossRef] 15. Berer, M.; Pinter, G.; Feuchter, M. Fracture mechanical analysis of two commercial polyoxymethylene homopolymer resins. J. Appl. Polym. Sci. 2014,131, 1–15. [CrossRef] 16. Furmanski, J.; Pruitt, L.A. Peak stress intensity dictates fatigue crack propagation in UHMWPE. Polymer 2007 ,48, 3512–3519. [CrossRef] 17. Favier, V.; Giroud, T.; Strijko, E.; Hiver, J.; G’Sell, C.; Hellinckx, S.; Goldberg, A. Slow crack propagation in polyethylene under fatigue at controlled stress intensity. Polymer 2002,43, 1375–1382. [CrossRef] 18. Pruitt, L.; Sreekanth, P.R.; Badgayan, N.; Sahoo, S. Fatigue of Polymers. In Reference Module in Materials Science and Materials Engineering; Elsevier BV: Amsterdam, The Netherlands, 2017; pp. 1–15. 19. Harris, J.S.; Ward, I.M. Fatigue-crack propagation in vinyl urethane polymers. J. Mater. Sci. 1973,8, 1655–1665. [CrossRef] 20. Kanters, M.J.; Kurokawa, T.; Govaert, L.E. Competition between plasticity-controlled and crack-growth controlled failure in static and cyclic fatigue of thermoplastic polymer systems. Polym. Test. 2016,50, 101–110. [CrossRef] 21. Kanters, M.J.; Stolk, J.; Govaert, L.E. Direct comparison of the compliance method with optical tracking of fatigue crack propagation in polymers. Polym. Test. 2015,46, 98–107. [CrossRef] 22. Arbeiter, F.; Spoerk, M.; Wiener, J.; Gosch, A.; Pinter, G. Fracture mechanical characterization and lifetime estimation of nearhomogeneous components produced by fused filament fabrication. Polym. Test. 2018,66, 105–113. [CrossRef] 23. Hutaˇr, P.; Ševˇcík, M.; Náhlík, L.; Pinter, G.; Frank, A.; Mitev, I. A numerical methodology for lifetime estimation of HDPE pressure pipes. Eng. Fract. Mech. 2011,78, 3049–3058. [CrossRef] 24. Frank, A.; Arbeiter, F.J.; Berger, I.J.; Hutaˇr, P.; Náhlík, L.; Pinter, G. Fracture Mechanics Lifetime Prediction of Polyethylene Pipes. J. Pipeline Syst. Eng. Pr. 2019,10, 04018030. [CrossRef] 25. Frank, A.; Hutaˇr, P.; Pinter, G. Numerical Assessment of PE 80 and PE 100 Pipe Lifetime Based on Paris-Erdogan Equation. Macromol. Symp. 2012,311, 112–121. [CrossRef] 26. Benhamena, A.; Bouiadjra, B.B.; Amrouche, A.; Mesmacque, G.; Benseddiq, N.; Benguediab, M. Three finite element analysis of semi-elliptical crack in high density poly-ethylene pipe subjected to internal pressure. Mater. Des. 2010 ,31, 3038–3043. [CrossRef] 27. Arbeiter, F.; Trávníˇcek, L.; Petersmann, S.; Dlhý, P.; Spoerk, M.; Pinter, G.; Hutaˇr, P. Damage tolerance-based methodology for fatigue lifetime estimation of a structural component produced by material extrusion-based additive manufacturing. Addit. Manuf. 2020,36, 101730. [CrossRef] 28. Puigoriol-Forcada, J.M.; Alsina, A.; Salazar-Martín, A.G.; Gomez-Gras, G.; Pérez, M.A. Flexural fatigue properties of polycarbonate fused-deposition modelling specimens. Mater. Des. 2018,155, 414–421. [CrossRef] 29. Jones, R.; Kinloch, A.; Michopoulos, J.; Brunner, A.; Phan, N. Delamination growth in polymer-matrix fibre composites and the use of fracture mechanics data for material characterisation and life prediction. Compos. Struct. 2017,180, 316–333. [CrossRef] 30. Evans, J.W.; Sinha, K. Applications of fracture mechanics to quantitative accelerated life testing of plastic encapsulated microelectronics. Microelectron. Reliab. 2018,80, 317–327. [CrossRef] 31. Gosch, A.; Berer, M.; Hutaˇr, P.; Slávik, O.; Vojtek, T.; Arbeiter, F.J.; Pinter, G. Mixed Mode I/III fatigue fracture characterization of Polyoxymethylene. Int. J. Fatigue 2020,130, 105269. [CrossRef] 32. Gosch, A.; Arbeiter, F.J.; Berer, M.; Vojtek, T.; Hutaˇr, P.; Pinter, G. Fatigue characterization of polyethylene under mixed mode I/III conditions. Int. J. Fatigue 2021,145, 106084. [CrossRef] 33. Schrader, P.; Gosch, A.; Berer, M.; Marzi, S. Fracture of Thin-Walled Polyoxymethylene Bulk Specimens in Modes I and III. Materials 2020,13, 5096. [CrossRef] 34. ANSYS Help Release 2020 R2; ANSYS Inc.: Canonsburg, PA, USA, 2020. 35. Shih, C.F.; Moran, B.; Nakamura, T. Energy release rate along a three-dimensional crack front in a thermally stressed body. Int. J. Fract. 1986,30, 79–102. [CrossRef] 36. Ingraffea, A.R.; Manu, C. Stress-intensity factor computation in three dimensions with quarter-point elements. Int. J. Numer. Methods Eng. 1980,15, 1427–1445. [CrossRef] Materials 2021,14, 2656 24 of 24 37. Lind, J.; Lindholm, P.; Qin, J.; Kassman, R. Friction and wear studies of some peek materials. Tribologia 2015,33, 20–28. 38. Chaudri, A.M.; Suvanto, M.; Pakkanen, T.T. Non-lubricated friction of polybutylene terephthalate (PBT) sliding against polyoxymethylene (POM). Wear 2015,342–343, 189–197. [CrossRef] 39. Chen, J.; Cao, Y.; Li, H. Investigation of the friction and wear behaviors of polyoxymethylene/linear low-density polyethylene/ethylene-acrylic-acid blends. Wear 2006,260, 1342–1348. [CrossRef] 40. Hoskins, T.; Dearn, K.; Chen, Y.; Kukureka, S. The wear of PEEK in rolling–sliding contact—Simulation of polymer gear applications. Wear 2014,309, 35–42. [CrossRef] 41. Chen, Y.T.; Liu, K.X. Crack propagation in viscoplastic polymers: Heat generation in near-tip zone and viscoplastic cohesive model. Appl. Phys. Lett. 2015,106, 061908. [CrossRef]