# Approximation of the Integrals of the Gaussian Distribution of Asperity Heights in the Greenwood-Tripp Contact Model of Two Rough Surfaces Revisited

Authors: Radosław Jedynak; Jacek Gilewicz
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Some examples of application of Padé approximant techniques to approximate the integrals in question and comparison with previous results, essentially the recent Panayi-Schock and Green-wood results, show the efficiency of this simple natural approximation.


---

## Body

## 1. Introduction

Modelling of the contact between moving rough surfaces allows a better understanding of friction and wear mechanisms, which can be used in engineering solutions. This issue has been examined using a number of approaches. The statistical type of a contact model is still the most popular model used in rough surfaces contact. This statement does not mean that it is the best solution of contact rough surfaces. Instead of using the complete roughness data, only probability density function is used. This function means the probability of the asperity with the height betweenh and h+dh. The first well-known statistical model was introduced by Greenwood and Williamson [1] (GW). They joined a statistical process with a classical Hertzian contact to deal with the rough surfaces contact. They adopted the following assumptions: (1) the asperity height distribution is Gaussian, (2) asperity contact is modelled by the Hertzian spherical contact theory, (3) the asperity tip radius is assumed constant, and (4) adhesion contact between asperities is ignored. A rough surface was described by three parameters: (1) standard deviation of asperity height distribution, (2) average asperity summit radius of curvature, and (3) areal asperity density. This model has been widely accepted and developed by numerous researchers. The main reason of its popularity is its simplicity and it’s predictions are in accordance with the carried out experiments. Some interesting review articles in the frame of dry friction are written by Bhushan [2, 3], Buczkowski and Kleiber [4], and Jedynak and Sulek [5].The model proposed by Greenwood and Tripp [6] (GT) extended the GW model to contact between two rough surfaces. They showed that the contact between two rough surfaces can be modelled by a contact between an equivalent single rough surface and a flat one. The equivalent rough surface is characterized by an asperity curvature and the peak-height distribution of the equivalent surface. They used the simple formula for a standard deviation of the statistical distribution as the square root of the sum of squares of the standard deviations of asperity height distributions on the two surfaces. The GT model gained large popularity in the field of elastohydrodynamic analysis, and all the following papers refer to [6]. The most frequently cited equations given by the GT model are for the following asperity contact area:
(1)Ae=π2(ξκσ)2AF2(λ)
and load carried by the following asperities:
(2)We=8215π(ξκσ)2σκE*AF5/2(λ),
where ξκσ—roughness parameter, A—nominal contact area,—Stribeck oil film parameter, first defined by Stribeck [3, 6] as λ=h/σ, E*—effective elastic modulus, andF2, F5/2(λ)—statistical functions introduced to match the assumed Gaussian distribution of asperities.A parameter,λ (ratio of film thickness to the composite roughness of the contiguous surfaces), is used to ascertain any boundary contributions, which occur because of asperity interactions. It is a very important parameter for EHD analysis because interruptions in a coherent film may occur when λ<3 [7].At this point, it is worth noting that GT equations are only valid for surfaces which have a Gaussian distribution of asperities. A lot of real engineering surfaces do not meet these requirements. For example, a typical analysis of the compression ring of a piston against the cylinder bore is not correct if we use a Gaussian distribution. Spencer et al. [8] showed that the bore usually has a liner insert, which is crosshatch honed, making a plateau geometry. This geometry provides a smooth surface to allow for hydrodynamic film buildup between the piston rings and the cylinder liner surface. On the basis of these facts we can definitely say that the mentioned surfaces do not have a Gaussian distribution of asperities ringcylinder.On the other hand, cam-tappet contact is a good example for using the mentioned distribution. Both surfaces are well known to follow Gaussian topography. This kind of conjunction is the most loaded contact in the valve train mechanism. Jedynak and Sulek [5] showed that some kinds of galvanic coatings have also a Gaussian distribution of asperities.Many surfaces do not follow at least a symmetric Gaussian-type distribution. Some researchers assume that the rubber-type material has a symmetric Gaussian-type distribution, but it is against the actually presented results, for example, Prokopovich et al. [9]. This paper clearly shows skew-symmetric Gaussian-type distribution for rubber seals. The mechanism of friction generation in the seal conjunction is complex, arising from adhesion of rubber in contact with the moving interface, viscous action of a thin film of fluid, and deformation of seal asperities. The authors developed a friction model based on the said mechanisms and then performed tests at both nano- and component level scales. Because the topography of the examined surfaces was nonsymmetric Gaussian, they used an alternative approach for adhesive friction, which is described in [3, 10]. They found that adhesion is a dominant mechanism of friction.The GT model is still being superseded by better models. Recently, Chong et al. [11] presented a friction model at fractal level, combining the effects of asperity adhesion, elastoplastic deformation, and a thin adsorbed boundary active molecular layer. This model forms the basis for the development of a generic multispecies lubricant-random roughness surface methodology. They stated that traditional boundary friction models do not apply to randomly distributed topology.Statistical functionsFn(h) which are used in extended GT and GW models have the following form:
(3)Fn(h)=12π∫h∞(s-h)ne-s2/2ds,n=12,1,32,2,52.
This function can only be solved using numerical methods and that is why we discuss the most efficient approximations formulas for F5/2(h) in the next chapter. They are used by many researchers in elastohydrodynamic analysis, In this case F5/2 is a function of λ. Then we propose a new formula which is based on Padé approximation method.It is important to emphasize that the approximations discussed next can be only applicable to surfaces which have a Gaussian distribution of asperities but not for all engineering surfaces which actually exist.


## 2. Previous Approximations

The GT model required the numerical integration of the Gaussian distribution. Greenwood and Tripp [6] provided tabulated results for this integral over its effective range: 9 points on the interval [0,4], 0≤Fn<1 for all n, the last digits given in [6] corresponded to 10-5 and from h=4 these values were as follows:
(4)F3/2(4)=F2(4)=F5/2(4)=0.00000,Fn(h>4)=0.
These values were recently used by a number of authors [7, 12, 13] to approximate the integrals Fn by simpler functions. Unfortunately, all these authors used old tabulated values, which were recently corrected [14] as shown in Table 1 for F5/2(h).Table 1
Computed values ofF5/2(h) used to build the approximation formulas.
h
Old values [6]
Improved values [14]
0
.616 64
.616 634
.5
.240 40
.240402
1
.080 56
.0805623
1.5
.022 86
.022 8631
2
.005 42
.00542371
2.5
.001 06
.001063 69
3
.000 17
.0001709
3.5
.000 02
.000 022 31
4
.000 00
.000 002 354
4.5
—
.000 000 2
5
—
.000 000 0
⋮
⋮
⋮
6.5
10
-
12The following models of approximation ofF5/2 based on the coefficient fitting to the tabulated values were proposed:(i)
Hu et al. power law approximation [13]:
(5)F5/2(C)(h)=α(4.0-h)βh<4.0,F5/2(h≥4)=0,(ii)
Teodorescu et al. 5-degree polynomial approximation [7]:
(6)F5/2(A)(h)=a0+a1h+⋯+a5h5h<4.0,F5/2(h≥4)=0,(iii)
Panayi and Schock exponential form [12]:
(7)F5/2(P)(h)=e(-.48347-1.65118h+.06043h2)/(1-.19256h+.01610h2)-.00000350796forh<4.0,F5/2(P)(h≥4)=0,(iv)
Greenwood quartic fit tolnFn [14]:
(8)F5/2(G)(h)=.616634×e.00121735h4-.02342695h3-.2696729h2-1.7433954h.The paper [14] contains the analysis by Greenwood of the paper [12], the corrected tabulated values of F5/2, and the model (7), which remains better than the Panayi-Schock model even after the correction of the coefficients in (6) by the new values of F5/2. Greenwood proposed also in [14] to consider the following modified function:
(9)Gn(h)=eh2/2Fn(h)
which is more adapted to the numerical treatment, because it is slower decreasing than Fn, as shown in Figure 1 in the case of n=5/2. Figure 1 shows both function F5/2(h) and G5/2(h) for h∈[0,4].Figure 1
Graphs ofF5/2(h) and G5/2(h) for h∈[0,4].


## 3. Padé Approximation (PA) ofF5/2 or of G5/2

The authors of the present paper are surprised that no Padé approximation was tried to approximateF5/2 or G5/2 or to approximate the integrated function. It is an automatic and low-cost method leading frequently to the surprisingly good results.Letf be an analytic function at N real different points having the following power expansions:
(10)∑k=0∞ck(xj)(x-xj)k,j=1,…,N.
In practical situations, we know only a few first coefficients of each expansion, and then we have to deal with the limited information characterized by the following truncated power series:
(11)∑k=0pj-1‍ck(xj)(x-xj)k+O((x-xj)jp),j=1,…,N.
The N-point Padé approximant (NPA) to f, if it exists, is a rational function Pm/Qn=[m/n] denoted, if it is needed, as follows:
(12)[mn]x1x2⋯xNp1p2⋯pN(x)=a0+a1x+⋯+amxm1+b1x+⋯+bnxn,m+n+1=p=p1+p2+⋯+pN
satisfying the following relations:
(13)f(x)-[mn](x)=O((x-xj)pj),j=1,2,…N.
Each pj represents the number of coefficients ck(xj) of expansion actually used for the computation of NPA. In the present problem, all pj=1, and the construction of NPA [m/n] is equivalent to the construction of the rational interpolation Pm/Qn. The classical Padé approximant (PA) is a one-point PA computed from the development of the function f at the origin x=0. It is defined by the linear system obtained by the definition (13) with right-hand side O(xm+n+1).At the first time, let us estimate crudely the functionF5/2 by an exponential function in the interval [0,4]:
(14)F5/2(h)=αe-βh;α=.616634;β=-14ln.000002354.616634.
Of course this simple interpolation over two nodes badly reproduces F5/2 on this interval; however, it shows that our Padé approximant must decrease on [0,4] (8) and must become rapidly zero after. That is, the degree of its denominator must be sensibly greater than those of numerator.To obtain the optimal results using the Padé approximation method, one must select the “best” PA[m/n]; it means the “best” degrees m and n. Many practical algorithms of choice of the best PA were proposed in [15]. For the present problem, we select two following algorithms: ρ-method and the method of valleys in c-table, both based on the analysis of the coefficients of the expansion of the considered function f(x)=∑‍cjxj. The ρ-method analyze the sequence {ρn=|cn/cn+1|}. The nonmonotonicity of the few first terms of {ρn}, if the presumed PA is such that m>n, indicates the difference m-n. The c-table is the table of determinants of matrices of linear systems defining the denominators of the classical (one-point) PA [m/n]:
(15)Cnm=|cmcm-1⋯cm-n+1cm+1cm⋯⋯⋮⋮⋮⋮cm+n-1⋯⋯cm|.
At the same time, the dominant error of {f-[m/n]} near the origin is proportional to Cnm. The minima of each antidiagonal in the table of |Cnm|, of the entries builded with the same information, indicate the paradiagonal with fixed m-n of the best PA.By the integration of parts of (3), one obtains the following recurrence formulae:
(16)Fn+1(h)=nFn-1(h)-hFn(h).
The following formulas for the derivatives are also useful to obtain the development in power series of Fn:
(17)Fn′(h)=-hFn(h)-Fn+1(h)=-nFn-1(h).
Now we can compute the expansion in the power series of F5/2 near the origin, where a≔F1/2(0)=.411089 and b≔F3/2(0)=.43002:
(18)F5/2(h)=(3a2)-(5b2)h+(15a8)h2-(5b8)h3+(5a64)h4+(b64)h5-(a256)h6-(5b5376)h7+⋯=.616634-1.07505h+.770793h2-.268762h3+.0321164h4+.00671906h5-.00160582h6-.000399944h7+.000100364h8+.0000249965h9+⋯,
and the power series expansion of the function G5/2:
(19)G5/2(h)=eh2/2F5/2(h)=.61663422-1.07505h+1.07911h2-.8062875h3+.49459203h4-.262043h5+.12364801h6-.0530326h7+.0209760h8-.0077339h9+⋯.
The ρ-method consists of the analysis of the regularity of the sequence {ρn=|cn/cn+1|}. The beginning of the sequence {ρn} corresponding to (19) is regular:
(20)ρ0,ρ1…ρ7=.57,.996,1.34,1.63,1.89,2.12,2.33,2.53.
This indicates that the rational approximation {m/n} corresponding to this series is such that m≤n. More significant result is observed studying the c-table corresponding to (19), where the notation .1-5 means 0.1×10-5 (see Table 2).Table 2
m
n
0
1
2
3
4
5
6
7
8
0
1
.61
.38
.23
.145
.89-1
.55-1
.34-1
.21-1
1
1
−1.08
.49
−.12
.15-1
−.47-3
−.72-4
.104-5
.98-6
2
1
1.08
.30
.29-1
.11-2
.14-4
.10-6
.2-8
3
1
−.81
.12
−.45-2
.53-4
−.19-6
.16-9
4
1
.49
.33-1
.48-3
.18-5
.2-8
5
1
−.26
.75-2
−.38-4
.42-7
6
1
.12
.14-2
.24-5
7
1
−.053
.22-3
8
1
.021The minima are located on the positions0/3, 1/4, 2/5 and 3/6. The best PA is then as follows:
(21)[36]G(h)=(h3.616634+.0587436h+.0186433h2-.0000525721h3)×(1.+1.83868h+1.48582h2+.680187h3+.187771h4+.0300186h5+.0022127h6)-1.
In Table 3 we compare the exact values of F5/2(h) with those obtained by multiplying [3/5]G and [3/6]G by e-(h2/2) and with the values calculated by the Greenwood approximation [14] and the Panayi approximation [12].Table 3
Approximation ofF5/2 via one point (h=0) Padé approximants of G5/2.
h
exactF5/2(h)
e
-
h
2
/
2
[
3
/
5
]
G
(
h
)
e
-
h
2
/
2
[
3
/
6
]
G
(
h
)
Greenwood
Panayi
0.
.616634
.616634
.616634
.616634
.616636
.5
.240402
.240402
.240402
.240401
.240395
1.
.0805623
.0805623
.0805623
.0805597
.0805637
1.5
.0228631
.0228627
.022863
.0228629
.0228449
2.
.00542371
.00542304
.00542366
.00542412
.00541762
2.5
.00106369
.00106315
.00106365
.00106379
.00106965
3.
.000170873
.000170616
.000170852
.000170864
.000175778
3.5
.0000223124
.0000222326
.0000223059
.000022307
.0000225912
4.
2.35338-6
2.33623-6
2.35197-6
2.35373-6
1.17725-9
The bold digits are exact.Clearly the[3/6] PA leads to the best values, but it is computed with 10 coefficients (information) on G5/2 at h=0. The results obtained with [3/5] PA (9 information) are better in the range 0÷2.5 than the Greenwood results (obtained also with 9 information). It is natural, because the PA calculated from the development of a function near the origin gives good results in the neighborhood of the origin. In the following, we will compare the N-point PA [3/5]h1h2⋯h911⋯1 and [3/6]h1h2⋯h921⋯1 with other results. The NPA [3/6], corresponding to our best choice, is easily calculated with 9 values of G5/2 and one additional coefficient of development (19):
(22)[35]G(h)=(h3.616634+.10073h+.00818188h2-.000161362h3)×(1.+1.90693h+1.58667h2+.740254h3+.19751h4+.0298511h5)-1.
Table 4 contains the values of the approximations also in the points located between the nodes, because the NPA are, of course, exact in the nodes, and we must analyze other points to compare our approximations with those of Greenwood and Panayi. The advantage to determine the values of F5/2(h) via e-h2/2[3/5]G(h) appears evident. In fact, the last values have in the intermediate points all digits exact and it is not the case of the values of the NPA [3/5]F(h) determined directly by the values given in [14] (Table 1):
(23)[35]F(h)=(h3.616634-.457287h+.1141502h2-.0095798h3)×(1.+1.3788066h-.997978h2+2.988216h3-1.651924h4+.5588187h5)-1.
Figures 2 and 3 illustrate this fact clearly.Table 4
9-point Padé approximation ofF5/2 obtained directly and also via approximation of G5/2. Comparison with Greenwood and Panayi results.
h
exactF5/2(h)
[
3
/
5
]
F
(
h
)
e
-
h
2
/
2
[
3
/
5
]
G
(
h
)
Greenwood
Panayi
0.
.616634
.616634
.616634
.616634
.616636
.25
.391978
.3849
.391977
.39198
.391941
.5
.240402
.240402
.240402
.240401
.240395
.75
.141962
.142463
.141962
.141959
.141971
1.
.0805623
.0805623
.0805623
.0805597
.0805637
1.25
.0438561
.0438004
.0438561
.0438548
.0438441
1.5
.0228631
.0228631
.0228631
.0228629
.0228449
1.75
.0113963
.0114051
.0113963
.0113967
.0113818
2.
.00542371
.00542371
.00542371
.00542412
.00541762
2.25
.00246124
.00245933
.00246124
.0024615
.00246296
2.5
.00106369
.00106369
.00106369
.00106379
.00106965
2.75
.00043732
.000437939
.00043732
.000437338
.000443886
3.
.000170873
.000170873
.000170873
.000170864
.000175778
3.25
.0000633931
.0000630551
.0000633931
.0000633825
.0000658778
3.5
.0000223124
.0000223124
.0000223124
.000022307
.0000225912
3.75
7.44495-6
.7.44495-6
7.44495-6
7.44359-6
6.11165-6
4.
2.35338-6
2.35338-6
2.35338-6
2.35373-6
1.17725-9
standard dev.:
stand. dev.:
stand. dev.:
9.6-18
9.3-7
7.3-6
The bold digits are exact.Figure 2
Graph of the error of Greenwood (7) approximation: F5/2(h)-F5/2(G)(h).Figure 3
Graph of the error of our approximation:F5/2(h)-e-(h2/2)[3/5]G(h).Here, we see that the simple NPA[3/5], that is, an NPA builded with the same information that Greenwood or Panayi approximations, is the best. The comparison of the standard deviations is significant. Other advantage consists in its simple calculation. The following figures representing the error of our approximations illustrate the sharpness of our results.It is possible to try other approximations using the Padé techniques. One can use other interpolation nodes forh>4 with values F5/2(h)=0 and in particular the point at infinity F5/2(∞)=0. For instance, if we compute 11 points Padé approximant [m/n](m+n+1=11), we can choose [0/10],[1/9],[2/8],[3/7],[4/6],… where the behavior at infinity is, respectively:
(24)α1h10,α2h8,α3h6,α4h4,α5h2,…,
where αi can be crudely estimated at h=4 as follows:
(25)α1=.000002354×410=2.36;…;α5=.000002354×42.
The numerical tests show that simple use of these additional nodes does not improve the approximation. The effect of appearing oscillations is negative.The Padé approximant method allows, if it is needed, to obtain the two-sided estimates of the approximated function. This consists of construction of one NPA[m/n] from K=m+n+1 given coefficients and of the second one, [m-1/n] or [m/n-1], with K-1 coefficients. In this case, if the first NPA is greater than approximated function between two nodes, the second is less than this function. Both NPA bound the approximated function from top and from below [16].


## 4. Conclusion

The efficiency of one-point as well asN-point Padé approximation method applied to the special tribology problem is presented and compared with other methods. This method starts from the choice of the best Padé approximant and is easily feasible. Our results confirm an ingenious idea of Greenwood proposing to approximate a modified function Gn, more smooth that the function Fn (see Figure 1), and then return to the calculus of Fn, which allows to reduce considerably the numerical errors.



---

*Source: 459280-2013-05-23.xml*