Modelling Friction and Abrasive Wear of Elastomers

1. Introduction

Tribological modelling of elastomer parts has high importance in engineering practice due to their widespread industrial use. Generally, the experimental investigation of the wear behaviour is time consuming and expensive, which led to the development of numerical techniques.

In order to properly model the friction and wear of elastomer materials, it is essential to understand some specific mechanical and tribological features. Mechanical characteristics of elastomer materials are different from those of metals and other relatively hard materials. Elastomers can bear large deformations, the stress-strain curves are non-linear, time and temperature dependent. When the elastomer is deformed, a part of the strain energy is dissipated. This feature makes the elastomer parts fulfil damping functions. The contact mechanical aspects of elastomers are also important for modelling. When an elastomer part is pressed against a counter-surface, the contact area would be larger than that estimated by the Hertzian theory, because of the effect of adhesion [1]. For the same reason, when the elastomer is contacting a rough surface, the real area of contact would be bigger than it is common for harder materials, as the elastomer tends to fill the valleys of the rough surface.

Friction of elastomers is usually considered as a joint effect of the friction on the surface by adhesion (and fluid film shearing in lubricated friction) and the internal friction generated by the hysteresis of the deformation by the counterpart [2]. Since the hysteretic contribution is caused by the visco-elastic properties of the material, it depends on the temperature and on the exciting frequency. When the elastomer is sliding on a hard, rough counter-surface, the asperities of the rough surface repeatedly deform the elastomer and so cause hysteresis loss. Thus the hysteresis contribution is dependent on the sliding speed. The hysteresis caused by the surface roughness of the counterpart is usually referred as micro-hysteresis. Hysteresis can occur on the macro level as well, where it can be a significant part of the friction.

When designing elastomer parts for tribological applications, modelling friction is inevitable in order to determine the contact conditions, which usually have an important effect on the functionality of the part. Wear, on the other hand, has less importance in short term use (an average passenger car tire can run up to 50 000 km before the need to change). However, in most cases, the reason for the elastomer parts to be replaced is the wear. In long term or when the lubrication is not ideal wear can dramatically change the geometry of the part and thus the nature of the contact. Therefore wear modelling is needed to estimate the product lifespan and to schedule maintenance operations. Beside experimental tests, which are cost and time consuming, numerical wear simulation techniques were developed in the last decades and became more and more commonly used. Although there are wear simulation methods incorporating various numerical methods, such as the boundary element method [3] or the discrete element method [4], the most successful and popular approaches deal with the finite element method, since it is a general method for mechanical stress analysis.

The most popular FE based wear simulation method was developed by Pödra and Andersson [5]. In this iterative method, firstly the contact pressure distribution is determined by FEM and the nodal wear increments are calculated from the wear equation of Archard [6]. Then the contacting nodes are moved with respect to the nodal wear values. Finally the FE contact calculation is carried out again with the modified mesh and the cycle is repeated according to the simulated wear process. This method was improved by several researchers (e.g. [7,8]) for various applications, different geometries and materials, such as metals, ceramics, polymers and composites. Kónya and Váradi in [9] improved the method in order to consider heat generation and time dependent material properties in the wear simulation.

Even though this chapter shows methods for modelling and simulation of elastomers, some experimental work is inevitable for the calculations. The parameters of the material models are usually determined based on experimental data. In my simulations I used the material model described by Pálfi et al. in [10]. It is a Mooney-Rivlin model with an attached generalized Maxwell model. This material model was constructed based on stress relaxation measurements. The tribological properties, like the coefficient of friction and the specific wear rate also require experimental tests. The basic mechanical and tribological properties of the investigated materials are summarised in [11] and in [12]. The coefficient of friction depends on many factors, yet the lubrication has one of the most important effects. Figure 1 shows the result of a starvation test, which was carried out at the Institute for Composite Materials, Kaiserslautern, Germany on a shaft-on-plate test rig (Dr. Tillwich GmbH, Horb-Ahldorf, Germany). In this configuration the rubber sheet was pressed against the rotating steel shaft. The normal and the tangential force were measured online. In this experiment lubricant (Hydraulan 407-1, BASF, Ludwigshafen, Germany) was applied only at the beginning of the test. Later, as the amount of lubricant decreased between the rubber and the steel shaft, the nature of the lubrication changed. In the beginning, a mostly constant coefficient of friction with a value of µ = 0.06 lets us assume soft EHD lubrication regime (the green section in the figure) [13]. In the next two phases, the coefficient of friction is increasing, which means more and more solid body contacts as the lubricant decreases (blue and orange). In the last phase, the coefficient of friction has an approximately constant value that is higher than one, which implies technically dry friction (red).

The presented finite element (FE) models were generated with MSC.Marc Mentat software package; the calculations were done by the MSC.Marc Solver. The developed wear simulation algorithms operate the MSC.Marc model and result files, however they are independent of the FE environment and can easily be adapted to other FE systems.

Figure 1.

The coefficient of friction with starving lubrication (v = 250 mm/s F_n = 28 N)

2. Modelling internal friction of elastomers

During deformation of elastomer materials, because of the visco-elastic material behaviour, a part of the applied strain energy is transformed to heat as a result of internal friction, hysteresis [14]. The amount of this energy loss can be characterized by the loss factor (tan(δ)), which is the ratio of the loss and the storage moduli of the material. When repeated loads are present, the hysteresis contribution to the friction is more significant. In case of sliding friction between elastomer and a rough rigid counter-surface, the elastomer is subjected to repeated, cyclic deformation by the asperities of the rough surface. This kind of internal heat generation can lead to fatigue wear [15].

2.1. Modelling micro-hysteresis

In a series of numerical calculations the micro-hysteresis was investigated. First the hysteretic friction caused by a single asperity was studied. The asperity was modelled by a rigid cylinder with a radius of 250 µm (Figure 2). The asperity was pressed and rubbed against a rubber block with varying sliding speed. In this model the surface friction was neglected, so only the friction caused by the hysteresis was considered.

Figure 2.

The model of one cylindrical asperity sliding on rubber

Figure 3 shows the results of the numerical simulations. The coefficient of friction was calculated as the ratio of the tangential and the normal reaction forces. It can be seen that the velocity dependence of the hysteretic coefficient of friction is in agreement with the measured loss factor curve.

Figure 3.

The hysteretic coefficient of friction vs. the exciting frequency and the corresponding loss factor curve

2.2. Modelling macro-hysteresis

A 3D finite element model was created in order to model a Pin-on-Plate test described in reference [12]. A 4 mm thick rubber sheet was modelled by 7410 3D brick elements (Figure 4). The hemispherical steel slider (r = 5 mm) was modelled as a rigid body. The rubber sheet was secured on a steel substrate, which was modelled as a rigid body too. The assumption of rigid bodies is reasonable if one considers the fact that the steel has four orders of magnitude higher elastic modulus than the rubber specimen.

Figure 4.

finite element model of a Pin-on-Plate configuration (the elements are blue; the surfaces considered to be rigid are magenta). Adaptive local remeshing was applied in the contact area.

The length of the elements varied between 0.5 and 2 mm. The hemispherical counter surface was moving along a circular track of 33 mm diameter at various sliding velocities. F_n = 100 N normal load acted on the substrate, which pressed the rubber sheet against the slider.

No displacement was possible between the substrate and the specimen; the coefficient of friction specified between the specimen and the counter surface was µ = 0.06 in order to take into consideration the friction developing on the lubricated rubbing surfaces (See Figure 1).

A series of calculations were run on this FE model with different sliding velocities (in a range between 3 and 1036 mm/s).

After the calculations, the friction corresponding to each velocity value was determined from the finite element result files. Changes in the coefficient of friction are shown in Figure 5. The figure also indicates the corresponding measured results.

A slight decrease in coefficient of friction can be observed both in the calculation and in the test results. Note that these calculations did not consider the changing viscosity of the lubricant. Any changing of the friction can be attributed to the changing of the internal friction, since the surface friction was modelled by a constant prescribed coefficient of friction. The decreasing friction can be explained by the decreasing penetration depth of the slider and the decreasing contact area. The higher the sliding velocity, the lower the penetration depth, because of the visco-elastic behaviour of the material (the rubber seems harder at higher strain rates). This causes decreasing in the excited volume, and the hysteresis part of the friction depends highly on the excited volume. Of course the hysteresis – as it was already mentioned – depends on the exciting frequency as well, however in the investigated sliding velocity range the resultant frequency (approximately from 0.3 Hz to 103.6 Hz) does not change the tan(δ) and the coefficient of friction significantly.

Figure 6 shows the contact area and the penetration depth of the counter surface for the sliding velocities examined. It can be established that the counter surface moving at a higher velocity can be pressed into the rubber surface to a lesser extent; thereby the contact area will also be smaller. There is a 0.25 mm difference between the penetration values corresponding to the lowest and the highest velocity.

Figure 5.

Coefficient of friction as a function of sliding velocity

Figure 6.

Contact area and penetration depth of the counter surface in function of sliding velocity

3. Modelling wear of elastomers

Three different wear simulation approaches are presented in the followings. The first one is the method of moving the contact nodes according to the nodal wear increments as proposed by Pödra and Andersson in [5]. It will be shown, that this method is suitable for modelling wear, which is much smaller than the elements in the FE model. The second and the third methods can model relatively large wear by global remeshing and by deactivation of certain elements, respectively.

3.1. Wear simulation by moving the nodes

The wear simulation method of Pödra and Andersson was adapted to model the wear process of an elastomer part considering its non-linear material properties. The simulation process is shown in Figure 7. In the modelled configuration a rubber sheet (thickness: 2 mm) was pressed against a rotating steel cylinder (diameter: 10 mm) with a normal force of 16 N. The rubber part was modelled in 2D assuming plain strain conditions (Figure 8). The rotating steel cylinder was modelled as a rigid body. A coefficient of friction with a value of 1 was prescribed between the rubber and the cylinder, according to the dry wear test results.

After the FE contact calculation, the wear increments were determined based on the stress distribution. The Δh_i nodal wear of the i-th node in the contact area by Archard [16] is defined as:

$$\Delta h_i=W_s{p}_{i}v\Delta t$$E1

where

W_s is the specific wear rate,

p_i is the nodal contact pressure at node i,

v is the sliding velocity,

Δt is the time increment.

In the next step of the simulation the contacting nodes were moved with respect to the nodal wear increments. Then the cycle was restarted with a new contact calculation using the modified geometry.

Figure 9 shows the result of the wear simulation after one minute of sliding. The profile of the worn specimen is also plotted in the figure. It can be seen that the wear is not symmetric. The simulation and the experimental results are in good agreement both by the means of the amount of wear and the worn shape. However, in longer term simulations and with finer FE mesh, this method failed to accurately estimate the experimental results.

Figure 7.

The flowchart of the wear simulation

Figure 8.

The finite element mesh of the shaft-on-plate configuration

Figure 9.

The calculated and measured worn profile of the rubber sheet (t = 60 s, F_N = 16 N, ω = 50 rad/s)

Despite the widespread use of the Pödra-Andersson method, the wear simulation is highly limited; only the top layer of elements can be worn. This limitation does not affect the usability of the method in case of relatively hard materials such as metals or ceramics, since minimal wear can have a great effect on the contact area and pressure distribution and thus on the performance of these parts (e.g. bearings). In case of rubber parts, wear in the magnitude of some μm does not change the contact conditions significantly. The wear needed for a rubber part to malfunction is much greater than those of metal parts; therefore the wear of the top layer is insufficient in rubber applications. To increase the volume of the wear to be modelled, the size of the elements can be increased; however it is not recommended, since it decreases the accuracy of the calculation. In the followings two methods are proposed to simulate wear regardless of the element size. It will be shown that these methods can model the material loss due to wear even in a scale comparable with the size of the part.

3.2. Wear simulation by global remeshing

In order to model wear that is larger than the elements of the FE mesh, a wear simulation procedure was developed using global remeshing. To demonstrate this method a wear simulation study of a reciprocating sliding seal (Figure 10) is presented. In the investigated application a rubber seal is coupled with an aluminium rod (diameter: ø22.2 mm). Since this method integrates the FE mesh generation, the geometric model is the starting point of the wear simulation. The seal was modelled assuming axisymmetry.

The counter-surface was modelled as rigid body, since its elastic modulus is some orders of magnitude larger than the elastomer material of the seal. The seal was mounted axially by fitting it in the housing, which was considered ideally rigid. The rod was moving with alternating motion at a speed of 20 mm/s with 9 mm amplitude. The contour of the seal section was modelled by line segments (Figure 11). A fine approximation in the vicinity of the lip as well as the ridges was required; in the region of the ridges and the lip the average length of the line segments was 10 µm.

The seal section was meshed by the built-in automatic meshing procedure using three-node axisymmetric triangular finite elements (Figure 12). In the FE calculations the interference fit between the seal and the rod was also taken into consideration: the unloaded inner diameter of the seal was 1.8 mm smaller than the diameter of the rod. This caused the seal to be stretched in the first load step. In the next step the alternating motion was modelled considering the friction during the wear simulation. According to the working conditions, in the inward strokes of the rod a pressure load of p = 2 MPa was applied on the right side of the seal, modelling the sealed pressure. There was no pressure applied in the outward strokes. Between the seal and the rod a prescribed coefficient of friction with a value of µ = 0.1 was defined in order to model the lubricated friction.

In the wear simulation process, the FE mesh was created based on the line segments describing the geometry. The automatic mesher was set to create elements of different sizes, so the elements in the vicinity of the contact zone were smaller than those deep inside the material. In the initial mesh approximately 4000 elements were distributed in a way that the regions, that are important for the wear calculations, were modelled more accurately with the small elements. Even so the total number of the elements remained in a reasonably low range, so the CPU time for one calculation was short enough to handle the iterative simulation process.

Figure 10.

Section view of the investigated reciprocating seal

Figure 11.

The geometry of the seal consisting of line segments and points

Figure 12.

Deformed shape of the seal in the first outward stroke (without pressure). The curves represent the rigid bodies (the rod and the housing). The element size varies from 10 to 200 µm

The flowchart of the simulation can be seen in Figure 13. In the first step the FE mesh was created based on the initial geometry. After adjusting the simulation parameters (coefficient of friction, sliding velocity, material properties, applied pressure, time increment, etc.) the contact calculation was run. In the first stroke the rod was moving rightward. The wear calculation was performed using Equation 1. The moving of the contacting nodes according to the nodal wear increments was the next step, however in a slightly different way. The surface nodes were attached to the contour points of the seal, so if a point was moved the attached node would also move with respect to the nodal wear value of the attached node. After moving the points, the whole FE mesh was deleted. Based on the new geometry the automatic mesher created the new elements and the cycle ran over again. The direction of the rod motion was changed in each cycle, so one simulation cycle represented one outward (left) or one inward (right) stroke. The applied pressure was also changed stroke by stroke, as for the inward strokes the working pressure of 2 MPa was applied, while in the outward strokes the working pressure was turned off in the simulation process. 200 simulation cycles were calculated in the frame of this study.

Figure 14 shows the deformed shape of the seal in the first outward stroke, when the seal was not pressurized. The resulting contact pressure comes from the stretching of the seal on the rod. It can be seen that the contact area is reduced: only five ridges on the left side, one on the right side and the lip are in contact with the rod. The maximum of the contact pressure can be observed at the lip. The values of the contact pressure are generally low, since there is no working pressure applied to press the seal surface to the rod.

Figure 13.

The flowchart of the wear simulation process

The calculated worn profiles of the seal are shown in Figure 15 after 0, 100 and 200 cycles of the simulation. It can be seen that the left side ridges of the seal and the lip edge wear first. Figure 15 also shows the contact pressure distributions in the investigated wear phases. Note that the FE mesh is different in each state, since the model was remeshed in every cycle; however the element sizes are distributed identically. It should be mentioned that the total simulated wear depth reaches 30 µm in some regions (e.g. at the lip), which is three times higher than the size of the elements in the contact area. This wear depth is so large that the Pödra-Andersson method can not be applied here.

Figure 14.

Deformed shape of the seal in the first outward stroke (without applied pressure) and its contact pressure distribution.

Figure 15.

The calculated worn profiles after 0, 100 and 200 cycles of the simulation in inward strokes (p = 2 MPa) and the corresponding contact pressure distributions along the axis

Based on the change of contact pressure distribution the followings can be established. The area of contact increased, while the maximum values of the contact pressure decreased. This phenomenon can also be seen in Figure 16. for the contact pressure distribution of the lip over the wear process. One can see that the seal at the beginning has about three times higher contact pressure than later. The pressure reduction is caused by the wear process, which also increases the contact area. It is remarkable that at the left side of the lip the contact pressure is increasing at first, and then slightly decreasing. It can be explained by the wear of the lip. As the lip has the highest contact pressure in the beginning of the simulation, it will suffer the most severe wear in the first cycles. As the lip wears, the contact pressure distribution becomes flatter, which causes the slight increase in the lower pressure region. Later this mostly uniform contact pressure remains, only the average pressure decreases.

Figure 16.

The change of the contact area and the contact pressure distribution at the lip over the wear process (after 0, 50 and 100 simulation cycles)

3.3. Wear simulation by deactivation of elements

In hydraulic systems the seals are well lubricated in most cases, so those parts show minimal wear and they can function for years without any sign of damage. However, in some special cases the rubber seals can work under mixed or boundary lubrication conditions. It is suspected that sometimes (typically after the machine is stopped for a longer period) the seal can squeeze the lubricant out of the space between the surfaces and some parts come to direct contact. If it happens, adhesion can occur, which significantly increases the frictional force and the wear mechanism would change completely. The different wear mechanisms of rubbers are examined by several studies e.g. in [17, 18]. In general it can be said that dry sliding wear of rubber can be traced back to the initiation and propagation of cracks in the rubber material. The modelling of this type of wear requires advanced simulation techniques.

A popular approach for modelling material detachment is the deactivation or deleting elements from the FE mesh (sometimes also referred as element death) based on the characteristic damage criterion of the given material. It is commonly used for simulation of manufacturing processes e.g. in [19, 20]. The principle of wear – from the simulation point of view – is similar to the manufacturing processes. Authors in [21] modelled the particle detachment process based on the strain history of the elements in dry sliding.

The wear process of an experimental test of the sliding seal described above was modelled. The seal was cut and a section of it was straightened and fixed in the clamping. The aluminium rod was then pressed and rubbed against the surface of the seal. The detailed description of the test can be found in [22]. In the series of wear tests, different lubrication conditions were studied, and it was found that when the amount of lubricant was decreased, the seal surface worn rather intensively. The surface of the seal was inspected with a MicroProf white light profilometry device (WLP, Fries Research & Technology GmbH., Bergisch Gladbach, Germany) before and after the test, so the wear of the seal can be evaluated and compared with the simulation results (See Figure 17). One can see that the most wear occurs in the symmetry plane of the seal, where the contact pressure was maximal; this region shows a rough surface texture, which is probably caused by the wear particles that were not removed from the contact region, since there are peaks that are higher than the original surface. In order to avoid the false results caused by the wear debris, a plane with slightly less but still significant wear was investigated, which was clear of wear debris (indicated in Figure 17).

Figure 17.

Worn surface of the seal

The previous wear simulation methods were not suitable to model the compound wear mechanism, which was caused by boundary lubrication. Beside the wear simulation method of moving the nodes, the deactivation of elements was integrated in the simulation procedure. The flowchart of the simulation can be seen in Figure 18.

Figure 18.

Flowchart of the simulation [22]

The test can be divided into repeated cycles of the inward (right) and outward (left) strokes. Since the stress distributions in the inward and outward strokes are different, the wear calculations were performed four times in every cycle (see Figure 19).

Starting with the initial geometry a frictional contact calculation is carried out assuming inward stroke. Based on the results, the normal stress of each element in the model was compared to the predefined critical stress limit, the tensile strength, which was 19.0 MPa for the material of the seal. If the tensile stress of an element was greater than the critical value, the element became deactivated. The damage check was done in every element in the following way. If σ_i > σ_crit (where σ_i is the normal stress in the i-th element, and σ_crit is the ultimate tensile strength of the material), the element i became deactivated.

Figure 19.

The motion of the rod (continuous red) and the points of the wear calculations in a cycle (blue dots)

Wear calculation was applied to determine the nodal wear of the nodes in the contact area using Equation (1). Once the nodal wear values were known, the nodes, which were in contact, were moved to their new position as in the Pödra-Andersson technique.

The seal section was modelled by four-node 2D finite elements assuming plain strain condition. The housing, the clamping, which fixed the seal to the housing, and the counterpart, was considered ideally rigid. In the first load step of the FE calculation, the seal was mounted in the holder, and the normal load was applied to the counterpart, according to the wear test. After this, the rod started its motion. There was a prescribed coefficient of friction between the seal and the rod with a value of µ_lub = 0.3, representing the boundary lubricated friction; and another between the seal and the holder parts with a value of µ_dry = 1, since the friction condition was dry in this case. Values of coefficient of friction in case of lubricated and dry friction were determined based on the results in [12] and [11], respectively. Note that the lubrication was taken into consideration only by the lower value of coefficient of friction. Other effects – like the lubricant separating the surfaces reducing real contact and wear – were neglected. The FE mesh with the ideally rigid parts and with the applied prescribed coefficients of friction is shown in Figure 20.

The calculated worn profiles are shown in Table 1. It can be seen that the ridges of the seal and the lip edge wear first, which is in accordance with the experimental findings. One can see that during the wear simulation the wear occurs not only in the top layer, but in several layers of elements. The maximal wear depth is about ten times bigger than the typical element edge length in the contact area.

Figure 20.

FE model of the seal section before (left) and after mounting in the holder (right)

The height of the elements in contact gets reduced due to wear and over time the contacting elements become more and more distorted and can even turn inside out, which can make the calculations unable to run. To avoid this and to maintain the stability of the simulation process, the following procedure was implemented. A filter was applied in the simulation algorithm, which checked the distortion of the elements before each FE calculation. Distortion in the FE software is determined at each node of the element as the cosine of the angle between the element edges at the nodes. The distortion of the element is defined as the maximum of these cosines. If the normal vector to the element changes sign from one node to the next (inside-out element), the distortion is one. If the distortion value of an element exceeded a critical value (0.99) before the load was applied, the element got deactivated. With this method the inside-out elements and the elements that were too flattened due to the wear, were became deactivated.

The wear process of one ridge can be tracked in the successive images of Figure 21. One can see the effect of the wear algorithm on the height of the elements. The dark shaded elements were deactivated by the damage check module of the simulation.

The measured and the calculated worn profiles of the seal are compared in Figure 22. The calculated worn profiles are mostly identical with the measured worn surfaces. However, in the vicinity of the lip, the simulation over predicts the wear. It is suspected that the deviation is caused by the assumption of linear wear theory in the wear calculations, while in reality wear is a non-linear function of contact pressure. This simplification nevertheless appears to be suitable for the regions with less contact pressure.

No. of cycles	Sliding time [min]	Calculated worn profile
0	0
600	10
1200	20
3600	60

Table 1.

Evolution of the calculated worn profile (The deactivated elements are not plotted.)

Table 1.

Figure 21.

The wear process of one ridge of the seal. The dark elements are deactivated by the damage check algorithm. The images show the left ridge of the seal model after 0, 15, 30, 45 and 60 minutes of sliding, respectively.

Figure 22.

The measured and calculated worn profile of the reciprocating seal

The developed wear simulation algorithm was able to model the wear process even in the case, when the amount of wear was much larger than the size of the elements used in the FE calculations.

4. Conclusion

Friction of elastomers has a hysteretic part, which is caused by the energy dissipation in the visco-elastic material. The hysteretic friction depends mainly on the sliding velocity and the geometry. It was shown that the effect of sliding velocity on the coefficient of friction can be modelled numerically both in macro and in micro level. Furthermore, the penetration depth of the counter surface and the size of the contact area can be modelled as functions of sliding velocity.

Three different wear simulation techniques were presented for different purposes. The first method can be useful for modelling relatively small amount of wear. The second can model wear even bigger than the elements in the model by global remeshing. The third method deactivates the damaged elements, thus the effect of the surface rupture of the elastomer part can be modelled.