Hertzian and adhesive plane models of contact of two inhomogeneous elastic bodies

Abstract

Previous study of contact of power-law graded materials concerned the contact of a rigid body (punch) with an elastic inhomogeneous foundation whose inhomogeneity is characterised by the Young modulus varying with depth as a power function. This paper models Hertzian and adhesive contact of two elastic inhomogeneous power-law graded bodies with different exponents. The problem is governed by an integral equation with two different power kernels. A nonstandard method of Gegenbauer orthogonal polynomials for its solution is proposed. It leads to an infinite system of linear algebraic equations of a special structure. The integral representations of the system coefficients are evaluated, and the properties of the system are studied. It is shown that if the exponents coincide, the infinite system admits a simple exact solution that corresponds to the case when the Young moduli are different but the exponents are the same. Formulas for the length of the contact zone, the pressure distribution and the surface normal displacements of the contacting bodies are obtained in the form convenient for computations. Effects of the mismatch in the Young moduli exponents are studied. A comparative analysis of the Hertzian and adhesive contact models clarifies the effects of the surface energy density on the contact pressure, the contact zone size and the profile of the contacting bodies outside the contact area.

1 Introduction

Interest in contact problems of interaction of bodies with elastic inhomogeneous foundations was originated in the forties of the previous century when civil engineers started taking into account the inhomogeneity properties of soil foundations. For the last thirty years, when novel functionally graded materials (FGMs) were designed and the necessity of the study of their properties arose [28], this interest became even stronger.

One of the most interesting classes of FGMs comprises inhomogeneous materials whose modulus of elasticity E varies with depth according to the power-law, $E(z)=E_\alpha z^\alpha$ . The first approximate solution of a contact problem of an axisymmetric foundation with the modulus of elasticity $E(z)=E_\alpha z^\alpha$ and subjected to a point force P applied to the boundary was obtained [18] in the form

(1.1) \begin{equation}\sigma_z=\frac{A P z^A}{2\pi R^{A+2}}, \quad\sigma_r=\frac{A P z^{A-2}r^2}{2\pi R^{A+2}},\quad \sigma_\theta=0,\quad \tau_{rz}=\frac{A P z^{A-1}r}{2\pi R^{A+2}},\end{equation}

where $R=\sqrt{r^2+z^2}$ . This solution satisfies the equilibrium equations for any values of the constant A. However, in general, the compatibility conditions for the strains are not met. It was found [18] that in only two cases, (1) $A=\alpha+3$ , $\nu=(2+\alpha)^{-1}$ and (2) $A=\alpha+2$ , $\nu=(1+\alpha)^{-1}$ , where $\nu$ is the Poisson ratio, the compatibility conditions are fulfilled. Also, in these particular cases it is possible to recover the normal displacement in the interior of the body by explicitly integrating the strain $e_z=(\sigma_z-\nu\sigma_r)(E_\alpha z^{\alpha})^{-1}$ . Upon passing to the limit $z\to 0$ in the resulting formula for the displacement, this gives the normal displacement on the surface of the half-space, $w(x,y,0)=\alpha P/{\pi E_\alpha} r^{-\alpha-1}$ , where $r=\sqrt{x^2+y^2}.$ Based on the solution obtained in these cases, Klein [18] suggested to extrapolate the formula for the displacement w(x,0), valid in only these two particular cases, to the general case when the Poisson ratio $\nu$ and the exponent $\alpha$ are not connected by any relation. Lekhnitskii [20] considered the plane problem of a wedge with a variable modulus of elasticity. On applying the separation of variables method to the equilibrium equations, he obtained an exact formula for the stress $\sigma_r$ in a half-plane $\{|x|<\infty, y>0\}$ when $E=E_\alpha y^\alpha$ for any constant Poisson ratio $\nu$ . By separating the variables in the equation for the Airy function Rostovtsev [27] not only revisited the Lekhnitskii formula for the stress but also obtained the exact representation for the normal displacement in the cases of concentrated and distributed normal load applied to the boundary.

A majority of results on plane and axisymmetric contact problems of power-law graded materials concern the indentation of a rigid two-dimensional or axisymmetric stamp into a half-plane or a half-space. In the case of a single contact zone, the plane problem reduces to the integral equation

(1.2) \begin{equation}\gamma_0\int_{-b}^b\frac{p(\xi)d\xi}{|x-\xi|^\alpha}=\delta-f(x), \quad -b<x<b,\quad 0<\alpha<1,\end{equation}

where $\delta$ is the indentation of the stamp, the function f(x) describes the stamp profile, and $\gamma_0$ is a function of $\alpha$ . The solution of this equation in the class of functions admitting integrable singularities at the endpoints $\pm b$ exists and unique. It can be constructed by a variety of methods including the method of Abelian integrals (see for example, [8, 1]), the method of dual integral equations, the Wiener-Hopf method and the method of orthogonal polynomials. The solution of this integral equation by the last method is presented in Section 5 of this paper. Popov [24] considered a more advanced case of this plane problem when there are two separate contact zones. He reduced the problem to two separately solvable equations with the Weber-Schafheitlin kernel and solved them approximately by the method of the Jacobi polynomials. The first exact solutions [19, 21] to the axisymmetric case were obtained by the method of dual integral equations under the assumption of the frictionless contact of a stamp and a power-law graded foundation. The same problem was later solved [22] by the Wiener-Hopf method. The method of Abelian operators was applied [25] to derive an exact solution to the axisymmetric problem of non-slipping adhesive contact of a punch with a power-law graded elastic half-space.

During the last twenty-five years, plane and axisymmetric contact problems of a stamp and a half-plane and a half-space with the Young modulus $E=E_\alpha z^\alpha$ have become the subject of interest due to modelling of micro- and nano-indentation processes arising in nanotechnology and therefore the necessity of characterisation of mechanical properties of a variety of biological materials with sizes approaching molecular or atomic dimension [12, 3]. In [9], the Rostovtsev formulas [27] were revisited and normal, sliding and rolling contact models were addressed. 2d and axisymmetric models of adhesive contact of a rigid stamp and a power-law graded semi-infinite body characterised by the Young modulus $E=E_0(y/c_0)^k$ were considered in [6] and [7]. It was found that the pull-off force attains its maximum at a certain critical value of $k\in(0,1)$ . The method of orthogonal polynomials [25] was employed [12] to study a model of non-slipping adhesive contact of a stamp and a power-law graded elastic half-space. Another axisymmetric adhesive contact model [15] assumes that the contact circular area is split into a non-slipping circular zone and an annulus where the tangential traction equals zero, while the normal traction is tensile and equals a prescribed constant. The penetration depth and contact radius in the case of axisymmetric contact of a stamp and a power-law graded half-space were calculated [14] by exploiting a correspondence between these quantities and the ones associated with a Winkler foundation. All papers discussed in this paragraph adopted the Johnson-Kendall-Roberts (JKR) adhesive model ([16, 17]) to examine plane and axisymmetric contact of a rigid punch with a half-plane and half-space, respectively, when the Young modulus of the foundation varies with depth according to a power-law. The feature of the JKR model is that it admits integrable singularities of the contact pressure at the endpoints and determines the contact zone radius from the condition of minimum of the total energy. The total energy $U_{total}$ is defined to be a sum of the elastic strain energy $U_e$ and the loss of surface energy $U_s$ .

There have been relatively limited efforts in studying Hertzian and adhesive contact of two elastic bodies whose Young moduli are power-law functions of depth. The axisymmetric Hertzian model of contact of two bodies whose Young moduli are expressed through the same power-law function, $E_1(z)=e_1 z^\alpha$ and $E_2(z)=e_2 (-z)^\alpha$ , was considered in [26]. In this work, the surface effects are taken into account according to the Shtayerman model [29]. Power-law kernels arise [13] in the problem of computing equilibrium measures for problems with attractive-repulsive kernels of the form $K(x-y)=\alpha^{-1}|x-y|^\alpha-\beta^{-1}|x-y|^\beta$ . For this problem, they proposed a numerical method of recursively generated banded and approximately banded operators acting on expansions in ultraspherical polynomial bases. To the authors’ knowledge, neither the two-dimensional nor the axisymmetric problem of Hertzian or JKR adhesive contact of two elastic bodies with different Young moduli, $E_1(z)=e_1 z^{\alpha_1}$ and $E_2(z)=e_2 (-z)^{\alpha_2}$ , has been considered in the literature.

In this paper, we aim to analyse the plane contact problem of two different power-law graded bodies. In Section 2, we formulate the problem and reduce it to the integral equation with two kernels of the form

(1.3) \begin{equation}\int_{-b}^{b}\left(\frac{\gamma_1}{|x-\xi|^{\alpha_1}}+\frac{\gamma_2}{|x-\xi|^{\alpha_2}}\right)p(\xi)d\xi=g(x), \quad -b<x<b.\end{equation}

This equation may be interpreted as a full integral equation with a single power kernel $|x-\xi|^{-\alpha_1}$ with the second kernel serving as a regular part [8]. However, the method of Abelian operators, when applied, leads to a Fredholm integral equation whose kernel is a chain of singular integrals and does not produce the solution in the form convenient for numerical purposes.

In Section 3, we describe the method of solution that expands the unknown function p(bt) in terms of the Gegenbauer polynomials $C_n^{\alpha_1/2}(t)$ with weight $(1-t^2)^{(\alpha_1-1)/2}(t)$ and reduces the task of finding the expansion coefficients to solution of an infinite system of linear algebraic coefficients with coefficients represented by integrals possessing the polynomials $C_n^{\alpha_1/2}(t)$ and $C_m^{\alpha_2/2}(t)$ . We manage to evaluate these integrals. The coefficients have certain remarkable properties which substantially simplify the system. We also show that in the limit case $\alpha_2\to\alpha_1$ , the solution of the infinite system can be derived explicitly, and it coincides with the solution of the contact problem of two bodies with different power-law Young moduli and the same exponent, $E_1=e_1y^{\alpha}$ and $E_2=e_2 (-y)^{\alpha}$ .

In Section 4, we derive formulas for the length of the contact zone, the parameter $\delta$ , the pressure distribution and the normal displacement on the surface outside the contact zone in the form convenient for computations. We emphasise that all the formulas except for the displacement are free of integrals. We also discuss the results of numerical tests.

In Section 5, we derive a closed-form solution of the problem of Hertzian contact of two bodies whose moduli of elasticity have the same exponents $\alpha_1=\alpha_2=\alpha$ but different factors $e_1$ and $e_2$ . We obtain exact formulas not only for the contact zone length and the pressure but also for the normal displacement outside the contact area.

In Section 6, we analyse the JKR model for both cases, when $\alpha_1=\alpha_2$ and $\alpha_1>\alpha_2$ . In both cases, we compute the elastic strain energy and the total energy. In the former case, we obtain a transcendental equation for the contact zone half-length b and show that it is possible to pass to the limit as $\alpha_j\to 0$ . In the case, $\alpha_1>\alpha_2$ we derive the equation for b approximately by computing the derivative of the strain energy numerically. We show that in both cases, the solution to the JKR model coincides with the solution to the Hertzian model when the surface energy half-density $\gamma_s\to 0$ .

In Appendix A, we discuss the limiting case $\alpha_j\to 0$ , $j=1,2$ . We compute an integral required for the displacements of surface points outside the contact zone in Appendix B. In Appendix C, we derive an exact solution of the integral equation (1.3) in a semi-infinite interval. In the case $\alpha_2=-2n$ , $n=0,1,\ldots$ , and $\alpha_1\in(0,1)$ , the integral equation (1.3) admits an exact solution. This is addressed in Appendix D.

2 Formulation

The problem of interest is the one of modelling of two-dimensional contact of two inhomogeneous elastic bodies, $B_1$ and $B_2$ (Figure 1(a)). The lower surface of the upper body $B_1$ and the upper surface of the lower body $B_2$ are described by curves $y=f_1(x)$ and $y=-f_2(x)$ . The functions $f_1(x)$ and $f_2(x)$ are even, continuously differentiable and share the tangent line $y=0$ at the point $x=0, y=0$ , the origin of the Cartesian coordinates (x,y), that is $f_1(0)=f_2(0)=0$ and $f^{\prime}_1(0)=f^{\prime}_2(0)=0$ . The Poisson ratios $\nu_1$ and $\nu_2$ of the bodies are constant, while the Young moduli vary according to a power law and equal $E_1(y)=e_1y^{\alpha_1}$ and $E_2(y)=e_2(-y)^{\alpha_2}$ , where $e_1$ and $e_2$ are positive parameters whose dimensions in the SI system are $N\, m^{-2-\alpha_1} $ and $N\, m^{-2-\alpha_2}$ , respectively, and $0<\alpha_j<1$ , $j=1,2$ . The bodies are subjected to compression by forces applied to the bodies parallel to the y-axis with the resultant force P balanced by the contact pressure p(x) arising in the contact area $(-b,b)$ , and the parameter b is unknown a priori. We also assume that the curve $y=f_1(x)$ is convex upward, while the second curve $y=-f_2(x)$ is either convex downward or flat or at least locally convex upward. To proceed with the contact modelling, we make the following Hertzian assumptions:

• the contact area is significantly less than the bodies’ sizes,

• the friction is absent, and the only nonzero traction component is $\sigma_y=-p(x)$ , where p(x) is the normal pressure distribution,

• the normal and tangential elastic displacements in the contact area are significantly smaller than the contact zone length.

Figure 1. Contact of two different power-law graded elastic bodies with Young moduli $E_1(y)=e_1y^{\alpha_1}$ (body $B_1$ ) and $E_2(y)=e_2(-y)^{\alpha_2}$ (body $B_2$ ). (a): Hertzian model and (b): adhesive JKR model.

Following Shtayerman (1949), we write the vertical displacements of any two points $A_1\in B_1$ and $A_2\in B_2$ which, as a result of compression, become the same point, a point A. These displacements are $f_1(x-u_1)+v_1-\delta_1$ and $-f_2(x+u_2)-v_2+\delta_2$ . Here, $(u_1,v_1)$ and $(-u_2,-v_2)$ are the elastic displacements of the points $A_1$ and $A_2$ , and the constants $\delta_1$ and $\delta_2$ are forward displacements of distant points. Approximating $f_1(x-u_1)\approx f_1(x)$ and $f_2(x+u_2)\approx f_2(x)$ , we can write at the point of contact A

(2.1) \begin{equation}v_1+v_2=\delta-f_1(x)-f_2(x), \quad -b<x < b, \quad \delta=\delta_1+\delta_2.\end{equation}

The parameter $\delta$ is to be determined a posteriori from the condition

(2.2) \begin{equation}\int_{-b}^b p(x)dx=P.\end{equation}

We next use the Rostovtsev relation (Rostovtsev, 1964, p. 747) between the normal displacement and the contact pressure for a half-plane to write down the displacements $v_1$ and $v_2$ in the contact area

(2.3) \begin{equation}v_j(x)=\frac{\theta_j}{\alpha_j}\int_{-b}^b\frac{p(\xi)d\xi}{|x-\xi|^{\alpha_j}}, \quad -b<x<b, \quad j=1,2.\end{equation}

Here,

(2.4) \begin{align} \theta_j&= \frac{C_j(1-\nu_j^2)q_j}{(\alpha_j+1)e_j}\sin\frac{\pi q_j}{2},\quad q_j=\sqrt{(1+\alpha_j)\left(1-\frac{\alpha_j\nu_j}{1-\nu_j}\right)}, \nonumber\\[3pt] C_j&=\frac{2^{\alpha_j+1}}{\pi\Gamma(\alpha_j+2)}\Gamma\left(\frac{\alpha_j}{2}-\frac{q_j}{2}+\frac32\right)\Gamma\left(\frac{\alpha_j}{2}+\frac{q_j}{2}+\frac32\right). \end{align}

Substituting the integral representations of the displacements $v_j$ into the condition (2.1), we derive the governing integral equation for the contact pressure distribution p(x)

(2.5) \begin{equation}\int_{-b}^b\left(\frac{\theta_1}{\alpha_1|x-\xi|^{\alpha_1}}+\frac{\theta_2}{\alpha_2|x-\xi|^{\alpha_2}}\right)p(\xi)d\xi=\delta-f_1(x)-f_2(x), \quad -b<x<b.\end{equation}

3 Solution of the integral equation

To solve the integral equation (2.5), it will be convenient to rewrite it in the interval $(-1,1)$

(3.1) \begin{equation}\int_{-1}^1\left(\frac{A_1}{|t-\tau|^{\alpha_1}}+\frac{A_2}{|t-\tau|^{\alpha_2}}\right)p(b\tau)d\tau=\delta-f(bt), \quad -1<t<1,\end{equation}

where

(3.2) \begin{equation}A_j=\frac{\theta_jb^{1-\alpha_j}}{\alpha_j}.\end{equation}

The right-hand side of equation (3.1) possesses the unknown parameter $\delta$ . To eliminate it from the equation, we represent the function p(bt) as

(3.3) \begin{equation}p(bt)=\phi^{(1)}(t)+\delta\phi^{(2)}(t)\end{equation}

and deduce

(3.4) \begin{equation}\int_{-1}^1\left(\frac{A_1}{|t-\tau|^{\alpha_1}}+\frac{A_2}{|t-\tau|^{\alpha_2}}\right)\phi^{(j)}(\tau)d\tau=g^{(j)}(t), \quad -1<t<1, \quad j=1,2,\end{equation}

where $g^{(1)}(t)=-f(bt)$ , $g^{(2)}(t)=1$ . The equilibrium condition (2.2) expresses the unknown parameter $\delta$ through the solutions $\phi_1$ and $\phi_2$ of the equations (3.4) which share the kernel and have different right-hand sides. We have

(3.5) \begin{equation}\delta=\left(\frac{P}{b}-\int_{-1}^1\phi^{(1)}(\tau)d\tau\right)\left(\int_{-1}^1\phi^{(2)}(\tau)d\tau\right)^{-1}.\end{equation}

3.1 Infinite system of algebraic equations

Without loss of generality, we assume further that $\alpha_1>\alpha_2$ and denote

(3.6) \begin{equation}\beta_{n}(\alpha)=\frac{\pi(\alpha)_n}{n!\cos\frac{\pi\alpha}{2}}, \quad n=0,1,\ldots,\end{equation}

where $(\alpha)_n=\alpha(\alpha+1)\ldots(\alpha+n-1)$ is the factorial symbol. Owing to the spectral relation for the Gegenbauer polynomials [23]

(3.7) \begin{equation}\int_{-1}^1 \frac{C_n^{\alpha/2}(\tau)d\tau}{|t-\tau|^{\alpha}(1-\tau^2)^{(1-\alpha)/2}}=\beta_{n}(\alpha)C_n^{\alpha/2}(t), \quad -1<t<1,\quad 0<\alpha<1,\end{equation}

and the orthogonality property of these polynomials

(3.8) \begin{equation}\int_{-1}^1 C_n^{\alpha/2}(t) C_m^{\alpha/2}(t)(1-t^2)^{(\alpha-1)/2}dt=h_n(\alpha) \delta_{mn}, \quad m,n=0,1,\ldots,\end{equation}

we seek the solution in the series form

(3.9) \begin{equation}\phi^{(j)}(t)=(1-t^2)^{(\alpha_1-1)/2}\sum_{n=0}^\infty\Phi_{n}^{(j)} C_n^{\alpha_1/2}(t), \quad -1<t<1, \quad j=1,2.\end{equation}

Here, $\Phi_n^{(j)}$ are unknown coefficients, $\delta_{mn}$ is the Kronecker symbol, and

(3.10) \begin{equation}h_n(\alpha)=\frac{\pi 2^{1-\alpha}\Gamma(n+\alpha)}{n!(n+\frac{\alpha}{2})\Gamma^2(\frac{\alpha}{2})}.\end{equation}

In the integral equations of a rigid stamp indented into an inhomogeneous power-law graded half-plane or Hertzian contact of two bodies with $\alpha_1=\alpha_2$ , there is only one power-law kernel. In these particular cases, the series coefficients can be derived explicitly by substituting the expansion (3.9) into the integral equation and taking into account the spectral relation (3.7) and the orthogonality property (3.8). In contrast to this, when $\alpha_1\ne \alpha_2$ , we have the second term in the kernel, and, in general, the series coefficients cannot be found exactly. On substituting (3.9) into (3.4), we have

(3.11) \begin{equation}A_1\sum_{n=0}^\infty\beta_n(\alpha_1)\Phi_{n}^{(j)}C_n^{\alpha_1/2}(t)+A_2\sum_{n=0}^\infty\Phi_{n}^{(j)}\int_{-1}^1\frac{G_n(\tau)(1-t^2)^{(\alpha_2-1)/2}d\tau}{|t-\tau|^{\alpha_2}}=g^{(j)}(t), \quad -1<t<1,\end{equation}

where $G_n(\tau)=C_n^{\alpha_1/2}(\tau)(1-t^2)^{(\alpha_1-\alpha_2)/2}$ . Since $\alpha_1>\alpha_2$ , we may expand the functions $G_n(\tau)$ in terms of the Genebauer polynomials $C_m^{\alpha_2/2}(\tau)$

(3.12) \begin{equation}G_n(\tau)=\sum_{m=0}^\infty G_m^{(n)} C_m^{\alpha_2/2}(\tau), \quad -1<\tau<1.\end{equation}

According to the orthogonality relation (3.8), the coefficients of the expansion are found to be

(3.13) \begin{equation}G_m^{(n)}=\frac{H_m^{(n)}}{h_m(\alpha_2)}, \quad H_m^{(n)}=\int_{-1}^1C_n^{\alpha_1/2}(\tau)C_m^{\alpha_2/2}(\tau)(1-\tau^2)^{(\alpha_1-1)/2}d\tau.\end{equation}

Notice that $H_m^{(n)}=0$ if $m<n$ . Indeed, the degree-m polynomial $C_m^{\alpha_2/2}(\tau)$ is a linear combination of the monomials $1,\tau, \ldots, \tau^m$ or, equivalently, a linear combination of the Gegenbauer polynomials $C_0^{\alpha_1/2)}(\tau), C_1^{\alpha_1/2}(\tau),\ldots, C_m^{\alpha_1/2)}(\tau)$ , and by the orthogonality relation (3.8) $H_m^{(n)}=0$ provided $m<n$ . Now, if we substitute the series (3.12) back to equation (3.11), use the spectral relation (3.8) for the Gegenbauer polynomials $C_m^{\alpha_2/2}(\tau)$ and change the order of summation, we find

(3.14) \begin{equation}A_1\sum_{n=0}^\infty\beta_n(\alpha_1)\Phi_n^{(j)}C_n^{\alpha_1/2}(t)+A_2\sum_{n=0}^\infty\Psi_{n}^{(j)}\beta_n(\alpha_2)C_n^{\alpha_2/2}(t)=g^{(j)}(t), \quad -1<t<1,\end{equation}

where

(3.15) \begin{equation}\Psi_{n}^{(j)}=\sum_{m=0}^n G_n^{(m)}\Phi_{m}^{(j)}.\end{equation}

The equation (3.14) can be recast by using the orthogonality relation (3.8) and written as an infinite system of algebraic equations. We have

(3.16) \begin{equation}A_1\beta_n(\alpha_1)h_n(\alpha_1)\Phi_{n}^{(j)}+A_2\sum_{m=0}^\infty\Psi_{m}^{(j)}\beta_m(\alpha_2)H_m^{(n)}=g^{(j)}_n, \quad n=0,1,\ldots,\end{equation}

where

(3.17) \begin{equation}g_n^{(j)}=\int_{-1}^1 g_j(t)C_n^{\alpha_1/2}(t)(1-t^2)^{(\alpha_1-1)/2}dt.\end{equation}

It is possible to simplify the system derived. On changing the order of summation in the series in the system (3.16) and using the relations (3.13) and (3.15), we obtain

(3.18) \begin{equation}\sum_{m=0}^\infty \Psi_{m}^{(j)}\beta_m(\alpha_2)H_m^{(n)}=\sum_{m=0}^\infty L_{nm}\Phi_{m}^{(j)},\end{equation}

where

(3.19) \begin{equation}L_{nm}=\sum_{k=\max(m,n)}^\infty \frac{H_k^{(n)}H_k^{(m)}\beta_k(\alpha_2)}{h_k(\alpha_2)}, \quad m,n=0,1,\ldots,\end{equation}

and therefore, the system has the form

(3.20) \begin{equation}A_1\beta_n(\alpha_1)h_n(\alpha_1)\Phi_{n}^{(j)}+A_2\sum_{m=0}^\infty L_{nm}\Phi_{m}^{(j)}=g_{n}^{(j)}, \quad n=0,1,\ldots.\end{equation}

3.2 Evaluation of the integrals $\boldsymbol{H}_{\boldsymbol{k}}^{(\boldsymbol{m})}$

We remind that $H_k^{(m)}=0$ if $k<m$ . To compute the integrals (3.13) when $k\ge m$ , we use the formula

(3.21) \begin{align} &\int_{-1}^1(1-x)^\alpha(1+x)^{\nu-1/2}C_m^\mu(x)C_n^\nu(x)dx =\frac{2^{\alpha+\nu+\frac12}\Gamma(\alpha+1)\Gamma(\nu+\frac12)\Gamma(\nu-\alpha+n-\frac12)}{m! n! \Gamma(\nu-\alpha-\frac12)\Gamma(\nu+\alpha+n+\frac32)} \nonumber\\[3pt] &\quad \times\frac{\Gamma(m+2\mu)\Gamma(n+2\nu)}{\Gamma(2\mu)\Gamma(2\nu)}{}_4 F_3\left(\begin{array}{cccc}-m, & m+2\mu,& \alpha+1,&\alpha-\nu+\frac32\\\mu+\frac12, &\nu+\alpha+n+\frac32,&\alpha-\nu-n+\frac32;&1\\\end{array}\right),\end{align}

where $\mathop{\rm Re}\nolimits\alpha>-1$ , $\mathop{\rm Re}\nolimits\nu>-\frac12$ and ${}_4 F_3$ is the generalised hypergeometric function. This relation can be derived from the general formula 16.4(20) [4] for the Jacobi polynomials. Notice that the corresponding formulas for the Gegenbauer polynomials 16.3(16) [4] and 7.314(7) [11] have the same error: instead of $\Gamma(\nu+\alpha+n+\frac32)$ in the right-hand side in (3.21), they write $\Gamma(\nu-\alpha+n+\frac32)$ .

On adjusting the relation (3.21) to our case when $k\ge m,$ we have

(3.22) \begin{equation}H_k^{\left(m\right)}=\dfrac{2\sqrt{\pi}\left(-1\right)^m \Gamma\left(\dfrac{\alpha_1+1}{2}\right)\left(\alpha_2\right)_k}{m!\Gamma\left(\dfrac{\alpha_1}{2}\right)\left(m+\alpha_1\right)}\Sigma,\end{equation}

where

(3.23) \begin{equation}\Sigma=\sum_{l=m}^k\dfrac{\left(-1\right)^l\left(\alpha_2+k\right)_l \left(\dfrac{\alpha_1+1}{2}\right)_l}{\left(k-l\right)!\left(l-m\right)!\left(\alpha_1+m+1\right)_l\left(\dfrac{\alpha_2+1}{2}\right)_l}.\end{equation}

This sum can be evaluated and the formula for $H_k^{\left(m\right)}$ simplified. We make the substitution $l-m=i$ , use the property of the factorial symbols

(3.24) \begin{equation}(a)_{m+i}=(a+m)_i(a)_m, \quad (k-n)!=\frac{(-1)^n k!}{(-k)_n},\quad k\ge n,\end{equation}

and express the sum $\Sigma$ through the function ${}_3 F_2$

(3.25) \begin{equation}\Sigma=\frac{(-1)^m(\frac{\alpha_1+1}{2})_m(\alpha_2+k)_m}{(k-m)!(\alpha_1+m+1)_m(\frac{\alpha_2+1}{2})_m}{}_3 F_2\left(\begin{array}{ccc}-k+m, & \alpha_2+k+m, & \frac{\alpha_1+1}{2}+m\\\frac{\alpha_2+1}{2}+m, & \alpha_1+2m+1; & 1\\\end{array}\right).\end{equation}

For the generalised hypergeometric function ${}_3 F_2$ in the right-hand side, we can employ Whipple’s formula [31]

(3.26) \begin{equation}{}_3 F_2\left(\begin{array}{ccc}a, & b, & c\\ \dfrac{a+b+1}{2}, & 2c; & 1\end{array}\right)=\dfrac{\sqrt{\pi}\Gamma\left(c+\dfrac12\right)\Gamma\left(\dfrac{a+b+1}{2}\right)\Gamma\left(\dfrac{1-a-b}{2}+c\right)}{\Gamma\left(\dfrac{a+1}{2}\right)\Gamma\left(\dfrac{b+1}{2}\right)\Gamma\left(\dfrac{1-a}{2}+c\right)\Gamma\left(\dfrac{1-b}{2}+c\right)}\end{equation}

and obtain the following representation for $\Sigma$ :

(3.27) \begin{align} \Sigma&=\dfrac{\left(-1\right)^m\left(\dfrac{\alpha_1+1}{2}\right)_m\left(\alpha_2+k\right)_m}{\left(k-m\right)!\left(\alpha_1+m+1\right)_m\left(\dfrac{\alpha_2+1}{2}\right)_m} \nonumber\\[3pt] &\quad \times\dfrac{\sqrt{\pi}\Gamma\left(\dfrac{\alpha_1}{2}+m+1\right)\Gamma\left(\dfrac{\alpha_2+1}{2}+m\right)\Gamma\left(\dfrac{\alpha_1-\alpha_2}{2}+1\right)}{\Gamma\left(\dfrac{m-k+1}{2}\right)\Gamma\left(\dfrac{\alpha_1+m+k}{2}+1\right)\Gamma\left(\dfrac{\alpha_2+m+k+1}{2}\right)\Gamma\left(\dfrac{\alpha_1-\alpha_2+m-k}{2}+1\right)}.\end{align}

This formula implies that $\Sigma=0$ and therefore $H_k^{(m)}=0$ if $k=m+1+2l$ , $l=0,1,\ldots$ . In the case when $k-m$ is even, $k=m+2l$ , $l=0,1,\ldots$ , we substitute (3.27) into (3.22) and find

(3.28) \begin{align} H_{m+2l}^{\left(m\right)}&=\dfrac{2\sin\dfrac{\pi\left(\alpha_2-\alpha_1\right)}{2}\Gamma\left(\dfrac{\alpha_1+1}{2}\right)\left(\dfrac{\alpha_1+1}{2}\right)_m\left(\dfrac{\alpha_1}{2}\right)_{m+1}\Gamma\left(\dfrac{\alpha_1-\alpha_2}{2}+1\right)\Gamma\left(\dfrac{\alpha_2+1}{2}\right)}{\pi m!\Gamma\left(\alpha_2\right)\left(m+\alpha_1\right)\left(\alpha_1+m+1\right)_m} \nonumber\\[3pt] &\times\dfrac{\Gamma\left(2l+2m+\alpha_2\right)\Gamma\left(l+\dfrac12\right)\Gamma\left(\dfrac{\alpha_2-\alpha_1}{2}+l\right)}{\left(2l\right)!\Gamma\left(\dfrac{\alpha_1}{2}+m+l+1\right)\Gamma\left(\dfrac{\alpha_2+1}{2}+m+l\right)}.\end{align}

On exploiting further the properties of the $\Gamma$ -function, it is possible to recast formula (3.28) as

(3.29) \begin{equation}H_{m+2l}^{\left(m\right)}=\dfrac{\sqrt{\pi}\Gamma\left(\dfrac{\alpha_1+1}{2}\right)}{\Gamma\left(\dfrac{\alpha_1}{2}+1\right)}\dfrac{\left(\alpha_1\right)_m\left(\alpha_2/2\right)_m}{m!\left(\dfrac{\alpha_1}{2}+1\right)_m}\dfrac{\left(\dfrac{\alpha_2-\alpha_1}{2}\right)_l\left(\dfrac{\alpha_2}{2}+m\right)_l}{l!\left(\dfrac{\alpha_1}{2}+m+1\right)_l}.\end{equation}

This formula is simpler and convenient for analysis of the coefficients asymptotics as $l\to\infty$ . Taking into account the asymptotic relation

(3.30) \begin{equation}\frac{\Gamma(z+a)}{\Gamma(z+b)}\sim z^{a-b}, \quad z\to \infty, \quad |\arg z|<\frac{\pi}{2},\end{equation}

we derive

(3.31) \begin{equation}H_{m+2l}^{(m)}\sim C_m l^{\alpha_2-\alpha_1-2}, \quad l\to \infty,\end{equation}

where $C_m$ are constants independent of l.

Having computed the coefficients $H_k^{(m)}$ we consider now two cases, $n=0,1,\ldots,m-1$ and $n=m,m+1,\ldots$ and evaluate the coefficients $H_k^{(n)}$ . In the former case according to formula (3.19) and since $H_k^{(m)}=0$ if $k=m+2l+1$ , $l=0,1,\ldots$ , we need to evaluate $H_k^{(n)}$ for $k=m+2l$ only. On replacing m by n and k by $m+2l$ in (3.22) and (3.27), we should have $H_{m+2l}^{(n)}=0$ , if $n-m$ is odd and $l=0,1,\ldots.$ Otherwise, if $n-m$ is even,

(3.32) \begin{align}H_{m+2l}^{\left(n\right)}&=\dfrac{2^{\alpha_2-\alpha_1}\sqrt{\pi}\Gamma\!\left(\dfrac{\alpha_2+1}{2}\right)\Gamma\!\left(\alpha_1+n\right)}{\Gamma\!\left(\dfrac{\alpha_1}{2}\right)\Gamma\!\left(\alpha_2\right)\Gamma\!\left(\dfrac{\alpha_2-\alpha_1}{2}\right)n!} \nonumber\\[3pt] &\quad \dfrac{\Gamma\!\left(\dfrac{\alpha_2+m+n}{2}+l\right)\Gamma\!\left(\dfrac{m-n+\alpha_2-\alpha_1}{2}+l\right)}{\Gamma\!\left(\dfrac{m-n}{2}+l+1\right)\Gamma\!\left(\dfrac{\alpha_1+m+n}{2}+l+1\right)}, \quad l=0,1,\ldots,\end{align}

and their asymptotics for large l is the same as for $H_{m+2l}^{(m)}$ . We have

(3.33) \begin{equation}H_{m+2l}^{(n)}\sim C^{\prime}_{mn} l^{\alpha_2-\alpha_1-2}, \quad l\to \infty,\end{equation}

where $C_{mn}$ are constants. We also give another, more convenient for numerical purposes, representation of the coefficients $H_{m+2l}^{(n)}$ when $n-m$ is even

(3.34) \begin{align}H_{m+2l}^{\left(n\right)}&=\dfrac{\sqrt{\pi}\Gamma\left(\dfrac{\alpha_1+1}{2}\right)}{\Gamma\left(\dfrac{\alpha_1}{2}+1\right)}\dfrac{\left(\alpha_1\right)_n}{n!}\dfrac{\left(\dfrac{\alpha_2}{2}\right)_{\left(m+n\right)/2}}{\left(\dfrac{\alpha_1}{2}+1\right)_{\left(m+n\right)/2}}\dfrac{\left(\dfrac{\alpha_2-\alpha_1}{2}\right)_{\left(m-n\right)/2}}{\left(\dfrac{m-n}{2}\right)!}\nonumber\\[3pt] &\quad \times \dfrac{\left(\dfrac{\alpha_2+m+n}{2}\right)_{\!l} \left(\dfrac{m-n+\alpha_2-\alpha_1}{2}\right)_{\!l}}{\left(\dfrac{m-n}{2}+1\right)_{\!l}\left(\dfrac{\alpha_1+m+n}{2}+1\right)_{\!l}}.\end{align}

3.3 Solution of the infinite system

By introducing new notations, we rewrite the system (3.20) in the canonical form

(3.35) \begin{equation}\Phi_{n}^{(j)}+\gamma\sum_{m=0}^\infty R_{nm}\Phi_{m}^{(j)}=d_{n}^{(j)}, \quad n=0,1,\ldots, \quad j=1,2,\end{equation}

where

(3.36) \begin{equation}\gamma=\frac{A_2}{A_1}, \quad R_{nm}=\frac{L_{nm}}{\beta_n(\alpha_1)h_n(\alpha_1)}, \quad d_{n}^{(j)}=\frac{g_{n}^{(j)}}{A_1\beta_n(\alpha_1)h_n(\alpha_1)}.\end{equation}

Owing to the fact that $H_{m+2l}^{(n)}=0$ , if $n-m$ is odd and $l=0,1,\ldots.$ , from formula (3.19) we deduce that $L_{nm}=0$ and therefore $R_{nm}=0$ if $n-m$ is odd. We have also derived that $H_{k}^{(m)}=0$ if $k=m+2l+1$ and $l=0,1,\ldots$ . This brings us to the following formulas for the coefficients $L_{nm}$ when $m-n$ is even:

(3.37) \begin{align} L_{nm}=\sum_{l=0}^\infty H_{m+2l}^{(m)}H_{m+2l}^{(n)}\Delta_{m+2l}, \quad n=0,1,\ldots,m-1, \nonumber \\[2pt] L_{nm}=\sum_{l=0}^\infty H_{n+2l}^{(n)}H_{n+2l}^{(m)}\Delta_{n+2l}, \quad n=m,m+1,\ldots,\end{align}

where

(3.38) \begin{equation}\Delta_k=\frac{\Gamma\!\left(\frac{\alpha_2}{2}\right)\Gamma\left(\frac{1-\alpha_2}{2}\right)\left(k+\frac{\alpha_2}{2}\right)}{\sqrt{\pi}},\end{equation}

$H_{n+2l}^{(m)}$ is obtained by interchanging n and m in (3.32), while $H_{n+2l}^{(n)}$ will coincide with (3.29) if m is replaced by n. To sum up, for all $n,m=0,1,\ldots$ , $L_{nm}=L_{mn}\ne 0$ if $n-m$ is even and $L_{nm}=0$ otherwise.

We remark that owing to the asymptotic relations (3.31) and (3.33) and formula (3.38), the coefficients in the series (3.37) behave for large l as $l^{2(\alpha_2-\alpha_1)-3}$ $(\alpha_1>\alpha_2)$ , and therefore, the series representations (3.37) for the coefficients $L_{nm}$ rapidly converge.

On passing to the limit $\alpha_2\to \alpha_1,$ we can show that the matrix of the infinite system is diagonal and the system admits an exact solution that coincides with that associated with the contact problem of two bodies with the same exponent $\alpha_1=\alpha_2$ . Indeed, when $\alpha_1=\alpha_2$ from (3.29) and (3.34) we deduce that in either case, $l>0$ or $n\ne m$ , the coefficients $H_{n+2l}^{(m)}$ and $H_{m+2l}^{(n)}$ are equal to zero, and the only nonzero coefficients are $H_n^{(n)}$ . They are given by

(3.39) \begin{equation}H_n^{(n)}=\frac{\sqrt{\pi}\Gamma(\frac{\alpha_1+1}{2})(\alpha_1)_n}{\Gamma(\frac{\alpha_1}{2})(\frac{\alpha_1}{2}+n)n!}.\end{equation}

This gives a simple formula for the coefficients $L_{nm}$ . It is $L_{nm}=[H_n^{(n)}]^2\Delta_n\delta_{nm}$ , and from (3.36), $R_{nm} =\delta_{nm}$ . The system (3.35) has a diagonal matrix, and the coefficients $\Phi_{n}^{(j)}=(1+\gamma)^{-1}d_n^{(j)}$ are the same as those obtained by solving the integral equation (3.4) when $\alpha_1=\alpha_2$ on using the standard method of orthogonal polynomials.

In the general case, when $0<\alpha_2<\alpha_1<1$ , the infinite system (3.35) does not admit an exact solution. Its approximate solution is found by the reduction method. The off-diagonal elements of the matrix of the system $\delta_{mn}+\gamma R_{mn}$ rapidly decay, and the numerical method demonstrates a rapid convergence.

The right-hand sides of the system (3.35) are represented by the integrals (3.17). The integral $g_n^{(2)}$ is evaluated immediately, $g_n^{(2)}=\Gamma_0\delta_{n0}$ , where

(3.40) \begin{equation}\Gamma_0=\frac{\sqrt{\pi}\Gamma(\frac{\alpha_1+1}{2})}{\Gamma(\frac{\alpha_1}{2}+1)}.\end{equation}

The other integral $g_n^{(1)}$ can be computed explicitly if we know the coefficients $a_k$ of the expansion of the function f(bt) in terms of the Gegenbauer polynomials

(3.41) \begin{equation} f(bt)=\sum_{k=0}^\infty a_k C_k^{\alpha_1/2}(t).\end{equation}

These coefficients are always computed exactly if the function f(bt) is a polynomial. Otherwise, we can employ either its approximate polynomial representation or use the corresponding Gauss’ quadrature formula. In the polynomial case, when all the coefficients $a_k=0$ , $k>N$ , we apply the orthogonality property (3.8) to find $g_n^{(1)}=-a_n h_n(\alpha_1)$ , $n=0,1,\ldots,N$ , and $g_n^{(1)}=0$ , $n>N$ .

4 Solution of the contact problem

4.1 Parameter $\boldsymbol\delta$ , the contact zone $(-{\boldsymbol{b}},{\boldsymbol{b}})$ , the contact pressure p(x), and the normal displacements ${\boldsymbol{v}}_{\boldsymbol{j}}$

After the system (3.35) for the two right-hand sides, $d_n^{(1)}$ and $d_n^{(2)}$ have been solved, and the values of the coefficients $\Phi_n^{(1)}$ and $\Phi_n^{(1)}$ have been found; we write down the series representations (3.9) of the solutions $\phi^{(1)}(t)$ and $\phi^{(2)}(t)$ of the integral equations (3.4). On substituting these series into (3.5), we can express the unknown parameter $\delta$ through the coefficients $\Phi_0^{(1)}$ and $\Phi_0^{(2)}$

(4.1) \begin{equation}\delta=\frac{P/b- \Phi_0^{(1)}\Gamma_0}{\Phi_0^{(2)}\Gamma_0}.\end{equation}

On having this parameter, we can write down the contact pressure as

(4.2) \begin{equation}p(x)=\phi^{(1)}\left(\frac{x}{b}\right)+\delta \phi^{(2)}\left(\frac{x}{b}\right).\end{equation}

Notice that the parameter $\gamma=\alpha_1\theta_2(\alpha_2\theta_1)^{-1}b^{\alpha_1-\alpha_2}$ and the right-hand sides of the system (3.35) depend on the unknown parameter b. That is why the contact pressure also depends on this parameter. Because of the smoothness of the bodies’ profiles, the contact pressure has to be bounded at the points $x=\pm b$ , $y=0$ . Owing to the representations (3.9), this implies that the contact pressure vanishes at these points,

(4.3) \begin{equation}\lim_{t\to 1}[\phi^{(1)}(t)+\delta\phi^{(2)}(t)]=0.\end{equation}

Equivalently, this reads

(4.4) \begin{equation}\sum_{n=0}^\infty \frac{(\alpha_1)_n}{n!}\left(\Phi_n^{(1)}+\frac{P/b-\Phi_0^{(1)}\Gamma_0}{\Phi_0^{(2)}\Gamma_0}\Phi_n^{(2)}\right)=0.\end{equation}

This is a transcendental equation with respect to the parameter b. On having solved this equation, we can determine the parameter $\delta$ and the contact pressure by formulas (4.1) and (4.2), respectively.

The final quantities we wish to determine are the displacements $v_j(x)$ of the surface points outside the contact zone. We assume that the curvatures of the surfaces of interest are sufficiently small. Since formula (2.3) for the normal displacement is valid not only in the contact area but also outside, we can write

(4.5) \begin{equation}v_j(tx)=A_j\int_{-1}^1\frac{p(b\tau)d\tau}{|\tau-t|^{\alpha_j}}, \quad |t|>1, \quad j=1,2.\end{equation}

Using formula (3.3) and substituting the series representations (3.9) into (4.5), we write the displacements as follows:

(4.6) \begin{equation}v_j(x)=A_j\sum_{n=0}^\infty[\Phi_n^{(1)}+\delta\Phi_n^{(2)}] I_{n}\left(\frac{x}{b};\alpha_j\right), \quad |t|>1,\end{equation}

where

(4.7) \begin{equation}I_{n}(t;\alpha_j)=\int_{-1}^1\frac{(1-\tau^2)^{(\alpha_1-1)/2} C_n^{\alpha_1/2}(\tau)d\tau}{|\tau-t|^{\alpha_j}}.\end{equation}

Series representations of this integral are derived in Appendix B. Since the function f(x) is even, all the coefficients $\Phi_{2m+1}^{(j)}=0$ , $m=0,1,\ldots,$ $j=1,2$ , and therefore

(4.8) \begin{equation}v_j(x)=A_j\sum_{n=0}^\infty[\Phi_{2n}^{(1)}+\delta\Phi_{2n}^{(2)}] I_{2n}\left(\frac{x}{b};\alpha_j\right), \quad |t|>1.\end{equation}

On differentiating these functions, we find out that the derivatives $v_j^{\prime}(x)$ are bounded at the points $x=\pm b$ if and only if the condition (4.4) is satisfied. In other words, if the contact zone parameter b is fixed by solving the transcendental equation (4.4), then not only the pressure p(x) vanishes at the endpoints but also the profiles of the contacting bodies are smooth at the endpoints.

4.2 Numerical results

The functions $f_1(x)$ and $f_2(x)$ have to be continuously differentiable and satisfy the conditions $f_j(0)=f^{\prime}_j(0)=0$ , $j=0,1$ . In the symmetric case, when both of the functions are even, in a neighbourhood of the point $x=0$ ,

(4.9) \begin{equation}f_j(x)=\frac{f_j^{\prime\prime}(0)}{2}x^2+\frac{f^{(IV)}(0)}{24}x^4+\ldots, \quad j=1,2.\end{equation}

For numerical tests, we confine ourselves to two polynomial cases of the function $f(x)=f_1(x)+f_2(x)$ . They are

1. (1) $f(x)=Q_0x^2$ , $Q_0>0$ , and

2. (2) $f(x)=Q_0x^2+Q_1x^4$ .

Case (1) occurs when one of the bodies has a parabolic profile, while the second one is either flat or also has a parabolic profile. In case (2), the profiles of the bodies are described by the polynomials $f_j(x)=c_{0j} x^2+c_{1j}x^4$ with some real coefficients $c_{0j}$ and $c_{1j}$ chosen such that $Q_0=c_{01}+c_{02}\ge 0$ and $Q_1=c_{11}+c_{12}> 0$ .

In case (1), we express the function f(x) through the degree-0 and 2 Gegenbauer polynomials to have

(4.10) \begin{equation}f(bt)=\frac{b^2Q_0}{\alpha_1(\frac{\alpha_1}{2}+1)}\left[\frac{\alpha_1}{2}C_0^{\alpha_1/2}(t)+C_2^{\alpha_1/2}(t)\right].\end{equation}

The orthogonality property (3.8) yields $g_n^{(1)}=0$ for all n unless $n=0$ or $n=2$ . In these cases,

(4.11) \begin{equation}g_0^{(1)}=-\frac{b^2Q_0\Gamma_0}{\alpha_1+2}, \quad g_2^{(1)}=-\frac{b^2Q_0\sqrt{\pi}\Gamma(\frac{\alpha_1+3}{2})}{(\frac{\alpha_1}{2}+1)(\frac{\alpha_1}{2}+2)\Gamma(\frac{\alpha_1}{2})},\end{equation}

where $\Gamma_0$ is given by (3.40).

In case (2), the corresponding representation of the function f(bt) has the form

(4.12) \begin{align} f(bt)& =\frac{24Q_1 b^4 C_4^{\alpha_1/2}(t)}{\alpha_1(\alpha+2)(\alpha_1+4)(\alpha_1+6)} \nonumber\\[3pt] &\quad +\frac{2b^2C_2^{\alpha_1/2}(t)}{\alpha_1(\alpha_1+2)}\left(Q_0+\frac{6Q_1b^2}{\alpha_1+6}\right)+\frac{b^2C_0^{\alpha_1/2}(t)}{\alpha_1+2}\left(Q_0+\frac{3Q_1b^2}{\alpha_1+4}\right),\end{align}

Except for $g_0^{(1)}$ , $g_2^{(1)}$ and $g_4^{(1)}$ , all the terms $g_n^{(1)}$ equal 0. The nonzero terms are given by

(4.13) \begin{align} g_0^{(1)}&=-\frac{b^2\Gamma_0}{\alpha_1+2}\left(Q_0+\frac{3Q_1b^2}{\alpha_1+4}\right), \nonumber \\[3pt] g_2^{(1)}&= -\frac{b^2\sqrt{\pi}\alpha_1\Gamma\!\left(\frac{\alpha_1+3}{2}\right)}{2\Gamma\!\left(\frac{\alpha_1}{2}+3\right)}\left(Q_0+\frac{6Q_1b^2}{\alpha_1+6}\right), \nonumber \\[3pt] g_4^{(1)}&= -\frac{b^4\sqrt{\pi}\alpha_1\left(\alpha_1+2\right)\Gamma\!\left(\frac{\alpha_1+5}{2}\right)}{4\Gamma\!\left(\frac{\alpha_1}{2}+5\right)}Q_1.\end{align}

For the numerical tests to be discussed, we choose the resultant force and the Poisson ratios to be $P=1$ , $\nu_1=\nu_2=0.3$ , the resultant moment to be zero and the function $f(x)=f_1(x)+f_2(x)$ to be even. This choice gives rise to a solution symmetric with respect to the y-axis. Figure 2 presents the half-length b of the contact zone for different values of the exponents $\alpha_1$ and $\alpha_2$ in the Young moduli of the bodies, $E_1=e_1y^{\alpha_1}$ and $E_2=e_2(-y)^{\alpha_2}$ when $e_1=e_2=1$ and (a) $f(x)=x^2$ and (b) $f(x)=x^4$ . It is seen that when $\alpha_1$ is fixed and $\alpha_2$ increases in the interval $(0,\alpha_1)$ , the contact zone length is also increasing. The same is true in the case when $\alpha_2$ is fixed and $\alpha_1$ grows. On comparing the results presented in Figure 2(a) and (b), we see that when the curvatures of the contacting bodies profiles are decreasing, the contact zone length is increasing.

Figure 2. The half-length b of the contact zone $(-b,b)$ versus the parameter $\alpha_2\in(0,\alpha_1)$ for $\alpha_1=0.5$ , $\alpha_1=0.7$ and $\alpha_1=0.9$ when (a) $f(x)=x^2$ ( $Q_0=1$ , $Q_1=0$ ) and (b) $f(x)=x^4$ ( $Q_0=0$ , $Q_1=1$ ).

Curves in Figure 3 give a clear demonstration of the dependence of the length of the contact zone upon one of the factors $e_1$ and $e_2$ while the second one is kept fixed. The parameters for this diagram are chosen as $e_2=1$ , $f(x)=x^2$ , $\alpha_2=0.3$ , and $\alpha_1$ is equal to either $0.5$ , $0.7$ , or $0.9$ .

Figure 3. The half-length b of the contact zone $(-b,b)$ versus the parameter $e_1\in(0,5)$ for $\alpha_1=0.5$ , $\alpha_1=0.7$ and $\alpha_1=0.9$ when $e_2=1$ , $\alpha_2=0.3$ , $f(x)=x^2$ .

Figure 4 shows how the parameter $\delta$ depends on the second body exponent $\alpha_2\in(0,\alpha_1)$ when the exponent $\alpha_1$ is fixed and chosen to have the values $0.5$ , $0.9$ or $0.95$ . The other parameters are $e_1=e_2=1$ , and the function $f(x)=x^2$ . It is seen that the parameter $\delta$ increases as $\alpha_2\to 0$ and also when $\alpha_1\to 1$ and $\alpha_2\to \alpha_1$ .

Figure 4. The parameter $\delta$ versus the exponent $\alpha_2\in(0,\alpha_1)$ for $\alpha_1=0.5$ , $\alpha_1=0.9$ and $\alpha_1=0.95$ when $f(x)=x^2$ .

The results of calculations of the pressure distribution p(x) are shown in Figures 5 and 6. In both cases, $e_1=e_2=1$ , and $f(x)=x^2$ . In Figure 5, the smaller exponent $\alpha_2$ is fixed as $\alpha_2=0.1$ , while $\alpha_1$ is equal to either $0.3$ , $0.7$ or $0.9$ . The corresponding values of the half-length b of the contact zone are computed to be $1.17365$ , $1.43214$ and $ 1.92390$ . The contact pressure p(x) vanishes at the endpoints $\pm b$ of the contact zone and attains its maximum at the origin. As the parameter $\alpha_1$ is increasing, the pressure maximum is decreasing. When the bigger exponent $\alpha_1$ is fixed (in Figure 6, $\alpha_1=0.95$ ), while the smaller exponent varies in the interval $(0,\alpha_1)$ , the variation of the pressure distribution p(x) for a fixed x is not large (Figure 6).

Figure 5. The contact pressure p(x), $x\in[0,b]$ for $\alpha_2=0.1$ when $\alpha_1=0.3$ , $\alpha_1=0.7$ , and $\alpha_1=0.9$ in the case $f(x)=x^2$ .

Figure 6. The contact pressure p(x), $x\in[0,b]$ for $\alpha_1=0.95$ when $\alpha_2=0.9$ , $\alpha_2=0.8$ , and $\alpha_2=0.1$ in the case $f(x)=x^2$ .

The normal elastic displacements $u_y(x,0)=v_1(x)$ and $u_y(x,0)=-v_2(x)$ of the upper and lower elastic bodies outside the contact zone are shown in Figure 7 (the displacements of the lower body $B_2$ are demonstrated by broken curves). As before, $e_1=e_2=1$ , $f(x)=x^2$ , and the functions $v_1(x)$ and $v_2(x)$ are even. For computations, we choose $\alpha_1$ to be either $0.5$ , $0.7$ or $0.9$ , while $\alpha_2=\alpha_1/2$ . Both displacements attain their maximum at the points $\pm b$ . It has been numerically verified that

(4.14) \begin{equation}\lim_{x\to\pm b^\pm} [v_1(x)+v_2(x)]=\delta-f_1(b)-f_2(b)\end{equation}

that is consistent with the boundary condition (2.1). The corresponding values of the half-length of the contact zone are $b= 1.28951$ for $\alpha_1=0.5$ , $b= 1.46450$ when $\alpha_1=0.7$ , and $b=1.94635$ in the case $\alpha_1=0.9$ .

Figure 7. The normal displacements $v_1(x)$ (solid curves) and $v_2(x)$ (broken curves), $x\in[-10,-b]$ for $\alpha_2=\alpha_1/2$ when $\alpha_1=0.5$ , $\alpha_1=0.7$ , and $\alpha_1=0.9$ in the case $f(x)=x^2$ .

Figure 8. The case $\alpha_1=\alpha_2=\alpha$ when $P=1$ , $f(x)=x^2$ , $e_1=e_2=1$ , $\nu_1=\nu_2=0.3$ . (a): the contact zone half-length b versus $\alpha\in(0,1)$ , (b): the parameter $\delta$ versus $\alpha$ , (c): pressure p(x) for $\alpha=0.3$ , (d): the displacement $u_y(x,0)=-v_2(x)$ for $x<-b$ and $u_y(x,0)=(x^2-\delta)/2$ for $x\in(-b,0)$ when $\alpha=0.3$ .

5 Hertzian contact of two power-law graded bodies when $\boldsymbol{E}_{\boldsymbol{1}}=\boldsymbol{e}_{\boldsymbol{1}}\boldsymbol{y}^{\boldsymbol\alpha}$ and $\boldsymbol{E}_{\boldsymbol{2}}=\boldsymbol{e}_{\boldsymbol{2}} (-\boldsymbol{y})^{\boldsymbol\alpha}$

Assume that the contacting bodies $B_1$ and $B_2$ have the Young moduli $E_1=e_1 y^\alpha$ and $E_1=e_2 (-y)^\alpha$ . The governing equation (2.5) with two kernels reduces to

(5.1) \begin{equation}A\int_{-1}^1\frac{p(b\tau)d\tau}{|\tau-t|^{\alpha}}=\delta-f_1(bt)-f_2(bt), \quad -1<t<1,\end{equation}

where $A=\alpha^{-1}(\theta_1+\theta_2)b^{1-\alpha}$ , and $\theta_j$ are defined by (2.4) with $\alpha_1=\alpha_2=\alpha$ . The pressure distribution has to be an even function, and we represent the solution in the form

(5.2) \begin{equation}p(bt)=(1-t^2)^{(\alpha-1)/2}\sum_{n=0}^\infty\Phi_{2n}C_{2n}^{\alpha/2}(t), \quad -1<t<1.\end{equation}

On substituting this function into (5.1), using the spectral relation (3.7) and orthogonality property (3.8) we find the coefficients $\Phi_{2n}$

(5.3) \begin{equation}\Phi_{2n}=\frac{g_{2n}^{(1)}+\delta\Gamma_0\delta_{n0}}{A\beta_{2n}(\alpha)h_{2n}(\alpha)},\end{equation}

where

(5.4) \begin{equation}g_{2n}^{(1)}=-\int_{1}^1 C_{2n}^{\alpha/2}(t)(1-t^2)^{(\alpha-1)/2}f(bt)dt\end{equation}

and $f(x)=f_1(x)+f_2(x)$ . To determine the parameter $\delta$ , we satisfy the equilibrium condition (2.2) and obtain

(5.5) \begin{equation}\delta=\frac{1}{\Gamma_0}\left(-g_0^{(1)}+\frac{\pi AP}{b\cos\frac{\pi\alpha}{2}}\right).\end{equation}

The half-length b of the contact zone is the positive root of the following transcendental equation (our numerical tests reveal that such a root is unique):

(5.6) \begin{equation}\sum_{n=0}^\infty \frac{(\alpha)_{2n}}{(2n)!}\Phi_{2n}=0\end{equation}

that reduces to

(5.7) \begin{equation}\frac{\delta\Gamma_0\alpha}{2\Gamma(\alpha)}+\sum_{n=0}^\infty \frac{g_{2n}^{(1)}(2n)!(2n+\frac{\alpha}{2})}{\Gamma(\alpha+2n)}=0.\end{equation}

The normal surface displacements of the bodies $B_1$ and $B_2$ are expressed through the integral $I_{2n}(x/b;\alpha)$

(5.8) \begin{equation}v_j(x)=A_j\sum_{n=0}^\infty \Phi_{2n} I_{2n}\left(\frac{x}{b};\alpha\right), \quad |x|>b,\quad j=1,2,\end{equation}

where

(5.9) \begin{equation}A_j=\frac{\theta_jb^{1-\alpha}}{\alpha}, \quad I_{2n}(t;\alpha)=\int_{-1}^1\frac{(1-\tau^2)^{(\alpha-1)/2} C_{2n}^{\alpha/2}(\tau)d\tau}{|\tau-t|^{\alpha}}.\end{equation}

This integral is a particular case $\alpha_j=\alpha$ of the integral $I_n(t;\alpha_j)$ evaluated in Appendix B and given by (B4) and (B5). As in the case $\alpha_1>\alpha_2$ , the displacements $v_j$ and their first derivative are bounded as $x\to \pm b^\pm$ , and the contacting surfaces are smooth at the endpoints. For numerical purposes, the Gauss quadrature order-N formula can also be employed

(5.10) \begin{equation}v_j(x)=\frac{\pi A_j}{N}\sum_{i=1}^N\frac{p(bx_i)}{|x/b-x_i|^\alpha}\sin\frac{(2i-1)\pi}{2N}, \quad x_i=\cos\frac{(2i-1)\pi}{2N},\quad |x|>b.\end{equation}

The numerical tests show that in the case of Hertzian contact, when the pressure vanishes at the endpoints, this approximation is in good agreement with the exact formulas (B4) and (B5).

Consider the particular case $f(x)=Q_0x^2$ . Owing to the fact that $g_n^{(1)}=0$ for all n except for $g_0^{(1)}$ and $g_2^{(1)}$ and employing formulas (4.11) for these nonzero terms, we specify the formulas for the parameter $\delta$ and find explicitly the half-length of the contact zone

(5.11) \begin{equation}\delta=\frac{Q_0b^2}{\alpha+2}+\frac{\pi AP}{\Gamma_0 b\cos\frac{\pi\alpha}{2}},\quad b=\left(\frac{\Gamma(2+\frac{\alpha}{2})\Gamma(\frac{1-\alpha}{2})(\theta_1+\theta_2)P}{\sqrt{\pi} Q_0}\right)^{\frac{1}{\alpha+2}}.\end{equation}

For this parabolic case, we also compute the pressure distribution and the normal displacement. From (5.2) we have

(5.12) \begin{equation}p(x)=\left(1-\frac{x^2}{b^2}\right)^{(\alpha-1)/2}\left[\frac{P}{\Gamma_0 b}+\frac{2b^2Q_0\cos\frac{\pi\alpha}{2}}{\pi A\alpha(\alpha+1)}\left(\frac{1}{\alpha+2}-\frac{x^2}{b^2}\right)\ \right].\end{equation}

On substituting the expression (5.11) for b into formula (5.12), we arrive at

(5.13) \begin{equation}p(x)=\frac{P\Gamma(\frac{\alpha}{2}+2)}{\sqrt{\pi}b\Gamma(\frac{\alpha+3}{2})}\left(1-\frac{x^2}{b^2}\right)^{(\alpha+1)/2}.\end{equation}

In the particular case, when one of the bodies is a rigid punch, this formula coincides with the corresponding expression of the pressure distribution derived in [9]. When $\alpha\to 0$ , the contact zone half-length b and the pressure distribution p(x) tend to $b_0$ and $p_0(x)$ which represent the contact zone half-length and the contact pressure, respectively, in the case when both bodies are isotropic elastic bodies and whose contact is governed by equation (A3)

(5.14) \begin{align} \lim_{\alpha\to 0} b&=b_0, \quad b_0=\sqrt{\frac{(\theta^\circ_1+\theta_2^\circ)P}{Q_0}}, \nonumber \\[3pt] \lim_{\alpha\to 0} p(x)&=p_0(x), \quad p_0(x)=\frac{2P}{\pi b_0^2}\sqrt{b_0^2-x^2},\end{align}

where $\theta_j^\circ$ are given by (A2). This expression coincides with the contact pressure associated with the elastic isotropic case of the Hertzian model and obtained by solving equation (A3) ([29], Chapter II, (23)).

We determine the normal displacements of surface points $u_y(x,0)=-v_2(x)$ of the lower body $B_2$ outside the contact zone when $B_2$ is a half-plane that is when $f_2(x)=0$ , while the upper body $B_1$ has the profile $y=Q_0x^2$ . As before, the Young moduli of both bodies have the same exponent and may have different factors $e_1$ and $e_2$ . Since the only nonzero coefficients are $\Phi_0$ and $\Phi_2$ , we transform formula (5.8) to the form

(5.15) \begin{equation}v_2(bt)=\theta_2 b^{1-\alpha}[\Phi_0 \tilde I_0(t;\alpha)+\Phi_2 \tilde I_2(t;\alpha)],\end{equation}

where

(5.16) \begin{equation}\Phi_0=\frac{P}{b\Gamma_0}, \quad \Phi_2=-\frac{4b^{\alpha+1}Q_0\cos\frac{\pi\alpha}{2}}{\pi\alpha(\alpha+1)(\alpha+2)(\theta_1+\theta_2)}, \quad\tilde I_n(t;\alpha)=\frac{1}{\alpha} I_n(t;\alpha).\end{equation}

From (B4),

(5.17) \begin{align} \tilde I_0(t;\alpha)=\frac{\pi}{\cos\frac{\pi\alpha}{2}}\left[\frac{1}{\alpha}-\frac{\Gamma(\frac{\alpha+1}{2})(-\frac{t+1}{2})^{(1-\alpha)/2}}{\Gamma(\alpha+1)\Gamma(\frac{3-\alpha}{2})}F\left(\frac{1-\alpha}{2},\frac{1+\alpha}{2},\frac{3-\alpha}{2};\frac{t+1}{2}\right)\right], \nonumber\\[3pt] \tilde I_2(t;\alpha)=\frac{\pi\alpha(\alpha+1)}{2\cos\frac{\pi\alpha}{2}}\left[\frac{\alpha+1}{2}F\left(-2,\alpha+2,\frac{\alpha+1}{2};\frac{t+1}{2}\right)\right. \nonumber\\[3pt] \left.-\frac{\Gamma(\frac{\alpha+1}{2})(-\frac{t+1}{2})^{(1-\alpha)/2}}{\Gamma(\alpha+1)\Gamma(\frac{3-\alpha}{2})}F\left(-\frac{3+\alpha}{2},\frac{5+\alpha}{2},\frac{3-\alpha}{2};\frac{t+1}{2}\right)\right],\quad -3<t<-1.\end{align}

For $t<-3,$ the integrals $\tilde I_n(t;\alpha)=\alpha^{-1} I_n(t;\alpha)$ ( $n=0,2$ ) are obtained from formula (B5). It is easy to see from formulas (5.15) to (5.17) that the displacements $v_j(t;\alpha)$ become infinite when $\alpha\to 0$ , and the limit transition $\alpha\to 0$ for the displacements outside the contact zone is impossible.

Figure 8 shows the results of computations in the case when the Young moduli of the contacting bodies are the same, $E_j(z)=e_j z^\alpha$ and $e_1=e_2$ . We choose $P=1$ , $f(x)=x^2$ , $e_1=e_2=1$ , and $\nu_1=\nu_2=0.3$ . Figure 8(a) and (b) demonstrate the variation of the half-length b and the parameter $\delta$ with the exponent $\alpha\in(0,1)$ . As $\alpha$ grows, the contact zone becomes larger. As in the case $\alpha_1\ne \alpha_2$ , as $\alpha\to 0$ , the parameter $\delta\to\infty$ . It also grows as $\alpha\to 1$ . In Figure 8(c) and (d), we present sample curves for the contact pressure for $x\in[0,b]$ and the normal displacement $u_y(x,0)=-v_2(x)$ of the surface of the lower body (a half-plane) when $f_1(x) =x^2$ , $f_2(x)=0$ for $x<-b$ and $u_y(x,0)=(x^2-\delta)/2$ for $x\in(-b,0)$ . In both Figure 8(c) and (d), $\alpha=0.3$ (in this case $b= 1.22072$ ).

6 Surface energy model

In this section following the JKR model [16, 17], we aim to take into account the effect of adhesive forces (Figure 1(b)) and study their impact on the contact zone size, the contact pressure and the normal displacement. In the two-dimensional case, the loss of surface energy is given by $U_s=-2\gamma_s b$ , where $\gamma_s$ is the work of adhesion (a half-density of the surface energy). The elastic strain energy is expressed through the normal displacement $v(x)=\delta-f(x)$ , and the contact pressure p(x) as

(6.1) \begin{equation}U_e=\frac12\int_{-b}^b p(x) v(x)dx,\end{equation}

and the total energy defined by $U_{total}=U_e-2\gamma_s b$ is a function of the contact zone half-length b. In contrary to the Hertzian model, the JKR model admits singularities of the contact pressure at the endpoints. Also, the parameter b is defined not from the condition that quenches the pressure singularities but from the condition of minimum of the total energy that is

(6.2) \begin{equation}\frac{d U_e}{db} -2\gamma_s=0.\end{equation}

6.1 Case $\boldsymbol\alpha_{\boldsymbol{1}}=\boldsymbol\alpha_{\boldsymbol{2}}=\boldsymbol\alpha$

We consider parabolic profiles of the contacting bodies, $f(x)=Q_0x^2$ . In this case, the pressure p(x) is given by (5.12) and has order $(\alpha-1)/2$ power singularities at the endpoints, while the parameter b is free. Since the resultant force P has to be balanced by the contact pressure p(x), we satisfy the condition (2.2) and define $\delta$ by formula (5.11). To evaluate the integral (6.1), we write the pressure and displacement in the following equivalent form:

(6.3) \begin{align} p(bt)=(1-t^2)^{(\alpha-1)/2}[\Phi_0 C_0^{\alpha/2}(t)+\Phi_2 C_2^{\alpha/2}(t)], \nonumber \\[3pt] v(bt)=\delta C_0^{\alpha/2}(t)-\frac{Q_0 b^2}{\alpha(\alpha+2)}[\alpha C_0^{\alpha/2}(t)+2 C_2^{\alpha/2}(t)],\end{align}

where $\Phi_0$ and $\Phi_2$ are determined in (5.16). Using the orthogonality property (3.8) of the Gegenbauer polynomials, we obtain

(6.4) \begin{equation}U_e=\frac{P^2(\theta_1+\theta_2)\Gamma(\frac{\alpha}{2}+1)\Gamma(\frac{1-\alpha}{2})b^{-\alpha}}{2\alpha\sqrt{\pi}}+\frac{\sqrt{\pi} Q_0^2 b^{\alpha+4}}{2(\theta_1+\theta_2)(\alpha+2)\Gamma(\frac{1-\alpha}{2})\Gamma(\frac{\alpha}{2}+3)}.\end{equation}

The derivative of the elastic strain energy in (6.2) can be evaluated exactly, and we arrive at the following transcendental equation with respect to the parameter b:

(6.5) \begin{equation}\frac{\sqrt{\pi} Q_0^2 b^{2\alpha+4}}{(\theta_1+\theta_2)(\alpha+2)\Gamma(\frac{1-\alpha}{2})\Gamma(\frac{\alpha}{2}+2)}-2\gamma_s b^{\alpha+1}-\frac{P^2(\theta_1+\theta_2)\Gamma(\frac{\alpha}{2}+1)\Gamma(\frac{1-\alpha}{2})}{2\sqrt{\pi}}=0.\end{equation}

Passing to the limit $\alpha\to 0$ and keeping $\gamma_s\ge 0,$ we reduce the transcendental equation to the quartic equation

(6.6) \begin{equation}\frac{Q_0^2 b^{4}}{2(\theta^\circ_1+\theta^\circ_2)}-2\gamma_s b-\frac{P^2(\theta^\circ_1+\theta^\circ_2)}{2}=0,\end{equation}

and when, in addition, $\gamma_s\to 0$ , we obtain the classical formula (5.15) for the value of b in the case of Hertzian contact of two elastic isotropic bodies.

Passing to the limit $\gamma_s\to 0$ in equation (6.5) and keeping $\alpha\in(0,1),$ we arrive at the equation with respect to b that admits an exact solution; it coincides with the value of b in the Hertzian model given by (5.11).

If we assume that $\nu_1=\nu_2=0.5$ , by passing to the limit $\alpha\to 1$ we derive from (6.5) the following cubic equation with respect to $b^2$ for two Gibson solids [10]:

(6.7) \begin{equation}\frac{8}{27}Q_0^2 b^6-2\gamma_s\left(\frac{1}{e_1}+\frac{1}{e_2}\right)b^2-\frac{3P^2}{8}\left(\frac{1}{e_1}+\frac{1}{e_2}\right)^2=0.\end{equation}

Notice that the transcendental equation (6.5) is different from the corresponding equations obtained in [9] and [6]. These authors split the solution into two parts, the first one gives the solution for a parabolic punch and the second one corresponds to the model of a flat punch. The discrepancies between the transcendental equations in [9] and [6] and equation (6.5) are caused by their disregard for the fact that the displacement $\delta$ is a function of the contact zone half-length (the contact radius) b. This explains why the limit transition $\alpha\to 0$ is impossible in the solutions [9] and [6]. We emphasise that the contact radius and the stresses recovered from the 2d contact model of two power-law graded bodies or a stamp and a power-law graded body in the limit $\alpha\to 0$ must give the corresponding quantities of the 2d contact model of homogeneous bodies governed by the integral equation with a logarithmic kernel (A3). Otherwise, the solution of the model of power-law graded bodies is incorrect.

A number of tests have been conducted to ascertain the impact of the surface energy density $2\gamma_s$ and the exponent $\alpha$ on the contact zone size 2b, the pressure distribution and the normal displacement. The curves in Figure 9(a) exhibit an increase of the contact zone size with the parameter $\gamma_s$ . It is seen from Figure 9(b) that the half-length b rapidly increases when the exponent $\alpha$ approaches 1. The $P-b$ curves in Figure 9(c) demonstrate the rate of growth of the half-length b when the total force P grows. From Figure 9(d), it is seen that when the Young moduli are $E_j=e_j |y|^\alpha$ , $e_2=1$ and the factor $e_1$ grows, the parameter b first rapidly decreases, and then, its rate of decrease is insignificant.

Figure 9. JKR model: the contact zone half-length b in the case $\alpha_1=\alpha_2=\alpha$ , $f(x)=x^2$ , $e_2=1$ , $\nu_1=\nu_2=0.3$ . (a): b versus the surface energy density $\gamma_s$ , (b): b versus the exponent $\alpha$ , (c): b versus the normal force P when $\gamma_s=1$ , (d): b versus the parameter $e_1$ when $\gamma_s=1$ .

Contact pressure curves computed according to the Hertz and JKR models are portrayed in Figure 10(a). Since the pressure p(x) is an even function, the curves demonstrate that the pressure vanishes at the endpoints $x=\pm b$ in the former model. In the JKR model, the contact stress is compressive for $-b_*<x<b_*$ and tensile at the edge zones $(-b,-b_*)$ and $(b_*,b)$ . The numerical tests show that a growth of the surface energy density $2\gamma_s$ shrinks the central zone, where the stress is compressive, and enlarges the zone, where the stress is tensile.

Figure 10. (a): the contact pressure p(x) and the normal displacement $u_y(x,0)=-v_2(x)$ for $x<-b$ and $u_y(x,0)=(x^2-\delta)/2$ for $x\in(-b,0)$ for the Hertzian ( $\gamma_s=0$ ) and JKR ( $\gamma_s=1$ ) models when $\alpha=0.5$ , $P=1$ , $f(x)=x^2$ , $e_1=e_2=1$ , $\nu_1=\nu_2=0.3$ .

As in the Hertzian contact model, we analyse the normal displacements $u_y(x,0)=-v_2(x)$ for the JKR model outside of the contact zone when the upper body has a parabolic profile, $f_1(x)=Q_0 x^2$ , and the lower body is a half-plane, $f_2(x)=0$ . The displacement $v_2(x)$ is given by the same formula (5.15). However, since the pressure p(x) does not vanish at the endpoints, has power singularities and is described by formula (5.12), the derivative of the displacement (5.15) tends to infinity as $x\to\pm b$ . A sample curve of the displacement $u_y(x,0)=-v_2(x)$ for $x<-b$ and $u_y(x,0)=(x^2-\delta)/2$ for $x\in(-b,0)$ ( $Q_0=1$ ) is shown in Figure 10(b). It is seen that in the case of Hertzian contact ( $\gamma_s=0),$ the contact surface is smooth near the contact zone endpoints, while in the case of the JKR model, due to the adhesion forces a part of the surface of the flat body $B_2$ is attracted to the interface, and the tangent lines to the surfaces of the contacting bodies at the endpoints have different slopes.

6.2 Case $\boldsymbol\alpha_{\boldsymbol{1}}>\boldsymbol\alpha_{\boldsymbol{2}}$

The direct method for computing the elastic strain energy described in the previous section can be generalised to the case when the contacting bodies have different exponents and as before, $\alpha_1>\alpha_2$ . On expanding the normal displacement $v(x)=\delta-f(x)$ $(-b<x<b)$ through the Gegenbauer polynomials of even order, we have

(6.8) \begin{equation}v(bt)=\delta-\sum_{k=0}^\infty a_{2k} C_{2k}^{\alpha_1/2}(t), \quad -1<t<1,\end{equation}

substituting it together with the contact pressure

(6.9) \begin{equation}p(bt)=(1-t^2)^{(\alpha_1-1)/2}\sum_{n=0}^\infty[\Phi_{2n}^{(1)}+\delta\Phi_{2n}^{(2)}]C_{2n}^{\alpha_1/2}(t), \quad -1<t<1,\end{equation}

into formula (6.1) and using the orthogonality of the polynomials we derive the series representation of the elastic strain energy

(6.10) \begin{equation}U_e=\frac{b\delta}{2}(\Phi_0^{(1)}+\delta\Phi_0^{(2)})\Gamma_0-\frac{b}{2}\sum_{n=0}^\infty(\Phi_{2n}^{(1)}+\delta\Phi_{2n}^{(2)})a_{2n} h_{2n}(\alpha_1).\end{equation}

As before, we simplify the formula for the parabolic case, $f(x)=Q_0x^2$ . In this case, $a_n=0$ unless $n=0$ or $n=2$ ,

(6.11) \begin{equation}a_0=\frac{b^2 Q_0}{\alpha_1+2}, \quad a_2=\frac{2 b^2 Q_0}{\alpha_1(\alpha_1+2)},\end{equation}

and ultimately, we have

(6.12) \begin{equation}U_e=\frac{b\Gamma_0}{2}(\Phi_0^{(1)}+\delta\Phi_0^{(2)})\left(\delta-\frac{b^2Q_0}{\alpha_1+2}\right)-\frac{2b^3Q_0\sqrt{\pi}\Gamma(\frac{\alpha_1+3}{2})}{(\alpha_1+2)(\alpha_1+4)\Gamma(\frac{\alpha_1}{2})}(\Phi_{2}^{(1)}+\delta\Phi_{2}^{(2)}).\end{equation}

The minimum of the total energy is attained if the contact zone half-length b solves the transcendental equation

(6.13) \begin{equation}\frac{dU_e}{db}-2\gamma_s=0.\end{equation}

Explicit differentiation is impossible for the coefficients $\Phi_0^{(j)}$ and $\Phi_2^{(j)}$ being a part of the solution to the infinite system (3.35), and there is no way to explicitly separate b from the unknowns of the infinite system. Approximately, equation (6.13) can be written as

(6.14) \begin{equation} \frac{U_e(b+\varepsilon)-U_e(b)}{\varepsilon}-2\gamma_s\approx 0,\end{equation}

where $\varepsilon$ is a small and positive.

The variation of the half-length of the contact zone b with the half-density $\gamma_s$ of the surface energy for three values of the exponent $\alpha$ is portrayed in Figure 11(a). It has been calculated by the method of orthogonal polynomials presented in Section 3. The difference between the scheme for the Hertzian and JKR models is only in the way how the parameter b is fixed. In the Hertzian model, it solves the transcendental equation (4.4) that guarantees that the pressure vanishes at the endpoints, while in the JKR model, it is defined from the approximate equation (6.14), the condition of minimum of the total energy. For computations, $\varepsilon$ is accepted to be $10^{-4}$ , and the differences between the results for $\varepsilon=10^{-3}, 10^{-4}, 10^{-5}$ are not significant. For example, for $\alpha_1=0.5$ , $\alpha_2=0.25$ , $e_1=e_2=1$ , $P=1$ , $Q_0=1$ , $\nu_1=\nu_2=0.3$ , and $\gamma_s=1,$ we have $b=1.97621$ if $\varepsilon=10^{-3}$ , $b=1.97666$ if $\varepsilon=10^{-4}$ , and $b=1.97670$ if $\varepsilon=10^{-5}$ . It turns out that as $\gamma_s\to 0$ , the pressure distribution vanishes at the endpoints, and the parameter b associated with the JKR model tends to the one for the Hertzian model. This is consistent with the homogeneous case $\alpha_1=\alpha_2$ considered in Section 6.1 where we managed to prove this analytically and with the result [2] obtained for an axisymmetric model of contact of a rigid stamp and a power-law graded half-space.

Figure 11. JKR model for the case of different exponents: $E_1=e_1 y^{\alpha_1}$ and $E_2=e_2 (-y)^{\alpha_2}$ , when $e_1=e_2=1$ , $P=1$ , $Q_0=1$ , $\nu_1=\nu_2=0.3$ . (a): the half-length b versus the half-density $\gamma_s$ of the surface energy for $\alpha_1=0.5, 0.7, 0.9$ and $\alpha_2=\alpha_1/2$ . (b): the contact pressure p(x) for $x\in(0,b)$ for $\gamma_s=0, 1, 5$ when $\alpha_1=0.9$ and $\alpha_2=0.5$ .

The pressure distribution p(x) is shown in Figure 11(b) for three values of the parameter $\gamma_s$ when $\alpha_1=0.9$ and $\alpha_2=0.5$ . As $\gamma_s\to 0$ , the contact pressure vanishes at the endpoints and coincides with the pressure found from the Hertzian model. When $\gamma_s>0$ , similarly to the case $\alpha_1=\alpha_2$ , the contact zone enlarges and the contact stress becomes tensile at two edge zones $(-b,-b_*)$ and $(b_*,b)$ and tends to $\infty$ as $x\to \pm b$ .

7 Conclusions

We analysed two plane problems, the Hertzian and JKR models, of frictionless contact of two inhomogeneous elastic bodies with distinct moduli of elasticity $E_1(y)=e_1 y^{\alpha_1}$ and $E_2(y)=e_2 (-y)^{\alpha_2}$ with $0<\alpha_2\le\alpha_1<1$ . On employing the Rostovtsev representation of the normal displacement in the contact zone through the pressure distribution, we showed that the model is governed by an integral equation with two different power kernels. For its solution, a novel method of Gegenbauer orthogonal polynomials was proposed. We reduced the integral equation to an infinite system of linear algebraic equations whose coefficients, after some transformations, become integral free. It was demonstrated that when $\alpha_2\to\alpha_1$ , the infinite system is decoupled, and its exact solution coincides with the one obtained by direct solution of the integral equation with one power kernel. When $\alpha_2\ne\alpha_1$ and the contact zone is semi-infinite, the integral equation was solved exactly by the Wiener-Hopf method. We also showed that the integral equation in a finite interval admits an exact solution for any $\alpha_1\in(0,1)$ when the second parameter $\alpha_2$ is either 0 or any negative even number.

We found a rigid body displacement $\delta$ (the total displacement of distant points of the bodies) from the equilibrium condition that balances the normal total force and the contact pressure. The length of the contact zone is determined from a transcendental equation that guarantees that the pressure vanishes at the endpoints in the Hertzian model and the total energy attains its minimum in the JKR model. The pressure distribution is found in a series form, and the coefficients of the expansion are determined from an infinite system of the second kind solved by the reduction method. The numerical tests implemented revealed rapid convergence of the method for all admissible values of the model parameters. By employing the method of Mellin’s convolution integrals and the theory of residues, we computed the normal displacements of surface points outside the contact zone. In the Hertzian model, the profile of the contacting surfaces at the endpoints is smooth, while in the JKR model, the derivative of the normal displacement is infinite, and a part of the contacting surfaces is attracted by adhesion forces to the interface. In contrary to the Hertzian model, the pressure distribution does not vanish at the endpoints, it tends to $-\infty$ , and there are two edge zones in the contact area where the contact stress is tensile.

Our numerical results showed that the parameter $\delta$ , the contact zone length, the contact pressure and the elastic displacement significantly depend on variation of the bigger parameter $\alpha_1$ and only slightly vary with the second, smaller, parameter $\alpha_2$ . When the exponent $\alpha_1$ is growing, the contact zone is also growing. In the case when the two exponents $\alpha_1$ and $\alpha_2$ are the same, $\alpha_1=\alpha_2=\alpha$ , we obtained the contact zone length, the parameter $\delta$ , the contact pressure and the normal displacements of the surface points exactly. By passing to the limit $\alpha\to 0$ , we showed that the result coincides with the classical solution of the problem of Hertzian contact of two isotropic elastic bodies.

For the JKR model, we found out that when the half-density of surface energy $\gamma_s\to 0$ , the contact zone length, pressure and normal displacement tend to the corresponding quantities associated with the Hertzian model. In both cases, $\alpha_1>\alpha_2$ and $\alpha_1=\alpha_2$ , the transcendental equation for the contact zone length admits passing to the limit $\alpha_j\to 0$ . This is possible not only for the contact zone length but also for the contact pressure. However, the normal displacement derived for both Hertzian and JKR models when $\alpha_1\ge \alpha_2>0$ become infinite when $\alpha_1\to 0$ . This is due to the presence of $\alpha_1^{-1}$ in the formula for the displacement $\delta$ of distant points of the contacting bodies.

The method we presented admits generalisations and modifications in different directions including the Hertzian and JKR axisymmetric contact models of two power-law graded bodies.