A Computational Model for Analysing the Wear Behaviour of Polymer Gears

1. Introduction

Polymer gears are used widely in many engineering applications, such as office appliances, mechatronic devices, household facilities, computer and laboratory equipment, medical instruments, etc. [1,2,3,4,5,6,7]. In these applications, the contact of two mechanical elements (i.e. gear flanks) is concentrated on a small contact area, resulting in higher contact stresses, which often exceed the yield stress of the material. Due to repeated contact loading during the exploitation, the local cyclic plastic deformations may lead to surface damage, usually termed surface wear [8,9,10,11].

In recent years, many wear models relating to rolling/sliding contacting mechanical elements were proposed. However, Archard’s model [12] is widely used in engineering applications. The model requires information about the coefficient of wear, relative sliding distance and load and has been applied to gear wear prediction by many researchers, who have derived equations that connect the wear volume or wear depth with the operating conditions and material properties of the gear tooth. Flodin et al. [13,14] derived a load distribution equation to determine the transferred load between the spur and helical gear pairs, using Herztzian contact theory to calculate the contact stress and the single point tracking method for relative sliding distance. Later, this model was extended to analyse the wear behaviour of double helical gears [15], worm gears [16], internal gears [17], hypoid bevel gears [18,19,20] and planetary gears [21,22].

Discretising the Archard model allows for a more detailed analysis of the wear behaviour of materials, especially where wear varies across the material surface. Põdra et al. [23] used the finite element method for simulating sliding wear between the rotation disk and pin. For sliding-rolling applications, such as gears, Lundvall et al. [24,25] derived time increment equations to predict wear using the nonsmooth Newton method. Employing the finite element method opens more possibilities to improve accuracy, which is why Bajpai et al. [26] combined the finite element method and the Archard model, using the iterative method to consider contact changes due to wear. Because mechanical elements like gears, bearings etc., are dynamically loaded components, it is also useful to consider time-varying sliding distance and time-varying contact load in the wear analysis to make results more accurate [27,28,29,30,31,32].

In recent years, polymer gears have become more widely used due to their good properties such as lubrication-free operation, corrosion resistance, lightweight, arbitrary geometry and low cost of manufacture. They are mostly found in gears within small appliances, where they can be combined with a variety of other materials [33]. As they are under cyclic load, they are dynamically loaded components. As the gears roll and slide between the tooth flanks of the gears, erosion and wear occur on the surfaces over time, which can be determined with experiments, analytically or numerically. The analytical calculation of polymer gears can be found in the VDI 2736 standard (Association of German Engineers) [34]. This covers the basic calculation and the presentation of the properties of polymer gears. Among those things, the wear depth of the gears can be calculated using this standard. Although this standard can be used to obtain certain values, it is limited to constant operating conditions, which, unlike gears from metallic materials, can change rapidly due to their poorer mechanical properties. As mentioned, the finite element method provides more accurate results and allows rapid change of material properties and operating conditions. Lin et al. [35] investigated the interaction between contact loads and tooth wear in POM (poly-oxy-methylene) and nylon tooth pairs, using the finite element method to analyse tooth stiffness and a quasi-linear iterative method to calculate wear depth. The result shows that peak load increases as the number of cycles reaches a threshold limit.

In the framework of the presented study, a computational model was built using open-source finite element numerical calculation software PrePoMax to predict wear behaviour based on Archard’s wear model. This allows the calculation of the wear depth in each loading cycle with constant mesh updating. The computational results that are obtained are compared with the analytical results according to the VDI 2736 guidelines.

2. Materials and Methods

2.1. Computational modelling

In this study, the PrePoMax open-source pre- and post-processor for the finite element method (FEM) [36] was used to analyse the wear behaviour of rolling/sliding contacting mechanical elements (i.e. mating gear flanks). The solver used by PrePoMax is the open-source CalculiX FEM solver [37] that does not support the wear analysis by itself, so the wear analysis was implemented into the PrePoMax. The linear Archard wear model [38] for predicting the volume of worn material was chosen for this purpose. The model was originally developed based on experimental results of cyclic wear of metal specimens under dry contact conditions. The model is represented by the following equation:

$$V=k·\frac{F}{H}·s$$

(1)

where $V$ is the volume of worn material, $k$ is the dimensionless wear coefficient, $F$ is the normal contact force, $H$ is the surface hardness in contact, and $s$ represents the sliding distance. To account for changes in the geometry of the contact surfaces during wear cycles, the aforementioned model, which neglects the geometry change, can be divided into multiple steps when used in numerical analyses using the FEM. In a multi-step numerical procedure, where the depth of wear needs to be calculated for each wear cycle [39], Eq. (1) can be expressed as follows:

$$h_i=\int \frac{k}{H}·p_i·\mathrm{d}s$$

(2)

where $h_i$ represents the depth of wear in the i-th wear cycle, $p_i$ denotes the normal contact pressure, and $\mathrm{d}s$ represents the differential sliding distance. The total wear depth $h$ over $n$ wear cycles is then computed by the Eq. (3), summing up the wear depths of individual cycles. Such an approach enables the wear analysis to be independent of the model geometry, boundary conditions, and loading.

$$h={\sum }_{i=1}^n{h}_i$$

(3)

The multi-step method implemented in PrePoMax is based on the following assumptions: (i) the entire wear process can be divided into multiple repetitive wear cycles, (ii) wear occurs only on the softer surface in the contact pair (slave surface), (iii) the change in contact surface geometry during one wear cycle is negligible and (iv) the dimensionless wear coefficient is constant.

For each wear cycle, a finite element model specific to that cycle needs to be prepared. The model must accurately capture the evolving contact conditions during the cycle, requiring a nonlinear contact analysis solved in multiple increments. For an accurate result, the number of increments must be large enough to result in a smooth movement of the contact pressure over the finite element mesh of the contact surfaces.

The wear calculation for the cycle begins after the nonlinear contact analysis is completed. All wear parameters are computed only at those slave surface finite element nodes that were in active contact with the master surface during the wear cycle. The normal contact pressure $p_{i,j}$ at the node $j$ and the nodal relative sliding distance $∆s_{i,j}$ are determined from the resulting contact fields. The wear depth at node $h_{i,j}$ is then computed using the expression:

$$h_{i,j}=\frac{k_j}{H_j}·p_{i,j}·∆s_{i,j}$$

(4)

where the material properties kj and Hj are taken at node j. Processing of each wear cycle concludes with the calculation of the change in the contact surface geometry due to wear, determined by the wear displacements at the nodes calculated from the wear depth. The node wear displacements due to wear in each cycle ${\stackrel{⃑}{h}}_{i,j}$ are calculated by a scalar product:

$${\stackrel{⃑}{h}}_{i,j}=-h_{i,j}·{\stackrel{⃑}{n}}_j$$

(5)

where where ${\stackrel{⃑}{n}}_j$ represents the node’s outward normal vector. The node’s outward normal vector is determined as an average value of the normal vectors of the FEs sharing the node $j$ (Figure 1). At the start of the next wear cycle, the wear displacements computed in the previous cycle are used to update the mesh. Thus, the simulation of the new wear cycle is carried out on a mesh where the wear from the previous cycles is considered

Applying the wear displacements to the mesh changes only the position of the surface nodes and thus deforms the shape of the finite elements, making all surface elements on the wear surface a little thinner. If the accumulated wear displacement is large enough, this might lead to unusable finite element shapes. This is especially problematic when a dense mesh is used on the contact surfaces to obtain an accurate contact pressure field since even a small wear displacement will deform a small finite element beyond being usable. To overcome this problem, local or global remeshing of the domain can be applied, or a method for mesh redistribution can be used. In PrePoMax software, a boundary displacement method (BDM) [40] for mesh redistribution was introduced. This is a two-step procedure where the surface nodes are first moved, and then new positions of the internal nodes are computed. From many different approaches to the BDM, the positions of the internal nodes were determined based on the solution of the linear elastic solid body problem. This was the most convenient approach since it can be solved with the FEM method on the existing mesh. The method was integrated into PrePoMax by an additional BDM simulation step that is added after each wear cycle. In this simulation step, the mesh from the previous wear step is taken, uniform elastic properties are assigned to all elements, and a prescribed boundary condition is applied to all nodes on the surface. In the nodes where the wear occurs, the value of the prescribed displacement vector is equal to the wear displacement ${\stackrel{⃑}{h}}_{i,j}$ and equal to $\stackrel{⃑}{0}$ for all other nodes. After the BDM simulation step completes, the computed nodal displacements for all nodes from the BDM step are used to deform the mesh before the next wear step is started.

Running multiple wear steps one after the other and deforming the mesh after each step can lead to the occurrence of high-frequency oscillations (roughness) of the mesh (see Figure 1), which results in an irregular contact pressure field. This causes an irregular wear displacement field, which additionally enhances the roughness of the contact surfaces and renders the results of the wear analysis unusable. A mesh smoothing algorithm on the wear displacement field was introduced to remove the high-frequency oscillations from the mesh. The smoothing algorithm is based on the Laplacian mesh smoothing procedure [41], where the new position of the node ${\stackrel{⃑}{x}}_j$ is computed from the positions of its neighbouring nodes ${\stackrel{⃑}{x}}_k$ as:

$${\stackrel{⃑}{x}}_j=\frac{1}{m}{\sum }_{k=1}^m{\stackrel{⃑}{x}}_k$$

(6)

where $m$ equals the number of node neighbours. Mesh smoothing is applied after the wear displacement field is computed and before the BDM step is prepared. The effect of the smoothing step can be increased by running the smoothing algorithm multiple times, which can be defined by the user. All these features make wear analyses applicable and robust for a large variety of problems, and the user can prepare the wear analysis without additional coding using only the built-in user interface. This makes the implementation usable for many engineers since all these features are already available in the latest release version of the PrePoMax v1.4.0, available online [42].

2.2. Geometrical and material parameters of the analysed gear pair

The developed computational model has been evaluated on the spur gear pair, where the pinion made of POM has meshed with a support gear made of steel. The basic geometrical and material parameters of the analysed gear pair are presented in Table 1.

2.3. FE model of gear pair

The geometry of the gear pair is shown in Figure 2, where the pinion is shown in grey, and the support gear is shown in green. Here, only one tooth of the gear and three teeth of the pinion are presented. It was due to the smaller stiffness of the pinion material (POM) if compared to the gear material (steel). The analysis of one load cycle was split into three simulation steps, as shown in Figure 3: static step (contact initialisation), wear step (computation of wear) and BDM step (mesh redistribution). In the simplified model presented in Figure 2, the reference points RP1 and RP2 have been defined in both axes of the gear pair. Here, RP1 was connected rigidly to the gear shaft, with the surfaces of the gear marked in red. Additionally, at this reference point, a torque of T = 1.6 Nm was applied acting to the pinion tooth, and a boundary condition was imposed, which allowed the pinion to move only about its axis. The point RP2 was connected rigidly to the pinion surfaces marked in blue. At this reference point, a boundary condition was prescribed first in the static step, which prevented displacement and rotation in all directions of the coordinate system. As sliding was required for the wear calculation procedure, the boundary condition in RP1 in the second step (Wear step) required a rotation around the pinion axis of ρ = 0.285 rad. The angle of rotation was prescribed in such a way that the gear only reached the position where the gear started to return to the already worn surface.

The computational model implemented into the PrePoMax software (see Section 2.1) was then used to analyse the wear behaviour of mating gear flanks. The numerical simulation of gear meshing consisted of three separate simulation steps, which followed one after the other. This simulation then repeats to account for the selected number of load cycles. Simulating each load cycle in multi-thousand load cycle analysis would require too much computational effort. At the same time, during one load cycle, the geometry of the model does not change significantly; thus, multiple load cycles can be merged into one simulation cycle that updates the mesh.

Archard’s model implemented in the PrePoMax software is limited to a 3D wear calculation, as it calculates the volume of wear. However, as the pinion and gear analysed in this study do not vary in thickness, it would be reasonable to simplify the model to 2D, but this was not possible. Therefore, a reduced thickness of 0.25 mm was assumed to reduce the computational time, and consequently, the load (torque) has also been reduced proportionally. Furthermore, an additional boundary condition has been added to both gears to prevent the gears from moving in the thickness direction. To best approximate the 2D model, a finite element mesh with prismatic finite elements was created using the boundary layer function, as shown in Figure 4. The overall model was meshed with a global mesh size of 1 mm. However, for better accuracy, the model around the contact was meshed with a local mesh size of 0.05 mm, which resulted in 3300 parabolic prismatic finite elements. The model used a linear overclosure-pressure contact relation with stiffness K = 105 N/mm³ and coefficient of friction μ = 0.25 mm.

Two computational models have been considered in the numerical simulation. In the first model (Model 1), a wear calculation after 10,000 load cycles was determined. First, for the purpose of comparison with VDI 2736, the 10,000 load cycles were calculated in one simulation cycle, which means that the geometry was not updated during the calculation. The purpose of this computational model was to compare the results with the VDI 2736 calculation, which does not consider the geometry change. In the second model (Model 2), 10,000 load cycles were used again. However, if compared to the first model, the number of cycles was divided into ten simulation cycles, meaning that after every 1,000 load cycles, the geometry was updated. The wear results were thus transferred directly to the finite element mesh, and the wear was recalculated on the updated geometry.

2.4. Analytical approach according to VDI 2736

Following the calculation according to VDI 2736, a wear behaviour can be obtained for the tooth flank profile of the polymer gear. The calculation is position-dependent, as the wear is different at each point due to the involute shape of the tooth. Here, the wear starts to decrease from the tip of the tooth up to the pitch point. Due to the involute shape of the tooth, sliding dominates at the starting engagement point, then the tooth rolls at the pitch circle, and later, towards the root of the tooth, sliding increases gradually again, and rolling decreases. Therefore, in the pitch point region, it can be observed that the wear is theoretically zero. Considering this fact, the local wear of the gear flank can be obtained as follows:

$$W_{local}=\frac{F_n}{b_w}\cdot N_L\cdot \zeta \cdot k_w$$

(7)

where b_w is a common face width, N_L is the number of load cycles, ζ is the specific local sliding of gear flanks, k_w is the wear coefficient, and F_n is the normal force acting on the tooth flank

$$F_n=\frac{2\cdot T_d}{d_w\cdot \cos {\alpha }_t}$$

(8)

where T_d is the nominal torque, d_w is the pitch circle diameter, and α_t is the pressure angle in the transverse section. It follows from Eq. (8) that the normal force F_n is a constant throughout the engagement line of gear flanks. However, its tangential and radial components change, resulting in a combination of sliding and rolling of gear flanks. When sliding is involved, a specific local sliding ζ appears, which should be considered when determining the local wear using Eg. (7). Figure 5 shows the sliding conditions in an arbitrary mating point Y before pitch point C. The circular velocities in an arbitrary mating point Y (v_y1= r_y1⋅ω₁ and v_y2= r_y2⋅ω₂), which are acting perpendicular to the lines O₁Y and O₂Y, can be distributed in the directions of common normal (v_ny1 and v_ny2) and common tangent (v_ty1 and v_ty2). Here, the normal velocity components must be the same, i.e. v_ny1 = v_ny2. Furthermore, the tangential components v_ty1 and v_ty2 can be obtained from the triangles O₁T₁Y and O₂T₂Y and the hatched triangles.

The sliding velocities of gear flanks in an arbitrary mating point Y (v_gy1 for a pinion and v_gy2 for a gear) are defined as a difference between the tangential velocity components v_ty1 and v_ty2 (v_gy1 = v_ty1 − v_ty2; v_gy2 = v_ty2 − v_ty1). From triangles O₁T₁C and O₂T₂C and with consideration O₁C = d_w1/2 and O₂C = d_w2/2, it follows:

$$v_{\mathrm{gy}1}=-\hspace{0.33em}\left({\mathsf{\omega }}_1+{\mathsf{\omega }}_2\right)\cdot \overline{\mathrm{YC}};v_{\mathrm{gy}2}=+\hspace{0.33em}\left({\mathsf{\omega }}_1+{\mathsf{\omega }}_2\right)\cdot \overline{\mathrm{YC}}$$

(9)

Because the sum (ω₁ + ω₂) is a constant value, it is clear from Eq. (9) that the sliding velocity v_gY only depends on the distance of point Y from the pitch point C. When an arbitrary mating point Y coincides with pitch point C, the sliding velocity equals zero.

Based on the sliding velocity v_gy1 in an arbitrary mating point Y, the specific local sliding of gear flanks ζ is defined as [44]:

$$\zeta =\frac{v_{\mathrm{gy}1}}{v_{\mathrm{ty}1}}=\frac{\left(1+\frac{z_1}{z_2}\right)\overline{YC}}{r_{w1}·\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_t+\overline{YC}}$$

(10)

Following the calculation procedure according to the VDI 2736, the wear behaviour of gear flanks is position-dependent due to different sliding velocities along the engagement line of mating gears. It is clear that the sliding velocities and, consequently, the wear of gear flanks are the highest at the starting and ending engagement points and then decrease toward pitch point C.

3. Results and discussion

Figure 6 shows the wear of the analysed gear made of POM after 10,000 loading cycles considering Model 1, where the geometry change due to wear of the gear flank was not considered in the numerical simulation. It is clear that the depth of the wear is the highest at the starting and ending engagement points, which corresponds to the theoretical background as described in Section 2.4. Similar conclusions can also be made for Model 2, which considered the geometry change during the numerical simulations (see Figure 7).

Figure 8 shows the comparison between the initial and worn geometry of the analysed POM gear after 10,000 loading cycles (a scaling factor of 3 was considered). It is clear that the main differences appear at the starting and ending mating points, where the wear of the gear flank is the highest.

Figure 9 shows the comparison between computational results (for Models 1 and 2) and the standardised approach according to VDI 2736. It should be noted that the VDI 2736 is based on the wear properties of the materials obtained from a rotating pin-bearing disc test, whereas the final equation for the calculation of the wear depth is based solely on the involute profile of the analysed gear. Based on the obtained computational results using the proposed computational model and analytical results according to VDI 2736, it can be concluded that both approaches give a similar wear behaviour of the analysed POM gear (the maximum wear appears at the starting and ending mating points). It is also evident that the computational results according to Model 1 (geometry changes due to wear of the gear flank were not considered) better describe the wear behaviour if compared to VDI 2736. However, it should be pointed out that, in real applications of the gear drive, the geometry of gear flanks constantly changes due to wear, which may affect the subsequent wear phenomenon of gear flanks significantly. By adding sensible updates to the geometry, we can get closer to more realistic wear values. Furthermore, it is possible to observe the change in wear compared to the model without a geometry update. In Figure 7, the green line shows the path of the wear curve of the Model 2. In the computational model, oscillations of the wear depth values between adjacent nodes have occurred. Wear is dependent directly on the contact pressure field, which means that if the contact pressure field is smooth, then the calculated wear field is smooth, and vice versa. Therefore, a smoothing function was applied to the contact pressure field, which acts as an average filter, averaging the results between adjacent nodes. For that reason, the obtained wear curves are corrected with dotted lines in Figure 9. The values at pitch point C are slightly higher, but this is not only due to smoothing but also due to the new tooth geometry, which is no longer ideal because of the tooth deformation. This results in sliding instead of rolling. It is possible to observe a marked change at the root of the tooth. The wear starts to decrease with the number of cycles also due to the change in geometry. Relief apparently occurs at the root of the tooth as the maximum pressure is reduced or redistributed.

4. Conclusions

In the framework of this study, the comprehensive computational model to analyse the wear behaviour of polymer gears was developed and evaluated on the spur gear pair, where the pinion made of POM has meshed with a support gear made of steel. A proposed computational model was built using PrePoMax finite element numerical calculation software considering Archard’s wear model to predict the wear behaviour of analysed contacting surfaces. Based on the obtained computational results and their comparison to the analytical results according to the VDI 2736 guideline, the following conclusions could be made:

• The proposed computational model could be used to simulate the wear behaviour of contacting mechanical elements like gears, bearings, etc.
• The use of the boundary displacement method (BDM) in the framework of PrePoMax software provides information on the geometry change that alters the operating conditions in subsequent load cycles.
• Due to the contact problem and the formation of contact pressure peaks in the finite element mesh, it is necessary to consider mesh smoothing in the model to allow a smooth distribution of wear over the surface. This avoids additional convergence problems in the use of BDM, but it does result in the averaging of values in locations where the differences in wear between adjacent element nodes should actually be larger (such as pitch point) and is not entirely correct. The validity of the results should, therefore, be checked with a model without smoothing.
• With the computational model using a multi-step geometry update, more accurate results were obtained, which show a reduced depth of wear at the root of the tooth. In the pitch point region, a non-zero value of the wear depth appears as the number of cycles increases. The tooth deflection and the new tooth flank geometry have a major impact on this.
• The main advantage of the model, if compared to the standardised procedure according to the VDI 2736 guidelines, is the geometry updating after a certain number of loading cycles, which enables a more accurate prediction of wear behaviour under rolling/sliding loading conditions.
• In the future, a comprehensive 3D numerical model will be developed to analyse the interaction of meshing gears. This model aims to compare the results obtained from experimental testing with those derived from numerical simulations.
• Potential future advancements of the existing Archard’s wear model could include incorporating a contact stress-dependent wear coefficient, as demonstrated in prior research. Currently, the wear coefficient remains constant throughout the analysis. Additionally, the model could be enhanced by establishing a relationship between the wear parameter and wear depth in order to account for diverse surface improvement techniques.