Axisymmetric frictionless contact between an elastic layer thickness and a circular base under a rigid punch

Abstract

The study presented in this work deals with analytical methods for an axisymmetric problem of an elastic layer partially reposing on a rigid circular base, and is indented along the upper surface with a rigid punch. The contact between the medium and the base is smooth. This boundary value problem is transformed into a system of dual integral equations. In contrast to the classical approach consisting in resolving the corresponding Fredholm equation of the second kind, the latter equations are obtained from an infinite algebraic system of simultaneous equations, where the particular case $h\to \infty$ is verified. The results of this system are also obtained numerically. The normal displacement, normal stress, and the stress singularity factor are given analytically and shown graphically with discussion. By comparison with those predicted by the finite element method, the accuracy of the numerical method is approved.

Introduction

The problems of the theory of elasticity are more significant in mechanics in particular those concerning contact. Indeed, the concentration of the efforts and the tiredness of the material often occurs in the zone of contact between structural components. This is why the zone of contact is classified critical zone. Due to its potential applicability to several constructions of practical importance, the contact problems in solid mechanics involving an elastic layer sitting on an elastic half-plane, elastic and rigid foundations have been extensively researched. Foundation grillages, pavements of highways and airfields, railways ballast, and ball and roller bearings are some of the application areas of contact mechanics. The classical theory of contact mechanics started with the work of Hertz. ¹ The axisymmetric problem of the surface indentation of an isotropic elastic half-space by a smooth rigid circular punch or indentor was first examined by Boussinesq. ² Lebedev and Ufliand ³ considered the punch problem for an elastic layer resting on a rigid foundation. The contact problem of the elastic medium substrate body has been investigated by Dhaliwal. ⁴ Hayes et al. ⁵ obtained a theoretical solution to the punch problem of axisymmetric indentation of an infinite elastic layer bonded to a rigid half-space as a model for the layered geometry of cartilage and subchondral bone. An axisymmetric contact problem for an elastic medium on a rigid foundation with a cylindrical hole has been studied by Dhaliwal and Singh. ⁶ The contact problem of a half-space pressed onto a rigid base with a cylindrical protrusion or pit has been treated by Shibuya et al. ⁷ Their second paper ⁸ deals with the contact stresses of an elastic plate pressed onto a rigid base with a cylindrical protrusion or pit. In the work, the layer thickness on a rigid base with a cylindrical protrusion or pit and is indented by a rigid stamp has been analyzed by Hara et al. ⁹ Kebli et al. ¹⁰ solved the problem of elastostatic deformation of the elastic medium. The thermoelastic contact problem has been studied by Guerrache and Kebli. ¹¹ The solution for an axisymmetric torsion problem was presented by Kebli et al. ¹²

Various programs based on the finite element method FEM have been created in parallel with developments in computer technology. The FEM is applied to the calculation of contact stresses between two or elastic bodies with friction on the contact surface by Ohte. ¹³ In the paper, the FEM for solving a contact problem of elasticity theory has been examined by Barlam. ¹⁴

The half-space is compressed by a rigid cylindrical surface has been examined by Komvopoulos, ¹⁵ the contact stresses, and displacements in the contact area are determined using a numerical method. An analytical method and a FEM are used to solve a receding contact problem for two elastic layers supported with a Winkler foundation by Oner et al. ¹⁶ Yaylaci ¹⁷ discussed a comparison between numerical and analytical solutions for the receding contact problem. Their second work ¹⁸ deals with the contact problem of an elastic layer resting on a rigid foundation, resolved by FEM.

The frictionless contact problem for an elastic medium on a rigid circular base can be reduced to a system of dual integral equations by using a Hankel transformation, and the Boussinesq stress functions, which is transformed to solving an infinite system of simultaneous algebraic equations through the Gegenbauer’s formula by expressing the normal stress under the indenter as an appropriate Chebyshev series. It seems there is a similar method to that in works^8,10 that is convenient to obtain the numerical results. The results analytical are presented with a discussion to illustrate the validity and accuracy of the proposed analytical method with this problem given by a half-space case treated by Shibuya et al. ⁷ An FEM based approach to this problem is also developed. The results have shown a graphical form that the FEM is in good agreement with the analytical method.

Governing equation of the problem and its solution

Problem description

Figure 1 shows the geometry and coordinate system studied. It uses a cylindrical coordinate system $(r,\theta ,z)$ , $u$ and $w$ are represented the components of displacement in the plane $r$ and $z$ , respectively. Shear modulus $G$ and Poisson’s ratio $\nu$ define the material’s characteristics. The stress tensor’s elements are denoted as ${\sigma }_r,{\sigma }_z,{\sigma }_{\theta },{\tau }_{\mathrm{rz}},{\tau }_{r\theta }$ , and ${\tau }_{z\theta }$ . $h$ is the thickness of the isotropic elastic medium, partially reposing on a smooth rigid circular base of radius $a$ . As we see, the layer tends to be infinite. An indented displacement $\delta$ is pressed on the plane $z=h$ using a rigid punch. It is assumed to be small. Displacement and stress fields must check the boundary conditions

$$w(r,0)=0,r\le a,$$

(1a)

$${\sigma }_z(r,0)=0,r>a,$$

(1b)

$${\tau }_{\mathrm{rz}}(r,0)=0,r\ge 0,$$

(1c)

$$w(r,h)=-\begin{array}{c}\delta \hfill \end{array},r\ge 0,$$

(1d)

$${\tau }_{\mathrm{rz}}(r,h)=0,r\ge 0.$$

(1e)

Figure 1.

Coordinates and configuration.

All stress components vanish at infinity.

Analytical solutions

To satisfy the field equations of the linear theory of elasticity and in the absence of the body force fields, the displacements and stresses in a medium without torsion can conveniently be expressed in terms of two harmonic stress Boussinesq ${\phi }_0(r,z)$ and ${\phi }_3(r,z)$ through the following relations ⁸

$$2\mathrm{Gu}(r,z)=\frac{\partial {\phi }_0(r,z)}{\partial r}+z\frac{\partial {\phi }_3(r,z)}{\partial r},$$

(2)

$$2\mathrm{Gw}(r,z)=\frac{\partial {\phi }_0(r,z)}{\partial z}+z\frac{\partial {\phi }_3(r,z)}{\partial z}-\chi {\phi }_3(r,z),$$

(3)

$$v(r,z)=0,$$

(4)

$${\sigma }_r(r,z)=\frac{{\partial }^2{\phi }_0(r,z)}{\partial r^2}+z\frac{{\partial }^2{\phi }_3(r,z)}{\partial r^2}-2\nu \frac{\partial {\phi }_3(r,z)}{\partial z},$$

(5)

$${\sigma }_{\theta }(r,z)=\frac{\partial {\phi }_0(r,z)}{r\partial r}+z\frac{\partial {\phi }_3(r,z)}{r\partial r}-2\nu \frac{\partial {\phi }_3(r,z)}{\partial z},$$

(6)

$${\sigma }_z(r,z)=\frac{{\partial }^2{\phi }_0(r,z)}{\partial z^2}+z\frac{{\partial }^2{\phi }_3(r,z)}{\partial z^2}-2(1-\nu )\frac{\partial {\phi }_3(r,z)}{\partial z},$$

(7)

$${\tau }_{\mathrm{rz}}(r,z)=\frac{\partial }{\partial r}[\frac{\partial {\phi }_0(r,z)}{\partial z}+z\frac{\partial {\phi }_3(r,z)}{\partial z}-(1-2\nu ){\phi }_3(r,z)],$$

(8)

$${\tau }_{r\theta }(r,z)={\tau }_{\theta z}(r,z)=0.$$

(9)

$\chi =3-4\nu$ denote the Kosolov’s constant. The stress functions ${\phi }_0(r,z)$ and ${\phi }_3(r,z)$ meets the following equation

$$\{\begin{array}{c}{\nabla }^2{\phi }_0(r,z)={\nabla }^2{\phi }_3(r,z)=0,\hfill \\ \multicolumn{1}{c}{{\nabla }^2\equiv \frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}+\frac{{\partial }^2}{\partial z^2}.}\end{array}$$

(10)

${\nabla }^2$ is Laplace’s operator. To solve the problem, it is useful to use a Hankel integral transform to develop equation (3), which is described as follows ¹⁹

$$H_n[f\left(\lambda \right),r]={\int }_0^{\infty }\lambda \begin{array}{c}f\left(\lambda \right)J_n\left(\lambda r\right)d\lambda \hfill \end{array}.$$

(11)

$H_n$ is the Hankel operator of order n, $J_n$ is the first kind of order n of the Bessel function. The integral representation ${\phi }_0(r,z)$ and ${\phi }_3(r,z)$ can be chosen to satisfy equation (3) ⁹

$$\begin{array}{c}{\phi }_0(r,z)=-\begin{array}{c}\frac{\nu \begin{array}{c}G\hfill \end{array}}{h}\hfill \end{array}\delta (2z-r^2)\hfill \\ +{\int }_0^{\infty }[A\left(\lambda \right)\cos h\lambda z+B\left(\lambda \right)\sin h\lambda z]J_0\left(\lambda r\right)d\lambda ,\hfill \end{array}$$

(12)

$${\phi }_3(r,z)=\frac{G}{h}\delta z+{\int }_0^{\infty }[C\left(\lambda \right)\sin h\lambda z+D\left(\lambda \right)\cos h\lambda z]J_0\left(\lambda r\right)d\lambda .$$

(13)

$A$ , $B$ , $C$ , and $D$ are arbitrary functions of $\lambda$ . Substituting equations (12) and (13) into equations (2) to (9), using boundary conditions (1c), (1d), and (1e), allows us to derive the expression

$$\lambda A\left(\lambda \right)=\frac{\lambda \begin{array}{c}h\hfill \end{array}-(1-2\nu )\sin h\lambda h\cos h\lambda h}{\sin h^2\lambda h}D\left(\lambda \right),$$

(14)

$$\lambda B\left(\lambda \right)=(1-2\nu )\begin{array}{c}D\left(\lambda \right)\hfill \end{array},$$

(15)

$$C\left(\lambda \right)=-\cot h\lambda h\begin{array}{c}D\left(\lambda \right)\hfill \end{array}.$$

(16)

The $w$ and ${\sigma }_z$ of equations (2) to (9) can be chosen using functions Boussinesq ${\phi }_0(r,z)$ and ${\phi }_3(r,z)$ of equations (2) to (9). By the boundary conditions (1a) and (1b), we find that the function D(λ) must meet the system of dual integral equations in terms of the Hankel transform as follows

$$w(r,0)=\eta {\int }_0^{\infty }D\left(\lambda \right)J_0\left(\lambda r\right)d\lambda =0,r\le a,$$

(17)

$${\sigma }_z(r,0)+\frac{E}{h}\delta =0\Rightarrow {\int }_0^{\infty }\lambda \begin{array}{c}q\left(\lambda \right)D\left(\lambda \right)J_0\left(\lambda r\right)d\lambda \hfill \end{array}=0,r>a.$$

(18)

$E$ is defined by Young’s modulus

$$\eta =\frac{-1+\nu }{G},$$

(19)

$$q\left(\lambda \right)=\frac{\lambda h+\sin h\lambda h\cos h\lambda h}{\sin h^2\lambda h}.$$

(20)

Assuming ²⁰ (formula 6.522.2), the integral formula

$${\int }_0^{\infty }\begin{array}{c}\lambda \hfill \end{array}{M}_n\left(\begin{array}{c}\lambda \hfill \end{array}x\right)J_0\left(\begin{array}{c}\lambda \hfill \end{array}r\right)\begin{array}{c}d\lambda \hfill \end{array}=\{\begin{array}{cc}\frac{2}{\pi r}\frac{T_{2n+1}(r/x)}{\sqrt{x^2-r^2}},\hfill & r<x,\hfill \\ \multicolumn{1}{c}{0,}& r>x.\hfill \end{array}$$

(21)

$T_n$ is the first kind of order n of Chebyshev’s polynomials. Here, putting

$$M_n\left(\begin{array}{c}\lambda \hfill \end{array}x\right)=J_{n+\frac{1}{2}}\left(\frac{\begin{array}{c}\lambda \hfill \end{array}x}{2}\right)J_{-(n+\frac{1}{2})}\left(\frac{\begin{array}{c}\lambda \hfill \end{array}x}{2}\right),\begin{array}{c}(n=0,1,2...)\hfill \end{array}.$$

(22)

Permits us to solve the dual integral equations (17) and (18), using the following expression

$$q\left(\lambda \right)D\left(\lambda \right)=\frac{\delta }{\eta }{\sum }_{n=0}^{\infty }{\alpha }_n{M}_n\left(\begin{array}{c}\lambda \hfill \end{array}a\right).$$

(23)

Substituting from equation (23) into equation (17), and making Gegenbauer’s formula ²⁰ (formula 8.531.1), we get

$$\begin{array}{c}\multicolumn{1}{c}{J_0\left(\begin{array}{c}\lambda \hfill \end{array}r\right)={\sum }_{m=0}^{\infty }(2-{\delta }_{0m})X_m\left(\begin{array}{c}\lambda \hfill \end{array}x\right)\cos \left(m\varphi \right),(\varphi =\arcsin \left(\frac{r}{x}\right)).}\end{array}$$

(24)

${\delta }_{0m}$ being the Kronecker delta function ${\delta }_{\mathrm{nm}}=\{\begin{array}{cc}1,\hfill & m=n,\hfill \\ \multicolumn{1}{c}{0,}& m\ne n.\hfill \end{array}$

$$X_m\left(\begin{array}{c}\lambda \hfill \end{array}x\right)=J_m^2\left(\frac{\begin{array}{c}\lambda \hfill \end{array}x}{2}\right),\begin{array}{c}(m=0,1,2...)\hfill \end{array}.$$

(25)

we obtained

$$\begin{array}{c}{\sum }_{n=0}^{\infty }{a}_n{\sum }_{n=0}^{\infty }(2-{\delta }_{0m})\cos m\varphi \hfill \\ {\int }_0^{\infty }{q}^{-1}\left(\lambda \right)M_n\left(\begin{array}{c}\lambda \hfill \end{array}a\right)X_m\left(\begin{array}{c}\lambda \hfill \end{array}a\right)d\lambda =0,r\le a.\hfill \end{array}$$

(26)

where

$$a_n=\delta \begin{array}{c}{\alpha }_n\hfill \end{array},\begin{array}{c}(n=0,1,2...)\hfill \end{array}.$$

(27)

Equation (26) must hold for an arbitrary value of $\varphi$ . The $a_n$ coefficients are carried from equation (28), as a result of matching the coefficients on both sides of equation (26).

$${\sum }_{n=0}^{\infty }{a}_n{\int }_0^{\infty }{q}^{-1}\left(\lambda \right)M_n\left(\begin{array}{c}\lambda \hfill \end{array}a\right)X_m\left(\begin{array}{c}\lambda \hfill \end{array}a\right)d\lambda =0,\begin{array}{c}(m=0,1,2...)\hfill \end{array}.$$

(28)

In matrix form, equation (28) becomes

$${\sum }_{n=0}^{\infty }{a}_n{A}_{\mathrm{mn}}=0,\begin{array}{c}(m=0,1,2...)\hfill \end{array}.$$

(29)

where

$$A_{\mathrm{mn}}={\int }_0^{\infty }{q}^{-1}\left(\lambda \right)M_n\left(\begin{array}{c}\lambda \hfill \end{array}a\right)X_m\left(\begin{array}{c}\lambda \hfill \end{array}a\right)d\lambda .$$

(30)

We can derive the results for the particular case. The result obtained in the problem ⁷ is easily obtained by taking $h\to \infty$ , the infinite integrals of the system equation (29), we find that equation

$${\sum }_{n=0}^{\infty }{a}_n{\int }_0^{\infty }{M}_n\left(\lambda a\right)\begin{array}{c}{X}_m\left(\lambda a\right)d\lambda \hfill \end{array}=0,\begin{array}{c}(m=0,1,2...)\hfill \end{array}.$$

(31)

We write the system equation (29) in dimensionless form, which allows us to facilitate their simplification and to generalize the results. By setting the following changes of variables

$$t=\lambda \begin{array}{c}a\hfill \end{array},H=\frac{h}{a}.$$

(32)

Then

$${\sum }_{n=0}^{\infty }{a}_n\prime \begin{array}{c}{A}_{\mathrm{mn}}\prime =0\hfill \end{array}.$$

(33)

where

$$a_n\prime =a_n,$$

(34)

$$A_{\mathrm{mn}}\prime ={\int }_0^{\infty }{q}^{-1}\left(t\right)M_n\left(t\right)X_m\left(t\right)\mathrm{dt},$$

(35)

and

$$q^{-1}\left(t\right)=\frac{\sin h^2\mathrm{Ht}}{\mathrm{Ht}+\cos \mathrm{hHt}\sin \mathrm{hHt}},$$

(36)

$$M_n\left(t\right)=J_{n+\frac{1}{2}}\left(\frac{t}{2}\right)J_{-(n+\frac{1}{2})}\left(\frac{t}{2}\right),$$

(37)

$$X_m\left(t\right)=J_m^2\left(\frac{t}{2}\right).$$

(38)

Numerical results

To estimate the unknown coefficients $a_n\prime$ , it is necessary to solve the infinite set of simultaneous equation (33), which involves the multiplication of Bessel functions in the integrand. Therefore, we follow the steps described in the work of Kebli et al. ¹⁰ The numerical values of coefficients $a_n\prime$ with values for the various values of the parameter $H$ are presented in the following Table 1. A good convergent solution can be achieved using eight terms.

Table 1.

Coefficient values $a_n\prime$ .

$a_n\prime$
n	$H=0.75$	$H=1$	$H=1.5$
0	0.000290751991865906	0.001088955486405	0.00106534756817198
1	−0.000360749489995402	−0.001526440370283	−0.00117263625004920
2	0.000104751189188284	0.000579796833755	0.00024352913575814
3	−0.000019911540295819	−0.000117589292941	−0.00005248678208980
4	−0.000000681653785655	0.000007418581197	−0.00000957747694767
5	−0.000001436598849622	−0.000003519173837	−0.00000759248737481
6	−0.000000862916750635	−0.000001808580884	−0.00000471927687651
7	−0.000000514438425896	−0.000001082331256	−0.00000280506434363
8	−0.000000208367030836	−0.000000410528142	−0.00000117807920084
9	−0.000000012790110256	0.000000000819269	−0.00000011117948242

Expressions for the physical quantities

The nondimensional normal displacement in the plane $z=0$ can be obtained by solving the following equation

$$\begin{array}{c}{w}^{*}(r,0)=\frac{w(r,0)}{\delta }={\sum }_{n=0}^{\infty }{a}_n\{{\int }_0^{\infty }[q^{-1}\left(\lambda \right)-1]M_n\left(\lambda \begin{array}{c}a\hfill \end{array}\right)J_0\left(\lambda r\right)d\lambda \hfill \\ +{\int }_0^{\infty }[M_n\left(\lambda \begin{array}{c}a\hfill \end{array}\right)-\frac{2}{\lambda \pi a}\sin \lambda a]J_0\left(\lambda r\right)\mathrm{dr}\hfill \\ +\frac{2}{\pi a}[H(a-r)\frac{\pi }{2}+H(r-a){\sin }^{-1}\frac{a}{r}]\}.\hfill \end{array}$$

(39)

where

$H$ denotes Heaviside’s step function $H(x-r)=\{\begin{array}{cc}1,\hfill & r<x,\hfill \\ \multicolumn{1}{c}{0,}& r>x.\hfill \end{array}$

Therefore, the nondimensional normal stress in the contact region between the elastic medium and the rigid base can be expressed in appropriate Chebyshev series, with unknown coefficients $a_n$ as follows as

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }^{*}(r,0)=\frac{h}{E\delta }[{\sigma }_z(r,0)+1]}\\ =\frac{2}{\pi }H(a-r){\sum }_{n=0}^{\infty }{a}_n\frac{T_{2n+1}(r/a)}{r\sqrt{a^2-r^2}}.\hfill \end{array}$$

(40)

The stress state around the plane of contact base can be described by the stress singularity factor defined as ¹⁰

$$S=\underset{r\to a^{-}}{\lim }\sqrt{2\pi (r-a)}\begin{array}{c}{\sigma }^{*}(r,0)\hfill \end{array}.$$

(41)

The substitution of equation (40) into equation (41), we obtained

$$S=\frac{2}{a\sqrt{\pi }}{\sum }_{n=0}^{\infty }{a}_n.$$

(42)

The behavior of nondimensional normal displacement $w^{*}(\rho ,0)$ where $\rho =\frac{r}{a}$ , on the region $\rho \ge 1$ for different values $H$ is given in Figure 2. We remarked that the value $w^{*}(\rho ,0)$ is decreasing with an increase $H$ . The nondimensional normal stress ${\sigma }^{*}(\rho ,0)$ can be seen in Figure 3. It is shown, the values ${\sigma }^{*}(\rho ,0)$ get their maximum at the center of the rigid base with increased values for the parameter $H$ . The slope ${\sigma }^{*}(\rho ,0)$ tends to infinity at $\rho =1$ . The nondimensional stress singularity factor $S$ around the region of contact is graphically illustrated in Figure 4. From the graph lines, it is clear that the $S$ gives a large value with decreasing layer thickness. It is important to note that, the layer thickness and the base radius affect the distribution of $S$ values

Figure 2.

Variation of $w^{*}(\rho ,0)$ with different values of $H$ .

Figure 3.

Variation of ${\sigma }^{*}(\rho ,0)$ with different values of $H$ .

Figure 4.

Variation of $S$ with different values of $H$ .

Finite element simulation modeling

In the modeling carried out, an axisymmetric finite element method Figure 5 is constructed based on the joint model of Figure 1. The method numerical solution requires as an input some material and geometrical properties. The layer is assumed isotropic and elastic, geometric are taken as $L=2\begin{array}{c}500\hfill \end{array}\mathrm{mm}$ in $x$ the direction, $h=20\begin{array}{c}\mathrm{mm}\hfill \end{array}$ in $y$ direction, and the radius of the base as $a=2\begin{array}{c}5\hfill \end{array}\mathrm{mm}$ . It should be noted that the length of the elastic medium is very large than the radius of the base. In the analysis, material properties are taken as Young’s modulus of elasticity $\begin{array}{c}21\times {10}^4\hfill \end{array}\mathrm{Mpa}$ , $0.29$ which is the value of Poisson’s ratio. The boundary conditions applied to correspond to the blockages in translation along $x$ and $y$ of the rigid base, and to the imposed displacement $\delta =\begin{array}{c}\begin{array}{c}\begin{array}{c}-\hfill \end{array}\hfill \end{array}0.05\hfill \end{array}\mathrm{mm}$ of the free vertical surface of the layer in the $y$ direction, Figure 5(a). The FEM mesh is based on the $\mathrm{CPS}4R$ type, and of size $(0.2\begin{array}{c}5\hfill \end{array}\times 0.2\begin{array}{c}5\hfill \end{array})\mathrm{mm}$ , Figure 5(b). The refined mesh is about of $800000$ elements, and the time step is $\begin{array}{c}0.05\hfill \end{array}$  s. Our convergence analysis was conducted on the maximum normal stress and the vertical displacement. We based ourselves on the calculation of the relative error of these two quantities by using an increasingly refined mesh. The value of the relative error is obtained compared to the value of normal stress, and the vertical displacement is calculated with the preceding mesh, coarser. In both cases, the relative error decreases at the end of the refinement of the model and remains below $0,2$ % for the vertical displacement and $2$ % for the case of vertical stress, which affirms the convergence of the two quantities.

Figure 5.

Boundary condition applied (a) and mesh of model (b).

The contour of normal stress distribution on the layer thickness obtained from the numerical method is plotted to the $y$ -axis in Figures 6 and 7. It is obvious from Figures 6 and 7 that the normal stress attains infinite value at the edges of the rigid base. It is noted that the results obtained are a good agreement with the analytical results. Figures 8 and 9, show the contour plot of Von Mises stress distribution. We observe a very high concentration of stresses around the edges of the rigid base for the Von Mises stress and it is canceled while moving away from the base. The simulation results obtained of vertical displacement are given in Figures 10 and 11. The latter represents the displacement behavior through the elastic medium. The value of vertical displacement is maximum in the vicinity of the rigid base, then it decreases while moving away from the vertical axis, and that is due to the displacement applied. We found good agreement between previously obtained vertical displacement analysis results.

Figure 6.

Contour of normal stress distribution on the layer thickness.

Figure 7.

Normal stress of distribution along the rigid base surface.

Figure 8.

Contour of Von Mises stress distribution of a layer thickness.

Figure 9.

Von Mises stress distribution at the rigid base surface.

Figure 10.

Contour of vertical displacement distribution on the layer thickness.

Figure 11.

Distribution of vertical displacement in the plane $z=0$ .

Conclusion

In the present work, the frictionless axisymmetric contact problem for an elastic layer thickness partially reposing on a rigid base is considered. The problem is solved by using a boundary value, the theory of elasticity, and the integral transformation technique. The finite element method FEM is performed for this problem. Following is a summary of some findings provided from mathematical formulations and the results:

• − Using the Boussinesq stress functions and the Hankel integral transformations, the studied axisymmetric frictionless contact problem of an elastic medium is reduced to a system of dual integral equations.

• − Instead of the traditional Fredholm equation approach, the latter equations are converted into an infinite algebraic system of simultaneous equations. Only the problem’s geometrical parameters are dependent on this system. This allows numerical results of the infinite algebraic system to be obtained for various values of the radius of the rigid base and the thickness layer.

• − The thickness layer and the base radius effect on the distribution of normal displacement, normal stress, and stress singularity factor are clarified by numerical computations and graphical results.

• − The methods of analysis and the result of this work are consistent with the result obtained by Shibuya et al. ⁷

• − The computational approach based on the finite element method is validated by the analysis described.