Numerical Simulation of 3D Rough Surfaces and Analysis of Interfacial Contact Characteristics

Yang Wang

1　Introduction

In recent decades,the study about the contact characteristics of rough surfaces on micro-scale especially in terms of three dimensional attracted great attention because it has been aware of a key factor determining the static and dynamic thermal performance of machine systems.It is proved that the numerical simulation of rough surfaces is not only a more efficient method to generate rough surfaces with the inputs determined by different height and spatial statistical parameters[Chilamakuri and Bhushan(1998)],but also an optimal method to capture the contact characteristics comprised of the load-deformation relationship,the contact pressure distribution,the real contact area and the interfacial loading/unloading characteristics of rough interfaces.

Thus far,some studies have been accomplished on generating various Gaussian and non-Gaussian rough surfaces and analyzing the contact characteristics of rough surfaces by utilizing different methods.For instance,Patir(1978)used the linear transformation of random sequence and proposed a method for generating rough surfaces with prescriptive ACFs.His scheme,however,can be impractical for its time consuming in solving a series of non-linear equations.Wu(2000)and Newland(1984)proposed numerical methods for generating Gaussian rough surfaces with applying FFT method and the prescriptive power spectral density(PSD)or ACF.Reizer(2011)simulated measured rough surfaces with Wu’s and New land’s FFT methods respectively,and then conducted comparisons between the statistical characteristics of the measured and simulated rough surfaces.Watson,King,Spedding and Stout(1979)and Watson and Soedding(1982)applied time series models to simulate Gaussian and non-Gaussian engineering surfaces.Seong and Peterka(1997)and Suresh Kumar and Stathopoulos(1999)used exponential autoregressive peak generation model(EARPG)and exponential peak generation model(EPG)to find the phases for generating non-Gaussian line pro files.However,the formulas for the random phases only apply to exponential random variables,so only rough surfaces with certain skewness and kurtosis values can be generated[Wu(2004)].According to[Hill,Hill and Holder(1976);Nagahara(2003)],non-Gaussian translator systems have been developed to generate non-Gaussian rough surfaces.Be-sides,Wu(2004)de fined one-dimensional phase sequences to generate updated white noise and non-Gaussian height sequences,which provides a basis for the generation of the 2D phase sequences in this study to make the average ACFs of the simulated rough surfaces coincide well with the input ones.As for the speci fic simulation process,Greenwood and Williams(1966)simplified the process by using the frictionless contact of an elastic hem is phere and a rigid flat plane to stochastically model an entire contacting surface with a postulated Gaussian height distribution.The aforementioned GW model assumed that the asperities of the surfaces deformed independently and the substrate below the asperities did not deform at all.Supplementing the GW model,many elasto-plastic asperity models have been devised[Chang,Etsion and Bogy(1987);Majumdar and Bhushan(1991);Jackson and Green(2006)].Meanwhile,the asperity interactions were considered in[Zhao and Chang(2001);Gao,Bower,Kim,Lev and Cheng(2006);Ciavarella,Greenwood and Paggi(2008);Buczkowski and Kleiber(2009)]where the asperities were simpli fied geometrically for normal contact and the contact of two rough surfaces was also simplified to the contact between one rough surface and one smooth rigid plane.Although these above models have been proven feasible,most of them were proposed based on many assumptions and simplifications of the shape,magnitude and height distribution of asperities,which can lim it their applications and affect the solution accuracies to some extent.Thus,to solve this problem,Jackson and Green(2005)proposed a 2D axisymmetric FE model of an elastic-perfectly plastic hemisphere in contact with a rigid flat surface.Admittedly,FEM is an effective way to overcome the defects of other analytical methods,but it calls a great challenge to control mesh quality and solution ef ficiency.Sahoo and Ghosh(2007)and Hyun,Pei,Molinari and Robbins(2004)simulated the interface of one rigid smooth plane and one self-affine rough surface based on fractal theory.However,there are disputes concerning whether the fractal theory can adapt to generate any engineering rough surface.Meanwhile,the above simplification of the interface may not reveal the true contact characteristics perfectly.

In this study,a simple and practical scheme is proposed to accurately generate Gaussian and non-Gaussian rough surfaces and comprehensively analyze the contact characteristics of rough surfaces.In the first section,rather than traditional measuring methods where the collected data are inaccurate and redundancy due to the measuring equipment and the choice of measuring area,the statistical characteristics consisted of amplitude and spatial characterization parameters are utilized here to effectively capture the detailed surface topography data[Wu(2000);Wu(2004);Bakolas(2003)].Next,the FFT method is adopted in the simulation process to improve the calculation efficiency since it is relatively faster and more convenient when compared with time series tool as a common used method.Be-sides,another calculation strategy is also proposed for generating non-Gaussian sequences directly from white noise sequences with the input statistical parameters and ACFs.More than that,this strategy is not only able to generate rough surfaces with any skewness and kurtosis in the whole skewness-kurtosis plane,but also narrow down the errors generated during the filter process by eliminating this complicated procedure,which is a great improvement over the previous methods.In addition,after the generation of non-Gaussian rough surfaces,array restructure and rearrangement methods,substituting for the conventional intricate linear transformation and filter process aimed to remain the satisfied first four moments unchanged,are introduced to ensure that the generated Gaussian and non-Gaussian rough surfaces have almost the same ACFs with the input ones.Last but not the least,based on the simulated models,the quantitative study,more than the qualitative study,is also carried on through FEM about the contact characteristics,such as load-deformation relationship,the contact pressure distribution,the real contact area and the interfacial loading/unloading characteristics,to reveal the superiority of the proposed models as well as lay a solid foundation for further study on surface topography optimization.

2　Numerical simulation of rough surfaces

2.1　Statistical characteristics

The surface topographies of rough interfaces on micro-scale are so irregular that can only be well described by the statistical characteristics of which the parameters are calculated with references to the height distribution,probability density function(PDF),power spectral density(PSD)or the autocorrelation function(ACF).

2.1.1Amplitude characterization

Height distribution function and the corresponding PDF play an important part in deciding most of the statistical height parameters of rough surfaces.In most cases,Gaussian distribution is assumed for simplification,but in fact,most machined rough surfaces are non-Gaussian[Watson,King,Spedding and Stout(1979)].For a surface described by a continuous univariate PDFf(z),thevthcentral momentuvcould be expressed as follows:

where E is the expectation operator;ˉzis the average height of the rough surface.Generally speaking,ˉzis defined to be zero at the beginning of numerical simulation,and finally modified to the input mean valueμby addingμto the height sequencez.

According to Eq.(3)and Eq.(4),the first central momentu1is zero;the second one is equal to the varianceσ2,whereσis the standard deviation of the height sequencez;Whenv≥3,thevthcentral momentuvcan also be standardized through being divided byespecially,the skewness and kurtosis are defined to be the 3rdand 4thstandard moments,respectively.

For a rough surface described by a discrete height sequencez(m,n),thevthcentral momentuvis as follows:

To sum up,the first four standard moments can be described as follows:

wheremandnstand for the numbers of row and column,respectively.Whenzhas a mean of zero and a variance of 1,SkzandKuzare equivalent to the 3rdand 4thoriginal moments,respectively.

2.1.2Spatial characterization

The statistical parameters of spatial characteristics are mainly determined by ACF or PSD,and the PSD can be obtained from the ACF function with fast Fourier transform(FFT)method.For a continuous height functionz(x,y),the normalized ACF can be described as follows:

whereLxandLyare the lengths of rough surface in x and y directions,respectively;τxandτyare the distance in x and y directions between any two points in sample area withis defined as a discrete sequence with spacing ofand whenmandnapproach∞,the corresponding ACF can be expressed as Eq.(10),especiallyR(0,0)=σ2.

wherek=0,1,2,...,m-1;l=0,1,2,...,n-1.

Nowadays,PSD has become a widely accepted method for characterizing the spatial frequencies of rough surfaces which can provide spectrum features of spatial frequencies and contains the information about spatial frequencies from 0 to∞(frequencies 0 and∞stand for the surfaces of in finite flat and ones with in finite small asperities,respectively).In this study,the PSD can be derived from the height information(continuous height distribution functionz(x,y)or discrete height sequencez(m,n))or the ACF of rough surfaces with FFT method as Eq.(11)or Eq.(12).

wherer=0,1,...,m-1,s=0,1,...,n-1;fxandfystand for the frequencies in x and y directions with

2.2　Simulations of Gaussian rough surfaces

As shown in Fig.1,the simulation of Gaussian rough surfaces is the process of producing 2D phase sequences from white noise sequences and then updating for generating Gaussian rough surfaces.In order to make the target Gaussian rough surfaces satisfy the input ACF,μ,σ,SkzandKuz,the discrete method for the input ACFs and the definition of phase sequences as well as the computation of the PSD constants of white noise sequences are presented in detail.The concrete steps for the simulation of Gaussian rough surfaces are listed as follows:

Figure 1:Numerical simulation of Gaussian rough surfaces.

1.The discrete form of ACFR(m,n)is extracted from the given ACFf(x,y).(1)The input ACF should be de fined at first.For example,the following exponential form of ACFf(x,y)is widely quoted.

whereσis the standard deviation of height sequence;βxandβystand for the autocorrelation lengths inxandydirections,respectively;φstands for the prescriptive orientation of surface texture.

(2)An equal discrete spacing inxandydirections needs to be speci fied(for example,and then,discretize the given ACF into sequenceR(m+1,n+1)within the range ofand

wheremandnare numbers of surface points inxandydirections;g=-m/2,-m/2+1,...,m/2;h=-n/2,-n/2+1,...,n/2.

(3)One of the repeated rowsm/2+1 orm/2+2,and the repeated columnsn/2+1orn/2+2 in the sequenceR(m+1,n+1)is required to be deleted to generate the sequenceR(m,n)of the same size with the target rough surface.

2.The PSD and TF can be obtained fromR(m,n).

(1)FFT method is applied here to get the PSDP(I,J)according to Eq.(12).

(2)Since the PSD of white noise is a constantC,assumingC=1,the transfer functionH(m,n)can be obtained:

3.The phase sequence Φ(m,n)can be produced from the white noise sequenceη(m,n)which needs to be updated with inverse fast Fourier transform(IFFT)method.

(1)The white noise generator randn(m,n)is used to obtain the white noise sequenceη(m,n).

(2)The phase sequence Φ(m,n)is generated from the sequenceη(m,n)with Eq.(16).

where all the elements in the phase sequence Φ range from 0 to 2π,and the computation of Φ(m,n)can be conducted with FFT method to improve the efficiency.(3)The white noise sequenceη(m,n)is updated from Φ(m,n)with IFFT method:

(4)FFT method is then introduced to the updated the white noise sequenceη(m,n).

4.The Gaussian sequencez0(m,n)is generated from the dot-product ofAandHwith IFFT method.

5.Repeat steps 3～4,if the skewness and kurtosis values ofz0approach 0 and 3,respectively,otherwise,repeat the step 2.

6.The Gaussian sequencez0is required to be scaled and translated to obtain the target Gaussian height sequence.

(1)According to Eq.(2),the meanμ0and standard deviationσ0ofz0can be calculated.

(2)The sequencez0is scaled towhose standard deviation is equal toσ.At the same time,the PSD(equal to the constant C)of the updated white noise sequence as well asHcan be updated:and

(3)The height sequencez1is translated in the height direction:zG=z1-μ0+μ,to obtain the target Gaussian sequencezGwith the mean ofμand the standard deviation ofσ.

The Gaussian rough surfaces produced with aforementioned simulation method are able to meet the requirements of the height statistical characteristics and ACFs,and moreover,can provide the spatial distribution rules for the corresponding non-Gaussian ones investigated in the follow ing section..

2.3　Simulations of non-Gaussian rough surfaces

There are usually two steps in the simulation of target non-Gaussian rough surfaces including the non-Gaussian transformation and the linear transformation[Wu(2004);Hill,Hill and Holder(1976);Nagahara(2003);Hu and Tonder(1992);Bakolas(2003)]Non-Gaussian height sequences are produced through non-Gaussiantransformation with the transitional skewnessSkηand kurtosisKuηdetermined from Eq.(20),and subsequently filtered to satisfy the input first four moments and ACFs.

where SkzandKuzare the input skewness and kurtosis values;θi=h(k,l)whereh(k,l),working as the filter coefficient,can be obtained from the transfer function(TF),andk=1,2,...,m,l=1,2,...,nas well asi=(k-1)m+l.Because the calculation of the summation ofitems is too time consuming,Bakolas(2003)simplified the above computation of the kurtosis from Eq.(21)to Eq.(22)

AlthoughSkzandKuzare within the whole skewness-kurtosis plane1≥0),the linear transformation as shown in Eq.(20)may enableSkηandKuηto deviate from the corresponding skewness-kurtosis plane:So the method for generating non-Gaussian rough surfaces cannot cover the whole skewness-kurtosis plane.

In this end,an optimal non-Gaussian translator system is introduced in this section to transform varied white noise sequences directly into non-Gaussian height sequences with the accurate input first four moments as shown in Fig.2.The transformed sequences,however,are completely disorganized and that’s why the reconfiguration is an indispensable procedure for generating the target non-Gaussian height sequences with the prescriptive ACF according to the spatial distribution characteristics of the corresponding Gaussian rough surfaces.

2.3.1Choice of non-Gaussian translator systems

Nowadays,two non-Gaussian transformation methods,Johnson’s and Pearson’s translator systems,are accepted by most academics.Both of them can achieve the aim of generating the sequences of height distributions covering the whole skewness-kurtosis plane:but only one of them will be utilized here to generate non-Gaussian sequences in the terms of transformation accuracy and efficiency.Johnson system transforms a random Gaussian sequence to a non-Gaussian sequence with given mean,standard deviation,skewness and kurtosis by using three types of fitting methods,SB,SUand SL.

The lognormal distribution SL:

The bounded distribution SB:

Figure 2:Numerical simulation of non-Gaussian rough surfaces.

The unbounded distribution SU:

whereηis the white noise sequence andz2is the non-Gaussian sequence;γandδare shape parameters;λandξare proportional coef ficient and position parameter,respectively.The value of all the parameters includingγ,δ,λandξcan be determined by the given skewness and kurtosis.Through the transformation,any non-Gaussian rough surface with given skewness and kurtosis can be obtained.

Johnson’s translator system tends to use theSLmethod when neitherSUnorSBcan converge to a constant,which may give rise to the reduction of transformation accuracy.Pearson’s non-Gaussian translator system of frequency curves uses a probability functionp(x)satisfying the follow ing differential equation Eq.(26):

wherea,c0,c1andc2are constants that can be obtained with given standard de-viation,skewness,kurtosis and the probability function derived from the above differential equation.

Pearson’s system introduces some widely used distributions including gamma,beta,F-distribution,t-distribution and so on to ensure that their first four moments coincide well with the prescriptive ones.However,its utility is limited because the PDF to be used must be chosen from the aforementioned distribution functions[Bakolas(2003)].In this study,Johnson’s translator system is given the priority to be considered first during the simulation of non-Gaussian rough surfaces.Nevertheless,when SUand SBin Johnson’s translator system are unable to converge,Pearson’s translator system will take the place to improve the accuracy and efficiency of non-Gaussian transformation.

2.3.2Non-Gaussian transformation

Johnson’s or Pearson’s translator system can directly transform Gaussian or white noise sequences to non-Gaussian sequences,respectively,with given values ofμ,σ,SkzandKuz.In this study,the complicated filter process is eliminated and height sequences with any inputSkzandKuzin the whole skewness-kurtosis plane(Kuz-can be obtained by carrying out the follow ing steps.

1.The white noise generator randn(m,n)is utilized to obtain the white noise sequenceη(m,n)and transform it to a non-Gaussian sequencez2(m,n)with Johnson’s translator system in terms of the input first four moments;If SBor SUare unable to converge,use the Pearson’s translator system to obtain the non-Gaussian sequencez2(m,n),instead.

2.The skewness and kurtosis of the sequencez2(m,n)is calculated according to Eq.(7)and Eq.(8)and repeat steps 1～3 if the accuracy cannot meet the satisfaction.

3.The sequencez2(m,n)needs to be scaled and translated withz3=σ·z2/σ1-μ1+μ(σ1andμ1stand for the standard deviation and the mean of the sequencez2)to update the non-Gaussian sequence to satisfy the given first four moments.

As mentioned above,non-Gaussian height sequences can be transformed from white noise sequences conforming strictly to the input first four moments.The accuracies of the skewness and kurtosis can be improved by the circulation of transformation process and the mean and standard deviation can be strictly guaranteed by the translating and scaling methods.After the aforementioned non-Gaussian transformation,the non-Gaussian height sequences with given height statistical characteristics will be generated.

2.3.3Restructure and rearrangement of the height sequences

The above section generates the non-Gaussian sequencez3(m,n)whose first four moments accurately approach those of the target height sequence.Afterwards,the restructure and rearrangement methods of height sequences are utilized to make the target non-Gaussian height sequence satisfy the given ACFs inxandydirections on the basis of the spatial distribution of the corresponding Gaussian rough surfaces.The following steps can be adopted:

1.The simulated Gaussian sequence zG(m,n)and the non-Gaussian sequencez3(m,n)need to be restructured into sequencesq(m·n,1)andQ(m·n,1),respectively.

2.The sorting and indexing tools in Matlab are introduced to arrange the sequenceq(m·n,1)and Q(m·n,1)then,rearrange the sequence Q(m·n,1)in terms of the sorting index of the sequenceq(m·n,1).

3.The rearranged sequence Q(m·n,1)is required to get restructured into the target non-Gaussian height sequencezN(m,n)which has almost the same ACFs with those of the corresponding Gaussian height sequencezG.

The proposed restructure and rearrangement methods,substituting for the conventional complicated linear transformation and filter process to remain the satisfied first four moments unchanged,will make the generated Gaussian and non-Gaussian rough surfaces have nearly the same ACFs with the input ones.So the non-Gaussian rough surfaces with the prescriptive first four moments and ACFs could be accurately simulated.

2.4　Simulation cases

Both of Gaussian and non-Gaussian rough surfaces have isotropic and anisotropic,geometrical and physical characteristics.As shown in Fig.3 and Fig.4,the simulated isotropic Gaussian and non-Gaussian rough surfaces as well as the anisotropic ones are obtained.It is obvious that isotropic rough surfaces have the same autocorrelation lengths with each other while anisotropic rough surfaces have different autocorrelation lengths inxandydirections.The simulation time of GS 1-1,NGS 1-2,GS 2-1,NGS 2-2,GS 3-1,NGS 3-2,GS 4-1 and NGS 4-2 are 0.7,0.8,0.5,1.3,2.5,3.3,2.6 and 3.3 seconds,respectively.Especially,it has been noticed that decreasing the discrete point count or the accuracy requirements of the skewness and kurtosis values means less time to simulate the rough surfaces with the same statistical parameters.

Figure 3:Isotropic Gaussian and non-Gaussian rough surfaces.

Figure 4:Anisotropic Gaussian and non-Gaussian rough surfaces.

2.4.1Statistical height characteristics

The height distribution characteristics of the first group rough surfaces(GS 1-1 and NGS 1-2)are shown in Fig.5 where the normal distribution plot indicates that NGS 1-2 has a clear deviation from the linear normal distribution.Similarly,the height distributions of the other simulated non-Gaussian rough surfaces can be obtained as shown in Fig.6 where it is obvious that the probability distribution of the rough surfaces with negative skewness skew to negative side of the normalized height,and vice versa.

Figure 5:Height distribution characteristics of GS 1-1 and NGS 1-2.

As common used roughness parameters,the pro file arithmetic meanRaand the root-mean-square(rms)roughnessRq(especially,Rq=σ)are calculated.As shown in Fig.3 and Fig.4,the ratio ofRatoRqof Gaussian rough surfaces approximately equals 0.8 while that of non-Gaussian rough surfaces deviates from 0.8.As for the accuracies of the first four moments,μandσcan be kept consistent with the output ones by the scaling and translating methods and the errors of the skewness and kurtosis values listed in Table 1 can be reduced sharply by the elimination of filter process.

Figure 6:Height distributions of other non-Gaussian surfaces.

Table 1:Errors of the skewness and kurtosis values.

The results show that the skewness and kurtosis errors are less than 0.03 and 0.05,respectively.Comapred with the results in[Reizer(2011)],the accuracy of the skewness and kurtosis values is greatly improved.Moreover,with more time given to simulation,the accuracy of the skewness and kurtosis values could be further improved.

2.4.2ACFs of rough surfaces

The normalized average ACFs of the simulated rough surfaces inxandydirections are extracted from the spatial characteristics based on the method with reference to[Wu(2004)].As shown in Fig.7 and Fig.8,the ACFs of the non-Gaussian rough surfaces coincide well with those of the corresponding Gaussian rough surfaces.What’s more,it can also be found that whenm·?x/βx≥20 orn·?y/βy≥20,the ACFs of simulated rough surfaces inxorydirection agree with the corresponding theoretical ones.

Figure 7:ACFs of isotropic group surfaces in x or y directions.

2.4.3Level-spacing distribution of asperities

Figure 8:ACFs of anisotropic group surfaces in x and y directions.

Eight-point method is introduced in this section to find the asperity peaks and valleys of rough surfaces.Specifically speaking,a pointKi,jin a rough surface can be found surrounded by eight points as shown in Fig.9.Once the height coordinate of pointKi,jis larger or smaller than those of the eight surrounding points,the pointKi,jcan be taken as an asperity peak or valley.Thus,based on the eight-point method,all of the asperity peaks and valleys of simulated non-Gaussian rough surfaces could be found out as shown in Fig.10 where"+"and"o"stand for the peaks and the valleys,respectively.Meanwhile,the spacing characteristics of asperity peaks inxandydirections are extracted as shown in Table 2.Moreover,the spacing distribution of the asperity peaks and valleys can be adjusted by the discrete spacing(?xand ?y)and the corresponding autocorrelation length inxandydirections(βxandβy).

Table 2:Spacing characteristics of asperity peaks.

2.4.4Texture direction of rough surface

The proposed scheme is suitable for generating Gaussian and non-Gaussian rough surfaces with isotropic and linear anisotropic surface texture.Among the anisotropic rough surfaces,as shown in Fig.4,GS 3-1 and NGS 3-2 withφ=0 andβx<βypresent longitudinal surface texture,while GS 4-1 and NGS 4-2 withφ=0 andβx>βyshow obvious transverse surface texture.If the orientation angleφof the rough surface NGS 4-2 is changed from 0 to 45?or 120?,considering that the orientation angleφequals the rotation angle ofzaxis,the corresponding surfaces will be generated as shown in Fig.11.

Figure 9:Eight-point method for peaks and valleys.

Figure 10:Asperity peaks"+"and valleys"o"of NGS 1-2.

Using the aforementioned simulation method,Gaussian and non-Gaussian rough surfaces can be generated with any skewness and kurtosis values in the whole skewness-kurtosis planewhich help lay a solid foundation for the following section about the study of the contact characteristics of two non-Gussian surfaceas.

Figure 11:Straight surface texture with different directions.

3　Analysis of contact characteristics of rough interfaces

3.1　Establishment of contact models

Two pairs of Gaussian height sequences generated with the above mentioned methods are utilized to establish two contact models with the application of the corresponding spacing of?xand ?yin ANSYS.In the two contact models,?x=?y=1μm,and the bottom and up contact rough surfaces have the same standard deviation,i.e.,σ1=σ2,hereσ1=σ2=1.6 in 1#contact andσ1=σ2=0.8 in 2#contact.During the modeling,one of the contact bodies should be turned upside down and translated to be in contact with the other one because of the different normal directions of the two contact rough surfaces.

Besides,some detailed information for modeling is listed as follows:

1.The contact models are established from bottom to up,i.e.,from points to rough surfaces to contact bodies which are meshed with Solid45 hexahedral grids.

2.All the contact bodies are supposed to adopt the same material with an elastic modulus of 71GPa and a Poisson’s ratio of 0.33;

3.Von Mises plasticity criterion and multi-kinematic hardening rule are introduced with taking the yield strength(σs=497Mpa)into consideration;

4.For convenience,the stress-strain curve is de fined exactly the same with that in[Oskouei,Keikhosravy and Soutis(2009)].

The constructed contact models are shown in Fig.12.During the establishment,conta173 and targe170 elements are utilized to model the contact pairs.Afterwards,Solid45 hexahedral elements are adopted to mesh all the contact bodies.In this study,each contact rough surface is simulated with 128×128 points,so the number of solid element which is about 129 thousand(127×127×4×2)will be within the calculation ability of common computers and FE software.Moreover,reasonable constraints are applied to improve the solution convergence.Specifically speaking,normal displacement d inzdirection,instead of the normal force or pressure,is applied at the top surface of the substrate to adjust the clamping force while the bottom surface of the substrate is fixed.Therefore,only from the reaction force or the contact force at the interface,the corresponding equivalent normal force can be derived.

Figure 12:Establishment of the contact model.

3.2　Resu lts and discussion

Using the above established contact models,this section is aimed at revealing the contact characteristics comprised of the load-deformation relationship,the contact pressure distribution,the real contact area and the loading/unloading characteristics at different load levels based on the finite element analysis of 1#and 2#contacts.

3.2.1Load-deformation relationship

First of all,P0is defined as the nominal contact pressure which is equivalent to the quotient of the corresponding contact forceFat the interface to the nominal contact areaA0,that isP0=F/A0.Then,the elastic-plastic deformation law of rough interface can be elucidated by applying the normal displacementdlittle by little from a small value.Meanwhile,through the post-processing of ANSYS,Von M ises equivalent stressσeand the corresponding deformation can also be derived.As shown in Fig.13,some asperities show the tendency to yield even under a very small load(P0=3.44 MPa),and with the increment ofP0from 3.44 MPa to 141.23 MPa,the asperities yield increasingly from locally to wholly.

Figure 13:Von Mises equivalent stress of the bottom contact body in 1#contact.

Since the deformation of the thin substrate is so small that could be ignored,the normal deformation of the contact interface is considered to be equal to the normal displacementd.From the data of the normal deformation and the equivalent loadP0,the load-deformation curve can be obtained as shown in Fig.14.It’s obvious that under the circumstance of the same normal deformation,theP0in 1#contact will be less than that in 2#contact.The slope of the deformation-load curve,enbodying the contact stiffness of the interfaces,increases withP0ord,and it changes slowly untilP0approaches 75MPa.Moreover,under the same load ofP0,the rougher the contact interface,the smaller the slope of the deformation-load curve.Thus,it is verified that the improvement of the sufaces quality and the increment of the assembly load can be beneficial to enhance the contact stiffness.

Figure 14:Deformation-load curve of different interfaces.

3.2.2Interface pressure distribution and real contact area

The interface pressure distribution can also be extracted by the post-processing of ANSYS.As shown in Fig.15,the maximum of contact pressurePis up to 3GPa even under a small normal load ofP0=3.44MPa,and it hardly change with the variation ofP0from 3.44MPa to 141.23MPa.

The real contact area of the interface is defined as the area with positive contact pressure which canalso be obtained from the post-processing of ANSYS.Asshown in Fig.16,the ratio ofAr/A0of real contact presents a nonlinear increment with the rise of the normal displacementdwhile a nearly linearly one with the rise of the nominal pressureP0.WhenP0<100Mpa,the ration ofAr/A0is less than 10%,and it is still less than 15%even under a load ofP0=180MPa.Besides,the ratio of Ar/A0can also be affected by the roughness of the interface,for example,the ratio ofAr/A0in 1#contact is smaller than that in 2#contact.

3.2.3Loading/unloading characteristics of the contact interfaces

Based on the above analysis of contact characteristics,this section is aimed to reveal the characteristics of interfaces during the loading/unloading processes in the aspect ofP0andAr/A0as shown in Fig.17 and Fig.18.

Figure 15:Interface pressure distribution under different load in 1#contact.

Figure 16:variation of Ar/A0 with the normal load.

As shown in Fig.17,the interfacial normalized deformationd/σ′is introduced,whereWhend/σ′≤1.5,the two loading curves coincide well with each other whileP0of 2#contact is smaller than that of 1#contact whend/σ′>1.5.Judging from the change of the curve slope,the conclusion that the interface stiffness can be considerably increased by the process of loading and unloading would be reached due to the fact that the slope of unloading curve is obviously larger than that of loading one.Since the characteristics of interfaces during loading and unloading processes in this study is analogous to the deformation-load characteristics of the gaskets in[Murali Krishna,Shunmugam and Siva Prasad(2007)],the normal stress,strain and other mechanical behaviors of the interfaces on micro-scale can also be simulated by adopting the gasket elements with a certain thickness to achieve the cross-scale coupling between micro-scales and macro-scales.

Figure 17:The relationship between d/σ′and P0 during loading and unloading process.

Figure 18:The relationship between Ar/A0 and P0 during loading and unloading process.

As shown in Fig.18,the ratio ofAr/A0of 1#and 2#contact show a linear increment with the rise ofP0during loading and an apparent non-linear relationship during unloading,which is similar with the results from[Pei,Hyun,Molinari and Robbins(2005)],and thus,veri fies the accuracy of the analysis in this study.

4　Conclusions

In this study,a novel scheme is proposed for numerically generating Gaussian and non-Gaussian rough surfaces with isotropic and linear anisotropic surface texture.Afterwards,the contact characteristics including the load-deformation relationship,the contact pressure distribution,the real contact area and the loading/unloading characteristics at different load levels are analyzed.Finally,several conclusions can be drawn as follows:

1.Three-dimensional rough surfaces with any skewness and kurtosis in the whole plane can be numerically generated without complicated linear transformation and filter process of non-Gaussian rough surfaces.Meanwhile,the errors of the skewness and kurtosis obtained from the simulated rough surfaces are less than 0.03 and 0.05,respectively,which is a great improvement compared with the outcome of previous methods

2.Phase sequences are produced from white noises sequences with higher calculation efficiency due to the application of FFT method,and the first four moments of generated Gaussian and non-Gaussian rough surfaces are controlled to closely approach the input ones through the numerical operations like translating,scaling and reconfiguring.

3.The proposed method for generating Gaussian and non-Gaussian rough surfaces would provide massive and accurate surface topography data with different statistical parameters for the study on the contact characteristics of rough surfaces.

4.The simulated models are established by FEM with a high-quality meshing method of hexahedral elements and reasonable constraints which give rise to the improvement of solution accuracy and efficiency.Meanwhile,the contact characteristics including load-deformation relationship,the contact pressure distribution,the real contact area and the interfacial loading/unloading characteristics of rough interfaces are analyzed,which could illuminate the further study of surface topography optimization.

Acknowledgement:This project is supported by the National Natural Science Foundation of China[51405377],the Natural Science Foundation of Hunan Province,China[2015JJ2054],and the Research Foundation of Education Bureau of Hunan Province,China[14C0435].

New land,D.E.(1984):An Introduction to Random Vibration and Spectral Analysis.Longman.

Bakolas,V.(2003):Numerical generation of arbitrarily oriented non-Gaussian three-dimensional rough surfaces.Wear,vol.254,pp.546-554.

Buczkowski,R.;K leiber,M.(2009):Statistical Models of Rough Surfaces for Finite Element 3D-Contact Analysis.Arch.Comput.Method.E.,vol.16,pp.399-424.

Chang,W.R.;Etsion,I.;Bogy,D.B.(1987):An Elastic-Plastic Model for the Contact of Rough Surfaces.J.Tribol.-T ASME,vol.109,pp.257-263.

Chilamakuri,S.K.;Bhushan,B.(1998):Contact analysis of non-Gaussian random surfaces.P.I.Mech.Eng.J-J.Eng.,vol.212,pp.19-32.

Ciavarella,M.;Greenwood,J.A.;Paggi,M.(2008):Inclusion of"Interaction"in the Greenwood and Williamson Contact Theory.Wear,vol.265,pp.729-734.

Gao,Y.F.;Bower,A.F.;Kim,K.S.;Lev,L.;Cheng,Y.T.(2006):The Behavior of an Elastic-Perfectly Plastic Sinusoidal Surface under Contact Loading.Wear,vol.261,pp.145-154.

Greenwood,J.A.;William s,J.B.(1966):Contact of Nom inally Flat Surfaces.Proc.R.Soc.London,Ser.A,vol.295,pp.300-319.

Hill,I.D.;Hill,R.;Holder,R.L.(1976):Fitting Johnson curves by moments.Appl.Stat.,vol.25,pp.180-189.

Nagahara,Y.(2003):Non-Gaussian filter and smoother based on the Pearson distribution system.J.Time Ser.Anal.,vol.24,pp.721-738

Hu,Y.Z.;Tonder,K.(1992):Simulation of 3-D random rough-surface by 2-D digital- filter and Fourier analysis.Int.J.Mach.Tool.Manu.,vol.32,pp.83-90.

Hyun,S.;Pei,L.;M olinari,J.F.;Robbins,M ark O.(2004):Finite-Element Analysis of Contact between Elastic Self-Af fine Surfaces.Phys.Rev.E,vol.70,pp.1-12.

Jackson,R.L.;Green,I.(2005):A Finite Element Study of Elasto-Plastic Hem ispherical Contact against a Rigid Flat.J.Tribol.-T ASME,vol.127,pp.343-354.

Jackson,R.L.;G reen,I.(2006):A Statistical Model of Elasto-Plastic Asperity Contact between Rough Surfaces.Tribol.Int.,vol.39,pp.906-914.

Majum dar,A.;Bhushan,B.(1991):Fractal Model of Elastic-Plastic Contact between Rough Surfaces.J.Tribol.-T ASME,vol.113,pp.1-11.

Murali Krishna,M.;Shunmugam,M.S.;Siva Prasad,N.(2007):A study on the sealing performance of bolted flange joints with gaskets using finite element analysis.Int.J.Pres.Ves.Pip.,vol.84,pp.349-357.

Oskouei,R.H.;Keikhosravy,M.;Soutis,C.(2009):Estimating Clamping Pressure Distribution and Stiffnessin Aircraft Bolted Jointsby Finite-Element Analysis.P.I.Mech.Eng.G-J.Aer.,vol.223,pp.863-871.

Patir,N.(1978):A numerical procedure for random generation of rough surfaces.Wear,vol.47,pp.263-277

Pei,L.;Hyun,S.;M olinari,J.F.;Robbins,M.O.(2005):Finite element modeling of elasto-plastic contact between rough surfaces.J.Mech.Phys.Solids,vol.53,no.11,pp.2385-2409.

Reizer,R.(2011):Simulation of 3D Gaussian surface topography.Wear,vol.271,pp.539-543.

Sahoo,P.;Ghosh,N.(2007):Finite Element Contact Analysis of Fractal Surfaces.J.Phys.D:Appl.Phys.,vol.40,pp.4245-4252.

Seong,S.H.;Peterka,J.A.(1997):Computer simulation of non-Gaussian multiple w ind pressure time series.J.Wind Eng.Ind.Aerod.,vol.72,pp.95-105.

Suresh Kumar,K.;Stathopoulos,T.(1999):Synthesis of non-Gaussian w ind pressure time series on low building roofs.Eng.Struct.,vol.21,pp.1086-1100.

Watson,W.;King,T.G.;Spedding,T.A.;Stout,K.J.(1979):The machined surface—time series modeling.Wear,vol.57,pp.195-205.

Watson,W.;Soedding,T.A.(1982):The time series modelling of non-Gaussian engineering processes.Wear,vol.83,pp.215-231.

Wu,J.J.(2000):Simulation of rough surfaces with FFT.Tribol.Int.,vol.33,pp.47-58.

Wu,J.J.(2004):Simulation of non-Gaussian surfaces with FFT.Tribol.Int.,vol.37,pp.339-346.

Yu,N.;Polycarpou,A.A.(2004):Combining and Contacting of Two Rough Surfaces with Asymmetric Distribution of Asperity Heights.J.Tribol.-T ASME,vol.126,pp.225-232.

Zhao,Y.W.;Chang,L.(2001):A Model of Asperity Interactions in Elastic-Plastic Contact of Rough Surfaces.J.Tribol.-T ASME,vol.123,pp.857-864.