PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Int J Eng Sci. Author manuscript; available in PMC 2010 November 1.
Published in final edited form as:
Int J Eng Sci. 2009 November 1; 47(11): 1274–1283.
doi:  10.1016/j.ijengsci.2008.08.005
PMCID: PMC2808696
NIHMSID: NIHMS106455

Indentation of an elastic half space with material properties varying with depth

Abstract

We consider the effect of an elastic modulus that decreases with depth on the load-displacement relation for indentation of a graded half space by a rigid indenter. A closed-form approximation incorporating features of the plate on an elastic substrate and the Hertzian contact theory is compared with finite element results for the case of a uniform stiff layer on a homogeneous substrate. Some general results are presented for the case where the grading has inverse power-law form and the effects of truncation to a finite surface value are investigated numerically. Finally, a more practical error-function grading is considered. In all cases, the load-displacement relation is closer to linear than in the homogeneous case. We conclude that the experimental data can be used to determine parameters in a predetermined form of grading, but that comparative insensitivity to the exact form of the grading would make it difficult to distincguish experimentally between different models based on indentation experiments alone.

Keywords: layered materials, indentation, inhomogeneous material, contact problems, functionally graded material

1 Introduction

Modern technological developments have led to considerable interest in materials whose properties vary with depth [1], [2]. Applications range from barrier coatings designed to resist extreme environments [3] to the development of custom layers on polymers in the development of biodiagnostic devices [4]. The variation in properties may either be in the form of a continuous function of depth as in a graded material, or as one or more discrete layers on the surface of a homogeneous substrate. The resulting layers can be extremely thin, which makes it very difficult to perform direct measurements of their mechanical properties. However, it has long been recognized that indentation tests can be used to extract more mechanical properties than merely the hardness [5]. In particular, the atomic force microscope (AFM) and the nano-indenter permit this technique to be extended to examine materials layered in the nanometer range, for which other measurement techniques would be impractical [6]. Techniques and analyses for extracting the elastic properties of homogeneous materials by indentation are well established [7], and many modifications to indentation analyses have been proposed to account for the effects of layered materials [8], [9], [10].

In the present paper, attention is focused on the effects of a substrate that is significantly more compliant than the surface. Elementary considerations suggest that the load-displacement relation should reflect elastic properties in a region of material comparable in thickness to the linear dimensions of the contact area, so that by examining this curve over a range of indentation depths, we might hope to obtain information about the way in which the elastic properties are graded with depth. However, if the elastic modulus decreases substantially with depth, the system will behave somewhat like a plate on an elastic foundation and the indentation load-displacement relation will be influenced by properties of the substrate even when the contact area is extremely small. In this paper, we shall explore the extent to which the experimentally measured load-displacement relation can be used to deduce information about the grading of elastic modulus in such cases.

The topic is introduced by considering limiting analyses for the mechanics of indentation for a discrete stiff layer on a compliant substrate, as would be appropriate for a metal layer on a polymer. This is followed by a general study of how the indentation behavior is affected if the surface is graded, with a continuous increase in compliance away from the surface, as would be appropriate for a polymer with a surface stiffened by a chemical reaction such as oxidation.

2 The uniform layer

The simplest situation is that in which a uniform layer of thickness h and Young's modulus Ef is bonded to a uniform half space of a more compliant material with modulus Es. If Ef [dbl greater-than sign] Es and the layer is in some sense thin, it is natural to expect the behaviour of this system to be well-approximated by an elastic plate of stiffness

D=Efh312(1νf2)
(1)

supported on an elastic half space representing the substrate. Timoshenko and Woinowsky-Krieger [11] give the solution for the case where a concentrated normal force P is applied to the surface of such a plate, but the interface between the plate and the half space is frictionless. The surface displacement at a distance r from the force is

w(r)=Pl22πD0J0(λρ)dλ(1+λ3),
(2)

where

ρ=rlandl=2D(1νs2)Es3.
(3)

If the plate is instead bonded to the half space, this will tend to reduce the tangential displacement at the surface of the substrate and hence increase its effective modulus. In the Appendix we consider the extreme case where the in-plane stiffness of the layer is sufficient to prevent any tangential displacement at the interface and show that l is then modified to

l=(D(1+νs)(34νs)2Es(1νs))1/3.
(4)

This is identical to (3:ii) for the case νs = 0.5 and is only 10% lower when νs = 0. For a plate of finite in-plane stiffness, the effective value of l will be intermediate between these quite tight bounds. Substututing for D from (1), we can define the dimensionless parameter

l^lh=(Ef(1+νs)(34νs)24Es(1νs)(1νf2))1/3,
(5)

which is a measure of the modulus mismatch between the layer and the substrate.

2.1 Local displacement and curvature

Ol'shanskii [12] gives a series solution to the integral (2) in the form

w(r)=Pl212Dm=0(1)m{43[ρ6m(3m)!(3m)!ρ6m+4(3m+2)!(3m+2)!]+3ρ6m+5Γ(3m+7/2)Γ(3m+7/2)+6ρ6m+2π(3m+1)!(3m+1)![ln(ρ)ψ(3m+2)]},
(6)

where ψ is Euler's psi function ([13] §8.36). For sufficiently small values of ρ (i.e. points sufficiently near the load), (6) can be approximated by the first three terms

w(r)=Pl212D{43+6ρ2π[ln(ρ)+C1]},
(7)

where we have used the result ψ(2) = 1 − C from [13]:§8.365.4 and C = 0.577215… is Euler's constant. Numerical calculations show that (7) is a good approximation to (6) in the range 0 ≤ ρ < 0.5.

The Laplacian of (7),

2w=d2wdr2+1rdwdr=2PπD{ln(η)+C},
(8)

provides a measure of the curvature of the initially plane plate surface. This curvature is unbounded as r → 0, showing that a finite contact area must be established for any rigid indenter of finite radius R.

2.2 Indentation by a rigid sphere

If the plate is indented by a rigid sphere of radius R, the resulting contact pressure will comprise a line load of P/2πa per unit length distributed around a ring of radius a. To demonstrate this, we use (8) as a Green's function to determine the curvature due to such a ring load. The constant term will clearly contribute a constant curvature at all points on the surface and it is a well-known result of two-dimensional potential theory that with a logarithmic Green's function, the potential due to the uniform ring source is equal to that of the corresponding point source at the origin for r > a and is constant for r < a. Superposing these results, we obtain

2w=2PπD{ln(rl)+C};r>a
(9)
=2PπD{ln(al)+C};0r<a.
(10)

It follows that the profile of the plate surface due to the ring force will conform exactly to that of the indenting sphere and leave a positive gap outside r = a as required if

2PπD{ln(al)+C}=2Roral=exp(CπP^),
(11)

where we define the dimensionless load

P^PRD.
(12)

Once the radius of the ring is found from this equation, the displacement at the center (and hence the indentation of the punch) is readily found by superposition using (7) as a Green's function. We obtain

w(0)=Pl212D{43+6a2πl2[ln(al)+C1]}.
(13)

Using (11) to eliminate a and introducing the dimensionless parameters (5, 12), we can also write this result in the form

w(0)=P^l^2h26R{233π[1+πP^]exp(2C2πP^)}.
(14)

2.3 The Hertzian solution

At small values of the dimensionless load P, the contact radius (11) is extremely small and the load-displacement relation (14) is virtually linear. However, in this range, the actual contact area and contact pressure distribution will be more accurately defined by the Hertzian theory for a body composed entirely of the layer material. In this limit, the Hertzian contact radius is

aH=(3(1νf2)PR4Ef)1/3oraHh=12(P^2)1/3
(15)

[14], using (12, 1). The Hertzian indentation depth d is given by

dH=aH2R=h24R(P^2)2/3
(16)

[14], using (15, 12). These results can be regarded as a correction that can be ‘patched in’ to the plate solution (14) at small values of P. In particular, the rigid body motion of the indenter d will then be the sum of the plate displacement w(0) and the Hertzian compliance dH, giving

d^dRh2=P^[l^233l^22π[1+πP^]exp(2C2πP^)+14(14P^)1/3].
(17)

This superposition is strictly meaningful only when the Hertzian radius aH [double less-than sign] h and hence P [double less-than sign] 1. However, for larger values of P, the last (Hertzian) term in (17) becomes small compared with the remaining terms and hence this equation might be expected to give reasonable results as long as the contact radius is sufficiently small to justify the use of the three-term approximation (7).

2.4 Dimensional considerations

Notice that the radius of the indenter R appears only as a linear multiplier in equation (17) and can be subsumed into the dimensionless indentation depth d. Thus, the dimensionless load-displacement relation depends only on the single dimensionless material parameter l and in particular is independent of the dimensionless ratio R/h. This reduction in parameter dependence also occurs for the exact elasticity solution, provided that the indenter is replaced by the equivalent paraboloid so that the normal displacement in the contact region is

uz(r,0)=dr22R.
(18)

If we normalize the length dimensions with h and the displacements with h2/R, defining

u^z=uzRh2;d^=dRh2;r^=rh,

equation (18) takes the form

u^z(r^,0)=d^r^22

and the resulting boundary-value problem in these variables clearly depends only on the dimensionless contact radius a/h and on the dimensionless material parameters l, νf, νs. Also P of equation (12) is an appropriate dimensionless indentation force resulting from this normalization.

2.5 Finite element results

To evaluate the approximations involved in the plate solution, an axisymmetric finite element solution of the problem was developed. Rectangular elements of side 0.1h were used in the layer, with coarser meshing distant from the contact region. The total depth of the modeled region is 6250h, and the maximum radius is 2500h. Fig.1 compares the predictions of equation (17) with the finite element results for l = 3. Agreement is excellent up to P ≈ 3.5, d ≈ 6, at which load equation (11:ii) predicts a/l ≈ 0.7. We also note that the contact pressure distribution exhibits a peak that moves towards the outer radius of the contact region in this load range.

Fig. 1
Finite element results (circles) and theoretical prediction (solid line) from equation (17) for indentation of a uniform layer bonded to an elastic half space.

3 Continuously graded material properties

The uniform layer solution is appropriate to situations in which a layer of a stiffer material is deposited on a more flexible substrate. However, in many cases, stiff surface layers are generated by chemical changes at the surface of an initially homogeneous material [15] and it is then reasonable to expect the modulus to be a continuous function of depth. If the modulus at the surface is much higher than that of the substrate, the qualitative behaviour might still be expected to resemble that of a plate on an elastic substrate. However equations (14, 17) require us to determine an equivalent stiffness D for the ‘layer’ and this in turn requires a decision as to the point where the layer is considered to stop and the substrate to start. Unfortunately, this essentially arbitrary choice has a substantial effect on the calculated value of D, because it involves the second moment of the modulus distribution across the layer thickness.

3.1 Power-law grading

Information about the nature of the grading can be obtained from the shape of the experimental load-displacement curve. As background to this, it is instructive to consider the case where the grading has power law form and the rigid indenter has a power law profile. We assume the most general form of anisotropic Hooke's law

σij=cijklukxl,
(19)

with the modulus defined by

cijkl=x3λBCijkl,
(20)

where Cijkl is a dimensionless tensor and B is a dimensional constant, so that all the components of the elasticity tensor follow the same grading. It should be remarked that this idealized grading function implies unrealistic behaviour at x3 = 0 and x3 → ∞ for all values of λ except the uniform case λ = 0. In particular, for λ < 0, the modulus tends to zero as x3 → ∞. However, we shall demonstrate later in §3.2 that more realistic grading functions with a power law central segment but finite limiting values exhibit behaviour close to that of the idealized case. We consider the case where the half space x3 ≥ 0 is indented by a frictionless rigid punch with a power-law conical profile. A formal statement of the boundary value problem is

g(x1,x2)u3(x1,x2,0)d+rβf(θ)0
(21)
p(x1,x2)σ33(x1,x2,0)0
(22)
σ31(x1,x2,0)=σ32(x1,x2,0)=0
(23)
p(x1,x2)g(x1,x2)=0,
(24)

where d is the indentation, g is the gap between the indenter and the surface of the half space and the function f(θ) describes a representative cross section of the punch in polar coordinates x1 = r cos θ, r sin θ = x2.

This problem has no intrinsic length scale and hence the solution for all loads must be similar to each other. This fact can be exposed by defining the dimensionless parameters

ξi=d1/βxi;ρ=d1/βr;Uj=ujd;Sij=d(1λβ)/βσijB,

in terms of which equations (19, 2124) take the form

Sij=Cijklξ3λUkξl
(25)
γ(ξ1,ξ2)U3(ξ1,ξ2,0)1+ρβf(θ)0
(26)
S33(ξ1,ξ2,0)0
(27)
S(ξ1,ξ2,0)=S(ξ1,ξ2,0)=0
(28)
S(ξ1,ξ2)γ(ξ1,ξ2)=0,
(29)

This problem is clearly independent of the punch indentation d. Once it is solved, the total force required to produce a given indentation can be obtained by considering the equilibrium of the region 0 < x3 < d1h, giving

P=σ33(x1,x2,d1/βh)dx1dx2=d(1+λ+β)/βBS33(ξ1,ξ2,h)dξ1dξ2.
(30)

The integral is independent of d and hence the load-displacement relation must have the form

Pdα,
(31)

where

α=(1+λ+β)β.
(32)

Similar considerations show that the contact area retains the same shape at all values of P and that a typical linear dimension a of the contact area increases as

ad1/βP1/(1+λ+β).

For the case of a paraboloidal punch, β = 2 and we have P ~ d(λ+3)/2, a ~ d1/2 ~ P1/(λ+3). These results agree with the classical Hertzian analysis [14] in the homogeneous case λ = 0 and with the solution of the isotropic axisymmetric indentation problem due to Giannakopoulos and Suresh [16]. The latter authors claim that their solution applies only to the case where λ > 0 — i.e. when the modulus increases with depth — but there seems to be no basis to this restriction. In the special case where the modulus increases linearly with depth (λ = 1), the half space mimics a Winkler elastic foundation [18]. This result was first remarked by Gibson [19].

Notice that the non-linearity of the load-displacement relation is reduced as λ becomes more negative and a completely linear relation is obtained for all power-law punch profiles (all β) when λ = −1. A physical explanation for this latter result is provided by the point force solution of Giannakopoulos and Suresh [17]. They show that the normal surface displacement uz(r, 0) at a distance r from the point of application of the force varies with r−λ−1. Thus, when λ → −1, the normal surface displacement becomes independent of r. In other words, a point force causes the entire surface to deflect downwards by the same distance (though this distance tends to zero in the limit). It follows that a rigid punch of any shape (not necessarily power law or axisymmetric) will make only point contact with the surface and in fact all such punches will have the same force-displacement relation.

3.2 More realistic grading

The preceding results suggest that some information about the nature of the grading can be obtained by plotting the experimental load-displacement data on a log-log scale and approximating the resulting curve by a straight line. Of course, any power law grading other than uniform (λ = 0) implies an unrealistic zero or infinite modulus at the surface and at infinite depth, so it is important to estimate how the truncation of equation (20) at large and small values of x3 might influence the resulting curve.

This question was investigated using a finite element approximation to the problem. In view of the localized loading, a graded mesh was used, with the smallest elements near the surface being 100 times thinner than the coarsest elements distant from the surface. A piecewise constant modulus distribution was generated from the power law function, which implies that the layer adjacent to the surface has the (finite) modulus appropriate to its mid-point. The material was assumed to be isotropic and incompressible (ν = 0.5). No explicit truncation was imposed at large depths, but of course the finite size of the model implies a form of truncation. This effect was explored by changing the total depth of the model and the results showed less than 1% change in the load-displacement curve with a change of a factor of five in the total depth.

Fig.2(a) shows the load-displacement relation for indentation of a power-law graded half space with λ = −0.5 by a rigid paraboloidal punch whose tip radius is 1250 times smaller than the total depth of the modeled region. A good straight line approximation can be fitted on the logarithmic scale. Results were obtained for several different values of λ in the range −1 < λ < 1 and the resulting slopes are compared with the theoretical prediction (32) in Fig.2(b). We conclude that truncation due to discretization has only a small effect on the results, and hence the ‘unphysical’ values of modulus implied by a strict power law do not preclude the use of these results in cases of approximately power-law grading. Incidentally, these finite element results can be regarded as a ‘practical’ confirmation that the results of Giannakopoulos and Suresh [16] can indeed be extended to negative values of λ > −1.

Fig. 2
(a) Load-displacement relation for the indentation of a half space with modulus proportional to x31/2. (b) Slope of the logarithmic load-displacement relation as a function of the modulus exponent λ. Finite element results are represented ...

To separate the effects of discretization and power law truncation, we next consider a half space with the explicitly truncated modulus distribution

E(x3)=E0[1+λ2λ(x3h)2];0x3<h=2E02λ(x3h)λ;x3>h.
(33)

where E0 is the modulus at the surface and h is the thickness of the region in which the power-law distribution is truncated. This distribution preserves continuity of modulus and gradient at x3 = h. The indentation problem now ceases to be self-similar because of the length parameter h, but dimensional considerations similar to those in §2.4 show that for indentation by a paraboloidal punch, the dimensionless indentation d must depend only on the dimensionless load

P=PRE0h3.
(34)

As in the case of the uniform layer, we anticipate that the load-displacement relation will be dominated by local ‘Hertzian’ stresses corresponding to the surface value of the modulus for P~ [double less-than sign] 1, but that the substrate will play an increasingly important rôle as P~ increases. This suggests a relation P~ ~ d3/2 for P~ [double less-than sign] 1 and P~ ~ d5/4 for P~ [dbl greater-than sign] 1, based on equation (32).

Fig.3 shows finite element results for the case λ = −0.5. The two limiting slopes 32 and 54 are indicated by the solid lines and the results confirm that a transition between these limiting behaviours occurs predominantly in the range 0.03 < d < 1. In particular, the load-displacement relation is dominated by the substrate power-law grading when

Fig. 3
Load-displacement relation for the truncated power-law grading of equation (33) with λ = −0.5.
d^>1ordh>hR.

Notice that since the dimensionless indentation contains the radius R of the indenter, this places restrictions on the ability of the measurement to probe the modulus very close to the surface. For example, if the AFM tip radius is 20 nm and the displacement d is measured with an accuracy of ±8 nm, a modulus distribution of the form (33) would give results essentially indistinguishable from pure power-law grading if h < 15 nm.

3.3 Error function grading

If the stiff surface layer is generated by a chemical reaction, it is likely that the resulting modulus will be some function of the concentration of a reactant and the period to which the local material is exposed to the reactant at an appropriate temperature. This in turn is likely to be determined by an appropriate diffusion equation, suggesting the possibility of error-function grading

E(x3)=Es+(E0Es)erfc(x3h),
(35)

where E0 is the modulus at the surface x3 = 0 and h is a characteristic thickness dimension for the stiffened layer.

Fig.4 shows the dimensionless load PR/Esh3 as a function of d for various values of the modulus ratio E0/Es. Notice that the substrate modulus Es is used here in normalizing the load P to permit us to present disparate curves on the same plot. A striking feature of these results is the almost linear form of the load-displacement relation, particularly at large modulus ratios. For example, at E0/Es = 103, the best straight line fit to the corresponding logarithmic plot has a slope of about 1.04. By contrast, the power law form of equations (20, 33) give a strictly linear load-displacement relation only in the limit λ → −1.

Fig. 4
Load-displacement relation for the error function grading of equation (35) and various values of E0/Es.

3.4 Effect of residual stress

All of the theoretical and finite element results reported in this paper are based on the assumption that the graded half space is stree-free in the unloaded state. However, if the graded modulus results from a chemical or physical change in the surface layers, it is possible for this same process to generate a state of residual stress, which could have a significant effect on the load-displacement relation. In particular, if a thin surface layer is much stiffer than the substrate and is in a state of residual biaxial tension, these membrane stresses will tend to stiffen the apparent load-displacement relation.

To assess the importance of this effect, finite element calculations were performed for a uniform layer in a state of uniform biaxial tensile strain ε11 = ε22 = ε0, ε33 = 0 and the results were compared with those reported in Figure 1. The load-displacement curve was essentially unaffected by residual strains of 10−3 or less, but a strain of ε0 = 0.005 increases P by about 30% in the range 0 < d < 20. Thus, if the layer formation process is likely to produce significant residual tensile strains, it is important that these be taken into account in interpreting the indentation data.

4 Conclusions

The classical Hertzian analysis for the indentation of a homogeneous elastic half space by a paraboloidal punch predicts that the load will vary with indentation to the power 3/2. The various results presented in this paper all show that when the modulus of the half space decreases with depth, this relation becomes closer to linear. Physically, this effect might be explained by noting that the stiffening effect of an increased contact area is to some extent offset by the fact that the resulting contact stress field penetrates to a greater depth where the modulus is lower. Also, when the surface layers are much stiffer than the substrate, a significant part of the indentation results from ‘plate’ deformation of this stiff layer.

For power-law grading, the load-displacement relation is itself a power law of the form of equation (32). However, this implies an unphysical infinite modulus at the surface and if this condition is relaxed by defining a thin layer in which the modulus passes to a finite surface value, the relation transitions to the Hertzian at small values of dimensionless indentation d.

In principle, one might expect that the information contained in a sufficiently accurate experimental measurement of the load-displacement relation could be used to define an inverse problem for the determination of the modulus as a function of depth. However, the results obtained here suggest that this inverse problem, even if well-posed, would be ill-conditioned in regard to experimental variance. For example, the load-displacement relation for error-function grading (Fig.4) with E0/Es = 1000 is very close to that obtained with power-law grading with λ = −0.92, and the difference between these two dissimilar forms of grading would very probably lie within the range of experimental uncertainty.

We conclude that indentation measurements can be used to estimate parameters in a predetermined grading function for the modulus such as a truncated power law or an error function, but that more detailed information about the grading would require the indentation results to be supplemented by additional experimental information, such as for example the surface wrinkling wavelength obtained in a compression text [20].

Acknowledgments

This work was supported in part by NIH (EB003793-01). The authors also wish to thank Professor N.Triantafyllidis for numerous useful discussions on this and related topics.

Appendix

If the interface between the layer and the substrate is frictionless, the substrate is loaded only by normal tractions and the elastic field can be defined in terms of a single harmonic function [var phi], using Solution F of [21]. In particular, the relevant surface tractions and displacements are

uz(r,0)=2(1νs2)Esφz(r,0);σzz(r,0)=2φz2(r,0)
(36)

At the opposite extreme, if the substrate is bonded to a plate with sufficient in-plane stiffness to prevent all tangential surface displacement, the most general stress state in the substrate can be defined in terms of the harmonic function ω of Solution B [21], for which

uz(r,0)=(34νs)(1+νs)Esω(r,0);σzz(r,0)=2(1νs)ωz(r,0).
(37)

The expressions for uz(r, 0) in (36, 37) can be made identical by the substitution

ω=2(1νs)(34νs).
(38)

Making the same substitution in the expressions for σzz(r, 0), we obtain results of identical form except for the multiplying constants. We conclude that for a given value of contact traction σzz(r, 0), the displacement in the frictionless case will exceed that in the radially restrained case in the constant ratio

4(1νs)2(34νs).
(39)

This also implies that the radial restraint is equivalent to an increase in the apparent elastic modulus Es in the ratio (39) relative to the frictionless case. In particular, Timoshenko's solution (2) remains valid if l is replaced by

l=(D(1+νs)(34νs)2Es(1νs))1/3.
(40)

References

1. Mian MA, Spencer AJM. J Mech Phys Solids. 1998;46:2283–2295.
2. Selvadurai APS. Int J Num Anal Meth Geomechanics. 1996;20:351–364.
3. Evans AG, Mumm DR, Hutchinson JW, Meier GH, Pettit FS. Prog Matls Sci. 2001;46:505–553.
4. Zhu X, Mills KL, Peters PR, Bahng JH, Liu EH, Shim J, Naruse K, Csete ME, Thouless MD, Takayama S. Nature Matls. 2005;4:403–406. [PubMed]
5. Loubet JL, Georges JM, Marchesini O, Meille G. ASME J Trib. 1984;106:45–48.
6. Doerner MF, Nix WD. J Matls Res. 1986;1:601–609.
7. Oliver WC, Pharr GM. J Matls Res. 1992;7:1564–1583.
8. Gao H, Chiu CH, Lee J. Int J Solids Struct. 1992;29:2471–2492.
9. Gao YF, Xu HT, Oliver WC, Pharr GM. J Mech Phys Solids. 2008;56:402–416.
10. Xu ZH, Rowcliffe D. Thin Solid Films. 2004;447-448:399–405.
11. Timoshenko SP, Woinowsky-Krieger S. Theory of Plates and Shells. 2nd. McGraw-Hill; New York: 1969.
12. Ol'shanskii VP. PMM. 1987;51:681–683.
13. Gradshteyn IS, Ryzhik IM. Tables of Integrals, Series and Products. Academic Press; New York: 1980.
14. Johnson KL. Contact Mechanics. Cambridge University Press; Cambridge: 1985.
15. Mills KL, Zhu XY, Takayama SC, Thouless MD. J Matls Res. 2008;23:37–48. [PMC free article] [PubMed]
16. Giannakopoulos AE, Suresh S. Int J Solids Struct. 1997;34:2393–2428.
17. Giannakopoulos AE, Suresh S. Int J Solids Struct. 1997;34:2357–2392.
18. Calladine CR, Greenwood JA. Quart J Mech Appl Math. 1978;31:507–529.
19. Gibson RE. Géotechnique. 1967;17:58–67.
20. Lee D, Triantafyllidis N, Barber JR, Thouless MD. J Mech Phys Solids. 2008;56:858–868. [PMC free article] [PubMed]
21. Barber JR. Elasticity. 2nd. Kluwer; Dordrecht: 2002.