Subsurface stresses in an elliptical Hertzian contact

Abstract

The traditional solution for the stresses below an elliptical Hertzian contact expresses the results in terms of incomplete Legendre elliptic integrals, so are necessarily based on the length of the semi-major axis a and the axis ratio k. The result is to produce completely different equations for the stresses in the x and y directions; and although these equations are now well-known, their derivation from the fundamental, symmetric, integrals is far from simple. When instead Carlson elliptic integrals are used, they immediately match the fundamental integrals, allowing the equations for the stresses to treat the two semi-axes equally, and so providing a single equation where two were needed before. The numerical evaluation of the Carlson integrals is simple and rapid, so the result is that more convenient answers are obtained more conveniently. A bonus is that the temptation to record the depth of the critical stresses as a fraction of the length of the semi-major axis is removed. Thomas and Hoersch’s method of finding all the stresses along the axis of symmetry has been extended to determine the full set of stresses in a principal plane. The stress patterns are displayed, and a comparison between the answers for the planes of the major and minor semiaxes is made. The results are unchanged from those found from equations given by Sackfield and Hills, but not previously evaluated. The present equations are simpler, not only in the simpler elliptic integrals, but also for the “tail” of elementary functions.

Introduction

The analysis of elliptical Hertzian contacts is most conveniently done using Love’s potential function ¹ $\phi \equiv {\int }_{\gamma }^{\infty }(1-\frac{x^2}{a^2+s}-\frac{y^2}{b^2+s})\mathrm{ds}$ , where $\gamma$ is the positive root of $\frac{x^2}{a^2+\lambda }+\frac{y^2}{b^2+\lambda }+\frac{z^2}{\lambda }=0$

Thus, from Johnson’s Contact Mechanics, ² the displacements due to a Hertzian pressure distribution $p=p_0{[1-x^2/a^2-y^2/b^2]}^{1/2}$ are $w=\frac{ab{p}_0}{2E^{\prime }}{\int }_{\gamma }^{\infty }\left(1-\frac{x^2}{a^2+s}-\frac{y^2}{b^2+s}\right)\frac{ds}{{[s(a^2+s)(b^2+s)]}^{1/2}}$ $\frac{\mathrm{ds}}{{[s(a^2+s)(b^2+s)]}^{1/2}}$ ,

Within the loaded region we require $w=x^2/2R_1+y^2/2R_2$ , so that (since within the loaded region $\gamma =0$ ), we have

$$\begin{array}{c}\multicolumn{1}{c}{\frac{1}{2R_1}=\frac{\mathrm{ab}{p}_0}{2E\prime }{\int }_0^{\infty }\frac{\mathrm{ds}}{s^{1/2}{(a^2+s)}^{3/2}{(b^2+s)}^{1/2}}\mathrm{and}}\\ \multicolumn{1}{c}{\frac{1}{2R_2}=\frac{\mathrm{ab}{p}_0}{2E\prime }{\int }_0^{\infty }\frac{\mathrm{ds}}{s^{1/2}{(a^2+s)}^{1/2}{(b^2+s)}^{3/2}}.}\end{array}$$

Those with considerable mathematical dexterity can transform these integrals into Legendre elliptic integrals ( $K(e)-E(e)$ and $E(e)-(1-e^2)K(e)$ (where $e^2=1-b^2/a^2$ )); but an easier alternative is to use Carlson elliptic integrals, defined as

$$\begin{array}{c}\multicolumn{1}{c}{R_F(p,q,r)=\frac{1}{2}{\int }_0^{\infty }\frac{\mathrm{dt}}{{(t+p)}^{1/2}{(t+q)}^{1/2}{(t+r)}^{1/2}}}\\ \multicolumn{1}{c}{R_D^{•}(p,q;r)\equiv \frac{1}{2}{\int }_0^{\infty }\frac{\mathrm{dt}}{{(t+p)}^{1/2}{(t+q)}^{1/2}{(t+r)}^{3/2}}}\end{array}$$

(Carlson ³ or Greenwood ⁴ (Note that all three parameters in $R_F$ are interchangeable, but only the first two in $R_D^{•}$ . ${\mathrm{The}\mathrm{spot}}^{•}$ indicates that this is one-third of Carlson’s $R_D$ )). Now no mathematical dexterity is required to see that

$$\frac{1}{2R_1}=\frac{\mathrm{ab}{p}_0}{E\prime }{R}_D^{•}(0,b^2;a^2);\frac{1}{2R_2}=\frac{\mathrm{ab}{p}_0}{E\prime }{R}_D^{•}(0,a^2;b^2).$$

(Equally readily, the displacement is found to be $\delta =\frac{\mathrm{ab}{p}_0}{E\prime }{R}_F(0,a^2,b^2)$ .)

Thus, using the natural affinity between Love’s potential function and Carlson elliptic integrals, we have replaced a pair of strongly differing equations by a pair of simple equations. Actually, we have done better than that: we have replaced the pair by asingleequation. For the curvature in the x-direction, the square of the semi-axis in the x-direction is the third argument: for the curvature in the y-direction, the third argument is the square of the semi-axis in the y-direction. We have a single equation, readily obtained, with no mention of major and minor axes.

True, there are no tables of Carlson integrals, and it appears that the two functions K(e) and E(e) of a single variable have been replaced by two functions of two variables. But here the Carlson scaling laws come in. A simple change of integration variable (substitute $s=\mathrm{ct}$ ) shows that $(c^{1/2})R_F(\mathrm{cp},\mathrm{cq},\mathrm{cr})=R_F(p,q,r)$ and $(c^{3/2})R_D^{•}(\mathrm{cp},\mathrm{cq};\mathrm{cr})=R_D^{•}(p,q;r)$ , so that one argument can readily be reduced to 1. Thus the complete Legendre integrals of a single variable can be replaced by Carlson integrals with two passive arguments (0, 1) and calculated as functions of the single remaining variable. The incomplete Legendre integrals, functions of $e$ (more often denoted by $k$ or $m\equiv k^2$ ) and the modular angle could be replaced by Carlson integrals with just two active arguments (i.e. differing from unity). In the computer age, this will probably never be done: a very simple, very rapid algorithm for calculating Carlson integrals is available. It is as simple as the recommended algorithm for the complete Legendre integrals, and since the same algorithm is used for all arguments, is distinctly preferable to the usual method for incomplete Legendre integrals.

A previous article ⁴ shows the advantages to be gained by the use of Carlson integrals in the analysis of the surface properties of Hertzian contacts: the determination of the contact geometry as cited above, the determination of the surface stresses under normal pressures as discussed by Hertz, and those caused by tangential loads as studied by Mindlin and others. Here we consider to what extent the analysis of the subsurface stresses could equally benefit.

The prospective advantages are of course the same: (1) Natural affinity between Love’s potentials and Carlson integrals (2) No need to distinguish between major and minor semi-axes: to assume $a\ge b$ and focus on the ratio $b/a$ : (3) Treating the two directions equally promotes a beneficent mind set: (4) The convenience of ending with a single equation instead of two, usually very different, ones. And to the author, perhaps alone, the advantage of dimensioned equations, providing an immediate check for errors: others may not need this!

The convenience of Carlson integrals over complete Legendre integrals is perhaps minor, the real gain is in the easier derivation: but see Discussion for what happens when we encounter incomplete Legendre integrals.

Analysis of Elliptical Hertzian Contact

The basis of the solution is spelt out by Love, ¹ and taken on from there by Thomas and Hoersch, ⁵ [“T&H”] and it is hard to improve on their presentation, up to the point where they are reduced to converting their answers into Legendre integrals, and so forcing their admirable, symmetric, equations into the inevitable contorted form.

Love introduced two potential functions ϕ and χ, defined as

$$\begin{array}{c}\multicolumn{1}{c}{\phi ={\int }_{\gamma }^{\infty }[1-\frac{x^2}{a^2+t}-\frac{y^2}{b^2+t}-\frac{z^2}{t}]}\\ \multicolumn{1}{c}{·\frac{\mathrm{dt}}{\sqrt{(a^2+t)(b^2+t)(t)}}\mathrm{and}\chi ={\int }_z^{\infty }\phi \mathrm{dz},}\end{array}$$

(1)

where $\gamma$ is the largest root of

$$f\equiv \frac{x^2}{a^2+\gamma }+\frac{y^2}{b^2+\gamma }+\frac{z^2}{\gamma }-1=0$$

(2)

Love gives equations only for the displacements: it seems that T&H were the first to continue and find the strains and hence obtain the basic equations for the stresses (see also Johnson ² or Hills et al. ⁶ (HNS)).

(A factor $(\mathbf{\pi }\mathbf{/}2\mathbf{)}\mathit{ab}{\mathit{p}}_0$ will subsequently be inserted for a Hertzian pressure distribution.)

$$2\pi {\sigma }_x=2\nu \frac{\partial \phi }{\partial z}-z\frac{{\partial }^2\phi }{\partial x^2}-(1-2\nu )\frac{{\partial }^2\chi }{\partial x^2}$$

(3a)

$$2\pi {\sigma }_y=2\nu \frac{\partial \phi }{\partial z}-z\frac{{\partial }^2\phi }{\partial y^2}-(1-2\nu )\frac{{\partial }^2\chi }{\partial y^2}$$

(3b)

$$2\pi {\sigma }_z=\frac{\partial \phi }{\partial z}-z\frac{{\partial }^2\phi }{\partial z^2}$$

(3c)

Subsequently Fessler and Ollerton ⁷

$$2\pi {\tau }_{\mathrm{xz}}=-z\frac{{\partial }^2\phi }{\partial x\partial z}2\pi {\tau }_{\mathrm{yz}}=-z\frac{{\partial }^2\phi }{\partial y\partial z}$$

(3d)

The equation for ${\tau }_{\mathrm{xy}}$ , which we do not need, is given by Johnson ² and Hills et al. ⁶ added

Note the complete symmetry between the pairs $({\sigma }_x,{\tau }_{\mathrm{xz}})$ and $({\sigma }_y,{\tau }_{\mathrm{yz}})$ : this is lost in the final answers when Legendre integrals are used, but is preserved when we use Carlson integrals.

Since ϕ and χ are both harmonic, and $\phi =\frac{\partial \chi }{\partial z}$ so that $\frac{{\partial }^2\chi }{\partial z^2}=\frac{\partial \phi }{\partial z}$ , then from equations (3)

$$\begin{array}{c}\multicolumn{1}{c}{2\pi ({\sigma }_x+{\sigma }_y+{\sigma }_z)=(1+4\nu )\frac{\partial \phi }{\partial z}+(1-2\nu )\frac{\partial \phi }{\partial z}}\\ \multicolumn{1}{c}{=2(1+\nu )\frac{\partial \phi }{\partial z}}\end{array}$$

(4)

We note that $({\sigma }_{\mathrm{xx}}+{\sigma }_{\mathrm{yy}})$ is independent of the second potential function $\chi$ . Indeed, not all the stresses involve χ, or elliptic integrals; for $({\sigma }_z,{\tau }_{\mathrm{xz}},{\tau }_{\mathrm{yz}})$ are all elementary functions, independent of ν.

Extracts from Thomas and Hoersch’s analysis

We introduce a spurious $c$ , which is identically zero, in order to make clearer the permutations to obtain $\frac{\partial }{\partial y},\frac{\partial }{\partial z}$ from $\frac{\partial }{\partial x}$ : thus we write the potential (1) as

$$\begin{array}{c}\multicolumn{1}{c}{\phi ={\int }_{\gamma }^{\infty }[1-\frac{x^2}{a^2+t}-\frac{y^2}{b^2+t}-\frac{z^2}{c^2+t}]·}\\ \multicolumn{1}{c}{\frac{\mathrm{dt}}{\sqrt{(a^2+t)(b^2+t)(c^2+t)}}}\end{array}$$

(1’)

Differentiating $\phi$ is straightforward, but lengthy. The derivatives of both integrand and lower limit are needed, (though not for the first derivative, for there the integrand vanishes) and T&H show that

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial \gamma }{\partial x}=\frac{2x}{a^2+\gamma }·{(\frac{x^2}{{(a^2+\gamma )}^2}+\frac{y^2}{{(b^2+\gamma )}^2}+\frac{z^2}{{(c^2+\gamma )}^2})}^{-1}}\\ \multicolumn{1}{c}{=\frac{\partial f}{\partial x}/H\mathrm{where}}\end{array}$$

(5)

$$\begin{array}{c}\multicolumn{1}{c}{f=\frac{x^2}{a^2+\gamma }+\frac{y^2}{b^2+\gamma }+\frac{z^2}{c^2+\gamma }-1\mathrm{and}}\\ \multicolumn{1}{c}{H\equiv \frac{x^2}{{(a^2+\gamma )}^2}+\frac{y^2}{{(b^2+\gamma )}^2}+\frac{z^2}{{(c^2+\gamma )}^2}.}\end{array}$$

$$\mathrm{Thus},\frac{\partial \phi }{\partial x}={\int }_{\gamma }^{\infty }-\frac{2x}{a^2+t}·\frac{\mathrm{dt}}{\sqrt{(a^2+t)(b^2+t)(c^2+t)}}$$

and permuting,

$$\frac{\partial \phi }{\partial z}={\int }_{\gamma }^{\infty }-\frac{2z}{c^2+t}·\frac{\mathrm{dt}}{\sqrt{(a^2+t)(b^2+t)(c^2+t)}}$$

(6)

This is readily seen to be $\frac{\partial \phi }{\partial z}=-4z{R}_D^{•}(a^2+\gamma ,b^2+\gamma ;c^2+\gamma )$ .

Hence, for our Hertzian contact,

$$\begin{array}{c}\multicolumn{1}{c}{({\sigma }_x+{\sigma }_y+{\sigma }_z)/p_0=-\mathrm{ab}(1+\nu )·z{R}_D^{•}}\\ \multicolumn{1}{c}{(a^2+\gamma ,b^2+\gamma ;\gamma )}\end{array}$$

The second derivatives involve both terms: thus

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2\phi }{\partial x^2}=-2{\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{\sqrt{{(a^2+t)}^3(b^2+t)(c^2+t)}}}\\ \multicolumn{1}{c}{+\frac{1}{\sqrt{{(a^2+\gamma )}^5(b^2+\gamma )(c^2+\gamma )}}·\frac{4x^2}{H}}\end{array}$$

(7)

which, by writing $p\equiv a^2+\gamma$ ; $q\equiv b^2+\gamma$ ; $r\equiv c^2+\gamma (\equiv \gamma )$ , becomes

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2\phi }{\partial x^2}=-2{\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{\sqrt{{(a^2+t)}^3(b^2+t)(c^2+t)}}}\\ \multicolumn{1}{c}{+\frac{1}{\sqrt{{(p)}^5(q)(r)}}·\frac{4x^2}{H}}\end{array}$$

(7a)

Using Carlson integrals this is $\frac{{\partial }^2\phi }{\partial x^2}=-4R_D^{•}(b^2+\gamma ,c^2+\gamma ;a^2+\gamma )+\frac{1}{\sqrt{{(p)}^5(q)(r)}}·\frac{4x^2}{H}$

Permuting the symbols gives

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2\phi }{\partial z^2}=-2{\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{\sqrt{(a^2+t)(b^2+t){(c^2+t)}^3}}}\\ \multicolumn{1}{c}{+\frac{1}{\sqrt{(p)(q){(r)}^5}}·\frac{4z^2}{H}}\end{array}$$

(7b)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{or}\frac{{\partial }^2\phi }{\partial z^2}=-4R_D^{•}(a^2+\gamma ,b^2+\gamma ;c^2+\gamma )}\\ \multicolumn{1}{c}{+\frac{1}{\sqrt{(p)(q){(r)}^5}}·\frac{4z^2}{H}}\end{array}$$

Note that $z\frac{{\partial }^2\phi }{\partial z^2}=\frac{\partial \phi }{\partial z}+\frac{1}{\sqrt{\mathrm{pq}{r}^5}}·\frac{4z^3}{H}$ .

For $\frac{{\partial }^2\phi }{\partial x\partial z}$ the derivative of the integrand vanishes, leaving (using (5)) just

$$\frac{{\partial }^2\phi }{\partial x\partial z}=+\frac{2x}{\sqrt{p^{3}qr}}·\frac{2z}{r}·\frac{1}{H}=\frac{4\mathrm{xz}}{\sqrt{p^{3}q{r}^3}}·\frac{1}{H}$$

(7c)

From the basic equation (3) it will be clear that for an incompressible solid $(\nu =1/2)$ these are all we need: $2\pi {\sigma }_x=2\nu \frac{\partial \phi }{\partial z}-z\frac{{\partial }^2\phi }{\partial x^2}-(1-2\nu )\frac{{\partial }^2\chi }{\partial x^2}$ reduces to $2\pi {\sigma }_x=\frac{\partial \phi }{\partial z}-z\frac{{\partial }^2\phi }{\partial x^2}$ , while for all values of $\nu$ , $2\pi {\sigma }_z=\frac{\partial \phi }{\partial z}-z\frac{{\partial }^2\phi }{\partial z^2}$ (and equally from (3d) we can obtain ${\tau }_{\mathrm{xz}},{\tau }_{\mathrm{yz}}$ ).

Otherwise we need the second potential χ. T&H begin by showing

$$\frac{{\partial }^2\chi }{\partial x^2}=z\frac{{\partial }^2\phi }{\partial x^2}+F_x(z)\mathrm{where}{F}_x(z)\equiv {\int }_z^{\infty }z\frac{{\partial }^3\phi }{\partial x^2\partial z}\mathrm{dz}$$

(8)

and we need to differentiate (7a): $\frac{{\partial }^2\phi }{\partial x^2}=-2{\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{\sqrt{{(a^2+t)}^3(b^2+t)(c^2+t)}}+\frac{1}{\sqrt{{(p)}^5(q)(r)}}·\frac{4x^2}{H}$ .

The differentiation is practicable only when the second term is absent: so for $x=0$ ( only ) we have

$$\frac{{\partial }^3\phi }{\partial z\partial x^2}=\frac{2}{\sqrt{{(a^2+\gamma )}^3(b^2+\gamma )(c^2+\gamma )}}\frac{\partial \gamma }{\partial z}\mathrm{so}\mathrm{that}$$

$$F_x(z)=2{\int }_{{\gamma }_1}^{\infty }\frac{zd\gamma }{\sqrt{{(a^2+\gamma )}^3(b^2+\gamma )(c^2+\gamma )}}$$

(9)

where the $z$ in the numerator is itself $z(\gamma )$ and makes the integration difficult! But on the axis we have simply $\gamma =z^2$ : so substituting (and dropping $c^2$ ) we get:

$$F_x(z)=2{\int }_{{\gamma }_1}^{\infty }\frac{d\gamma }{\sqrt{{(a^2+\gamma )}^3(b^2+\gamma )}}=2I(b;a|{\gamma }_1)$$

and T & H show that the integral $I(b;a|\gamma )$ equals $\frac{2}{a^2-b^2}[1-\sqrt{\frac{b^2+\gamma }{a^2+\gamma }}]$ as is readily checked by differentiation. T&H stop there, with an equation that appears to fail for a circular contact, but a little algebra gives

$${\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{\sqrt{{(a^2+t)}^3(b^2+t)(c^2+t)}}\equiv 2R_D^{•}(q,r;p)$$

(10)

However T&H needed numerical answers for their stresses, and in 1930 the only elliptic integrals available were Legendre’s: so they transformed (spoilt?) their beautiful, symmetric, integrals by substituting $t=a^2{\omega }^2,\gamma =a^2{\varsigma }^2$ and writing $k=b/a$ , so that, for example, the $\frac{{\partial }^2\phi }{\partial x^2}$ integral (7a)

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{\sqrt{{(a^2+t)}^3(b^2+t)(c^2+t)}}\equiv 2R_D^{•}(q,r;p)}\\ \multicolumn{1}{c}{\mathrm{became}\frac{1}{a^3}{\int }_{\varsigma }^{\infty }\frac{\mathrm{dt}}{\sqrt{{(1+{\omega }^2)}^3(k^2+{\omega }^2)}},}\end{array}$$

which they identified as $\frac{1}{a^3}·\frac{1}{{k\prime }^2}[F(\vartheta ,k\prime )-E(\vartheta ,k\prime )]$ (where ${k^{\prime }}^2=1-k^2$ and $\varsigma =\cot \vartheta$ . (Inconsiderate of Love to use for his potential the standard symbol for the modular angle)).

But the $\frac{{\partial }^2\phi }{\partial y^2}$ integral (7b), which we know from symmetry to be ${\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{\sqrt{(a^2+t){(b^2+t)}^3(c^2+t)}}\equiv 2R_D^{•}(p,r;q)$ ,

was found to become $\frac{1}{a^3}\frac{1}{{k\prime }^2}[E(\vartheta k\prime )/k^2-F(\vartheta ,k\prime )]$   +   algebraic terms (see Discussion).

The Carlson forms needed no changes of variable or careful mathematics, and retain the natural symmetry.

Stresses on the axis of symmetry

The stresses are conventionally given as the sum of the answer for $\nu =0$ and a term proportional to $\nu$ , as $\frac{{\sigma }_x}{p_0}={\Omega }_x+\nu \Omega {\prime }_x$ , so we have (after restoring the factor $(\pi /2)\mathrm{ab}{p}_0$ )

$$\begin{array}{c}\multicolumn{1}{c}{{\Omega }_x=-(\mathrm{ab}/2)[z\frac{{\partial }^2\phi }{\partial x^2}+0.5F_x(z)];}\\ \multicolumn{1}{c}{\mathrm{and}\Omega {\prime }_x=(\mathrm{ab}/2)[\frac{\partial \phi }{\partial z}+z\frac{{\partial }^2\phi }{\partial x^2}+F_x(z)].}\end{array}$$

Using Carlson integrals, $\frac{\partial \phi }{\partial z}=-4z{R}_D^{•}(p,q;r)$ and $\frac{{\partial }^2\phi }{\partial x^2}=-4R_D^{•}(q,r;p)$ ,

(on the axis the algebraic term disappears). and also, from (9), (10), $F_x(z)=I(b;a|\gamma )\equiv \frac{2}{\sqrt{p}[\sqrt{p}+\sqrt{q}]}$ .

Then since on the axis $r\equiv \gamma =z^2$ , the stresses are

$${\Omega }_x=2\mathrm{ab}z{R}_D^{•}(q,\gamma ;p)-\frac{\mathrm{ab}}{\sqrt{p}[\sqrt{p}+\sqrt{q})}$$

(11a)

$$\begin{array}{c}\multicolumn{1}{c}{\Omega {\prime }_x=-(2\mathrm{ab}z)\{R_D^{•}(p,q;\gamma )+R_D^{•}(q,\gamma ;p)\}}\\ \multicolumn{1}{c}{+\frac{2\mathrm{ab}}{\sqrt{p}[\sqrt{p}+\sqrt{q})}}\end{array}$$

(11b)

There is a problem here, for at $z=0$ , $R_D^{•}(p,q;z^2)$ is infinite (though the product $z·R_D^{•}(p,q;z^2)$ is finite): but we evade this by using the identity

$$R_D^{•}(p,q;r)+R_D^{•}(q,r;p)+R_D^{•}(r,p;q)=\frac{1}{\sqrt{\mathrm{pqr}}}$$

(12)

(proved in Greenwood ⁴ (Appendix A3); quoted (in a different notation) in Hills et al. ⁶ (10.52)).

Then equation (11b) can be converted into

$$\begin{array}{c}\multicolumn{1}{c}{\Omega {\prime }_x=-(2\mathrm{ab})\{\frac{z}{\sqrt{\mathrm{pq}\gamma }}-z{R}_D^{•}(p,\gamma ;q)\}}\\ \multicolumn{1}{c}{+\frac{2\mathrm{ab}}{\sqrt{p}(\sqrt{p}+\sqrt{q})}}\end{array}$$

or combining the first and last terms (using $\gamma =z^2$ )

$$\Omega {\prime }_x=(2\mathrm{ab}z)R_D^{•}(q,\gamma ;p)\}-\frac{2\mathrm{ab}}{\sqrt{q}[\sqrt{p}+\sqrt{q})}$$

(11c)

To bring out that a simple permutation is all that is needed to give the stress in the b-direction, the two results (11a and 11c) are displayed as Table 1 (We have used the permuted $F_y(z)=\frac{4}{\sqrt{q}[\sqrt{p}+\sqrt{q}]}$ ).

Table 1.

Stresses on the axis of symmetry.

	$[{\Omega }_x//{\Omega }_y]/\mathrm{ab}$	$[\Omega {\prime }_x//\Omega {\prime }_y]/\mathrm{ab}$
“a dirn”	$2z{R}_D^{•}(q,\gamma ;p)-\frac{1}{\sqrt{p}[\sqrt{p}+\sqrt{q})}$	$(2z)R_D^{•}(p,\gamma ;q)\}-\frac{2}{\sqrt{q}[\sqrt{p}+\sqrt{q})}$
“b dirn”	$2z{R}_D^{•}(p,\gamma ;q)-\frac{1}{\sqrt{q}[\sqrt{p}+\sqrt{q})}$	$(2z)R_D^{•}(q,\gamma ;p)\}-\frac{2}{\sqrt{p}[\sqrt{p}+\sqrt{q})}$

${\sigma }_z/p_0=-\mathrm{ab}/\sqrt{\mathrm{pq}}$ ; where $p=a^2+z^2,q=b^2+z^2,\gamma =z^2.$

It must be emphasised that these are not the two distinct equations for the stresses along the directions of the major and minor axes: they are the stresses along the “a” and “b” directions, either of which might be the major semi-axis of the ellipse: there is really only a single equation .

Verification

Thomas and Hoersch, ² with the handicap (in 1930) of working from tables of incomplete Legendre functions at $5^o$ intervals, tabulate the stresses on the axis for $\nu =0.25$ and a number mildly elliptical geometries (appropriate for wheel on rail contacts). Their answers are related to the principal curvatures at the contact, but have been converted to values of $\sigma /p_0$ by taking the value they give for ${\sigma }_z$ at $z=0$ to be $p_0$ . For their moderate case, $k=\sin ({35}^0)=0.57358$ there is almost perfect agreement with the present calculations (occasional $\pm 1$ in the fourth decimal). Figure 1 shows the (less good) comparison for $k=\sin ({20}^0)=0.34202$ , the markers being the T&H scaled values: the lines the present ones. It is clear that the present equations using Carlson integrals are completely equivalent to the T&H equations using Legendre integrals.

$$[b/a=0.342:\nu =0.25].$$

Figure 1.

Comparison of present answers with results tabulated by Thomas and Hoersch.

It is well-known that for a circular contact the value of the critical stress on the axis depends on Poisson’s ratio, and T&H’s equations have always allowed an investigation of how this result transfers to elliptical contacts, but seems not to have been done. It is very readily accomplished using the equations above, and Figure 2 shows the result.

Figure 2.

Maximum von Mises stress. Poisson’s ratio has a significant effect.

Note that we plot the values of the greatest von Mises stress on the axis, defined as $\sqrt{0.5·[{({\sigma }_z-{\sigma }_x)}^2+{({\sigma }_x-{\sigma }_y)}^2+{({\sigma }_y-{\sigma }_z)}^2]}$ . The maximum is for a circular contact, and for $\nu =0.3$ is $0.62p_0$ . Johnson ⁶ quotes a critical stress of $0.31p_0$ , but this is for the maximum principal shear stress, here equal to $({\sigma }_z-{\sigma }_r)/2$ . For a circular contact where ${\sigma }_x={\sigma }_y$ this is indeed just half our von Mises stress.

Progress so far

A single, simple equation has been derived from which both the stresses on the axis ${\sigma }_x$ and ${\sigma }_y$ can be found by permuting the variables. These replace the pair of more complicated equations derived by T&H. T&H regard converting their symmetric integrals into Legendre incomplete elliptic integrals as “the authors’ innovation”: and comment “The direct derivation of these results is rather tedious” (it is not given). Perhaps I should write “the present innovation is converting the integrals into Carlson elliptic integrals” and “The direct derivation of these results is quick and straightforward”!

The full stress distribution on a principal plane

The general equations still hold, so combining (3a) with (8), we get

$$2\pi {\sigma }_x=2\nu \frac{\partial \phi }{\partial z}-2(1-\nu )z\frac{{\partial }^2\phi }{\partial x^2}-(1-2\nu )F_x(z)$$

(12)

where $F_x(z)\equiv {\int }_z^{\infty }z\frac{{\partial }^3\phi }{\partial z\partial x^2}\mathrm{dz}=2{\int }_{{\gamma }_1}^{\infty }\frac{zd\gamma }{\sqrt{{(a^2+\gamma )}^3(b^2+\gamma )(c^2+\gamma )}}$ provided x = 0.

In describing the evaluation of $F_x(z)$ for points on the axis, the actual significance of $F_x(z)$ has been glossed over. It is a contribution to ${\sigma }_x$ from points in the y-plane $(x=0)$ . For points on the x-plane $(y=0)$ →the plane containing the major axis →we need to use $F_y(z)$ , and this will only give us ${\sigma }_y$ , not the in-plane stress ${\sigma }_x$ . However since equation (4) for ${\sigma }_x+{\sigma }_y+{\sigma }_z$ does not involve $\chi$ , (nor does ${\sigma }_z$ ) this is good enough: so oddly, to find the in-plane stresses on a principal plane, we find the out-of-plane stress!

The switch from $F_x(z)$ to $F_y(z)$ simply requires the interchange of $a$ and $b$ , but the evaluation (of either) will now be harder, for again the integrand contains $z$ both explicitly and implicitly as $\gamma (z)$ . But the $\gamma (z)$ relation is no longer the simple $\gamma =z^2$ : instead (for the $y=0$ plane) it is $x^2/(a^2+\gamma )+z^2/\gamma =1$ . However the reverse $z(\gamma )$ relation is usable, and converting the original integration wrt $z$ into integration wrt $\gamma$ becomes possible. The somewhat involved procedure is described in the Appendix: we note here that the forward $\gamma (z)$ relation – finding $\gamma$ as the positive root of a quadratic – is still needed, but now only to determine the lower limit of integration: it does not enter the integrand. Recalling equation (9) for $F_x(z)$ but using its dual to benefit from y = 0 on the x-plane and so obtain $F_y(z)$ , we have

$${\Omega }_y=-(\mathrm{ab}/2)[-4z{R}_D^{•}(p,\gamma ;q)+0.5F_y(z)];$$

(11a)

$$\Omega {\prime }_y=(\mathrm{ab}/2)[-4z{R}_D^{•}(p,q;\gamma )-4z{R}_D^{•}(p,\gamma ;q)+F_y(z)].$$

(11b).

From the Appendix, we have

$$F_y(z)=\frac{2}{(a^2-b^2)}·[2\{\sqrt{\frac{d^2+\gamma }{b^2+\gamma }}-1\}+x^{2}J]$$

(12)

where

$$\mathrm{for}\mathit{a}>\mathit{b},J=\frac{2}{xe}[{\tanh }^{-1}\left\{\frac{e}{x}\right\}-{\tanh }^{-1}\left\{\frac{e}{x}\sqrt{\frac{d^2+\gamma }{b^2+\gamma }}\right\}]$$

(A5b)

$$\mathit{b}>\mathit{a}J=\frac{2}{xe}[{\tan }^{-1}\left\{\frac{e}{x}\right\}-{\tan }^{-1}\left\{\frac{e}{x}\sqrt{\frac{d^2+\gamma }{b^2+\gamma }}\right\}]$$

(A6)

where $r\equiv \sqrt{|a^2-b^2|}$ (but see Appendix for d = b, i.e. for $x^2=a^2-b^2$ )

To avoid the confusion of denoting the in-plane principal stress as ${\sigma }_x$ when on the $y=0$ plane but as ${\sigma }_y$ for the $x=0$ plane, we shall (as for a circular contact) use ${\sigma }_r$ for the radial, in-plane, stress and ${\sigma }_{\theta }$ for the out-of-plane stress. Thus, equations (11) are for ${\Omega }_{\theta },\Omega {\prime }_{\theta }$ and give the hoop stress ${\sigma }_{\theta }$

The sum of the principal stresses is given by equation (4), which with the $(\pi /2)\mathrm{ab}{p}_0$ factor included becomes (using (6)):

$$\begin{array}{c}\multicolumn{1}{c}{({\sigma }_r+{\sigma }_{\theta }+{\sigma }_z)/p_0=\frac{\mathrm{ab}(1+\nu )}{4}\frac{\partial \phi }{\partial z}=}\\ \multicolumn{1}{c}{=-\mathrm{ab}(1+\nu )·z{R}_D^{•}(p,q;\gamma )}\end{array}$$

(Note that on the surface within the contact area, the product $\mathrm{ab}·z{R}_D^{•}(p,q;\gamma )$ must be replaced by its limiting value $4\sqrt{1-x^2/a^2}$ (the mathematical limit, not just copied from Hertz!)).

Then, since for the normal stress ${\sigma }_z$ we have (7c) $\frac{{\sigma }_z}{p_0}=-\frac{ab{z}^3}{\sqrt{pq{\gamma }^5}}\cdot \frac{1}{H}$ we at last obtain the radial stress ${\sigma }_r$ .

Figure 3 shows a typical hoop stress distribution. It is indeed much the same as for a circular hertz contact. In the region shown, the hoop stress is entirely compressive, and almost entirely confined to the region below the contact area. The greatest values are of course on the axis, and the greatest value is (for $\nu =0.3$ ) close to $0.9p_0$ . But as is known for a circular contact, and as found by Thomas and Hoersch (and here: see Figure 1), at greater depths very small (hoop) tensile stresses are found. Note a problem of identification here: on the axis $x=y=0$ , which is the hoop stress?].

Figure 3.

Hoop stress distribution on the principal plane containing the major axis.

Figure 4 shows the matching radial stress distribution. Some differences from the hoop stress distribution are clear. The radial stresses extend substantially further down than the hoop stresses; but this we have already seen from Figure 1 showing the stresses along the axis (for a similar geometry but $\nu =0.25$ ). For $b/a=0.3,\nu =0.3$ , we find that at $z=b$ on the axis, ${\sigma }_r=0.187p_0$ while the hoop stress has fallen to ${\sigma }_{\theta }=0.063p_0$ , so this is quite different from a circular contact, where on the axis, the two must be equal.

Figure 4.

Radial stresses on the principal plane containing the major axis.

In contrast, Figure 5 shows that the surface stresses do behave as for a circular contact, with the hoop stress not falling to zero at the contact edge, while the radial stress has become tensile there; and as for a circular contact becoming equal and opposite outside the contact. Figure 4 brings out that the tensile region is not purely at the surface, extending to perhaps $0.2b$ (though whether this is a sensible measure, here on the a-axis, is not clear for comparison, we need to examine the stresses on the principal plane along the minor axis.

Figure 5.

Hoop stresses in the principal plane of the minor axis.

Not surprisingly, significant hoop stresses exist much further beyond the edge of the contact area, although Figure 7 shows that actually on the surface, the behaviour is the same as for a circular contact. The radial stresses (Figure 6), somewhat to my surprise, show the same tendency more strongly.

Figure 6.

Radial stresses in the plane of the minor axis.

Figure 7.

Surface stresses along major and minor axes for $k=1/3$ .

Figure 7 compares the surface stresses along the major and minor axes. We see that they are what one might expect from a knowledge of the pattern for a circular contact. Here they decay from the values at the origin $-{\sigma }_x/p_0=a/(a+b)+2\nu b/(a+b);$ $-{\sigma }_y/p_0=b/(a+b)+2\nu a/(a+b)$ to the known values at the ends of the axis ² (CM 3.68a), followed by a state of pure shear ( ${\sigma }_r=-{\sigma }_{\theta }$ ). The surface stresses outside the contact area decay, but the variation $(1-2\nu )ab/3r^2$ as for a circular contact is found only when $b\approx a$ or $r$ is large: otherwise we have (for $x>a$ ) $\frac{{\sigma }_{\theta }}{p_0}=\frac{\mathrm{ab}(1-2\nu )}{a^2-b^2}[\{(x/e){\tanh }^{-1}(e/x)\}-1]$ or $\frac{{\sigma }_{\theta }}{p_0}=\frac{\mathrm{ab}(1-2\nu )}{a^2-b^2}[\{(x/e){\tan }^{-1}(e/x)\}-1]$ (where $e=\sqrt{|a^2-b^2|}$ ) depending on whether $b<a$ or $b>a$ .

Note that in order to display the comparison, the “major” and “minor” axis results were obtained by setting $a=1$ and taking “b” successively as 1/3 and 3: the different decay rates of the stresses along the two axes is partly due to this scaling. Note also that the equations (3.68a,b) in Hills et al. ⁶ are incorrect: the arguments of the ${\tanh }^{-1}$ and ${\tan }^{-1}$ terms should have been $\mathrm{ae}/x$ and $b^2/\mathrm{aey}$ , not as printed.

Discussion

Equations for the full stress distribution are not new: indeed Sackfield and Hills, ⁸ “S&H,” give equations for the complete subsurface stress distribution, not only, as here, for the stresses on the principal planes.

But they give no derivation. (The outline in HNS ⁶ (chapter 12) suggests that χ was obtained by field point integration over the ellipse as described by Barber, ⁹ but while this could certainly be used to obtain $\phi$ , it is not clear how it could be used when the kernel $1/R$ is replaced by $\ln (R+z)$ .) S&H give a simpler form for the stresses on the principal plane of the major axis, but do not evaluate their equations, possibly dissuaded by the Legendre incomplete integrals. It is easy to sympathise.

For the S&H solution (and ours) depends on the Thomas & Hoersch integrals $I_1,I_2,I_3$ , declared by S&H to be elliptic integrals, as indeed they are. Thus

$$\begin{array}{c}\multicolumn{1}{c}{I_2={\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{{(a^2+t)}^{1/2}{(b^2+t)}^{3/2}{t}^{1/2}}}\\ \multicolumn{1}{c}{=2R_D^{•}(\gamma ,a^2+\gamma ;b^2+\gamma );}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{I_3={\int }_{\gamma }^{\infty }\frac{\mathrm{dt}}{{(a^2+t)}^{1/2}{(b^2+t)}^{1/2}{t}^{3/2}}}\\ \multicolumn{1}{c}{=2R_D^{•}(a^2+\gamma ,b^2+\gamma ;\gamma )}\end{array}$$

and $I_1$ is just the third permutation of the same three parameters. The Carlson forms above are immediately found (just set $t=w+\gamma$ ) from the integrals in the form in which they first appear in T&H’s analysis.

But $I_1,I_2,I_3$ are not “Legendre elliptic integrals”: while $I_1$ is indeed just the difference of two Legendre integrals, $I_2$ and $I_3$ are each the sum of one or two Legendre integral integrals and algebraic terms. To derive them T&H begin by setting $t=a^2{\varpi }^2$ : then (dropping a factor $2a^{-3}$ ) $I_2={\int }_s^{\infty }\frac{d\varpi }{{(1+{\varpi }^2)}^{1/2}{(k^2+{\varpi }^2)}^{3/2}}$ ; $I_3={\int }_s^{\infty }\frac{d\varpi }{{(1+{\varpi }^2)}^{1/2}{(k^2+{\varpi }^2)}^{1/2}{\varpi }^2}$ (the relation to the Hertz problem is now lost). They then find (laboriously)

$$I_2=\frac{1}{{k^{\prime }}^2}\left[\frac{E(\phi ,k^{\prime })}{k^2}-F(\phi ,k^{\prime })-\frac{1}{k^{2}s}\sqrt{\left(\frac{k^2+s^2}{1+s^2}\right)}+\frac{1}{s}\sqrt{\left(\frac{1+s^2}{k^2+s^2}\right)}\right]$$

$$I_3=\frac{1}{k^{2}s}\sqrt{\frac{k^2+s^2}{1+s^2}}-\frac{E(\phi ,k\prime )}{k^2}$$

two apparently unrelated equations, which cannot immediately be checked.

Unfortunately, the elliptic integrals are only a part of the solution, and are accompanied by a “tail”: our $F_y(z)$ or S&H’s $I_4$ , which are complicated “elementary” quantities. It is not clear how the very different expressions are related, but numerical tests (but using our forms for $I_2,I_3$ ) establish that the answers are the same (A different form of the Sackfield and Hills equations for the full stress distribution is given by Hills et al., ⁶ from which the tail for the stresses on the $y=0$ plane may be abstracted. But our attempt to evaluate this gave impossible answers). So equations for the subsurface stresses have been available since 1983. What is new, is that the derivation of our equations is given; and that for the first time the stresses have been evaluated and plotted, so that the (minor) distortion of the stresses under a circular contact by the ellipticity may be seen and evaluated. The undoubted easing of the calculation of the elliptic integrals is somewhat offset by the existence of the “tails,” although we would claim that our tail is the simpler.

And indeed, with a computer program capable of complex arithmetic, the need to distinguish between $a>b$ and $b>a$ finally disappears: we found (using MATLAB) that the single equation (A5) does produce the same answers as the alternative two real equations (A6) and (A5b).

Conclusion

The simple, direct, link between the equations for the stresses underneath an elliptical Hertzian contact and (readily calculable) Carlson elliptic integrals has been used to obtain equations for the stresses in the principal planes, while avoiding the distinction between major and minor axes. The stresses along the axis of symmetry, well known from the work of Thomas and Hoersch, ⁵ are now given by a single (simpler) equation, in which “a” and “b” appear on an equal basis. Away from the axis of symmetry, the equations for the stresses involve both elliptic integrals and a long “tail” of complicated elementary functions. The considerable gain from using the immediately obvious, neatly matched, Carlson integrals instead of the rather messy transformation into non-matching incomplete Legendre integrals is diluted by the existence of the tails, though we would argue that ours is the simpler. The final numerical answers are the same as those obtained from the equations given in 1983 by Sackfield and Hills, now for the first time evaluated: but our equations are complete with a derivation.

Maps of the principal stresses ${\sigma }_r$ and ${\sigma }_{\theta }$ for a typical elliptical contact ( $k=1/3$ ) are presented. These are much as might be guessed from a knowledge of the stresses under a circular contact, but some differences are noted, particularly when the stresses in the plane of the minor axis are compared to those in the plane of the major axis.

Sadly, the interchangeability of a and b in the elliptic integrals is not maintained in the tails.

The tail for the stresses in the plane of the major principal axis contains a term in arctanh (or a logarithm): the corresponding term in the tail for stresses in the minor principal plane is an arctan. But this was there already in the equations for the surface stresses, and so unavoidable.