A Simple Semi-Analytical Method for Solving Axisymmetric Contact Problems Involving Bonded and Unbonded Layers of Arbitrary Thickness

Abstract

In the present work, a recently extended version of the method of dimensionality reduction (MDR) for layered elastic media is applied for the first time using a semi-analytical approach. It is based on a priori knowledge of the cylindrical flat punch solution which is determined numerically using the boundary element method (BEM). We consider arbitrary indenters of revolution producing a circular area of contact with bonded and unbonded layers of arbitrary thickness. The proposed method reduces the contact solution to the numerically efficient evaluation of simple one-dimensional integrals. We further show that the solution of JKR-adhesive contacts with layers and contacts with linear-viscoelastic layers is straightforward using the well-known MDR formalisms. A specific focus has been devoted to study the thickness effect in different application examples. Comparisons with the literature and finite element simulations show very good agreement with the proposed method.

1. Introduction

Contact solutions based on the half-space assumption, i.e., requiring the contact area to be much smaller than the dimensions of the contacting bodies, are usually simple and often available in closed-form analytical expressions. Therefore, even for contact problems with significant influence of finite geometry, this requirement is often disregarded and the half-space-based contact solutions are used as rough estimates. However, if the body dimensions are of the order of the contact dimensions or smaller, finite geometry must necessarily be considered. For very thin layers, again simple solutions are found within the thin layer assumption [1,2]. Unfortunately, for many contact problems neither asymptotic solution is applicable. A prominent example are layers of articular cartilage in biological joints. The cartilage layers are bonded to significantly stiffer bone material and the contact radius is usually of the order of the thickness or greater [3]. It has been shown that small variations in the thickness of cartilage layers have a pronounced effect on their viscoelastic response [4,5]. For adhesive contacts, the layer thickness directly influences the adhesive strength [6], but can also cause adhesive instabilities [7].

In fully computational methods such as the finite element method, finite geometry is readily considered. However, the calculation times are too high for comprehensive parameter studies and common effects such as adhesion are not implemented natively. Li et al. [8] developed a boundary element method (BEM) for layered half-spaces which is advantageous in terms of numerical complexity and allows for a simple implementation of JKR-type adhesion [9]. Further simplification and reduction in computational cost is possible if only axisymmetric contacts are considered. Then, the solution of the contact problem of a rigid indenter and an elastic layer can be reduced to the solution of one-dimensional Fredholm integral Equations [10,11,12,13]. Recently, Argatov et al. [14] proposed an even simpler approach; they extended the well-known method of dimensionality reduction (MDR) to layered elastic media. Within the MDR, which was introduced by Popov and Heß [15], a very simple equivalent one-dimensional contact problem of a rigid plane profile and a one-dimensional Winkler foundation is solved. It exactly reproduces the solution of three-dimensional contact problems including adhesion [16], viscoelastic half-spaces [17] and functionally graded half-spaces [18]. The proposed MDR extension to layered media is based on a priori knowledge of the contact solution of a rigid cylindrical flat punch with an elastic layer for the whole parameter space. If this flat punch solution is available, the elastic solution for arbitrary axisymmetric profiles is simple and the consideration of adhesion and viscoelasticity follows the well-known rules of the MDR framework.

In the present paper, the MDR for layered media as proposed by Argatov et al. [14] is employed for the first time for a single layer, which is either bonded to a rigid substrate or slides frictionless on the substrate (unbonded case). Using a semi-analytical approach, the required flat punch solution is provided in numerical form, determined by a BEM framework. By investigating the thickness effect in different application examples including adhesive layers and viscoelastic layers, the efficiency and variability of the proposed method are shown. The obtained numerical BEM data needed to solve the MDR integrals are provided in the Supplementary Materials.

2. MDR Framework

We consider the frictionless axisymmetric indentation of a linear elastic layer of thickness h with a rigid indenter with profile $f\left(r\right)$ and a compact contact area as shown in Figure 1A. The layer rests on a rigid substrate and is either completely bonded to the substrate or slides frictionless on the substrate (unbonded case). The surface displacement within the contact area at given indentation depth $\delta$ is given by

$$w\left(r\right)=\delta -f\left(r\right).$$

(1)

In the framework of the method of dimensionality reduction (MDR), the much simpler equivalent contact problem of a one-dimensional bedding of independent springs and a one-dimensional indenter with profile $g(x;h)$ shown in Figure 1B is used to find the complete contact solution. The one-dimensional spring displacement within the contact area is

$$w_{1\mathrm{D}}(x;h)=\delta -g(x;h)$$

(2)

and, thus, if $g(x;h)$ is continuous at $x=a$, the MDR profile of the nonadhesive contact is simply the relation of indentation depth and contact radius:

$$\delta =g(a;h).$$

(3)

The normal force is found by integrating the spring forces,

$$F_{\mathrm{N}}(a;h)=2\underset{0}{\overset{a}{\int }}{c}_{\mathrm{N}}(x;h)w_{1\mathrm{D}}(x;h)\mathrm{d}x,$$

(4)

with the spring stiffness $c_{\mathrm{N}}$.

By taking into account the universality of the contact stiffness,

$$k_{\mathrm{N}}:=\frac{\mathrm{d}{F}_{\mathrm{N}}}{\mathrm{d}\delta }=\frac{\mathrm{d}{F}_{\mathrm{N}}\left(a\right)}{\mathrm{d}a}\frac{\mathrm{d}a\left(\delta \right)}{\mathrm{d}\delta },$$

(5)

the relation of spring stiffness and contact stiffness can be found from Equations (3) and (4),

$$c_{\mathrm{N}}(x;h)=\frac{1}{2}{k}_{\mathrm{N}}^{\prime }(x;h).$$

(6)

In contrast to the well-known half-space MDR, the spring stiffness $c_{\mathrm{N}}(x;h)$ for layered bodies is generally variable. According to Betti’s reciprocal theorem, the normal force can also be written using the pressure distribution resulting from the unit indentation by a rigid cylindrical flat punch $p^{*}(r,a/h)$,

$$F_{\mathrm{N}}=2\pi {\int }_0^{a}w\left(r\right)p^{*}(r,a/h)r\mathrm{d}r.$$

(7)

By equating Equations (4) and (7) and differentiating with respect to a, the transformation that gives the MDR profile $g\left(x\right)$ based on a given profile $f\left(r\right)$ is found:

$$g(x;h)=\frac{\pi }{c_{\mathrm{N}}(x;h)}\frac{\mathrm{d}}{\mathrm{d}x}\left\{\underset{0}{\overset{x}{\int }}{p}^{*}(r;x/h)f\left(r\right)r\mathrm{d}r\right\}.$$

(8)

In the context of the MDR, the transformation in Equation (8) was formulated by Argatov et al. [14] based on the previous work by Efimov et al. [19]. It is well known that the contact pressure can be obtained from the superposition of incremental flat punches indentations. Using the previous definitions, it can be obtained from the one-dimensional spring displacement,

$$p(r;a/h)=-\underset{r}{\overset{a}{\int }}{p}^{*}(r;x/h)w_{1\mathrm{D}}^{\prime }(x;h)\mathrm{d}x+p^{*}(r;a/h)w_{1\mathrm{D}}(a;h),$$

(9)

where the second part vanishes if the one-dimensional profile is continuous at $x=a$. Let us now introduce the following dimensionless variables that are more convenient for the numerical evaluation:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \xi =\frac{x}{h},{\xi }_a=\frac{a}{h},\alpha =\frac{r}{h},\overline{f}\left(\alpha \right)=\frac{f\left(r\right)}{h},{\overline{w}}_{1\mathrm{D}}\left(\xi \right)=\frac{w_{1\mathrm{D}}(x;h)}{h},{\overline{k}}_{\mathrm{N}}\left(\xi \right)=\frac{k_{\mathrm{N}}(x;h)}{2E^{*}x},\hfill \\ \hfill & {\overline{c}}_{\mathrm{N}}\left(\xi \right)=\frac{c_{\mathrm{N}}(x;h)}{E^{*}}={\overline{k}}_{\mathrm{N}}\left(\xi \right)+\xi {\overline{k}}_{\mathrm{N}}^{\prime }\left(\xi \right),\overline{p}(\alpha /\xi ;\xi )=\frac{\pi x}{2E^{*}{\overline{k}}_{\mathrm{N}}(x;h)}{p}^{*}(r/x;x/h),\hfill \end{array}\end{array}$$

(10)

with the elastic modulus $E^{*}=E/(1-{\nu }^2)$. Using these normalizations, the transformations (8) and (9) are rewritten as

$$g\left(\xi \right)=\frac{2h}{{\overline{c}}_{\mathrm{N}}\left(\xi \right)}\frac{\mathrm{d}}{\mathrm{d}\xi }\left\{\frac{{\overline{k}}_{\mathrm{N}}\left(\xi \right)}{\xi }\underset{0}{\overset{\xi }{\int }}\overline{p}(\alpha /\xi ;\xi )\overline{f}\left(\alpha \right)\alpha \mathrm{d}\alpha \right\},$$

(11)

and

$$p(\alpha ;{\xi }_a)=\frac{2E^{*}}{\pi }\left\{-\underset{\alpha }{\overset{{\xi }_a}{\int }}\frac{{\overline{k}}_{\mathrm{N}}\left(\xi \right)}{\xi }\overline{p}(\alpha /\xi ;\xi ){\overline{w}}_{1\mathrm{D}}^{\prime }\left(\xi \right)\mathrm{d}\xi +\frac{{\overline{k}}_{\mathrm{N}}\left({\xi }_a\right)}{{\xi }_a}\overline{p}(\alpha /{\xi }_a;{\xi }_a){\overline{w}}_{1\mathrm{D}}\left({\xi }_a\right)\right\},$$

(12)

respectively. For the solution of the contact problem using Equations (4), (11) and (12), only the normalized contact stiffness ${\overline{k}}_{\mathrm{N}}(a/h)$ and the normalized pressure distributions under the cylindrical flat punch $\overline{p}(r/a;a/h)$ are needed. It should be noted that for a layer that is bonded to the rigid substrate, these normalized quantities also depend on Poisson’s ratio $\nu$. We omitted this dependency here for brevity.

For the sake of completeness, it should be mentioned that the displacement outside of the contact zone $w(r,a/h)$ can be determined using a similar superposition,

$$w(r,a/h)=-\underset{0}{\overset{a}{\int }}{w}^{*}(r;x/h)w_{1\mathrm{D}}^{\prime }(x;h)\mathrm{d}x+w^{*}(r;a/h)w_{1\mathrm{D}}(a;h),$$

(13)

where $w^{*}(r,a/h)$ is the displacement outside of the contact zone for an unit indentation with a cylindrical flat punch.

2.1. Asymptotic Cases

In this section, the asymptotic cases of very thick and very thin layers are discussed. Within these assumptions, the given MDR relations may be reduced to very simple expressions that allow for the derivation of closed-form solutions for many axisymmetric indenter profiles $f\left(r\right)$. The derivations of the contact stiffness and normalized flat punch pressure distribution are, for example, described by Barber [20] in more detail.

2.1.1. The Elastic Half-Space

For the asymptotic case of a thick layer with $a/h\to 0$, the solution is independent of the boundary condition between layer and rigid substrate. With the contact stiffness of the elastic half-space $k_{\mathrm{N}}=2E^{*}a$ and the normalized pressure distribution under a flat punch

$$\overline{p}(r/a)=\frac{1}{2}{\left(1-{\left(\frac{r}{a}\right)}^2\right)}^{-1/2},$$

(14)

the bedding stiffness is constant $c_{\mathrm{N}}=E^{*}$ and the well-known MDR relations are recovered from (11) and (12):

$$g\left(x\right)=x\underset{0}{\overset{x}{\int }}\frac{f^{\prime }\left(r\right)\mathrm{d}r}{\sqrt{x^2-r^2}},p(r;a)=\frac{E^{*}}{\pi }\left\{-\underset{r}{\overset{x}{\int }}\frac{w_{1\mathrm{D}}^{\prime }\left(x\right)\mathrm{d}x}{\sqrt{x^2-r^2}}+\frac{w_{1\mathrm{D}}\left(a\right)}{\sqrt{a^2-r^2}}\right\}.$$

(15)

2.1.2. The Compressible Bonded Thin Layer

In the case of a bonded thin layer with $\nu <0.5$ and $a\gg h$, the in-plane strains will be negligible. Hence, the layer itself acts as a three-dimensional linear Winkler foundation. Then, the relation of indentation depth and contact radius a is given by the profile, $\delta =f\left(a\right)$, the contact pressure under a flat punch is constant and a simple analysis of the elasticity equations yields the contact stiffness,

$$k_{\mathrm{N}}=\pi \frac{E(1-\nu )}{(1+\nu )(1-2\nu )}\frac{a^2}{h}\Rightarrow {\overline{k}}_{\mathrm{N}}=\frac{\pi }{2}\frac{{(1-\nu )}^2}{(1-2\nu )}\frac{a}{h}.$$

(16)

Thus, the MDR relations for this case are very simple:

$${\overline{c}}_{\mathrm{N}}=\pi \frac{{(1-\nu )}^2}{(1-2\nu )}\frac{x}{h},g\left(x\right)=f\left(x\right),p\left(r\right)=\frac{E(1-\nu )}{(1+\nu )(1-2\nu )}\frac{w\left(r\right)}{h}$$

(17)

2.1.3. The Incompressible Bonded Thin Layer

Due to the incompressibility ($\nu =0.5$) and the geometric confinement, the material of an indented thin incompressible layer with $a\gg h$ flows laterally and bulges upwards at the contact edge. Barber [20] used a quadratic approximation for the in-plane displacement to derive the solution under a spherical indenter. The according derivation for the cylindrical flat punch yields the contact stiffness

$$k_{\mathrm{N}}=\frac{\pi E{a}^4}{8h^3}\Rightarrow {\overline{k}}_{\mathrm{N}}=\frac{3\pi }{64}{\left(\frac{a}{h}\right)}^3,$$

(18)

and the contact pressure

$$\overline{p}(r/a)=2\left(1-{\left(\frac{r}{a}\right)}^2\right).$$

(19)

Thus, the resulting MDR relations for this case are

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\overline{c}}_{\mathrm{N}}& =\frac{3\pi }{16}{\left(\frac{x}{h}\right)}^3,g\left(x\right)=\frac{2}{x^2}\underset{0}{\overset{x}{\int }}f\left(r\right)r\mathrm{d}r,\hfill \\ \hfill p(r;a)& =\frac{E}{4h^3}\left\{-\underset{r}{\overset{a}{\int }}(x^2-r^2)w_{1\mathrm{D}}^{\prime }\left(x\right)\mathrm{d}x+(a^2-r^2)w_{1\mathrm{D}}\left(a\right)\right\}.\hfill \end{array}\end{array}$$

(20)

2.1.4. The Unbonded Thin Layer

In the case of an unbonded thin layer with $a\gg h$, the plane section remains as plane and the displacement is just a function of radial coordinate. As for the compressible thin layer, the layer itself acts as a three-dimensional linear Winkler foundation but with a lower stiffness due to the less confined situation,

$$k_{\mathrm{N}}=\pi E^{*}\frac{a^2}{h}\Rightarrow {\overline{k}}_{\mathrm{N}}=\frac{\pi }{2}\frac{a}{h}.$$

(21)

The MDR relations are

$${\overline{c}}_{\mathrm{N}}=\pi \frac{x}{h},g\left(x\right)=f\left(x\right),p\left(r\right)=E^{*}\frac{w\left(r\right)}{h}.$$

(22)

It is worth noting that the above equations also hold for an unbonded incompressible layer, because material can flow laterally.

2.2. Contact of Two Coated Rigid Bodies with Curved Surfaces

The MDR formalism described above is not limited to the case of a rigid indenter in contact with a flat layer, but can also be extended to the contact of two coated rigid bodies with different curvatures $R_1$ and $R_2$ where coating thickness h and material are equal (see Figure 2). Then, an equivalent contact problem is that of a rigid indenter and a flat layer with the same thickness h, where the curvature $R^{*}$ of the indenter is found under the condition that the total displacement in the contact area remains equal,

$$\frac{1}{R^{*}}=\frac{1}{R_1}+\frac{1}{R_2}.$$

(23)

The condition of equal contact stiffness leads to the effective elastic modulus of the flat layer of the equivalent problem,

$$E^{*}=\frac{E}{2(1-{\nu }^2)}.$$

(24)

These relations are, of course, the well-known relations from the Hertzian contact theory of two curved bodies with $a/h\to 0$. Due to the more complicated dependency on Poisson’s ratio for the bonded layer, the equivalent problem is only given for equal materials. If the thickness of the layers differs slightly, a good MDR approximation is found for an equivalent flat layer with thickness

$$h^{*}=\frac{h_1+h_2}{2}.$$

(25)

For compressible bonded thin layers and unbonded thin layers, Equation (25) is "exact" within the thin layer assumption.

3. Numerical Solution of Flat Punch Indentations Using FFT-Based BEM

3.1. FFT-Based BEM for Bonded and Unbonded Elastic Layers

To obtain a numerical flat punch solution needed for the MDR transformations in Section 2, we reformulate the BEM formulation by Li et al. [8] which solves the problem of an elastic layer on an elastic half-space for a single bonded or unbonded layer. In the BEM formalism, the pressure–displacement relation is solved using direct and inverse Fast Fourier Transforms (FFT),

$$\mathit{u}=\mathrm{IFFT}\left[\mathit{K}·\mathrm{FFT}\left(\mathit{p}\right)\right],$$

(26)

with the compliance matrix in the Fourier space $\mathit{K}$, which is the Fourier transform of the discretized fundamental solution: $\mathit{K}=\mathrm{FFT}\left({\mathit{U}}_{\mathit{0}}\right)$. As shown by Li et al., the compliance matrix can be directly found in the Fourier space,

$$\mathit{K}=\frac{2}{k{E}^{*}}\chi \left(kh\right),$$

(27)

with the modulus $E^{*}=E/(1-{\nu }^2)$ and the variable k which is the absolute value of the wave vector $\mathit{k}$. Completely analogous to the work by Li et al., the dimensionless kernel $\chi \left(kh\right)$ for the unbonded layer is obtained by restraining the out-of-plane displacement at the layer–substrate interface and allowing for frictionless in-plane motion,

$$u_z(z=h)=0\mathrm{and}{\tau }_{\{xz,yz\}}(z=h)=0.$$

(28)

The resulting kernel function of the unbonded layer reads:

$$\chi \left(kh\right)=\frac{2{sinh}^2\left(kh\right)}{2kh+sinh\left(2kh\right)}.$$

(29)

For the bonded layer with fixed displacement at the layer–substrate interface,

$$u_{\{x,y,z\}}(z=h)=0,$$

(30)

the kernel has an explicit dependency on Poisson’s ratio,

$$\chi \left(kh\right)=\frac{2\kappa sinh\left(2kh\right)-4kh}{1+4k^2{h}^2+{\kappa }^2+2\kappa cosh\left(2kh\right)},\mathrm{with}\kappa =3-4\nu .$$

(31)

It should be noted that the above kernels have been found by several other authors [21,22]. With the compliances in Equations (29) and (31), the normal stress at the surface due to a sinusoidal surface displacement is given by

$$\sigma (x,y)=\frac{E^{*}k}{2\chi \left(kh\right)}{u}_{0}sin\left(kx\right)$$

(32)

The described BEM framework with the derived kernels allows for an efficient contact solution for arbitrary indenter shapes. This also includes JKR-adhesive problems using the mesh-size dependent criterion introduced by Pohrt and Popov [9].

3.2. Contact Stiffness and Pressure under a Cylindrical Flat Punch

We now consider the indentation of a layer by a cylindrical flat punch using the BEM formalism described in Section 3.1. The resulting contact stiffness (see Equation (6)) normalized with the that of the elastic half-space, $2E^{*}a$, is shown in Figure 3 in terms of the ratio of contact radius to layer thickness $a/h$ for different Poisson’s ratios, while Figure 4 shows the normalized contact pressure under the cylindrical flat punch indented into a bonded, incompressible layer for different $a/h$. In case of the unbonded layer, the normalized stiffness does not depend on Poisson’s ratio, as shown in Section 3.1. We also included the asymptotic thin layer solutions described in Section 2.1.2, Section 2.1.3 and Section 2.1.4 in Figure 3 and Figure 4 and the numerical solution of the governing Fredholm integral equations by Hayes et al. [11] in Figure 3. As expected, the normalized contact stiffness increases with the confinement. The confinement is increased by the geometric ratio of contact radius to layer thickness $a/h$, restricting displacements at the layer–substrate interface (the bonded case) and, for the bonded case, the degree of incompressibility. The normalized pressure distribution exhibits higher values in the middle of the indenter for increasing $a/h$ and lower at the contact edge. In very confined situations, the singularity at the contact edge disappears. The asymptotic behavior for small and large $a/h$ is met in all cases for contact stiffness and the contact pressure and the agreement of with the data by Hayes et al. [11] appears to be very satisfactory. Taking into account the asymptotic solutions summarized in

Table 1. Normalized contact stiffness and relation of indentation depth and contact radius $\delta \left(a\right)$ for a paraboloid with radius of curvature R.

	Norm. Contact Stiffness ${\overline{\mathit{k}}}_{\mathit{N}}$	${\mathit{\delta }}^{\mathbf{P}}\left(\mathit{a}\right)$ for a Paraboloid
Half-space	1	${\displaystyle \frac{a^2}{R}}$
Bonded thin layer ($\nu <0.5$)	${\displaystyle \frac{\pi }{2}}{\displaystyle \frac{{(1-\nu )}^2}{(1-2\nu )}}{\displaystyle \frac{a}{h}}$	${\displaystyle \frac{a^2}{2R}}$
Bonded thin layer ($\nu =0.5$)	${\displaystyle \frac{3\pi }{64}}{\left({\displaystyle \frac{a}{h}}\right)}^3$	${\displaystyle \frac{a^2}{4R}}$
Unbonded thin layer	${\displaystyle \frac{\pi }{2}}{\displaystyle \frac{a}{h}}$	${\displaystyle \frac{a^2}{2R}}$
General layer	$\kappa (a/h,\nu ){\left({\displaystyle \frac{a}{h}}\right)}^{m(a/h,\nu )}$	${\displaystyle \frac{a^2}{(1+m(a/h,\nu \left)\right)R}}$

, an evident representation of the contact stiffness is given by

$${\overline{k}}_{\mathrm{N}}(a/h,\nu )=\kappa (a/h,\nu ){\left({\displaystyle \frac{a}{h}}\right)}^{m(a/h,\nu )},$$

(33)

where the exponent m is calculated from the relation of indentation depth and contact radius of a paraboloid with radius of curvature R at the same ratio of contact radius to layer thickness $a/h$,

$${\delta }^{\mathrm{P}}(a;a/h)=\frac{a^2}{(1+m(a/h,\nu \left)\right)R}.$$

(34)

To obtain relation (34) for all $a/h$, we can simply evaluate Equation (11) for a given parabolic profile $f\left(r\right)=r^2/\left(2R\right)$ and layer thickness h. The dimensionless parameter $\kappa$ is then calculated from Equation (33). Figure 5 shows parameters m and $\kappa$ in terms of the ratio $a/h$. Conveniently, both parameters are bounded and tend toward the asymptotic solutions with $m=0$ for the half-space, $m=1$ for the bonded compressible thin layer as well as the unbonded thin layer and $m=3$ for the bonded incompressible thin layer. The agreement with the results of Hayes et al. [11] is very good. In the ranges $a/h<0.1$ and $a/h>50$, the asymptotic solutions of half-space and thin layer, respectively, are generally good approximations. Bonded nearly incompressible and incompressible layers are slower in the asymptotic behavior. It is interesting to note that the unbonded layer agrees almost exactly with the more artificial case of a bonded layer with $\nu =0.0$. In fact, the evaluation of the BEM kernel functions of the unbonded layer (29) and the bonded layer (31) with $\nu =0.0$ gives almost the same numerical results for the whole parameter space.

Finally, the needed flat punch solution has been obtained in numerical form for the whole parameter space. With that, we show in the following section that the solution of any axisymmetric indenter shape including adhesion or viscoelasticity is straightforward using the MDR relations presented in Section 2. The numerical BEM solutions needed for evaluation of the MDR relations, that is, the normalized contact stiffness and the contact pressure distribution under a flat punch, are uploaded as Supplementary Materials to this article. Here, we discretized the $a/h$ space using 400 points with logarithmic spacing between $a/h=0.001$ and $a/h=1000$ and the normalized radial coordinate $r/a$ is discretized with approximately 500 equally spaced points.

4. Case Studies

In this section, we demonstrate the capability of the MDR formulation presented in Section 2 based on the numerical flat punch solutions presented in Section 3.2 by analyzing exemplary cases. An advantage of the MDR is that the non-adhesive elastic formulation can be easily expanded for JKR-adhesive problems and linear viscoelasticity. Here, the implementation follows the well-known rules for the half-space, which are summarized, e.g., by Willert [23]. The evaluation of the one-dimensional integrals in Equations (11) and (12) is realized by numerically efficient matrix–vector products. It is sufficient to set up the matrix kernel, a given combination of normalized contact stiffness and flat punch pressure, once for a given boundary condition (bonded or unbonded) and Poisson’s ratio (in the bonded case), since it does not depend on other simulation parameters (profile shape, adhesive properties, viscoelastic properties). The numerical BEM solution described in Section 3.2 needed to perform the following MDR calculations is uploaded as Supplementary Materials to this article.

4.1. Indentation by a Truncated Cone

Consider the indentation of a layer with thickness h by a rigid truncated cone with the profile

$$\begin{array}{c}\hfill \begin{array}{c}\hfill f\left(r\right)=\left\{\begin{array}{cc}0,\hfill & \forall r<b\hfill \\ (r-b)tan\left(\theta \right),\hfill & \forall r\ge b,\hfill \end{array}\right.\end{array}\end{array}$$

(35)

depicted in Figure 6. In the general case, the contact solution can only be obtained numerically as shown in Section 2. However, the asymptotic solutions can be derived explicitly using the relations in Section 2.1. The resulting solution for the elastic half-space can be found, for example, in [24]. For the bonded compressible thin layer and the unbonded thin layer, the layer acts as a three-dimensional Winkler foundation and the one-dimensional MDR profile and contact pressure are directly given by Equations (17) and (22). The resulting normal force is

$$F_{\mathrm{N}}=\tilde{E}\frac{\pi a^3}{3h}tan\theta \left(1-{\left(\frac{b}{a}\right)}^3\right),$$

(36)

where $\tilde{E}=E(1-\nu )/\left((1+\nu )(1-2\nu )\right)$ for the compressible thin layer and $\tilde{E}=E^{*}$ for the unbonded thin layer. In the case of a bonded incompressible layer, the asymptotic thin layer relations (20) yield the one-dimensional MDR-profile,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill g\left(x\right)=\left\{\begin{array}{cc}0\hfill & \forall x<b\hfill \\ {\displaystyle \frac{1}{3}}\left(2x+b\right){\left(1-{\displaystyle \frac{b}{x}}\right)}^{2}tan\theta \hfill & \forall x\ge b,\hfill \end{array}\right.\end{array}\end{array}$$

(37)

the normal force,

$$F_{\mathrm{N}}=\frac{E\pi a^2}{120}{\left(\frac{a}{h}\right)}^{3}tan\theta \left(2-5{\left(\frac{b}{a}\right)}^3+3{\left(\frac{b}{a}\right)}^5\right),$$

(38)

and the contact pressure,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill p\left(r\right)=\frac{Etan\theta }{36}{\left({\displaystyle \frac{a}{h}}\right)}^3\left\{\begin{array}{cc}2-{\displaystyle \frac{b^3}{a^3}}\left[6log\left({\displaystyle \frac{a}{b}}\right)+2\right]-3\left[2-3{\displaystyle \frac{b}{a}}+{\displaystyle \frac{b^3}{a^3}}\right]{\left({\displaystyle \frac{r}{a}}\right)}^2\hfill & \forall r<b\hfill \\ 2-{\displaystyle \frac{b^3}{a^3}}\left[6log\left({\displaystyle \frac{a}{r}}\right)-3\right]-3\left[2+{\displaystyle \frac{b^3}{a^3}}\right]{\left({\displaystyle \frac{r}{a}}\right)}^2+4{\left({\displaystyle \frac{r}{a}}\right)}^3\hfill & \forall r\ge b.\hfill \end{array}\right.\end{array}\end{array}$$

(39)

The solutions for the one-dimensional MDR profile and the contact pressure are shown for different confinement ratios with $b/a=0.5$ in Figure 7. The MDR profile is normalized with the profile height at the contact radius. For increasing confinement ratio $a/h$, the profile is flattened. Since the MDR profile is simply the relation of indentation depth and contact radius, this means that the contact radius of a layer is always larger then that of a respective half-space at the same indentation depth. For the bonded layer with high Poisson’s ratios and $a/h\ge 2$, the indentation depth is smaller than the profile height at $r=a$, $g\left(a\right)/f\left(a\right)<1$, meaning that the material outside of the contact zone bulges upwards, $w(r>a)<0$. For increasing confinement ratios $a/h$, the contact pressure singularity at the discontinuity of the profile at $r=b$ weakens and disappears for $a/h>5$. At the same time, the normalized stress increases in the center and decreases toward the contact edge. The asymptotic solutions are reached in all cases. However, for the bonded case with high Poisson’s ratios, agreement with the thin layer solution is achieved only for much larger $a/h$ than for the unbonded case.

4.2. Adhesive Detachment of a Cylindrical Flat Punch from an Elastic Layer

Generalized JKR-adhesion can be mapped exactly in the MDR space as shown by Heß [16]. For this, we assume that the springs adhere to the one-dimensional indenter. In the equilibrium state, the outermost spring have an maximum elongation of

$$\Delta l\left(a\right)=\sqrt{\frac{2\pi \Delta \gamma a}{c_{\mathrm{N}}(a;h)}},$$

(40)

where $\Delta \gamma$ is the work of adhesion [18]. Thus, the indentation depth for the adhesive contact is

$$\delta =g(a;h)-\Delta l(a;h).$$

(41)

With that, the normal force of the adhesive problem is again determined by Equation (4) with the one-dimensional spring displacements in Equation (2). The determination of the critical values for complete detachment follows from a simple stability analysis as shown by Heß [18]. For the following analysis, we normalize the maximum spring elongation in Equation (40) using the definitions in Equation (10),

$$\overline{\Delta l}\left(\xi \right)=\Delta l(a;h)\sqrt{\frac{E^{*}}{2\pi \Delta \gamma h}}=\sqrt{\frac{\xi }{\overline{c_{\mathrm{N}}}\left(\xi \right)}}.$$

(42)

Figure 8 shows the normalized spring elongation for the whole parameter range including the asymptotic solutions for half-space and thin layer which follow directly from Section 2.1. It is interesting to note that, for the unbonded case and the bonded case with $\nu \le 0.3$, the normalized spring elongation is monotonically increased with $a/h$, but for the bonded layer with Poisson ratio’s $\nu \ge 0.4$, it has a maximum in the range of $a/h=0.46-0.92$ ($a/h=0.46$ for the incompressible case). We show in the following that this has some important consequences for the detachment of axisymmetric indenters from elastic layers.

For the simple case of the adhesive detachment of a cylindrical flat punch of radius $a^{*}$ from an elastic layer, the normal force before detachment is given by

$$F_{\mathrm{N}}=k_{\mathrm{N}}(a^{*};h)\delta$$

(43)

and detachment begins when ${\delta }_{\mathrm{c}1}=-\Delta l(a^{*};h)$ is reached. If

$$\frac{\partial \Delta l(a;h)}{\partial a}{|}_{a=a^{*}}<0,$$

(44)

in a displacement-controlled trial, the contact area decreases stably; otherwise, unstable detachment occurs. Thus, for $a^{*}<0.46h$, detachment always occurs abruptly when ${\delta }_{\mathrm{c}1}=-\Delta l(a^{*};h)$ is reached. For the bonded layer with Poisson ratio’s $\nu \ge 0.4$ and large $a^{*}/h$, the contact area in a displacement-controlled trial decreases until

$${\delta }_{\mathrm{c}2}=-\sqrt{\frac{2\pi \Delta \gamma h}{E^{*}}}\max \left(\overline{\Delta l}(a/h)\right).$$

(45)

is reached and detaches abruptly afterwards. For bonded incompressible layers, this detachment behavior has already been found and experimentally confirmed by Webber et al. [7] (they found the value $a=0.45h$). The force displacement relation for the adhesive contact of a cylindrical flat punch and a incompressible layer is shown in Figure 9 for a flat punch with radius $a^{*}$ on layers of different thickness h. By decreasing the layer thickness, the linear force–displacement relation before detachment becomes a steeper slope due to the increased contact stiffness. Although the critical displacement ${\delta }_{c1}$ decreases for higher $a^{*}/h$, the adhesion force ($\left|\min \right(F_{\mathrm{N}}\left)\right|$) increases. As discussed, stable detachment occurs if ${\delta }_{c1}<{\delta }_{c2}$ (dotted lines in Figure 9). The absolute value of the normal force during stable detachment,

$$F_{\mathrm{N}}=-k_{\mathrm{N}}(a;h)\Delta l(a;h),$$

(46)

is always smaller than the adhesion force $F_{\mathrm{A}}=k_{\mathrm{N}}(a^{*};h)\Delta l(a^{*};h)$. Thus, in a load-controlled trial, detachment always occurs abruptly.

It should be noted that if the high confinement ratio is high, adhesive instabilities are expected during the detachment of a flat punch from a bonded layer. In this case, the contact area loses its axisymmetric and compact shape and, hence, the MDR is not applicable. Webber et al. [7] argue that the detachment mechanisms change when the pressure distribution under the flat punch has a second maximum in the center of the punch. According to Figure 4, the second maximum develops for $a/h\ge 2$ for the bonded incompressible layer.

4.3. Quasi-Static Impact with an Incompressible Viscoelastic Layer

In many cases, it is possible to solve contact problems involving viscoelastic bodies based on the solution of the elastic problem (correspondence principle, see e.g., [23]). Argatov and Popov [25] have shown that the MDR can be used to solve rebound indentation problems involving incompressible viscoelastic half-spaces, where the contact radius has a single maximum. One common option to include viscoelastic behavior is the use of rheological models that consist of linear elastic and linear viscous elements. A large range of different linear viscoelastic behaviors may be modeled in this way. In the one-dimensional MDR picture, the independent springs are simply replaced by the chosen rheological model as shown in Figure 10. To consider viscoelastic layers, the dimensionless scaling factor ${\overline{c}}_{\mathrm{N}}(x/h)$ that accounts for the geometric confinement has to be adopted from the purely elastic bedding. Using the simple example of a Kelvin–Voigt model, a parallel combination of spring and a dashpot as shown in Figure 10, the contact force for a viscoelastic incompressible layer may be calculated similar to Equation (4),

$$\begin{array}{cc}\hfill F_{\mathrm{N}}& =8\underset{0}{\overset{a}{\int }}{\overline{c}}_{\mathrm{N}}(x/h)\left[G{w}_{1\mathrm{D}}(x;h)+\eta \dot{\delta }\right]\mathrm{d}x\hfill \\ \hfill & =8G\underset{0}{\overset{a}{\int }}{\overline{k}}_{\mathrm{N}}(x/h)x{g}^{\prime }(x;h)\mathrm{d}x+8\eta a{\overline{k}}_{\mathrm{N}}(a/h)\dot{\delta },\hfill \end{array}$$

(47)

with the indentation velocity $\dot{\delta }$ and the shear modulus $G=E^{*}/4$ (incompressible).

One of the simplest cases is the quasi-static impact of a rigid cylindrical flat punch with an incompressible viscoelastic layer whose response is modeled using the Kelvin–Voigt model. Within the assumption of quasi-stationarity, inertia effects such as wave propagation are neglected. The differential equation of the impact of a flat punch with mass M and radius $a^{\mathrm{FP}}$ takes the following form:

$$M\ddot{\delta }+8{\overline{k}}_{\mathrm{N}}(a^{\mathrm{FP}}/h)\eta a^{\mathrm{FP}}\dot{\delta }+8{\overline{k}}_{\mathrm{N}}(a^{\mathrm{FP}}/h)G{a}^{\mathrm{FP}}\delta =0.$$

(48)

The easily found solution for the coefficient of restitution (COR, ratio of rebound velocity to initial velocity) is for example given by Argatov [17],

$$ϵ=exp\left[-\frac{2D}{\sqrt{1-D^2}}arctan\left(\frac{\sqrt{1-D^2}}{D}\right)\right],D<1.$$

(49)

The expressions for higher damping ratios, $D\ge 1$, are listed by Willert [23]. For the viscoelastic layer, the damping ratio according to Equation (48) is

$$D=\eta \sqrt{\frac{2a^{\mathrm{FP}}{\overline{k}}_{\mathrm{N}}(a^{\mathrm{FP}}/h)}{MG}}=D_{\mathrm{HS}}\sqrt{{\overline{k}}_{\mathrm{N}}(a^{\mathrm{FP}}/h)},$$

(50)

where $D_{\mathrm{HS}}$ is the damping ratio for the viscoelastic half-space. We can interpret the damping ratio as the ratio of relaxation time of the Kelvin–Voigt element $\tau =\eta /G$ and the impact duration in the purely elastic case. For this simple impact case, the enlarged stiffness due to the finite geometry reduces the impact duration of the elastic problem and thus increases the damping ratio. Since Equation (49) is a monotonic decreasing function of the damping ratio, the COR decreases with the confinement ratio. Figure 11 shows the COR in terms of the damping ratio for the elastic half-space and the confinement ratio for bonded and unbonded layers. For $a/h\le 0.1$, the half-space solution is a good approximation. Due to the much larger dimensionless contact stiffness for bonded layers with $a/h>1$, the COR has a significant dependence on the boundary condition at the interface.

A less simplistic impact problem is that of a rigid sphere with an incompressible viscoelastic layer that is modeled with the well-known standard model shown in Figure 12, where usually $G_{\infty }\ll G_1$. If the confinement factor ${\overline{c}}_{\mathrm{N}}(x/h)$ is taken into account (as for the Kelvin–Voigt model described above), the simple explicit time integration procedure described by Willert et al. [26] for the impact with a viscoelastic half-space can be applied. The only difference is the spatially varying factor for the reaction force of a single element at position $x_i=i\Delta x$,

$$f_i=4\Delta x{\overline{c}}_{\mathrm{N}}(x_i/h)\left[G_{\infty }{w}_i+\eta ({\dot{w}}_i-{\dot{\overline{w}}}_i)\right],$$

(51)

where ${\overline{w}}_i$ denotes the displacement of the inner point that is determined from the equilibrium condition,

$$G_1{\overline{w}}_i+\eta ({\dot{\overline{w}}}_i-{\dot{w}}_i)=0.$$

(52)

For the impact of a spherical indenter with a viscoelastic half-space, the ratio $G_{\infty }/G_1$ and the dimensionless variable

$${\chi }^{\mathrm{HS}}=\eta {\left(\frac{R{v}_0}{M^2{G}_{\infty }^3}\right)}^{1/5}$$

(53)

are defining [23]. For the asymptotic thin layer cases described in Section 2.1, we can easily find the defining dimensionless variables for incompressible layers,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\chi }^{\mathrm{TL}}& =\eta {\left(\frac{R{v}_0}{hM{G}_{\infty }^2}\right)}^{1/3},\mathrm{for}\mathrm{the}\mathrm{unbonded}\mathrm{thin}\mathrm{layer},\hfill \\ \hfill {\chi }^{\mathrm{TL}}& =\eta {\left(\frac{R^2{v}_0^2}{h^{3}M{G}_{\infty }^3}\right)}^{1/4},\mathrm{for}\mathrm{the}\mathrm{bonded}\mathrm{thin}\mathrm{layer}.\hfill \end{array}\end{array}$$

(54)

It is impossible to find similar simple defining variables outside of the asymptotic ranges, if the contact area changes during the impact. Thus, the influence of the impact parameters on the finite geometry is more complex for the sphere.

In Figure 12, we used the parameters $G_{\infty }/G_1=0.\overline{11}$ and ${\chi }^{\mathrm{HS}}\approx 2.65$ and varied the layer thickness. Then, it can be shown numerically that impact velocity $v_0$ and radius R of the indenter are free parameters. The other non-free parameters were chosen as $\tau =\eta /G_1=1$ ms, $G_1\approx 7.5$ MPa, $M=20$ g and layer thicknesses $h=\{1\mathrm{mm},2\mathrm{mm},3\mathrm{mm}\}$. Figure 12A shows typical hysteresis curves of normal force in terms of indentation depth, both normalized with the maximum value of the elastic half-space with $G=G_{\infty }+G_1$. Figure 12B shows the confinement ratio $a/h$ and Figure 12C shows the normalized velocity $v/v_0$, both in terms of the normalized impact time. For the elastic layer and, with the given parameter combination, also for the viscoelastic layer, the impact duration decreases with decreasing layer thickness. The maximum confinement ratio of approximately 1 is found for the 1 mm layer and here the maximum normal force is twice as high compared to the viscoelastic half-space. The contact radius of the elastic and viscoelastic layer differs almost only during the retraction phase, which is shorter then the indentation phase for the viscoelastic layer. The COR at the given parameter variation can be read from Figure 12C. It is approximately 0.45 for the viscoelastic half-space and 0.6 for the 1 mm thick viscoelastic layer, and thus increases with decreasing layer thickness. However, for other parameter combinations of $G_{\infty }/G_1$ and ${\chi }^{\mathrm{HS}}$, opposing trends can be obtained.

Since the validity of the MDR for rebound problems with viscoelastic media has been formally proven only for the half-space geometry [25], we also simulated the quasi-static impact problem of sphere and viscoelastic layer with the finite element method (FEM) using ABAQUS Explicit. For the material behavior, we chose the elastic parameters in combination with a Prony series corresponding to the standard medium mentioned above. Figure 12 shows that FEM results (dashed lines) and MDR results (solid lines) agree almost exactly. It should be noted that the BEM formulation described in Section 3.1 could also be used to solve the viscoelastic impact problem by applying the correspondence principle. For the impact with a viscoelastic half-space, BEM and MDR have been shown to agree except for negligible numerical errors, with MDR having an enormous advantage in terms of computational time [26]. In this newly proposed MDR implementation, where the elastic flat punch solution obtained with BEM is the basis for the MDR solution, FEM is used as a completely independent method for comparison. The results indicate that the MDR in the present form is indeed valid for rebound problems with viscoelastic layers. The MDR algorithm for the viscoelastic impacts in Figure 12, implemented in Python on a regular laptop (Intel i7 processor), took less than one second for a single impact. In comparison, the FEM implementation is poorly parameterizable and much slower (several minutes of computational time), so the proposed MDR routine is very useful for comprehensive parameter studies.

5. Conclusions

In this paper, a simple, versatile and numerically efficient semi-analytical solution method for axisymmetric contact problems involving bonded and unbonded layers is presented and applied for exemplary cases. The method represents a recently published generalized version of the method of dimensionality reduction (MDR) [14], which is based on the pre-required contact solution of a cylindrical flat punch. In the present paper, this flat punch solution is obtained numerically using the boundary element method (BEM) for a single bonded or unbonded layer. It is provided as Supplementary Materials to this paper. With this flat punch solution, the contact solution for arbitrary axisymmtric indenter shapes and arbitrary layer thickness is reduced to the numerical evaluation of one-dimensional integrals. Furthermore, we show that solely based on this provided flat punch solution, the adhesive contact with elastic layers and the influence of viscoelasticity can easily and efficiently be investigated using the well-established MDR formalisms. In the application examples, qualitative and quantitative effects due to the finite geometry are elaborated in extensive parameters studies. Comparisons with the literature and finite element simulations show very good agreement. For the asymptotic case of a bonded incompressible thin layer, closed-form MDR expressions for arbitrary indenter shapes are presented for the first time.

Due to the possibility of rapid parameter studies, the proposed semi-analytical method will be very helpful to study thickness effects in various applications such as contacts with cartilage layers. The stress tensor and hydrostatic pressure gradient inside the layer can be determined using the superposition idea as was performed for the half-space geometry [27]. It is easy to adapt the proposed MDR procedure for other interesting contact problems including bodies with multiple layers or functionally graded layers. The extension for slightly non-axisymmetric contacts involving elastic layers applying the principles of the work by Popov [28] is planned for the future.