FE simulation of the hysteretic friction considering the surface topography

Present study investigates the hysteretic rubber friction in case of ideally rigid, rough surfaces sliding on a rubber plate. To characterize the topography of the rough surface stylus and AFM measurements were performed. Twoand threedimensional models were created in order to predict the hysteretic rubber friction by finite element analysis. The characteristic wavelengths and amplitudes of the engineering surface studied were determined and simple, two-dimensional roughness models with different length scales were constructed. In 3D, NURBS surface was used to model the surface roughness. In order to describe the non-linear and viscoelastic material behaviour of the observed EPDM rubber, two-parameter Mooney-Rivlin material law combined with 40-term generalized Maxwell-model has been used. The effect of sliding speed, temperature and geometric parameters of the roughness model were also investigated.


Introduction
Sliding friction of rubber is involved in many technical applications such as seals, wiper blades, tires etc.The frictional resistance occurring when an elastomer part is pressed and rubbed against a rough counter-surface is usually explained by the joint effect of the adhesion, abrasion (crack propagation), hysteresis and other physical processes.Present study focuses on the hysteretic contribution only.
During deformation of elastomer materials, due to the viscoelastic material behaviour, a part of the strain energy is transformed to heat as a result of the hysteresis [1].The amount of this energy loss depends on the loss factor (tan(δ)) -which is the ratio of the loss and the storage moduli of the material -and on the excited volume [1,2].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 [1,2].The reliability of tribological models 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 technique to characterize the local features of topographies based on the identification of asperities and scratches, while the other is the "global" surface characterization method using complex mathematical tools, such as power spectral density (PSD) analysis.Beside 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.
Generally the studies on friction phenomena are focusing on the measuring of the complex friction behaviour.However, it is difficult to distinguish the hysteretic part of the rubber friction in experimental way, because the other sources of the friction (adhesion and abrasion in dry case, fluid film shearing and boundary lubrication in lubricated case) can not be separated during the measurement.
It is important to know the factors, which are affecting the hysteretic friction, because hysteresis can cause fatigue wear and due to the generated heat, it changes the mechanical and tribological properties of the rubber [3].
In the literature several papers study the hysteretic friction experimentally or in analytical way.Grosch [4] found that there are two sources of the friction when a rubber slides on a dry, rough counter-surface, one is the adhesion, the other is the dissipation inside the rubber, the hysteresis.
In the papers of Persson [1] and Klüppel [2] the authors investigated the energy dissipation via the internal friction of the rubber sliding on a hard, rough surface.Hysteretic friction was determined and 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).
Our aim in this paper is to determine the hysteresis part of the friction when a rubber part is sliding on a rough metal surface, and to investigate the effect of temperature and sliding velocity using FE technique.

Characterization of the surface roughness
Surface roughness measurements were done using a Mahr Perthometer Concept type stylus instrument and an Atomic Force Microscope (AFM) on the plunger of a brake cylinder.Measured topographies can be seen in Fig. 1.Three different methods were used to characterise topographies and find dominant wavelength and amplitude pairs: 3D parameter based technique, slicing method and power spectral density analysis.Details of the methods can be found in [5]- [7].
The 3D surface parameters can be seen in Table 1.The cause of the differences in the parameters is the different size of the measured area.Based on Fig. 1 it can be observed, that the AFM result contains only some dominant topographic elements, so the parameter S z of the AFM measurement (1.888 µm) is in correlation with the amplitude of dominant elements.Of course -as S z of the stylus measurement shows -higher elements also appear on the surface.
It is important that the average slope is in the same range in both measurements.Definition of S dq summarizes the average slope in direction x and y. x and y are the horizontal coordinates while z is the height coordinate.Because the plunger surface has no well defined orientation, it has similar slope in both directions.Thus S dq measured should be halved to get real results, in the sliding direction.So the average slope in the sliding direction is in the range of 1.8-2.4˚.
The asperity analysis can be carried out only in case of stylus measurements, because the AFM topography does not contain enough asperities to make statistical analysis.In the current analysis it was supposed that the area of an asperity is larger than 100 µm 2 .Altogether 65 asperities were found on the scanned area.The height distribution of the peak points, the peak angle distribution in the direction of minor axis and the distribution of the ratio of major and minor axis can be seen in Fig. 2. Most of the peak points are in the range of 0.3 -0.6 µm, and their peak angle is 177˚, which means that the slope is about 1.5˚.The peak value of the axis ratio is 1.7, so if we suppose that every asperity has the same size, and their area is a rectangle with a ratio of edges of 1.7, the distance among the peak points is 95 µm in one direction, and 161 µm in the other.
The PSD of micro-and nano-surfaces can be seen in Figs. 3  and 4.Many local maximum points can be found, that means the surface studied contains many wavelength components.It is important that the largest length scale (1900 nm) in the nano-scale result extends to the micro scale, so the connection is insured.
The surface characterization results are as follows: • The dominant wavelengths (λ) are in the range of 400 nm -250 µm.
• The average slope of the surface (α) is about 1.5-5 o in both the micro-and the nano-scale.
• Amplitude of the different wavelengths (A) can be calculated as:

Material model
The investigated elastomer material was an EPDM rubber filled with carbon black.To describe the material behaviour of the rubber, Mooney-Rivlin material law having two parameters combined with a 40-term generalized Maxwell model was used.Parameters of the Mooney-Rivlin material law were determined by the equations proposed in [8], so the two parameters were c 10 =406.66MPa and c 01 =101.66MPa.These parameters define the glassy modulus of the rubber.
Parameters of the viscoelastic material model were determined based on a DMTA test by applying the method presented in [9].Results of the DMTA tests are storage/loss modulus vs. frequency isotherms which make the construction of master curve at different reference temperatures possible by shifting them horizontally, according to the time-temperature equivalence principle.The master curve is able to describe the behaviour of the material in a broad frequency range at a given reference temperature.As a next step, based on the master curve, a 40-term generalized Maxwell-model was constructed by fitting the material model to the storage modulus master curve using the software ViscoData [10].Following this, the Maxwellparameters were manually adjusted to the tan(δ) master curve, using a trial-and-error technique, to improve the agreement between the measured and the simulated tan(δ) curve.For more detail see [11].Simulations were carried out at temperatures of T= -50, 25, and 150 o C. Temperature-dependency of the material behaviour was taken into consideration by constructing master curves at these temperatures.The simulated storage modulus and the simulated loss factor vs. frequency master curves fitted to the loss factor at three different temperatures are shown in Fig. 6 and Fig. 7.

FE models
To predict the hysteretic friction between an ideally smooth rubber and a rough, rigid counter part, series of FE analysis were performed.

2D FE models with regular surface topography
Analyses were done to study the hysteresis component of rubber friction when a rigid surface with a sinusoidal roughness profile is pressed against a rubber block with pressure (p) and sliding on it with sliding speed (v).The surface profile had a wavelength of λ=200 µm and amplitude of A=1 µm, which is in agreement with the measured surface roughness.
The rubber block with a height of 2 mm and width of 4 mm was modelled with 4000, 2D plain strain finite elements where the nodes on the bottom of the rubber block were fixed.Between the contacting surfaces there was no prescribed coefficient of friction.

Effect of normal load
The rigid surface was pressed into the rubber block with three different pressure values: p=0.05, 0.5 and 5 MPa.The sliding speed was set to v=1 mm/s.Table ?? shows the calculated values of the coefficient of friction as a function of pressure.The In this set of calculations, the configuration described in the previous chapter was used with a pressure of p=0.5 MPa.The effect of the roughness parameters was studied by using sinusoidal roughness profiles with different wavelengths and amplitudes.The basic sine wave, as it was mentioned above, had a wavelength of λ=200 µm and amplitude of A=1 µm.In the next step, this sine wave was replaced with one that had a wavelength of 200 µm and an amplitude of 2 µm, and then one with a wavelength of 400 µm and an amplitude of 1 µm (referred as 1×200, 2×200 and 1×400, respectively).The sliding speed was v=1 mm/s in each case.
The calculated values of the coefficient of friction as the functions of time can be seen in Fig. 7.One can see that the two times larger amplitude produces about 5 times higher hysteretic friction, while the effect of the wavelength is not so significant.

Effect of the combination of different length scales
In case of roughness modelling, the question of length scale is unavoidable.Roughness profiles are often modelled with a set 2D FE models were created to study the effect of roughness models with different scales.Three different roughness models were studied as shown in Fig. 8.
A 0.1 mm high and 0.2 mm wide section of the rubber part was modelled using 4000 2D plain strain elements.In the first case, one sine wave (referred as 'A') with a wavelength of 200 µm and amplitude of 1 µm and without smaller length scale asperities formed the counter-surface.In all the FE simulations performed in this section, the counter-surface was pressed against the rubber block with a constant penetration depth of 0.002 mm, and moved with a sliding speed of 0.1 mm/s.In the next step, the sine wave with a wavelength of 20 µm and amplitude of 0.1 µm was superposed on the sine wave 'A' (referred as 'B').Then, in the third one (referred as 'C'), a sine wave with a wavelength of 2 µm and amplitude of 0.01 µm was superposed on counter-surface 'B' (Fig. 8).The real area of contact can be observed based on the FE contact pressure distribution in Fig. 9.
It can be seen that the rubber is not able to fill out the valleys of the rough counter-surface.The equivalent von Mises stress plots can be seen in Fig. 10.The values of the real area of contact were also investigated.The ratios of the real and the nominal contact area and the predicted hysteretic coefficient of friction can be seen in Fig. 11.Note that the calculated area of contact in the case of roughness model 'A', equals the nominal area of contact.In fact this was considered as the nominal area of contact for cases 'B' and 'C'.It can be concluded that the smaller wavelength and amplitude of the sine waves that are taken into account produces smaller real area of contact, while the same penetration depth was applied.
It can be seen that, in the case of the 'B', the hysteretic coefficient of friction is more than two times higher than in the case 'A'.In case 'C', the value of coefficient of friction is even higher, but the increase, in this case, is not considerable.

Effect of the temperature
The configuration described in chapter 4.1.3was used with the sine wave 'A'.The material model properties were changed in order to model the effect of the temperature on the hysteretic friction.At 150˚C, the Mooney-Rivlin parameters are the same as at 25˚C (see section 3).Only the relaxation times of the Maxwell-model are changed according to the time-temperature superposition.At -50˚C, both the Mooney-Rivlin parameters and the relaxation times are changed relative to the ones applied at 25˚C in order to obtain non-zero loss factor above 1E+4 Hz.
The curves of the coefficient of friction as a function of time can be seen in Fig. 12.The results show that, in case of 150 ˚C, the values of the coefficient of friction are slightly higher than at 25˚C.In the case of the lowest temperature, the coefficient of friction is more than three times higher than at room temperature.This can be explained by the fact that, at the excitation frequency (0.5 Hz), the loss factor is much higher at -50˚C than at the other two temperatures, while the excited volume is the same.The excitation frequency can be evaluated as the ratio of the width of the contact area and the sliding speed.

3D FE models with measured surface topography
In order to investigate the effect of the sliding velocity and temperature on the hysteretic friction when a rubber block sliding on a hard rough substrate a 3D FE model was constructed.
The microscopic surface roughness of the hard counter-surface was modelled using a rigid surface in the FE model, which was described by a NURBS surface, where the measured microtopography (Fig. 1.a) was taken, and 100×100 points were selected from sampling area.The distance among the points was 10 µm.
To reduce the CPU time only a small segment of the rubber block (30 µm × 500 µm × 1000 µm) was modelled (Fig. 13), which was set up by using 8 node incompressible Hermann type elements [8].The rigid counter-surface slides in negative zdirection.Lateral walls of the rubber block that are parallel to the x − y plane were fixed in direction z.The upper plane of the rubber block was also fixed, i.e. the nodes in this plane were constrained in x-, y-and z-direction.The modelled rubber segment was high enough (1000 µm) to avoid the effect of the fixed upper part to the contact zone.To reduce the problem size glued mesh was used in the model.The average element size in the contact zone was 2 µm.The sliding contact between the counter-surface and the rubber was modelled as follows.Firstly, the rigid surface was pressed with a pressure of p=1 MPa against the rubber block using incremental technique meanwhile the relative tangential velocity between contacting bodies was zero.The second step was the horizontal motion with the specified sliding speeds (v=1, 10 and 100 mm/s).Time curves of the pressure load and the horizontal motion are shown in Fig. 14.
Fig. 15 shows the contact state and the distribution of the normal stress σ y at three different temperatures.The areas in contact are visualised by white colour while the black areas are not in contact.
It can be seen in the figures that as the temperature rises the contact area becomes larger and larger.Comparing the normal stress and contact state figures one can estimate which part of the rough surface is penetrating into the rubber block.Where the normal stresses have high negative values, the asperities are The sliding speed has a smaller effect on the contact area but it can be seen that with increasing sliding speed the ratio of the real and nominal contact area is slightly decreasing.It can be seen in Fig. 16 that as the temperature is raising the rough surface can penetrate deeper and deeper in the rubber and the contact area becomes larger.The reason for the higher penetration depth at higher temperatures is the lower storage modulus.It is important to characterize the penetration depth, which is proportional to It can be established on the basis of Fig. 17 that the highest hysteretic frictional coefficient is 0.006 obtained at T = −50 o C and v=1 mm/s.The reason can be that, the rubber can not fill out the valleys of the rough surface at higher sliding velocities and there through the excited volume becomes smaller.

Conclusions
2D FE models have been developed in order to predict the hysteretic friction contribution of different harmonic components of the measured surface roughness from nano-to microscale.3D FE model has also been evolved to predict the real contact area and the hysteretic friction force by taking into account the real, measured surface roughness.
The hysteretic component of the friction is affected by many parameters, such as ratio of the amplitude and wavelength, slid- ing speed, temperature, normal load etc.It was shown that in case of roughness modeling the combination of different length scales is an essential issue.Considering the correct scale(s) is inevitable in order to get the right results.One can see that the apparent or nominal contact area can be much larger than the real contact area in the case of superposing the sine waves of the different length scales.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.
It can be also concluded that the calculated hysteretic friction is very small in case of rough surfaces used in engineering applications, but for example it can be higher in case of vehicle tyres when they are rolling/sliding on the rough asphalt.

Fig. 7 .
Fig. 7. Calculated hysteretic friction vs. time curves in case of different amplitudes and wavelengths

Fig. 9 .Fig. 10 .
Fig. 9. Vertical normal stress (σ y ) distribution in case of roughness model A, B and C

Fig. 11 .
Fig. 11.Ratios of the real and the nominal contact area and the predicted hysteretic coefficients of friction for the surface roughness model A, B and C

Fig. 12 .
Fig. 12. Coefficient of friction vs. time curves at three different temperatures (counter-surface A)