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

Authors: Marilena Glovnea; Sergiu Spinu; Emanuel Diaconescu
Journal: Advances in Tribology
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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 high-order bounding surface.For revolution surfaces, Shtaerman [2] proofed that a rigid punch of equation z∝r2n, 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-r2/a2, 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, P2n, to the Hertz square root generates a normal surface displacement which, within the circle, is proportional to P2n.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].Later on, these authors [6] found that circular contacts bounded by even monomial surfaces of order 2m lead to a nonuniform pressure distribution, having a local minimum in contact center, a circumferential maximum near contact boundary, and zero value in contour points. This pressure distribution was improved by imposing that the difference between current and central pressures p(r)-p0 possesses a multiple root in the origin of multiplicity order 2m-2. Because all odd order derivatives vanish in the origin, this requirement is satisfied if all even order derivatives of pressure, up to order 2m-2, are zero for r=0. These conditions yield the required central pressure plateau.This paper derives a similar procedure for elliptic elastic contacts between high-order surfaces aiming to get a flat central pressure. A generalized Hertz pressure of ordern is chosen to this end. This is defined as the product between typical Hertz square root H and an even polynomial of degree n in a-2x2+b-2y2, Pn(a-2x2+b-2y2), and of order 2n in x and y. First, a recurrence relation is derived for the coefficients of polynomial Pn which yield pressure distributions possessing a flat central plateau. Then the paper derives the surfaces of equivalent punches which generate the proposed pressure distribution when the punch is pressed against an elastic half-space. Finally, equations are derived for contact parameters, namely, contact ellipse half-axes, maximum central pressure, and normal approach or indentation of half-space.


## 2. 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:
(1)p(x,y)=p01-x2a2-y2b2∑i=0nci(x2a2+y2b2)i,c0=1.
The sum is an even polynomial of degree n in a-2x2+b-2y2, and of degree 2n in x and y. This pressure is expressed in terms of elliptic parameter ρ as follows:
(2)p(ρ)=p01-ρ2∑i=0nciρ2ic0=1.
When applied over an elastic half-space boundary, generalized Hertz pressure generates following polynomial normal displacement:
(3)w(x,y)=1πηp0b∑i=0n[ci∑k=0i(-1)kCikIk]
in which:
(4)Ik=π(2k+1)!!2k(k+1)!∑j=0k+1[(-1)jCk+1j∑l=0jC2j2l(xa)2j-2l(ya)2lIjl].Ijℓ being simple integrals depending on ellipse eccentricity:
(5)Ijl=Ijl(e)=∫0π/2sin2lγ⋅cos⁡2j-2l⁡γ(1-e2sin2γ)j+1/2dγ.
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 forceQ, 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:
(6)z(x,y)=δ-1πηp0b∑i=0n[ci∑k=0i(-1)kCikIk]=a1x2n+2+a2x2ny2+⋅⋅⋅+aN-1x2+aNy2δ 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 ai the coefficients of punch surface.Equation (6), in which Ik 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.


## 3. 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ρ:
(7)p(ρ)=p01-ρ2(1+2ρ2),forn=1,p(ρ)=p01-ρ2(1+43ρ2+83ρ4),forn=2,p(ρ)=p01-ρ2(1+65ρ2+85ρ4+155ρ6),forn=3.
Axial pressure profiles predicted by (7) exhibit a non-uniform pressure central plateau, which is unsatisfactory for both electric current density distribution and contact strength. The degree of nonuniformity increases with surface order.


## 4. 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 p0 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 ci 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:
(8)∂2p∂ρ2|ρ=0=2c1-1.
The derivative in (8) vanishes if the coefficient c1 is:
(9)c1=12.
In a similar way, the forth derivative of pressure vanishes in the origin if coefficient c2 has the following value:
(10)c2=38,
whereas
(11)c3=516,c4=35128,c5=63256
make the sixth, eighth and tenth derivatives of pressure to vanish in the origin, respectively.The above values of coefficientsci yield simply following general recurrence relation:
(12)c0=1,ci=2i-12ici-1.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:
(13)p0(ρ)=p01-ρ2,(14)p1(ρ)=p01-ρ2(1+12ρ2),(15)p2(ρ)=p01-ρ2(1+12ρ2+38ρ4),(16)p3(ρ)=p01-ρ2(1+12ρ2+38ρ4+516ρ6),(17)p4(ρ)=p01-ρ2(1+12ρ2+38ρ4+516ρ6+35128ρ8),(18)p5(ρ)=p01-ρ2(1+12ρ2+38ρ4+516ρ6+35128ρ8+63256ρ10).These pressure profiles exhibit a well-defined central pressure plateau surrounded by a monotonous decrease to zero in contour points. Moreover, asn 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.


## 5. Punch Surface and Contact Parameters

Once pressure distribution is known, the coefficientsai 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 Ik given in [6].


### 5.1. Fourth-Order Surfaces

In the caseI1=(3π/4)[I00-2(x2/a2)I10-2(y2/a2)I11+(x4/a4)I20+6(x2/a2)(y2/a2)I21+(y4/a4)I22], which means fourth order surfaces, (6) takes following form:
(19)a1x4+a2x2y2+a3y4+a4x2+a5y2=δ-1πηp0b[(1+c1)I0-c1I1].
As c1=1/2, (19) becomes
(20)a1x4+a2x2y2+a3y4+a4x2+a5y2=δ-12πηp0b(3I0-I1).
According to [6], the integrals I0 and I1 are:
(21)I0=π(I00-x2a2I10-y2a2I11),I1=3π4[I00-2x2a2I10-2y2a2I11+x4a4I20+6x2a2y2a2I21+y4a4I22].
The integrals Ijℓ are expressed in terms of complete elliptical integrals of first and second kind in [6]. By using these expressions of Ijℓ, coefficient identification in (20) leads to following equations of normal approach and surface coefficients:
(22)δ=98ηp0bK,(23)a1=ηp0b8e4a4[(2+e2)K-2(1+e2)E],(24)a2=3ηp0b4e4a2b2[(2-e2)E-2(1-e2)K],(25)a3=ηp0b8e4b4[(2-5e2+3e4)K-2(1-2e2)E],(26)a4=3ηp0b4e2a2(K-E),(27)a5=3ηp0b4e2b2[E-(1-e2)K].
Coefficients a1⋯a5 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 a1 and a3. The division of (23) and (25), member by member, yields the following transcendental equation having the eccentricity as unknown:
(28)a1a3=(1-e2)2(2+e2)K(e)-2(1+e2)E(e)(3e4-5e2+2)K(e)-2(1-2e2)E(e).
The eccentricity e is the solution of (28) which is solved numerically. Once knowing the eccentricity, remaining coefficients a2, a4, and a5 result from (24), (26), and (27). Therefore, surface polynomial possesses only two independent coefficients, in this case a1 and a3, which yield ellipse eccentricity. The remaining coefficients must have specified values given by (24), (26), (27).Force balance equation onz direction
(29)Q=∫∫Ap(x,y)dxdy
yields following relation between central pressure and normal load Q:
(30)p0=54Qπab.By substituting (30) for central pressure in (23) and (25), one finds following expressions of ellipse half-axes:
(31)a=na⋅54ηQa15,na=18πe4[(2+e2)K-2(1+e2)E]5,b=nb⋅54ηQa35,nb=β8πe4[(2-5e2+3e4)K-2(1-2e2)E]5.Maximum contact pressure, which occurs now in contact center, yields by substituting (31) for ellipse half-axes in (30):
(32)p0=npπ⋅5343Q3η2a1a35,np=1nanb.Finally, (22) leads to following expression of normal approach:
(33)δ=nδ2⋅5444η4Q4a15,nδ=94Kπna.
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 to2n+2.Punch surface is described by following equation:(34)z(x,y)=a1x4+a2x2y2+a3y4+a4x2+a5y2,
where the coefficients a1 and a3 are arbitrary imposed and a2, a4, a5 given result from (24), (26), (27).


### 5.2. Sixth-Order Surfaces

As in previous case, the coefficientsai 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:
(35)a1x6+a2x4y2+a3x2y4+a4y6+a5x4+a6x2y2+a7y4+a8x2+a9y2=δ-1πηp0b[(1+c1+c2)I0-(c1+2c2)I1+c2I2].
As c1=1/2 and c2=3/8, (35) becomes
(36)a1x6+a2x4y2+a3x2y4+a4y6+a5x4+a6x2y2+a7y4+a8x2+a9y2=δ-1πηp0b(158I0-54I1+38I2).
As established in [6], the integral I2 is given by following expression:
(37)I2=5π8(I00-3x2a2I10-3y2a2I11+3x4a4I205π8+18x2y2a4I21+3y4a4I22-x6a6I305π8-15x4y2a6I31-15x2y4a6I32-y6a6I33).The substitution of expressions (21), (37) for I0, I1, and I2, respectively, in (36) and identification of coefficients yields following results:
(38)δ=7564ηp0bK,(39)a1=ηp0b64e6a6[(8+3e2+4e4)K-(8+7e2+8e4)E],(40)a2=15ηp0b64e6a4b2[(8-3e2-2e4)E-(8-7e2-e4)K],(41)a3=15ηp0b64e6a2b4[(8-17e2+9e4)K-(8-13e2+3e4)E],(42)a4=ηp0b64e6b6[(8-23e2+23e4)E-(8-27e2+34e4-15e6)K],(43)a5=5ηp0b64e4a4[(2+e2)K-2(1+e2)E],(44)a6=15ηp0b32e4a2b2[(2-e2)E-2(1-e2)K],(45)a7=5ηp0b64e4b4[(2-5e2+3e4)K-2(1-2e2)E],(46)a8=45ηp0b64e2a2(K-E),(47)a9=45ηp0b64e2b2[E-(1-e2)K].All coefficientsa1⋯a9 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 a1 and a4. Dividing member by member equations (39) and (42), following transcendental equation having the eccentricity as unknown results:
(48)a1a4=(1-e2)3(4e4+3e2+8)K-(8e4+7e2+8)E(23e4-23e2+8)E-(8-15e6+34e4-27e2).The solution of (48) is found numerically. Once knowing the eccentricity, all remaining coefficients a2, a3, and a5 to a9 result from (40), (41) and (43)–(47). Therefore, surface polynomial possesses only two independent coefficients, in this case a1 and a4, which define ellipse eccentricity. The remaining coefficients must have specific values given by above-mentioned equations.Force balance equation onz direction yields now following relation:
(49)p0=76Qπab.By substituting (49) for central pressure in (39), (42), one finds following expressions of contact ellipse half-axes:
(50)a=na⋅76ηQa17,na=126πe6[(4e4+3e2+8)K-(8e4+7e2+8)E]7,b=nb⋅76ηQa47,nb=β26πe6[(23e4-23e2+8)E-(8-15e6+34e4-27e2)]7.Central or maximum contact pressure yields by substituting (50), for ellipse half-axes in (49):
(51)p0=npπ7565Q5η2a1a47;np=1nanb.Finally, (38) leads to following normal approach:
(52)δ=nδ27666η6Q6a17,nδ=7532Kπna.Surface equation is now straightforward:(53)z(x,y)=a1x6+a2x4y2+a3x2y4+a4y6+a5x4+a6x2y2+a7y4+a8x2+a9y2,
where, excepting the coefficients a1 and a4 which can be chosen arbitrarily, all coefficients ai result from (40), (41) and (43)–(47).


### 5.3. Eighth-Order Surfaces

As in previous cases, the coefficientsai 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 (6) becomes
(54)a1x8+a2x6y2+a3x4y4+a4x2y6+a5y8+a6x6+a7x4y2+a8x2y4+a9y6+a10x4+a11x2y2+a12y4+a13x2+⋯+a14y2=δ-1πηp0b×[(1+c1+c2+c3)I0-(c1+2c2+3c3)I1+(c2+3c3)I2-c3I3],
where I3 is given by following relation [6]:
(55)I3=35π64×(I00-4x2a2I10-4y2a2I11+6x4a4I20+36x2y2a4I21+6y4a4I22-4x6a6I30-60x4y2a6I31-60x2y4a6I32-4y6a6I33+x8a8I40+28x6y2a8I41+70x4y4a8I42+28x2y6a8I43+y8a8I44).Because coefficientsci are now c1=1/2, c2=3/8 and c3=5/16, (54) takes a simpler form, as follows:
(56)a1x8+a2x6y2+a3x4y4+a4x2y6+a5y8+a6x6+a7x4y2+a8x2y4+a9y6+a10x4+a11x2y2+a12y4+a13x2+⋯+a14y2=δ-116πηp0b(35I0-35I1+21I2-5I3).The substitution of expressions (21), (37), (55) for I0, I1, I2, and I3 respectively, in (56) and then identification of resulting coefficients yields following relations between surface coefficients ai and contact parameters a, b, e, p0, and δ:
(57)δ=12251024ηp0bK,(58)a1=53072ηp0ba8e8[(48+16e2+17e4+24e6)K53072ηp0ba8e8-(48+40e2+40e4+48e6)E],(59)a2=35768ηp0ba6b2e8[(48-16e2-9e4-8e6)E35768ηp0ba6b2e8-(48-40e2-4e4-4e6)K],(60)a3=175512ηp0ba4b4e8[(16-32e2+15e4+e6)K175512ηp0ba4b4e8-(16-24e2+4e4+2e6)E],(61)a4=35768ηp0ba2b6e8[(48-128e2+103e4-15e6)E35768ηp0ba2b6e8-(48-152e2+164e4-60e6)K],(62)a5=53072ηp0bb8e8[(48-208e2+353e4-298e6+105e8)K53072ηp0bb8e8-(48-184e2+264e4-176e6)E],(63)a6=7768ηp0ba6e6[(8+3e2+4e4)K-(8+7e2+8e4)E],(64)a7=35256ηp0ba4b2e6[(8-3e2-2e4)E-(8-7e2-e4)K],(65)a8=35256ηp0ba2b4e6[(8-17e2+9e4)K-(8-13e2+3e4)E],(66)a9=7768ηp0bb6e6×[(8-23e2+23e4)E-(8-27e2+34e4-15e6)K],(67)a10=35512ηp0ba4e4[(2+e2)K-2(1+e2)E],(68)a11=105256ηp0ba2b2e4[(2-e2)E-2(1-e2)K],(69)a12=35512ηp0bb4e4×[(2-5e2+3e4)K-2(1-2e2)E],(70)a13=175256ηp0ba2e2(K-E),(71)a14=175256ηp0bb2e2[E-(1-e2)K].All coefficientsa1⋯a14 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 a1 and a5. Dividing member by member (59)–(61) and (42), following transcendental equation having the eccentricity as unknown results:(72)a1a5=(1-e2)4[(48+16e2+17e4+24e6)K-(48+40e2+40e4+48e6)E][(48-208e2+313e4-188e6+35e8)K-(48-184e2+264e4-176e6)E].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 a1 and a5, which define ellipse eccentricity. The remaining coefficients must have specified values given by above-mentioned equations.Force balance equation onz direction yields now following correlation between central pressure and applied normal load:
(73)p0=98Qπab.Following equations for half-axesa and b result from (58), (62), and (73):(74)a=na⋅98ηQa19,(75)na=12⋅56πe8[(48+16e2+17e4+24e6)K-(48+40e2+40e4+48e6)E]9,(76)b=nb⋅98ηQa59,nb=12⋅5β6πe8[(48-208e2+313e4-188e6+35e8)K-(48-184e2+264e4-176e6)E]9.Central or maximum contact pressure yields by substituting (75) and (76) for ellipse half-axes in (73):
(77)p0=npπ⋅9787Q7η2a1a59;np=1nanbFinally, (57) leads to following normal approach:
(78)δ=nδ2⋅9888η8Q8a19;nδ=1225512Kπna.
Surface equation is now straightforward:
(79)z(x,y)=a1x8+a2x6y2+a3x4y4+a4x2y6+a5y8+a6x6+a7x4y2+a8x2y4+a9y6+a10x4+a11x2y2+a12y4+a13x2+a14y2,
where, excepting the coefficients a1 and a5 which can be chosen arbitrarily, all coefficients ai result from (59) to (61) and (63) to (71).


## 5.1. Fourth-Order Surfaces

In the caseI1=(3π/4)[I00-2(x2/a2)I10-2(y2/a2)I11+(x4/a4)I20+6(x2/a2)(y2/a2)I21+(y4/a4)I22], which means fourth order surfaces, (6) takes following form:
(19)a1x4+a2x2y2+a3y4+a4x2+a5y2=δ-1πηp0b[(1+c1)I0-c1I1].
As c1=1/2, (19) becomes
(20)a1x4+a2x2y2+a3y4+a4x2+a5y2=δ-12πηp0b(3I0-I1).
According to [6], the integrals I0 and I1 are:
(21)I0=π(I00-x2a2I10-y2a2I11),I1=3π4[I00-2x2a2I10-2y2a2I11+x4a4I20+6x2a2y2a2I21+y4a4I22].
The integrals Ijℓ are expressed in terms of complete elliptical integrals of first and second kind in [6]. By using these expressions of Ijℓ, coefficient identification in (20) leads to following equations of normal approach and surface coefficients:
(22)δ=98ηp0bK,(23)a1=ηp0b8e4a4[(2+e2)K-2(1+e2)E],(24)a2=3ηp0b4e4a2b2[(2-e2)E-2(1-e2)K],(25)a3=ηp0b8e4b4[(2-5e2+3e4)K-2(1-2e2)E],(26)a4=3ηp0b4e2a2(K-E),(27)a5=3ηp0b4e2b2[E-(1-e2)K].
Coefficients a1⋯a5 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 a1 and a3. The division of (23) and (25), member by member, yields the following transcendental equation having the eccentricity as unknown:
(28)a1a3=(1-e2)2(2+e2)K(e)-2(1+e2)E(e)(3e4-5e2+2)K(e)-2(1-2e2)E(e).
The eccentricity e is the solution of (28) which is solved numerically. Once knowing the eccentricity, remaining coefficients a2, a4, and a5 result from (24), (26), and (27). Therefore, surface polynomial possesses only two independent coefficients, in this case a1 and a3, which yield ellipse eccentricity. The remaining coefficients must have specified values given by (24), (26), (27).Force balance equation onz direction
(29)Q=∫∫Ap(x,y)dxdy
yields following relation between central pressure and normal load Q:
(30)p0=54Qπab.By substituting (30) for central pressure in (23) and (25), one finds following expressions of ellipse half-axes:
(31)a=na⋅54ηQa15,na=18πe4[(2+e2)K-2(1+e2)E]5,b=nb⋅54ηQa35,nb=β8πe4[(2-5e2+3e4)K-2(1-2e2)E]5.Maximum contact pressure, which occurs now in contact center, yields by substituting (31) for ellipse half-axes in (30):
(32)p0=npπ⋅5343Q3η2a1a35,np=1nanb.Finally, (22) leads to following expression of normal approach:
(33)δ=nδ2⋅5444η4Q4a15,nδ=94Kπna.
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 to2n+2.Punch surface is described by following equation:(34)z(x,y)=a1x4+a2x2y2+a3y4+a4x2+a5y2,
where the coefficients a1 and a3 are arbitrary imposed and a2, a4, a5 given result from (24), (26), (27).


## 5.2. Sixth-Order Surfaces

As in previous case, the coefficientsai 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:
(35)a1x6+a2x4y2+a3x2y4+a4y6+a5x4+a6x2y2+a7y4+a8x2+a9y2=δ-1πηp0b[(1+c1+c2)I0-(c1+2c2)I1+c2I2].
As c1=1/2 and c2=3/8, (35) becomes
(36)a1x6+a2x4y2+a3x2y4+a4y6+a5x4+a6x2y2+a7y4+a8x2+a9y2=δ-1πηp0b(158I0-54I1+38I2).
As established in [6], the integral I2 is given by following expression:
(37)I2=5π8(I00-3x2a2I10-3y2a2I11+3x4a4I205π8+18x2y2a4I21+3y4a4I22-x6a6I305π8-15x4y2a6I31-15x2y4a6I32-y6a6I33).The substitution of expressions (21), (37) for I0, I1, and I2, respectively, in (36) and identification of coefficients yields following results:
(38)δ=7564ηp0bK,(39)a1=ηp0b64e6a6[(8+3e2+4e4)K-(8+7e2+8e4)E],(40)a2=15ηp0b64e6a4b2[(8-3e2-2e4)E-(8-7e2-e4)K],(41)a3=15ηp0b64e6a2b4[(8-17e2+9e4)K-(8-13e2+3e4)E],(42)a4=ηp0b64e6b6[(8-23e2+23e4)E-(8-27e2+34e4-15e6)K],(43)a5=5ηp0b64e4a4[(2+e2)K-2(1+e2)E],(44)a6=15ηp0b32e4a2b2[(2-e2)E-2(1-e2)K],(45)a7=5ηp0b64e4b4[(2-5e2+3e4)K-2(1-2e2)E],(46)a8=45ηp0b64e2a2(K-E),(47)a9=45ηp0b64e2b2[E-(1-e2)K].All coefficientsa1⋯a9 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 a1 and a4. Dividing member by member equations (39) and (42), following transcendental equation having the eccentricity as unknown results:
(48)a1a4=(1-e2)3(4e4+3e2+8)K-(8e4+7e2+8)E(23e4-23e2+8)E-(8-15e6+34e4-27e2).The solution of (48) is found numerically. Once knowing the eccentricity, all remaining coefficients a2, a3, and a5 to a9 result from (40), (41) and (43)–(47). Therefore, surface polynomial possesses only two independent coefficients, in this case a1 and a4, which define ellipse eccentricity. The remaining coefficients must have specific values given by above-mentioned equations.Force balance equation onz direction yields now following relation:
(49)p0=76Qπab.By substituting (49) for central pressure in (39), (42), one finds following expressions of contact ellipse half-axes:
(50)a=na⋅76ηQa17,na=126πe6[(4e4+3e2+8)K-(8e4+7e2+8)E]7,b=nb⋅76ηQa47,nb=β26πe6[(23e4-23e2+8)E-(8-15e6+34e4-27e2)]7.Central or maximum contact pressure yields by substituting (50), for ellipse half-axes in (49):
(51)p0=npπ7565Q5η2a1a47;np=1nanb.Finally, (38) leads to following normal approach:
(52)δ=nδ27666η6Q6a17,nδ=7532Kπna.Surface equation is now straightforward:(53)z(x,y)=a1x6+a2x4y2+a3x2y4+a4y6+a5x4+a6x2y2+a7y4+a8x2+a9y2,
where, excepting the coefficients a1 and a4 which can be chosen arbitrarily, all coefficients ai result from (40), (41) and (43)–(47).


## 5.3. Eighth-Order Surfaces

As in previous cases, the coefficientsai 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 (6) becomes
(54)a1x8+a2x6y2+a3x4y4+a4x2y6+a5y8+a6x6+a7x4y2+a8x2y4+a9y6+a10x4+a11x2y2+a12y4+a13x2+⋯+a14y2=δ-1πηp0b×[(1+c1+c2+c3)I0-(c1+2c2+3c3)I1+(c2+3c3)I2-c3I3],
where I3 is given by following relation [6]:
(55)I3=35π64×(I00-4x2a2I10-4y2a2I11+6x4a4I20+36x2y2a4I21+6y4a4I22-4x6a6I30-60x4y2a6I31-60x2y4a6I32-4y6a6I33+x8a8I40+28x6y2a8I41+70x4y4a8I42+28x2y6a8I43+y8a8I44).Because coefficientsci are now c1=1/2, c2=3/8 and c3=5/16, (54) takes a simpler form, as follows:
(56)a1x8+a2x6y2+a3x4y4+a4x2y6+a5y8+a6x6+a7x4y2+a8x2y4+a9y6+a10x4+a11x2y2+a12y4+a13x2+⋯+a14y2=δ-116πηp0b(35I0-35I1+21I2-5I3).The substitution of expressions (21), (37), (55) for I0, I1, I2, and I3 respectively, in (56) and then identification of resulting coefficients yields following relations between surface coefficients ai and contact parameters a, b, e, p0, and δ:
(57)δ=12251024ηp0bK,(58)a1=53072ηp0ba8e8[(48+16e2+17e4+24e6)K53072ηp0ba8e8-(48+40e2+40e4+48e6)E],(59)a2=35768ηp0ba6b2e8[(48-16e2-9e4-8e6)E35768ηp0ba6b2e8-(48-40e2-4e4-4e6)K],(60)a3=175512ηp0ba4b4e8[(16-32e2+15e4+e6)K175512ηp0ba4b4e8-(16-24e2+4e4+2e6)E],(61)a4=35768ηp0ba2b6e8[(48-128e2+103e4-15e6)E35768ηp0ba2b6e8-(48-152e2+164e4-60e6)K],(62)a5=53072ηp0bb8e8[(48-208e2+353e4-298e6+105e8)K53072ηp0bb8e8-(48-184e2+264e4-176e6)E],(63)a6=7768ηp0ba6e6[(8+3e2+4e4)K-(8+7e2+8e4)E],(64)a7=35256ηp0ba4b2e6[(8-3e2-2e4)E-(8-7e2-e4)K],(65)a8=35256ηp0ba2b4e6[(8-17e2+9e4)K-(8-13e2+3e4)E],(66)a9=7768ηp0bb6e6×[(8-23e2+23e4)E-(8-27e2+34e4-15e6)K],(67)a10=35512ηp0ba4e4[(2+e2)K-2(1+e2)E],(68)a11=105256ηp0ba2b2e4[(2-e2)E-2(1-e2)K],(69)a12=35512ηp0bb4e4×[(2-5e2+3e4)K-2(1-2e2)E],(70)a13=175256ηp0ba2e2(K-E),(71)a14=175256ηp0bb2e2[E-(1-e2)K].All coefficientsa1⋯a14 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 a1 and a5. Dividing member by member (59)–(61) and (42), following transcendental equation having the eccentricity as unknown results:(72)a1a5=(1-e2)4[(48+16e2+17e4+24e6)K-(48+40e2+40e4+48e6)E][(48-208e2+313e4-188e6+35e8)K-(48-184e2+264e4-176e6)E].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 a1 and a5, which define ellipse eccentricity. The remaining coefficients must have specified values given by above-mentioned equations.Force balance equation onz direction yields now following correlation between central pressure and applied normal load:
(73)p0=98Qπab.Following equations for half-axesa and b result from (58), (62), and (73):(74)a=na⋅98ηQa19,(75)na=12⋅56πe8[(48+16e2+17e4+24e6)K-(48+40e2+40e4+48e6)E]9,(76)b=nb⋅98ηQa59,nb=12⋅5β6πe8[(48-208e2+313e4-188e6+35e8)K-(48-184e2+264e4-176e6)E]9.Central or maximum contact pressure yields by substituting (75) and (76) for ellipse half-axes in (73):
(77)p0=npπ⋅9787Q7η2a1a59;np=1nanbFinally, (57) leads to following normal approach:
(78)δ=nδ2⋅9888η8Q8a19;nδ=1225512Kπna.
Surface equation is now straightforward:
(79)z(x,y)=a1x8+a2x6y2+a3x4y4+a4x2y6+a5y8+a6x6+a7x4y2+a8x2y4+a9y6+a10x4+a11x2y2+a12y4+a13x2+a14y2,
where, excepting the coefficients a1 and a5 which can be chosen arbitrarily, all coefficients ai result from (59) to (61) and (63) to (71).


## 6. Discussion

An improved pressure distribution between elastic bodies bounded by high-order symmetrical surfaces is a generalized Hertz pressure in which coefficientsci result from recurrence (12).Maximum pressurep0, which occurs now in contact center, relates to average pressure Q/(πab) by following equation:
(80)p0=2n+32n+2Qπab.Naturally, ifn=0, (80) yields the well-known maximum Hertz pressure. Equations (22), (38), (57) lead to following recurrence relation for δi:
(81)δ1=98ηp0bK,δi=(2i+1)23⋅2n+1δi-1.Contact parametersa, b, δ and p0 are expressed as roots of order 2n+3 as follows:
(82)a=na⋅CnηQa12n+3,b=nb⋅CnηQan+22n+3,δ=nδ2⋅Cn2n+2η2n+2Q2n+2a12n+3,p0=npπ⋅Cn2n+1Q2n+1η2a1an+22n+3,
where Cn is a coefficient depending on n:
(83)Cn=2n+32n+2,
and na, nb, nδ, and np are coefficients which depend both on degree n and contact ellipse eccentricity. Once na is found, the remaining coefficients result from following equations:
(84)nb=β⋅an+2a12n+3na,np=1nanb,nδ1=94Kπna,nδi=(2i+1)23⋅2n+1nδi-1.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 ai yielding the equivalent punch surface need to be recomputed every time a perturbation is introduced in the input (i.e., in load, contact half-axes, or surface degree), as they depend explicitly on central pressure and on contact ellipse eccentricity, formulas (12) for ci hold in all cases.


## 7. 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 toz 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 ofO(N2) when using direct multiplication—summation can be reduced significantly to O(NlogN) if convolution is performed in the frequency domain, as element-wise 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 106 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 high-order 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 loadQ=1kN. 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 ai can be computed explicitly, using corresponding closed-form relations advanced in this work.Table 1
Formulas used to validate the predictions of the numerical program.
Central
Rigid-body
Pressure
pressure,p0
approach,δ
distribution
Forth order surface
(30)
(22)
(14)
Sixth order surface
(49)
(38)
(15)
Eighth order surface
(73)
(57)
(16)To assure that the discretization error is reduced to an acceptable level, a256×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:(85)ε=Δ∑i,j∈A|pij-pijn|Q,
where subscripts i,j are used to index the grid cells in the numerical model, Δ is the elementary patch area, pijn,i,andj∈A denote the nodal pressures predicted by the numerical program, and pij 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 Tables2, 3, and 4. The most important errors are found when assessing contact half-axis 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.Table 2
Contact parameters in case of fourth-order surface.
Contact major half-axis,
Central pressure,
Rigid-body approach,
ε, (%)
a, (m)
p
0, (Pa)
δ, (m)
Analytical
1
e
-
3 (imposed)
7.9577
e
+
8
8.3660
e
-
6
5.4011
e
-
2
Numerical
9.9844
e
-
4
7.9567
e
+
8
8.3659
e
-
6Table 3
Contact parameters in case of sixth-order surface.
Contact major half-axis,
Central pressure,
Rigid-body approach,
ε, (%)
a, (m)
p
0, (Pa)
δ, (m)
Analytical
1
e
-
3 (imposed)
7.4272
e
+
8
8.1336
e
-
6
6.0842
e
-
2
Numerical
9.9844
e
-
4
7.4268
e
+
8
8.1335
e
-
6Table 4
Contact parameters in case of eighth-order surface.
Contact major half-axis,
Central pressure,
Rigid-body approach,
ε
,(%)
a, (m)
p
0, (Pa)
δ, (m)
Analytical
1
e
-
3 (imposed)
7.1620
e
+
8
8.0065
e
-
6
6.4630
e
-
2
Numerical
9.9844
e
-
4
7.1629
e
+
8
8.0068
e
-
6Figure 1
Pressure profiles in the planey=0.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.


## 8. 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 coefficientsci 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 half-space are nonhomogeneous and contain all even power terms up to2n+2. Analytical expressions of surface coefficients involve η, p0, 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 a1 and an+2 as the independent coefficients yielding contact eccentricity. All remaining parameters result from equations proposed in the paper.Analytical expressions of contact half-axesa and b, maximum pressure p0, 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 whenn=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.



---

*Source: 832859-2012-08-16.xml*