VIEWS: 21 PAGES: 9 POSTED ON: 5/26/2012 Public Domain
NATIONAL THE ANNALS OF UNIVERSITY 39 TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 THE STUDY OF NON-HERTZIAN CONCENTRATED CONTACTS BY A GC-DFFT TECHNIQUE Spiridon Creţu1, Eduard Antaluca1 , Ovidiu Creţu2 1 Iasi Technical University, Romania, 2 Tribo-Consulting, Olympia, USA spcretu@mec.tuiasi.ro ABSTRACT In the last two decades the title problem has been solved by discretizing the virtual contact domain as a finite union of rectangles with non-overlapping interiors, both the surface pressures and tangential tractions being taking constant over each rectangle. On the other hand, to analyze the contact between rough surfaces a number of grid points greater than 105 are needed, so that, new numerical methods have been evolved to speed up the computational process. In this respect, the paper uses the Conjugate Gradient Method (CGM) to solve the set of linear equations for unknown pressures. The Discrete Convolution Fast Fourier Transform (DC-FFT) is further involved in each iteration of CGM to multiply both the pressure vector and the direction vector with the matrix of the influence coefficients. The final comparison points out the superiority of the CG+DC-FFT method, as long as, the computing time are considered as the comparison criterion. The internal stresses field has been further obtained for the general case when both, the normal pressure distribution and the corresponding tangential friction traction were present on the contact area. KEYWORDS: concentrated contacts, non-Hertzian, conjugate gradient, discrete convolution, fast Fourier transform. 1. INTRODUCTION amount of the computing time for N greater than 5000. To achieve the maximum reduction in size and To solve a non-conforming concentrated contact increased efficiency, Creţu [1996, 2002], used a non- with smooth surfaces the classic approach uniform mesh, the finest mesh being in the critical zo- approximates the surfaces around contact region as ne where the pressure concentration occurs (Fig. 1.). quadratic. In this assumption is possible to use the direct equations supplied by Hertz theory of the elastic contact. Unfortunately, most of the machine elements that operate under severe contact condition, such as rolling bearings, gears, cams, variable drives, usually achieve a non-hertzian contact. Some numerical approaches of the problem had been developed during the last decades: Hartnett [1977], Webster and Sayles [1986]. Both, the pressure distribution and the contact area calculation are based on Boussinesq integral equation, Johnson [1985]. The virtual contact area is divided into rectangles; the Fig. 1 Stress raisers at the ends of the real line deformation in a particular rectangle is simply the contacts. sum of the effects of the pressure of each rectangle. The system of equation formed by using the Furthermore the working surfaces of machine geometric equation of elastic contact and the force elements are always rough. Asperities and the balance equation is solved by Gauss elimination interaction among asperities significantly affect the method. Because the computation time for Gauss contact performances. The detailed behaviour of elimination is O (N3), this approach needs a huge rough surfaces in contact may only be determined by 40 NATIONAL THE ANNALS OF UNIVERSITY TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 numerical solution for samples of digitised real pressure acting at point j is expressed by using the surfaces. The Nyquist sampling theorem requires that influence coefficient Cij: the product of the sampling interval and the maximum 1 −ν I 1 − ν II 2 2 possible frequency of the waveform of the real u z (i , j ) = π ⋅E + π ⋅E ⋅ C ij ⋅ p j (1) surface should be smaller than 0.5. To account for all I II the roughness components that are meaningful in friction and wear a sufficiently long sample length must be considered. Consequently, the number of data in the roughness sample will be extremely large. In this respect, modern equipment for 3-D surface topography analysis typically produces data sets about 500 by 500 points. There are a number of approaches to contact analysis of real rough surfaces. In the last decade, new methods for the numerical solution of rough contact problems have emerged that are orders of Fig, 2 Non-uniform mesh (arithmetic progression). magnitude faster that the conventional methods based on Gauss elimination. The influence coefficient Cij is given by the Brand and Lubrecht [1990] developed multilevel- expression: multisummation (MLMS) technique, that was used [ Cij = f (xl , yl ) + f (xr , yr ) − f (xl , yr ) − f (xr , yl ) (2) ] further by Lubrecht and Ioannides [1991] in their where: study on dry contact, but the applicability of their iteration scheme to rough contact problems appears to xl = x − ∆x j , xr = x + ∆x j , x = xi − x j (3) be limited. yl = y − ∆y j , yr = y + ∆y j , y = yi − y j Ren and Lee [1993] used their Moving Grid Method (MGM) to simulate the 3-D contact of rough ( ) f (x, y) = x ⋅ ln y + x2 + y2 + y ⋅ ln x + x2 + y2 ( ) surfaces. The MGM reduces the required computer The total displacement at i is therefore expressed by: memory size but has no advantage over the N 1 1 − ν I 2 1 − ν II 2 N ⋅ (C ⋅ p ) (4) conventional matrix inversion method in reducing the z ∑ u (i) = u (i, j) = ⋅ j =1 z + π EI EII ∑ j =1 ij j computing time. Tian and Bhushan [1996] developed a method Considering the geometric form of the elastic contact based on a variational principle so that the equation, the following system of linear algebraic convergence of the numerical solution is guaranteed. equation is developed: But the computing time of the method is greatly 1 −ν I 2 2 N ∑ (C ) 1 1 − ν II dependent on the number of initial contact points and ⋅ + ⋅ ij ⋅ p j (5) π E I E II proved considerably time-consuming when the j =1 contact of a sphere on a plane was simulated. = δ − h (i) A number of researchers applied Fast Fourier Equation (5) represents a set of N simultaneous Transform (FFT) to contact problems. The 2-D line equations for unknown pressures p. In matrix form contact problem was examined by Ju and Farris the equations become: [1996], while 3-D point contact problem was studied C⋅ p=B (6) by Nogi and Kato [1997], Ai and Sawamiphakdi [1999], Polonsky and Keer [1999, 2000]. The equilibrium condition gives the last equation: N 2. INFLUENCE COEFFICIENTS ∑ (p j =1 j ⋅ dx j ⋅ dy j = F ) (7) A plane rectangular domain larger than the The domain Dc of the contact points represents a expected contact area is divided into rectangular discrete analogue of the actual area of contact. But Dc elements, the pressure over each element being is not known in advance and needs to be found. In this considered constant. purpose two supplementary restrictions have to be To improve accuracy without increasing considered: computing time both, the length and width of the pi > 0, ∀ I ∈ Dc (8a) rectangular patches may be established non-uniform, (Fig. 2), being finest in the most interested zone, pi = 0, ∀ I ∉ Dc (8b) where the stress concentrator acts. The effect of an uniform pressure distribution 3. GAUSS ELIMINATION METHOD acting on a rectangular area has been investigated by Love and the necessary equations are given in classic The algebraic system represented by (6) and (7) is contact books as Johnson [1985]. The normal solved by Gauss elimination method (GEM). It is deflection of particular point i due to constant NATIONAL THE ANNALS OF UNIVERSITY 41 TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 found that values of pr near the periphery of the virtual contact domain, Dv, are negative which implies that a tensile traction is needed to maintain the contact at these points. A new iteration is necessary to consider these points having the uniform pressure pr=0. Repeated iterations converge to a set of values pi which are positive or zero. The boundary between the region with pi >0 and the region with pr =0 defines the contact area, to the accuracy of the mesh size. A cylindrical roller with circular profile has been chosen to compare the results obtained with this method with those obtained with the more time consuming finite element method (FEM). The contact geometry and the finite element results are those given by J. de Mul [1986], figure 3 and Table1. In both cases the number of points was 63x31 = 1953, the non-uniform mesh being obtained by using arithmetic progressions, (Fig. 4). The comparison presented in figure 5 points out the superiority of the non-uniform mesh upon the uniform one. In fact the non-uniform mesh gave the same variation as obtained with FEM, but the uniform mesh was not Fig. 3 The geometry of the bodies in contact able to catch the real stress concentration value. (Jan de Mull, [1985]). Table 1 The contact geometry ( Jan de Mull,[1985]) Dw Lw R1 R2 ax Dci [mm] [mm] [mm] [mm] [mm] [mm] 15 16 1114 1.006 6.994 58.5 Fig. 4 3-D contact pressure distributions obtained by using GEM with non-uniform mesh (63x31 = 1953points). Fig. 5 2-D contact pressure distributions obtained with uniform and progressive grids by using GEM. 42 NATIONAL THE ANNALS OF UNIVERSITY TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 Fig. 6 2-D contact pressure variations along roller centerline for uniform grid using CG M. 4. CONJUGATE GRADIENT METHOD, for further iteration with respect to the contact (CGM) area; - ii the force balance equation is enforced Roughness samples with very large number of during each iteration for the contact pressure. data points need to be used to analyze rough contact. The rigid deflection δ is not obtained in the Consequently, solution of a rough contact problem by iteration process but it may easily be determined, as the Gauss elimination method, may take days, even long as the pressure distribution was obtained. The on a modern computer with more than 2GHz working comparative results are depicted in figure 3 for three frequency. different mesh grids: 63x31=1,953 grid points; Several iterative methods as Jacobi, Gauss-Seidel 255x31=7,905 grid points and 512x31=15,872 grid and Conjugate Gradient are available to solve the points. The necessity of a fine mesh to catch the stress system (3) and (4). From these, the Conjugate concentration is clearly pointed out. Gradient Method (CGM) has important advantages: - there exists a rigorous mathematical proof of 5. THE FAST FOURIER TECHNIQUE the method convergence; APPLIED TO CONTACT PROBLEMS - it offers a very high rate of convergence (super-linear); Although, the computation of the elastic - it requires very modest storage capacity that is displacements by using the influence coefficients extremely advantageous when very large considerably reduces the computational work, still systems of equations are involved. remains the time consuming computation represented The CGM is an iterative method that generates a by the multiplication operations between the matrix C sequence of approximations of the solution starting with the pressure vector p (Eq. 4), as well as, the from an arbitrary initial approximation. The direction vector t (Eq. 9). These multiplication recurrence formulae of the CGM are: operations require O (N2) operations, and if N is r ⋅t larger than 10,000 the computing time increases at p I +1 = p I + I I ⋅ t I (9a) unacceptable values for current analyses. rI ⋅ c ⋅ t I The linear elasticity has a convolution nature that rI ⋅ t I is clearly pointed out in the Boussinesq’s equation: rI +1 = rI − ⋅ C ⋅ tI (9b) rI ⋅ c ⋅ t I 1−ν 2 p(ξ,η) rI +1 ⋅ rI +1 t I +1 = rI +1 + ⋅ tI (9c) uz (x, y,0) = ∫∫ πE D (x −ξ)2 + ( y −η)2 dξdη (10) rI ⋅ rI This observation directs the computational where i = 0,1,2..; rI and tI are vectors of N elements; algorithm to the application of the Fourier Transform P0 is an arbitrary start vector and t0=r0=B - C·P0. technique to solve the corresponding convolution Polonsky and Keer [1999, 2000], designed a fast products. However, the direct application of FFT iteration scheme based on the conjugate gradient me- approach to non-periodic problems, such as the thod (CGM) and a multi-level multi-summation me- concentrated contacts are, introduces an error known thod (MLMSM) to solve the concentrated non- as the periodicity error. To reduce the periodicity hertzian contact problem .The iteration process error, the virtual contact domain has to be extended proposed by Polonsky and Keer and adopted into the far beyond the actual contact domain that can have current work has some distinctive particularities as the adversary effect to reduce the efficiency of the follows: FFT implementation into the contact numerical - i the contact area is established in the course algorithm. of the pressure iteration, so that, there is no need NATIONAL THE ANNALS OF UNIVERSITY 43 TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 The FFT-based methods fall into two categories: When the DC-FFT technique is used to solve a -i the continuous convolution and fast Fourier convolution product the number of operations is transform (CC-FFT) method, that uses the reduced to O (Nw·log2NW), where frequency response function (FRF) and needs a Nw=(mfft*mfft)*(N1*N2) is the number of mesh computational domain four to eight times as large points involved in the windows of DC-FFT as the target domain; operation. -ii the discrete convolution and fast Fourier For the contact conditions mentioned in figure 3, transform (DC-FFT) that uses either frequency res- the pressure distribution and actual contact area ponse function (FRF) or influence coefficients (IC). obtained with CGM+DC-FFT are presented in figure In order to increase the efficiency of the 7. numerical algorithm developed, a dedicated real Because the extension of the virtual domain discrete fast Fourier transform routine for 3-D contact might severely reduce the computational efficiency of problems has been developed and incorporated into DC-FFT, a comparative study has been developed, to the code. The conversion of the linear convolution validate both, the accuracy and computational speed into a cyclic convolution appeared as compulsory superiority. preliminary treatments. In this respect, three techniques have been used: 6. COMPARATIVE NUMERICAL i - the extension of the initial mesh of the RESULTS contact area in each direction by a certain factor denoted mfft. A comparative numerical study has been carried ii - the wrap-around order of the influence out to point out the main goals of the methods coefficients incorporated in the presented. The comparison criteria have been as (mfft*N1)*(mfft*N2) extended mesh. follows: the finesse of the grid, the computing time, iii - the pressure vector p is maintained in the and the relative error. original domain (0<i<N1, 0<j<N2) and is zero- The comparative results are presented in figure 8 padded in all other zones of the extended domain. and summarised in Table 2. Fig. 7 3-D pressure distribution and actual contact area obtained by using CGM+DC_FFT (512x512 =262,144 grid points). 44 NATIONAL THE ANNALS OF UNIVERSITY TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 Fig. 8 2-D pressure distributions on centerline near stress concentrator for FEM, CGM and CGM+FFT methods. Table 2. Summarized comparative results Evaluation criteria Nr. of grid ∆xmi, Computing Method points ∆ymin , time Error [µm] [sec] 15 ! Progressive grid 440 GAUSS 63x31 30 (GE) (1,953) 240 Uniform grid 320 50 63x31 240 5 10-3 (1,953) 50 CONJUGATE GRADIENT 127x31 119 23 10-3 (CG) (3,937) 50 255x31 60 95 10-3 (7,905) 50 256x32 60 <1 10-3 (8,192) 50 CONJUGATE GRADIENT + 128x128 120 3 10-3 DISCRETE FAST FOURIER (16,384) 12 TRANSFORM 256x256 60 12 10-3 (CG + DFFT) (65,536) 6 512x512 30 ! 50 10-3 (262,144) 3 NATIONAL THE ANNALS OF UNIVERSITY 45 TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 7. SUB-SURFACES STRESSES clearly appeared, as long as, the computing time was taken as comparison criterion. 7.1. Normal loading 5. Further, the internal stresses field has been numerically obtained for the general case when Having established a set of solutions for the both, the normal pressure distribution, and the uniform pressure steps pj, it is straightforward to corresponding tangential friction traction were compute the sub-surface stresses due to each step, and present. from the superposition define the full sub-surface 6. The superposition of tangential friction traction stresses. The DC-FFT technique has been once again over the normal loading increases the maximum involved in any convolution products computation. value of von Mises stress and moves its position closer to the surface. 7.2. Normal load and surface shear traction REFERENCES The authors also considered the sub-surface 1. Ai X., Sawamiphadi R., 1999, ”Solving Elastic Contact stress distribution for a surface shear traction with Between Rough Surfaces as an Unconstrained Strain Energy various friction coefficients. The relevant equations Minimization by Using CGM and FFT Techniques”, Trans. of the are from Hill, 1993. Setting qj = µpj and adding these ASME, Journal of Tribology, vol. 121, pp 639-637. 2. Brandt A., Lubrecht A. A., 1990, ”Multilevel Matrix stresses to those produced by the pure normal loading Multiplication and Fast Solution of Integral Equations”, Journal of allows the definition of the sub-surface stresses due to Computational Physics, vol. 90, pp. 348-370. a normal and tangential loads to be established for a 3. Bryant G., Keer L. M.,1985, “ “ sliding surface interface with a given coefficient of 4. Creţu S., 1996, ”Initial Plastic Deformation of Cylindrical Roller Generatrix: Stress Distribution and Fatigue Life Tests”, Acta friction µ. Tribologica, 4, No. 1-2, pp. 1-6 The DC-FFT technique has been once again 5. Creţu S., 2002, Mecanica Contactului, vol. I, Ed. ”Gheorghe involved in any convolution products computation. Asachi”, Iaşi. A comparison of the present method with those 6. Hartnett M. J., 1979, ”The Analysis of Contact Stresses in Rolling Element Bearings”, Trans. of the ASME, Journal of of existing analytical stress solutions for two smooth Tribology, vol. 101, pp 105-109. spheres in contact with friction traction has been done 7. Hill D. A., Nowell D., Sackfield A., 1993, Mechanics of to verify the sub-surfaces numerical scheme. Elastic Contacts, Butterworth,Oxford In figure 6 the von Mises stresses are plotted as 8. Johnson K. L., 1985, Contact Mechanics, Cambridge University Press. functions of the depth along the x=y =0 line for three 9. Ju Y., Farris T. N., 1996 ”Spectral Analysis of Two- different coefficients of friction values. The dimensional Contact Problems”, Trans. of the ASME, Journal of distributions presented in figure 9 agree entirely with Tribology, vol. 118, pp 320-328. corresponding analytical distributions obtained by 10. Lee G., Keer L. M., 1994,” “ 11. Liu S., Wang Q., 2002, ”Studying Contact Stress Fields Bryant and Keer, (1985), as well as with those Caused by Surface Tractions with a Discrete Convolution and Fast obtained by Lee and Keer, (1994). Fourier Transform Algorithm”, Trans. of the ASME, Journal of The von Mises stresses plotted in the yoz plane as Tribology, vol. 124, pp 36-45. equal values curves are depicted in figures 10…13, 12. de Mul J. M., Kalker J. J., Fredriksson B., 1986, ”The Contact between Arbitrarily Curved Bodies of Finite Dimensions”, for the cases fy =0 , fy =0.1 and fy =0.25 respectively. Trans. of the ASME, Journal of Tribology, vol. 108, pp 140-148. 13. Nogi T., Kato T., 1997, ”Influence of a Hard Surface Layer 8. CONCLUSIONS on the Limit of Elastic Contact – part I: Analysis Using a Real Surface Model”, Trans. of the ASME, Journal of Tribology, vol. 119, pp 493-500. 1. For the case of non-Hertzian concentrated 14. Polonsky I. A., Keer L. M., 1999, ”A Numerical Method For contacts three numerical approaches have been Solving Contact Problems Based On The Multilevel Multisumation developed to obtain the actual contact area and and Conjugate Gradient Techniques”, Wear, 231, pp. 206-219. 15. Polonsky I. A., Keer L. M., 2000, ”Fast Methods For the pressure distribution on it. Solving Rough Contact Problems: A Comparative Study”, Trans. 2. When the purpose is to point out the position and of the ASME, Journal of Tribology, vol. 122, pp 36-41. the value of the stress concentration, the using of 16. Ren N., Lee S. C., 1993, “Contact simulation of 3- a non-uniform grid offers the possibility to dimensional rough surfaces using moving grid methods”, Trans. of the ASME, Journal of Tribology, vol. 115, pp 597-601. maintain the number of grid points lower than 17. Sayles R. S.,1996, “Basic Principles of Rough Surface 5000. In this case the Gauss Elimination Method Contact Analysis Using Numerical Methods”, Tribology is recommended to be used. International, vol. 29, No. 8, pp639-650. 3. When the purpose is the study of the roughness 18. Tian X., Bushan B., 1996, ”A Numerical Three Dimensional Model For Rough Surfaces by Variational Principle”, Trans. of the influence, a uniform mesh is evidently needed ASME, Journal of Tribology, vol. 118, pp 33-42. and the number of the grid points must be as 19. Webster M. N., Sayles R. S., 1986, ”A Numerical Model for great as 2.5·105. The fast numerical algorithms Elastic Frictionless Contact of Real Rough Surfaces”, Trans. of the are essential to be involved in order to reduce the ASME, Journal of Tribology, vol. 108, pp 314-320. computing time by three or even four magnitude order. 4. The superiority of the Conjugate Gradient Method coupled with Fast Fourier Technique 46 NATIONAL THE ANNALS OF UNIVERSITY TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 Fig. 9 Von Mises stresses along the x = y =0 line for three different coefficients of friction. Fig. 10 Equal values von Mises plots for pure normal loading. NATIONAL THE ANNALS OF UNIVERSITY 47 TRIBOLOGY “DUNĂREA DE JOS“ OF GALAŢI CONFERENCE FASCICLE VIII, TRIBOLOGY 24-26 September 2003 2003 ISSN 1221-4590 Fig. 11 Equal values von Mises plots for normal and tangential friction loading, (fy =0.1). Fig. 12 Equal values von Mises plots for normal and tangential friction loading, ( fy =0.25).