CC BY
69
УДК 539.62
Investigation of the dry normal contact between fractal rough surfaces using the reduction method, comparison to 3D simulations
R. Pohrt, V.L. Popov
Berlin University of Technology, Berlin, 10623, Germany
The reduction method was applied to solving the elastic contact problem of fractal rough surfaces having different fractal dimensions. These surfaces are of particular interest, as contact spots influence each other on large and small scales, depending on the power spectrum. Results were compared to previous boundary-element method simulations that had been conducted with 3D surfaces, taking into account all interdependences between contact spots. Excellent agreement in the contact stiffness was found. Even though no interactions are allowed in the reduced model, the same power-law behavior can be observed as in the 3D case, due to the transformation of the power spectrum. Theoretical considerations about the power-law for different fractal dimensions are given. Keywords: fractal surfaces, contact problems, fractal rough surface
1. Introduction
The roughness of surfaces has a great influence on many physical phenomena, such as friction, wear, sealing, adhesion, as well as electrical and thermal conductivity [1-4]. At the same time, it is one of the main reasons preventing the broad usage of numerical simulation methods in tribo-logy — in contrast to many other contemporary areas in physics and engineering. Bowden and Tabor [2] were the first to realize the importance of the surface roughness of bodies in contact. Because of the roughness, the real contact area between the two bodies is typically orders of magnitude smaller than the apparent contact area. The works of Bowden and Tabor triggered an entirely new line of theory for contact mechanics regarding rough surfaces in the 50's and 60's, with basic contributions by Archard [3] and Greenwood and Williamson [5]. The main result of these examinations was that the contact area between rough elastic surfaces is approximately proportional to the normal force. At the same time, it was realized that in many tribological problems, roughness on different spatial scales plays an important role in determining contact properties, friction, and wear. This understanding led in the last years to extensive studies of contact properties of bodies with rough surfaces [4, 6-8]. The main focus of these studies was put on the determination of the real contact area for self-affine surfaces with relevant fractal dimensions.
2. Direct numerical approaches
Solving the general contact problem is not a trivial thing to do. The fact that surfaces are rough, makes it a discrete problem with regions either being in contact or not and different boundary conditions applying in both cases [9]. As our goal is to include roughness on many length scales, the finite element method is not the best option because grids would have to be very fine in all spatial dimensions. Even when only discretizing the surface itself, as in the frame of the boundary element method, there still is the problem of all grid points influencing each other. One must deal with full (not sparse) matrices which rapidly increase in size [10]. In the right column of Fig. 1, typical shapes of the surface deformation are shown. There is a considerable deformation outside of the contacting areas, showing how these areas influence each other. The use of multigrid techniques, such as the multilevel integration introduced by Brandt and Lub-recht [11] and the specially adapted CG-method by Polonsky and Keer [12], facilitated the solution of this problem dramatically.
3. Scaling and dimensionless variables
In [13], as well as in the present paper, we use a set of dimensionless expressions for the contact stiffness as well as for the normal force. The motivation for these definitions will be elaborated here.
© Pohrt R., Popov V.L., 2012
The maximum stiffness that can be obtained with the finite square sample A0 is when 100% contact area is given. In this case, the stiffness is
ks = 1.012 • 2E y Ao/ n = 1.1419 • E*^.
It only depends on the contact configuration, e.g., the distribution of contact patches and free surface areas. Notably, it does not depend on the stress distribution within the contact, nor on the roughness. Any contact configuration on A0 with less than 100% contact area will lead to k < ks, but it can be found that k approaches ks quickly with an increase in contact area. We can now define a dimensionless stiffness related to the maximum value:
k =-k * --- (1)
1.1419 • E JA
The most remarkable property of the contact stiffness is
its proportionality to the spatial length of the whole system
L = ^/A^, not to the area.
The normal force FN for a given indentation d can be
d
expressed as FN = J k(z) dz. Thus, FN scales in the same
o *
manner as k with respect to A0 and E . In addition, it also
scales linearly with respect to the surface height z. The best choice for a dimensionless force is, thus,
FN =
FN
(2)
E \/A0'
where h is the root mean square roughness of the surface topography.
4. Surface description
The interactions of different contact spots dramatically depend on the surface characteristics. Contact patches may be small but close to each other or they might be concentrated in a specific area and rarely occur further away. The proportion of smaller and greater wavelength amplitudes determines this characteristic. Therefore, it is reasonable to choose fractal rough surfaces for the investigation. These surfaces have a spectral density obeying a power law [7]:
-2( H+1)
C2D(q) = const
%
(3)
The parameter H ranges from 0 to 1 and makes the surface look either wavy or fuzzy. By choosing different values of
Fig. 1. Graphical representations of fractal surfaces in 1D and 3D. For 4 values of the fractal dimension Df = 2,7/3,8/3, 3 (a-d), a typical 3D topography is shown in the second column. One specific value of the dimensionless normal force FN has been applied to these surfaces to give a plot of the real area of contact in the third column. In the last column, one can see a graphical representation of the surface deflection for the specific sample of a contact configuration. For comparison, an equivalent 1D rough line is shown on the left hand side that has been generated according to (7). All illustrations shall be seen as examples, because phases in 3D as well as in 1D are random
H, one can see the transition from dominant long range interactions to shorter, almost equally distributed contact spots. In the latter case only a few adjacent spots are strongly affected by the local pressure. Figure 1 shows some realizations of such surfaces and examples of their associated contact area shapes for a specific value of the normal force. Another advantage of these surfaces is the fact that they have been examined by various researchers in the past [4, 6, 8, 14]. As our goal is to compare 3D results to the 1D line, only isotropic surfaces will be considered.
The authors recently published extensive results for the 3D contact of fractal rough 3D surfaces, including the contact stiffness, as well as the real contact area for a series of fractal dimensions [13]. The basic result was that the linear dependence of the contact area on the normal force [7, 8, 15]
A = kFn/(E*Vh) (4)
was confirmed and a new power law dependence of the contact stiffness k = dFn/dd was found:
/ \a
■ = Z
E*jA0 1 E*hjA0
(5)
In the present paper, all 1D results will be compared to the exact 3D datasets obtained from boundary element method.
5. Reduction method
An approach that is totally different from the classical numerical methods was first introduced by Popov and Psa-khie [1]. They showed that the initial 3D contact can be mapped to a 1-dimensional line, consisting of a set of independent springs. In other words, the original three-dimensional body can be formally substituted by an equivalent one-dimensional "rough line" with an elastic foundation having a normal stiffness per length equal to
k = E *. (6)
To be equivalent to the initial three-dimensional body, the one-dimensional spectral density C1D(q) must be defined according to:
CiD(q ) = nqC2D(| q |). (7)
This choice ensures that the one-dimensional rough line has the same rms roughness (as well as the same moments and moments of derivatives) of the surface profile as the original three-dimensional line. In the case of self-affine rough surfaces, the 3D and 1D bodies will have the same Hurst exponent. The validity has been proven analytically for singular contacts of any bodies of revolution [16]. The reduction method has been elaborated in many publications and will not be discussed in detail here. It shall be noted that Geike [17] also investigated the distribution and curvature of the surface's roughness peaks for very specific, narrow-banded surfaces and found good agreement. The reduction method has also been applied to rough surfaces
Fig. 2. Interpretation of an indentation of the rough line. The overall stiffness is the sum of springs in contact k = E*/n ^ L., and the normal force is the shaded area FN = E A^c (- d), where xc are the positions of the springs in contact and A is the grid spacing
[14], but up until recently [13], no 3D data for comparison was available.
The contact stiffness of a given configuration can be easily found in the 1D case. As all springs individually act in parallel, the overall stiffness is simply the total length of the springs in contact, multiplied with the specific stiffness.
The stiffness has a natural saturation value both in 1D and 3D. In the latter case, it is obtained by forcing the imaginary indenter so far into the opposite half-space that the surface roughness can be neglected and only the macroscopic shape dominates the stiffness. In the above-mentioned simulation, a square indenter shape was chosen, so that the stiffness approached the theoretical value of ks = 1.1419E*^JA0 for very high normal forces [15]. Interestingly, the saturated stiffness value can be found long before all grid points have come into contact. In fact, only a relatively small real contact area is needed.
In the 1D model, the saturated stiffness is only obtained when all springs have entered into contact.
6. Results
In the one-dimensional model, the stiffness as a function of the normal force has been obtained for 6 values of the fractal dimension Df = 2.0, 2.2, 2.4, 2.6, 2.8, 3.0 by averaging 500 surface realizations. Results are shown in Fig. 3. For low to medium normal forces, a power-law behavior can be observed in the logarithmic plot, with an exponent varying from 0.5 to 0.75. For very high forces, a saturation value is reached, when all springs have come into contact. These results agree very well with the results from the direct 3D simulation, see Fig. 4.
Furthermore, the 1D approach allows for the investigation of highly discretized surfaces. It, thereby, shows that the power-law can correctly describe the contact stiffness for very low forces as well. These low forces could not be represented in 3D due to numerical limitations. In both models, the power-law behavior was observed down to the lowest forces that were representable on the grid.
10°
10"
10"
10"
: 1D yy/// 3D plot limits-
10"
10"
10°
Fig. 3. Dimensionless contact stiffness k1D as a function of the applied dimensionless load F1D for 6 different values of the fractal dimension Df = 2.0, 2.2, 2.4, 2.6, 2.8, 3.0 following the indicated arrow. Results are obtained using the 1D reduction method. k1D =dF/dd(EL) FD = F (E *hL)-1
7. Analytical considerations
In both investigations, the power a in Eq. (5) varied between 0.5 and 0.8 for fractal dimensions of the surface varying from 2 to 3. These values were obtained by fitting the resulting dependencies with a logarithmic least-square approximation. In the 3D boundary-element calculations, these dependencies were found to be the average stiffness value from an ensemble of 60 surfaces, having the same fractal power spectrum, but different random phases. The dependency of the exponent a as a function of the fractal dimension Df of the surface was observed to be a ^ ^ 0.2567 Df. As the 1D reduction method allows for simulations to be done in almost no time, statistics were taken for 500 surfaces. Furthermore the increased grid resolutions allowed for a much greater range in the normal force applied. Therefore, the exponent a could be deduced with great accuracy from the 1D simulation data.
We now want to estimate the dependence of the exponent a in (5) on the fractal dimension analytically.
Fig. 4. Dimensionless contact stiffness k3D over the dimensionless load F3D taken from [13]. Results are obtained using the 3D boundary element method. Fractal dimensions are the same as in Fig. 3. k3D =3F/dd x x (1.1419E*y[A0)-1, F,D = FEhyflO)-1
in z:
Let us assume a square indenter with a fractal rough surface h and an edge length L pressed into the elastic half space with a given indentation depth d. The contact configuration will be the distribution of contact areas and areas of non-contact:
A : [- L 2, L 2f a {'contact','no contact'}. (8)
Now, consider the same indenter magnified by a factor C in x and y, while scaled by CH
h'(x, y) = CHh(x/ C, y/C), L = CL.
We will get the same scaled contact configuration if the indentation d is scaled in the same way as h:
d ' = CHd,
A' ( x, y ) = A( x/C, y/C ).
The contact stiffness is only a function of A(x, y) and it scales linearly with the magnification, resulting in
k = Ck. (11)
Considering k = dF/dd, we get
(9)
(10)
( C
CHC
(12)
xa
Ff = Cl+HF .
Entering (9), (11) and (12) in (5) leads to C =
C
underlining again the validity of this equation.
The fractal dimension of the indenter is directly related to the Hurst exponent H = 3 - Df, which graphically describes the nature of a fractal rough (or self-affine) surface. A magnified version of the surface cannot be distinguished from the original and a different magnification is allowed for the axis perpendicular to the surface (z) with respect to x and y. The transformed surface h = CHh(x/C, y/C) "looks" exactly like the original one, z = h(x, y) [7].
This means that if the original force F is applied to the scaled surface h and we assume that
A'(x, y):
[A(x, y), V(x, y) e [-L2, L/2]2, I'no contact' else,
(13)
then, of course, the stiffness will still be k. Entering this
i 1 \a
into (5) this time leads to
a = -
1
1 C
1
KCHCj
resulting in
H +1
(14)
Equation (13) is a very strict assumption, but it can be seen that in order to only obtain the same stiffness value, the new contact configuration does not have to be exactly the same but only "look" the same and its position within the larger [-Lr/2, Ly2]2 is not of importance. It is only vital that all contact patches appear in a region with spatial extension of approximately [0, L] . This assumption fails for very low values of H when contact spots are widely spread over the entire indenter.
As can be seen from Fig. 5, the analytical prediction
a =-=- matches the values obtained from si-
H +1 4 - Df
+ 3D data - - 1D data / / / / /
-— 0.2567 Df * V / / x +
- —- 1/(4 - Df)
i i i i I I
2.0 2.2 2.4 2.6 2.8 D,
Fig. 5. The exponent a from Eq. (5) plotted over the fractal dimension. We show the results for 3D and 1D calculations and two approximations, the latter one being deducted from analytical considerations. The 1D curve in this figure was obtained by evaluating 100 values of Df with 500 rough lines in each case
mulation extremely well for Df < 2.6. For higher fractal dimensions, a small discrepancy can be observed.
8. Conclusion
The 1D reduction method was used to study the contact stiffness between rough surfaces. In order to account for short and long range interactions, fractal surfaces were chosen and results were compared to 3D boundary-element calculations. The results from the 1D model clearly confirm the power law behavior and suggest its validity for very low forces, which could not be simulated in 3D. As the reduction method uses a set of independent springs, no numerical difficulties arise. In the initial 3D boundary-element simulations, highly optimized multilevel-methods were applied but calculation time was several weeks for a
grid resolution of 2049x2049. In contrast, the 1D model
20
allowed for highly discretized samples of 2 points to be processed within minutes. Furthermore, a simple explana-
tion for the power-law behavior is given, supporting not only the power-law itself but also predicting a value for the exponent.
References
1. Popov V.L., Psakhie S.G. Numerical simulation methods in tribology // Tribol. Int. - 2007. - V. 40. - P. 916-923.
2. Bowden F.P., Tabor D. The Friction and Lubrication of Solids. - Oxford:
Clarendon Press, 1986. - 424 p.
3. Archard J.F. Elastic deformation and the laws of friction // Proc. R. Soc. A. - 1957. - V. 243. - P. 190-205.
4. Persson B.N.J. Elastoplastic contact between randomly rough surfaces // Phys. Rev. Lett. - 2001. - V. 87. - No. 11. - P. 116101.
5. Greenwood J.A., Williamson J.B.P. Contact of nominally flat surfaces // Proc. R. Soc. A. - 1966. - V. 295. - P. 300-319.
6. Hyun S., Pei L., Molinari J.-F., Robbins M.O. Finite-element analysis of contact between elastic self-affine surfaces // Phys. Rev. E. - 2004. -V. 70. - P. 026117.
7. Persson B.N.J. Contact mechanics for randomly rough surfaces // Surf.
Sci. Rep. - 2006. - V. 61. - P. 201-227.
8. Hyun S., Robbins M.O. Elastic contact between rough surfaces: Effect of roughness at large and small wavelengths // Tribol. Int. - 2007. -V. 40. - P. 1413-1422.
9. Leidner M. Kontaktphysikalische Simulation von Schichtsystemen: PhD Thesis. - Technische Universität Darmstadt, 2009.
10. Venner C.H., Lubrecht A.A. Multilevel Methods in Lubrication. -Amsterdam: Elsevier, 2000. - 379 p.
11. Brandt A., Lubrecht A.A. Multilevel matrix multiplication and fast solution of integral equations // J. Comput. Phys. - 1990. - V. 2. -P. 348-370.
12. Polonsky I., Keer L. A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques // Wear. - 1999. - V. 231. - P. 206-219.
13. Pohrt R., Popov V.L. Normal contact stiffness of elastic solids with fractal rough surfaces // Phys. Rev. Lett. - 2012. - V. 108. - P. 104301.
14. Popov V.L., Filippov A.E. Force of friction between fractal rough surface and elastomer // Tech. Phys. Lett. - 2010. - V. 36. - No. 9. -P. 525-527.
15. Popov V.L.Contact Mechanics and Friction. Physical Principles and Applications. - Berlin: Springer-Verlag, 2010. - 362 p.
16. Heß M. Über die exakte Abbildung ausgewählter dreidimensionaler Kontakte auf Systeme mit niedrigerer räumlicher Dimension. -Göttingen: Cuvillier-Verlag, 2011. - 172 p.
17. Geike T. Theoretische Grundlagen eines schnellen Berechnungsverfahrens für den Kontakt rauer Oberflächen. - Berlin: TU Berlin, 2007.
Поступила в редакцию 10.06.2012 г.
Сведения об авторах
Pohrt Roman, Dipl.-Ing., Berlin University of Technology, Germany, [email protected]
Попов Валентин Леонидович, д.ф.-м.н., проф. Берлинского технического университета, [email protected]
