Improved Pressure Distribution in Elliptic Elastic Contacts between High-Order Surfaces

The improvement of mechanical contacts or microcontacts seeks a nearly uniform current density over most of contact area. When microtopography is homogeneous, this aim is achieved if nominal shape of contacting surfaces yields a nearly uniform central pressure which decreases monotonously to zero in contour points. These authors derived recently this shape for circular contacts by employing high-order surfaces. This paper extends this result to elliptical contacts. Some results are used to this end, derived for elliptical elastic contacts between high-order surfaces. As homogeneous high order surfaces lead to a highly nonuniform pressure distribution, central pressure is flattened by making the first derivatives of pressure vanish in contact center. Then, the contacts between fourth, sixth, and eighth, order surfaces are analyzed and recurrence relations for pressure distribution and contact parameters are proposed.


Introduction
Both high-load-carrying capacity of the contact and the avoidance of pressure and stress risers require a contact pressure distribution as even as possible over most of contact area.Shape improvement of circular contacts uses a pressure distribution made of a central plateau surrounded by a narrow annulus of monotonous decrease to zero [1].Such a pressure distribution yields if the equivalent rigid punch is bounded by polynomial surface of a higher order than two [1].Hertz theory fails when dealing with such surfaces.Therefore, a new Hertz kind of theory is needed for highorder bounding surface.
For revolution surfaces, Shtaerman [2] proofed that a rigid punch of equation z ∝ r 2n , pressed against an elastic half-space, generates a pressure expressed by the product between an even order polynomial in radius r, of degree 2n− 2, and typical Hertz square root √ 1 − r 2 /a 2 , a being the outer contact radius.Klubin, and later on Popov, quoted in [3], found that a pressure given by the ratio of a Legendre polynomial of order 2n, P 2n , to the Hertz square root generates a normal surface displacement which, within the circle, is proportional to P 2n .
For elliptical contact domains, Shtaerman [2] showed that a pressure written as the ratio between an even order polynomial of degree n in 1 − ρ 2 and typical Hertz square root gives rise to a surface potential expressed by an even polynomial of degree 2n in x and y when applied to an elastic half-space.This potential leads to an explicit integral expression of normal displacement in the points of bounding plane.In 1953, Galin [3] proved the following theorem: if a punch of front surface described by a polynomial of degree n is pressed against an elastic half-space over a domain bounded by an ellipse of semiaxes a and b, the contact pressure can be written as the ratio between another polynomial, of the same degree n, and Hertz square root.Gladwell [4] supplied an alternative proof of Galin's theorem for a general anisotropic half-space.For transversely isotropic half-spaces, Gladwell found a polynomial for the displacement, expressed by definite integrals, when the pressure is the ratio between associated Legendre function and typical Hertz square root.The problem of normal indentation of an elastic half-space by a rigid frictionless axisymmetric punch described by a fractional power series of radial coordinate was analyzed by Borodich [5].

Basic Equations
The solution for elliptic elastic contacts is based on results establishing the class of high-order surfaces leading to an elliptical contact area and the correlation between pressure distribution and equations of bounding surfaces [6].A generalized Hertz pressure of order n has following form: The sum is an even polynomial of degree n in a −2 x 2 + b −2 y 2 , and of degree 2n in x and y.This pressure is expressed in terms of elliptic parameter ρ as follows: When applied over an elastic half-space boundary, generalized Hertz pressure generates following polynomial normal displacement: in which: x a (4) I j being simple integrals depending on ellipse eccentricity: Equation ( 5) leads to combinations of complete elliptical integrals of first and second kind.Generalized Hertz pressure reduces to Hertz pressure when n = 0.If the pressure is generated by a rigid punch pressed against the elastic half-space by a normal force Q, the deformed half-space boundary coincides with front surface of the rigid punch in the points of contact area.According to interference equation, front surface of the punch generating the pressure given by one of (1), ( 2) is expressed by following even polynomial of order 2n + 2: δ being normal approach or half-space penetration, N = (n + 1)(n + 4)/2 maximum number of terms in an even polynomial of order 2n + 2, and a i the coefficients of punch surface.Equation (6), in which I k is given by (4), defines an explicit one to one correspondence between generalized Hertz pressure and punch surface.A generalized Hertz pressure of order n yields an even polynomial punch surface of degree 2n + 2 with respect to co-ordinates x and y, with no free term.

Contact of Homogeneous Surfaces
Applied to the contact of homogeneous even order surfaces, the above equations lead following pressure distributions expressed in terms of elliptic parameter ρ: Axial pressure profiles predicted by ( 7) exhibit a nonuniform pressure central plateau, which is unsatisfactory for both electric current density distribution and contact strength.The degree of nonuniformity increases with surface order.

Flat Central Pressure
All pressure distributions given by ( 7) possess an average negative concavity in central region.This requires important pressure maxima in peripheral zone to get a vanishing pressure in contour points.Intuitively, one feels that central pressure dimple decreases when adding lower order terms in surface equation, especially second order terms, which generate a Hertz-like pressure component compensating central crater.Of course, this means that contacting surfaces are not any longer homogeneous.Mathematically, the avoidance of peripheral pressure maxima requires central concavity of pressure distribution be zero or the difference between current generalized Hertz pressure p(ρ) and central pressure p 0 possess a multiple root in the origin of multiplicity order 2n.This means that the first 2n derivatives of pressure must be zero for ρ = 0.All odd order pressure derivatives contain the ρ as factor and consequently they vanish in the origin.Therefore, outer pressure maxima vanish if all even order derivatives of pressure up to order 2n take zero value when ρ = 0.These conditions yield the coefficients c i in (1) or (2) of pressure distribution which assure a flat central plateau.It is more convenient to use (2) because this implies fewer calculations.For instance, the second derivative of (2) with respect to elliptic parameter ρ takes in the origin following expression: The derivative in (8) vanishes if the coefficient c 1 is: In a similar way, the forth derivative of pressure vanishes in the origin if coefficient c 2 has the following value: whereas make the sixth, eighth and tenth derivatives of pressure to vanish in the origin, respectively.The above values of coefficients c i yield simply following general recurrence relation: Equation ( 12) allows writing the equation of pressure distribution of any desired order, which possesses a flat central plateau surrounded by a peripheral zone of monotonously decreasing pressure.Several such pressure distributions are: (18) These pressure profiles exhibit a well-defined central pressure plateau surrounded by a monotonous decrease to zero in contour points.Moreover, as n increases, the central flat region of pressure distribution extends, whereas maximum pressure decreases.This means that either maximum pressure decreases for the same load or load carrying capacity increases for the same maximum pressure.It is thus obvious that newly proposed pressure distributions are superior with respect to classical Hertz pressure.

Punch Surface and Contact Parameters
Once pressure distribution is known, the coefficients a i of polynomial surface of equivalent punch result by coefficient identification in (6).Several surfaces leading to flat central pressure are derived below for n = 1, n = 2 and n = 3 by using the expressions of integrals I k given in [6].

Fourth-Order Surfaces. In the case
, which means fourth order surfaces, (6) takes following form: As 4

Advances in Tribology
According to [6], the integrals I 0 and I 1 are: The integrals I j are expressed in terms of complete elliptical integrals of first and second kind in [6].By using these expressions of I j , coefficient identification in (20) leads to following equations of normal approach and surface coefficients: Coefficients a 1 • • • a 5 depend all on the eccentricity of contact ellipse.This eccentricity is found by aid of two of these coefficients, arbitrarily chosen.Following the procedure derived in [6] for homogeneous surfaces, it is convenient to find eccentricity e by involving a 1 and a 3 .The division of ( 23) and ( 25), member by member, yields the following transcendental equation having the eccentricity as unknown: The eccentricity e is the solution of (28) which is solved numerically.Once knowing the eccentricity, remaining coefficients a 2 , a 4 , and a 5 result from (24), (26), and (27).Therefore, surface polynomial possesses only two independent coefficients, in this case a 1 and a 3 , which yield ellipse eccentricity.The remaining coefficients must have specified values given by ( 24), ( 26), ( 27).Force balance equation on z direction yields following relation between central pressure and normal load Q: By substituting (30) for central pressure in ( 23) and ( 25), one finds following expressions of ellipse half-axes: Maximum contact pressure, which occurs now in contact center, yields by substituting (31) for ellipse half-axes in (30): Finally, (22) leads to following expression of normal approach: In all above equations SI system is used.As guessed initially, the surfaces of equivalent punch leading to a flat central pressure are non-homogeneous.They contain all even powers of coordinates, up to 2n + 2.
Punch surface is described by following equation: where the coefficients a 1 and a 3 are arbitrary imposed and a 2 , a 4 , a 5 given result from (24), ( 26), (27).

Sixth-Order Surfaces.
As in previous case, the coefficients a i of polynomial surface of equivalent punch leading to flat central pressure result by coefficient identification in (6).Now n = 2, which means sixth order surface, and (6) takes following form: Advances in Tribology 5 As c 1 = 1/2 and c 2 = 3/8, (35) becomes (36) As established in [6], the integral I 2 is given by following expression: (37) The substitution of expressions ( 21), (37) for I 0 , I 1 , and I 2 , respectively, in (36) and identification of coefficients yields following results: All coefficients a 1 • • • a 9 depend on the eccentricity of contact ellipse.As in previous case, this eccentricity is found by aid of two of these coefficients, arbitrary chosen.Following the procedure derived in [6] for homogeneous surfaces, it is convenient to find eccentricity e by involving a 1 and a 4 .Dividing member by member equations (39) and (42), following transcendental equation having the eccentricity as unknown results: The solution of (48) is found numerically.Once knowing the eccentricity, all remaining coefficients a 2 , a 3 , and a 5 to a 9 result from (40), ( 41) and ( 43)-(47).Therefore, surface polynomial possesses only two independent coefficients, in this case a 1 and a 4 , which define ellipse eccentricity.The remaining coefficients must have specific values given by abovementioned equations.
Force balance equation on z direction yields now following relation: By substituting (49) for central pressure in (39), (42), one finds following expressions of contact ellipse half-axes: Central or maximum contact pressure yields by substituting (50), for ellipse half-axes in (49): Finally, (38) leads to following normal approach: Surface equation is now straightforward: where, excepting the coefficients a 1 and a 4 which can be chosen arbitrarily, all coefficients a i result from (40), ( 41) and ( 43)-(47).

Eighth-Order Surfaces.
As in previous cases, the coefficients a i of polynomial surface of equivalent punch leading to flat central pressure result by coefficient identification in (6).This time n = 3, meaning an eighth order surface, and where I 3 is given by following relation [6]: Because coefficients c i are now c 1 = 1/2, c 2 = 3/8 and c 3 = 5/16, (54) takes a simpler form, as follows: a 1 x 8 +a 2 x 6 y 2 +a 3 x 4 y 4 +a 4 x 2 y 6 +a 5 y 8 +a 6 x 6 +a 7 x 4 y 2 +a 8 x 2 y 4 +a 9 y 6 +a 10 x 4 +a 11 x 2 y 2 +a 12 y 4 +a 13 The substitution of expressions ( 21), (37), (55) for I 0 , I 1 , I 2 , and I 3 respectively, in (56) and then identification of resulting coefficients yields following relations between surface coefficients a i and contact parameters a, b, e, p 0 , and δ: All coefficients a 1 • • • a 14 depend on the eccentricity of contact ellipse.As in previous cases, this eccentricity is found by aid of two of these coefficients, arbitrary chosen.Following the procedure derived in [6] for homogeneous surfaces, it is convenient to find eccentricity e by involving a 1 and a 5 .Dividing member by member (59)-( 61) and (42), following transcendental equation having the eccentricity as unknown results: The solution of (72) is found numerically.Once knowing the eccentricity, the remaining coefficients result from (40), ( 41) and ( 63)-(71).Therefore, surface polynomial possesses only two independent coefficients, in this case a 1 and a 5 , which define ellipse eccentricity.The remaining coefficients must have specified values given by above-mentioned equations.Force balance equation on z direction yields now following correlation between central pressure and applied normal load: Following equations for half-axes a and b result from (58), (62), and (73): (75) Central or maximum contact pressure yields by substituting (75) and ( 76) for ellipse half-axes in (73): Finally, (57) leads to following normal approach: Surface equation is now straightforward: z x, y = a 1 x 8 +a 2 x 6 y 2 +a 3 x 4 y 4 +a 4 x 2 y 6 +a 5 y 8 +a 6 x 6 +a 7 x 4 y 2 +a 8 x 2 y 4 +a 9 y 6 +a 10 x 4 +a 11 x 2 y 2 +a 12 y 4 +a 13 x 2 +a 14 y 2 , (78) where, excepting the coefficients a 1 and a 5 which can be chosen arbitrarily, all coefficients a i result from (59) to ( 61) and ( 63) to (71).

Discussion
An improved pressure distribution between elastic bodies bounded by high-order symmetrical surfaces is a generalized Hertz pressure in which coefficients c i result from recurrence (12).Maximum pressure p 0 , which occurs now in contact center, relates to average pressure Q/(πab) by following equation: Naturally, if n = 0, (80) yields the well-known maximum Hertz pressure.Equations ( 22), ( 38), (57) lead to following recurrence relation for δ i :

Advances in Tribology
Contact parameters a, b, δ and p 0 are expressed as roots of order 2n + 3 as follows: where C n is a coefficient depending on n: and n a , n b , n δ , and n p are coefficients which depend both on degree n and contact ellipse eccentricity.Once n a is found, the remaining coefficients result from following equations: ( As shown in [7], the newly advanced framework can also be used to get a nearly flat maximum pressure along the line contact.The severe end effects in the contact between parallel cylinders of different lengths can be attenuated by various methods such as partial crowning, hollow ended roller or generatrix profiling, the latter usually based on the logarithmic profile proposed by Lundberg [8].Reusner [9] evidenced the advantages of a new special logarithmic profile used by SKF, while Teutsch and Sauer [10] advanced a fast method for roller-race contact analysis in roller bearings, based on a theoretical and implicit load-deflection relationship. In the framework reviewed and extended herein, maximum extension of flat central pressure results if pressure coefficients for a given polynomial order are adequately chosen.Following this idea, a different approach to end effect attenuation in elastic finite length line contact between revolution bodies was advanced in [7], by allowing the roller generatrix to be a polynomial yielding a nearly flat maximum pressure along most of the contact length.A special choice of polynomial coefficients results in a contact width constant along most of the contact length, which can be interpreted as a contact area in a modified line contact.The central region of contact width increases with the degree of the surface polynomial.Eventually, the generatrix of the equivalent rigid roller, found numerically by aid of interference equation, proves to be a high-order even polynomial.For imposed load, contact area extents and degree of the polynomial, the procedure advanced in [7] establish a convenient pressure distribution in line contacts, which in its turn yields to a symmetrical polynomial roller generatrix and a central roller cross radius.While coefficients a i yielding the equivalent punch surface need to be recomputed every time a perturbation is introduced in the input (i.e., in load, contact halfaxes, or surface degree), as they depend explicitly on central pressure and on contact ellipse eccentricity, formulas (12) for c i hold in all cases.

Numerical Validation
The frictionless contact of elastic bodies assimilated to elastic half-spaces can be simulated numerically for an irregular initial clearance using the well-known algorithm based on the conjugate gradient (CG) method advanced by Polonsky and Keer, [11].The method is fast because the rate of convergence of CG is superlinear and robust because there is mathematical proof of convergence for the CG when the system matrix is symmetrical and positive definite.
The system arising from digitization of geometrical condition of deformation with respect to z direction is essentially nonlinear due to presence of rigid-body approach.Early attempts to linearize it resulted in an additional outer loop, in which normal approach was iterated with respect to static force equilibrium.Another difficulty stems from the fact that the contact area, which determines the size of the system having the nodal pressures as unknowns, is also a priori unknown, and a trial-and-error technique is required.Different techniques to overcome these difficulties are overviewed and benchmarked by Allwood, [12].
The algorithm developed by the authors is based on the work reported in [11].Linearization is achieved by assessing, in every iteration, estimates for the rigid-body translations and rotations, the latter being related to bending moments transmitted through conformal contacts [13].These estimates are numerically derived through the least square method, as the best fit of an overdetermined system of equations, assembled from equations corresponding to grid cells included in the current contact area.In order to force the solution to verify the static equilibrium equations, a correction of pressure is imposed after every iteration in the CG algorithm.The size of the system, that is, the contact area, is also adjusted at every CG cycle, according to boundary conditions expressed in terms of pressure and gap.Although the algorithm is essentially based on one level of iteration only, a reset of descent directions in the CG minimization process is required every time the contact area changes.
On every iteration of the CG, two convolution type products must be calculated.The order of computation of O(N 2 ) when using direct multiplication-summation can be reduced significantly to O(N log N) if convolution is performed in the frequency domain, as elementwise multiplication.The errors introduced by problem periodization, implicitly assumed when the spatial series of pressure and influence coefficients are transferred to frequency domain via fast Fourier transform are avoided with a minimal additional computational cost using the Discrete Convolution Fast Fourier Transform (DCFFT) technique advanced by Liu et al. [14].The computer code obtained using the CG combined with the DCFFT technique is fast enough to solve grids of up to 10 6 points in the contact area in a reasonable amount of time, allowing for simulation of deterministic rough contact scenarios.
An implementation of this algorithm is used herein to validate the framework for the smooth contact between highorder surfaces.The initial contact clearance fully authorizes the use of half-space assumption, allowing for the use of appropriate Green functions when computing free-surface deflections.
All simulations are performed with a normal load Q = 1 kN.The data used for validation was generated by imposing a major half-axis of contact ellipse of 1 mm and an aspect ratio of contact ellipse β = 0.5, using appropriate formulas advanced in this work (see Table 1).With these parameters fixed (but otherwise arbitrarily chosen), contact geometry was readily available as input for the numerical program, as all terms a i can be computed explicitly, using corresponding closed-form relations advanced in this work.
To assure that the discretization error is reduced to an acceptable level, a 256 × 256 uniformly spaced rectangular grid, usually assumed for smooth contact scenarios, was imposed in a surface domain exceeding with 20% symmetrically the contact axes predicted by the analytical framework.The imposed precision for pressure convergence was fixed at ε 0 = 256 −3/2 , according to numerical experimentations reported in [11].
The relative error of pressure distribution was calculated using the following formula: where subscripts i, j are used to index the grid cells in the numerical model, Δ is the elementary patch area, p n i j , i, and j ∈ A denote the nodal pressures predicted by the numerical program, and p i j discrete pressures computed using the newly advanced closed-form relations, at coordinates matching grid cells control points.
The predictions of the numerical program are found to agree well with the closed-form relations advanced in the theoretical framework, as shown in Tables 2, 3, and 4. The most important errors are found when assessing contact halfaxis and can be attributed to grid resolution, as the numerical prediction for this parameter can only vary with integers of grid steps.The analytical and numerical solutions for pressure profiles along contact major half-axis are depicted in Figure 1.Dimensionless pressure is defined as ratio to Hertzian pressure corresponding to the imposed load and contact area.
As shown in [15], a central plateau of uniform pressure can also be found in the elastic-plastic nonconforming contact, when the residual term of displacement and of subsurface stresses becomes significant.The residual print due to permanent deformation of the surface, acting together with the hardening of the elastic-plastic material, decrease the central pressure computed according to the elastic model.The elastic-plastic material responds to loading by developing residual stresses, which decrease stresses induced by pressure, as to oppose further plastic yielding.In the framework proposed herein, the flattened pressure plateau results in purely elastic conditions, and is related only to the fine-tuning of the surface profile.It is clear that the newly advanced indenter will accommodate larger loadings prior to yield inception compared to the quadratic one, when the same contact area is established.

Conclusions
Pressure distribution in elastic elliptic contacts between high-order homogeneous surfaces is nonuniform, possessing a local minimum in contact center, maxima in peripheral region, and zero value in contour points, thus reducing the load carrying capacity of the mechanical contact.A more advantageous pressure distribution must possess a flat central region surrounded by a peripheral monotonous decrease to zero.Such pressure distribution yields from a generalized Hertz pressure if coefficients c i result from a recurrence equation.As the degree of generalized Hertz pressure increases, the central plateau of pressure extends and maximum pressure decreases, giving overall a better pressure distribution.
The surfaces of equivalent rigid punch generating these favorable pressures when pressed against an elastic halfspace are nonhomogeneous and contain all even power terms up to 2n + 2. Analytical expressions of surface coefficients involve η, p 0 , a, b, and algebraic combinations of complete elliptical integrals of first and second kind.Only two of these coefficients can be independent to get an elliptical contact area.Although these can be arbitrarily chosen, it is convenient to take a 1 and a n+2 as the independent coefficients yielding contact eccentricity.All remaining parameters result from equations proposed in the paper.
Analytical expressions of contact half-axes a and b, maximum pressure p 0 , and normal approach δ are derived by similarity to Hertz equations for fourth, sixth, and eighth order surfaces.These are expressed as roots of order 2n + 3 of involved parameters and they can be applied directly to calculate contact elements.
General or recurrence equations are established for main contact parameters, which reduce to the results derived in here when n = 1, n = 2, and n = 3.The newly advanced formulas are verified well by a numerical program for the elastic contact with known, but otherwise arbitrarily distributed initial clearance, giving confidence in the proposed framework.

Table 1 :
Formulas used to validate the predictions of the numerical program.

Table 2 :
Contact parameters in case of fourth-order surface.

Table 3 :
Contact parameters in case of sixth-order surface.

Table 4 :
Contact parameters in case of eighth-order surface.