FE Prediction of Hysteretic Component of Rubber Friction

The hysteretic part of the friction coefficient for rubber sliding on an ideal rigid, rough surface has been investigated by FE technique. The FE models were created by using two different FE softwares, ABAQUS and MSC.MARC. The surface roughness has been considered by using two different sine waves having a wavelength of 100 μm and 11.11 μm, as well as their superposition. Parameters of the viscoelastic material models of the rubber were gained, firstly from a fit to the measured storage modulus, secondly from a fit to the measured loss factor master curve of the rubber. The effect of viscoelastic material models, comparing 10-term and 40-term generalized Maxwell models was also considered together with the temperature effect between −50 and 150◦C. According to the results, both postprocessing methods, namely, the reaction force and the energy-based approach, show very similar coefficients of friction. The 40-term Maxwell model fitted to both the storage modulus and loss factor curve provided the most realistic results. The tendency of the FE results has been explained by semianalytical theory.


Introduction
It is important to characterize the friction behaviour of rubber components sliding on rough counterpart in a number of applications, such as seals, wiper blades, and tires. Many papers found in the literature deal with improving the tribological properties of sliding pairs and decreasing the friction. In the absence of lubricant when contacting surfaces are dry and clean, the rubber friction is mainly due to adhesion (especially at low sliding speeds) and hysteresis [1][2][3][4]. The hysteretic friction comes into being when the rubber is subjected to cyclic deformation by the macroand/or microgeometry (surface roughness) of the hard, rough substrate [2,4]. When a rubber component slides on a hard, rough substrate, the surface asperities of the substrate exert oscillating forces on the rubber surface leading to energy "dissipation" via the internal friction of the rubber. In the most engineering applications, rubber/metal sliding pairs are lubricated in order to decrease the friction force arising in dry case. In presence of lubricant, rubber friction is due to hysteretic losses in the rubber, boundary lubrication, and fluid friction. The lubrication decreases the contribution of nanoroughness to hysteretic friction because lubricant fills out the nanovalleys, which makes the penetration of the rubber impossible into these regions. In the case of fluid friction, friction force comes from the shearing of a continuous, relatively thick fluid film. At the same time, in regions where a very thin lubricating film separates contacting surfaces, the friction force comes primarily from the shearing of this thin boundary layer. In a recent paper, Kozma [5] explains the friction force measured in lubricated rubber/metal sliding pair with the boundary lubrication theory.
Many articles in the literature deal with the sliding behaviour of rubber in order to develop models for the prediction of hysteretic friction force. Pioneering work of Grosch [1] has shown that rubber friction arises, first of all, from interfacial adhesion and hysteresis losses in the rubber itself. Furthermore, it has been found that both adhesion and deformation losses are directly related to the viscoelastic properties of the rubber.
The hysteretic friction force is caused by the macrogeometry when a hard, smooth ball rolls or slides on a rubber substrate. In case of rolling, hysteresis in the rubber is the primary source of the frictional resistance. Felhos et al. [6] have studied the rolling resistance to be arisen when a steel ball rolls on a rubber plate. Besides the experimental measurements, the rolling process was also simulated by finite element method. The viscoelastic material behavior of the rubber was described by a 15-term generalized Maxwellmodel fitted to the results of dynamic mechanical thermal analysis (DMTA). Results of the experimental tests and of the simulation were compared, and a good agreement was found between them.
In [7], authors have investigated the friction behaviour in the case of steel ball sliding on rubber plate by FE technique. 40-term Maxwell-model was used in their study at three different temperatures. The hysteretic friction predicted was in good agreement with the experimental results obtained at room temperature. In lubricated case, the effect of mixed friction was taken into account by a prescribed coefficient of friction. [8] has investigated hysteretic friction at asperity level in the case of unlubricated steel/rubber sliding pair.
When there is no hysteretic friction caused by the macrogeometry, hysteretic friction is due to the surface microtopography (surface-roughness) of the hard counterpart. Persson et al. [2] has developed a friction model motivated by rubber physics to predict the surface roughnessinduced hysteretic friction for rubber/road track sliding pairs. Although all length scales of the surface roughness contribute to the hysteretic friction, Persson [3] concludes that rubber friction on rough surfaces is mainly due to the viscoelastic deformations of rubber on length scales larger than 1 μm. In [4], Klüppel and Heinrich presented a theory for rubber friction on rough surfaces, which relates friction force to the dissipated energy of the rubber. The model presented was in agreement with the measurements of Grosch [1]. In [9], an experimental study has been performed to evaluate the predictive capability of the models of Persson [2] and Klüppel/Heinrich [4]. Sliding friction coefficient of rubber/road asphalt contact pairs was measured under dry conditions and at a contact pressure of 0.4 MPa. The measured data were then compared against the predictions of the two friction theories. The main conclusion was that, above 10 mm/s, both models overestimate the measured friction force. As a possible cause, the flash temperature effect not included in the models was mentioned. The authors in [10] have applied the Persson model for lubricated metal/rubber sliding pairs, where the metal counterpart has an apparently smooth surface. The theory applied emphasizes that the whole length scale of the surface roughness (from nano-to microlevel) has to be taken into account at the prediction of the hysteretic friction force. Considering the coefficient of friction values presented, it must be noted that the theory resulted in surprisingly high friction force even if the countersurface is apparently smooth. The measured friction force has been attributed to the hysteretic losses only although, very likely, the insufficient lubrication (boundary lubrication) and the rubber wear contributed also to the high coefficients of friction measured.
The aim of this study is to predict the microtopographyinduced hysteretic component of rubber friction by using the finite element method. In order to increase the confidence in the predicted results, calculations were performed by two different commercial FE softwares (ABAQUS and MSC.Marc). The FE-based modeling approach used here is the same as in [6] where the hysteretic friction was predicted on macrolevel by using a steel sphere and a rubber plate. In [6], a very good agreement was found between the FE prediction and the measurement, which proved the applicability of the FE modeling approach used here. Through the FE modeling of hysteretic friction of apparently smooth surfaces, authors take an attempt to investigate (1) how an additional harmonic component of the surface roughness superimposed to a larger one influences the hysteretic friction, (2) how it is worth to determine the Maxwell parameters for the description of viscoelastic nature of the rubber at different temperatures, (3) the relation between the reaction force based and the energy-based derivation of the friction force.

Description of the Problem Investigated
In order to examine the hysteretic portion of the coefficient of friction when smooth rubber block slides on a hard, rough surface, two-dimensional FE models were constructed to make the interpretation of results easier. The microscopic surface roughness of the hard countersurface was modelled by two different sine waves having a wavelength of 100 μm and 11.11 μm (countersurface A and B) and by their combination (countersurface A+B) where the countersurface B was superimposed on countersurface A. Figure 1 shows the three different surface roughness models used in the present study. In order to emphasize the difference among the models, we depicted them in a single graph by using different scale in the horizontal and vertical directions. Horizontal and vertical peak-to-peak distances that is, the wavelength (λ) and the double of the amplitude (2 · A) of the countersurface A and B can be seen in Table 1. The excitation frequency with which the asperities excite a surface layer of the rubber can be estimated as the ratio of the sliding speed and the wavelength. The sliding speed was considered to be 10 mm/s. 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 micro-to nanolevel) to the hysteretic friction.

FE Model
Due to the periodicity in the surface roughness, it is sufficient to model a small, repetitive segment of the rubber. Naturally,   for such a small model, the condition of repetition has to be applied as boundary condition. Due to the repetitive symmetry, the length of the 2D FE model in sliding direction was the characteristic wavelength of the surface roughness models. According to this, the length of the FE models for countersurface A, B, and A+B was 100 μm, 11.11 μm, and 100 μm, respectively. The length of the rigid countersurface is larger than the wavelength of the irregularity. Its length was specified to insure continuous contact with the rubber during the simulation. It means that the rigid countersurface is composed of identical irregularities having a wavelength of 100 or 11.11 μm, thus many irregularities make a contact with the rubber during the simulation. The modelled rubber segment had a large enough height (500 μm) to ensure that the load application zone (upper part in Figure 2), which was defined to remain perfectly horizontal during the simulation, is not able to affect the contact zone (lower part in Figure 2). In the contact zone between the rubber and the rigid counterpart, the FE mesh consisted of elements having a size of 0.25 μm in x-direction. The element size of 0.25 μm was obtained from sensitivity analysis. The rubber has been modelled as a deformable large strain viscoelastic material, being treated as a plane strain solid consisting of CPE4H elements [11] and QUAD80 elements [12]. In Abaqus, the hard, rough countersurface has been modelled by means of "stiff " elements, for which plane strain quadratic solid elements (CPE8R), with mechanical properties corresponding to a linear elastic material several orders of magnitude stiffer than rubber, have been used. In respect of Abaqus, authors used version 6.3 (released in 2002) because newer version was not available at the institutions  involved in the current research. According to the authors' experiences, such a modelling strategy approximates the solution, provided by a rigid body approach, very accurately. In MSC.Marc, the rough countersurface has been modelled as an analytical rigid surface. The sliding contact between the rubber and the countersurface was modelled as follows. Firstly, the rubber segment was compressed with a pressure of p = 1 MPa against the countersurface by using incremental techniques. During the indentation, the relative tangential velocity between contacting bodies was zero (see Figure 3). The second step was the acceleration of the rubber block up to a sliding speed of 10 mm/s. During this, the countersurface was fixed in vertical direction and the nodes on the lateral walls of the rubber were forced to move identically (see Figure 2). 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 steps.
In the course of steps 2 and 3 (acceleration and motion with constant speed), nodes on the lateral walls can move horizontally in order to allow shear of the model, but both in the same way to assure that the small segment represents a finite rubber sample that is composed of identical sections. With this purpose, constraint equations have been defined for the corresponding nodes on both lateral walls (see Figure 2), which imposes the same horizontal displacement upon the node and its equivalent on the contrary lateral wall. In the same way, the same condition has been imposed 4 Advances in Tribology for the vertical displacement of the nodes at both lateral walls. Figure 2 shows the deformed shape of the FE model, according to the boundary conditions discussed above. Time curves of the pressure load and the horizontal motion of the rubber segment are shown in Figure 3.

Viscoelastic Material Model.
To model the mechanical behaviour of structural elements made of rubber and rubberlike materials, first of all, an appropriate material model is needed, which is able to describe their static as well as dynamic material behaviour simultaneously. Thus the description of the material model can be divided into two parts: one is responsible for the nonlinear σ-ε function, and another is responsible for the time-and temperaturedependent material properties. To model the nonlinear stress-strain curve of the rubber, the commonly used Neo-Hooke and Mooney-Rivlin material laws were employed (nonlinear springs), while to model the time-dependency of the material behaviour an n-term Maxwell-model was used. The short-time (instantaneous) behavior of the generalized Maxwell-model applied is specified by an energy density function (W 0 ) that is given, in our case, by a two-term Mooney-Rivlin material model. Thus our largestrain viscoelastic model shows nonlinearity between stresses and strains that, at the same time, depends on time. The instantaneous strain energy density is distributed among the branches of our spring-dashpot model by assuming that the instantaneous strain energy in the ith Maxwell element can be calculated as where e i is the nondimensional energy parameter of the ith Maxwell element and W 0 is the instantaneous strain energy defined by the Mooney-Rivlin model. The interrelation among the energy parameters can be written as where e ∞ is the energy parameter of the spring-that determines the relaxed response of the rubber-connected parallel to the Maxwell elements and N is the number of the Maxwell branches. In other words, the material behaviour of the rubber was modelled by a generalized Maxwellmodel ( Figure 4) in which the nonlinear spring and the n-term Maxwell-model are assembled in a parallel form.
The generalized Maxwell-model can be a reasonable choice as it describes the measured material behaviour (see later) properly and is available in the most commercial finite element software packages. Parameters of the material model were determined from a fit to DMTA measurements by applying the method presented in [6]. For more details on the nonlinear, viscoelastic material model used, see [12][13][14][15].
Results of the DMTA measurements are storage/loss modulus versus frequency isotherms (from T = −90 • C to T = +150 • C with a temperature increment of 5 • C), which makes the construction of master curve at different reference temperatures possible. In the course of the DMTA test, excitation frequency varied between 1 and 100 Hz. The master curve considering a reference temperature describes the behaviour of the material in a broad frequency range at the given reference temperature. According to the timetemperature superposition principle, the master curve can be constructed by horizontal shifting of the measured isotherms. Naturally, the isotherm belonging to the reference temperature remains unchanged and only the rest of the isotherms are shifted horizontally.
After constructing the master curve, the next step was to create the generalized Maxwell-model for the finite element models. Parameters of the generalized Maxwell-model were determined by fitting the material model to the storage Advances in Tribology 5 Non-linear spring modulus master curve. Curve fitting was performed by software ViscoData [16]. Two variants of an EPDM rubber (EPDM v1 and EPDM v2) containing different amount of carbon black were investigated in this study. In the case of EPDM v1, the parameters of a 10-term generalized Maxwell-model were gained from a fit to the storage modulus master curve constructed at temperature T = −45 • C. The glass transition temperature of the EPDM rubber studied was −46.3 • C. The parameter identification can be performed by different ways, for example, by the software Abaqus itself (see the section "Calibration the Prony series parameters" in [11]). At this point, it must be noted that most of the commercial FE software have limited capability for "calibrating" n-term Maxwell model. In most cases, the number of terms is limited anywhere between ten and fifteen. The reason why a 10-term Maxwell-model was used is that such a material model is able to describe the variation of E in a wide frequency range and can be calibrated by without over passing the capability of the applied, commercial FE software. The research work summarized here has been carried out by two research groups. One of these groups used Abaqus, while the other used MSC.Marc because these softwares were available at the research institutes. The group using Abaqus did not investigate the effect of Maxwell-terms on the predicted hysteretic friction. They concentrated their attention on the comparison of energy-based and reaction force-based computation of hysteretic friction only. The other group focused its attention on the effect of Maxwell parameters on the predicted coefficient of friction and the development of a semianalytical model. Naturally, both the Abaqus and MSC.Marc are able to handle high-order Maxwell-models in the computations; however, they have strongly limited capability for model calibration. This is the reason why the software ViscoData was used at the calibration of 40term Maxwell-models. As non-linear spring, the Neo-Hooke material model with a constant of c 10 = 0.968 MPa (E ∞ = 9 MPa, where E ∞ is the relaxed modulus of the rubber) was used in the generalized Maxwell-model. In all the simulations with EPDMv1, the temperature dependency of the material behaviour was defined by the WLF transformation using universal constants C 1 = 17.4, 1E − 08 1E − 04 1E + 00 1E + 04 1E + 08 1E + 12 1E + 16 This transformation evaluates the material behaviour at a temperature other than the reference temperature by shifting the master curve constructed at the reference temperature horizontally. Consequently, this transformation influences the relaxation times of the material model only. Figures 5 and 6 show the E versus frequency and the tan(δ) versus frequency curves of the 10-term generalized Maxwell model at different temperatures. The storage modulus curve (E ) shows a physically realistic behavior, while the tan(δ) curve provides unrealistic behaviour due to the high oscillation. Accurate numerical prediction for the hysteretic friction cannot be anticipated as long as the material model is not able to describe with satisfactory accuracy both the storage modulus and the loss factor curve. Since surface roughness can be approximated as the superposition of sine waves with different amplitudes and wavelengths, this has to be insured in a sufficiently wide frequency range that is called interesting frequency range. Additionally, any variation in temperature shifts the frequency-dependent material behaviour towards higher or lower frequencies, which results in situations where different parts of the material model curves fall into the interesting frequency range. Therefore, the parameter identification of the material model has to be performed carefully.
Basically, both the loss factor (tan(δ)) and the storage modulus (E ) influence the hysteretic component of rubber friction. The former affects directly the magnitude of the energy dissipation, while the latter has an effect on the excited volume. The actual value of the hysteretic friction is determined by their joint effect. In order to reduce the oscillation of the fitted loss factor curve, the number of terms has to be increased [13]. This is why, in the case of EPDM v2, more sophisticated 40-term Maxwell models were constructed.
By using the common practice, as a first step, the Maxwell parameters were determined from a fit to the storage modulus master curve. In this case, the agreement between  the measured frequency-dependent storage modulus master curve and the Maxwell model is very good (see Figure 7). Contrary to this, the Maxwell model provides much smaller tan(δ) values in the interesting frequency range than the measurement (see Figure 8). Consequently, FE prediction using this material model will underestimate the hysteretic friction force. The reason why authors performed manual modification is that the software ViscoData is not able to fit the Maxwell model to the measured tan(δ) master curve or the measured E and E curves simultaneously. To improve the agreement with the tan(δ) master curve, as a second step, the Maxwell parameters were manually adjusted using a trialand-error technique [13]. The aim of this modification is to obtain parameters with which both the storage modulus and the loss factor curves can be described with acceptable accuracy in the interesting frequency range that is, close to the characteristic excitation frequencies. This method can be regarded as an optimization of the material model for the measured material behaviour. Naturally, such a manual modification of the Maxwell-parameters influences not only the tan(δ) but also the E curve of the material model. In other words, such a modification results in better agreement between the frequency-dependent tan(δ) curves and weaker agreement between frequency-dependent E curves of the measurement and the Maxwell model than in case of fitting to the measured storage modulus curve only (see Figures 7  and 8).
For the sake of clarity, this means that the estimation of the excited volume is more realistic using Maxwell model fitted to the storage modulus, while the estimation of the energy dissipation is more realistic in the case of modified Maxwell model (Maxwell model fitted to tan(δ)).
In case of EPDM v2, the Maxwell-model was combined with a two-parameter Mooney-Rivlin model that corresponds to a non-linear spring. Its parameters were as follows: c 10     master curve. The Mooney-Rivlin parameters change due to the manual modification. The nondimensional energy parameters (e i ) of the Maxwell model fitted to the loss factor (tan(δ)) master curve were obtained by manual modification from the energy parameters of the Maxwell model fitted to the storage modulus master curve. Due to this modification, the sum of the energy parameters of the Maxwell elements ( 40 i=1 e i ) has been changed. Consequently, if we used the same glassy Mooney parameters in both cases, we would obtain an accurate relaxed modulus at the Maxwell model fitted to E , while at the Maxwell model fitted to tan(δ), the relaxed modulus would be considerably underestimated in the interesting frequency range. To minimize this underestimation that is, to approach E the real values, the glassy Mooney parameters have been changed. Thus we obtain higher glassy modulus than the measured one but the agreement, within the interesting frequency range, between the modeled and measured E values will be better than in the case of same Mooney parameters.
1E − 08 1E − 04 1E + 00 1E + 04 1E + 08 1E + 16 1E + 12  FE models with 40-term Maxwell-model were solved by MSC.Marc only. According to [12], the Mooney-Rivlin parameters specified correspond to the glassy modulus of the rubber. Temperature-dependency of the material behaviour was taken into consideration by constructing master curve at temperatures of T = −50, 25, and 150 • C. Figures 9 and 10 show frequency-dependent E and tan(δ) curves of 40-term Maxwell models fitted to the storage modulus as well as to the loss factor master curve at three different temperatures. In Figure 9, one can observe the difference in storage modulus between the two different fitting techniques. In case of fitting to the loss factor at 25 and 150 • C, the storage modulus is underestimated below a certain frequency, while it is overestimated above that. At −50 • C, the interesting frequency range falls right from the maximum of the tan(δ) curve; therefore, the Maxwell parameters differ from the ones applied at 25 and 150 • C in case of fitting to the measured loss factor curve. Due to this, at −50 • C, the storage modulus curves of the two types of fitting run together up to 10,000 Hz, while it becomes overestimated above this frequency.
One can also see some difference in case of loss factor curves in Figure 10. It can be concluded that at 25 • C under 10,000 Hz and at 150 • C under 1E + 8 Hz, the loss factor becomes underestimated in case of material model fitted to the storage modulus. At −50 • C, the loss factor becomes also underestimated above 10,000 Hz using fitting to the storage modulus master curve.
It can be concluded that, within the observed frequency range, the material model fitted to the measured loss factor master curve describes more accurately the measured behaviour.

Calculation Based on Reaction
Forces. The simplest way to determine the hysteretic friction coefficient is based on the horizontal and vertical reaction forces μ Force . Since the magnitude of friction force can be calculated as the sum of horizontal reaction forces, the coefficient of friction was defined as the ratio of the total horizontal and vertical reaction forces:

Calculation Based on Energies.
The hysteretic friction force can also be calculated from the dissipated energy. If a zero coefficient of friction is prescribed at the contact interface μ interface , then the only source of energy dissipation is the hysteresis in the rubber. Thus, the dissipated energy can be interpreted as the work of the hysteretic friction force. If the contribution of lower-scale roughness or physical processes, other than hysteresis (adhesion, boundary friction, etc.), to the friction force has to be included in the present model, then a nonzero coefficient of friction at the contact interface (μ interface = 0) is needed. In that case, energy dissipation is due to hysteresis and interfacial friction. According to this, it has been considered that there are mainly three energetic contributions which are assumed to be involved in the friction process: elastic energy (ALLSE-coming from elastic deformations), viscous dissipation (ALLCD-coming from viscoelastic dissipation in the rubber), and "frictional" energy (ALLFD-coming from the friction coefficient value included in the FE model). The contribution to the total friction coefficient from each energy can be calculated and it is assumed that, in the present FE model, the sum of coefficients of friction coming from these three energetic  contributions corresponds approximately to the total friction coefficient value μ Energy :

FE Results
In the followings, authors focus their attention on the interpretation of results obtained during motion with constant speed (step 3). Table 2 shows the average steady-state coefficients of friction (COF) for temperatures of T = −45, −35, and 25 • C. The difference between coefficients of friction μ Force and μ Energy is less than 5% in most of the cases; therefore, both techniques can be applied. As we assumed the same average contact pressure at both length-scales studied (countersurface A and B) the COF calculated for countersurface A+B may not be considered as the sum of those values. Comparing the two FE softwares, it can be concluded that in the most cases the differences in COFs are less than 26%. Due to the oscillation in the tan(δ) curve, the 10-term Maxwell-model is inadequate for COF predictions by FE technique. In order to avoid the high oscillation of the loss factor, the number of Maxwell terms in the material model has to be increased. For this purpose, 40-term Maxwell models were constructed. Table 3 shows the coefficient of frictions coming from the hysteresis in case of three different roughness models and three different temperatures considering the two different 40-term Maxwell models, in the case of zero coefficient of friction at the contact interface (μ interface = 0).

40-Term Generalized Maxwell-Model.
A semianalytical theory has been developed to interpret our FE predictions for the hysteretic component of rubber friction. The semi-analytical theory developed is based on qualitative discussion of hysteretic friction. Similar analyses can be found, among others, in [17,18].
The COF coming from hysteresis can be expressed as the ratio of the friction force and the normal force. The friction force can be calculated from the energy loss during sliding with a distance of λ (wavelength) The energy loss during a single cycle in a unit volume of a prismatic rod under cyclic compression is W unit = σ 0 · ε 0 · π · sin(δ). Thus as an approximation the hysteretic COF can be written as where V excited is the excited rubber volume. Substituting the excited volume by V excited = h · k · λ · 1, where k is a proportionality factor (see Figure 11), which is used to define the depth of the excited layer, and using Hooke's law σ 0 = E * · ε 0 , the COF is The normal force can be calculated from the pressure on the given surface: The correlation between Advances in Tribology 9 2·a 1 λ h Figure 11: Illustration of a sine wave pressed against a rubber block. the complex and the storage modulus is E * = E / cos(δ). By inserting this in (8), the hysteretic COF will be that is By assuming that the strain amplitude is proportional to the ratio (k · h)/λ, the hysteretic COF will be or So, in case of a given countersurface (A or B), μ hysteresis will be proportional to the product of E / p, tan(δ), and (h/λ) 3 . By knowing the storage modulus, the loss factor, and the wavelength of the countersurface, as well as the penetration depth calculated by FE technique the μ 0 (semianalytical prediction) can be calculated for both countersurfaces. These μ 0 values can be used to determine how the hysteretic COFs relate to each other at various temperatures and different countersurfaces. Obviously, not only the hysteretic component of COF can be determined this way, but also the hysteretic friction force as It can be established on the basis of the COFs Table 3 that for both 40-term material models, the COF calculated at countersurface A is the highest at 25 • C while the lowest at −50 • C. On the basis of Figures 9 and 10, it can be concluded that the product E ·tan(δ) decreases as the temperature increases, while the ratio (h/λ) 3 increases by the temperature increase. As a resultant of these two effects, μ hysteresis will first rise as the temperature increases, then, having reached a peak value, it begins to fall. The same tendency can be observed in case of semi-analytical prediction.
In the case of countersurface B, a similar tendency can be observed; however, in case of a modified 40-term Maxwell model, the highest COF can be observed at 150 • C, in accordance with the semi-analytical prediction. Equation (12) shows that hysteretic friction depends not only on E and tan(δ), but it cubically depends on the ratio of penetration depth (h), and the wavelength (λ) of the rough surface as well as on the proportionality factor (k). Obviously, the ratio h/λ is not identical in case of countersurfaces A and B.
In case of countersurface A+B, it can be established for all three temperatures and all three material models that hysteresis will be greater than at countersurfaces A or B separately.
Also in case of surface A+B, the greatest hysteresis loss can be calculated at 25 • C. At −50 • C, nearly identical hysteretic COF can be calculated in case of 40-term material models. At 25 • C and 150 • C, the greatest COF is generated at countersurface A+B, followed by countersurface B and then by countersurface A.
As regards countersurfaces A and B and the two 40term material models, Figure 12 shows the COF values normalized by the highest COF in case of FEM and semianalytical prediction. In case of semi-analytical prediction, the diagrams indicate μ 0 /μ 0 max , while μ Force /μ Force max in case of FEM.
It can be established that-taking the three temperatures into consideration-semi-analytical prediction is not suitable for reproducing numerical results but, in tendency, it follows FEM results. It can be seen in Figures 12(a)-12(d), that the COF proportions calculated by semi-analytical prediction are lower than ones calculated by FEM, except the highest COF.
In order to verify our FE predictions for hysteretic friction rubber sliding, experiments on lubricated or contaminated sinusoidal rough surfaces would be needed. As it is well known even a very thin adhered layer of lubricant or contaminant is able to eliminate the adhesion between contacting surfaces. Hence the hysteretic component of rubber friction can be studied separately from adhesion. In spite of the great importance of microhysteresis, only a few experimental studies devoted to this topic are available in the literature. As authors have no test results for verification purposes, the ones of [19] are used to prove the reliability of FE predictions. Rana et al. [19] have developed an experimental reciprocating rig to study the contact conditions of an elastomeric seal with a hard surface at room temperature. The stroke length, stroke speed, and the load on the rectangular seal were varied. The counter surfaces included smooth glass (R a = 0.0075 μm) and steel having a mean surface roughness of R a = 1.4 μm. The tests were carried out at a mean surface pressure of p ≈ 3 MPa in dry and lubricated conditions. The stroke length was 5 mm, while the reciprocating frequency varied between 1.5 and 11 Hz. The velocity of the slider changed with time in a sinusoidal manner since rotary motion was converted into linear. A reciprocation frequency of 1.5 Hz implied an average stroking speed of 15 mm/s with a peak velocity of 47 mm/s. They concluded that the reciprocating with a rough steel plate increases the friction due to the interaction of hard steel asperities with soft rubber ones. Coefficients of friction computed based on the measured friction force values of [19] can be seen in Figure 13. Measured COFs show clearly that the higher surface roughness of steel causes higher sliding friction due to the increased hysteretic friction. The increase in rubber friction is about 0.09, which is comparable with our predictions for a sine profile having an amplitude of 1 μm. It is worth to mention that as the sliding speed is low neither the hydrodynamic effect nor the frictional heat generation play a role in the measured sliding behavior of rubber.

Conclusion
(1) An 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.
(2) It can be concluded that the 10-term generalized Maxwell-model is not able to describe the complex material behaviour of the rubber accurately due to the high oscillation in the tan(δ) curve. The effect of oscillation in the tan(δ) curve on the hysteretic friction is the most critical at room temperature and at temperatures above that. In order to avoid this oscillation, a Maxwell-model with 40-term was proposed.
(3) It can also be concluded that around the observed frequency the 40-term generalized Maxwell-model fitted to the storage modulus produces significantly smaller COF than the 40-term generalized Maxwell-model fitted to the measured loss factor curve. Therefore, an FE prediction using a Maxwell-model fitted to the storage modulus master curve underestimates the hysteretic component of rubber friction (see Figure 10).
(4) In order to predict the hysteretic friction accurately, the material model applied has to be able to describe with sufficient accuracy both the measured storage modulus and loss factor master curves. A Maxwell-model fitted, as usual, to the measured storage modulus master curve only may Glass, 1.6 Hz, R a = 0.0075 micron Coefficient of friction (-) Figure 13: Coefficients of friction of lubricated rubber/glass and rubber/steel sliding pairs measured at room temperature. (Based on the results of [19]). Oil viscosity = 65.9 mm 2 /s. result in an inaccurate prediction of the hysteretic friction force especially at room temperature and above it.
(5) It is important to mention that not only the loss factor influences the hysteretic component of rubber friction, but the storage modulus too, which has an impact on the excited volume. 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. (6) 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.

Nomenclature λ:
W avelength A: Amplitude u: Displacement in x-direction v: Displacement in y-direction p: Pressure v x : Sliding speed in x direction T: T emperature E i : Elastic modulus of the ith spring η i : Relaxation time of the ith Maxwell-element c 10 , c 01 : Mooney-Rivlin parameters C 1 , C 2 : WLF constants a T : shift factor T ref : reference temperature E : storagemodulus tan(δ): loss factor f : frequency μ Force : hysteretic friction coefficient calculated by the reaction forces of the FE model F friction : friction force F normal : normal force μ interface : prescribed coefficient of friction at the contact interface μ Energy : coefficients of friction calculated by energies μ pred : coefficient of friction predicted by semi-analytical model W: energyloss W unit : energy loss per unit volume V excited : excited volume h: pe netratio nd epth k: proportionality factor E * : complexmodulus σ 0 : stressamplitude ε 0 : strainamplitude V Original : original volume COF: coefficient of friction.