Analytic solution for the indentation of a transversely isotropic elastic layer bonded to a rigid foundation

D. H. Cortés¹ and J. J. García²

¹ Dto. de Ing. Mecánica ,Universidad del Valle,Cali, Colombia
dhcortes@petecuy.univalle.edu.co
² Dto. de Ing. Civil, Universidad del Valle,Cali, Colombia
josejgar@mafalda.univalle.edu.co

Abstract — The analytical solution for a transversely isotropic linear elastic layer bonded to a rigid foundation and indented by a rigid cylinder or sphere was developed. This solution follows procedures used by others to solve contact problems of linear anisotropic materials. The solution can be used to find the stress distribution and the displacement field in anisotropic layers like articular cartilage. The solution was used to compare stresses and displacements in articular cartilage assuming two sets of engineering properties with different degrees of anisotropy. The results may support current research about the relation between impact loading on the articular cartilage and the development of osteoarthrosis.

Keywords — Indentation. Transversely Isotropic. Articular Cartilage. Osteoarthrosis.

I. INTRODUCTION

Articular cartilage is the material covering the end of the bones inside the synovial joints. It has excelent mechanical properties to transmit loads and to allow relative movements without significant wear (Mow et al., 1980). In some cases, articular cartilage begins to deteriorate and the underlying bone grows until direct contact bone to bone is produced inside the joint causing inmobilization and severe pain. This desease, known as osteoathrosis, can lead to considerable hospital stays and extended periods of lost work days (Mackenzie et al., 1988).

Significant research has been undertaken to understand the ethiology of osteoartrosis (Ewers et al. 1998). Some authors have suggested that osteoarthrosis is due to impact loading to the joint and in vivo experiments in animals have been used to study this correlation (Newberry et al., 1998). The most widely used mechanical model for the articular cartilage is the biphasic (Mow et al., 1980), which considers the tissue composed of a solid phase and a fluid. Under impact loading and equilibrium, biphasic cartilage can be analyzed as an equivalent elastic layer (García et al., 1998). It has been shown that an isotropic model for articular cartilage is unable to simulate the response in indentation experiments (Mow et al., 1989). On the other hand, a transversely isotropic model in which the Young's modulus in the plane of the cartilage (plane of isotropy) is higher than that in the direction of the loading, provides a good fit to the experimental curves (Cohen et al., 1993; García et al., 2000).

In situ indentation tests with spherical and cylindrical indenters have been widely used to determine the elastic properties of articular cartilage (Töyräs et al., 2001), which are necessary to asses the condition of the tissue in animal experiments and to undertake finite element analysis of the joints. If a transversely isotropic model is adopted for the articular cartilage, there is no analytical elastic solution for this layer firmly bonded to the rigid foundation. In this study, the analytical solution for the indentation of a transversely isotropic elastic layer bonded to a rigid foundation was developed based on the general equations of the anisotropic elasticity presented by Lekhnitskii (1981) and procedures followed by others to solve linear elastic contact problems (Sakamoto et al., 1991). This solution may help to analyze articular cartilage under impact loading and to develop procedures to determine their elastic constants from in situ tests.

II. METHODS

A. Problem formulation

The model to be solved consists of a transversely isotropic linear elastic layer bonded to a rigid foundation and indented by a cylindrical or spherical punch. Isotropic planes are perpendicular to z-axis (Fig. 1), which is the axial symmetry axis. A cylindrical coordinate system was used with the origin located at the intersection between the rigid foundation and the symmetry axis. The material is characterized by five elastic constants, which are the Young's modulus and Poisson's ratio (E, Q) in the isotropy plane and the Young's modulus, shear's modulus and Poisson's ratio out of the isotropy plane (E', Q' and G'). These engineering constants can also be related to the elasticity coefficients C₁₁, C₃₃, C₄₄, C₁₃ and C₁₂ used by Lekhnitskii (1981) as shown in Appendix A.

This model can be used to represent the mechanical behavior of articular cartilage, firmly attached to the underlying subchondral bone which is much more rigid than the cartilage.

Figure 1. Problems to solve

Boundary conditions are as follows. Normal stresses on stresses can be calculated using Eqs. (7)-(12) as follows the surface outside of the contact area as well as shear stresses on the entire surface are zero. Displacements inside of the contact area are known and displacements in the bottom plane are zero since the layer is bonded to a rigid foundation. These conditions can be written as follows:

(w)_z=h = -ε₀,+ f(r) (0 £ r £ r₀), (1)

(σ_z)_z=h= 0, (r > r₀) (2)

(τ_rz)_z=h = 0, (3)

(w_z)_z=0 = 0, (4)

(u)_z=0 = 0, (5)

where w(r,z) and u(r,z) are the displacement functions in the z and r directions, respectively. Letters σ and τ are used for normal and shear stresses respectively. The r function f(r) = 0 applies for the cylindrical indenter and f(r) = r²/2R applies for the spherical indenter, where R is the sphere radius.

Linear elasticity equations can be formulated as proposed by Lechnitskii (1981) for axial symmetry and anisotropic materials. two functions f₁(r,z) and f₂(r,z) are defined which must satisfy equilibrium and the compatibility equation as follows:

. (6)

The linear elastic solution consists of finding functions f₁(r,z) and f₂(r,z) which satisfy Eq. (6) and boundary conditions (1)-(5). Then, displacements and stresses can be calculated using Eqs. (7)-(12) as follows

	(7)
	(8)
	(9)
	(10)
	(11)
	(12)

where:

and ; for i = 1, 2

and μ_i are the roots of the equation:

C₁₁C₄₄ μ² + { C₁₃(2C₄₄+C₁₃) - C₁₁C₃₃ }μ + C₃₃C₄₄ = 0

Other stress components and the circumferential displacements are zero due to axial symmetry.

B. Solution

To be able to satisfy the boundary conditions for the layer bonded to the rigid foundation, the functions f_i are assume to be:

(13)

where i = 1,2, A_i(λ) and B_i(λ) are unknown functions of O and J₀(λ) is the zeroth-order Bessel function of the first kind. In addition

Functions f_i defined in (13) are a variation of those used for the indentation of a transversely isotropic layer not bonded to the rigid foundation (Sakamoto et al., 1991).

Substitution of equation (13) into equations (3)-(5) leads to:

B2(λ) = - B1(λ),	(14)
	(15)
	(16)

in which A₂(λ), B₁(λ) and B₂(λ) are defined in terms of A₁(λ) which is still an unknown function. The main part of the procedure, which is quite tedious, consists of using the equations presented above to find function A₁(λ), which is expressed in terms of an infinite series of coefficients a_n and b_n. This procedure is thoroughly explained in Appendix B. The numerical procedure used to calculate the coefficients a_n and b_n given a set of materials properties is explained by Hara et al. (1990).

C. Parametric Analysis

A parametric analysis was undertaken assuming two sets of properties with different degree of anisotropy. In set 1 the ratio between the elasticity moduli was 20 and set 2 represents approximately an isotropic material, with a unit ratio between elasticity moduli (Table 1). For each set, displacements and stresses were calculated for several values of z/h, by Eqs. 8, 11, 27, 31, 39 and 37. In addition, the finite element program ALGOR (Algor Inc., Pittsburgh, USA) was used to compare some results of this analytic solution.

Table 1. Elastic properties used in the analysis

III. RESULTS

The a_n and b_n coefficients decrease rapidly (Tables 2 and 3) and function A₁(λ) can be accurately evaluated using the first seven terms of the series.

Table 2. Values of coefficients a_n and b_n for set 1.

Table 3. Values of coefficients a_n and b_n for set 2.

Nondimensional axial displacement ( = w/ε₀ ) for the second set is higher at the axis of symmetry and lower at a distance of 1.5 times the indenter radius (Fig. 2, 3).

Figure 2. Nondimensional displacements for spherical indenter

Figure 3. Nondimensional displacements for cylindrical indenter

Nondimensional axial stress ( = σ* r₀²/ P ) for the second set is higher at the axis of symmetry. The difference at the axis of symmetry is substantial, especially at the center of the layer, in which the stress for set 2 is approximately 2.5 times the stress for set 1 (Fig. 4, 5).

Figure 4. Nondimensional axial stress distribution for the spherical indenter

Figure 5. Nondimensional normal stresses for cylindrical indenter

Good correlation is observed between pressure distributions obtained with the Hertz's solution and this solution for material properties of Set 2, aspect ratio of 10 and a spherical indenter (Fig. 6).

Comparison of vertical displacements with those obtained with the finite element method also shows good agreement for the plane indenter (Fig. 7). Small differences are observed between the analytical and computational axial stress distributions (Fig. 8).

Figure 6. Pressure distributions under contact area of spherical indenter

Figure 7. Nondimensional axial displacements for cylindrical indenter, material Set 1

Figure 8. Nondimensional axial stresses for cylindrical drical indenter indenter, material Set 2

III. CONCLUSIONS

A new analytical linear elastic solution has been developed for the indentation, with a flat ended cylinder or a sphere, of a transversely isotropic layer bonded to a rigid foundation. The solution follows the procedure presented by Sakamoto et al. (1991) for the indentation with a flat ended cylinder of a transversely isotropic layer resting without friction on a rigid foundation and has been extended to consider the case of a spherical indenter. Unlike the solution developed by Sakamoto et al. (1991), which allows relative horizontal displacements at the rigid foundation, the solution presented here considers the layer completely bonded to the foundation, which is an appropriate assumption for the articular cartilage, firmly bonded to the subchondral bone.

This new solution can be employed to find relations between the indentation force and the elastic properties of the layer, which may be helpful to find material properties of articular cartilage under indentation tests. The solution also allows finding the complete distribution of stresses and displacements in the layer. These results could be used to validate finite element contact analysis of anisotropy layers.

This is a new solution and there are no reported results in the literature to compare those obtained in this study. For the nearly isotropic case, good correlation was observed between pressure distribution in the contact area obtained with this solution and that given by the Hertz. For an anisotropic material, vertical displacements were almost equal to those obtained with the finite element code Algor (Algor Inc., Pittsburgh, USA) and small differences were observed in the stress distributions. These differences can be understood since the precision of the finite element displacement method is reduced for stresses.

Results of the analysis showed smoother displacement and stress distributions for the more anisotropic layer, i.e. displacement and axial stress were lower below the contact area and descended more gradually towards the radial direction. Experimental studies (Wang et al., 2003) have shown that the degree of anisotropy of the articular cartilage is similar to that of set 1 of this analysis. This suggests that under equilibrium or impact loading, when the biphasic cartilage can be considered as an equivalent elastic layer, the anisotropy of the cartilage plays an important role to redistribute loads and to reduce stresses in the solid phase.

The solution presented here is valid for an elastic layer under infinitesimal deformation. Significant errors may be produced if this solution is applied to articular cartilage when large deformations and viscoelastic effects are significant.

Appendix A. Relation between flexibility coefficients and engineering constants

Appendix B. Detailed Procedure to Develop Stress and Displacement Equations.

In order to develop an expresion for A₁(λ) the following procedure must be undertaken. Substitution of Eqs. (14), (15) and (16) into Eqs. (1) and (2) leads to:

	(17)
	(18)

where y(λ) and f(λ) are defined as:

and

Next, contact stresses are assumed to be (Sakamoto et al., 1991):

(19)

where T_2n+1 are Tchebycheff polynomials and a_n' are unknown coefficients. This distribution contains a singularity equal to that of an isotropic layer indented by a cylindrical punch (Hayes et al., 1972). By using the Hankel inversion of Eq. (18) by means of Eq. (19) and the following identity of Bessel functions:

, (20)

where:

, n = 0,1,2, (21)

we obtain:

. (22)

Substitution of Eq. (22) into Eq. (17) leads to:

, (23)

where g(λ)=y(λ)/f(λ). Now, if the Gegenbauer's equation:

, (24)

(where ) is substituted into Eq. (23), it is obtained:

		, (25)
where		(26)

and δ_0m is the Kronecker's delta.

By matching the coefficients of cos(mθ) in both sides of Eq. (25) the equation (27) is obtained:

(27)

which represents an infinite system of simultaneous equations which are used to determinate the a_n coefficients. Then, the a_n' terms the function A₁(λ) are determined by Eqs. (26) and (22) respectively. Numerical solution of Eq. (27) was undertaken using the method of Hara, et al., (1990).

Contact stresses σ_z are calculated substituting the Eq. (26) into Eq. (19) as follows:

, (28)

where H(r₀ - r) is the Heaviside's step function.

In order to calculate w(r,z), Eq. (22) is sustituted into Eq. (17).

. (29)

The above integrals were calculated following the procedure developed by Sakamoto et al. (1991), which consists of calculating the limits when λ tends to infinite of λZ_n(λ) and g(λ), and rearranging the equation (29) as shown bellow.

, (30)

where

. (31)

Substitution of Eqs. (30), (31) into Eq. (29) leads to Eq. (32):

(32)

So far, the above equations are valid for a cylindrical indenter. To find equivalent equations for spherical indenter we assume that:

, (33)

where the a_n coefficients are the same of equation (27) and the b_n are unknown coefficients which can be found as follows.

Substitution of Eq. (32) into Eq. (22) leads to:

. (34)

Next, substitution of Eqs. (34), (24) and f(r) = r²/2R into Eq. (17) yields to Eq. (35):

. (35)

Next, by matching the coefficients of cos(mθ) on both sides of last equation we obtain:

, (36)

from which the b_n coefficients can be calculated as explained before for equation (27). In order to obtain contact stresses and displacements, Eq. (33) is substituted into Eq. (19) and the contact stress at r = r₀ is made equal zero, it is obtained the following entity

. (37)

Using Eqs. (37) and (33), the contact stresses and dispacements can be written as follows:

	(38)
	(39)

Then, using Eqs. (30) and (31) we can write the Eq. (39) as:

(40)

With the above equations the analytical solution for spherical indenter is completed. Now, the indentation force is obtained for both indenters by integrating (Eq. (28) and (38))

The indentation forces for the cylindrical and spherical indenters are respectively:

	(41)
	(42)

Acknowledgment
The authors are grateful to Colciencias to provide support to undertake this project.

REFERENCES
1. Cohen, B., T.R. Gradner and G.A. Ateshian, "The influence of transverse isotropy on cartilage indentation behavior. A study on the human humeral head", Trans. Aorthop. Res. Soc., 18, 185 (1993).         [ Links ]
2. Ewers, B.J., W.N. Nedwberry, J.J. García and R. Haut, "Alterations of the mechanical properties of bone underlying cartilage in a traumatized joint", Proc. 44^th Annual Meeting Orthopaedic Research Society, (1998).         [ Links ]
3. García, J.J., N.J. Altiero and R.C. Haut, "An approach for the stress analysis of transversely isotropic bi-phasic cartilage under impact load", Journal of Biomechanical Engineering, 120, 608-613, (1998).         [ Links ]
4. García, J.J., N.J. Altiero and R.C. Haut, "Estimation of in situ elastic properties of biphasic cartilage based on a transversely isotropic hypo-elastic model", Journal of Biomechanical Engineering, 122, 1-8, (2000).         [ Links ]
5. Hara, T., T. Akiyama, T. Shibuya and T. Koisumi, "An axisymmetric contact problem of an elastic layer on a rigid base with a cylindrical hole", JSME Int. Journal, 33, 461, (1990).         [ Links ]
6. Hayes W.C., L.M. Keer, G. Herrmann and L.F. Mockros, "A mathematical analysis for indentation tests of articular cartilage", J. Biomechanics, 5, 541-551, (1972).         [ Links ]
7. Lekhnitskii, S.G.: Theory of elasticity of an anisotropic body, MIR Publishers, Moscow (1981).         [ Links ]
8. MacKenzie, E.J., J.H. Siegel and S. Shapiro, "Functional recovery and medical cost of trauma: an analysis by the type and severity of injury", Journal of Trauma, 28, 281-298, (1988).         [ Links ]
9. Mow, V.C., S.C. Kuei, W.M. Lai and C.G. Armstrong, "Biphasic creep and stress relaxation of articular cartilage: Theory and experiment", Journal of Biomechanics, 102, 73-84, (1980).         [ Links ]
10. Mow, V.C., M.C. Gibbs, W.M. Lai, W.B. Zhu and K.A. Athanasiou, "Biphasic indentation of articular cartilage. II. A numerical algorithm and an experimental study", Journal of Biomechanics, 22, 853-861, (1989).         [ Links ]
11. Newberry, W.N., J.J. García, C.D. Mackenzie, C.E. Decamp and R.C. Haut, "Analysis of acute mechanical insult in an animal model of posttraumatic osteoarthrosis", Journal of Biomechanical Engineering, 120, 704-709, (1998).         [ Links ]
12. Sakamoto, M., T. Hara, T. Shibuya and T. Koizumi, "Indentation by a circular rigid punch of an transversely isotropic layer on a rigid foundation", JSME international Journal, Serie I, 34, 130-134, (1991).         [ Links ]
13. Töiräs, J., T. Lyyra-Laitinen, R. Niinimäki, M.T. Lindgren, M.T. Nieminen, I. Kirivanta and J.S. Jurvelin, "Estimation of the Young's modulus of articular cartilage using an arthroscopic indentation instrument and ultrasonic measurement of tissue tichness", Journal of Biomechanics, 34, 251-256, (2001).         [ Links ]
14. Wang, C.C-B., N.O. Chahine., C.T. Hung and G.A. Ateshian, "Optical determination of anisotropic properties of bovine articular cartilage in compression", Journal of Biomechanics, 36, 339-353, (2003).         [ Links ]