Hysteretic friction of rubber in tribological tests

The goal of this paper is to summarize different investigations of FE modelling of rubber hysteretic friction which was produced in the last 4 years in the KRISTAL project (Knowledgebased Radical Innovation Surfacing for Tribology and Advanced Lubrication, EU Project Reference). First, the characterization of the viscoelastic material properties of an EPDM 75 rubber by 40 term generalized Maxwell model fitted to the loss factor is presented. After it, the validation of the material model using it in FE model’s for simulating different tribological tests are produced. The DMTA test, after micro hardness tests are modelled. Then, a ball on plate tribological test was simulated by FEA and to check the reliability of the material model the test and FE results were compared. Finally, rubber sliding on a rough surface was modelled and the hysteretic friction was determined by 2D and 3D FE models.


Introduction
During the design process of different mechanical engineering applications wherein components are sliding on each other it is important to know the mechanical behaviour and tribological properties of the components. In many technical applications one of the sliding surfaces is made of rubber or a rubber-like material (such as seals contacting with a metal or polymer surface, wiper blades-wind screen and tires-road contact etc.). In case of clean surfaces and in the absence of lubricant when contacting surfaces are dry the main source of rubber friction is the adhesion and hysteresis [1][2][3]. The joint effect of them generates frictional resistance when an elastomer part is pressed and rubbed against a rough surface.If the rubber slides on a rough surface which is rigid compared to the rubber, the asperities of the surface repeatedly deform the surface of the rubber. Due to the viscoelastic material behaviour of rubber, during deformation, part of the strain energy is transformed to heat as a result of the hysteresis [3]. Generally, rubber components operate under lubricated conditions, to separate rubbing surfaces thus reducing the impact of adhesion, as well as by filling up the valleys with the lubricant to reduce the exciting effect of asperities. The shearing of the fluid film also acts as a source of friction. In general, the fluid film is extremely thin, therefore many times boundary lubrication occurs.
Numbers of articles deal with the analytic and semi-analytic determination of friction force induced by from hysteresis. Grosch [3] performs pioneering work by revealing that the two main sources of the friction force are adhesion between the surfaces and the energy loss generated in the material, that is, hysteresis. In addition, he stated that both physical phenomena are closely connected to internal friction (hysteresis) of the rubber. In the papers of Persson and Klüppel [4,5] they investigated the rubber when it is sliding on a hard, rough substrate and energy dissipated inside the rubber due to the surface asperities of the substrate exert oscillating forces on the rubber surface. Hysteretic friction was determined and the results were compared with the experimental data of Grosch. Persson stated that rubber friction on rough surfaces in presence of lubricant is mainly due to the viscoelastic deformations of rubber (hysteresis).
In the literature only a few paper deal with the hysteretic friction prediction by the finite element method. For instance, in [6,7] the hysteretic friction at asperity level was studied by FE technique.Nevertheless, the reliability of hysteretic friction prediction caused by rough surfaces is in close connection with the effectiveness of topographical analyses. Beside the parameter based techniques, nowadays two dominant research trends can be observed. One is the "global" surface characterization method using complex mathematical tools, such as power spectral density (PSD) analysis, while the other is the technique to characterize the local features of topographies based on the identification of asperities and scratches. Due to the development of characterization techniques remarkable improvement has occurred in the field of measurement techniques: for example, the atomic force microscopy (AFM) extents the surface characterization to nano-scale [8].

Experimental
The mechanical behaviour of rubber-like materials is principally characterized by a non-linear stress-strain curve and time and temperature dependency.

DMTA measurement and results
In order to characterize the viscoelastic nature of the observed peroxide cured EPDM rubber with a ≈ 65 phr carbon black content (CB N347), IRH 75 (referred further on as EPDM 75) DMTA measurements were performed. The DMTA measurements were carried out with equipment GABO Eplexor 100N (Ahlden,Germany) at TU Kaiserslautern in the IVW (Institut für Verbundwerkstoffe GmbH). During the test the complex modulus (E*) of the material, its storage (E') and loss (E") parts, as well as its mechanical loss factor (tan(δ)) were recorded. Later, the master curve of the observed material can be constructed from the DMTA test results to characterize the viscoelastic properties of the rubber. Details of the test are written in the [9].

Micro hardness test and results
Micro hardness tests were performed on the EPDM 75 material using a dynamic ultra micro hardness testing set (Shimadzutype DUH 202) on the 1 mm thick rubber samples. Due to the measurement three different load modes were used Fig. 1 [10]. The micro hardness test procedure begins with the surface detection. During this procedure the diamond pyramid approaches and touches the investigated surface until the contact force reaches the 2 mN value. The loading and unloading speeds were the same 1.4 mN/sec in each loading. In mode-1 and mode-3 after reaching the maximum loads, 2 sec holding time was applied, the unloading starts after the holding time. In mode-3 the unloading was also followed by 2 sec holding time. The details of the micro hardness test can be seen in [10].
The results of the micro hardness test (and the later described FE results) can be seen in Fig. 6. For the first case (mode-1) a 2 sec long holding time was located between the uploading and downloading phase Fig. 6a while the material undergoes a creep process. After the load is completely removed it can be seen that the investigated material undergoes retarded deformations due its viscoelastic nature and the surface of the material doesn't recover to its original position. In the case of mode-2 micro hardness tests hysteresis and hardening can be observed Fig. 6b. In mode-3 Fig. 6c. Mullins-effect like process can observed, superimposed with creep during the 2 sec holding times.

Ball on plate test and results
The ball on plate tests were performed on a CSM INSTRU-MENTS microtribometer [11] at TEKNIKER. In the ball on plate configuration test a steel ball of 2 mm diameter was pressed against to the EPDM rubber plate of 10 mm×4 mm×2 mm. Then the ball was forced to perform reciprocating motion with an amplitude of A=0.3 mm at different, sinusoidally varying speeds. The normal force was F n = 100 mN while the maximums of the sinusoidally varying sliding speeds were v max = 0.1, 1, and 10 mm/s. During tests, the frictional and normal forces were measured under lubricated conditions.

40-term generalized Maxwell model
Three different 40-term generalized Maxwell models fitted to the loss factor (see later), considering three temperatures (T = −50, 20 and 150 • C), were constructed based on the DMTA measurement results.

Construction of the master curve using the WLF-Arrhenius equation
The temperature vs. time equivalence developed for polymers can be described by the Williams-Landel-Ferry (WLF) equation [12][13][14][15]. The WLF equation can be used for determining the shift factors compared to the curve associated with the respective reference temperature, thus a master curve for the corresponding material can be constructed.As we are going to study the material behaviour at three different reference temperatures, master curves associated (T = −50, 20 and 150 • C) were constructed, by horizontal shifting the measured curves, in case of DMTA tests. In Fig. 2 the curves represent the master curves pertaining to the three reference temperatures, respectively, and indicate the storage (continuous line) and the loss modulus (dashed line) as a function of excitation frequency. Fig. 2 show the course of the ready master curves. It can be seen that the value of the loss modulus (E") does not increase continuously as the frequency rises, but after reaching a maximum value it begins to decrease and approaches a plateau value [16]. The loss factor (tan(δ)) master curve is derived as the ratio of the loss (E") and storage (E') modules Fig. 3. Due to this it has also a global maximum value near, to the maximum of the loss modulus. Based on the constructed master curves three different 40term generalized Maxwell models were created for the three different reference temperature fitted them for the loss factor (by using the ViscoData software [17]).The details of the generalized Maxwell model fitting can be seen in [9,10,[18][19][20].The non-linear stress-strain behaviour of the rubber was modelled with the frequently used Mooney-Rivlin model applying two parameters [21] with values of c 10 = 406.66 MPa and C 01 = 106, 66 MPa (E 0 = 3050 MPa) The construction of the generalized Maxwell model in detail can be found in [10,[18][19][20].

FE simulation of the micro hardness test
In the following chapter the presented FE models were created by MSC.MARC 2007(r1) FEA software.

The FE model of the DMTA test
The material model described was used in a simulated DMTA test by FEA to model the dynamic behaviour and dissipated energy in the rubber. Geometry of the FE model was created according to the sample size in the DMTA measurement considering 1/8 model due to the symmetry conditions. To describe the rubber sample 8-node incompressible Herrmann elements were used [21]. The applied loads were ε s = 0.06581% in static case and ε= 0.01% in dynamic case, while the observed frequency range was f = 10 −5 − 10 15 Hz.
from the material model to the FE modelling The result of the simulation was the shift ( t ) between the stress and the strain responses. The phase shift (δ) loss can be calculated from the shift by [16]: The FE results for the loss factor follow the generalized Maxwell models at every frequency In the 40-term generalized Maxwell model fitted to the loss factor by trial and error technique, the loss factor shows an acceptable correspondence with the measured loss factor curve, but we have to note that due to the smaller underestimation of the storage modulus the estimation of the dissipated energy will be slightly higher than the real value. The FE model of the DMTA test and their results can be found in detail in [9].

FE model of the micro hardness test
In order to verify the generalized Maxwell model, FE models were prepared to simulate the micro hardness test. The scheme of the micro hardness test can be seen in Fig. 4. In the FE simulations the surface touching process is also considered and involved as in the measurement. The FE mesh of the micro hardness tests considering axisymmetric model is shown in Fig. 5.
To model the rubber sample 4-node incompressible Hermann type elements [21], while to model the diamond indenter a rigid surface were used. Due to the axisymemetric FE model the inclination angle of the cone had to be changed to 19 • to represent the same contact surface as in case of the diamond pyramid in the measurement.
The results of the FE calculation can be seen in Fig. 6. It can be concluded that both fitted Maxwell models can describe the material behaviour of rubber in case of a micro hardness test. The FE model predicts larger penetration depth compared to the measurement. But the creep of the material during the 2 sec holding time was simulated successfully. The reason can be found in the penetration frequency. The penetration frequency can be evaluated from the penetration speed and the Although the fitted Maxwell model can describe adequately the upload phase in case of mode-1, and can follow the relaxation phase of the load curve (when the force kept at a constant level) but in case of download phase the material model shows a higher 'stiffness' considering the calculated load at the diamond indenter. This effect shows an error in the penetration of the indenter at the end of each load cycles, that is cumulated. That is why one can see handsome difference between the measured and calculated penetration in repetitive cycle type of loads Unfortunately this material model underestimates the storage modulus because of the manual modification [20] and this is the reason why this model can predict too large penetration. The micro hardness test and the FE model results can be found detailed in [10].

FE model of the ball on plate test
In the FE model of the ball on plate configuration to model the rubber plate incompressible Herrmann type elements (υ=0.5), while to model the steel ball an ideally rigid component was used. Only half of the rubber plate was modelled using symmetry condition, so the nodes in the symmetry plane were fixed in z-direction Fig. 7. The bottom of the rubber plate was fixed, i.e. the nodes on this plane were constrained in x-, y-and zdirection. At first, the steel ball was pressed (in negative ydirection) into the rubber plate with the force specified, that was built up linearly. Afterwards, the ball was drawn horizontally (in negative x-direction) at sinusoidally varying speed according to the measurement, at an amplitude of A = 0.3 mm, through six cycles. The results of the simulation were the friction force arisen from the hysteretic friction and the calculated coefficient of friction.
Two models were produced. One of the models was made to verify the FE model, so that a coefficient of friction -preset by an iterative method -was taken into consideration to character-  ize fluid film shearing and boundary lubrication. This prescribed coefficient of friction had to be set in a way that its value -superimposing with the coefficient of friction caused by hysteresis, as calculated by the FE software -should be equal to the measured value of the coefficient of friction pertaining to sliding. By setting coefficient of friction of µ= 0 in the other model, calculations were performed to determine the hysteresis component of the coefficient of friction. Based on preliminary calculations the change of the coefficients of friction calculated can be considered as "steady state" by cycle 6.In Fig. 8, the black lines represent the coefficient of friction, taking into consideration the fluid film shearing, boundary lubrication and the hysteresis required during iterative calculations. The separate rhombics and square symbols show the values of the measured coefficient of friction for two subsequent measurements; the circles connected by a continuous grey line represent the results of the calculation at µ= 0. By this iterative approach we could identify the portion of the coefficient of friction due to hysteresis as well as fluid film shearing and boundary lubrication. Fig. 8 show that, comparing to the measured coefficient of friction the FE iterative calculation can approximate the friction behaviour of the sliding pair.
To see the impact of the sliding speed on the predicted hysteresis the sliding speed range was extended 0.01-100 mm/s, thus the excited frequency range was between 10 −2 Hz and 10 +2 Hz in the FE simulation. Fig. 9 show the hysteretic friction and the penetration depth of the ball in this extended sliding speed range. The highest hysteresis can be observed at 150 • C at each speed, due to a low elastic modulus and a significant loss factor. The highest hysteresis was calculated at 1 mm/s. At this sliding speed (exciting frequency about 1 Hz) the peak of the hysteresis is corresponds to the small peak of the little bit oscillating loss factor (see Fig. 3. At 20 • C the same peak of the loss factor can cause very similar effect but to a smaller extent.
In order to understand the results at −50 • C, relative to the ones at room temperature two different tendencies should be followed. The loss factor is greater by nearly one order of magnitude and E' is greater by more that one order of magnitude.
The loss factor increases nearly proportionally the friction force, while the larger E' highly reduces the volume of material subjected to deformation. This volume is proportional to the penetration depth. In the FE studies, at first, the ball was pressed into the rubber sample. At −50 • C the depth volumes is lower nearly by one order of magnitude compare to the penetration depth at room temperature, yielding to much smaller volume deformed, affected by the loss factor, which is greater. The joint effects of these two different tendencies characterize the hysteretic friction. It can be stated that the material model models the loss factor accurately but underestimates the measured storage modulus. In the case of increasing frequencies, the degree of underestimation gradually decreases. As a result of the underestimated elastic modulus, the results in Fig. 9 for room temperature and 150 • C can be considered as upper estimates as regards the COF. More in detail about this investigation can be found in [17,18].

Hysteretic friction prediction in case of rubber sliding on a hard rough surface
In order to examine the hysteretic friction when smooth rubber block slides on a hard, rough surface, 2D and 3D FE models were constructed. Firstly, to characterize the surface roughness of a plunger of a brake cylinder surface roughness measurements were done using a Mahr Perthometer Concept type stylus instrument and an Atomic Force Microscope (AFM).
The PSD (Power Spectral Density) analysis, that decomposes the surface roughness into harmonic components, makes a connection between the real, measured surface roughness and the surface roughness model applied. With the help of the PSD analysis of a measured rough surface, one can determine the characteristic wavelengths and amplitudes of its harmonic components. Using the harmonic components of a measured rough surface one by one or their combination one can examine the contribution of each component of the surface roughness (from microto nano-level) to the hysteretic friction. The surface roughness measurement and analysis is described in [8]. The results of the surface characterization show that the dominant wavelengths (λ) are in the range of 400 nm −250 µm.
In the case of the 2D FE model an assumption was made. The microscopic surface roughness of the rough surface was modelled by two different sine waves having a wavelength of 100 µm and 11.11 µm (called them surface A and B) where the peak to peak distance (double of the amplitude of the sine waves) were 8 and 2 µm, respectively. Furthermore, by the combination of surface A and B (surface A+B) where the surface B was superimposed on surface A. The exciting frequency considering 10 mm/s sliding speed 100 Hz in case of surface A and 900 Hz in case of B surface. Due to the periodicity in the surface roughness it is sufficient to model a small, repetitive segment of the rubber. For this purpose repetitive symmetry was used. The rubber has been modelled by plane strain QUAD80 elements, while the rough surface has been modelled as an analytical rigid surface [21].
The sliding contact between the rubber and the rough surface was modelled as follows. Firstly, the rubber segment was compressed with a pressure of p= 1 MPa against the surface by using incremental techniques. During the indentation the relative tangential velocity between contacting bodies was zero. The second step was the acceleration of the rubber block up to a sliding speed of 10 mm/s. During this the rough surface was fixed in vertical direction and the nodes on the lateral walls of the rubber were forced to move identically (see Fig. 9). The third one was the horizontal motion with a constant speed of 10 mm/s, where boundary conditions were the same as in the previous load case. The pressure on the top wall was also applied during the second and third step.
It can be established on the basis of the coefficient of friction in Table. 1 that in the case of surface A the highest friction coefficient was obtained at room temperature while the lowest one at −50 • C. In the case of surface B, a similar tendency can be observed; however, the highest coefficient of friction can be observed at 150 • C. In case of surface A+B, it can be established for all three temperatures that hysteresis will be greater than at surfaces A or B separately. Also in case of surface A+B, the greatest hysteresis loss can be calculated at room temperature. The 2D model and their results can be found in detail in [20].

Conclusions
The time-and temperature-dependent material behaviour of the investigated EPDM 75 rubber was measured by DMTA equipment. Using the temperature-time equivalence principle, master curves were constructed from the measurement results considering three different temperatures (−50 • C, room temperature and 150 • C). After that, three different generalized Maxwell models were constructed. Using the Maxwell models, a series of FE models has been created to investigate the viscoelastic behaviour of the rubber in a broad frequency range.
Based on the results of the FE simulation of the DMTA test, it can be concluded that however the 40-term Maxwell model fitted to the loss factor master curve shows a smaller underestimation considering the storage modulus but describe the energy loss accurately in case of EPDM 75 rubber.For the verification of the material models an axisymmetric FE model was built to simulate micro hardness test for the investigated EPDM rubber. The 40-term generalized Maxwell model fitted to the loss factor predicts larger penetration depth compared to the measurement. But the creep of the material during the 2 sec holding time was simulated successfully.3D FE model has been developed for the Fig. 9. The coefficient of friction in function of the distance covered by the ball,at F=100 mN load and v=0.1 mm/s maximum speed prediction of hysteretic friction force in case of a reciprocating steel ball sliding on a rubber plate. The 40-term Maxwell model described the material behaviour at −50 • C, room temperature and 150 • C. At −50 • C two tendencies characterise the material behaviour. The loss factor (tan(δ)) is greater by nearly one order of magnitude, and at the same time the storage modulus (E') is higher, representing stiffer behaviour.
It can be concluded that the hysteretic friction is specified by the joint effect of the magnitude of the energy dissipation and the excited volume of the rubber on what the energy dissipation has an effect. In this manner, not only the loss factor has an impact on the hysteretic friction of the rubber, but also the storage modulus. In order to predict the hysteretic friction more accurately or to provide below and/or upper limit for the hysteretic coefficient of friction one has to construct a material model that describes both the storage modulus and the loss factor master curves measured with sufficient accuracy.2D FE model has been developed which is able to predict the contribution of any harmonic component of the measured surface roughness to the hysteretic friction force. Comparing the different wavelength values the hysteretic friction is higher if the smaller one is considered and we have the highest friction coefficient values if the smaller wavelength was superposed to the larger one. These tendencies are remarkable at room and higher temperatures.