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

                 BY A GC-DFFT TECHNIQUE
                        Spiridon Creţu1, Eduard Antaluca1 , Ovidiu Creţu2
                        Iasi Technical University, Romania, 2 Tribo-Consulting, Olympia, USA

                     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
                      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
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 
                                                                                                          ⋅ (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
                                                                                      ∑ (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
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

                                  Table 2. Summarized comparative results
                                                                           Evaluation criteria
                                                 Nr. of grid     ∆xmi,        Computing
                                                                ∆ymin ,       time               Error
                                                                [µm]             [sec]
                                                                    15 !
                         Progressive grid                                         440
         GAUSS                                      63x31           30
         (GE)                                      (1,953)         240
                          Uniform grid                                            320
                                                     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
      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
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).

To top