2003_Simulation of non-Gaussian surfaces with FFT

ÃTel.:+886-2-28239548.E-mail addresses:jjw5277@;jjwu@faculty.pccu.e-du.tw(J.J Wu).0301-679X/$-see front matter#2003Published by Elsevier Ltd. doi:10.1016/j.triboint.2003.11.005Therefore,statistical parameter,FFT and Johnson translator system will be discussed in this section.2.1.Statistical parametersSome statistical parameters are used to characterize a surface.For a surface,thefirst moment is the mean height,z which is generally removed before data pro-cessing and is therefore zero.z¼ð1À1zpðzÞd zwhere p(z)is the probability density function of surface height z.The second moment is the variance,r2which is the square of the RMS value or standard deviation.r2¼ð1À1z2pðzÞd zThe third moment is the skewness Sk,which represents asymmetric spread of the height distribution.Sk¼rÀ3ð1À1z3pðzÞd zThe fourth moment,K,represents the peakedness of the distribution.K¼rÀ4ð1À1z4pðzÞd zA Gaussian surface has zero skewness and a kurtosis of3,but a non-Gaussian surface has various Sk and K values.If a surface is described by a set of data f z r j r¼0;1;2;ÁÁÁ;NÀ1g,these four moments can be obtained byz¼1X NÀ1r¼0z rr2¼1NX NÀ1r¼0z2rSk¼1NX NÀ1r¼0z3rr¼1NX NÀ1r¼0z3r1NX NÀ1r¼0z2r!3=2K¼1NX NÀ1r¼0z4ir4¼1NX NÀ1r¼0z4r1NX NÀ1r¼0z2r!22.2.Fast Fourier TransformFast Fourier Transform(FFT)is a useful tool in generating a profile or a surface.Suppose a profile is described by a set of points{z r r¼0;1;2;...;NÀ1j}. Its FFT is defined asZ k¼1NX NÀ1r¼0z r expÀi2p krN340J.J Wu/Tribology International37(2004)339–346fork ¼0;1;2;...;N À1The IFFT of Z k is z r ¼XN À1k ¼0Z k expi 2p krNThe auto-correlation function of {z r }isR r ¼1X NÀ1s ¼0z s z s þrr ¼0;1;2;...;N À1where z r þs ¼z r þs ÀN for r þs >N .Its spectral densityisS k ¼1N X N À1r ¼0R r exp Ài 2p krN 2.3.Johnson’s system of translationNon-Gaussian random series is necessary for gener-ating a non-Gaussian surface.The non-Gaussian ran-dom series with different skewness and kurtosis can be generated by Johnson’s translator system [11].The Johnson system of frequency curves based on the method of moments provides a family of curves that can be used to generate an equation for the distri-bution for which the first four moments are known.In this system,there are three main types of frequency curves:1.The unbound system ðS U Þ:g 0¼c þd sinh À1g 00Àn k 2.The logarithm normal system ðS L Þ:g 0¼c þd log ðg 00ÀnkÞwhere g 00>n 3.The bounded system ðS B Þ:g 0¼c þd log ðg 00Ànn þk ÀgÞwhere n <g 00<n þk where g 0is Gaussian random series,g 00is the non-Gaussian random series with given skewness and kur-tosis,and c ,d ,f and k are constants to be determinedfor the first four given moments by using the method of moments.The schemes for generating such frequency curves was suggest by Hill et al.[11].3.Current methodsCurrently there are two methods using FFT to gen-erate non-Gaussian surfaces.They will be discussed in this section.But for convenience,we use one-dimen-sional profiles as examples.3.1.Hu and Tonder’s methodA digital filter technique was employed by Hu andTonder [8]to generate an output sequence of known auto-correlation function by a linear transformation system.Supposed f g r g is a series of random numbers,and the simulated rough profiles heights are taken as z r ¼X N À1s ¼0h s g s þrr ;s ¼0;1;2;...;ðN À1Þð1ÞIf f g r g is Gaussian random series,z r is a Gaussian pro-file.The Fourier transform of this equation is given by Z k ¼H k A kwhere A k is the FFT of g and H k is the transfer func-tion of the system given by H k ¼XN À1r ¼0h r exp Ài 2p krNk ;r ¼0;1;2;...;ðN À1Þð2ÞThe inverse Fourier transform of H k is the filter func-tion hh r ¼1N X N À1k ¼0H k exp Ài 2p krNð3ÞTheir spectral density has the following relation.S z ;k ¼j H k j 2S g ;kð4ÞSince S g ;k is the spectral density of random series,it is constant.The spectral density S k can be obtained from its auto-correlation function (ACF)R r .S k ¼1X N À1r ¼0R r exp Ài 2p kr ð5ÞFor the generation of non-Gaussian random sur-faces,the Gaussian input sequence is transformed to an input sequence with appropriate skewness Sk g and kur-tosis K g by using Johnson translator system of distri-bution.The skewness Sk z and kurtosis K z of the output sequence are related to those of the input sequence.Sk z ¼Pm À1r ¼0h 3rP m À1r ¼0h 2r3=2ðSk g Þð6ÞK z ¼Pm À1r ¼0h 4r K gþ6P m À2r ¼0P m À1p ¼r þ1h 2r h 2pP m À1r ¼0h 2r2ð7ÞJ.J Wu /Tribology International 37(2004)339–346341Given skewness Sk z,kurtosis K z and an auto-corre-lation function R k,the non-Gaussian profile can be obtained by the following procedure.1.Obtain the spectral density by Eq.(5).2.Obtain H k by Eq.(4).3.Obtain h r by Eq.(3).4.Find the skewness Sk g and kurtosis K g of a randomseries by Eqs.(6)and(7).5.Generate a non-Gaussian random series g by Johnsontranslator system.ing Eq.(1)to generate a non-Gaussian profile. However,there are some problems about this method. Theoretically,the spectral density and ACF of a random series areSðxÞ¼1ð8ÞRðsÞ¼1s¼00s¼0ð9ÞPractically,they are not so perfect.Fig.1shows the spectral density of a random series,which is generated by the software‘‘matlab’’.Fig.2shows the ACF.From Eq.(4),it is clear that the imperfection of random series will cause the spectral density and ACF of generated surfaces distorted.Wu[7]also found that there are mathematical mistakes in this method.Although Chilamankuri and Bhushan[12]found that the average ACF of profiles for generated surface is very close to the desired one.Wu[7]found that,in some cases,the ACF of generated surfaces are not so good as expected.Therefore,this method has its limita-tions.3.2.Non-Gaussian wind pressure time seriesNon-Gaussian wind pressure time series was investi-gated by Seong and Peterka[9],and also by Kumar and Stathopoulos[10].Although they investigated wind pressure time series,their methods can be used for generating non-Gaussian surfaces or profiles,too. Since the phase part does not affect the spectral ampli-tude characteristics,they found that IFFT is a good tool to generate the wind pressure time series.Using the same method,a profile can also be gener-ated by IFFT.z r¼X NÀ1k¼0ffiffiffiffiffiS kpexp i2p/kþkrMr¼0;1;2;...;ðNÀ1Þð10ÞThe phases are important for generating a non-Gaus-sian profile.In order tofind the phases,Seong and Peterka used an exponential autoregressive peak gener-ation model(EARPG),Y t¼aY tÀ1þ0with probability bE t with probability1Àb;0b1ð11ÞBut,Kumar and Stathopoulos used exponential peak generation model(EPG),Y t¼0with probability bE t with probability1Àb;0b1ð12Þwhere E t is an exponential distributed variable.Then,the phases of this series is/k¼tanÀ1ÀP nÀ1t¼0Y t sinð2p kt=nÞP nÀ1t¼0Y t cosð2p kt=nÞB BB@1C CC AThen,employing phases shift and using Eq.(10),the profile can be obtained.Adjusting the value of b(in EPG model)or the values of a and b(in EARPG model),the profile with given skewness and kurtosis can be generated.In this method,random phases are generated by Eqs.(11)or(12),which can generate exponential ran-dom variables only.Therefore,although this method can generate surfaces with correct spectral density or ACF,only surfaces with a certain types of skewness and kurtosis can be generated.This method also has its limitations.4.A new procedure for generating non-Gaussian surfaces4.1.One-dimensional methodFrom previous section,it is found that both current methods have their merits and their limitations.A new method,which adopts the merits of both methods is proposed.Since phase part does not affect the spectral density, Eq.(10)is adopted to generate non-Gaussian surfaces. Because Johnson translator system can form a system of distributions covering the whole skewness-kurtosis plane,it is also adopted in this new method.Similar to H uand Tonder’s method,the skewness and ku rtosis of random series can be obtained by Eqs.(6)and(7) Given skewness Sk z,kurtosis K z and given a spectral density or ACF,a new profile can be generated by the following procedure.1.First of all,find the discrete spectral density,whichwill be used to generate the profiles.If given a spec-tral density SðxÞ,S k¼S NÀk¼1SðxÞd x k¼0;1;2;...;N=2If ACF is given,Eq.(5)is used.2.Then,generate a Gaussian random series byg r¼X NÀ1k¼0exp i2p/kþkrMr¼0;1;2;...;ðNÀ1ÞThis random series g would have constant spectral density and skewness of zero and kurtosis of3.ing the skewness Sk z and kurtosis K z as the initialguess Skð1Þand Kð1Þ,find the skewness Sk g and thekurtosis K g of a non-Gaussian random series bySk g¼P mÀ1r¼0h2r3=2P mÀ1r¼0h3rSkð1ÞK g¼ðKð1ÞÞP mÀ1r¼0h2r2À6P mÀ2r¼0P mÀ1p¼rþ1h2rh2pP mÀ1r¼0h4ring the skewness Sk g and the kurtosis K g,gener-ate a new random series g0by Johnson’s distribution curve.This random series g0would have skewness of Sk g and kurtosis of K g.5.The new phase is obtained by/0k¼tanÀ1ÀP NÀ1t¼0g0r sinð2p kr=NÞP NÀ1r¼0g0r cosð2p kr=NÞB BB@1C CC Ae the new phase/0k to generate a new randomseries g00byg00r¼X NÀ1k¼0exp i2p/0kþkrMr¼0;1;2;...;ðNÀ1ÞThis random series g00would have constant spectral density.A new profile is generated byz r¼X NÀ1k¼0ffiffiffiffiffiS kpexp i2p/0kþkrMr¼0;1;2;...;ðNÀ1Þwhich is the convolution offfiffiffiffiffiS kpand g00.7.Check the skewness Skð2Þand kurtosis Kð2Þof thisnew profile.If this skewness and kurtosis do not meet our requirement,go to step3,and adjust Skð1Þand Kð1Þ,and follow the procedure3to7,till Skð2Þ%Sk z and Kð2Þ%K z.There are several numerical techniques can be used for step3to7,such as bisection method.This method will produce a profile with correct skewness,kurtosis and spectral density.In order to generate a profile with given ACF or spectral density,a random series with constant spectral density should be generated.In order to generate a profile with skewness Sk z and kurtosis K z,a random series with skewness Sk g and kurtosis K g should be generated.In step2,a random series g is generated by FFT.The spectral density of this random series g is constant.But the skewness and kurtosis of this random series g are not Sk g and K g.Hence,from step3to4, another random series g0is generated.The skewness of this random series g0is approximate Sk g,and itsJ.J Wu/Tribology International37(2004)339–346343kurtosis is approximate K g.But usually,the spectral density of this random series g0is not constant.Thus, we keep the phase of g0,but set the amplitude to be constant,and generate a new random series g00.The spectral density of this new random series g00is constant. Although its skewness and kurtosis are still not Sk g and K g,they are closer to Sk g and K g.Therefore, although the skewness and kurtosis of the generated profile are not Sk z and K z,they are close to Sk z and K z.It still seems very difficult tofind the adequate ran-dom series.Usually,the skewness and kurtosis of generated profile are different from the input skewness and kurtosis.That is,Skð1Þ¼Skð2Þand Kð1Þ¼Kð2Þ.However,we found that the generated profile is a function of input skewness Skð1Þand kurtosis Kð1Þ.It is also found that the larger is Skð1Þ, the larger is Skð2Þ;and so are Kð1Þand Kð2Þ.Thus,we can adjust Skð1Þand Kð1Þ,and obtain a profile with different S kð2Þand Kð2Þ.By using some numerical technique,we can find a profile with Skð2Þ%Sk z and Kð2Þ%K z.4.2.Two-dimensional methodThis method can be extended into two-dimension easily.The equations for generating surfaces becomesz p;q¼XMÀ1k¼0X NÀ1l¼0ffiffiffiffiffiffiffiS k;lpexp i2p/k;lþkpMþlqNp¼0;1;2;...;ðMÀ1Þq¼0;1;2;...;ðNÀ1Þð13ÞEqs.(6)and(7)becomeSk z¼P MÀ1r¼0P NÀ1s¼0h3r;sP MÀ1r¼0P NÀ1s¼0h2r;s3=2ðSk gÞð14ÞKu z¼P MÀ1r¼0P NÀ1s¼0h4r;sKu gþ6P MÀ2r¼0P MÀ1p¼rþ1P NÀ2s¼0P NÀ1q¼sþ1h2r;sh2p;qP MÀ1r¼0P NÀ1s¼0h2r;s2ð15ÞAll the other parts are the same.5.ExamplesIn this section,three non-Gaussian surfaces are simulated.First of all,a256Â256point-surface is simulated with sampling interval D x¼0:1l m and given spectral densitySðx x;x yÞ¼r2b2pðx2xþx2yþb2Þ3=2All the other statistical parameters are:r=1nm,b= 2l m,b=2.3/b,Sk=0.5,K=5.0.From Nayak[13], we know that the profile spectral density isSðxÞ¼r2bpðx2þb2ÞBy using our new method,a new surface is simu-lated.The surface is shown in Fig.3.Its height distri-bution is shown in Fig.4.The average spectral density of256profiles is shown in Fig.5.Its statistical para-meters are shown in Table1.All show good agreement between the target profile and the simulatedone.It is similar to generate a surface with given ACF.If we want to generate a 256Â256-point surface with given ACF,statistical parameters and sampling intervalas following:R ðk ;l Þ¼r 2exp f À2:3ðk 2þl 2b 2Þ12g ,r =1nm,b =2l m,Sk =0.5,K =5.0,sampling interval D x ¼0:1l mBy using the new method,a surface can be gener-ated.Its height distribution is shown in Fig.6.Its aver-age ACF of 256profiles is shown in Fig.7.Its statistical parameters are shown in Table 2.All show good agreement between the expected one and the simulated one.But not all the cases are so good as the previous examples.If a surface is generated with special skew-ness and kurtosis,the simulated surface may not be as good as expected.Now,a 256Â256point-surface with given ACF,statistical parameters and sampling intervalas following:R ðk ;l Þ¼r 2exp f À2:3ðk 2þl 2b 2Þ12g ,r =1nm,b =2l m,Sk =0,K =8.0,sampling interval D x ¼0:1l mIn this case,kurtosis is large,but skewness is zero.By using the new method,a surface is simulated.Its height distribution is shown in Fig.8.Its ACF is shown in Fig.9.Its statistical parameters are shown in Table 3.It shows good agreement between the expected one and the simulated one except the skewness.The new method can generate surfaces with better skewness,kurtosis,spectral density or ACF than other methods.But there is limitation for this new method.The new method cannot generate surfaces with every skewness and kurtosis,especially when the skewness and kurtosis are too large.It may be becausesuchTable 1Statistical parameters of case 1StdSkewness Kurtosis Target surface 1nmÀ0.5 5.0Simulation surface0.9746nmÀ0.50005.0000Table 2Statistical parameters of case 2StdSkewness Kurtosis Target surface 1nmÀ0.5 5.0Simulation surface0.9964nmÀ0.50005.0000J.J Wu /Tribology International 37(2004)339–346345surface does not exist.It may be because of the limi-tation for FFT.6.ConclusionA new procedure is proposed for generating non-Gaussian rough surfaces.This method based on FFT approach is set preserve the skewness,kurtosis together with auto-correlation function or spectral density of surfaces.Johnson translator system is used for the simu-lation of phase part.By iteration and adjusting the input skewness and kurtosis,the required non-Gaussian surfaces can be generated.Several numerical experi-ments show that this new procedure can generate the required non-Gaussian surfaces,as long as the skewness and kurtosis are not too large.In fact,our proposed method does not limit to Johnson’s translator system.If the probability of height distribution follows other systems,same procedure can be used,too.References[1]Shinozuka M,Jan C-M.Digital simulation of random processand its application.Journal of Sound and Vibration 1972;25(1):111–128.[2]Patir N.A numerical procedure for random generation of roughsurfaces.Wear 1978;47:263–77.[3]Watson W,Spedding TA.The time series modelling of non-gaussian engineering processes.Wear 1982;83:215–31.[4]Gu X,Huang Y.The modelling and simulation of a rough sur-face.Wear 1990;137:275–85.[5]YouSJ,Ehmann pu ter synthesis of three-dimensionsurfaces.Wear 1991;145:29–42.[6]Newland DE.An Introduction to Random Vibration and Spec-tral Analysis,2nd ed.London:Longman;1984.[7]Wu J-J.Simulation of rough surfaces with FFT.TribologyInternational 2000;33:47–58.[8]Hu YZ,Tonder K.Simulation of 3-D random rough surface by2-D digital filter and fourier analysis.International Journal of Machine Tools Manufacturing 1992;32(1/2):83–90.[9]Seong SH,Peterka puter simulation of non-Gaussianmultiple wind pressure time series.Journal of Wind Engineering and Industrial Aerodynamics 1997;72:95–105.[10]Kumar KS,Stathopoulos T.Synthesis of non-Gaussian windpressure time series on low building roofs.Engineering Structure 1999;21:1086–100.[11]Hill ID,Hill R,Holder RL.Fitting Johnson curves by moments.Applied Statistics 1976;25:149–76.[12]Chilamankuri SK,Bhushan B.Contact analysis of non-gaussianrandom surfaces.Proceedings of the Institution of Mechanical Engineers,Part.J 1998;212:19–32.[13]Nayak PR.Some aspects of surface roughness measurement.Wear1973;26:165–74.Table 3Statistical parameters of case 3StdSkewness Kurtosis Target surface 1nm8.0Simulation surface0.9964nmÀ0.14298.0000346J.J Wu /Tribology International 37(2004)339–346。