A numerical procedure for random generation of rough surfaces

Wear, 47 (1978) 263 - 277 0 Elsevier Sequoia S.A., Lausanne

263 - Printed

A NUMERICAL PROCEDURE ROUGH SURFACES

NADIR

in the Netherlands

FOR R~DOM

GENERATION

OF

PATIR

~eparf~e~t

of Me~~a~~cal Engineering, Northwestern

il~i~ersit~, Eua~ston, ~~lino~

(U.S.A.) (Received

May 19,1977)

Summary A numerical procedure for randomly generating any general tbreedimensional surface roughness with prescribed statistical properties is presented. Through the use of linear transformations on random matrices, this procedure is capable of generating Gaussian or non-Gaussian rough surfaces with any given surface autoco~elation function. As an example of the numerical procedure, a Gaussian surface having a predetermined exponential autocorrelation function is generated.

1. Introduction Although the effects of surface roughness on lubrication have been investigated during the past decade, numerical simulation has only recently been used. Archard et al. [l] used the digital form of profilometer readings of a profile for simulation of the contact of rough surfaces. Patir and Cheng [Z] derived an average Reynolds equation for rough surfaces based on a numerical solution of the Reynolds equation for a model bearing with a randomly generated surface roughness. It is considered that numerical simulation will be an important tool in the near future for analyzing the effects of surface roughness on tribological systems. Such a study would require the roughness heights of a surface, which can either be measured from a rough surface or generated numerically, to be in digital form. Although accurate measurement of a profile is relatively simple using a stylus profilometer, the measurement of a surface is much more complex. Surface me~uremen~ are done by taking a number of parallel profile measurements which requires an accurate relocation technique and an additional software requirement to align the profiles numerically. Randomly generating surface roughness by numerical means, however, is much simpler and offers certain advantages. In addition to eliminating all

264

the hardware requirements for surface measurement, it also eliminates the need to filter out the unwanted wavelengths from a measured surface. Furthermore any parametric study involving roughness requires surfaces with known statistical properties, and it is much more convenient to generate them numerically rather than to measure mechanically produced rough surfaces. This paper describes methods of generating random surface roughness with prescribed statistical properties as an aid to future research in roughness-related phenomena through simulation.

2. Autocorrelation

function

Most of the statistical parameters of a rough surface can be derived from two statistical functions: the frequency density function and the autocorrelation function (ACF). Therefore a convenient way to generate surfaces with known statistical parameters would be to generate surfaces having predetermined ACFs and frequency density functions. Since most engineering surfaces are approximately Gaussian, the frequency density of roughness heights can be assumed to be Gaussian. However, the algorithm will not put any restrictions on the ACF. To aid in the choice of a meaningful ACF, some of its properties are reviewed. The height z&y) of a rough surface may be considered as a twodimensional random variable. For convenience z is measured from the mean plane of the surface. Assuming the surface to be homogeneous (statistical properties are invariant with respect to a translation along the surface), the ACF of a surface is defined as

where E is the expectancy (averaging) operator and h, , h, are the delay lengths. R(O,O) is equal to u2 where u is the standard deviation of the roughness heights (r.m.s. roughness). The power spectral density (PSD) of a surface is defined as the Fourier transform of the ACF. Most of the vital statistics of a rough surface such as density and distribution of summit heights, mean summit curvature, mean surface gradient, r.m.s. surface gradient etc. can be related to the moments mpg of the PSD, which can be obtained through the formula

ap+q R&,X,) ahpx;

= iP+4m

P9

(2)

A,=0 hy=O

where i = m These relations are given by Nayak [3] for isotropic surfaces. For nonisotropic surfaces an extensive theory is not available at present, but such surfaces can be characterized by considering the properties of surface profiles in different directions.

265 The ACF of a profile

along a direction

0 is related to the surface ACF

where A, =

x* cos e

(4)

h, = he sin 0 Hence the ACFs of the x and y profiles are &(A,)

= R&,0)

&J(k)

= NO,&)

Since ient to use x = i Ax, y directions.

(5)

a digital form of surface roughness is sought, it is more convenan index notation. Let Zij denote the roughness anplitude at = j Ay, where Ax and Ay are the sampling intervals in the x and y Similarly, R,, is defined as

R pq = R@Ax,

4AY) = E(zijzi+p, j+q)

(6)

As the arguments of the ACF increase, R decreases to a small value, which is assumed here to be zero, so that a finite order autocorrelation matrix is obtained. Let n and m be two integers such that R,, is zero if p > n or q 2 m. This yields an n X m autocorrelation matrix. One can define the correlation length of a profile as the length at which the ACF becomes zero. Hence the correlation lengths of the x and y profiles are X: =nAx (7)

h;“,=mAy

It has been shown by Whitehouse and Archard [4] that many statistical parameters of a measured profile depend on the sampling interval. Choosing the sampling interval to be equal to the correlation length reveals the main structure of the roughness, while choosing it to be-equal to one-tenth of the correlation length reveals the short wavelength fine structure of the roughness. Therefore, n, m, Ax and Ay should be chosen carefully so as to generate the desired roughness structure; the number of calculations required for roughness generation increases significantly with increasing n and m. The dependence of the statistical properties on the sampling interval could also be seen from the following example. The square of the r.m.s. x slope ok, which is equal to mm, can be calculated from eqn. (2) by numerical differentiation (noting the symmetry of R about the origin):

(uy~m20=_;~ E

=_Rlo-~yR-lo A,=0 A,

=

,’

2

(R oo-RlO)

=o

(8)

266

3. Generating

Gaussian surfaces

It is possible to generate an N X M matrix of roughness amplitudes [Zij], having a Gaussian distribution of heights, and a given IZX m autocorrelation matrix [R,,] using linear transformations on random matrices. Using random number generators, one first generates an (N + n) X (M + m) matrix [Qij] whose components are independent identically distributed Gaussian random numbers with zero mean and unit standard deviation. The roughness heights are then obtained through the linear transformation

zij

=

i

2

k=l

=k17)i+k,

j+l

I=1

i = 1,2,...,N j=1,2 ,***, M

(9)

where akz are the coefficients to be determined so as to give the desired autocorrelation matrix. Since vii are independent and have unit variance, 1 if i=k,j=l

E(vijvkl)

=

I

Using this relation

R pq

=

along with the definition

n-P

m-q

x

c I=1

k=l

(10)

0 if otherwise of R,,

(eqn. (6)), we obtain

p =

0,l ,...,Iz - 1

q =

O,l,...,m

aklak+p,l+q

- 1

(11)

Equations (11) represent nm simultaneous non-linear equations for the determination of the coefficients akz. They can be solved by an iterative technique such as the Newton method. The details of the iterative scheme are summarized in the Appendix.

4. A numerical

example

To illustrate the numerical procedure, a rough surface having a Gaussian height distribution and an exponential ACF is generated. The exponential function has been used by many authors such as Whitehouse and Archard [4] and Peklenik [ 51 in modeling the ACF of a profile. This model can be extended to the ACF of non-isotropic surfaces by assuming that all the profiles on the surface have an exponential ACF with a decay constant depending on the orientation of the profile. A possible ACF for such a surface is

R(P,,P,) = 2 exp

]--2.3 i(2)’+($,‘i”‘l

where u is the r.m.s. roughness,

(12)

and p:, 0; are the lengths at which the ACFs

267

of the x and y profiles reduce to lengths have been defined as the by Peklenik [5]. The definition arbitrary; therefore this definition length. The ACF of a profile along

Re (PO) = u2exp

the 8 direction

is then

(-2.3 /$1)

where pi is the 0.10 correlation

6=

10% of their values at the origin. These correlation lengths of the x and y profiles of the correlation length is somewhat will be referred to as the 0.10 correlation

(13)

length of the 6 profile given by

P: 0: ((@Z

(14)

sin 0)2 + (0: cos 0)2}1’2

The non-isotropic properties of the surface can be seen from the polar coordinate representation of pi as suggested by Peklenik [ 51. For the assumed ACF this representation forms an ellipse with an ellipticity ratio n* (15) This ellipse can also be considered as the locus of all points whose heights have the same correlation with the height of a point located at the center of the ellipse. Therefore asperities should have roughly elliptical shapes with the given ellipticity ratio. The parameter y, which is the ratio of the x and y correlation lengths, shows the degree of non-isotropy of a rough surface. A value of y = 1 corresponds to an isotropic surface, while the limiting cases y = 0 and y = 00 correspond to one-dimensional transverse or longitudinal ridges. Having constructed an appropriate ACF model, we now consider the generation of a non-isotropic Gaussian surface having an ACF of the form of eqn. (12) with y = 2. The relevant parameters are chosen to be n=7

m= 4

Ax=Ay

0: = 6Ax

/3; = 3Ay

u= 1

(16)

Choosing u = 1 will produce normalized roughness amplitudes. To obtain roughness with a given r.m.s. value, it suffices to multiply the roughness amplitudes by the desired r.m.s. value. The roughness amplitude generated will be comparable with that of a real surface measured at sampling intervals Ax and Ay equal to l/6 and l/3 of the respective 0.10 correlation lengths. According to the convention of Whitehouse and Archard [ 41 this should reveal the intermediate structure of the roughness. In order to obtain a finite order autocorrelation matrix the exponential ACF is assumed to drop to zero after the 0.10 correlation lengths. Hence the digital form of the ACF is R,,=exp

[-2.3\(f)2+(i)2\1’2]

R,,=O

if pa7

or q>F17)

268

The coefficient matrix required to generate the desired roughness structure is obtained by solving the non-linear eqns. (11) using the Newton method outlined in the Appendix. The solution is given in Table 1. TABLE 1 Theoretical ACF and the coefficient

matrix obtained from a solution of eqns. (11)

Normalized ACF

0

1 2 3 4 5 6

0

1

2

3

1.000000 0.681292 0.464159 0.316228 0.215443 0.146780 0.100000

0.464159 0.423957 0.337750 0.250654 0.179740 0.126611 0.088289

0.215443 0.205502 0.179740 0.146780 0.114075 0.085667 0.062827

0.100000 0.096874 0.088289 0.076202 0.062827 0.049922 0.038529

3

4

0.191422 0.281895 0.309520 0.415534 0.172936 0.067346 0.162762

0.116858 0.258345 0.414711 0.134551 0.157780 0.128041 0.208701

Coefficients of transformation matrix 2

1 0.184575 0.126058 0.084033 0.099300 -0.149043 -0.053339 0.135184

0.156891 0.156822 0.071952 -0.136758 -0.031144 -0.043678 0.122985

The desired roughness amplitudes are then obtained using the algorithm described in Section 3. Four different surfaces, each with 1600 points (40 X 40 matrix), were generated. To compare the ACF of the generated surfaces with the expected ACF, the x and y profile ACFs of the generated surfaces are calculated by R po =

(18)

Ef~i+~,j) =

R04 = E(zijZi,j+s) =

1

N(M--q)

N M-Q x

izl

c j=l

z..z.

. ’

(19)

‘I’+’

where N = M = 40. Figure 1 shows the x and y profile ACFs of the generated surfaces which agree well with the desired ACFs. Since the ACF is fundamentally a random function, one cannot expect the ACF of every generated surface to be identical. Good agreement between the expected and observed ACFs, however, depends critically on how well the random number generator

0.0 -r 0

(4

1.0

a R

oq

4

(b) Fig. 1. x and y profile ACFs of the generated surfaces: ed surfaces; @ expected ACF.

x profile ACF of the generat-

generates a set of mutually independent identically distributed Gaussian random numbers (~8). Figure 2 shows a contour plot of a generated surface, with contour levels 0, +a, k 20,~3a, The non-isotropy of the roughness can be seen from longer dimensions of the asperities in the x direction.

Fig. 2. Contour plot of a generated surface. Contour levels: 0, c U, + 2~, c 3~.

Figure 3 shows typical x and y profiles taken from a generated surface. The x profiles are smoother and have longer asperity dimensions than the y profiles, owing to the higher correlation length of the x profiles.

Fig. 3. Typical x and y profiles of a generated surface: (a) two consecutive x profiles; (b) two consecutive y profiles.

271

Comparison of consecutive profiles also reveals the difference in the correlation functions of the x and y profiles. The two consecutive y profiles, separated by a distance Ax, are highly correlated as seen by the matching of the valleys and asperities at corresponding points of the two profiles. This is due to the fact that the correlation between corresponding points of the two profiles is Rio = 0.68. However, the two consecutive x profiles are not as highly correlated since Rol = 0.46 is smaller.

5. Generating

non-Gaussian

surfaces

While linear transformations on independent Gaussian random variables result in Gaussian variables, most non-Gaussian random variables do not have this property. The frequency density function of the roughness heights, obtained through a linear transformation on [TQ] , will be a function of the frequency density of the [nij] and the coefficients of the transformation matrix. Therefore one should choose the frequency density of the random numbers [gij ] and the transformation matrix [akl] such that the roughness heights [ZQ] have the desired frequency density and the desired autocorrelation matrix. This involves an additional set of non-linear equations relating the frequency densities of [no] and [zij] which will be coupled to eqns. (11). Since the frequency density of a random variable is uniquely determined by its moments, these additional equations can relate the moments of the two frequency densities. The effect of linear transformations on the frequency density function of a random variable can be conveniently analyzed by utilizing characteristic functions. The Fourier transform Lz of the frequency density function fz of a random variable 2 is called the characteristic function of 2 [6] : +L,(u)

=

I -cc

fz(z)eiU” dz

(20)

where u is a real parameter and i = dq Recalling the definition of expectation, the characteristic function can also be defined as the expected value of eiUZ : L,(u)

= E(eiUZ)

(21)

If the characteristic function of a random variable is known, the frequency density function can be obtained through the inverse relation fz(z) =

&

jw

LZ(u)edur

du

-00

An important property of the characteristic function is that the moments M,(k) of the frequency density of a random variable 2 can be obtained by

272

M,(k)

= E[Zk]

= (-i)kLhkk)(0)

(23)

whenever the hth derivative Likk’ exists at the origin in u [6] . Let qij be mutually independent identically distributed random numbers having a frequency density fi, and moments M,(h). If the characteristic function of qij is L,(U), it can be shown [6] that the characteristic function of the Zij obtained through the linear transformation (9) is

&Zt”) = ii k=l

ii

'&(a,@)

(24)

I=1

The moments M,(h) of the frequency density of the roughness heights Zij C~II then be related to the moments M,(k) of the nij using eqn. (23). These relations, along with eqns. (ll), should then be solved for the coefficients ckl of the transformation matrix and the moments M,(k) of the frequency density nij. Equations (ll), however, implicitly assume that the first two moments (mean and variance) of 17~ are 0 and 1 respectively. Equations (11) also impose the necessary conditions on ckl that the roughness heights will have zero mean and the required variance (r.m.s. roughness). Therefore the equations relating the moments of Zij and nij should only be included for the third and higher moments. After these non-linear equations are solved for the coefficients of the transformation matrix and the moments of the frequency density of vii, the roughness heights Zij with the desired frequency and correlation properties can be obtained through the use of eqns. (9). The random number generator for the nij should generate independent random numbers having a frequency density function defined by the moments obtained through the solution. As can be seen from the complexity of the formulation the generation of non-Gaussian surfaces requires tedious effort and is impractical. However, if one sacrifices the freedom of choosing any type of autocorrelation function, the generation of Gaussian or non-Gaussian surfaces becomes very simple, as will be illustrated in the next section.

6. A simple method In the preceding sections a numerical procedure capable of generating Gaussian or non-Gaussian rough surfaces having a given ACF was described. Since no restrictions were placed on the shape of the ACF, this procedure required the solution of a system of non-linear equations. In this section a simple ACF will be chosen which will lead to a constant coefficient matrix. Hence generation of surfaces having this ACF will be much simpler since it will not require the solution of non-linear equations.

273

Consider

a family of ACFs of the form

NL, hy) =

lh,l G

x:

IX,l G

h:

(25)

otherwise

0

where X: and hz will be defined as the correlation lengths of the x and y profiles (i.e. the lengths at which the profile correlation functions become zero). This surface ACF results in linear profile autocorrelation functions in the x and y directions: R,(X,)

=

(

u2

1 - y

1

(26)

R&p..(l-$)

(27)

The profiles in the other directions have a parabolic ACF which resembles a decaying exponential function. This ACF, which seems rather simplistic, is still a reasonable approximation for modeling many engineering surfaces. Inspection of profile ACFs obtained from measured profiles by Peklenik [ 5 ] shows that most of them can be approximated by a linear ACF, at least for delay lengths smaller than the correlation length. Furthermore most engineering surfaces have strong directional properties in the x and y directions owing to either machining processes or wear. Therefore the shapes of the profile ACFs in the other directions are less critical than those of the x and y profiles. Apart from being easy to generate, another advantage of this type of ACF is that the surface is characterized by three parameters u, h: and X;, and most surface statistics can be obtained in terms of these parameters by simple functions. This makes the use of such surfaces suitable for parametric simulation studies. The discrete form of the ACF (eqn. (25)) is R,,

=02(1-;) (1-i)

p
(28)

q
where

R pp=O cl

=nAx

G =mAy

ifpan

or q>rn

(29)

Choosing large n and m reveals the fine structure of the roughness since this means that the sampling interval is chosen as a small fraction of the correlation length. The ratio 7 of the x and y correlation lengths indicates the degree of non-isotropy of the surface:

(30)

The generation of Gaussian or non-Gaussian surfaces having an ACF of the form (28) is accomplished by the linear transformation i = 1,2,...N j = 1,2,...M

(31)

where qij are mutually independent identically distributed random numbers with zero mean and unit variance. The frequency density of the vii should be chosen so that the roughness heights Zij have the desired frequency density. If a Gaussian distribution of heights is desired, it suffices to generate the nij with a Gaussian distribution with zero mean and unit variance. For non-Gaussian distributions eqn. (24) reduces to a simple relationship between the characteristic functions of the Zij and nij : Lz(l.4) = L”,“(au)

(321

where (33) Hence the characteristic L,(u) = ,ynm

function

of the Qij is

0 ;

The frequency density of the vii required to generate roughness heights Zij with a frequency density fz and a characteristic function Lz can then be obtained through the relation

As an example, consider generation of roughness heights having a gamma distribution. The gamma distribution is asymmetric with a positive skewness, and therefore may be convenient in modeling roughness having skewed dis~ibutions. The frequency density and characteristic function of a gamma variable with zero mean, standard deviation o and skewness Sk are t71 (2 - c-Q~-~ exp {-(.a -- d)/b} z>d f&) = (36) bT(c)

L,(U)

= eiud(l

- iub)-”

(37)

where 4 c=-.._ 8,”

d=-bc

(331

275

Using eqn. (34) the characteristic = exp jiu(&)/

L,(u)

function

of the vii is obtained

as

(39)

(1 - i,f)E’“m

Comparison of this equation gamma variable with parameters

with eqn. (37) reveals that 7) is also a

b, =A _ SJz(nmY 2 a C

c’ = --

nm d’ =

= -

4 (40)

Sjfnm

_d =-_b’c’ anm

Hence the nij will have a gamma distribution variance and a skewness S; given by Sk = Sk(nm)l12

with zero mean, unit

(41)

7. Conclusion With complex problems involving random inputs, numerical simulation is a promising tool for analyzing roughness-related phenomena in tribology. An essential requirement for such simulation studies is the ability numerically to generate rough surfaces which have statistical properties similar to real surfaces. Most of the statistical properties of a rough surface can be derived from a knowledge of two statistical functions, the frequency density and the ACF. Hence a good algorithm should be able to generate surfaces having prescribed frequency density functions and ACFs. To generate surfaces with a given ACF, linear transformations on random matrices are utilized. First, an input matrix, whose components are independent identically distributed random numbers, is generated using random number generators. The roughness heights are then obtained through a linear transformation on the input matrix. The ACF of the generated surface depends on the coefficients of the transformation matrix. To obtain the desired ACF these coefficients have to be determined by solving a system of non-linear equations. The desired frequency density of the roughness heights is obtained by controlling the frequency density of the random numbers in the input matrix. For a Gaussian distribution of roughness heights it suffices to generate the input matrix with a Gaussian density function. The generation of roughness heights with non-Gaussian density functions is more complicated because the linear transformation required for the desired ACF also changes the frequency density of the input matrix. Hence, additional

276

relations are required to obtain the frequency density of the input matrix in terms of the desired frequency density of the roughness heights. Generation of Gaussian or non-Gaussian surfaces is greatly simplified if the ACF is assumed to be of a specific form which yields a constant coefficient matrix. The resulting surface has linear profile ACFs in the x and y directions, and is a reasonable approximation to real surfaces.

Acknowledgments The author thanks Professor Herbert S. Cheng for his assistance during the course of this work. The financial support for this work was provided by a research grant from the National Science Foundation under NSF Grant ENG74-21330.

Nomenclature “kl

b, c,d E

fz L.Z Mz(k)

mm

R

sk

k,P,: LAY

coefficients of the transformation matrix parameters of the gamma distribution expectancy operator frequency density function of variable 2 characteristic function of variable 2 kth moment of 2 moment of the power spectral density autocorrelation function skewness coefficient (normalized third central moment) roughness height 0.10 correlation lengths of x and y profiles ratio of x and y correlation lengths sampling intervals (distance between successive x and y points on the generated surface) a matrix of independent identically distributed random numbers correlation lengths of x and y profiles r.m.s. roughness (standard deviation of roughness heights) r.m.s. x slope (standard deviation of x slopes of roughness)

References 1 J. F. Archard, R. T. Hunt and R. A. Onions, Stylus profilometry and the analysis of the contact of rough surfaces, Proc. UITAM Symp. on the Mechanics of the Contact of Deformable Bodies, Delft Univ. Press, Delft, 1975. 2 N. Patir and H. S. Cheng, An average flow model for determining effects of three dimensional roughness on partial hydrodynamic lubrication, J. Lubr. Technol., submitted for publication, paper no. 77 - Lub-17 (1977). 3 P. R. Nayak, Random process model of rough surfaces, J. Lubr. Technol., F 93 (1971) 398. 4 D. J. Whitehouse and J. F. Archard, The properties of random surfaces of significance in their contact, Proc. R. Sot. London, Ser. A, 316 (1970) 97.

211 J. Peklenik, New developments in surface characterization and measurements by means of random process analysis, Proc. Inst. Mech. Eng., London, 182 (3K) (1967 - 68) 108. W. B. Davenport, Jr., Probability and Random Processes, McGraw-Hill, New York, 1970. N. L. Johnson and S. Kotz, Continuous Univariate Distributions-l, Houghton Mifflin, Boston, 1970.

Appendix Iterative solution by Newton’s method The non-linear system of eqns. (11) is solved iteratively relation a V+l = a” - [J”]-lf(a”)

using the

V = O,l,...

where a=

{a11a12.. . alma21.. . and

f=

~foofo~...fo.m-1flo...fn-1.m-l~ m--Q

n-P

f P9

=&

= I=1

and J is the Jacobian dfPll

_

akiak+p,

2+q

-

RP,

matrix having the components

”

J:: = - ai+p, j+q + aLp, j-q aayj where r=pm+q+l s=(i-l)m+j An initial approximation the formulae

to the coefficient

0 aij = SCij

where Cij

=

Ri-l.j-1 (n -i

f: i=l

+ l)(m -i

5 c; j=l

+ 1)

vector may be obtained

by