Biological Tissues Using Nanoindentation

Published Date: 02 Nov 2017

Introduction

In order to understand how different types of biological tissues perform their functions inside the human body or to develop artificial materials that could replace them it is important to know their mechanical properties. Biological tissues need to have certain mechanical properties in order to fulfill their role inside the human body. Most obvious example are bones. Main role of bones is to support the weight of the body. Diseases like osteogenesis imperfecta that influences the production of collagen, result in inferior mechanical properties of bones making them very brittle (Wikipedia-Osteogenesis imperfecta). Osteoarthritis is a progressive joint disease in which articular cartilage inside the joint degenerates, which results in inflammation, painful and limited motion of joints (Wikipedia-Osteoarthritis)(Lopez, Amrami et al. 2008). Articular cartilage is a very complex material composed of 65-80% of water and the rest 20-35% are chondrocytes, proteoglycans and collagen type II. The arrangement and density of these constituents varies throughout the thickness of articular cartilage (Lopez, Amrami et al. 2008). In order to understand the causes and progression of osteoarthritis it is necessary to understand mechanical properties of healthy and diseased articular cartilage.

There are many different types of mechanical tests that are used to measure certain mechanical properties, like tensile test, compression test, indentation, nanoindentation. Most test procedures, like tensile or compression tests, require special preparation of material samples. In order to test material samples using compression or tensile test, the material sample needs to be a certain size in order to attach it to the testing equipment. This can pose a problem when dealing with material samples like articular cartilage which is rather thin and has several layers with different mechanical properties.

Depending on the type of material response, different mechanical properties are measured. When dealing with elastic materials (most biological tissues exhibit some type of elasticity) Young’s elastic modulus is one of the most important mechanical properties to be measured (Fischer-Cripps and C. 2004; Ebenstein and Pruitt 2006) .

Nanoindentation is an (almost) noninvasive method for obtaining material properties of the specimen material from experimental readings of indenter load (P) and depth of penetration (h). Forces are usually in milinewton range with the resolution of a few nanonewtons. Depths of penetration are in the order of micrometers with a resolution of less than a nanometer(Fischer-Cripps and C. 2004). In most cases the material properties that are obtained by nanoindentation are elastic modulus (E) and hardness (H) (Oliver and Pharr 1992; Fischer-Cripps and C. 2004). Main advantage of this method is that no extensive preparation of the specimen is required, testing procedure is automated and equipment is easy to use (Oyen and Cook 2009).

Although equipment is easy to use, resulting data is not easily interpreted. Interpretation of data depends on the mechanical response of material specimen to nanoindentation (linear or non-linear response, time dependency). This is especially true for biological materials which usually exhibit significant time-dependent properties (viscoelastic or poroelastic) (Oyen and Cook 2009), since the standard analyses are derived for indentation of elasto-plastic materials which are inappropriate for viscoelastic or viscoelastic-plastic materials (Oyen 2006). Although the nanoindentation instruments have significantly evolved during the past decades, such tests still require significant skill in order to obtain useful data and interpret it in a correct way.

Nanoindentation tests are extremely sensitive to thermal expansion, mechanical vibration and acoustical noise. In order to minimize the impact of thermal expansion, specimen and indenter need to be in thermal equilibrium (Fischer-Cripps and C. 2004). Indenter is mounted on a shaft that needs to be very stiff and lightweight in order to minimize compliance. Indenter must be connected to the indenter shaft firmly and have a minimal compliance in order to get accurate readings. That is why indenters are made out of very hard materials, such as diamond. The downside of using diamond for making indenters is that it is very brittle and can easily be chipped. Indenters must be absolutely clean and free from any contaminants. Diamond indenters are cleaned by pressing them into a block of dense polystyrene.

There are two types of controlling the nanoindentation test equipment, load control and displacement control. In load controlled machines, the load imposed on the specimen is the input variable and the displacement at a certain load is being measured. In displacement controlled machines the input variable is displacement and the reaction force of the material sample is being measured. Most nanoindentation devices are load-controlled(Fischer-Cripps and C. 2004). The choice between load and displacement control is dependent mostly on the type of material that is being tested. When fractures occur in the material sample during nanoindentation testing, there will be a change in the slope of the P-h curve when load-controlled machines are used. With displacement controlled machines there will be no change in the slope of the P-h curve, so the onset of fractures cannot be detected with displacement controlled machines (Oyen and Cook 2009). When testing viscoelastic or poroelastic materials two types of test are usually being performed: creep and load relaxation tests. In creep tests the load P(t) is constant and the displacement h(t) changes with time. This requires the use of load controlled machines, because the constant value of load must be inputted. For load relaxation test, the displacement h(t) is constant and the load P(t) changes with time. In this case displacement controlled machines are more suited (Cheng, Xia et al. 1999).

Two main approaches for calculating material properties from nanoindentation tests are analytical and computational. Analytical methods aim to derive analytical expressions for determining mechanical properties of materials from load-displacement data obtained by nanoindentation tests. Analytical techniques are mostly used for simpler material responses like linear elasticity, hyperelasticity although there are analytical techniques developed for viscoelastic and poroelastic materials. Downside of analytical techniques is that sometimes analytical calculations are very complex and numerical methods must be employed in order to find the solution. Computational methods are used when there are no analytical solutions for certain indentation problems like the evolution of the plastic zone beneath the indenter (Fischer-Cripps and C. 2004) or for some advanced poroelastic theories like Fiber reinforced poroelastic model (FRPE) (Gupta, Lin et al. 2009). Finite elements method is the most often employed computational method. Computational methods are especially convenient for finding solutions to non-linear problems, since such problems are usually solved using iterative procedures(Fischer-Cripps and C. 2004).

The main goal of this literature review is to gather analytical solutions from different sources for nanoindentation tests of different material models that could be used for determining material properties of biological tissues. The paper is divided into several chapters. First chapter briefly explains different nanoindenter tip geometries. Subsequent chapters deal with different material models and analytical solutions for obtaining material properties from load-displacement data. First material model that will be described is the simplest one, linear elasticity. Next chapter deals with elastoplasticity, since in most nanoindentation tests elasticity is coupled with plasticity, which is reflected on the way the load-displacement data should be interpreted. Next is the hyperelastic material model, developed from nanoindentation tests performed on elastomers. After hyperelasticity comes the chapter on viscoelasticity, a material model suited for soft hydrated tissues. The last chapter deals with poroelasticity, which is a material model very similar to viscoelasticity, but more complex in its aim to explain the physical characteristics of materials that lead to phenomena like creep and load relaxation.

Nanoindenter tip geometry

There are several types of nanoindenter tip geometries being used in nanoindentation testing: cylindrical (flat) tips, spherical tips and conical (pyramidal) tips. Different tip geometries are used for studying different types of materials. Microstructure of the material will also dictate the size of the tip, since it has to be larger than the size of smallest constituent of the material (like the size of the cell, fiber or crystal that a certain material is made of) (Ebenstein and Pruitt 2006). For indentation of soft tissues, hydratation state of the material sample is crucial for obtaining correct results because most tissues exhibit poroelastic behavior which is greatly dependent on the amount of fluid absorbed by the tissue (Ebenstein and Pruitt 2006; Kurland, Drira et al. 2012).

Main advantage of cylindrical (flat) tips is that the resulting unloading curve remains linear because the cross-section area of flat indenter tips does not change with indentation depth. This also means that the contact area does not depend on the penetration depth because it can be easily calculated from the radius of the indenter (Fischer-Cripps and C. 2004). Main problem with flat (cylindrical) indenters is that the load is concentrated on the outer edge of the indenter (Choi and Shield 1981; Ebenstein and Pruitt 2006).

Spherical indenters are often used for indentation of compliant materials (soft tissues or polymers) because they minimize plastic deformation and stress concentrations (Ebenstein and Pruitt 2006). Unlike flat indenters, the area of contact changes with depth for spherical indenters, but it can be easily calculated from indentation depth and radius of the indenter (Fischer-Cripps and C. 2004).

Several types of pyramidal tips are in use today: Berkovich, Vickers and Knoop. Main difference between these indenters is the face angle (or half angle) of the tip. Berkovich indenter is a three-sided pyramid indenter with a half angle of 65.27° (Fischer-Cripps and C. 2004). The original Berkovich indenter had a face angle of 65.0333° in order to have the same area to depth ratio as the Vickers indenter. Another important geometry parameter is the equivalent cone angle. This is the angle of the tip of a conical indenter for which the conical indenter will have the same area to depth ratio as the chosen pyramidal indenter. For Berkovich indenter the equivalent cone angle is 70.296° (Fischer-Cripps and C. 2004). Vickers indenter is a four sided pyramid with a face angle of 68° and the equivalent cone angle is the same as for the Berkovich indenter, 70.296° (Fischer-Cripps and C. 2004). Cube corner indenters are more acute then Vickers and Berkovich. The face angle is 35.26° and the equivalent cone angle is 42.278° (Fischer-Cripps and C. 2004). The Knoop indenter is similar to Vickers, but has unequal lengths of edges which results in an impression with uneven diagonals (longer diagonal is approximately seven times longer than the shorter one). The face angles are 86.25° and 65°. The equivalent cone angle for Knoop indenter is 77.64° (Knoop, Peters et al. 1939; Fischer-Cripps and C. 2004).

Knowledge of the exact indenter tip geometry is crucial, because it is very hard to determine the surface area of the impression left by the indenter on the surface of material specimen, so the surface area is calculated using the depth of impression which can be established from the load-displacement data (P-h curve) (Oliver and Pharr 1992; Fischer-Cripps and C. 2004). Surface area of the impression on the material specimen is important because it is used for calculating material properties such as hardness (H), stiffness (S) or Young’s modulus (E). Indenter area function (or shape function) shows the cross-sectional area of the indenter as a function of the distance from the tip of the indenter (Oliver and Pharr 1992; Oliver and Pharr 2004).

Depending on the type of material we wish to examine, a suitable indenter tip must be chosen. Different types of indenter tip geometries are suitable for studying different types of material responses (linearly elastic, visco-elastic, elasto-plastic). The shape of load-displacement curve depends not only on the shape of the loading function, but also on the choice of indenter tip (Oyen and Cook 2009).

In most cases the displacement is limited to around 10% of the material sample thickness in order to avoid the influence of the substrate material sample is resting on (Doerner and Nix 1986; Fischer-Cripps and C. 2004; Chen and Diebels 2012). Chen and Diebels (Chen and Diebels 2012) also report that if the stiffness of the substrate is much larger than the stiffness of the material sample (100 times stiffer), there will be no influence of the substrate on the force displacement data.

Analysis of linear elastic materials

Linear elasticity is the simplest type of material response. The relation between stress and strain is governed by a constant elastic modulus (E), which is the most important material property that needs to be obtained from nanoindentation tests. Elastic deformation is time independent, rate independent, it has no threshold and it is reversible (Oyen and Cook 2009). Analysis of nanoindentation data is based on Hertz theory of elastic contact. This theory relates force to depth of penetration (displacement) for different combinations of indenter and material geometries. For the purpose of nanoindentation tests, the contact is always between a certain geometry of the indenter and a flat surface (material specimen is usually flat ) (Fischer-Cripps and C. 2004). The slope of force-displacement curve is called stiffness and is defined as (Oliver and Pharr 1992):

Cylindrical indenters

Based on Hertz theory of elastic contact, the relation between the force and the depth of penetration:

(1)

where E* is the combined elastic modulus of the indenter and material specimen, a is the radius of cylinder and h is the depth of penetration and P is the loading force. The combined elastic modulus E* is defined as (Fischer-Cripps and C. 2004):

ν’ and E’ apply to the indenter properties and ν and E apply to the material specimen.

By differentiation of force over elastic displacement, the slope of loading curve is obtained:

(2)

It can be seen from equation (2) that for cylindrical indenters, elastic modulus can be obtained from the slope of the loading curve using only the radius of the indenter. The elastic modulus can also be calculated using the contact area (A) (Fischer-Cripps and C. 2004):

(3)

Conical indenters

Analytical solution for elastic contact between a conical indenter and elastic half-space was derived by Hertz (Fischer-Cripps and C. 2004):

(4)

where a is the radius of circle of contact, E is the Young’s modulus and α is the cone half- angle.

Spherical indenters

The force - displacement relationship is derived from Hertz equation for elastic contact between a sphere and elastic half space:

(5)

Where he is the elastic deformation, R is the radius of the indenter and E* is the combined elastic modulus.

Using eq. (5) we arrive to the expression for stiffness (slope of the loading curve) (Fischer-Cripps and C. 2004):

(6)

Analysis of elasto-plastic deformations

In real nanoindentation tests it is often hard to produce only elastic deformation. Usually, elastic deformations are coupled with plastic deformation, so the resulting load-displacement curve displays elasto-plastic material response. In order to calculate the Young’s modulus from such load-displacement curve, only the slope of the unloading curve must be used. This is because during loading a combined elasto-plastic deformation occurs, while unloading is purely elastic (D.Tabor 1951; Oliver and Pharr 1992). Analytical solutions are the same as for the purely elastic case, the only difference is that the unloading part of the force-displacement curve must be used.

Except for the cylindrical indenters, the unloading curve for other types of indenters is not linear because the surface area of contact changes with depth (Doerner and Nix 1986; Oliver and Pharr 1992). Most authors consider that the initial part of unloading curve is linear(Doerner and Nix 1986; Fischer-Cripps and C. 2004), although some authors disagree with that (Oliver and Pharr 1992). Doerner and Nix (Doerner and Nix 1986) base their claims on the fact that in the initial part of unloading the contact area between the material sample and the indenter remains constant for all types of indenter geometries, so the initial part of the unloading curve must be linear. Oliver and Pharr (Oliver and Pharr 1992) conducted a number of indentation experiments and found that the initial portion of the unloading curve is almost never linear except for the case of cylindrical indenters. Nevertheless, the initial slope of the unloading curve is used to calculate the elastic modulus and is considered to be linear.

Cylindrical indenters

Depth of penetration is made of two parts: hr which is the depth of residual impression (plastic deformation) and he which is the elastic displacement during unloading. By differentiation of force over elastic displacement, the slope of unloading curve is obtained:

(2)

The elastic modulus can also be calculated using the contact area in the same way as for the purely elastic indentation (Fischer-Cripps and C. 2004):

(3)

Figure 1 Load-displacement curve for cylindrical indenters (Fischer-Cripps and C. 2004)

Conical indenters

Pyramidal indenters are considered to be conical since the pyramidal geometry has a small effect on the final result (Pharr, Oliver et al. 1992; Fischer-Cripps and C. 2004). The initial part of the unloading curve made with Berkovich pyramidal indenter is linear for a wide range of materials (Doerner and Nix 1986). This means that equations for cylindrical indenters can be used to calculate the surface of contact area for conical indenters from depth measurements. When the material specimen is being unloaded using a conical indenter, the contact radius remains constant until the surface of the material specimen stops conforming to the shape of the indenter. This results in linear initial part of the unloading curve.

The derivation of analytical solution starts with expression for contact between a conical indenter and elastic half space derived by Hertz (Fischer-Cripps and C. 2004):

(4)

where a is the radius of circle of contact, E is the Young’s modulus and α is the cone semi- angle.

Figure 2. shows the P-h curve for conical indenters. It can be seen that the initial part of the unloading curve is linear, when the contact area between the indenter and material specimen is constant. If the indenter was cylindrical, the unloading curve would follow the BD path. In reality when using conical indenters, unloading curve follows the BC path which is non-linear. Using the P-h curve and the linear unloading approximation, the depth of residual impression (hrc ) can be obtained:

(7)

For Vickers and Berkovich indenters, the relationship of projected area Ap of indentation and the depth of residual impression hc is:

(8)

Using the equation for the slope of unloading curve for cylindrical indenters and equation (6), we arrive to the expression for the slope of the initial part of unloading curve for conical indenters (Fischer-Cripps and C. 2004):

(9)

Figure 2 Load-displacement curve for conical indenters (Fischer-Cripps and C. 2004)

The cylindrical method is a good approximation for conical indenters only when the initial part of unloading curve is linear and hrc<<hmax. This is valid for materials that do not show a significant elastic recovery (E/Y ratio is large, where E is the elastic modulus and Y is the yield stress). For highly elastic materials where hrc>>hmax, initial part of the unloading curve is also non-linear, cylindrical method approximation will lead to large errors (Fischer-Cripps and C. 2004). In order to solve the problem of non-linear unloading curve, Oliver and Pharr (Oliver and Pharr 1992; Oliver and Pharr 2004) started with load-displacement relationship for simple punch geometries derived by Sneddon (I.Sneddon 1965):

(10)

Where α and m are constants related to the geometry of the indenter, h is the maximum depth and hf is the depth of residual plastic impression (Figure 3). Exponent m has a value of m=1 for flat cylindrical indenters, m=2 for conical indenters and m=1.5 for spherical indenters.

Figure 3 Indentation depths(Oliver and Pharr 1992)

Based on eq. (10), Oliver and Pharr derived a general equation for the slope of the P-h curve (elastic unloading stiffness) that applies to any axisymmetric indenter and takes into account non-linear property of the elastic unloading curve (Oliver and Pharr 2004).

(11)

Correct value of correction factor β is essential in order to obtain correct value of elastic modulus E (and hardness H when plastic material response is measured). For small deformations of an elastic material made by an axisymmetric indenter, the value of correction factor β is equal to 1. Small deformations can be achieved by using indenters with a half- angle close to 90°, which means that in most real cases the small deformation criterion will not be met. King (R.B.King 1987) used numerical methods to obtain the value of correction factor β for flat-ended indenters. He reported that for square based indenters β=1.012 and for triangular based indenters β=1.034. For instrumented indentations using Berkovich indenter (triangular pyramid) value of β=1.034 is used. Using a more precise method, different value of the correction factor β was found by Vlassak and Nix (Vlassak and Nix 1994) for flat-ended triangular indenter, β=1.058. Using simple elastic analysis procedures and assuming that pressure profile is perfectly flat, Hendrix found the values of correction factor β for Vickers indenter β=1.0055 and for Berkovich indenter β=1.0226 (Hendrix 1995). Using a full 3D finite element analysis for purely elastic materials, Larsson et al. (Larsson, Giannakopoulos et al. 1996) found that β is dependent on Poisson’s ratio:

(12)

Assuming that ν=0.3, the correction factor β=1.14. This value of correction factor β applies to indentation of flat elastic half-space and does not take into account distortion of the surface caused by plastic deformations.

Cheng and Cheng (Cheng and Cheng 1998; Cheng and Cheng 2004) calculated the values of β for a wide range of materials, both work-hardening and non-hardening, using finite elements method. They found that β was independent of both the Poisson ratio and E/Y (elastic modulus and yield strength ratio). For non-hardening materials they found the value of β=1.05 and for work-hardening materials β=1.085. There are slight inconsistencies in the results in the two papers by Cheng and Cheng, so the question remains which values are correct or are they correct at all (Oliver and Pharr 2004). Oliver and Pharr (Oliver and Pharr 2004) found that the value found by King (R.B.King 1987) β=1.034 produces good results with errors not larger than 6.5% in contact area and hardness. Oliver and Pharr (Oliver and Pharr 2004) concluded that the correction factor β falls into the range

and choose β=1.05 as an optimal choice with a potential for error

.

Spherical indenters

As mentioned before, Oliver and Pharr (Oliver and Pharr 1992; Oliver and Pharr 2004) noted that the unloading response for most materials is non-linear which led them to use a power law relationship of a cone rather than the linear load-displacement relationship of a flat cylindrical indenter (Fischer-Cripps and C. 2004). Field and Swain (Field and Swain 1993) applied elastic equations directly to the unloading force-displacement curve, not to the slope of a fit to the data for spherical indenter, but this method can be used for other types of indenters as well (Fischer-Cripps and C. 2004). The relationship between the force and displacement is derived from Hertz equation for elastic contact between a sphere and elastic half space:

(5)

where:

hmax is the maximal depth, while hr is the depth of the residual impression (depth of plastic deformation).

Using eq. (5) we arrive to the expression for stiffness (slope of the unloading curve) (Fischer-Cripps and C. 2004):

(6)

Figure 4 Load-displacement curve for spherical indenters(Fischer-Cripps and C. 2004)

Analysis of hyperelastic materials

Hyperelastic materials have a non-linear relationship between stress and strain. In linear elastic materials stress and strain are related through Young’s modulus E, while in hyperelastic materials they are related through strain energy function W (Giannakopoulos and Panagiotopoulos 2009; Chen and Diebels 2012).

There are many types of strain energy functions, but most researchers are concentrated on neo-Hooken and Mooney-Rivlin strain energy functions because of their simplicity (Chen and Diebels 2012). Neo-Hookean strain energy function is used to describe hyperelastic material properties at low strains (Marckmann and Verron), while Mooney-Rivlin strain energy function is used to describe hyperelastic material properties at moderate strain (Marckmann and Verron ; Giannakopoulos and Triantafyllou 2007). Neo-Hookean and Mooney-Rivlin models are too simple to describe the material behavior at large strains, so some authors looked into other hyperelastic models, such as Yeoh model (Chen and Diebels 2012). Neveretheless, Neo-Hookean and Mooney-Rivlin models proved to be adequate for use in nanoindentation tests, since strains remain small during nanoindentation testing (Busfield and Thomas 1999).

Cylindrical indenters

Main advantage of cylindrical indenters is the linear relation between load and displacement which makes it easy to determine the slope of the curve without knowing the initial point of contact (Gent and Yeoh 2006).

Choi and Shield (Choi and Shield 1981) derived an analytical expression for calculating the shear modulus from force-displacement data. Their formula shows the relation between the force and displacement for compressible and incompressible materials. The solution is made out of two parts, first order and second order load. For compressible materials it is not possible to find a solution because of singularities that occur at the edge of the cylinder, so the expression is only made out of first order load:

(13)

where G is the shear modulus, d is vertical displacement, α is the radius of the impression (which is equal to the radius of the cylindrical indenter) and ν is the Poisson ratio.

For incompressible materials the singularity at the edge of the indenter disappears and the expression for the force is made out of first and second order:

(14)

Conical indenters

Giannakopoulos and Panagiotopoulos derived relations between the P-D (force-displacement) response and mechanical properties for cone indenters (Giannakopoulos and Panagiotopoulos 2009). They used a non-linear analysis of the second order to derive closed form results that gives the relation between indentation depth, contact radius, half angle of the cone and the applied load P.

It should be emphasized that Giannakopoulos and Panagiotopoulos use the angle φ between the cone and material sample and the cone half angle is defined as (Figure 5):

Figure 5 Cone angle as defined by Giannakopoulos and Panagiotopoulos (Giannakopoulos and Panagiotopoulos 2009)

Although adhesion and friction might be important at certain cone angles, they were ignored in order to simplify the solution. The material is also considered to be incompressible (Poisson ratio is 0.5). In this paper three types of hyperelastic strain functions were analyzed: Mooney-Rivlin, Gent and one term Ogden.

The relation between the force and displacement was defined as:

(15)

(16)

where P is the force, D is the displacement, E is the Young’s modulus, G is the shear modulus and φ is the cone half angle.

Mooney-Rivlin strain function is defined as:

(17)

Where I1 and I2 are the invariants of the left Cauchy deformation tensor and c1 and c2 are constants, c1+c2>0. If c2=0, then the Mooney-Rivlin model reduces to Neo-Hookean.

It is not possible to calculate the value of c1 and c2 independently, but their sum can be obtained from the shear modulus at zero strain using the expression:

G=2(c1+c2) (18)

Wagner (Wagner 1994) showed that c1=c2 for small values of c1 (less than 0.1 MPa for natural rubber).

Using the finite element analysis, Giannakopoulos and Panagiotopoulos (Giannakopoulos and Panagiotopoulos 2009) concluded that the results obtained for all types of Mooney-Rivlin materials can be obtained by liner elastic analysis with more than 95% accuracy, but only when blunt cones are used. They also showed that the c2 constant has to be larger than zero, otherwise the indentation problem cannot be solved. As it was mentioned before, Bushfield and Thomas (Busfield and Thomas 1999) reported that this unusually good prediction of the linear elastic theory was the result of small strains during nanoindentation experiments. Giannakopoulos and Panagiotopoulos (Giannakopoulos and Panagiotopoulos 2009), however, report that the strains were of the order of 50%, in the range of validity of Mooney-Rivlin model. The answer lies in the cone angle, since other authors also report that sharper cone angles produce larger strains (Fischer-Cripps and C. 2004). It is not clear if these results are valid for pyramidal indenters such as Vickers, Berkovich and Knoop, but Giannakopoulos and Panagiotopoulos proved they are valid for blunt pyramids. Fischer and Cripps (Fischer-Cripps and C. 2004) also reported that the results obtained using cone and pyramid indenters could be analyzed using the same formulas because the ratio of the length of the diagonal or the radius of the circle of contact to the depth of indentation remains constant for increasing indenter load.

Spherical indenters

For contact analysis between a rigid spherical indenter and an incompressible and isotropic hyperelastic material, a general solution can be provided only in the case when the indentation is not very deep, if the ratio of the indentation radius and the radius of the indenter is very small (Giannakopoulos and Triantafyllou 2007):

There are a few papers that tried to provide analytical solutions to contact problems between a spherical indenter and a hyperelastic material sample. First one was done by Choi and Shield (Choi and Shield 1981). For incompressible rubber materials they derived the following terms for predicting the load P and the depth of indentation D:

(19)

(20)

where G is the shear modulus at zero strain, α is the contact radius (radius of the impression on the material sample) and R is the indenter radius. These equations predict linear elastic response for the contact problem of spherical indenter and hyperelastic material.

Second analysis was done by Sabin and Kaloni (Sabin and Kaloni 1983) who obtained the following equations for calculating the force P and the depth of indentation D:

(21)

(22)

Giannakopoulos and Triantafyllou started their analysis with strain energy function suggested by Murnaghan (Murnaghan 1951):

(23)

where α1, α2, α3, α4, α5 are constants and J1,J2,J3 are invariants of the Lagrange strain tensor.

In case of uniaxial response, eq. (22) reads:

(24)

Giannakopoulos and Triantafyllou declare that eq. (24) does not represent real material response, but a fictive response that is useful in indentation analysis. Using eq. (24), the following expressions for the force P and depth of penetration are derived:

(25)

(26)

It can be seen that eq. (26) is the same as eq. (22) derived by Sabin and Kaloni (Sabin and Kaloni 1983), but eq. (25) is different than eq. (21). The P-D response derived by Giannakopoulos and Triantafyllou results in a stiffer response than linear elasticity and Sabin and Kaloni solution, which can be seen on Figure 6.

Giannakopoulos and Triantafyllou (Giannakopoulos and Triantafyllou 2007) also derived force and displacement equations for Mooney-Rivlin hyperelastic response:

(27)

Constants c1 and c2 can be calculated from shear modulus at zero strain:

(28)

In order to avoid body forces, c1=c2. In this case the Mooney-Rivlin strain function reads:

(29)

In this case the second order elasticity solution is the same as linear incompressible elasticity solution:

(30)

(31)

Figure 6 Experimental results compared with linear elasticity and second order elasticity derived by Giannakopoulos and Triantafyllou (G&T) and Sabin and Kaloni (S&K)(Giannakopoulos and Triantafyllou 2007)

The main issue with the reverse problem of estimating material properties from experimental data is that constants c1 and c2 can be calculated only as a sum from zero strain shear modulus. This makes it impossible to capture the exact form of the Mooney-Rivlin strain function from load-displacement curve, except for the case c1=c2 (Giannakopoulos and Triantafyllou 2007).

Analysis of linearly viscoelastic materials

Viscous deformation is time dependent, irreversible and does not have an onset threshold. Viscous flow changes the apparent value of elastic modulus as a function of experimental time or frequency and is often considered together with elastic deformation. Most important mechanical properties that are used for describing viscoelastic material behavior are elastic modulus E, viscosity η and time constant τ (Oyen and Cook 2009). Time constant can be defined in two ways: empirical viscoelastic and poroelastic view. In viscoelastic view the time constant is a function of viscosity η and elastic modulus E (Oyen and Cook 2009):

In poroelastic view, time constant is defined by the characteristic length scale l and hydraulic permeability coefficient k from Darcy’s law (Oyen and Cook 2009):

Viscoelastic material response is time-dependent and comprised of instantaneous elasticity, delayed elasticity and viscous flow (Cheng, Xia et al. 2005).

There are several mathematical models that describe viscoelasticity: Maxwell model, Kelvin-Voight and Standard model. Most researchers use the Standard model to describe viscoelastic response of polymers since this model provides sufficiently accurate results without being too complicated (Cheng, Xia et al. 1999; Sakai 2002; Cheng, Xia et al. 2005; Oyen 2006). Standard model is made out of a spring and a dashpot connected in parallel and an additional spring connected in series with the parallel spring and dashpot.

In the correspondence approach, which is the most common way of analyzing viscoelastic material behavior in nanoindentation, viscoelastic operators are substituted for elastic constants (E.H.Lee and J.R.M.Radok 1960). Main goal of this method is to remove the effects of time-dependent material behavior on elastic modulus measurements (Oyen 2006). This is done by removing the time variable in the governing equations using Laplace transformations with respect to time (Cheng, Xia et al. 2005).

For measuring different responses of viscoelastic material (creep or load relaxation) different nanoindentaion control techniques must be used (load or displacement control). For measuring creep response, load control must be utilized since creep occurs when the applied load is constant (Mattice, Lau et al. 2006). For measuring load relaxation response, displacement control must be utilized since load relaxation occurs when displacement is constant (Cheng, Xia et al. 1999; Cheng, Xia et al. 2005).

When dealing with viscoelasticity, loading conditions are very important for validity of measurement results. Most researchers use step functions in load controlled indentation studies because they are analytically more convenient, but are impossible to implement in real experiments. Ramp correction factor is utilized in order to cross the bridge between the analytically convenient step loading and ramp loading used in real nanoindentation experiments (Oyen 2005).

Cylindrical indenters

Sakai (Sakai 2002) derived the solution for viscoelastic indentation using the standard linear solid model that consists of one Maxwell element (serial combination of a spring with modulus E0 and a dashpot with viscosity η) and a parallel spring with equilibrium modulus Ee. Sakai used a step loading function (Heaviside step function) and did not consider the ramp correction factor. The derivation started with the load-displacement relation for linearly elastic materials derived by Sneddon (I.Sneddon 1965):

(32)

where n=1 for flat indenters, E is the elastic modulus and C is given by:

(33)

where D is diameter of the indenter and ν is the Poisson’s ratio.

Using the Boltzmann hereditary integral, the following relation between the load and displacement has been derived:

(34)

where 0 ≤ t’ ≥ t, Y(t-t’) is the relaxation function and C is the factor defined in eq. (33) for flat cylindrical indenters .

The relaxation function is equal to:

(35)

where E0 is the elastic modulus responsible for initial elastic response (lower slope in Figure 7), while Ee is the long term elastic modulus responsible for the secondary elastic response (upper slope in Figure 7), t is the complete time of the indentation test and Ï„ is the relaxation time.

Figure 7 Normalized indentation load versus time (or depth) for a flat cylindrical indenter(Sakai 2002)

Sakai assumed a constant rate of penetration v0 which results in:

(36)

Using eq. (35) and eq. (36) and putting them in eq. (34), Sakai arrived to the final expression for indentation load:

(37)

where is the normalized penetration depth:

(38)

Cheng et al. (Cheng, Xia et al. 1999) derived their solution for viscoelastic nanoindentation using a three element model, where a Kelvin-Voight model with a parallel spring (E2, ν2) and a dashpot (η) is connected in series with another spring (E1, ν1), Figure 8. Three different loading conditions were assumed: a single step function (Heaviside step function), multiple step function and an arbitrary loading function. Cheng et al. also derived solutions for load control (creep test) and displacement control (relaxation test), unlike Sakai (Sakai 2002) who derived solutions only for displacement control. Nevertheless, both the creep and relaxation tests should yield the same solution. Only the solution for an arbitrary load and displacement histories will be presented here since that is the most general solution and there have been some concerns regarding the validity of using a Heavisde step function since it cannot be applied to real nanoindentation experiments (Oyen 2005).

Figure 8 Three element viscoelastic model (Cheng, Xia et al. 2005)

Using Boltzmann superposition principle, Cheng et al. (Cheng, Xia et al. 1999) derived the following expression for relaxation response:

(39)

where R is the radius of the indenter, Gr(t) is the relaxation modulus,

(40)

(41)

(42)

(43)

(44)

(45)

(46)

(47)

(48)

(49)

Where K is the bulk modulus and G is the shear modulus defined as:

(50)

(51)

And creep response:

(52)

where R is the radius of the indenter, Jc is the creep compliance

(53)

(54)

(55)

(56)

(57)

(58)

(59)

Conical indenters

Sakai (Sakai 2002) derived the solution for nanoindentation of viscoelastic materials using conical indenters in the same way as for the cylindrical indenters. The load P(t) is defined in the same way as in eq. (37). The only difference is in the value of factors n and C. For conical indenters n=2 and C:

(60)

where β is the angle between the material sample surface and the indenter, γ=π/2 for conical indenters and g=πcot2 β.

Spherical indenters

Sakai (Sakai 2002) derived the solution for nanoindentation of viscoelastic materials using spherical indenters in the same way as for the cylindrical and conical indenters. The load P(t) is defined in the same way as in eq. (37). Again, the only difference is in the value of factors n and C. For spherical indenters n=3/2 and C:

(61)

where γ=2 for spherical indenters

Oyen (Oyen 2005) uses the same three element viscoelastic model as Cheng et al. (Cheng, Xia et al. 1999; Cheng, Xia et al. 2005) and also calls it the standard model although the definition of the standard model is different in other sources (Wikipedia-Viscoelasticity). Although Oyen uses the same viscoelastic model as Cheng et al., the creep function is different in these two papers. Oyen derives only the solution for load control (creep test), since it is more convenient when ramp loading is taken into account. The solution provided by Oyen is valid only for incompressible materials (ν=0.5). Main difference between the work of Oyen (Oyen 2005) with the other two models described here (Sakai 2002; Cheng, Xia et al. 2005) is that Oyen uses ramp loading which is more realistic than step loading, since step loading cannot be achieved in real experiments. In this way the solution is made out of two parts, the ramp part and the hold part:

For 0 ≤ t ≥ tR – Ramp time

(62)

For t ≥ tR – Hold time

(63)

Where k is the slope of the ramp, R is the radius of the spherical indenter and J(t) is the material creep function:

(64)

The ramp correction factor:

(65)

Once the factors C0 and Ci have been obtained, the shear modulus is calculated:

(66)

Cheng et al. (Cheng, Xia et al. 2005) derived solutions for creep and relaxation responses using spherical nanoindeters for step loading and for arbitrary loading history. Only arbitrary loading history will be presented here because of afore mentioned issues concerning step loading. The viscoelastic material is modeled as a three element model (Kelvin-Voight element in series with another spring, Figure 8).

For relaxation test:

(67)

where P(t) is the load, R is the radius of the indenter, h(Ï„) is the displacement (depth of penetration), t is the time, Ï„ is the material time constant and Gr is the relaxation modulus defined as:

(68)

(69)

(70)

(71)

(72)

(73)

(74)

(75)

(76)

(77)

(78)

(79)

(80)

(81)

For creep test:

(82)

where h(t) is displacement, R is the indenter radius, P(Ï„) is the load, t is time, Ï„ is material time constant and Jc is the creep compliance:

(83)

(84)

(85)

(86)

(87)

(88)

Indentation experiments on polystyrene showed that data fit the formulas with less than 2% deviation. Values of viscosity η and elastic modulus E2 obtained using creep and relaxation tests are within 2% difference (Cheng, Xia et al. 2005).

Analysis of poroelastic materials

Poroelastic materials are composed of two phases, a solid phase and a fluid phase. Solid phase has interconnected pores which are filled with fluid. It is essential that the pores inside the solid matrix are connected because this enables the fluid to move through the solid matrix and flow in or out depending on the type of load imposed on the material. This interaction between the fluid and the solid phases gives rise to poroelastic material response which is similar to viscoelastic response since it also exhibits creep and stress relaxation, but is more complicated because it depends on the interaction of two different phases. In most cases the fluid is considered to be incompressible, Poisson ratio ν=0.5, while the solid phase is considered to be linearly elastic and isotropic (Oyen 2011).

There are three major approaches for derivation of the same basic equations of poroelasticity. The main difference between these approaches is the employed averaging process (S.C.Cowin 1999). First approach uses a representative volume element (RVE) for the averaging process. This approach originates in solid mechanics and appears in its early form in the Biot formulation of the effective medium approach (M.A.Biot 1941). The second approach is the mixture theory approach which is based on diffusion models. Main advantage of this approach becomes evident when a number of different fluid types are present in a relative motion (S.C.Cowin 1999). For the mixture theory the averaging is density weighted on the basis of the density of each fluid type in the mixture, while for the RVE it is averaged over a finite volume of the porous medium. For small strains the difference in these two methods is negligible (S.C.Cowin 1999). The third approach was developed by Burridge and Keller (Burridge and Keller 1981) and is called the homogenization derivation approach. Using a two space method of homogenization the dynamic form of the same basic set of equations was derived. In this way macroscopic equations which manage the behavior of the medium on the micro scale are derived.

All three theoretical approaches lead to the same set of equations at the continuum level, although they are focused on different phenomena in the poroelastic material response. The effective medium approach explains the nature of parameters associated with the solid phase. The mixture theory approach is used for averaging over different fluid phases. The homogenization approach explains the dynamical characteristics of the theory (wave propagation) (S.C.Cowin 1999).

There are almost no analytical methods for determining the material propersties of poroelastic materials using nanoindentation. Most researchers focus on numerical methods or FEM analysis (Warner, Taylor et al. 2001; S.Selvadurai 2004; Cao, Youn et al. 2006; Gupta, Lin et al. 2009). One analytical method that will be described here was developed by M. L. Oyen (Oyen 2011). It describes derivation of material properties from creep response of a poroelastic material (hydrated bone) using an impermeable spherical indenter.

The solid part of the poroelastic material is considered to be isotropic and linearly elastic. It is described using two material properties, in this case the shear modulus G and Poisson’s ratio ν. The fluid and fluid solid interaction are described using two additional parameters: the undrained Poisson’s ratio νu and Skempton pore pressure coefficient B. These four parameters completely describe the linear poroelastic solid, but other sets of four parameters could also be used (bulk modulus K instead of shear modulus G). The parameter that describes the flow through the porous elastic solid is called hydraulic permeability (κ) and it is utilized in Darcy’s law:

(89)

Where vi is the fluid velocity and p,i is the pore pressure gradient.

Hydraulic permeability can be calculated using intrinsic permeability (k) which is directly related to the square of a characteristic length scale, often taken as the pore size (Oyen 2011) and fluid viscosity (μ):

(90)

In order to describe the material creep function, viscoelastic hereditary integral is used, developed for viscoelastic material response (eq. 64):

(64)

Using the viscoelastic hereditary integral, the expression for displacement – time (h-t) response following the ramp loading at finite ramp time tR is derived:

(91)

Agbezuge and Deresiewicz (Agbezuge and Deresiewicz 1974) have presented the solution for three sets of different boundary conditions (indenter and half-space surface are permeable, half-space is impermeable, the indenter is impermeable and half-space surface is permeable) using normalized displacement (H*) and normalized time (T*).

(92)

(93)

Where G is the shear modulus, κ is the hydraulic permeability, t is time, a(t) is the contact radius, h(0) is the zero-time depth of penetration and h() is the final displacement following the creep exhaustion. Shear modulus and zero-time depth are known from the zero-time response:

(94)

The contact radius a(t) is equal to:

(95)

Where R is the indenter radius and h(t) is the depth of penetration.

In this case the half-space surface is considered to be permeable and the indenter is considered to be impermeable, so the normalized displacement is approximated by a sigmodial function:

(96)

Where A, P and T0 are empirical fitting parameters A=0.928, P=2.0837 and T0=0.772.

Poisson’s ratio of the solid skeleton is calculated from the final displacement h().

(97)

Linear elasticity

Cylindrical indenter

Conical indenter

Spherical indenter

Force control

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Hertz(1882)

Sneddon(1965)

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Displacement control

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Hertz(1882)

Sneddon(1965)

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Elasto-plasticity

Cylindrical indenter

Conical indenter

Spherical indenter

Force control

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Displacement control

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Pharr, Oliver, and Brotzen (1992)

Hertz(1882)

Sneddon(1965)

Hyperelasticity

Cylindrical indenter

Conical indenter

Spherical indenter

Force control

Choi and Shield (1981)

Giannakopoulos and Panagiotopoulos (2009)

Choi and Shield (1981)

Sabin and Kaloni (1983)

Giannakopoulos and Triantafyllou (2007)

Displacement control

Choi and Shield (1981)

Giannakopoulos and Panagiotopoulos (2009)

Choi and Shield (1981)

Sabin and Kaloni (1983)

Giannakopoulos and Triantafyllou (2007)

Viscoelasticity

Cylindrical indenter

Conical indenter

Spherical indenter

Force control

Cheng, Xia et al. (1999)

N/A

Cheng, Xia et al. (1999)

Oyen (2005)

For 0 ≤ t ≥ tR – Ramp time

For t ≥ tR – Hold time

Displacement

control

Sakai (2002)

Cheng, Xia et al. (1999)

Sakai (2002)

Sakai (2002)

Cheng, Xia et al. (1999)

Poroelasticity

Cylidrical indenter

Conical indenter

Spherical indenter

Force control

N/A

N/A

Oyen (2011)

Displacement

control

N/A

N/A

N/A

Table 1. Analytical solutions for all indenter types and material responses presented in this paper

Agbezuge, L. K. and H. Deresiewicz (1974). "On the indentation of a consolidating half-space." Israel J. of Tchnol. 12.

Burridge, R. and J. B. Keller (1981). "Poroelasticity equations derived from microstructure." J. Acoust. Soc. Am 70: 1140-1146.

Busfield, J. J. C. and A. G. Thomas (1999). "Indentation Tests On Elastomer Blocks." Rubber Chem. Technol. 75: 876-894.

Cao, L., I. Youn, et al. (2006). "Compressive properties of mouse articular cartilage determined in a novel micro-indentation test method and biphasic finite element model." J Biomech Eng 128(5): 766-771.

Chen, Z. and S. Diebels (2012). "Nanoindentation of hyperelastic polymer layers at finite deformation and parameter re-identification." Archive of Applied Mechanics 82(8): 1041-1056.

Cheng, L., X. Xia, et al. (2005). "Spherical-tip indentation of viscoelastic material." Mechanics of Materials 37(1): 213-226.

Cheng, L., X. Xia, et al. (1999). "Flat-Punch Indentation of Viscoelastic Material." Journal of Polymer Science: Part B: Polymer Physics 36.

Cheng, Y. and C. Cheng (1998). "Scaling approach to conical indentation in elastic-plastic solids with work hardening." Journal of Applied Physics 84(3): 1284.

Cheng, Y. and C. Cheng (2004). "Scaling, dimensional analysis, and indentation measurements." Materials Science and Engineering: R: Reports 44(4-5): 91-149.

Choi, I. and R. T. Shield (1981). "Second-order effects in problems for a class of elastic materials." Journal of Applied Mathematics and Physics 32.

D.Tabor (1951). "Hardness of solids." Oxford University Press.

Doerner, M. F. and W. D. Nix (1986). "A method for interpreting the data from depth-sensing indentation instruments." J. Mater. Res 1(4).

E.H.Lee and J.R.M.Radok (1960). "The Contact Problem for Viscoelastic Bodies." J.Appl.Mech. 27.

Ebenstein, D. M. and L. A. Pruitt (2006). "Nanoindentation of biological materials." Nano Today 1(3): 26-33.

Field, J. S. and M. V. Swain (1993). "A simple predictive model for spherical indentation." J. Mater. Res 8(2).

Fischer-Cripps and A. C. (2004). "Nanoindentation."

Gent, A. N. and O. H. Yeoh (2006). "Small Indentations Of Rubber Blocks Effect Of Size And Shape Of Block And Of Lateral Compression." Rubber Chem. Technol. 79: 674-693.

Giannakopoulos, A. E. and D. I. Panagiotopoulos (2009). "Conical indentation of incompressible rubber-like materials." International Journal of Solids and Structures 46(6): 1436-1447.

Giannakopoulos, A. E. and A. Triantafyllou (2007). "Spherical indentation of incompressible rubber-like materials." Journal of the Mechanics and Physics of Solids 55(6): 1196-1211.

Gupta, S., J. Lin, et al. (2009). "A fiber reinforced poroelastic model of nanoindentation of porcine costal cartilage: a combined experimental and finite element approach." J Mech Behav Biomed Mater 2(4): 326-337; discussion 337-328.

Hendrix, B. C. (1995). "The use of shape correction factors for elastic indentation measurements." J.Mater.Res. 10(2).

I.Sneddon (1965). "The Relation Between Load And Penetration In The Axisymmetric Boussinesq Problem For A Punch Of Arbitrary Profile." Int. J. Engng Sci 3: 47-57.

Knoop, F., C. C. Peters, et al. (1939). "A sensitive pyramidal diamond tool for indentation measurments." Journal of Research of the National Bureau of Standards 23.

Kurland, N. E., Z. Drira, et al. (2012). "Measurement of nanomechanical properties of biomolecules using atomic force microscopy." Micron 43(2-3): 116-128.

Larsson, P.-L., A. E. Giannakopoulos, et al. (1996). "Analysis Of Berkovich Indentation." Int. J. Solids Structures 33(2): 221-248.

Lopez, O., K. K. Amrami, et al. (2008). "Characterization of the dynamic shear properties of hyaline cartilage using high-frequency dynamic MR elastography." Magn Reson Med 59(2): 356-364.

This work evaluated the feasibility of dynamic MR Elastography (MRE) to quantify structural changes in bovine hyaline cartilage induced by selective enzymatic degradation. The ability of the technique to quantify the frequency-dependent response of normal cartilage to shear in the kilohertz range was also explored. Bovine cartilage plugs of 8 mm in diameter were used for this study. The shear stiffness (mu(s)) of each cartilage plug was measured before and after 16 hr of enzymatic treatments by dynamic MRE at 5000 Hz of shear excitation. Collagenase and trypsin were used to selectively affect the collagen and proteoglycans contents of the matrix. Additionally, normal cartilage plugs were tested by dynamic MRE at shear-excitations of 3000-7000 Hz. Measured micro(s) of cartilage plugs showed a significant decrease (-37%, P < 0.05) after collagenase treatment and a significant decrease (-28%, P < 0.05) after trypsin treatment. Furthermore, a near-linear increase (R(2) = 0.9141) in the speed of shear wave propagation with shear-excitation frequency was observed in cartilage, indicating that wave speed is dominated by viscoelastic effects. These experiments suggest that dynamic MRE can provide a sensitive quantitative tool to characterize the degradation process and viscoelastic behavior of cartilage.

M.A.Biot (1941). "General Theory Of Three-dimensional Consolidation." Journal of Applied Physics 12(2): 155-164.

Marckmann, G. and E. Verron "Comparison Of Hyperelastic Models For Rubber-like Materials." Rubber Chem. Technol. 79: 835-858.

Mattice, J. M., A. G. Lau, et al. (2006). "Spherical indentation load relaxation of soft biological tissues." Journal of Materials Research 21(08).

Murnaghan, F. D. (1951). "Finite Deformations of an Elastic Solid."

Oliver, W. C. and G. M. Pharr (1992). "An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments." Journal of Materials Research 7(06): 1564-1583.

Oliver, W. C. and G. M. Pharr (2004). "Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology." J. Mater. Res 19(1).

Oyen, M. L. (2005). "Spherical Indentation Creep Following Ramp Loading." Journal of Materials Research 20(08): 2094-2100.

Oyen, M. L. (2006). "Analytical techniques for indentation of viscoelastic materials." Philosophical Magazine 86(33-35): 5625-5641.

Oyen, M. L. (2011). "Poroelastic nanoindentation responses of hydrated bone." Journal of Materials Research 23(05): 1307-1314.

Oyen, M. L. and R. F. Cook (2009). "A practical guide for analysis of nanoindentation data." Journal of the Mechanical Behavior of Biomedical Materials 2(4): 396-407.

Pharr, G. M., W. C. Oliver, et al. (1992). "On the generality of the relationship among contact stiffness,contact area, and elastic modulus during indentation." J. Mater. Res 7(3).

R.B.King (1987). "Elastic Analysis Of Some Punch Problems For A Layered Medium." Int.J.Solid.Struct. 23(12): 1657-1664.

S.C.Cowin (1999). "Bone poroelasticity." J.Biomech. 32: 217-238.

S.Selvadurai, A. P. (2004). "Stationary damage modelling of poroelastic contact." International Journal of Solids and Structures 41(8): 2043-2064.

Sabin, G. C. W. and P. N. Kaloni (1983). "Contact Problem Of A Rigid Indentor In Second Order Elasticity Theory." Journal of Applied Mathematics and Physics 34.

Sakai, M. (2002). "Time-dependent viscoelastic relation between load and penetration for an axisymmetric indenter." Philosophical Magazine A 82(10): 1841-1849.

Vlassak, J. J. and W. D. Nix (1994). "Measuring The Elastic Properties Of Anisotropic Materials By Means Of Indentation Experiments." J. Mech. Phys. Solids, 42(8): 1223-1245.

Wagner, M. H. (1994). "The origin of the C2 term in rubber elasticitya)." Journal of Rheology 38(3): 655.

Warner, M. D., W. R. Taylor, et al. (2001). "Finite element biphasic indentation of cartilage: A comparison of experimental indenter and physiological contact geometries." Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 215(5): 487-496.