.v‘a.

_ . 2:5}- ‘

‘ ,H‘.‘ 2“
A3429 ' '.
vhmfik .35

r... s. .
. .... HIV: .mdir.
. ..- 2.. . .
. .. n. : I
. , mwrwf.
unmufluwwuhmw ..
3'2:

 

z
34.: »

.q.

'«za
, ‘u
:1‘ 5-? .
..
.J‘ .3
.
l
A“,
1‘”
.4

h
1}“
Infi- _ _

Jr‘
1
4 in
(

2.
up;
of?
...
m‘
I": : u
i.

g! . 7
3W
#2::
M :1
*3

Z .

i
451

. Wmvu4nir mauve o ! ...
«man... ....Ahzwum. d fimwhwrw. ..
yum... ...:m. . “bay...

5

.Vxnffllu
i .

£5“?

4)

In. 4,.»
I ...
2....Irvotkubv.
Va: 010 .l.
a I :0
4am...§i.. ..
$-u9fl4dud ......
(lupinfii
«Ii

.....l O
ll-Olf'lli
Q l

....
...”..me

00 l D..-‘v.l.|bo.vl.
. .11!‘
E .50.... .7kqurl.
11.0%) v1.1.0!
II). Eli)!“
.1... ..D gig.
.Illl); tie)...“
is I;
1.81.. (.11.:111
. (I

{it

...”.

 

 

 

 

4 .-...A. || 1.
§ .I.10L...l ..-...
,t‘vl *vttvllutltrw..l

— u
3.11.5..Arohn l. ..

ID I «LA-u.

 

 

  

....,§&..W$fi$ .mefiwmwfifi

A.
”. .-- ...... .1}... . . .
U .. ...... . .1.... 3.. .--: 3 r . .
M v . : .-. .. E . ..I. ... 3-..... . ..u..u..m..MA¢.Lar&MWwWF....AHv
'21:... -2: 2..
I {y
L r

 

1. hut-17$

Illllllllllllllllllllllllllllllllll

3 1293 01688 0621

 

 

 

 

 

 

 

 

This is to certify that the

dissertation entitled

A TRANSVERSELY ISOTROPIC HYPO-ELASTIC BIPHASIC

MODEL OF ARTICULAR CARTILAGE UNDER IMPACT LOADING

presented by

Jose Jaime Garcia

has been accepted towards fulfillment
of the requirements for

Ph . D . degree in Mechanic 3

 

 

{uojlmlo/qW flfifiw

Maj profersor

Date 6%4/08 .J
/ I

MS U is an Affirmative Action/Equal Opportunity Institution 0- 12771

 

 

 

 

 

LIBRARY

Michigan State
University

 

 

 

PLACE IN RETURN BOX
to remove this checkout from your record.
TO AVOID FINE-S return on or before date due.

 

I MTE DUE DATE DUE DATE DUE

«$231M: :'

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1’98 WM“

A TRANSVERSELY ISOTROPIC HYPO-ELASTIC BIPHASIC

MODEL OF ARTICULAR CARTILAGE UNDER IMPACT LOADING

By

Jose Jaime Garcia

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements for the degree of
DOCTOR OF PHILOSOPHY
Department of Materials Science and Mechanics

1998

ABSTRACT

A transversely isotropic hypo-elastic biphasic model is presented to describe the
mechanical behavior of articular cartilage under impact loading. This model allowed the
consideration that the Young’s modulus in the plane of the surface (plane of isotropy) is
higher than that in the direction of loading, which is consistent with the microstructural
organization of the tissue. Unlike the isotropic model, the transversely isotropic model
was capable of fitting high peak loads observed during the initial stage of rapid indentation
on the tissue. This made it possible to devise a method of estimating the infinitesimal
material properties from small strain indentation relaxation experiments. Infinitesimal
analyses using finite element models showed significant differences in the stress distribu-
tions between transversely isotropic and isotropic models for idealized contact problems.
Stress distributions in transversely isotropic models were more consistent with surface
injuries observed during impact experiments. A hypo-elastic model for the solid skeleton
was used to model the nonlinear behavior of cartilage under large deformations. Using this
model, the nonlinear elastic properties of rabbit retropatellar cartilage were determined
from indentation, relaxation experiments. Experimental and theoretical results revealed
the important role of interstitial fluid in supporting compressive stress during rapid impact
loading with large deformations of retropatellar cartilage. A significant strain-stiffening
effect was observed for rapid loading, limiting tissue strains during blunt insult to the
joint. The hypo-elastic properties were used to obtain shear stresses in a plane strain model
of the rabbit patella under impact loading. Peak shear stresses on the surface were consis-

tent with surface injuries observed in our rabbit model of post-traumatic osteoarthrosis.

DEDICATION

To my wife Olga Marina and my children Lina Maria and Jose Alejandro, whose
affection has been a continuous source of energy during the course of these studies. To my
mother Matilde and the memory of my father Camilo who taught me to be responsible in

all activities of my life.

iii

ACKNOWLEGMENTS

My profound appreciation to my advisors, Professors Nicholas J. Altiero and
Roger C. Haut, for their guidance and advice during numerous hours we spent together
discussing different aspects of this dissertation. Their assistance has been very valuable in
introducing me in the discipline of Biomechanics.

My gratitude to Professors D. Liu and C. DeCamp for their advice, comments and
support during my studies.

This study was supported by a grant from the Centers for Disease Control and Pre-
vention (R49/CC503607). Its contents are solely the responsibility of the author and do
not necessarily represent the official views of the Centers for Disease Control and Preven-
tion. My studies were also made possible by the support of the Universidad de Valle, Cali,

Colombia.

iv

TABLE OF CONTENTS

LIST OF TABLES .................................................... ix

LIST OF FIGURES ................................................... x

CHAPTER 1 INTRODUCTION ......................................... 1
1.1 Osteoarthrosis ............................................... l
1.2 Isotropic biphasic model for articular cartilage ..................... 3
1.3 Inconsistencies of the isotropic model ............................ 5
1.4 Large deformation under impact events ........................... 5
1.5 Proposed transversely isotropic model ............................ 6

CHAPTER 2 AN APPROACH FOR THE STRESS ANALYSIS OF

TRANSVERSELY ISOTROPIC BIPHASIC CARTILAGE

UNDER IMPACT LOAD ......................................... 9
2.1 Derivation of Effective Elastic Material Properties .................. 9
2.2 Numerical Methods ........................................... 16
2.2.1 Verification ............................................. 16
2.2.2 Stress Analysis .......................................... 17
2.3 Results ..................................................... 18
2.3.1 Verification ............................................. 18
2.3.2 Stress Analysis .......................................... 18
2.4 Discussion .................................................. 19

CHAPTER 3 A METHOD TO DETERMINE IN SITU MATERIAL

PROPERTIES OF TRANSVERSELY ISOTROPIC BIPHASIC

CARTILAGE .................................................. 28
3.1 Numerical Method ........................................... 28
3.1.1 Determination of candidate elastic constants ................... 28
3.1.2 Determination of candidate isotropic permeabilities ............. 32
3.1.3 Determination of actual elastic constants ...................... 33
3.1.4 Estimation of anisotropic permeabilities ....................... 34
3.1.5 Automation of the algorithm ................................ 35
3.2 Experimental Method ......................................... 37
3.3 Results ..................................................... 38
3.4 Discussion .................................................. 38

CHAPTER 4 ESTIMATION OF IN SITU ELASTIC PROPERTIES OF

BIPHASIC CARTILAGE BASED ON A TRANSVERSELY

ISOTROPIC HY PO-ELASTIC MODEL ............................ 50
4.1 Theoretical Analysis .......................................... 50
4.2 Methods .................................................... 5 7

4.2.1 Experimental ............................................ 57
4.2.2 Numerical .............................................. 58
4.3 Results ..................................................... 60
4.4 Discussion .................................................. 61

vi

CHAPTER 5 PLAN STRAIN ANALYSIS OF THE RABBIT

RETROPATELLAR CARTILAGE UNDER IMPACT LOADING ......... 70
5.1 Numerical method ............................................ 70
5.2 Results .................................................... 71
5.3 Discussion ................................................. 72
CHAPTER 6 CONCLUSIONS .......................................... 76

APPENDIX A DESCRIPTION OF THE PROGRAM T [PRORM USED TO

DETERMINE INFINITESIMAL TRANSVERSELY ISOTROPIC

PROPERTIES OF CART ILAGE ................................... 81
Al Instructions for the user of the program tiprop.m ................... 81
A. 1.1 Data provided by the user .................................. 81
A.l.2 Programs and data necessary to run the program ................ 81
A. l .3 Steps to run the program .................................. 82
A2 Description and list of the program .............................. 83
A3 Names and organization of resident data files ...................... 100
A4 Procedure used to calibrate the master curve ....................... 101

APPENDIX B SOME EQUATIONS FOR THE HYPO-ELASTIC AND
HYPERELASTIC MODELS ...................................... 104

3.1 Integration of the grade zero hypo-elastic equation for the simple

extension problem ............................................... 104

B2 Relations between the engineering constants and the coefficients

of the strain energy function ....................................... 105

vii

B.3 Analysis of uniaxial impact loading according to the
hyperelastic formulation .......................................... 106

B4 Description and list of notebooks used in the nonlinear analysis ....... 107

3.5 Subroutine hypela.f used to define the components of the
hypo-elastic tensor in the program MARC ............................ 124

viii

LIST OF TABLES

Table 2.1 Skeleton material properties from Cohen et al. (1993) ................ 23
Table 2.2 Maximum pressures and aspect ratios .............................. 24

Table 2.3 Maximum shear and normal hoop stresses in the solid skeleton

(solid volume fraction = 0.25) at the locations A, B and C (Figure 2.4) ..... 24

Table 3.1 Transversely isotropic properties of rabbit retropatellar cartilage ........ 43

LIST OF FIGURES

Figure 2.1 Stress decomposition in a biphasic element at time zero .............. 25

Figure 2.2 Maximum shear stress contours for (a) incompressible elastic

and (b) biphasic at time zero ....................................... 26

Figure 2.3 Contours of radial effective stresses for (a) incompressible

elastic and (b) biphasic at time zero ................................. 26

Figure 2.4 Contours of maximum shear stress for the 200 mm indenter
and (a) I, (b) TI and (c) I-TI models ................................. 27

Figure 3.1 Typical force-time relaxation curve from an indentation experiment ..... 43

Figure 3.2 Points in a typical fpe, feq plane for an isotropic biphasic layer of
0.5 mm. The arrows represent a typical experimental point for which

the Poisson’s ratio was equal to -0.1 ................................. 44

Figure 3.3 Points in a typical fpe, feq plane where 1311’: 5.0 MPa, .35.): 0.10,
and thickness = 0.5 mm. The indicated linearly interpolated point

corresponded to E253) = 1.4 MPa and 613 = 0.65 MPa ................... 45

Figure 3.4 Force-time relaxation theoretical curves obtained for two sets of
material properties for the (a) 1 mm diameter indenter and

(b) 1.5 mm diameter indenter ...................................... 46
Figure 3.5 Simple model to estimate anisotropic permeabilities ................. 47

Figure 3.6 Comparison of the master curves with the force-time curves

from Abaqus, 3:4 and or = 1 ....................................... 48

Figure 3.7 Comparison of the experimental curves with the curves from
Abaqus run with the constants and anisotropic permeabilities predicted
with the program tiprop.m ........................................ 49

Figure 4.1 Stress-Stretch curves from the uniaxial solution with a constant

coefficient a1=0.8 and different values for the coefficient

11 (0.020, 0.067 and 0.113) ........................................ 65

Figure 4.2 Stress-Stretch curves from the uniaxial solution with a constant
coefficient n = 0.1 13 and different values for the coefficient
a] (0.8, 2.0, 3.0) ................................................ 66

Figure 4.3 Experimental force-displacement curves obtained with the spherical
indenter of 2 mm diameter showing the typical degree of fit obtained
with the theoretical model ......................................... 67

Figure 4.4 Comparison of stress and lateral displacement between
the hypo-elastic and hyperelastic solutions for uniaxial loading ........... 68

Figure 4.5 Uniaxial stress-stretch curves obtained using the average properties
from this study for rabbit retro-patellar cartilage. Observe the substantial

contribution of water pressure to total stress in the rapid loading response. . .69

Figure 5.1 Plane strain model for the effective analysis of the
retropatellar cartilage ............................................ 74

Figure 5.2 Original and deformed shapes of the cartilage and contours of shear
stress for the (a) low severity and (b) high severity cases ................. 75

Figure A] Functions of the program tiprop.m ............................... 83

CHAPTER 1

INTRODUCTION

1.1 Osteoarthrosis

One of the main functions of a diarthrodial joint, such as the knee or shoulder, is to
facilitate body movement and locomotion (Mow et al., 1984). Articular cartilage is the
white material covering the ends of bones in diarthrodial joints. It is a highly efficient
bearing system which provides insignificant friction and presents negligible wear for the
entire life span of an individual. Articular cartilage is composed of a solid skeleton swol-
len by water. The solid skeleton, which accounts for approximately 25% of volume, is
formed by collagen fibers, proteoglycans and chondrocytes (Mow et al., 1991). The fluid
plays an important role in the capacity of cartilage to transmit loads, provide lubrication

and n'ansport nutrients and waste products (McCutchen, 1978; Maroudas, 1975).

Osteoarthrosis is a disease which affects cartilage and subchondral bone of diar-
throdial joints. This disease may account for a societal cost of $54 billion annually in both
treatment costs and lost work days (Atkinson et al., 1997). Osteoarthrosis is diagnosed
radiographically by a reduction of joint space, as articular cartilage is lost and bone spurs
appear (Hamerman, 1989). Mechanical insult has long been associated with the initiation
of this disease, but clinical correlations are extremely difficult (Kellgren and Lawrence.

1958; Davis et al., 1989). Factors leading to the progression of osteoarthrosis are not fully

understood. Brittberg et al. (1994) suspect that minor lesions disrupting either the cartilage
or the subchondral bone lead to osteoarthrosis. It has been hypothesized that softening of
the cartilage after blunt impact would increase stresses in the subchondral bone producing
microcracks and sclerosis that would lead to degeneration of the overlying cartilage
(Donohue et al., 1983; Newberry et al. 1997). Prevention of acute tissue damage from
blunt impact loading, or early diagnosis and intervention, may mitigate the disease process
(Newberry and Haut, 1996). However, a better understanding of the disease etiology is

necessary to solve this problem.

Animal studies have been conducted to better understand the relationship between
a single blunt impact and the pathogenesis of osteoarthrosis. Thompson et al. (1991) docu-
ment occult microcracks at the cartilage-bone interface and surface injuries following
blunt insult to the canine patellofemoral joint. Armstrong et al. (1985) reported separation
of the articular cartilage from the subchondral bone and reduction in the mechanical prop-
erties of the cartilage surface layer in blunt experiments on the patellofemoral joint of pig
specimens. In experimental studies from our laboratory on the rabbit patellofemoral joint
(Haut et al., 1995; Newberry et al., 1997, Li et al., 1995), surface injuries have been
reported on retropatellar cartilage without apparent injury in the subchondral bone or at
the bone-cartilage interface. In these experiments, impact rise times and total contact
times were similar to those reported by Armstrong et al. (1985) and Thompson et al.
(1991). Newberry et al. (1998) documented early stages of osteoarthrosis, characterized
by degenerative changes in cartilage and in subchondral bone of the patella, after a single

blunt impact to the patellofemoral rabbit joint.

1.2 Isotropic biphasic model for articular cartilage

A biphasic model has been developed to represent the mechanical behavior of
articular cartilage (Mow et al., 1980). In this model, cartilage is represented by two intrin-
sically incompressible phases: solid and fluid. Stresses in the tissue can be decomposed
into the pressure in the fluid and the stresses in the solid skeleton. Given the low penne-
ability of articular cartilage, the amount of water that can leave the tissue during an impact
event is negligible and, therefore, the external mechanical response of the tissue can be
represented as that of an incompressible elastic material (Eberhardt et al. 1990). The effec-
tive elastic properties of the isotropic, incompressible material were presented by Ann-

strong et al. (1984).

Numerous theoretical models for isotropic cartilage under short time loading indi-
cate that the highest shear stresses occur at the interface with bone and that compressive
normal stresses occur near the surface in the zone of contact (Armstrong, 1986; Eberhardt
et al., 1990; Eberhardt et al., 1991; Ateshian et al., 1994; Kelly and O’Connor, 1996; Wu
et al., 1996; Askew and Mow, 1978). It is only for small aspect ratios that shear stresses
are maximum on the contact surface (Atkinson et al., 1997) and tensile stresses develop on
the surface of the cartilage near the edge of the contact area (Eberhardt et al., 1991). Plane-
strain models of the rabbit patellofemoral joint under impact loading also show the highest
shear stresses at the bone-cartilage interface (Li et al., 1995; Garcia et al., 1997). However,

the experiments only describe surface injuries on retropatellar cartilage.

Both shear and tensile stresses have been associated with injuries in articular carti-

lage (Kelly and O’Connor, 1996; Silyn-Roberts and Broom, 1990). Recent impact experi-

ments undertaken by Atkinson et a1. (1997) with small radius spherical impactors, have
shown a strong correlation of surface injuries with high levels of shear stress. The radially-
oriented nature of these fissures also suggests that hoop tensile stresses might be responsi-
ble for these surface lesions. However, maximum principal stresses in the solid skeleton of
the cartilage, predicted with a computational isotropic model, were not well correlated
with the occurrence of surface fissures (Atkinson et al., 1997). We hypothesize that sur-
face lesions during blunt experiments are due to high stresses in the solid skeleton of the
cartilage. However, based on analyses of isotropic models, surface injuries in blunt exper-
iments on the patellofemoral joint of the rabbit are inconsistent with a stress-based failure

criterion, if one assumes the same strength for all depths of cartilage.

Several methods have been established to determine the mechanical properties of
articular cartilage. In the initial isotropic model presented by Mow et al. (1980), creep and
relaxation tests on confined specimens were used to determine the elastic constants of the
solid skeleton and the permeability of the tissue. Armstrong et al. (1984) later conducted
relaxation tests on unconfined specimens to determine the isotropic properties of articular
cartilage. Those tests yielded ratios of initial to equilibrium load higher than a maximum
value of 1.5 predicted by the isotropic biphasic model. The authors attributed the problem
of model fit to the fact that frictional effects which were thought to develop on the contact-
ing surface during these experiments were not included. Subsequent analyses by Spilker et
al. (1990), however, showed that high ratios of peak to equilibrium forces in unconfined
relaxation experiments can only be partially explained by friction on the load interface.
Using creep indentation experiments Mow et al. (1989) determined the in situ isotropic

properties of articular cartilage. The analysis, however, had difficulty fitting the early time

response of the cartilage. In creep indentation experiments on canine femurs Jurvelin et a1.
(1988) found that the initial and equilibrium loads could be fit only with negative values of
Poisson’s ratio using the isotropic, biphasic model. Indentation relaxation tests undertaken
in our laboratory on rabbit patellar articular cartilage show typically high ratios of peak to
equilibrium forces (Newberry et al., 1997), which also could only be explained with nega-

tive Poisson’s ratios using an isotropic theory.

1.3 Inconsistencies of the isotropic model

In summary, two inconsistencies arise from the use of an isotropic model for artic-
ular cartilage. First, stress analyses in patellofemoral contact problems using isotropic
models reveal that shear stresses are maximum at the interface with the bone, and that
principal stresses are compressive at the surface, in the area of contact. This distribution is
not consistent with experimental studies on the rabbit patellofemoral joint (Haut et a1,
1995; Newberry et al., 1997) which show surface fissures on the retropatellar cartilage in
the zone of contact without apparent injury at the bone-cartilage interface. And secondly,
the isotropic biphasic model is unable to fit the high ratios of peak to equilibrium load
observed in stress relaxation experiments on the cartilage tissue, unless a negative Pos-

son’s ratio is assumed which is physically unlikely.

1.4 Large deformations under impact events

Impact experiments conducted on the rabbit and human patellofemoral joint show
maximum contact pressures of 40-50 MPa (Haut et al., 1995; Newberry et al., 1997;

Atkinson et al., 1995). Based on these pressures, excessively large deformations are pro-

duced in finite element analyses of rabbit and human patellofemoral joints under impact
loading when one uses the isotropic elastic constants derived from small strain experi-
ments reported in the literature, e.g. Young’s modulus from 0.4 to 1.0 MPa (Athanasiou et

al., 1992) and shear modulus values of approximately 0.97 MPa (Newberry et al., 1997).

Articular cartilage is known to behave nonlinearly for large deformations (Mow et
al., 1984). Hyperelastic functions have previously been proposed to describe the nonlinear
behavior of the solid phase of isotropic (Holmes 1986; Kwan et al., 1990; Holmes and
Mow, 1990) and transversely isotropic (Almeida et al., 1995) articular cartilage. In addi-
tion, Ateshian et a1. (1997) employed the model proposed by Holmes and Mow (1990) to
determine isotropic nonlinear material properties of articular cartilage by fitting stress
relaxation curves from confined compression experiments. Almeida and Spilker (1997)
also recently presented a finite element code to model the nonlinear behavior of biphasic
soft tissues under more general loading conditions. To date, however, there has been no
attempt to study the nonlinear response of cartilage in the early stages of rapid indentation.
This information is crucial if one is to understand the state of tissue stress under impact
loading to joints, and the role of interstitial fluid pressurization and stresses in the solid

skeleton on injury mechanisms.

1.5 Proposed Transversely Isotropic Model

The microstructure of the cartilage. with a superficial zone formed by sheets of
tightly woven collagen fibrils (Clark, 1990) suggests a continuum model with a Young’s
modulus in the plane of the cartilage different than that in the direction perpendicular to

the surface. A transversely isotropic model captures this characteristic.

Cohen et al. (1993) have shown that a transversely isotropic model is able to pro-
vide excellent fits for the early response period of force-time curves obtained in creep
indentation experiments on articular cartilage. In addition, Donzelli and Spilker (1995)
recently reported that a transversely isotropic model shifts maximum stresses from the
interface with the subchondral bone to the surface of the cartilage where we see impact
induced surface lesions in our animal model. Based on these results, a transversely isotro-
pic biphasic model was adopted in this dissertation. Wave effects were not considered in
this study, based on the analyses by Anderson et a1. (1991) and Li (1995) who showed that
the amplitude of the waves developed in articular cartilage during impact events is not

large enough to appreciably alter the prevailing quasi-static stress distribution.

In Chapter 2, the stress distribution in the solid skeleton of transversely isotropic
articular cartilage under impact load is obtained by the analysis of an incompressible elas-
tic material equivalent to the biphasic cartilage at time zero. The elastic properties of the
transversely isotropic equivalent material are obtained using stress-strain relations and the
incompressibility condition. A subsequent stress analysis with a commercial finite element
code revealed peaks of maximum shear stress and tensile stress near the surface of carti-
lage in the contact zone. The results appear consistent with observed surface lesions in our
impact experiments using the rabbit, and help support the notion of a stress-based failure

criterion for injury to articular cartilage under blunt insult.

An efficient algorithm is presented in Chapter 3 to obtain transversely isotropic
biphasic properties of articular cartilage, based on small strain relaxation tests. These tests
are easier to undertake than low load creep tests with the majority of existing test equip-

ment, such as a servohydraulic testing machine. The estimated properties provided a good

fit for the entire range of the relaxation response curve, including the early time response
of the tissue. Transversely isotropic properties are reported for rabbit retropatellar carti-

lage.

In Chapter 4 a method is presented to determine the non-linear elastic properties of
biphasic cartilage based on a transversely isotropic hypo-elastic model. The elastic proper-
ties were estimated by fitting two force-displacement curves (in rapid loading and equilib-
rium) obtained from large deformation, indentation relaxation tests on cartilage using a
non-porous spherical indenter. The model was able to fit the material response for rapid
loading, and equilibrium indentation test data to approximately 50% strain. This material
model suggested even higher percentages of stress supported by the fluid phase of carti-
lage during impact loading versus that given earlier by small deformation theories for a

biphasic material.

In Chapter 5 a plane strain analysis of the rabbit patella under impact loading is
presented using the hypo-elastic properties reported in Chapter 4. The shear stress distri-
butions from this model were consistent with surface injuries observed during high energy

impact experiments on the patellofemoral joint.

CHAPTER 2

AN APPROACH FOR THE STRESS ANALYSIS OF TRANSVERSELY
ISOTROPIC BIPHASIC CARTILAGE UNDER IMPACT LOAD

The stress distribution in the solid skeleton of transversely isotropic contact mod-
els of articular cartilage was obtained by the analysis of an incompressible elastic material
equivalent to the biphasic cartilage at time zero. The elastic properties of the transversely
isotropic equivalent material were obtained using stress-strain relations and the incom-
pressibility condition. A subsequent stress analysis with a commercial finite element code
revealed peaks of maximum shear stress and tensile stress near the surface in the contact
zone for transversely isotropic models of articular cartilage. The results appear consistent
with observed surface lesions in our rabbit experiments, and help support the notion of a

stress-based failure criterion for lesions to articular cartilage under blunt insult.

2.1 Derivation of Effective Elastic Material Properties

From the biphasic theory (Mow etal., 1980), stresses in a biphasic cartilage can be

decomposed as follows

(8) (f)
Gij = Oij +Gi' (2.1)

. . f . .
where 6:13) are the stresses 1n the solrd skeleton and oil.) are the stresses 1n the flurd.

The stresses in the solid skeleton can be found in terms of the pressure in the fluid,

p, and the so-called effective stresses, ope (Terzaghi, 1943)

lj)’

om: —¢"’p8ij + 6(8) (2.2)

in which 9(5) is the volume fraction of the solid phase and Sij is the Kronecker delta.

Stress in the fluid is equal to the pressure averaged over the whole volume, i. e.

o‘j” .. —p¢“’ 5.1.9.3)

in which om is the volume fraction of the fluid.

If equations (2.2) and (2.3) are substituted into equation (2.1), the total stresses in

the tissue can be written as

— pfiij+ (SF-e) (2.4)
Since both the fluid and the material forming the solid skeleton are assumed to be

incompressible, the strains in the tissue are due only to the effective stresses.

For biphasic materials composed of particles, effective stresses are proportional to
the contact forces among particles (Scheidegger, 1972), and this is the reason why effec-
tive stresses are widely used in the field of soil mechanics to predict failure. For a continu-
ous solid skeleton, as in articular cartilage, stresses in the solid skeleton are proportional to
the forces transmitted by the solid skeleton averaged over the total volume (Scheidegger,
1972). Thus, skeleton stresses, rather than effective stresses, would appear to be a more

appropriate predictor of change in articular cartilage.

For an isotropic material the incompressibility condition is attained if Poisson’s

10

ratio is 0.5. Similar restrictions on the Poisson’s ratios can be derived for a transversely
isotropic material. If 1-2 is the plane of isotropy, the compliance matrix relates the stress

and strain tensors in their vectorial description as follows (Christensen, 1979):

an 1/511 —v12/E“ 4231/1333 0 0 0 (on

£22 —v12/Ell l/Ell —v31/E33 0 0 0 622

533 =H‘V31/E33 "V3r/E33 1/E33 0 0 0 033 (2.5)
2223 o 0 o 1/Gl3 0 0 023

2213 o o 0 0 1/(313 0 013
gem _ 0 0 o 0 o 2(1+v12)/E1U_012J

 

 

 

 

 

 

where cij and eij are the components of the stress and strain tensors, respectively. The five
constants are the Young’s moduli, E11, E22, the Poisson’s ratios v12, v31, and the shear

modulus 613.
For the compliance matrix to be positive definite it is required that
E11 >0, E33>0, G13>0 (2.6)

Viz 5 0 (2.7)

and

 

(2.8)

For infinitesimal strains, the mathematical condition for incompressibility is a null
trace of the strain tensor. By considering the special case of a hydrostatic stress, it can be
shown that the incompressibility condition is reached if

2(1—v12)+(1—4v31) _

0 (2.9
E 11 E33 )

 

ll

23

ll;

which reduces to v = 0.5 for an isotropic material.

Eliminating v12 between equations (2.8) and ( 2.9) yield

(v31— 0.5)2 s 0 (2.10)

from which it is concluded that the incompressibility condition can be attained only if v31

is 0.5. Substituting this value into equation (2.9), we obtain

so that, if the material is isotropic, it follows that v12 is also 0.5.

It has been shown that the instantaneous response of biphasic cartilage is equiva-
lent to that of an incompressible material (Armstrong et al., 1984; Ateshian et al., 1994).
In addition, Eberhardt et al. (1990) calculated for typical properties of cartilage, that the
total stress in the unconfined test at 200 ms is only 0.14 % different than that at time zero.
Therefore, considering that typical durations of blunt impact experiments are lower than

50 ms, cartilage was treated as an incompressible material.

We hypothesized that it is possible to find a transversely isotropic incompressible
elastic material equivalent to a transversely isotropic biphasic material at time equal zero
for a general stress field. In Figure 2.1, the general stress field, characterized by the normal

stresses 0'11, 0’22 and 033, is decomposed into the pressure in the fluid and the effective

stresses in the solid skeleton, following equation (2.4). The pressure in the fluid at time

zero is proportional to the trace of the stress tensor, okk,

where for an isotropic material the parameter or is equal to 1/3, and the fluid pressure is

equal to the mean stress. The decomposition applies only to normal stresses, since it is

assumed that the fluid sustains only pressure.

At time zero, the deformation of the incompressible tissue under the stress field (I)
of Figure 2.1 must be equal to the deformation of the solid skeleton under the effective
stress field (111) of Figure 2.1, if the fluid is assumed to be incompressible. Therefore, if
the elastic constants of the equivalent incompressible elastic tissue at time zero are

denoted by

and the elastic constants of the solid skeleton (which are equal to the elastic constants of

articular cartilage in the equilibrium state) are denoted by
( ) (S) (S) (8)
E181, E33, V319 V129 G13
then equating the strains in (I) and (III) obtained using equation (2.5) yields

(8) (S)
011 v12622 V31633 _ (Grr‘ackk)_vrz(°22—°‘Gkk) V31(633—a6kk)

 

 

 

_ __ _ __ _ (2.14)
E11 E11 E33 13(181) E281) Egg)

(8) (S)
6—22_v12011_v31033 = (Ozz'aokk)-vlz(oll‘a0kk) v31 (O33—u‘ckk) (215)
IE 13 I3 (8) (S) (s) O

11 11 33 E11 E11 E33

(3) (S)

533 _ V31611_ V31022 = (0'33 - (101(k) _ V31(611_a6kk) V31(G33_°‘Okk) (216)

E— E E (s) (s) (s)
33 33 33 E33 E33 E33
Note that the shear modulus, 613, is the same for the equivalent incompressible elastic

material and the solid skeleton.

13

This system must be satisfied for any combination of 011, 622 and 0'33 . Substitut-

ing the Poisson’s ratios of equation (2.13) and subtracting equation (2.15) from equation

(2.14), one obtains

Boll—[3022 = 0 (2.17)

where B is a parameter in terms of the elastic constants and or. For a general stress field,
equation (2.17) can be satisfied only if [3 is equal to zero. in which case the following rela-

tion is obtained

 

(S)
E E
E33= (1)- ” ” <1 (2.18)
After eliminating 633 between equations (2.15) and (2.16) one obtains
"011‘“:022 = 0 (2.19)

where n and § are parameters in terms of the elastic constants and 01. Under the condition
of equation (2.18), n is equal to i. Thus, for a general stress field, equation (2.19) can be

satisfied only if § is equal to zero.

Equating E, to zero under the condition of equation (2.18), the elastic modulus of
the equivalent incompressible material is found in terms of the elastic constants of the

skeleton as follows

E‘lj’[133—4v‘3" + 2(1—v‘s’)(E:,§’/E§:)]

 

Eu= (2.20)

1 + (1 — v12
If equation (2.20) is substituted into equation (2.18), the other equivalent modulus

is found to be

14

Egg)“ 4v§j’+2(1—v‘1;’)(E‘3§’/E‘lj’)]
E33 = v(S) (S) (s) (8)2 ' (2'21)
2(112)(E33/E11) 4"31

 

If equations (2.20) and (2.21) are substituted into equation ( 2.16), or can be found

in terms of the elastic constants of the skeleton and the stresses as

=(011+622)[1 v9- v§j’(E§j’/:§§’)]+033(112v‘;’)(E‘,i’/E§§’)

s s s s 331535 (2°22)

 

In the case of isotropy, equations (2.20) and (2.21) give an equivalent incompress-
ible Young’s modulus equal to three times the shear modulus of the skeleton. This same
result is obtained for the unconfined compression of cartilage by Armstrong et a1. (1984)
by solving the differential equations of the biphasic mixture. It is also noted that 01 reduces

to the well known value of 1/3 for the isotropic case.

It can be shown algebraically that the in-plane shear modulus 612 calculated using

the equivalent incompressible elastic properties defined by equations (2.13), (2.20) and

(2.21) is equal to that defined in terms of the properties of the skeleton, i.e.

(8)
E11 _ E11

G : —— _ .—
12 2(1+v12) 2(1+v(182))

(2.23)

In summary, it has been shown that the solution of any biphasic transversely iso-
tropic problem at time zero can be made equivalent to that of an incompressible trans-
versely isotropic elastic problem, where the material constants given by equations (2.13),
(2.20) and (2.21) are used along with the stress decomposition defined by the parameter or

of equation (2.22).

15

2.2 Numerical Methods

2.2.1 Verification

The equivalence between the incompressible transversely isotropic elastic material
and the biphasic transversely isotropic material at time zero was verified using the
poroelastic capabilities of the program Abaqus 5.5 (Hibbit, Karlsson & Sorensen, Inc.,
Pawtucket, RI, USA). Mow and Lai (1980) showed that the biphasic theory reduces to the
same poroelastic governing equations developed by Biot (1941) under the current assump-

tions.

The geometry of the problem was taken from that used by Ateshian et al. (1994).
This geometry is relevant because it can approximately reproduce the aspect ratio docu-
mented in contact analysis of diarthrodial joints. It consisted of a large, axisymmetric
cylindrical disk of cartilage 1 mm thick, indented by a frictionless sphere of 400 mm
radius. A 100 N load was applied to the spherical indenter and the bottom nodes of the
disk were fixed. The mesh for both analyses consisted of 640 eight-node quadrilateral ele-

ments, CAX8P for the poroelastic analysis and CAX8 for the elastic analysis.

The poroelastic analysis used the soil consolidation capabilities of Abaqus 5.5 with
the material properties of the solid skeleton as reported by Cohen et al. (1993), see Table
2.1. Material properties for the elastic analysis were then obtained from equations (2.13),

(2.20) and (2.21).

16

2.2.2 Stress analysis

A comparative stress analysis was performed using three axisymmetric models of
the articular cartilage. The geometry was the same as that used for the verification analy-

sis. The models considered were:
a. isotropic cartilage (I)
b. transversely isotropic cartilage (TI)

c. isotropic cartilage with a transversely isotropic superficial layer of thickness
0.125 mm (LTD. The thickness of the superficial layer in the I-TI model was chosen
because it is representative of the superficial tangential zone which is between 10% and

20% of the total thickness of the cartilage (Mow et al., 1991).

The equivalent elastic incompressible properties were calculated from the proper-
ties of the solid skeleton reported by Cohen et al. (1993) for both, isotropic and anisotropic
models. The mesh consisted of 640 eight—node quadrilateral elements. The spherical
indenters had radii of 20, 50, 200, and 400 mm. For each spherical indenter, the magnitude
of the applied load was chosen so that the maximum vertical displacement at the cartilage
surface was 0.1 mm for the (I-TI) model. This maximum displacement ensured that strain
levels remained small as was assumed in this analysis. Stresses in the solid skeleton were
calculated using equation (2.2). Both, shear and normal hoop stresses were calculated as

they may be the cause of initial fissuring at the surface of cartilage.

17

2.3 Results

2.3.1 Verification

The percentage difference between the results of the two analyses for the maxi-
mum shear stress (Figure 2.2) and the maximum vertical displacements were 2.0 % and
1.4 %, respectively. In addition, the percentage difference between the results of the two
analyses for the maximum effective stresses in the radial (Figure 2.3), axial, and hoop
directions were 0.81%, 2.28% and 0% respectively. Locations of maximum stresses were

the same in the two models.

2.3.2 Stress Analysis

The average of the maximum vertical displacement for the (1) model and (TI)
model were 0.107 mm and 0.085 mm, respectively. Based on these displacements and on
the aspect ratio data (Table 2.2) the isotropic model was the most compliant. The applied
load increased with indenter radius. There was a significant difference in the contours of
maximum shear among the three models (Figure 2.4). In the isotropic model, the location
of the maximum shear stress was at the interface with the subchondral bone, except for the
case of a 20 mm radius indenter, for which the maximum shear stress at the interface was
89 % of that at the center of contact (Table 2.3). In the anisotropic models, the maximum
shear was always highest at the surface, either in the center of contact (200, 50 and 20 mm
indenters) or at the edge of the area of contact (400 mm indenter). Furthermore, a two-
fold increase in the maximum shear stress value was observed in the anisotropic models

compared to the isotropic, except for the case of the 400 mm indenter in which maximum

18

shear stresses were of similar magnitude for all models. The maximum shear stress was
observed to be higher in the l-TI model versus the TI model. The ratio of maximum shear
suess to the maximum contact pressure increased with a decrease of the radius of the
spheres, reaching maxima of 50%, 26% and 25% for the l-TI, TI and I models respec-

tively.

Since the contact between the sphere and the cartilage was assumed to be friction-
less, the hoop stress on the surface of the solid skeleton was a principal stress, and equal to
the radial stress at the center of contact (Table 2.3). Hoop stress in the solid skeleton was
always compressive at the center of contact for the isotropic model. This stress became
tensile on the whole surface of contact in both anisotropic models, except for the 400 mm
indenter, where the hoop stress was negative at the center of contact. A small area of ten-
sile stress also existed at the edge of the contact area for the isotropic model. except for the

400 mm indenter.

2.4 Discussion

We hypothesized that stresses in the solid skeleton are responsible for fissures
observed in cartilage under blunt insult, as they are proportional to the forces transmitted
by the solid phase. Solid skeleton stress contours in transversely isotropic models of artic-
ular cartilage were shown to be significantly different from those in isotropic models.
Shear stress was maximum at the surface of the anisotropic models, in contrast to the iso-
tropic models, in which shear was in general maximum at the interface with the subchon-
dral bone. This location of maximum shear for the isotropic biphasic model was also

found by Ateshian et al. ( 1994) for the 400 mm radius indenter. Askew and Mow (1978)

19

showed a reduction, in the transversely isotropic models, of the normalized radial stress at
the center of the contacting area. Since the maximum shear stress is proportional to the
difference between the maximum and minimum principal stresses, the results in the analy-
sis of Askew and Mow (1978) would also indicate that the maximum shear at the center of
the contact area is higher in the transversely isotropic model than that of an isotropic
model. The compressive hoop stresses found at the center of contact for the isotropic
model were also consistent with results of numerous contact models of isotropic articular
cartilage under short time load (Armstrong, 1986; Eberhardt et al., 1991; Kelly and

O’Connor, 1996; Athesian et al., 1994; Wu et al., 1996).

Tensile hoop stresses in the solid skeleton were found at the surface in the center of
the contact zone for anisotropic models. In addition to high shear stresses, tensile stresses
have been suggested as a possible cause of surface injuries in articular cartilage under
impact load (Kelly and O’Connor, 1995; Eberhardt et al. 1990; Li et al., 1995). Recent
experiments undertaken by Atkinson et al. (1997) with spherical impactors on lateral tibial
plateaus of Flemish Giant rabbits show radial fissures in the center of the area of impact.
The presence of tensile hoop stresses, maximum at the center of the contact zone in the
solid skeleton of the anisotropic models of cartilage, may explain the radial fissures docu-
mented by Atkinson et al. (1997). In repeated impact experiments on distal ends of bovine
femora, Silyn-Roberts and Broom (1990) found fissures at the center of the impact zone in
23 of 32 specimens. After removing a thin layer of the surface, no surface damage was
induced, even after 150 cycles of impact. As removal of the surface makes the cartilage
more nearly isotropic, the observations of Silyn-Roberts and Broom (1990) would appear

to be consistent with our results since anisotropic models showed higher shear, and even

20

tensile stresses, on the surface of the solid skeleton as opposed to the isotropic model
which generally exhibited low shear stresses and only compressive stresses on the sur-

face.

The analysis of transversely isotropic models also showed, for high aspect ratios.
magnitudes of shear stress similar at the interface and the surface, and tensile stresses low
at the surface. For certain geometries, and depending on the anisotropy of the cartilage and
the solid volume fraction, injuries could be developed at the interface with bone instead of
the surface, as has been observed in many previous experiments (Thompson et a1. 1991;
Atkinson and Haut, 1995). In order to undertake a more detailed analysis of different
modes of failure observed in different species we would require additional information

regarding geometry, material properties and strength of the cartilage.

In conclusion, stress analysis of anisotropic models of articular cartilage seem to

offer better explanations for observed phenomena. Haut et al. (1995) reported an average

contact area over the lateral facet equal to 14.1 mm2 and an average thickness of 0.474
m. If the area is assumed to be circular, the aspect ratio is equal to 4.47, which falls
between the aspect ratios for the 50 mm and 200 mm diameter indenters. For both of these
cases, tensile hoop stresses and the highest shear stress were reported on the surface of the
cartilage. This could help explain the surface lesions observed during impact experiments

using the rabbit model.

A method has been proposed to study the stress field in the skeleton of a biphasic
transversely isotropic model of articular cartilage under short time load by analyzing an

equivalent incompressible elastic transversely isotropic model. The equivalent elastic

21

incompressible model can then be solved by standard analytical or numerical methods,
and the resulting stress field can be uncoupled into the fluid pressure and the stresses in the
solid skeleton. A standard finite element program can then be used to study the stress field

in a biphasic transversely isotropic model of cartilage at time zero.

The foregoing analysis was conducted using material properties for the human

humeral head for which the ratio of the in-plane Young’s modulus to the axial Young’s

modulus is high. This high value of 15:51) is needed to predict the initial response of bipha-

sic cartilage (Cohen et a1. 1993). The inability of the isotropic model to predict the early
response of cartilage has been previously reported (Armstrong et al., 1984; Mow et al.,
1989; Jurvelin et al. 1988). Based on the high ratios of the unrelaxed to the relaxed shear
moduli reported by Haut et al. (1995) for rabbit retropatellar cartilage one could presume
that only with a transversely isotropic model would it be possible to fit the early response
of the force-time relaxation curves. This suggests that the trends observed in the stress
analysis of transversely isotropic models may be applicable to the rabbit retropatellar car-

tilage. Nonetheless, further studies need to be done, to determine transversely isotropic

Properties for the rabbit retropatellar cartilage.

Another limitation of the analyses was the assumption of infinitesimal deforma-
tions. While no experimental deformation have been reported for impact experiments, the
comparative analysis would appear to be valid at least for the first part of an impact event
when deformations are still low. The stress analysis conducted here was also based on an
ideal geometry of a cylindrical disk of cartilage indented by a rigid sphere. However, it is

likely that the same trends would be found in diarthrodial joints with aspect ratios similar

22

to those studied in this analysis. While the transversely isotropic models considered here
are still rather crude approximations of cartilage, we believe they provide a significant
improvement over an isotropic model in the study of the mechanisms of blunt impact to

articular cartilage in a diarthrodial joint.

Table 2.1 Skeleton material properties from Cohen et al. (1993)

 

 

 

 

 

 

 

. Isotropic Transversely Isotropic
Material property Assumption Assumption
0.69 0.46
Egg) (MPa)
N/A 5.8
13:81) (MPa)
(s) 0018 0.0002
V12
(3) N/A 0.
V31
(33] N/A 0.37
Permeability (x10-15 m4/(N s) 3.0 5.1
Solid volume fraction 0.25 0.25

 

 

 

 

 

23

Table 2.2 Maximum pressures and aspect ratios

 

 

 

 

 

 

 

 

 

 

 

 

 

Maxium contact pressure Aspect ratio (radius of contact
Indenter .
size Force (MPa) area/thrckness)
(N)

(m) (I) (TI) (III) (I) (II) (I-TI)
4(1) 289.5 1.81 2.15 1.87 11.26 9.91 11.02
200 28.25 0.67 0.99 0.74 5.84 4.71 5.82

50 9.56 0.47 0.73 0.52 3.72 3.06 3.57
20 2.57 0.33 0.57 0.35 2.43 1.78 2.17

 

Table 2.3 Maximum shear and normal hoop stresses in the solid skeleton (solid

volume fraction = 0.25) at the locations A, B and C (Figure 2.4)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Indenter Shear (MPa) Hoop (MPa)
radius Model

(mm) A B C A B
400 11 0.23 0.24 0.21 0.15 +0.11
1-r1 0.24 0.25 0.21 -0.066 +0079
1 0.087 0.040 0.11 -0072 +0021
200 Ti 0.22 0.16 0.11 +0.15 +0.11
1-11 0.23 0.19 0.11 +0.28 +0093
1 0.083 0.0063 0.10 0044 +0035
50 TI 0.18 0.15 0.086 +0.14 +0.13
I-TI 0.22 0.19 0.10 +0.29 +0.14
1 0.081 0.050 0.072 -00010 +0036
20 TI 0.15 0.13 0.070 +0.13 +0.14
1-11 0.18 0.15 0.080 +0.34 +0.13

 

24

 

 

p = -(XO'kk

(II)
(I) p

 

 

 

 

 

 

\ on (III) “33 ' p

 

 

Figure 2.1 Stress decomposition in a biphasic element at time zero

25

0.145 MPa

Surface 0.148 MP3.

 

 

 

 

Center

0. 123 MPa
(a)

0.144 MPa Surface 0.151 MPa ‘1

 

 

.5-»—2I‘l \IA. \ ‘JLLHMS AL

 

 

Center

0.123 MFa

(b)

Figure 2.2 Maximum shear stress contours for (a) incompressible elastic and (b) biphasic

 

    

 

 

 

at time zero
“+0249 MPa Surface f -0 223 MPa
Q) gt, ,. .7
U ........ as
(a)
+0.247 MPa Surface -0 227 MPa
9 \“\\\ \“k 33/ /
E ' \ .
8 VMM

 

(b)

Figure 2.3 Contours of radial effective stresses for (a) incompressible elastic and (b)

biphasic at time zero

26

 

A (0.087 MPa) Surface B (0.040 MPa)

 

 

 

 

 

 

C (0.11MPa) 0‘
(a)

A (0-22 MP3) Surface B (0 16 MPa)

 

 

 

 

 

Center

 

 

C (0.11 MPa)

0))

A (0.23 MPa) Surface B (0.19 MPa)

 

 

 

Center

 

 

 

 

 

 

(C)

Figure 2.4 Contours of maximum shear stress for the 200 mm indenter and (a) I, (b) TI and

(c) I-TI models

27

CHAPTER 3

A METHOD TO DETERMINE IN SITU MATERIAL PROPERTIES OF
TRANSVERSELY ISOTROPIC BIPHASIC CARTILAGE

A method is presented to determine the material properties of biphasic cartilage
based on a transversely isotropic model. Given stress relaxation curves obtained from tests
performed with two cylindrical, non-porous, plane-ended indenters of different radius,
four of the elastic constants of the solid skeleton and two permeabilities (parallel and per-
pendicular to the plane of isotropy) can be obtained. The fifth elastic constant, the in-plane
Poisson’s ration, was assumed to be zero. The calculated transversely isotropic properties
provided a good fit for the entire range of the relaxation curves, including the early time
responses of the tissue. Transversely isotropic properties are reported for the rabbit retro-

patellar cartilage.
3.1 Numerical Method
3.1.1 Determination of candidate elastic constants

The first part of the force-time curve (Figure 3.1) obtained in an indentation relax-

ation experiment on articular cartilage is represented by the peak force (fpe) attained dur-

ing rapid indentation, and assumed to occur without a significant flow of fluid in the tissue.

During this stage the articular cartilage can be considered elastic and incompressible

28

(Eberhardt et al., 1990) with effective elastic constants that can be found in terms of the
elastic constants of the solid skeleton. The last part of the relaxation curve is the equilib-

rium response in which all the extemal force (feq) is supported by the solid skeleton. In

this stage the elastic properties of the cartilage are those of the solid skeleton alone.

If we assume isotropic cartilage, the effective Poisson’s ratio and Young’s modulus
are v = 0.5 and E = 1.5 E“) / (1+v(s)), respectively (Armstrong et al., 1984), where Em

and V“) are the properties of the solid skeleton. If we know the cartilage thickness we can

use the elastic solution of this problem (Hayes et al., 1972) to undertake a trial and error
procedure to find the set (E‘s), vm) that uniquely matches the experimental point (feq, fpe).

Another approach consists of calculating theoretical (feq, fpe) points corresponding to

(ES), vm) sets and then drawing those points in a cartesian coordinate system (fe fpe), as

q,
shown in Figure 3.2. Thus, given an experimental (feq, fpe) point, we can estimate by inter-
polation the set (E‘s), vm) corresponding to that point. One problem is that we do not
know in advance the thickness of the tissue layer, but plots can be generated for various
thicknesses. Then, depending on the actual thickness of the cartilage at a particular loca-
tion, one can interpolate between the plots to determine the elastic constants at a given test

location.
Typically in a relaxation test on articular cartilage an isotropic model can only
match an experimental (feq, fpe) point if the Poisson’s ratio is negative. A transversely iso-

tropic biphasic model, on the other hand, is characterized by five elastic constants for the

solid skeleton. These constants are the Young’s modulus E1? in the plane of isotropy (par-

29

allel to the surface of articular cartilage), the Young’s modulus normal to the plane of isot-

ropy (i.e. the axis of loading) 5‘32} , the shear modulus 613, the out-of-plane Poisson’s ratio

VS) , and the in-plane Poisson’s ratiov‘lsz) .

The effective properties E11, E33, V3] and v12 of transversely isotropic biphasic

articular cartilage can be given in terms of the properties of the solid skeleton by equations

(2.13), (2.20) and (2.21) presented in Chapter 2. The shear modulus G 13 of the effective

elastic material is equal to the shear modulus of the solid skeleton, since the fluid is

assumed inviscid and supports only hydrostatic loads.

Four elastic constants of the solid skeleton were determined in this numerical pro-
cedure. Parametric studies with a finite element model revealed that the out-of-plane Pois-
son’s ratio had a substantially greater influence on the equilibrium and peak forces than

did the in-plane Poisson’s ratio. In one parameter study, for example, it was assumed that

E‘lsl)=5.0 MPa, Egg) =1.5 MPa and G13=0.5 MPa. For a fixed 48341.1, changes in v‘lsz)
from 0.0 to 0.2, 0.3, 0.4 and 0.5, caused changes of 0.51%, 5.10%, 11.05% and18.4%,
respectively, in the peak force and changes of 0.84%, 1.39%, 1.95% and 3.06%, respec-
tively, in the equilibrium force. Conversely, for vi? fixed at 0.1, changes in (1‘3? from 0.0

to 0.2 caused changes in the peak and equilibrium forces of 16.4% and 22.6), respectively.
Similar effects were observed using other values for the elastic constants. Hence, since it
is not possible to determine the in-plane Poisson’s ratio for this type of experiment we

assumed a value of 0.0, based on the previous study by Cohen et a1. (5) which documented

a value of 0.0002 :1: 0.0005 (n=9).

30

In order to match a (f , fpe) point, we followed a procedure similar to that

described for isotropic cartilage. However, instead of a unique solution, a spectrum of

combinations of elastic constants could be found to match a given experimental (feq, fpe)

point. A computation scheme was therefore developed to determine this spectrum of com-

binations of constants. This scheme is as follows.

For given values ongsl) andvgsl), corresponding combinations of 13g? and G13

were plotted on (qu, f e) planes. For example, Figure 3.3 depicts such plot assuming E1”? =

5 MPa and vgj’ =0.10, and a cartilage thickness equal to 0.5 mm. Given an experimental

(feq, fpe) point, a linear interpolation scheme was used on each plane in order to determine

the elastic pairs (Egg), G13) associated to the point. All plots have been generated for

thicknesses of 0.5 mm and 1.0 mm because this range coincides with the thickness of rab-
bit retropatellar cartilage being studied in our laboratory. Other ranges of values could be
plotted in a similar fashion. As was done in the isotropic case, combinations of elastic con-
stants for each experiment were determined by linear interpolation based on the actual car-
tilage thickness at the test location. A finite element linear elastic model consisting of 560
8-noded axisymmetric elements (AN SYS 5.3, SAS IP Inc., PA, USA) was used to gener-

ate these plots.
Some trends were observed in the sets of elastic constant satisfying experimental

(f , fpe) points. First, only with relatively low values of the Poisson’s ratio vgj’ (less than

0.25) was it possible to match experimental (feq, fpe) points. Secondly, since Egg) was

31

largely determined by the equilibrium force, this material constant had no substantial vari-

ations for a given (qu, fpe). Finally, it was observed that for a given (feq, fpe) point high

values of 5131) implied low values of 613 and vice versa.

3.1.2 Determination of candidate isotropic permeabilities

Given these combinations of elastic constants, we then intended to fit the interme-
diate part of the indentation curve. The shape of this curve depends on the permeability

coefficient, the elastic constants and the geometry of the problem (McCutchen, 1982).

No analytical solution has been developed for the indentation problem of a trans-
versely isotropic biphasic material. Thus the finite element program Abaqus 5.5 (Hibbit,
Karlsson & Sorensen, Inc., RI, USA) was employed to generate force-time curves for
indentation of a material with poroelastic properties. An axisymmetric mesh of 328
CAX8P elements was used for the analysis. A displacement was applied at the top nodes
of the model in the region of contact. A free flow condition was assumed at the surface of
the cartilage outside the zone of contact with the indenter, and an impermeable condition

was assumed everywhere else along the boundary.

If we assume at this time that the permeability of the cartilage is isotropic, we can
input a given combination of elastic constants satisfying an experimental (feq, fpe) point,
and run the program for various values of the permeability coefficient until we find the
best fit of the experimental curve obtained with a 1 mm diameter indenter (which will be

subsequently labeled 1 mm curve). In this way, we can associate an isotropic permeability

(11%) to each combination of elastic constants. Now, with various combinations of elastic

32

constants and permeabilities, we can fit the entire 1 mm curve. It was observed that lower

14‘}, coefficients were associated with sets of elastic constants containing higher E‘lsl) val-

ues, and higher k‘ll}, coefficients were associated with sets of elastic constants containing

lower E131) values.
3.1.3 Determination of actual elastic constants

It was determined that different combinations of elastic properties and permeabili-
ties lead to nearly identical response curves, see, for example, Figure 3.4 (a) where curves
from two different sets are compared. Therefore, in order to determine a unique set of elas-
tic constants and permeabilities, we seek to fit another force-time curve obtained with a
1.5 mm plane-ended solid cylindrical indenter (which will be subsequently labeled 1.5 mm
curve). In Figure 3.4 (b) we can observe differences in the theoretical curves obtained for

the 1.5 mm indenter with combinations that equally fit the 1 mm curve.

The 1.5 mm curve was fit following a procedure similar to that used for the 1 mm
curve, i.e. given a combination of elastic constants, the program Abaqus was run for vari-
ous values of the permeability coefficient until the best fit of the experimental curve was

found. Based on this procedure one combination was selected that best fit the 1.5 mm

curve. Another permeability coefficient ( k‘li; ) was also obtained in this procedure.

All combinations of elastic constants satisfying exactly a (f , fpe) point for the 1
mm curve were able to match a (feq, fpe) point for the 1.5 mm curve to within 5%. It was

felt that an exact match could not be achieved because of the nature of the indentation tests

33

employed. Contact between the plane-ended surface of a cylindrical indenter and the sur-
face of rabbit retropatellar cartilage might not be perfect, considering the convex small
patella of the rabbit. Since the elastic solution for a perfectly plane cartilage layer indented
by a plane-ended cylinder exhibits a high contact pressure at the perimeter of the circle of
contact, the error in the prediction of an experimental force by using this solution would

be higher for large indenters.

3.1.4 Estimation of anisotropic permeabilities

It was observed that for a given set of properties that fit both indenter tests. kilo?)

was different that 1611}, . One possible explanation for this difference was that permeabili-

ties were anisotropic, with a permeability coefficient k1, associated with the plane of isot-
ropy, and another permeability coefficient k3, associated with the axis of loading. The
finite element analysis indicated that the direction of the flow was approximately radial
(direction 1) under the solid indenter and predominantly axial (direction 3) near the bor-
der. This meant that the fluid flow was governed mainly by the permeability k1 under the
solid indenter and by the permeability k3 near the tip. In order to estimate k1 and k3 a sim-
ple model was developed (Figure 3.5). Using Darcy’s law, and assuming a steady flow
governed only by the permeability k1 under the indenter and by the permeability k3 near

the border, anisotropic permeabilities were related to the average isotr0pic permeabilities

by using the formula presented by Terzaghi (1943)

_ -5. 1__§
"1.1“ (3.1)

141:, )9

34

__ = 11. 1_—fl
. k1+ k3 (3.2)

The geometrical parameters i and C were determined as follows. Elastic constants and iso-
tropic permeabilities were calculated as explained before for three experiments on rabbit
patella cartilage to be described later. In order to provide a higher sensitivity, three experi-

ments were chosen for which the effective isotropic permeabilities kg?) and kg]; were

significantly different. Then, the program Abaqus was run considering various sets of
anisotropic permeabilities k1 and k3 calculated with equations (3.1) and (3.2) for different
values of g and C. Then, the force-time curves from Abaqus where compared with those
obtained in experimentally. With the values i = 0.5, Q = 0.99, the calculated anisotropic

permeabilities produced the best fit of the experimental curves.

3.1.5 Automation of the algorithm

An application of the above method to routine experiments would be computation-
ally prohibitive. In order to automate the procedure we sought to define an analytical curve
to represent the curves obtained with Abaqus. After a trial and error procedure, the curve

defined by the expression

 

_ 1 — or 01
f“) — feq + (fpe — feq)[m + 1+ a - sinh(ct)] (3'3)

was shown to follow the force-time curves from Abaqus very well. In equation (2.3), or
and a are parameters dependent on the elastic constants and c is dependent on the isotropic

permeability coefficient. The force is equal to fpe at t = 0 and it is equal to feq when t tends

to infinity. Comparison of the master curve with the output from Abaqus is presented in

35

Figure 2.4 for one set of elastic constants (fixed or, a) and two permeabilities. Correlation

coefficients (r2) higher than 0.997 were obtained for all sets of force-time curves from

Abaqus calibrated with the master curve.

Considering than in this problem the effective length for fluid flow is proportional
to the radius of the indenter and not to the thickness of the cartilage (McCutchen, 1982).
the solution was rather insensitive to cartilage thickness. Thus the calibration of the output
of Abaqus with the master curve was undertaken considering a fixed thickness of 0.75
mm, which is approximately the media thickness for rabbit retropatellar cartilage used in

this particular study.

From all combinations of elastic properties satisfying the (fe fpe) point for the 1.0

(1’
curve, the constants ES) and G13 had the most variation between sets. For example, a typ-

ical point (feq = 0.18 N, fpe = 1.04 N) could be matched with combinations containing val-

ues of E331) from 2.0 to 5.0 MPa and associated values of G13 from 0.87 to 0.10 MPa,

respectively. Values of Hg?) and G13 outside these ranges would not match the (feq, fpe)

point for the 1.0 mm curve. The test undertaken with the 1.5 mm indenter was intended to

provide additional information for further refinement of material properties. Since the

shape of the 1.5 mm curve was mainly affected by the elastic constant E351) , no significant

differences were observed in the shape of the curve for material sets containing similar

values of E‘fl’ , e.g. E3? :30 and 3.2 MPa. Thus, only three material property sets were

considered form the combination that fit equally well the 1 mm curve. The three sets con-

36

tained the highest, the median and the lowest values of the elastic constant E351) , which, as

noted before, were associated with the lowest, the median and the highest values of G13

and with a relatively constant value of Egg). Lower values of Va? were also associated

with lower values of E181) . Even though this was not a precise determination of the mate-

rial set which ‘best’ fit the l.5 mm curve, it allowed a selection from among the sets of

elastic constants that fit equally well the 1 mm curve.

All details about the actual implementation of the procedure in the program

tiprop.m, written using the Matlab code can be found in Appendix A.

3.2 Experimental method

Indentation relaxation experiments were undertaken on retropatellar cartilage of
mature, Flemish Giant rabbits following a protocol described by Newberry et a1. (1997).
In brief, the patellae were excised, placed in phosphate buffered saline (pH 7.2) and
mounted in a fixture with the retropatellar surface facing an indenter probe that was
attached to a stepper micrometer (model M-168.30, PI Physik Instrumente, Auburn, MA).
First, a plane-ended, non-porous, cylindrical indenter of 1.0 mm diameter was pressed 0.1
mm into the cartilage within 60 milliseconds and maintained for 150 seconds at three loca-
tions (proximal to distal) along the center of the lateral facet. The indentation force from a
transducer (model JP25; Data Instrument, Acton, MA, U.S.A.) (25 lb [11.3 kg] capacity),
attached between the' probe and the stepper, was amplified (model SG71, Validyne,
Northridge, CA, U.S.A.) and collected at 1000 Hz during the first second and at 20 Hz

during the remainder period by a personal computer with an analog-to-digital board. Next,

37

after a period of 360 s, the same procedure was repeated in the same location with a plane-
ended, non-porous, cylindrical indenter of 1.5 mm diameter. Finally, tissue thickness was
measured with a penetrating needle probe. The elastic constants and permeabilities deter-

mined from indentations at the three locations were averaged.
3.3 Results

The procedure was applied to determine the elastic properties and permeabilities
of patellar cartilage (Table 3.1). The time for the tiprop.m program to perform each esti-
mation was between 20-30 8 on a Sun Ultra 1 workstation (Sun Microsystems, Inc., CA,
USA). To validate the predictions of tiprop.m program, the program Abaqus was run a
posteriori. The best fit predicted with the program Abaqus always corresponded to the best
fit predicted with the program tiprop.m. For this study on rabbit retropatellar cartilage, the

best fit was generally obtained with the set of elastic constant containing the maximum

value of E‘fl’. Comparison of the force-time curves from Abaqus and experiments indi-

cated a good fitting of both curves (Figure 3.7).

3.4 Discussion

An automated numerical procedure was developed to estimate four elastic con-
stants and two permeabilities for biphasic cartilage based on a transversely isotropic
model. “0th this procedure, implemented in a Matlab program, it was possible to deter-
mine one set of elastic constants and permeabilities that simulated the high ratios of peak
to equilibrium load typically found in relaxation indentation experiments on rabbit retro-

patellar articular cartilage. Moreover, the transversely isotropic model provided a excel-

38

lent curve fit of the force-time relaxation curve for the 1 mm and 1.5 mm diameter

indenters.

Determination of the elastic constants was based on matching peak and equilib-
rium forces in the relaxation experiment with a 1 mm diameter indenter. This required

generation of points in the (fpe, feq) space by running a finite element elastic model numer-

ous times considering different combinations of elastic constants and thicknesses. This
demanded a lot of computational time. On the other hand, after all necessary tables were
generated, the matching of the peak and equilibrium load from an actual experiment
required only a few seconds. It was found that peak and equilibrium forces could be
matched with various combinations of elastic constants. From these combinations, three
representative sets were chosen containing high, medium and low values of the in-plane
(transverse) Young’s modulus. The final elastic set was subsequently selected as the one

best matching the shapes of the force-time curves for the 1.5 mm diameter indenter.

A master curve was found that followed very closely the relaxation curves from an
Abaqus code modeling the indentation problem. The mathematical curve utilized three
parameters, two related to the elastic constants and one related to permeability for the best
fit. With the master curve it was possible to implement efficiently and automatically the

calculation of permeabilities and elastic constants.

Using a transversely isotrOpic model, Cohen et al. (1993) fit force-time creep
curves for porous indenters employing the methodology developed by Athanasiou et al.
(1991), which was based on a finite element analysis and a non-linear optimization tech-

nique. The high computational cost of the methodology has impeded its widespread use

39

(Athanasiou, 1991). On the other hand, the methodology proposed here could be used on a

routine basis, once the points in the (fpe, ) space were determined for various combina-

feq

tions of elastic constants and thicknesses of cartilage.

One limitation of the method is that the final set of elastic constants was chosen
from a discrete number of combinations (three) that fit equally well the 1 mm curve. The
estimation of anisotropic permeabilities was based on a very simple model which
described very roughly the flow of water inside the cartilage. A posteriori verifications
with the program Abaqus revealed that the force-time relaxation curves obtained with the

anisotropic permeabilities for the two indenters followed closely the experimental curves.

The Young’s modulus E251) (Table 3.1) for rabbit retropatellar cartilage (5.85 MPa)

compared very well to that reported by Cohen et al. (1993) for human humeral cartilage

(5.6 MPa). The shear modulus G13 (0.18 MPa) was between 0.14 MPa measured by

Khalsa et al. (1997) in the bovine femur and the equilibrium modulus of 0.23 MPa

reported by Zhu et al. (1986) for human patellar cartilage. The Poisson’s ratio value vffl’
(Table 3.1) was in the range of isotropic values reported in the literature (Jurvelin et al.,

1997; Whipple et al., 1985; Athanasiou et al., 1991). Only with small vgj’ values (less that

0.25) was it possible to satisfy both peak and equilibrium forces. This appeared to be due
to the high ratios of peak to equilibrium force (generally greater than 4.5) documented in

this study and often by others. Theoretically, the higher the Poisson’s ratio, the lower the

ratio of fpe to feq. For instance, if E331) = 5 .0 MPa, Egg) = 1.5 MPa, and G13 2 0.7, the ratio

of fpe to feq was 2.13 for vff,’ = 0.30, contrasting with a ratio of fpe to feq equal to 5.7 for

40

vgj’ = 0.10. This dependence could also be inferred from equation (2.21) which predicts

no relaxation for vgsl) =0.5, since the effective Young’s modulus in the axis of loading is

equal to the modulus of the solid skeleton. Numerical experiments also revealed that the
indentation forces were not sensitive to variations of the in-plane Poisson’s ratio, which

was assumed to be zero in this analysis based on previous data by Cohen et al. (1993).

The Young modulus E22) (Table 3.1) was slightly higher than the range of values

reported for isotropic cartilage in the literature (Athanasiou et al., 1991). The difference
may, in part, be due to the relatively short time (150 s) we used to measure the equilibrium
force. Longer duration tests may be expected to yield values approximately 10% lower
(Haut et al., 1995). Another part of the difference could be explained by the amount of
indentation (0.1 mm) used in this study, which produced an average axial strain of 17%

(considering the average thickness reported in Table I). Since the method developed here

was based on infinitesimal analysis. the modulus Egg) could have been overestimated. We

have developed an algorithm for a nonlinear, transversely isotropic biphasic model that

will be reported in Chapter 4.

The average value of the anisotropic permeabilities k1 and k3 (Table 3.1), i.e.

3.67x10'15 m4/(N.s), was similar to the isotropic permeability of 3.84x10‘15 m4/(N.s)
reported by Mow et a1. (1984) for rabbit patellar cartilage. Permeability k1 was signifi-
cantly higher than k3 (paired test, p<0.05). Since the permeability increases with compac-

tion of the tissue (Mow et al., 1984), and the permeability was assumed to be constant in

this model, the measured permeabilities can be considered apparent, as explained by Mow

41

et al. (1984).

The in-plane Poisson’s ratio in this study was assumed to be zero, equal to the
value reported by Cohen et al. (1993). Parametric studies indicated that peak and equilib-
rium forces for this problem were not sensitive to variations of the in-plane Poisson’s ratio
for a range of low values, as expected in a porous material like the solid skeleton of carti—
lage. Therefore, even if the in-plane Poisson’s ratio were as high as 0.3, the prediction of
the other material properties would still be valid. On the other hand, significant differences
in the prediction of the other material properties would be expected if the in-plane Pois-

son’s ratio actually approached 0.5. Hence this is a limitation of the current methodology.

While micromechanical analyses indicate that articular cartilage is anisotropic,
insufficient evidence exists to verify a transversely isotropic model. Nonetheless, the
transversely isotropic model is the simplest anisotropic model able to predict the high
ratios of peak to equilibrium forces observed in many indentation relaxation experiments
on articular cartilage, as suggested by Cohen et al. (1993). Future studies could be directed
towards more microstructurally complex models, but this may not be useful in the reduc-

tion of large volumes of experimental data. a major advantage of the current model.

42

Table 3.1 Thansversely isotropic properties of rabbit retropatellar cartilage"
kl k3
5‘1? 13;? G13 v‘fl’ (m4/(N s) (m4/(N s)
x1015 x1015

0.58 5.85 1.07 0.18 0.092 4.79 2.54

(0.059) (2.00) (0.30) (0.047) (0.029) (2.26) (1.45)

 

Thickness
(mm

 

 

 

 

 

 

 

 

 

* n = 10, mean (standard deviation)

 

 

 

 

 

 

 

2.5 I l 1 I l I I
2 f -
1'5 Peak force 4
A f
a PC
8
LL 1 .. -
o s - t. .
Equilibrium force feq
o l * l l l l l l
0 20 40 60 80 100 120 140 160
Time (s)

Figure 3.1 Typical force-time relaxation curve from an indentation experiment

43

 

 

Peak force, fpe (N)

 

 

 

 

 

.3
0.12 0.14 0 16 0.18 02 022 024 0.26

Equilibrium force, feq (N)

Figure 3.2 Points in a typical fpe, feq plane for an isotropic biphasic layer of 0.5 mm. The
arrows represent a typical experimental point for which the Poisson’s ratio was equal to -

0.1

44

 

 

Peak force, fpc (N)

 

 

 

 

 

V?
a
II .
E
Q) :
' 1 1 l 1 l l
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Equilibrium force, feq (N)

Figure 3.3 Points in a typical fpe, feq plane where 5‘11): 5.0 MPa, 8;: = 0.10, and thick-

ness = 0.5 mm. The indicated linearly interpolated point corresponded to Egg) = 1.4 MPa

and 613 = 0.65 MPa

45

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

l l T l 1
E33“) 013 V31“) 1‘10i . _
2 (MPa) (MPa) 10'15m4/(Ns)
8 1.11 0.63 0.050 2.9 1“
E 1.32 0.15 0.025 2.0
23‘ 5* r--i---+o
250 300 350 400 450 500
Tlme(s)
(a)
H I I I I I I I
. 511“) 533“) Gl3 V31“) k1.0i
1Q (MPa) (MPa) (MPa) 1015m4/(Ns)~
‘; 2.20 1.1 1 0.63 0.050 2.9
09 6—0. 4.40 1.32 0.15 0.025 2.0 ~

Force (N)

 

 

 

 

0 4 i i i i i i 1'
‘ 20 4o 60 so 100 120 140 160
(b)

Figure 3.4 Force-time relaxation theoretical curves obtained for two sets of material prop-

erties for the (a) 1 mm diameter indenter and (b) 1.5 mm diameter indenter

46

 

 

 

 

(1) Physical model

Surface

k31 1
'1
/

Interior —> ,

———”

(b) Anisotropic permeabilities

1 Surface

. 1(1.0 l 1
i ,+

Interior —> ,/

———’

(c) Average isotropic permeability

Figure 3.5 Simple model to estimate anisotropic permeabilities

47

 

 

0.55 I I 1 I I L r L 4

 

 

 

 

 

 

 

 

 

 

 

. Output from Abaqus
O 0 Master curve
0. .- . . . H
t
0.45
g 0.4
0
8 035
12
0, . .. ..f,
0.25 ..
0.2
0
Time (s)
(a) Permeability = 0.5x10'15 m4/(N s), c = 0.005
05g I I I I . e .
c) Output from Abaqus
d) 9 ’ l ’ O 0 Master curve

 

 

 

04Q5........ ...... . .. I.... .. .....'....... ...I_...........:.................................... -(

A
v

 

 

 

 

 

I“; 1 1V 1 1V 1 V
40 60 80 100 120 140 160

Time (s)

(b) Permeability = 25x1015 m4/(N s), c = 0.23
Figure 3.6 Comparison of the master curves with the force—time curves from Abaqus, a=4

andor=l

48

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.2 T l r 1 I I
' ‘ ' Experimental ;
, .. O 0 Curve from Abaqus _
A
5 _
53
C
[L ..
E '29
“a 22 43 in , 8‘0 .30 1.120 do 160
T1me(s)
(a) Curves for the 1 mm indenter
4 r l r r I I I
I ' Experimental
35’ " O 0 Curve from Abaqus
3............?....,......'............i............g............§............§............:~...........l
-.- . ...J
A
Z
v y —
Q
g .
m . . . . l
o- ' , ...L-
_0_5 l 1 1 l l 1 l
o 20 40 so so 100 120 140 160

 

 

Time (s)
(b) Curves for the 1.5 mm indenter

Figure 3.7 Comparison of typical experimental curves with the curves from Abaqus run

with the constants and anisotropic permeabilities predicted with the program tiprop.m

49

CHAPTER 4

ESTIMATION OF IN SITU ELASTIC PROPERTIES OF BIPHASIC CARTILAGE
BASED ON A TRANSVERSELY ISOTROPIC HYPO-ELASTIC MODEL

A method was developed to determine the nonlinear transversely isotropic proper-
ties of articular cartilage by fitting two force-displacement curves obtained from large
deformation indentation relaxation tests on cartilage using a non-porous spherical
indenter. The first curve consisted of the peak forces obtained at the initial stage of rapid
indentation versus various associated displacements and, the second curve consisted of the
equilibrium forces versus the associated displacements. We modeled the solid skeleton as
a transversely isotropic hypo-elastic material. A commercial finite element code was
employed to solve the problem of a layer indented by a rigid sphere. The components of
the hypo-elasticity tensor were made dependent on deformation according to the varia-
tions defined by a transversely isotropic hyperelastic formulation. Incompressibility was
assumed to model the initial stage of rapid loading. The method was applied to indentation
data from rabbit retropatellar cartilage to document the non-linear elastic properties for

impact modeling.

4.1 Theoretical Analysis

The theory of hypo-elastic material behavior was formulated by Truesdell (1955).

Later, Noll (1955) showed that hypo-elasticity is the most general form of elasticity.

50

Under a hypo-elastic model an objective stress rate 8; is dependent on the rate of defor-

mation dkl by the equation

Sij = Hijkldkl (4-1)

where Him is the hypo-elasticity tensor. An objective stress rate is invariant to rigid body

rotation. Many objective stress rates have been proposed, e.g. Truesdell, Zaremba-Jauman

and logarithmic (Xiao et al., 1997).

Bernstein (1960) presented the conditions under which the hypo-elastic equations
are equivalent to those of Cauchy elastic materials. This means that equation (4.1) can be
integrated in order to obtain relations between stress and deformation. For a grade zero
hypo-elastic model (for which the components of the hypo-elastic tensor are constant)

under simple extension, integration of equation (4.1) yields (see details in Appendix B.1)

o‘fl’ = Hlmmotln H11221n(}t2)+ H1133ln(7\.3) (4.2)
0‘2? = H11221n(?tl)+ H22221n(12)+ H22331not3) (4.3)
6‘3? = H11331n(x1)+ H2233ln(}\.2) + H3333ln0t3) (4.4)
where 6:?) is the normal component of the Cauchy stress in direction i, and 7&1, X2 and K3
are the stretches in directions 1, 2 and 3, respectively. The components GE?) represent

effective stress in the solid skeleton. The total stress in the tissue is obtained (Mow et al.,

1980) by adding the fluid pressure p, i.e

_ (e)
oij — —p8ij +0ij (4.5)

51

The solid skeleton was assumed to be transversely isotropic, with the plane of isot-
ropy (plane 1—2) parallel to the surface of the cartilage and perpendicular to the axis of

loading (axis 3). Under this symmetry H1133 = H2233 and H1111 = H2222. The engineering
constants of the solid skeleton are the Young’s modulus 5‘1? in the plane of isotropy. the

Young’s modulus Egg) normal to the plane of isotropy, the shear modulus G13, the out-of-

plane Poisson’s ratio vg? , and the in-plane Poisson’s ratio VS) .

During the equilibrium stage, pressure in the fluid is zero and the solid skeleton

supports all of the load. For the case of uniaxial loading in the 3 direction, k1 = 12 and

equation (4.2) can be used to obtain it] in terms of 13 as follows

_ H1133
1mm) _ —H11“+H11221n(7t3). (4.6)

 

Differentiation of equation (4.6) yields

 

 

 

H1133 = _7‘3d}‘1 (4 7)
H1111 +H1122 ladle
By definition, the out of plane Poisson’s ratio vgsl) is
(s) : _[‘17f_i] = H1133 (4 8)
31 d7V3 Mia->1 H1111+H1122 '
and equation (4.7) becomes
A d?»
(S) _ _ 3 1
v31 — K1978. (4.9)

52

Substitution of equation (4.6) into equation (4.4) yields

 

2
2H
(8) 1133
533 - 633 = [H3333“H1111+Hlm]ln("3)'

Differentiation of equation (4.10) yields

 

2
H _ 2H1133 = A dc533
3333 H1111+ H1122 3 (17‘3

So that the infinitesimal Young’s modulus in the direction of loading is given by

 

2

(s) _ [(1533] _ H 2H1133
— — - 3333-
13-91 H1111+H1122

(4.10)

(4.11)

(4.12)

(4.13)

By solving problems similar to the foregoing it is possible to obtain additional

equations relating the five independent components of the transversely isotropic zero

grade hypo-elasticity tensor and the engineering constants of the solid skeleton. Inversion

of those equations yields the relations

(s) (s)_ EMV <s>2
Elll:E33 E11V3l ]

H1111 = s 3 s2
<1le’)<E"<1vlz’>-2Eii’v§1’>

 

E(S)2[1_ v(182)]

(s) (s) (s)v (8)2
E33(l V12) 213“" 31

 

H3333 =

53

(4.14)

(4.15)

(S) (S) (S) E(S) (8)2
JE11[E33V12 E11V31 1

 

 

H 22 = (4.16)
(S) (S) Vls)
E11E33v 31
H1133 ‘ EMU v(1s))_ 2Emv (5)2 (4.17)
—. 11 V31
H1313 = G13 (4.18)

During rapid loading, articular cartilage can be considered as an incompressible

elastic material (Eberhardt et al.,1990). Under this condition
1131213 = 1. (4.19)

For the case of uniaxial rapid loading in the 3 direction, 711 = M, and equations

(4.2) and (4.5) can be used to determine the pressure necessary to ensure stress free lateral

surfaces as follows
P = (H1111+ H1122)1“0‘-1)+ H1133ln013). (4-20)

Replacement of this pressure into equation (4.5) and use of equation (4.19) yields

 

H +H
“‘1 2 “22]111013) 21533111013) (4.21)

0‘33 = [H3333 — 2“1133 +

where E33 represents the effective Young’s modulus of cartilage under rapid loading. It

can then be shown that replacement of equations (4.14)-(4.17) into equation (4.21) yields

E _E‘,;’[1 4v‘3j’+2(1—v‘;’)(E§;’/E‘fi’)] 422
33 ' v11s) 1s) 1s) (sl)2 (' )
21(— 23)(E3-/E11) 4v3

 

which is equal to equation (2.21) obtained in Chapter 2 using the assumption of infinitesi-

54

mal deformations.

A grade zero transversely isotropic hypo-elastic model is not capable of describing
the stiffening behavior of articular cartilage. Therefore, in order to incorporate this effect,

the components of Him were made dependent on deformation. Since hyperelastic formu-

lations have already been employed to model articular cartilage, our hypo-elastic model
was made consistent with a hyperelastic formulation for the cases of uniaxial load and
pure shear. In other words, we solved two uniaxial problems (in the axis of loading and in
the plane of isotropy) using a hyperelastic formulation and calculated, at each deformation

( )

level, the Young’s modulus Egg) and Poisson’s ratiov;1 using equations (4.9) and (4.13).

Similar equations were used for the calculation of tag? andv‘lsz) . The variation of G13 with

deformation was obtained by solving the out-of-plane pure shear problem. The hypo-elas-
tic tensor was then calculated using equations (4.14)-(4.18) at each deformation level.

Since 13:51), vgj’, Egg) , vgi’ and 013 are now dependent on deformation, they will be

called the engineering constants when they are referred to the undeformed configuration

and the engineering parameters when they are referred to other deformed configurations.

In a hyperelastic formulation the Cauchy stress is obtained from (Holmes, 1986)

aw
oi). = —p5ij+ (2/J)Fik(a—(:;)ij (4.23)

where Fik is the deformation gradient, Ckm is the right Cauchy-Green deformation tensor,

J is the determinant of the deformation gradient and W is a strain-energy function govem-

ing the behavior of the hyperelastic material. We chose the strain-energy function pro-

55

posed by Almeida et a1. (1997 ) and defined as follows:

wal, 5,13,14,15) = a0 Exp[a1(11— 3) + 11,112 — 3) + a3(I1 — 3 )2+ (4.24)
1,,(14 — 1) + 115(14 -1)2 + 116(11— 3)(14 — 1) + 117(15 — 1)]/1§

where 11, 12, 13, are the strain invariants of the right Cauchy-Green deformation tensor C,

a0, a1, a2, a3, a4, a5, a6, a7 and n are coefficients independent of deformation, I4 and 15 are

strain invariants defined as 14 = trace(M C), 15 = trace (M C2), and M is the structural ten-

sor defining the preferred direction of the material in the reference configuration (direction

3 in this case). To ensure a stress free undeformed configuration, it is required that n = a1 +
2 a2 and a7 = -a4 (Almeida et al. 1997). Therefore, only seven coefficients of the strain-

energy function are independent.

The five infinitesimal engineering constants of the solid skeleton were expressed in
terms of the seven coefficients of the strain-energy function through equations (B.3)—(B .7)
presented in Appendix B. Analyses of the uniaxial and simple shear problems according
to the hyperelastic formulation were possible using the program Mathematica (Wolfram

Research Inc., Champaign, IL).

Again, for the case of rapid loading, articular cartilage was assumed to be incom-
pressible. The effective infinitesimal engineering constants (i.e. those resulting from the
combined action of solid skeleton and water during rapid loading) were also related to the
coefficients of the strain-energy function as given by equations (B8) and (B9) of Appen-

dix B.

In order to fit the experimental force—displacement curves, two coefficients of the

strain-energy function were adjusted in our model. The other five coefficients could then

56

be found in terms of the infinitesimal elastic constants using equations (A.2.3)-(A.2.7). In
one parametric study, the uniaxial stress-stretch equilibrium and rapid loading curves were

plotted (Figures 1 and 2) for different values of the coefficients a1 and 11. These two param-

eters were varied in order to fit simultaneously the peak and equilibrium experimental

curves obtained with the spherical indenter.

4.2 Methods

4.2.1 Experimental

Indentation relaxation experiments were undertaken on cartilage of five fresh left
patellae of mature, Flemish Giant rabbits. The patellae were excised, placed in phosphate
buffered saline (pH 7.2) and mounted in a fixture with the retropatellar surface facing an
indenter probe that was attached to a stepper micrometer (model M-168.30, PI Physik
Instrumente, Auburn, MA). The following procedure was followed at each of four loca-
tions on each patellae (two along the center of each facet). First, a plane-ended, non-
porous, cylindrical indenter of 1.0 mm diameter was pressed 0.1 mm into the cartilage at
an average velocity of 2.5 mm/s and maintained for 300 seconds. Next, after a period of
360 s, the same procedure was repeated in the same location with a plane-ended, non-
porous, cylindrical indenter of 1.5 mm diameter. Then, a spherical, non-porous indenter of
2 mm diameter was pressed 0.06 mm into the cartilage at an average velocity of 2.5 mm/s
and maintained for 300 seconds. The spherical shape was chosen in order to avoid damage
of the tissue at large indentations. The procedure for the spherical indenter was subse-

quently repeated increasing the initial indentation to 0.12, 0.18, 0.24 and 0.30 mm with a

57

waiting period of 360 s between indentations. The indentation force from a transducer
(model JP25; Data Instrument, Acton, MA, U.S.A.) (25 lb [11.3 kg] capacity), attached
between the probe and the stepper, was amplified (model SG71, Validyne, Northridge,
CA, U.S.A.) and collected at 1000 Hz in the first second and at 20 Hz in the remaining
period by a personal computer with an analog-to-digital board. The initial maximum force
for each indentation was reported as the peak force and the force at 300 s was reported as
the equilibrium force. Recovery of the cartilage was monitored after each 360 s waiting
period by indenting manually 0.005 mm and observing the force developed in the indenter.

Finally, tissue thickness was measured at each location with a penetrating needle probe.

4.2.2 Numerical

A numerical comparison of the hyperelastic and the hypo-elastic responses for the
case of uniaxial stress was undertaken. The hyperelastic solution was obtained using the
program Mathematica (Wolfram Research Inc., Champaign, IL) by implementing equa-
tion (4.23) and using the coefficients a0 = 3.073, a1 = 0.8, a2 = -0.390, a3 = 0.245, a4

=0.379, a5=0.175, a6 = 0.0137, a7=-0.379 and n = 0.020. For this set of coefficients the

engineering constants of the solid skeleton were Egg) =1.07 MPa, E281) =5.04 MPa, = 0.187,

VS) =0.0, G13 = 0.19 MPa. As the body deformed, the engineering parameters Egg) and

vgj’ were calculated in Mathematica using equations (4.9) and (4.13) and similar equa-

tions were used to calculate Egsl) , vgsz) and G13. Parameters E11s1) , Egg) , vgsl) , VS) were

updated every 0.1 increment of stretch and G13 was updated every 0.1 increment of shear

deformation.The hypo-elastic solution was obtained with the finite element program

58

MARC (MARC Analysis Research Corporation, Palo Alto, CA). A mesh of 2 rectangular
axisymmetric elements was employed in the analysis. The components of the hypo-elastic
tensor were calculated with equations (4.14)-(4.18) by using the values of the engineering
parameters. The components of the hypo-elastic tensor were input into the subroutine
hypela.f which allows the definition of hypo-elastic material behavior in the program

MARC (1992).

The fitting procedure for each pair (equilibrium and rapid loading) of force-dis-
placement experimental curves was as follows. First, the infinitesimal elastic constants
were estimated at each location using the force-time relaxation curves obtained with the
plane indenters and the thickness of the cartilage following the procedure described in

Chapter 3. Next, an iterative procedure was initiated by guessing the coefficients a1 and n

of the strain-energy function. The other coefficients were obtained by solving equations
(B.3)-(B.7). The uniaxial responses for equilibrium and rapid loading as well as the pure
shear response were obtained using Mathematica according to the hyperelastic formula-
tion described in equation (4.23). These responses were then used to obtain the engineer-

ing parameters at different levels of deformation using equations (4.8) and (4.13) for Eff.)

and vgi’ respectively, and similar equations for E381),v(152) and 613. The uniaxial rapid

loading responses (in directions 1 and 3) were used to obtain the effective engineering
parameters E11, E33 (i.e. parameters of the cartilage under rapid loading). The Poisson’s
ratios during rapid loading were defined by the incompressibility condition as v31 = 0.49

and v12 = l - 0.5(E11/E33). The engineering parameters E11. E33, v12, E381) , 13g? , Vii.) and

59

vgj’ were updated every 0.1 change of stretch, and 013 every 0.1 change of shear defor-

mation. Thus, the components of the hypo-elastic tensor, calculated according to equa-
tions (4.14)-(4.18), were input into the subroutine hypela.f A mesh of 1 mm thickness was
generated with axisymmetric rectangular elements of 0.0625 x 0.0625 mm. Before each
curve fitting, this master mesh was sealed in order to attain the real thickness at the loca-
tion of testing. For the solution of the contact problem the spherical indenter was consid-
ered as a rigid body. For each iteration the program was run two times, one with the hypo-
elastic tensor corresponding to the equilibrium response and the other with the hypo-elas-
tic tensor corresponding to the rapid loading response. The theoretical force-displacement
curves were then saved and compared with those obtained experimentally. Based on this

comparison, changes were made in the coefficients a1 and n and the iteration was repeated

until the least square error was observed to be minimum. The infinitesimal elastic con-
stants and the coefficients of the strain energy function from the four locations in each

patella were averaged.

4.3 Results

A resistance force was developed when indenting manually after each 360 s wait-
ing period, which indicated that the tissue recovered its initial configuration. Differences
were noticed between the shape of equilibrium and rapid loading force-displacement

curves (Figure 4.3).

Good correlation was observed between the hypo-elastic and hyperelastic uniaxial

responses (Figure 4.4). Slight differences may be due to the discrete variations prescribed

for the hypo—elasticity tensor.

The values of the infinitesimal elastic constants of the solid skeleton estimated

from the indentations tests with plane-ended indenters were (five animals, mean (standard
deviation)): 5‘11): 3.37 (0.99)M1>a.r~:‘3§’ = 0.86 (0.14) MPa, G13 = 0.18 (0.03) MPa and

vgj’ =0.13 (0.04). Acceptable fitting were obtained (Figure 4.3) of the nonlinear force-dis-

placement curves with determination coefficients (r2) of 0.932 (0.017) and 0.966 (0.025)
for the equilibrium and rapid loading cases, respectively. The values of the coefficients n
and al, used to adjust the nonlinear behavior, were 0.58 (0.35) and 4.51 (2.14), respec-
tively. The other coefficients of the strain-energy function can be obtained with equations
(B.3)-(B.7). The theoretical uniaxial stress curves in direction 3 obtained with these mate-
rial properties are presented in Figure 4.5, where, for the rapid loading curve, the total
stress was decomposed into the pressure in the fluid and the stress in the solid skeleton.
For this uniaxial stress case the contribution of fluid pressure to total stress was 74.7% for

infinitesimal deformations and 92.0% at a stretch of 0.6.

4.4 Discussion

A method was presented to determine the non-linear material properties of trans-
versely isotropic biphasic cartilage based on a hypo-elastic model. With this model it was
possible to fit simultaneously the high peak force attained during rapid indentation as well
as the equilibrium response. In addition, since the experimental protocol was based on
indentation tests on intact cartilage, the material properties obtained in this study provided

a representation of the in situ mechanical behavior of the tissue. In other test protocols, the

61

mechanical behavior may be modified by disruptions in the tissue when cutting specimens
for testing. The method was implemented using commercial software (MARC, Mathemat-
ica), readily available to other investigators, and it was made consistent with a previously

proposed hyperelastic formulation.

A more pronounced non-linear behavior was observed in the experimental force-
displacement curves obtained for rapid loading compared to those obtained for equilib-
rium. A significant nonlinear behavior was also observed in a theoretical uniaxial curve for
rapid loading (Figure 4.5 (b)), which showed the substantial contribution of fluid pressure
to total stress. This was even more considerable at large displacements. This revealed the
important role of interstitial fluid in supporting compressive stress during rapid loading
and large deformations, which is consistent with theoretical analyses presented by Kelkar
and Ateshian (1995) using infinitesimal assumptions. This effect is also observed in the
testing of cartilage cylindrical plugs under confined compression stress-relaxation (Soltz
and Ateshian, 1997). A significant value of fluid pressure also implies that tensile stress
components might be present in the solid skeleton, as was found at the surface of articular
cartilage in infinitesimal deformation analyses of a transversely isotropic layer indented
by a rigid sphere reported in Chapter 2. These tensile stresses might be responsible for sur-
face fissures observed during impact experiments on the patellofemoral rabbit joint (Haut
et a1, 1995; Newberry et al., 1997). However, more theoretical and experimental studies
need to be undertaken using more realistic contact geometries in order to confirm the pres-

ence of tensile stresses during impact events.

The theoretical uniaxial stress curve (Figure 4.5 (b)) obtained using the average

properties from this study also showed a significant strain-stiffening behavior under rapid

62

loading. This will help limit unrealistically large deformations anticipated using infinitesi-
mal elastic properties in finite element analyses of joints under impact loading. Since the
hypo-elastic model was capable of representing the nonlinear behavior of the tissue under

rapid loading, it can now be used in our analytical studies of joints under blunt impact.

We are not aware of properties reported elsewhere for non-linear transversely iso-
tropic cartilage. The uniaxial stress-stretch curve for the equilibrium response (Figure
5(a)) obtained using the average properties form this study indicated a significant linear
behavior for stretches between 0.8 and 1, which is consistent with experiments by Mow et

al. (1980) in bovine patellofemoral cartilage. A hyperelastic formulation similar to that of

equation (4.24), with I; in the denominator, has been used for isotropic cartilage. Using

this isotropic formulation Ateshian et al. (1997) obtained a value of 11 equal to 0.35 (0.29)
for bovine shoulder and Holmes and Mow (1990) reported values of 11 equal to 0.761 and
1.105, for bovine and human specimens, respectively. These compare with 0.58(0.35) in

the current analysis using rabbit cartilage.

Incompressibility was assumed to model the initial, rapid indentation of the bipha-
sic cartilage. The study undertaken by Eberhardt et al. (1990) using the infinitesimal
unconfined solution obtained by Armstrong et al. (1984) reports differences in force
between the biphasic solution and the elastic solution of only 0.14% at 200 ms after the
application of load. To date in the literature no analysis has been undertaken considering
other geometries and larger displacements. However, since fluid flow is proportional to
pressure gradient, it would appear that the unconfined problem analyzed by Eberhardt

(1990) may represent a conservative estimate of the degree of incompressibility for short

63

duration loading, taking into account the more confined stress field existing in a layer
indented by a rigid sphere. In addition, considering that permeability is reduced at large
strains (Mow et al., 1984) the incompressibility assumption we have made for short times

would be further justified.

A limitation of the model proposed here was that during the curve-fitting proce-
dure a significant time was necessary to complete each iteration (between one and three
hours on a Sun Ultra 1 workstation (Sun Microsystems, Palo Alto, CA)). Part of the time
was used to run a a finite element contact model for the analysis of a transversely isotropic
hypo-elastic layer indented by a rigid sphere, and part was used to determine the compo-
nents of the hypo-elastic tensor using the program Mathematica. Changes or simplifica-
tions in the numerical method would need to be devised in the future in order to efficiently
analyze large amounts of experimental data. This model, however, is a necessary step to
understand the states of mechanical stress in cartilage of diarthrodial joints under impact

loading situations.

 

 

 

 

 

0.2} <-
0.4 - 0.6 0.8 p 1 1:2
Stretch t '
4;: -023
I u- “:2- t
I I I :/ '0.4 r
I z ” ’ / -0.6 ’-
" ,/ -0.8 L
l/ '1 r
./ ’1-2} Stress (MPa)
(a) Equilibrium curves
0.4 0.6 0.8 1.2
Streteh‘ A i A A‘ i "
-10.
F
-201
-30.
1 Stress (MPa)
(b) Rapid loading curves

Figure 4.1 Stress-Stretch curves from the uniaxial solution with a constant coefficient

a1=0.8 and different values for the coefficient 11 (0.020, 0.067 and 0.113)

65

0.4 ‘
Stretch

(14

1 _A_

412

-0.4 ,-

416

-0.8
-1

032:

vvv'vvv

v'1111 vv v'

 

-L2

(a) Equilibrium curves

(16

(18

EStress (MPa)

1.2

 

 

Stretch

  

-100

-120 .

v firv'vrv'vvv'

 

-l40

(b) Rapid loading curves

'-Stress (MPa)

Figure 4.2 Stress-Stretch curves from the uniaxial solution with a constant coefficient n =

0.113 and different values for the coefficient a1 (0.8, 2.0, 3.0)

 

0.7 4 I r I l

Hypo-elastic model

 

 

 

 

0.5-

 

- _ _ Experimental

05............._.......... ..................,................................f................_,

04.1-...“ I ..

2 ’
, . . l . ,
v 03-. ........- “. ............;..~..,.... ..;........... ...L. .. ........;. ..... ..-...i,......-..—(
o ' t T f i i
U 2
5-
O a u /
LL : : : ’ . : :
. L L/ t L l
f /
- 1’
1 I .
D /
011—. ...; J’. .. ~ _
’1
I C
I

 

 

 

00 0.05 0.1 0. 15 0.2 025 0.3 0. 35
Indentation (mm)

(a) Equilibrium curves

 

6 1 1 g t r

Hypo-elastic model ’

 

 

 

 

 

Force (N)

 

 

 

0 0.05 0.1 0.15 02 025 0.3 0.35

Indentation (mm)
(a) Rapid loading curves

Figure 4.3 Experimental force-displacement curves obtained with the spherical indenter of

2 mm diameter showing the typical degree of fit obtained with the theoretical model

67

 

 

 

 

 

 

 

 

 

° ‘ ‘ l T l
Hypo-elastic
...—.... Hyperelastic
4.05- ............................................... .1
-o,1- ~
’6?
m _0161— ......................... —
v:
E -o.2~ x x -
U)
_025_...,.. .. ......3
_o'3L- ............................................................................... _1
_035 I 1 1 l l
0,4 05 06 07 08 0.9
Stretch

(3) Stress versus stretch

 

 

 

 

 

 

 

 

0.08 t t t 1 I
' ' ' __ Hypo-elastic
0.0% , . 44—111 Hyperelastic q
A .
E 0%.. ......................... 1 ............................. r ............................ ...
... .
t: :
g 006).. ....... ......... .. .......: ..................................................... .1
o .
3i ;
8" 0.04» I d
73' W 4
g z
0702-. ..................... .................................................... _
001,. ........................... i ........................................................ .1
0 l I I J l
04 05 06 07 08 1
Stretch

(b) Lateral displacement versus stretch

Figure 4.4 Comparison of Cauchy stress and lateral displacement between the hypo-elastic

and hyperelastic solutions for uniaxial loading

68

‘V‘

0.6
0.4 '-

0.2}
0.6; 2J1? .- : ‘ ‘ 51.2‘ - 31,4‘

Stretch
-0.2

'VVVV‘V‘V

 

'vvv‘vvv

.12 Stress (MPa)

(a) Equilibrium response

0.6 0.7 0.8 0.9

A A A

    

vvvv'vrvv‘rvvvv

_ _ _ — Fluid pressure (MPa)

Total stress (MPa) :'
-6 L Stress (MPa)

 

(b) Rapid loading response

Figure 4.5 Uniaxial stress-stretch curves obtained using the average properties from this
study for rabbit retropatellar cartilage. Observe the substantial contribution of water pres-

sure to total stress in the rapid loading response

69

CHAPTER 5

PLANE STRAIN ANALYSIS OF THE RABBIT RETROPATELLAR
CARTILAGE UNDER INIPACT LOADING

The hypo-elastic material properties reported in Chapter 4 were used to analyze
the shear stress distribution in the rabbit retropatellar cartilage under impact loading. The
analysis was undertaken using the effective properties of the biphasic cartilage. Stress
results from this analysis were then compared with those from the linear elastic analyses
reported in Chapter 2 and a large deformation isotropic model reported elsewhere (Garcia

et al., 1997).

5.1 Numerical method

A plane-strain finite element model (Figure 5.1) was used to determine the shear
stress distribution in the rabbit retropatellar cartilage. This two-dimensional model is
appropriate (Li et al., 1995) considering the geometry of the joint and the type of pressure
distribution observed on the retropatellar surface during impact experiments. The model
was run on the commercial code MARC (MARC Analysis Research Corporation, Palo
Alto, CA). The bone was assumed to be homogeneous with an effective Young’s modulus
of 1.0 GPa and a Poisson’s ration of 0.3 (Brown and Vrahas, 1984). The mesh consisted of
700 four node quadrilateral elements for the main body of the cartilage, 201 triangular ele—

ments for the ends of the cartilage and 631 quadrilateral elements for the bone. A local

70

coordinate system was used in order to define the anisotropy of the cartilage according to
the transversely isotropic model proposed in Chapter 2, i.e. the plane of isotropy is parallel
to the surface of the cartilage. Material properties for the cartilage were defined in the sub-
routine hypela.f listed in Appendix B.5 according to the average effective properties
obtained in Chapter 4. The load was applied on the posterior surface of the patella, and
some nodes on the anterior surface were fixed. Input loads on the patella were obtained
experimentally by Newberry et al. (1997) by recording pressures from Fuji film at the
cross section of interest. The model was subjected to two load cases. The first case, termed
low severity, consisted of a pressure distribution equivalent to 70 N on a slice of unit thick-
ness. This pressure was recorded during low energy impact experiments (0.9 J) in which
no injuries were generated on the surface of the cartilage. The second load case, termed
high severity, consisted of a pressure distribution equivalent to a 200 N load on a slice of
unit thickness. This pressure was recorded during high energy impact experiments (6.0 J)
in which injuries were generated at the surface of the cartilage, inside of the zone of con-

tact on the lateral facet.

5.2 Results

Contours of maximum shear stress on the cartilage (Figure 5.2) showed peak val-
ues at the surface of the cartilage inside of the zone of contact. Local maximum values
were also observed at the interface with the subchondral bone near the border of the lateral
facet and near the surface in the medial facet. An important distortional deformation was
present at the end of the lateral facet. Engineering strains in the direction of loading equal

to 35% and 37% for the low and high severity cases, respectively (Figure 5.2) were

71

obtained in the zone of maximum contact pressure.

5.3 Discussion

The stress distribution obtained in this study, which showed peak shear stresses at
the surface of the cartilage in the zone of maximum contact pressure, is consistent with
surface injuries observed during high energy impact experiments on the joint (Newberry et
al., 1997). The values of maximum shear at the surface of 12.5 MPa and 36.1 MPa for the
low and high severity load cases, respectively, were higher than the 6.2 MPa and 12.9 MPa
values obtained using an isotropic model under finite deformations (Garcia et al., 1997).
These results are consistent with those of the infinitesimal analysis presented in Chapter 2,
which shows that peak shear stresses at the surface of transversely isotropic models are
between two and three times higher than those at the surface of isotropic models. The
shear stresses at the interface with the subchondral bone of 6.6 MPa and 13.2 MPa, for the
low and high severity cases, respectively, were similar to the 6.0 MPa and 12.9 MPa values
obtained using the isotropic model. This is also consistent with the infinitesimal analysis
of Chapter 2 which showed similar values of shear stress at the interface with the bone

between the isotropic and transversely isotropic models.

Analyses of in vitro experiments using a large deformation isotropic model by
Atkinson et al. (1998) indicated that a shear stress level of 5.5 MPa is an appropriate crite-
rion for predicting surface fissures. Had a transversely isotropic model been employed, a
significantly higher shear stress level could be anticipated. Based on the shear stress value
obtained in this study at the surface of the lateral facet for the low severity case (12.5 MPa)

and the values of shear stresses at the interface with the bone (13.2 MPa) and near the sur-

72

face of the medial facet (10.5 MPa) for the high severity case, we could speculate that a
shear stress higher than 13.2 MPa would be associated with injury of the cartilage, since
no fissures have been observed at those locations. However, more analytical and experi-

mental work have to be done in order to better elucidate this issue.

Deformations in the direction of loading were reasonable compared to unrealistic
levels obtained using elastic constants from low strain experiments. This appears to be due
to high pressures generated in the fluid during large deformations and rapid loading which
causes a significant strain-stiffening effect. On the other hand, since maximum contact
pressures were developed near the center of the lateral facet, fluid pressure was transmitted
to the end of the facet and caused important distortional deformations there (Figure 5.2).
The significant deformations of the cartilage near the end of the lateral facet also caused
large distortions in the elements of the mesh, which made the analysis fail when relatively
coarse meshes were used. The problem was alleviated using a finer mesh and triangular

elements at the end of the facet.

The analysis presented here was based on the effective properties of the biphasic
cartilage under rapid loading obtained in Chapter 4. Since the engineering parameters of
the solid skeleton are dependent on deformation, it appears that there is not a simple way
to decompose total stresses into fluid pressure and solid skeleton stresses, as was done in
Chapter 2 for the infinitesimal case. Thus, poroelastic capabilities of finite element codes
would have to be used in order to obtain stress distributions in the solid skeleton using the
hypoelastic model. While such an analysis is beyond the scope of this dissertation, it is an
important next step since, as was shown in the infinitesimal analysis of Chapter 2, tensile

stresses develop in the solid skeleton inside the zone of contact of transversely isotropic

73

models, and these stresses, as well as shear stresses, are believed to be responsible for
injuries observed at the surface of the cartilage. Taking into account the important contri-
bution of fluid pressure to support compressive loads on the cartilage under large deforma-
tions, it could be expected that tensile stresses are also present in our hypo-elastic models.
More experiments and analyses with poroelastic models need to be undertaken in order to

better determine the best predictor of fissuring at the surface of the cartilage.

Anterior

 

Posterior

Figure 5.1 Plane strain model for the effective analysis of the retropatellar cartilage

74

6.6 MPa/
. ‘J’
:1 . ,g‘f'fflx-d -’

_ .- _ - ""J'11/ 25 MPa

.
Medral 4-9 MPa Lateral
(a)

 

Medial 12.8 MP3 Lateral

Figure 5.2 Original and deformed shapes of the cartilage and contours of shear stress for
the (a) low severity and (b) high severity cases

75

CHAPTER 6

CONCLUSIONS

In this dissertation a transversely isotropic biphasic model was presented to
describe the mechanical behavior of articular cartilage under impact loading. Cartilage
was modeled as a biphasic material composed of a transversely isotropic hypo-elastic
solid skeleton and a fluid. This model allowed a Young’s modulus in the plane of the sur-
face (plane of isotropy) higher than that in the direction of loading. This was consistent
with the microstructural organization of the tissue, exhibiting a superficial layer of tightly
woven collagen fibers aligned with the surface. Unlike the isotropic model, the trans-
versely isotropic model was capable of fitting high peak loads observed during the initial
stage of rapid indentation on the tissue. This made it possible to estimate infinitesimal
material properties of articular cartilage consistent with the entire force-time curves
obtained in indentation relaxation experiments. Finite element analyses of contact models
showed significant differences in stress distributions between transversely isotropic and
isotropic material models. Stress distributions in the transversely isotropic model appeared
to be more consistent with the observed surface injuries in our rabbit model for the study
of blunt trauma. A hypo-elastic model for the solid skeleton was developed to handle the
nonlinear behavior of cartilage under large deformations. Using this model, nonlinear
properties were estimated for rabbit retropatellar cartilage. These nonlinear properties

were employed to determine shear stresses in the rabbit retropatellar cartilage under

76

impact loading to the patellofemoral joint. The analyses showed the largest shear stresses
on the surface of cartilage, which was consistent with the area where tissue injury has been
observed in high energy impact experiments on the patellofemoral joint.

In Chapter 2 a parametric stress analysis was undertaken, using a simple geometry
for the cartilage layer indented by a solid sphere, in order to compare the stress distribu-
tions in the solid skeleton of isotropic and transversely isotropic models under impact
load. Employing the assumption of incompressibility, valid for short time duration impact
events (Eberhardt, 1990), effective elastic constants were obtained for a material equiva-
lent to a linear, transversely isotropic biphasic cartilage under impact loading. Elastic
methods were then used to analyze contact models of articular cartilage and the stress in
the tissue was decomposed into fluid pressure and stresses in the solid skeleton. Since
these stresses are proportional to internal forces transmitted by the solid phase, they are
good candidates for the prediction of injuries during impact events. By varying the aspect
ratio (radius of contact / thickness of cartilage) in this geometry, conditions similar to
those obtained with more realistic models of diarthrodial joints could be predicted. Results
of the analyses showed significant differences in the stress distribution between the isotro-
pic and transversely isotropic models. In the isotrOpic model the maximum shear stress
was generally located at the interface with the subchondral bone, and in the anisotropic
models it was largest on the surface of the tissue. Hoop stress was always compressive at
the center of contact for the isotropic model, and tensile on the whole surface of contact
for the transversely isotropic models, except for high aspect ratios, in which case the hoop
stress was compressive at the center of contact. The stress distributions for the anisotropic

models were consistent with surface lesions observed in impact experiments on the rabbit

77

patellofemoral joint, without lesions at the interface with the subchondral bone. In addi-
tion, the presence of tensile hoop stresses at the surface of the solid skeleton may explain
the radial fissures documented by Atkinson et al. (1997), and the surface fissures observed
by Silyn-Roberts and Broom (1990) for impact experiments on articular cartilage.

In Chapter 3 a method was developed to determine the infinitesimal elastic con-
stants and permeabilities for the transversely isotropic biphasic cartilage. The method was
intended to offer material properties more consistent with high peak forces observed dur-
ing the early response of a small displacement indentation experiment. These high forces
could be fit with an isotropic model only if the Poisson’s ratio was assumed to be negative.
which is not a realistic assumption for articular cartilage. The method, implemented in the
commercial program Matlab, consisted of fitting two force-time relaxation curves
obtained with two plane-ended cylindrical indenters. Four elastic constants of the solid
skeleton were determined and the other elastic constant, the in-plane Poisson’s ratio was
assumed to be zero. The method allowed the prediction of the two permeabilities. Results
of the analysis, applied to tests on the rabbit retropatellar cartilage, indicated that the
Young’s modulus and the permeability coefficient in the plane of the surface were higher
than those in the direction of loading.

In Chapter 4 a transversely isotropic, hypo-elastic model for the solid skeleton was
adopted to consider the nonlinear behavior of cartilage. Nonlinear properties were deter-
mined by fitting two force-displacement curves (in rapid loading and at equilibrium)
obtained from large deformation, indentation relaxation tests on cartilage using a non-
porous spherical indenter. The method was implemented using commercial software,

readily available to other investigators and it was made consistent with a previously pro-

78

 

posed hyperelastic formulation. The model was able to fit simultaneously the high peak
force attained in a large indentation, as well as the equilibrium response. A more pro-
nounced nonlinear behavior was observed in the experimental force-displacement curves
obtained for rapid loading than for those obtained at equilibrium. A significant nonlinear
behavior was also observed in a theoretical uniaxial curve for rapid loading, indicating the
substantial contribution of fluid pressure to total stress. This effect was even more pro-
nounced at large displacements. This revealed the important role of interstitial fluid in sup-
porting compressive stress during rapid loading and large deformations, which is
consistent with theoretical analyses presented by Kelkar and Ateshian (1995) using infini-
tesimal assumptions and with experiments on cartilage cylindrical plugs under confined
compression stress-relaxation (Soltz and Ateshian, 1997). A significant fluid pressure
under impacting loads also implies that tensile stress components might be present in the
solid skeleton, as found at the surface in the infinitesimal analyses described in Chapter 2.
An important stress-stiffening effect was also observed for rapid loading. This lead to real-
istic deformations in finite element analyses of joints under blunt impact.

In Chapter 5 a finite element analysis is presented for the rabbit patella under
impact loading. This analysis revealed peak shear stress at the surface of the cartilage
inside the zone of contact. This is consistent with surface injuries observed during high
energy impact experiments on the patellofemoral joint. Qualitative results of this nonlinear
analyses were similar to those found in the infinitesimal study presented in Chapter 2,
showing shear stress values at the surface of the cartilage between two and three times
higher than those calculated using isotropic models.

In summary, a transversely isotropic hypo-elastic model appears to offer a more

79

realistic representation of articular cartilage under rapid loading. The material properties
obtained in this study can now be used in poroelastic finite element analyses of the rabbit
patellofemoral joint under impact load in order to determine the stress distributions in the
solid skeleton of the tissue. Analysis of this distribution is crucial to understand the state
of tissue stress under impact loading to joints and the role of interstitial fluid pressuriza-
tion and stresses in the solid skeleton to injury mechanisms. While in this dissertation fluid
pressure and stresses in the solid skeleton have been obtained for the infinitesimal case,
the decomposition has yet to be done for the nonlinear case. Furthermore, since the ulti-
mate goal of this research is to propose strategies to protect the occurrence of osteoarthro-
sis in humans, the analysis proposed here should also be applied to human cartilage
properties. This is a necessary step in order to undertake analytical studies of human joints
for the optimization of padding systems for their protection from injuries during blunt

insults.

80

APPENDIX A

DESCRIPTION OF THE PROGRAM TIPRORM USED TO DETERMINE
INFINITESIMAL TRANSVERSELY ISOTROPIC PROPERTIES OF
CARTILAGE

A.1 Instructions for the user of the program tiprop.m

A.1.1 Data provided by the user

In order to use the program the user must provide the following data:

1. Two files containing the force-time curves obtained in two relaxation indenta-
tion experiments undertaken with plane-ended solid cylinders of 1.0 mm and 1.5 mm
diameter. In these files, the time in seconds must appear in the first column and the force in
Newtons in the second column. For the program to be applicable, the indentation has to be
applied rapidly, in a period less than 200 ms. If the values of force contained in the files
described above are not in Newtons, the user must provide the factor necessary to convert

those forces into Newtons.

2. The thickness at the location of testing. This thickness must be between 0.5 mm

and 1.0 mm.

A.1.2 Programs and data necessary to run the program

In order to run the program, the following information must be available.

1. A version of the program Matlab (The MATH WORKS Inc, N atick, MA). The

81

current procedure was developed using the version 5.1.0.421 of Matlab running in a

ULTRA 1 Sun workstation (Sun microsystems, Palo Alto, CA).

2. Copy of the files tiprop. m, elasvlm, nfun52.m nfun102. m, shapelm, mamemim

and anfun2.m listed in A2.

3. Copy of the data files containing the equilibrium and peak forces obtained for

different combinations of elastic constants. These files are described in A3.
4. Copy of the six files containing the parameters a, or and c, which give the cali-

bration of the master curve of equation (3.3). These files are described in A4.

A.1.3 Steps to run the program

The following steps must be followed to run the program.

1. The file containing the force—time curve for the 1.0 mm diameter indenter must
be copied into another file called dv. The file for the 1.5 mm diameter indenter must be

copied into a file called dx.

2. If the forces are not in Newtons the scale factor must be changed in the line:
fur=d(:,2)*2.452; of the program tiprop.m, where 2.452 is the scale used in our laboratory,
which relates Volts and Newtons. In the case that the force is not in the second column, but

in the third column, this instruction has also to be modified to be: fur=d(:,3)*2.452.
3. The program matlab must be started by typing matlab.
4. After the program matlab is started the word tiprop must be typed.
5. The program then asks for the thickness of the cartilage, which must be input by

the users.

82

6. The program runs until material properties are listed.

A.2 Description and list of the program tipr0p.m

A description and a complete list is presented of all the functions contained in the

program tiprop.m. Figure A] represents the relations among all the functions composing

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

the program.

r shape] .m

nfun52.m
I mamemi.m

elasv1.rn

l 1 shape1.m

tiprop.m nfun102.m
I mamemi.m

 

 

 

 

 

 

 

anfun2.m

 

 

 

Figure A.l Functions of the program tiprom

tiprop.m. This is the main function. First, it loads two force-time indentation
curves for the 1.0 and 1.5 mm indenters and determines the associated peak and equilib-
rium forces. Then, for the 1.0 mm curve it calls the function elasvl.m which determines

three combinations of elastic constants that fit the peak and equilibrium force for the 1.0

83

mm indenter. Thus, it calls the function anfun2.m which determines the permeabilities
associated with each combination of elastic constants for the 1.0 mm and 1.5 mm curves.
Finally, it calculates the anisotropic permeabilities, according to equations (3.1) and (3.2)
and prints the results.

%%%%%%%%%%%%%%%%‘71%%%‘7c%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

% Function tiprop.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%.u%%%%%%%%%
% Array loaded in dv (1 mm) and dx (1.5 mm), time in the first column

% and force (V) in the third column
%

load dv;

d = dv;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop on the indentation curves, kk = l, 1.0 mm indenter
% kk = 2, 1.5 mm indenter
%
forkk= 1:2;
%
time =d(:,l);
nn = length(time):
rm =nn - 100;

3959

The force is assumed to be contained in the third column of the arrays

% in case that it is contained in the second colurrm, the number 3 in the following
% instruction have to be change to 2

% The number 2.452 represents the calibration between Newtons and Volts

%

fur = d(:,3)*2.452;

3933

Determine and cuts before the maximum value

M] = maX(fur);
time = time(i:nn);
time = time - time( 1);
%
% Makes maximum value zero before filtering
%
forc = fur(i:nn) - y;
%
%
%
[13.21] = butter(l.0.07);
forf = filter(b,a,forc);
%
% Readjust maximum value after filtering
%
forc = forc + y;

84

forf = forf + y;
ti = time(]:lOO);
ti = forc(1:100);
ffr = forf( 1:100);
fma = ffi(l);
%
nk = size(time)'.
ni = nk - 300;
ti = time(nimk);
ffi = forf(ni:nk);
%
% Readjust from 1 s for the 1.5 mm indenter
%
if kk = 2,
aa = 0; k = 1;
while aa = 0.
if time(k) > 1.0.
tls = time(k-l);
aa = 1;
n1 = k - 2;
end
k = k + 1;
end
time = time(nlmk);
forc = forc(nl:nk);
end
area(kk) = trapdtimeforc);
%
%
nk = size(ti)-100;
ni = nk - 20;
ave = ffi(ni:nk);
fmi = mean(ave);
fprintf(l,’ Peak force ............. %6.2f \n’,fma);
fprintf(1.’ Equilibrium force ...... %6.2f \n’.fmi);
%
% Calls the function elasvl.m which determines elastic constants that satisfy the peak and
% equilibrium force for the 1.0 mm indenter
%
if kk = 1,
[cl l,e33,g13,n31] = elasv 1(fma,fmi);
end.
%
% Calls the function anfun2.m which determines isotropic permeabilities for each curve
%
[per,nrm] = anfun2(fmi.fma,time,forc.kk);
%
perme(kk,l:3) = per;
nrma(kk,l:3) = nrrn;
load dx;
(1 = dx;
end.
%
for i=l:3;

85

%
%
%

%
%
%

sum(i) = nrma(l,i) + nrma(2,i);
end.
1le = min(sum);
nrma

Estimates anisotropic permeabilities according to equations (3.1) and (3.2)

kil = perme(1,i)
kilS = perme(2,i)
if kilS > kil,
a1 = 0.5;
be = 0.99;
else
a1 = 0.5;
be = 0.99;
end
kl = (be - al)/(-(1-be)/ki1 + (l-al)/ki15)
k3 =(l - al/be)/(l/kil - (al/be)/ki15)
ifkl <= 0.0 |k3 <= 0.0.
sum = nrma(l,:);
1le = min(sum);
nrma(l,i) = nrma( l,i)/area(1);
nrma(2,i) = nrma(2,i)/area(2);

Prints results
fprintf(l”*************$*ll!'IUI‘********$**********$*$**#**************\n\n’);
fprintf(l,’ Ell ................................ %7.2f MPa\n’ ,el l(i));

fprintf(l,’ E33 ................................ %7.2f MPa\n’,e33(i));

fprintf(l.’ G13 ................................ %7.2f MPa\n’.gl3(i));

fprintf(l,’ nu31 ............................... %7.3f\n’.1131(i));

fprintf(l,’ Isotropic permeability 1 mm ind ...%7.2f m4/(N.s)XE-15 \n’,kil);

fprintf(l,’ Estimated norm .................... %7.3f\n’.nrma( 1.i));

fprintf(l,’ Isotmpic permeability 1.5 mm ind .%7.2f m4/(N.s)XE-15 \n’.ki15);
fprintf(l,’ Estimated norm .................... %7.3f\n’.nrma(2,i));

fprintf(l” *Ill’lllkilllllilullIIUIHIHIUIHIUII***************IMMFIIIIF’F'IHIUII'NF****************\n’)

fprintf(l,’ Warning : Anisotropic permeabilities do not appear to fit both experimental curves accurately

\n ’);

%

%

else

nrma(l,i) = mma(l,i)/area(1);
nrma(2,i) = nrma(2,i)/area(2);

fprintf(l” ********I"!!!***Ill************************$******ilflkfinlfllfllnk****\n\n’)

.
9

fprintf(l,’ Ell ................................ %7.2f MPa\n’,el l(i));

fprintf(l,’ E33 ................................ %7.2f MPa\n',e33(i));

fprintf(l,’ G13 ................................ %7.2f MPa\n’,gl3(i));

fprintf(l.’ nu31 ............................... %7.31\n’.n31(i)):

fprintf( l,’ Permeability direction 1 .......... %7.2f m4/(N.s)XE-15 \n’,kl);
fprintf(l,‘ Estimated norm .................... %7.3t\n’,nrma( 1,i));

fprintf( 1,’ Permeability direction 3 .......... %7.2f m4/(N.s)XE-15 \n’,k3);
fprintf( 1.’ Estimated norm .................... %7.3t\n’.nrma(2,i));

fprintf(l,’***********************mamarntumnarumwmanarwrunrwmmwwrwmrm’).
9

86

%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elasvl.m. This function reads the thickness of the cartilage and then calls the func-
tions nfun52.m and nfun102.m, which search the values of peak and equilibrium forces in
the tables for the 0.5 mm and 1.0 mm thickness cases, respectively. Finally, it interpolates
linearly between the elastic constant found by nfirn52.m and nfun102.m. according to the

actual thickness of the cartilage at the location of testing.

%%%%%%%%%%%%%%%%%%%%%%%%%%‘7c%%’7r%’71~%’7o%%%%%’71 %‘71'7r%%%%%
%
% Function elasvl.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%‘7c-%%%%%%%%%%%%%%%
function[el,e3,g3,n3] = elasvl(fple,fe 1e)
%
% fple : peak force from 1 mm indentation curve
% fele :equilibrium force from 1 mm indentation curve
% th : thickness of the cartilage
%
th = input(’ Enter thickness of the cartilage :’);
fa = (th - 0.5)/0.5;
%
if th <= 0.5.
[el,e3,g3,nB] = nfrm52(fple.fele);
fprintf(l,’ Warning: Calculations will be made assuming thickness = 0.5 m \n’);
elseif th >= 1.0.
[el,e3,g3,n3] = nfunlO2(fp1e,fele);
fprintf(l,’ Warning: Calculations will be made assuming thickness = 1.0 m \n’);
else
[xl,y1,wl,zl] = nfun52(fple,fele);
[x2,y2,w2,22] = nfun102(fple,fele):
el = x1 + (x2 - xl)"‘fa;
e3 = yl + (y2 - y1)*fa;
g3 = wl + (w2 - wl)*fa:
n3 = z1+(22 - zl)*fa;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nfun52.m. This function performs a loop on the tables containing the data repre-

sented in plots like Figure 3.3 for the 0.5 mm thickness case. If a point (feq, fpe) is con-

tained in the table it interpolates linearly in order to determine the set (Egg), G13)

associated to the point. For the interpolation procedure the area is divided into triangles.

87

The function shapelm is called in order to generate the shape functions for a three-noded

triangle. After it searches in all the tables for the 0.5 mm thickness case, it then calls the

subroutine mamemi.m which determines three sets of elastic constants with associated

high, medium and low values of E‘lsl).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function nfun52.m

))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[zell,ze33,zgl3,zn31] = nfun52(fple,fele);

%
%
% fpt and fet are the experimental peak and equilibrium forces
%
%
fpt = prC;
fet = fele;
%
% Initialize residual
%
ka = 1;
rmin = 0;
%

% Loop on tables for t = 0.5 mm. ifil from 0 to 12 indicate the values ofell
% associated with the table, and ipoi indicates the value of n31 associated
% with the table. E.g. ifil=3, ipoi=4 indicates ell = 3MPa and n31 = (4-1)*0.05=
% 0.15
%
km = 0;
for ifll = 0:12.
ki = 0;
for ipoi = 1:7,
%
31 = ’r’;
ipk = (ipoi - 1)*5;
nu31 = ipk/100.0;
ell = ifil;
if ifrl == 0.
ell = 0.5;
end;
32 = num23tr(ifll);
s3 = num23tr(ipk);
s3 = strcat(’n’,s3,’.dat’);
fno = strcat(sl,s2,s3);
%
fid = fopen(fi10.'r');

88

nc = fscanf(frd,’%5d’,l); nr = fscanf(fid,'%5d,l’);
c = fscanf(fid,’%5f’,[nr.ncl);
e = fscanf(fid,’%5f’,[nr.nc]);
%
e33 = fscanf(fid,’%5f‘,[l,nc]);
g13 = fscanf(fid,’%5f’,[l,nr]);
fclose(fid);
err = 9999;
%
%
%
%
xp = fet;
w=¢c
fl = 0;
for it = l:(nr - l),
for ic = l:(nc - 1).

%
% Localize coordinates of element
%
x1 =c(ir ,ic );
x2 = e(ir ,ic+1);
x3 = c(ir+l,ic+1);
x4 = c(ir+l.ic );
%
w=akrcx
y2 = e(ir ,ic+l);
y3 = e(ir+l,ic+1);
y4 = e(ir+l,ic );
%
%
if ((xp >= x1) & (xp <= x3)) & ((yp >= yl) & (yp <= y3))
%
x5 = x4;
x6 = x2;
y5 = yl + (y2 - yl)*(x5 - x1)/(x2 - x1);
y6 = y4 + (y3 - y4)*(x6 - x4)/(x3 - x4);
%
if (xp >= x1) & (xp <= x5)
ylow = yl + (yS - yl)*(xp - xl)/(x5 - x1);
yhig = yl + (y4 - yl)*(xp - x1)/(x4 - x1);
if (yp >= ylow) & (yp <= yhig)
fl=1;row=ir;col=ic;
%
xe(1)= x1; ye(l) = yl;
xe(2) = x5; ye(2) = y5;
xe(3) = x4; ye(3) = y4;
%
d1 = sqrt((x5 - x1)"2 + (y5 - yl')"2);
d = sqrt((x2 - x1)"2 + (y2 - yl)"2);
%

un( 1) = e33(col); un(3) = e33(col);
un(2) = un(l) + (e33(col+l) - un(1))*dl/d;
e331 = shapel(xe,ye,un,xp,yp);

89

un(l) = g13(row); un(2) = g13(row);
un(3) = gl3(row+l);
gl3i = shapel(xe.ye.un.xp.yp);

%
%

%

end

elseif (xp > x5) & (xp <= x2)
ylow = y1+(y2 - y1)*(xp - xl)/(x2 - x1);
ymed = Y4 + (y2 - y4)*(xp - x4)/(x2 - x4);
yhig = y4 + (y3 - y4)*(xp - x4)/(x3 - x4);
if (YP >= ylow) & (yp <= ymed)

%

%

%

%
%

%

%

%

%
%

fl=2;row=ir;col=ic;
xe(1)= x5; ye(l) = y5;
xe(2) = x2; ye(2) = y2;
xe(3) = x4; ye(3) = y4;

d1 = sqrt((x5 - xl)"2 + (y5 - yl)"2);
d = sqrt((x2 - x1 )Az + (y2 - yl)"2);

un(2) = e33(col+l); un(3) = e33(col);
un(1)= e33(col) + (e33(col+l) - e33(col))*dl/d;
e33i = shapel(xe,ye,un.xp.yp);

un(l) = g13(row); un(2) = g13(row);
un(3) = gl3(row+l);
g13i = shape1(xe.ye.un.xp.yp);

elseif (yp > ymed) & (yp <= yhig)

fl=3;row=ir;col= ic;
xe(l) = x2; ye(l) = y2;
xe(2) = x6; ye(2) = y6;
xe(3) = x4; ye(3) = y4;

d1 = sqrt((x6 - x4)"2 + (y6 - y4)"2);
‘1 = sqrt((x3 - x4)"2 + (y3 - y4)"2);

un(l) = e33(col+l); un(3) = e33(col);
tm(2) = e33(col) + (e33(col+l) - e33(col))*d1/d;
e33i = shape1(xe,ye,un.xp.yp);

un(1)= g13(row); un(2) = gl3(row+l);
un(3) = gl3(row+l);
g13i = shapel(xe,ye.tm.xP.yp);

end

elseif (xp > x2) & (xp <= x3)
ylow = y2 + (y3 - y2)*(xp - x2)/(x3 - x2);
yhig = y4 + (y3 — y4)*(xp - x4’)/(x3 - x4);
if (yp >= ylow) & (yp <= yhig)

90

fl=4;row=ir;col= ic;
%
xe(l) = x2; ye(1)= y2;
xe(2) = x3; ye(2) = y3;
xe(3) = x6; ye(3) = y6;
%
d1 = sqrt((x6 - x4)"2 + (y6 - y4)"2);
d = sqrt((x4 - x3)"2 + (y4 - y3)"2);
%
tm(l) = e33(col+l); un(2) = e33(col+l);
un(3) = e33(col) + (e33(col+1) - e33(col))*dl/d;
e33i = shapel(xe,ye,un.xp,yp);
%
un(l) = g13(row); un(2) = gl3(row+l);
un(3) = gl3(row+l);
g13i = shapel(xe,ye,rm.xp.yp);
%
end
end
end
end
end
%
if fl ~= 0.
ifki == ,
km = km +1;
ie11(km) = ell;
end

h=fi+h
in31(km,ki) = ipoi;
ie33(km,ki) = e33i;
igl3(km,ki) = g13i;
r = ell;
if r > rmin,
rmin = r;
e3 = e33i;
g3 = g13i;
el = ell;
n3 = nu31;
end
%%%%%
end
%
clear c;e;e33;gl3;
%
% End of the loop on the files
%
end
end
iell
in31
ie33
igl3

91

[zell,ze33,zg13,zn31] = mamemi(iel l,ie33.ig13,in31)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nfim102.m. This function performs the same operations of nfun52.m for the 1.0
mm thickness tables.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function nfun102.m
O'>(E')(E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
frmction[zel 1,ze33,zgl3.zn31] = nfun102(fple,fele);
%
%
fpt = fple;
fet = fele;
%
% Initializes residual
%
ka = l;
rmin = 0;
%
km = 0;
for ifil = 1:20,
k1 = 0;
for ipoi = 1:7.
%
$1 = ’ lr’;
ipk = (ipoi - 1)"‘5;
nu31 = ipk/100.0;
ell = ifil;
52 = num25tr(ifrl);
s3 = num251r(ipk);
s3 =strcat(’n’,s3.’.dat");
fno = strcat(sl,sZ,s3):
%
fid = fopen(fno,’r’);
nc = fscanf(fid,’%5d’,l); nr = fscanf(fld,’%5d,l’);
c = fscanf(fid,’%5f’,[nr,nc]);
e = fscanf(fid,'%5f’.[nr.nCl);
%
e33 = fscanf(fid,’%5f’ ,[l,nc]);
gl3 = fscanf(fid,’ %5f’,[1,nr]);
fclose(fid);
err = 0;
%
%
%
%
xp = fet;
yp = fpt:

92

%
%
%

%

%
%

%

%

%

%

%

%

%
%

%

fl=0;
forir= l:(nr- l).
for ic = l:(nc - l).

Localize coordinates of element

x1 =c(ir .ic );
x2 = e(ir ,ic+1);
x3 = c(ir+l,ic+l);
x4 = e(ir-+1.10 );

yl = e(ir .ic );
y2 = 601' .iC+1);
y3 -.- e(ir+l,ic+l);
y4 = e(ir+l,ic );

if ((xp >= x1) & (xp <= x3)) & ((yp >= yl) & (yp <= y3))

x5 = x4;
x6 = x2;
y5 = yl + (y2 - yl)*(x5 - xl)/(x2 - x1);
y6 = y4 + (y3 - y4)*(x6 - x4)/(x3 - x4);

if(Xp>=Xl) & (xp <= x5)
ylow = yl + (y5 - yl)*(xp - xl)/(x5 - x1);
yhig = yl + (y4 - yl)*(xp - x1)/(x4 - x1);
if (yp >= yIOW) & (yp <= yhig)
fl: I;I‘OW=ir; COI=IC;

xe(1)= x1; ye(l) = yl;
xe(2) = x5; ye(2) = y5;
xe(3) = x4; ye(3) = y4;

d1 = sqrt((x5 - xl)"2 + (y5 - yl)"2):
d = sqrt((x2 - xl)"2 + (y2 - yl)"2);

un(l) = e33(col); un(3) = e33(col);
un(2) = un(l) + (e33(col+l) - un(1))*dl/d;
e331: shapel(xe,ye,un,xp,yp);

un(1)= g13(row); un(2) = g13(row);
1m(3) = gl3(row+l);
g13i = shapel(xe,ye,un,xp.yp);

end

elseif (xp > x5) & (xp <= x2)
ylow = y1+(y2 - yl)"‘(xp - xl)/(x2 - x1);
ymed = y4 + (y2 - y4)"‘(xp - x4)/(x2 - x4);
yhig = 1’4 + (Y3 - y4)*(xp - x4)/(x3 - x4);
if (yp >= ylow) & (yp <= ymed)
fl = 2; row = in C01 = ic;

93

%

%

%

%
%

%

%

%

%
%

xe(l) = x5; ye(1)= y5;
xe(2) = x2; ye(2) = y2;
xe(3) = x4; ye(3) = y4;

d1 = sqrt((x5 - x1)"2 + (y5 - yl)"2);
d = sqrt((x2 - xl)"2 + (y2 - yl)"2);

un(2) = e33(col+1); un(3) = e33(col);
un(l) = e33(col) + (e33(col+l) - e33(col))*d1/d;
e33i = shapel(xe.ye,un.xp.yp);

un(l) = g13(row); un(2) = g13(row);
un(3) = g13(row+1);
g13i = shapel(xe,ye.un.xp,yp);

elseif (yp > ymed) & (yp <= yhig)

fl=3;row=ir;col=ic;
xe(1)= x2; ye(1)= y2;
xe(2) = x6; ye(2) = y6;
xe(3) = x4; ye(3) = y4;

d1 = sqrt((x6 - x4)"2 + (y6 - y4)"2);
d = sqrt((x3 - x4)"2 + (y3 - y4)"2);

un(1)= e33(col+1); un(3) = e33(col);
un(2) = e33(col) + (e33(col+1) - e33(col))*dl/d;
e33i = shapel(xe.ye.un.xp.yp);

un(1)= g13(row); un(2) = gl3(row+l);
un(3) = gl3(row+l);
g13i = shapel(xe.ye.rm.xp.yp);

end

elseif (xp > x2) & (xp <= x3)
ylow = y2 + (y3 - y2)*(xp - x2)/(x3 - x2);
yhig = y4 + (y3 - y4)*(xp - x4)/(x3 - x4);
if (yp >= ylow) & (yp <= yhig)

%

%

%

%

fl=4;row=ir;col=ic;
xe(l) = x2; ye(1)= y2;
xe(2) = x3; ye(2) = y3;
xe(3) = x6; ye(3) = y6;

d1 = sqrt((x6 - x4)"2 + (y6 - y4)"2);
d = sqrt((x4 - x3)"2 + (y4 - y3)"2);

un(1)= e33(col+1); un(2) = e33(col+1);
un(3) = e33(col) + (e33(col+l) - e33(col))*dl/d;
e33i = shape1(xe,ye.un,xp.yp):

un(l) = g13(row); un(2) = gl3(row+l);

94

un(3) = gl3(row+l);
g13i = shapel(xe,ye,un,xp,yp):
%
end
end
end
end
end
%
if fl ~= 0.
if ki = 0.
km = km + 1;
iel l(km) = ifil;
end,
h=h+h
in31(kIII.ki) = ipoi;
ie33(km,ki) = e33i;
ig13(km.ki) = g13i;
r = e11;
if r > rmin,
rmin = r;
e3 = e33i;
g3 = g13i;
el = ell;
n3 = nu31;
end
%%%%%
end
%
clear c;e;e33;gl3;
%
% End of the loop on the files
%
end
end
iell
in31
ie33
igl3
[zell,ze33,zgl3,zn31] = mamemi(ie11,ie33,igl3,in3l)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

shapelan. This function generates the three shape functions of a three node trian-
gle.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%‘71%%%%%%%%%%%%%%%
%

% Function shapel.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

95

function [u] = shapel(xn,yn,un,x,y)
%
a = zeros(1.3);
b = zeros(l.3);
c = zeros( 1,3);
%
a(l.l) = xn(2)"'yn(3) - xn(3)*yn(2);
a( 1.2) = xn(3)*yn(l) - xn( l)*yn(3);
2103) = xn(l)*yn(2) - xn(2)*yn(l);
%
b(l.l) = yn(2) - yn(3);
1102) = yn(3) - yn(l);
b(l,3) = yn11) - yn12);
%
C(l.1)= m(3) - xn(2);
C(12) = xn(l) - mm:
00.3) = xn(2) - xn(1);
%
del = (210.1) + a(l.2) + a(l.3));
%
ni =(a(l,1)+ b(l,l)"'x + c(1,1)*y)/del;
nj = (2111.2) + b(l,2)*x + c(l.2)*y)/del;
nk = (a(l.3) + b(l,3)*x + c(1,3)*y)/del;
%
u = ni*tm(l,l) + nj*un(l,2) + nk*rrn(l.3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mamemi.m. This function determines the three sets of elastic constants with asso-

ciated high, medium and low values of E381) from all the combinations of elastic constants

matching the peak and equilibrium force for the 1 mm indenter.

%%%%%%%%%%%%% "c r "r1 'n 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function mamemi.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionlzel l,ze33,zg13,zn3l] = mamemi(iel l,ie33,igl3,in31);
%
%
a = size(in31);
at = 3(1);
nc = a(2);
nth = round((nr+l)/2);
%
iin(l) = l;
iin(2) = nrh;
iin(3) = m';
%
for i=l:3;
ii = iin(i);
nz = length(find(in31(ii,:)));

96

nz = round(nz/2);

zell(i) = iel l(ii);

ze33(i) = ie33(ii,nz);

zgl3(i) = ig13(ii.nz):

2.1131(i) = in31(ii.nz):

m31(i) = 0.05*(zn31(i) - 1):

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

anfim2.m. This function fits the experimental curves with the master curve defined

in equation (3.3) and determines the isotropic permeabilities associated with each curve.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function anfun2.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[perme,nrma] = anfun2(fmi,fma,time,forc,kk);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%‘71.%%%%%
%
pera = [0.5 2.0 5.0 10.0 18.0 25.0];
%%%%%%%%%%%%%%%%%%%‘71‘70‘76%%%%%%‘7o%%%%%
delf = fma - fmi;
%
% LooponthevaluesofEll,ki= l, min
% ki = 2, med
% ki = 3, max
%
s1 = num2str(kk);
s2 = strcat(sl,’.dat’);
%
for ki = 1:3;
%
if k1 = 1.
s1 = 'min’;
elseif ki = 2,
81 = ’med’;
else
31 = ’max';
end
%
% Open files with data for the coefficients of master function
%
fno = strcat(sl,s2);
frd = fopen(fno,’r’);
nper = fscanf(fid,’%5d’,l);
nran = fscanf(fid,’%5d',1);
dan = fscanf(frd,’%5f’.nran);
aaa = fscanf(frd,’%5f’,nran);
alfa = fscanf(fld,’%5f’.nran);
c3an = fscanf(fid,’%7f’,[nper,nran]);
c3an = c3an’;
fclose(fid);

97

%
global XD KD;
if kk = l,
ifdelf>= dan(l) & delf <= dan(2),
KD = l;
elseif delf >= dan(2) & delf <= dan(3),
KD = 2;
elseif delf >= dan(3) & delf <= dan(4).
KD = 3;
elseif delf > dan(4)
fprintf(l,’ Delf is not contained in the array \n');
KD = 3;
delf = dan(4);
else
fprintf(l,’ Delf is not contained in the array \n’);
KD = 1;
delf = dan(l);
end
if KD ~= 0,
nd = KD;
XD = (delf - dan(nd))/(dan(nd+1) - dan(nd));
end
end
%%%%%%%%%%%%%%%%%%%%%%%
if k = 2.
nd = KD;
delf = dan(nd) + XD*(dan(nd+l) — dan(nd));
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Loop on the
delforce ranges above and below delf
%
for ii=1 :2
a = aaa(nd);
a1 = alfa(nd);
%
aa = 0;
% c3 = 0.01;
c3 = c3an(nd.1);
ni = 9999.9;
del = 0.005;
pp = 1:
%
while aa == 0,
%
%
% Builds analytical curves into the array fan (f analytical), according to master function
%
fan = fmi + delP"((1-al)./cosh(c3*time) + al./( 1 + a*sinh(c3"‘time)));
%
% Plots both curves
%
% If a graphic representation of the fitting procedure is desired, this plot must be
% activated

98

% figure( 1)

% plot(time,fan,time.forc,’r’)
% grid

%

nn = length(fan):
f3 = fan - forc;
t2 = time;
13 = abs(f3);
norm = trapz(t2.f3);
%
%
if norm > ni & pp == 1.
c3 = c3 - 0.01;
del = 0.0002;
pp = 2:
norm = trim;
elseifnorm>ni&pp==2.
a = 1;
nr(ii) = norm;
end
nim = ni;
ni = norm;
c3 = c3 4» del;
%
end
% fprintf(l,’ Value of c3 which better fits the exp. curve : %6.2f \n‘,c3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Determination of the permeability, based on delf and c3
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Now calculates peme from the value marked by nk
%
if c3 < c3an(nd,l),
fprintf(’ c3 is not contained in the range R \n‘);
per(ii) = pera(l);
aa = 1;
fprintf(l.’ Wanting : Permeability can be overestimated \n’);
elseif c3 > c3an(nd,6),
fprintf(’ c3 is not contained in the range k \n’);
per(ii) = pera(6);

aa = 1;
fprintf(l,’ Wanting : Permeability can be underestimated \n’);
else
aa = 0;
k = 2;
while aa = 0.
if c3 < c3an(nd,k),

c3bef = c3an(nd.k- 1);

c3aft = c3an(nd.k);

per(ii) = pera(k-l) + (pera(k) - pera(k-l))*(c3 - c3bef)/(c3aft - c3bef)
aa = 1;

99

else
k = k + 1;
end
end
end
nd = nd + 1;
end % End of for ii
%
% Now interpolates between permel and perrne2
%
nd = nd - 1;
perme(ki) = per( 1) + (per(2) - per(1))*(delf - dan(nd-l))/(dan(nd) - dan(nd-1));
nrma(ki) = nr(l) + (nr(2) - nr( 1) )"'(delf - dan(nd-l))/(dan(nd) - dan(nd-l));
%
% Ends loop on max, med and min
%
end,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A.3 Names and organization of resident data files

Tables representing plots like Figure 3.3 were generated for values of E‘fl’ from 1

MPa through 12 MPa for the 0.5 mm thickness case and from 1 MPa through 20 MPa for

the 1.0 mm thickness case at increments of 1.0 MPa, and values of vgj’ from 0.0 through
0.3 at increments of 0.05. The name of each table for the 0.5 mm thickness case was of the

form rxxnyy.dat, where xxx was the value of Ea? associated to the table and yy was equal

to 100 times the value of vgj’ associated to the table. For example, the table r12n5.dat was

created using values of E‘f,’ = 12 MPa and vgj’ = 0.05. Names of the tables for the 1.0 mm

thickness case followed the same convention, but the number 1 was added in the first posi-

tion. For example, the table 1r6n25.dat was generated for a cartilage thickness of 1.0 mm

and values of E3? = 6 MPa and vgj’ = 0.25. In each table, the numeric values were orga-

nized as follows (this is the r10n5.dat table):

100

3
3

0.129 0.162 0.206
0.225 0.270 0.335
0.757 0.851 0.997

2.003 2.465 3.891
2.048 2.514 3.948
2.553 3.064 4.595

0.5 1.0 4.0
0.2 0.7 3.0

where the first two numbers (3 and 3) indicated the number of values of Egg) and 013 that

were considered in the analysis. The following nine numbers were the equilibrium forces
obtained with the different combinations of Eff; and G13, where each row contained equi-

librium forces for a constant E253) , and each column contained equilibrium forces for a

constant G13. The following nine numbers were the peak forces, ordered in the same way

as the equilibrium forces. The following three numbers represented the values of Egg) con-

sidered in the analyses, i.e. 0.5, 1.0 and 4.0 MPa and the last three numbers represented the

values of G13.

A.4 Procedure used to calibrate the master curve
Details of the generation of force-time master curves according to equation (3.3)
were as follows.

Numerical tests revealed that the relation between the parameter c and the perrne-
ability was highly dependent on the difference between the peak and the equilibrium

force, i.e. fpe - feq. Therefore, the relation between c and permeability was determined for

101

sets of elastic constants matching assumed values of fpe and feq, such that their difference
covered a wide spectrum (0.3, 0.5, 1.0 and 2.0 N). Corresponding values of the pairs fpe
and feq were chosen similar to those found in experiments on rabbit patellae to be

described later. These were (0.5, 0.2), (0.7, 0.2), (1.25, 0.25), (2.4, 0.4) N, respectively.

Three elastic sets with associated max, med and min values of E251) were deter-

mined for each assumed pair (fpe, feq). This resulted in 12 elastic sets.

The program Abaqus was run using each of the 12 elastic sets and six assumed
permeabilities 0.5, 2, 5, 10, 18 and 25 m4/(N.s)x10'15 covering a wide range of permeabil-

ities most likely to be found for cartilage. This produced 72 force-time curves.

For each of the 12 elastic sets, values of or and a were chosen by trial and error so
that the master curve had the best fit of the Abaqus force-time curves for all the permeabil-
ities.

With these values of 01 and a for each of the 12 elastic sets, each of the force-time

curves from Abaqus (one for each assumed permeability) was fit with the master curve by

changing the parameter c. This defined a tabular relation between permeability and c.

Values of the parameters a, 01 and c obtained in the calibration procedure were
saved in six files. The names max1.dat, medl .dat and min1.dat were given to the files con-

taining the parameters obtained by using the elastic sets with high, medium and low values

of Egsl) , respectively, for the 1 mm diameter indenter. The names max2.dat, med2.dat and

min2.dat were given to the files containing the parameters obtained using the elastic sets

with high, medium and low values of E381) , respectively, for the 1.5 mm diameter indenter.

102

In each table, the numeric values were organized as follows (this is the maxl.dat table):

6
4

0.30 0.50 1.0 2.2
4.0 12.0 20.0 30.0
1.0 1.05 1.02 1.02

0.005 0.0193 0.046 0.092 0.16 0.23
0.0029 0.011 0.028 0.052 0.1 0.14
0.0055 0.02 0.052 0.091 0.184 0.25
0.006 0.024 0.06 0.11 0.2 0.28

where the first number indicates that the calibration was undertaken for six values of per-

meabilities (0.5, 2.0, 5.0, 10.0, 18.0 and 25.0 111015 m4/(N s)). The second number indi-

cates that four values of the difference (fpe - feq) were used in the calibration. The first row
contains the values of the differences (fpe - feq) used in the calibration procedure. The sec-
ond row contains the values of a corresponding to each value of the difference (fpe - feq).
The third row contains the values of 01 corresponding to each value of the difference (fpe -
feq). The last four rows, one for each value of (fpe - feq), contain the parameter c, which is

related to the permeability.

103

APPENDIX B

SONIE EQUATIONS FOR THE HYPO-ELASTIC AND HYPERELASTIC
MODELS

B.l Integration of the grade zero hypo-elastic equation for the simple extension prob-
lem
For the simple extension problem and an accelerationless motion the velocity

components vi in direction i in terms of the spatial coordinates x1, x2 and x3 are given by

the expression

-k—xi B1
vi_ i1+kit (')

where no sum convention applies, ki is a constant and t is time. Differentiation of equa-
tions (B.1) respect to the spatial coordinates yields the rate of deformation dkl (equal to the
symmetric part of the velocity gradient) which is diagonal in this case. In addition, for this
simple extension case the Zaremba-Jauman stress rate is equal to the time derivative of the

Cauchy stress, and therefore, the component 1 of equation (4.1) becomes

do11 kl k3

_— 2
dt " “”1+klt+H11221+k2t+H11331+k3t' (Bl)

Which can be integrated respect to time to yield equation (4.2). The same procedure

applies for equations (4.3) and (4.4). Note that 1, = l + kit.

104

B.2 Relations between the engineering constants and the coefficients of the strain-
energy function

The following five independent theoretical problems were used to relate the engi-
neering constants to the coefficients of the strain-energy function. All of these calculations
were carried out by implementing equation (4.23) in the notebooks presented in Appendix

B.4 . Zero pressure was assumed in the fluid, which is the condition for equilibrium.

1. Confined stress in direction 3, i.e. 11 = 12 = 0. The derivative of the stress

respect to stretch evaluated at the undeformed configuration gives the equation

E(S)2[1_ v(13)]

E(38)(1_ v(18))_ 2E1?" (351)2

 

- _4a0(2(a3+a5+a6+a7)+n) (B.3)

where the term in the left represents the consolidation coefficient in the direction of load-
ing.

2. Equibiaxial loading in the plane of isotropy with zero strain in direction 3, i.e.
11:19, 13:1. The derivative of stress respect to stretch evaluated in the undeformed con-
figuration yields the equation

E:)E(lsl)

 

=2a0(—a1 + 8a3 + 3n) (BA)

3. Pure shear in the plane 1-2. The derivative of shear stress respect to shear defor-

mation evaluated in the undefonrred configuration yields the equation

= a0(al +n) (B5)

4. Pure shear in the plane 1-3. The derivative of shear stress respect to shear defor-

105

mation evaluated in the undeformed configuration yields the equation

Gl3 = a0(a1 + 2a7 + 11) (B6)

5. Triaxial loading assuming equal stretches, i.e. 11:12:13. The derivative of the

stress in direction 3 respect to stretch evaluated in the undeformed configuration yields the
equation

( ( ( ) ( ) ( )
E33)[E3s3)(1- V152 ) "' 2E1: V381]

(S) ( ) (s) (s)2
E33(1“’182)—2Err"31

 

= 4210(—al +2(3a3+a5+2a6+a7+n)) (37)

B.3 Analysis of uniaxial impact loading according to the hyperelastic formulation

For the case of uniaxial loading in the 3 direction, 11 = 712. In order to satisfy the

condition of free stress on the lateral surfaces it was required that the normal components
of stress in directions 1 and 2 obtained with equation (4.23) were equal to zero. This pro-
vided an equation to calculate pressure in terms of stretches. In addition, employing the

incompressibility condition of equation (4.19), the stretch in direction 1 was found in
terms of the stretch in direction 3 as 711 = 13"”. Finally, pressure was replaced in the

normal component in direction 3 to obtain the Cauchy stress in the direction of loading in

terms of the stretch 13.

A similar procedure was followed for the case of uniaxial loading in the 1 direc-

tion. For this case the lateral stretches M and 13 were not equal, and an additional equation

was solved in order to guarantee free stress lateral surfaces.

The infinitesimal effective elastic constants E11 and E33 were expressed in terms of

the coefficients of the strain-energy function through the equations (B8) and (B9) pre-

106

sented below.

E33 = a0(3a1+ 8(a5 + a7) + 3n) (B.8)

_ 2a0(a1+ n)(3al + 8(a5 + a7) + 3n)

E“ _ 2a1+4(a5+a7)+2n (8'9)

 

By replacing the coefficients of the strain-energy function in terms of the infinites-
imal elastic constants of the solid skeleton (as defined by equations (B-3)-(B-7)), equa-
tions (B8) and (B9) reduce to equations (2.20) and (2.21) which were obtained using

infinitesimal stress-strain relations.

B.4 Description and list of notebooks used in the nonlinear analyses

Notebook ti2.nb was used to analyze the uniaxial stress problem in direction 3 for

the equilibrium case discussed in 4.2.2.

(*"' Uniaxial behavior in the preferred direction **)

Clear[aO,al,a2,a3,a4,a5,a6,a7,a8,n,il,i2,i3,i4,15]

w0 = a0 Exp[
a1(i1-3)+a2(12-3)+a3(il-3)"2+a4(i4-1)+

a5 (i4 - 1)"2 + a6(i1 - 3)(i4 - 1) + a7 (15-1) + a8(iS-1)"2];

w= w0/(i3"n);

wl = D[w,il];

w2 = D[w,i2];

w3 = D[w,i3];

W4 = D[W,i4];

w5 = D[w,iS];

(u- Set up the analysis of the uniaxial problem **)

ii = {{1,0,0},{0,1,0},{0,0,1}};
f= {{sl,0.0}.{0.Sl.0}.{0,0,s3} };
ft= Tmnsposelfl;

X4 = {{0,0,0},{0,0,0}.{0,0.1}};

x5 = 283 "2 x4;

11 = s3"2 + 2 s1"2;

12 = 2s3"2 sl"2 + sl"4;

i3 = s3"2 sl"4;

i4 = s3"2;

15 = s3"4;

(" wc is the derivative of w respect to c for the uniaxial case in the

107

preffered direction **)

n=al +2a2;
84 = -287;
wc=w1ii+(i1ii-c)w2+(i2ii-ilc+c2)w3 +x4 w4+ x5 w5;

(** Now we obtain the cauchy stress t **)
j = PowerExpand[Sqrt[i3l]:
[ = (2/j)f.WC.f[;

(** Now solve s1 in terms of $3 by making the entry 1.1 equal to zero **)
Simplifyltltum
Out[137]=
2 a0
E\"<<-2> a7 <(<-1> + s3\"2)> + a5 <(<-l> + $3\"2)>\"2 +
31 <(<-3> + 2 sl\"2 + $3\"2)> +
a6 <(<-l> + s3\"2)> <(<-3> + 2 sl\"2 + s3\"2)> +
a3 <(<-3> + 2 sl\"2 + s3\"2)>\"2 +
a2 <(<-3> + sl\"4 + 2 sl\"2 s3\"2)> +
a7 ((<-l> + s3\"4)> + 88 <(<-l> + S3\"4)>\"2> Sl\"2 S3
<(Sl\"4 S3\"2)>\"«-l> - 81 - 2 82>
((81 <(<-l) + s1\"2)> +
sl\"2 <(<-6> a3 - a6 + 4 a3 s1\"2 + 2 a3 s3\"2 + a6 s3\"2)
> + a2 <(<-2> + s1\"4 + sl\"2 s3\"2)>)>>
taux =
1(al ((<-l> + s1\"2)> +
sl\"2 <(<-6> a3 - a6 + 4 a3 sl\"2 + 2 a3 s3\"2 + a6 s3\"2)
> + 82 <(<-2> + sl\"4 + sl\"2 S3\"2)>)>; >>
Solve[taux=0,sl];
res = %;
lat = $1 /. res[[4]];

(*"' Check that res[[4]] is one at s3->l **)

Simplify[lat /. s3->l]

Out[l42]=

(“m-81> - 82 + 4 83 + ”((81 + 3 82 + 4 83)>\"2>/<82 +
4 a3»/2>

PowerExpand[%]

\!<‘<2 a2 + 8 83/52 “<82 + 4 83>»

Simplify[%]

Out[ 144]:

1

PowerExpand[%];

("”" Calculate the Poisson’s ratio n3l(s) **)
Simplify[D[lat.s3] /. s3->l];
PowerExpand[%];

Simplify[%]

\!«-«a2 + 2 a3 + a6»/<al + 3 a2 + 4 a3>>>>
Simplifyltll3.3]ll;

taux = % /. sl->lat;

taux = %:

(*"' The uniaxial behavior is represented in the variable taux for one set of coefficients *"')
Simplify[taux /.

{80 -> 0.08,a1->7.733,a2 ->-3.567.a3->l.905,a4->3.958.a5->0.8904,

a6->-0.01 1 l,a7->-1.979,a8->0}];
plot[taux.{s3,0.6,l.4}

(l***********************************************************************************)

108

Notebook tiZe.nb was used to analyze the uniaxial stress problem in the direction 3

for the rapid loading case discussed in 4.2.2.

(""" Uniaxial behavior in the preferred direction, rapid loading case.
i.e. no change of volume **)

w0 = 80 Exp[
81(i1-3)+32(12-3)+a3(il-3)"2+a4(i4-1)+
85 (i4 - l)"2 +a8(i5-l)"2 + a6(il - 3)(i4 - 1)+ a7 (i5-l)];

w= w0/(i3"n);
wl = D[w,il];
w2 = D[w,i2];
w3 = D[w,i3];
w4 = D[w,i4];
w5 = D[w,i5];

(*"' Set up the analysis of the uniaxial problem **)

ii = {11.0.0}.{0.l.0}.{0.0.1l};
f= {{s1,0,0},{0,sl,0},{0,0,s3} };
ft = Transposem;

c = ft.f;

c2 = c.c;

x4 = {{0,0,0},{0,0,0},{0,0,l}};
x5 = 233 "2 x4;

i1= s3"2 + 2 sl"2;

i2 = 2s3"2 sl"2 + sl"4;

i3 = 33"2 sl"4;

i4 = s3"2;

15 = s3"4;

Solve[i3 = 1.31];

S] /. %[[4]];

81 = %;

(** we is the derivative of w respect to c for the uniaxial case in the

preferred direction **)

n=a1+2a2;
a4 = -2a7;
wc=wlii+(ilii-c)w2+(i2ii-i1c+c2)w3 +x4 w4+ x5 w5;

(** Now we obtain the cauchy stress 1 *"‘)
j: PowerExpand[Sqrt[i3]]

Out[36]=

l

In[37]:=

I = (2/j)f.WC.fI;

(H NOW, pressure p is equal to the entry 1,1 **)
p.= 01.1.11]:

Srnrplrfy1p]

109

Our[40]=
\!«l/s3\"2<(
2 a0 E\"<a2 <(<-3> + l/s3\"2 + 2 $3» -
2 a7 <(<-1> + s3\"2)> + a5 <(<-1) + s3\"2)>\"2 +
a1 <(<-3> + %3 + s3\"2)> +
86 <(<-l> + s3\"2)> <(<-3> + %3 + S3\"2)> +
a3 <(1-3> + %3 + s3\"2)>\"2 + a7 <(<-l> + s3\"4)> +
88 <(<-1> + s3\"4)>\"2> <(<-1> + s3)>
«(s3 <(<-a1> + a6 + a6 s3)> + a2 <(<-l> - s3 + s3\"2)> +

2 a3 <(<-2> + S3 + S3\"2)>)>)>»

(** Substract p to the entry 3.

3 to find the total stress on the equivalent material at time zero **)

“[3931] = 'P + “[3931];
1113.311 = Simplifyltll3.3111

Out[43]=
\!<<l/S3\"2<(
2 a0 E\"<a2 <(<-3> + l/s3\"2 + 2 s3)» -
2 87 <(<-l> + S3\"2)> + 85 <(<-l> + s3\"2)>\"2 +
81 <(<-3> + %3 + s3\"2)> +
36 <(<-l> + S3\"2)> <(<-3> + 2/83 + s3\"2)> +
83 ((<-3> + 2/83 + S3\"2)>\"2 + 87 <(<-l> + s3\"4)> +
88 <(<-l> + 33\"4)>\"2> <(<-l> + 83))
1(a2 «(l + 33 + s3\"2)> +
2 a3 <(<-l) + S3)>\"2 ((2 + 3 S3 + 3 S3\"2 + S3\A3)> +
$3 <(al - a6 + a1 s3 - a6 s3 + a1 s3\"2 - 2 a6 s3\"2 +
2 a5 s3\"3 + 2 a6 s3\"3 + 2 a7 s3\"3 + 2 35 s3\"4 +
2 a6 s3\"4 + 2 a7 s3\"4 + 4 a8 s3\"5 + 4 a8 s3\"6 +
4 a8 s3\"7 + 4 a8 s3\"8)>)>)».

(** Calculate the Young’s modulus for rapid loading expressed in equation 8.8 **)

D[t[[3,3]],s3] /. s3->1

2 a0 (3 al+3 a2+4 a5+4 a7+16 a8)
Simplify[%]

2 80 (3 al+3 a2+4 (a5+a7+4 a8))

(""'I tax contains the uniaxial stress response in direction 3 for rapid loading **)

tax = 1113.31];

(*‘Fllflhlllllllhlukll'lfllllhlfllli1"!***********************Ill$110108IF********$***************'lfllfllnll’ll’l'*******)

Notebook tiapp.nb was used to derive the right hand side of equations (B.3), (B .4)

and (B7).

(** Derivation of right hand side of equations B.3-B.6 **)

w0 = 80 Exp[
81 (i1-3)+a2(12-3)+a3(i1 -3)"2+a4(i4- 1)+
a5 (i4 - 1)"2 +a6(il - 3)(i4 - 1) + a7 (i5-l)];
w= w0/(i3"n);
wl = D[w,il];
w2 = D[w,i2];
W3 = D[w,i3];

110

w4 = D[w,i4]:
w5 = D[w,iS];

(** Set up the analysis of the simple extension problem **)

ii = {{1,0,0}.{0,l,0}.{0.0,l}};
f= {{sl,0,0},{0,32,0},{0.0,s3}};
ft = Transpose[f];

c = ftf;

c2 = c.c;

x4 = {10.0.0}.{0.0.0}.{0.0.1}l;
x5 = 2s3 "2 x4;

11 = sl"2+ 32"2 + s3"2;

12 = sl"2 52"2 + sl"2 s3“2 + s2"2 s3"2;
i3 = sl"2 s2"2 8392;

i4 = s3"2;

i5 = s3"4;

(""'I we is the derivative of w respect to c for the uniaxial case in the
preffered direction **)

82=(n-81)/2;
a4=-2a7;
wc=wl ii+(ilii-c)w2+(iZii-ilc+c2)w3+x4 w4+x5 w5;

(*"‘ Now we obtain the cauchy stress 1 **)

j = PowerExpand[Sqrt[i3]];

1= (2/j)f.wc .11;

(’"' Determines right hand side of equation 8.3 **)
s1 = s2;

D[t[[3,3l] /.s2->1,s3] /. s3->l;

Simplify[%]

4 a0 (2 a3+2 a5+2 a6+2 a7+n)
(** Determines the right hand side of equation 8.4 **)

D[t[[l,l]] I. s3->l, s2] /. s2->l;
Simplify[%]

-2 a0 (31-8 83-3 n)
(** Determines the right hand side of equation 8.7 *"‘)

Clear[s1,52,s3];

31 = s2;

53 = s1;

D[t[[3,3]],s3] /. s3->l;
Simplify[%]

4 a0 (-a1+2 (3 a3+aS+2 a6+a7+n))
(*****************II"!IIUIUIUIUIUIUII******************#***$***********************************)

Notebook ti4.nb to determine the right hand side of equation (B.6).

Clear[80,al ,a2,a3,a4,85.a6,a7,i 1 .12,i3,i4,i5]

lll

(*"' Analysis of the pure shear in the 2-3 plane **)

w0= a0 Exp[
a1 (il -3)+a2 (12-3)+a3 (il -3)"2+a4 (i4 -1)+
a5 (14 - 1)"2 + a6(il - 3)(i4 - 1) + a7 (15-1)];

W: w0/(i3"n);
wl = D[w,il];
W = D[WJZ];
w3 = D[w,i3];
w4 = D[w,i4];
w5 = D[w,iS];

(** Set up the analysis of pure shear problem in the 2-3 problem ”"0

ii = {{1,0,0},{0,1,0},{0,0,1}};

f: {I190v0}9{001k}9109091}};

ft= Transpose[f];

c = ft.f;

c2 = c.c;

x4 = {{0,0.0},{0,0,0},{0.0,1}};

x5 = {{0.0.0}.{0.0.c[[2.3ll1.1001132112611331“ };
Eigenvalues[c]

Out[19]=

\!<{ 1, 1/2 <(2 + k\"2 - k ‘14 + k\"2>)>,
1/2 «(2 + k\"2 + k ‘14 + k\"2>)> }.

m2 = %[[211

“(I/2 ((2 + k\"2 - k ‘4 + k\"2>)>>

m3=%%[[3l]

\!<l/Z «(2 + k\"2 + k “<4 «1» k\"2>)»

m1 = 1;

i1 = l + m2 + m3;

12=m1m2+m1m3+m2m3;

i3 = l;

x4.c

Out[32]=

\!<{ {0. 0, 0}. {0, 0, 0}, {0, k, 1+ k\"2} })

i4 = 20113.31];

x4.c2

Out[34]=

\!<{ {0, 0, 0}, {0, 0, 0}, {0, k + k <(l + k\"2)),

k\"2 + <(1 + k\"2)>\"2} )1

i5 = %[[3.311;

a2 = (n - a1)/2;

a4 = -2a7;

wc=wlii+(i1ii-c)w2+(i2ii-il c+c2)w3 +x4 w4+x5 w5;
(** Now we obtain the cauchy stress t **)

j = PowerExpand[Sqn[i3]];

t = (2/j)f.wc.ft;

(** Now determine the infinitesimal shear modulus **)

D[t[[2.3]l.k] /. k->0;
Simplify[%]

112

Out[45]=
a0 (al-r-2 a7+n)

(******************3101'#310101"?*ilt'lfllulI*Il‘*IIUII*III*iIUIUII’II*ilfllfll"1‘*1"!**********ll"!**********************)

Notebook ti5.nb was used to determine the right hand side of equation (B.5).

(*"‘ Analysis of the pure shear problem in the 1-2 plane **)

w0 = 80 Exp[
81 (i1 -3)+82(12-3)+a3 (il -3)"2+a4(i4- 1)+
85 (i4 -1)"2 + a6(il - 3)(i4 - 1) + a7 (i5-l)];
w= w0/(i3‘n);
wl = D[w,il];
w2 = D[w,i2];
w3 = D[w,i3];
w4 = D[w,i4];
w5 = D[w,iS];

(""" Set up the analysis of pure shear problem in the 2-3 problem **)

ii = {{1,0,0},{0,1,0},{0,0,1}};

f= {{1k,0},{0,1,0},{0,0.1}};

ft= Transpose[f];

c = ftf;

c2 = c.c;

x4 = {{0,0,0}.{0,0,0}.{0.0,1}}:

x5 = l{0.0.0}.{0.0.c{[2.3lll.{0.c[[3.211.20113.3ll} }:

Eigenvalues[c]

Out[19]=

\!<{ 1, 1/2 <(2 + k\"2 - k “<4 + k\"2>)),
1/2 (2 + k\"2 + k “<4 «I» k\"2>))}>

1112 = %[[211

“<14 <(2 + k\"2 - k “<4 + k\"2>)»

m3=%%[[3]]

\1<W <(2 + k\"2 + k ‘(4 + l(\"2>)>>

ml = l;

(*"' Strain invariants il, 12, 13, i4 and i5 **)

11 = 1+ m2 + m3;

12=m1m2+ml m3+m2m3;

i3 = l;

x4.c;

i4 = %[[3311;

x4.c2;

i5 = %[[3311;

(** we is the derivative of w respect to c **)

82 = (n - al)/2;
a4 = -2a7;

wc = wl ii 4» (i1 ii - c)w2 +(12 ii - i1 c + c2) w3 + x4 w4 + x5 wS;
(** Now we obtain the cauchy stress t *"')

I = (2/j)f.WC.ft;

C” Now determine the infinitesimal shear modulus *"‘)

D[t[[1.2]].k] /- k->0

113

Out[38]=

\!<2 <(aO a1+ 1/2 a0 <(<-a1> + n).)»
Simplify[%]

a0 (a1+n)

(*************************$****************#******#**$************************$*****)

Notebook ti1.nb was used to determine the left hand side of equations (B.3), (B.4)
and (B .7).

(** Infinitesimal strain-
stress relations for a transversely isotropic material **)
d ={ { l/ell,-n12/el 1,-n31/e33},{-n12/e1 1.1/e1 1,-n31/e33},{-n31/e33,—n3l/e33.
1/e33} };
Inverse[d];
d=Simplify[%];
(*"' The term [[33]] is equal to the left hand side of equation B.3 **)
6113.31]
Out[8]=
\l<<c33\"2 <(<-l> + n12)>»/<e33 <(<-l> + n12)> +
2 e11 n31\"2»
(*"' The left hand side of equation 8.4 is equal to the sum of the terms
[[1,1]] and [[1,2]] as follows **)
Simplify[dlll.l]l + 4111,2111
Out[10]=
\!<<e11e33»/<e33 - e33 n12 - 2 e11n31\"2»
(** The left hand side of equation 8.7 is equal to the sum of the terms
[13.11]. [13.21] and [13.31] *"‘)
Simplify[d[[3.l]] +d[[3.2]] + (“13.3111
\!<<e33 <(e33 <(<-1> + n12)> - 2 ell n31)
>x/<e33 <(<-l> + n12)> + 2 ell n31\"2>>

(*****************************¥*************$***$************************************)

Notebook prove3.nb used to prove that equations (B.8) and (B9) are equivalent to

equations (2.20) and (2.21)

("'1' Notebook prove3.nb **)
("'1' First, define the right hand side of equations (B.3)-(B.7) "l
otas = 4a0(-a1 + 2(3a3 + a5 +2a6+a7+n));

g12a = a0(al +n);

h3as = 4a0( 2(a3+ a5 + a6 + a7 )+n);

hlas = 2a0(-al + 8a3+3n);

h3 = e3s"2(1 - n12s)/(e3s(1-n12s)-2els n3ls"2);

(""" Define left hand side of equations (B.3)-(B.7) *"')
bl = els e3s l(e3s (l-nl2s) - 2els n315"2);

g12 = els/(2(1+ n12s));

ots = e3s(e3s(l-n125) + 2els n313)/(e3s (l - n12s) - 2els n3ls"2);

(** Define effective Young’s modulus in direction 3 according to equation B.8 **)
e3 = a0 (3a1 + 8(a5 + a7 )+3n);

(*** Solve equations (B.3)-(B.7)in order to find e333,

114

e1 Is and n12s in terms of the parameters of the strain energy function **)
Solve[gl2 = g12a, els];
els I. %;
els = %11111;
Solve[hlas=hl,e3s];
e3s /. %;
638 = %[[l]];
e3s = Simplify[e3s];
Solve[h3as=h3,nl2s];
n125 I. %;
11128 = %[[lll;
n123 = Simplify[nl2s];
Solve[otas==ots,n3 ls];
n31s /. %;
11318 =%[[lll;
(*"' Now input equation (2.20) for the effective Young’s modulus in direction 3 **)
for3 = Simplify[
e3s(1 - 4n313 + 2( l-n12s)(e3s/els))/(2(l-n12s)(e3s/els) - 4n31s"2)];
(*** Check that for3 is equal to e3 **)
Simplify[for3 - e3]
Out[32]=
0
(*** Now input equation (2.21) for the effective Young’s modulus in direction 1 **)
forl = Simplify[
els(l - 4n313 + 2(l-nl2s)(e3s/els))/(
1 + (l-nl2s"2)(e3s/els) - 2n31s(l + n125) - n3ls"2(els/e3s))];

forl=FullSimplify[%]
Out[35]=
\!<<aO «(81+ n)> <(3 81+ 8 <(a5 + a7)> + 3 n)>»/<al +

2 <(a5 + a7)» + n»

(*"‘ Define the Young’s modulus in direction 1 according to equation (8.9) ***)
el = 280(al + n)(3a1 + 8(a5 + a7 H3n)/(2a1 + 4(a5 + a7) + 2n);

("”' Check that forl is equal to e1 *"‘)
Simplify[forl - e1]

Out[38]=
0

(***$*************************II"!#81"!!!*ikllnlfllulflllll‘fll*IBIIUIUIUIUF’II********#***$*****$****)

Notebook tiappnb was used during the iterative procedure discussed in 4.2.2 in
order to, first, determine the coefficients of the strain-energy function according to equa-

tions (B.3)-(B.7) given the infinitesimal elastic constants and the coefficients a1 and 11,

second, analyze the equilibrium and rapid loading cases in direction 3 in order to define

the variations of BS, E33 and vgj’ according to equations (4.13), and third, analyze the

out-of-plane pure shear problem in order to define the variation of G13 with shear strain.

115

(***********************************************#************************************)

(*"' Notebook tiapp.nb **)

Clear[a0,al,a2.a3,a4,a5,a6,a7,n,h3,hl,g12,ha3.hal,na3l,gal3.ga12]
(*"‘ Left hand side of equations (A.4) - (A.7) **)
\!<<h3 =
<e33\"2 <(<-l> + n12)>»/<e33 <(<-l> + n12)> +

2 e11n31\"2>;»
\!<<h1 = <ell e33»/<e33 - e33 n12 - 2 ell n31"2>; >>
hot = e33(e33(-l + n12) - 2elln31)/(e33(-1+ n12) + 2ell n3l"2);
g12 = e11/(2(1+n12));

(*"‘ Define the value of the infinitesimal elastic constants **)
{e33 = 0.8,el l=2.7,n12=0.0.n3 1:0. l3.g 13:0.13};

("”" Right hand side of equations (A.4) - (A.7) **)

ha3 = 4aO( 2( a3 + 85 + a6 + a7 )+n);
hal = 2a0(-a1 + 8a3+ 3n);

haot = 4a0 (-a1 + 2(3a3+a5+2a6+a7+n)) ;
ga13 = a0(al + n + 237)

a0 (a1+2 a7+n)

gal2 = a0 (a1 + n);

(*******$**************$*******************************#**********************)

Clear[al,a2,a3.a4.85.a6.a7.n];

(** Solve for the coefficients of the strain—energy function. given a1 and n **)
{a1 = 5. n=0.7846};

Solve[gal2=g12.a0];
30 I. %[[1]];
80 = %;

Solve[gal3 == g13.a7];
a7 /. %[[1]];
87 = %;

Solve[hal = hl.a3];
a3 /. %[[1]];
a3 = %;

Solve[haot=hot.a6];
a6 /. %[[1]];
a6 = Simplify[%];

Solve[ha3==h3,aS];
35 /- %[[1]];

85 = %;
a8 =0;
82: (n-a1)/2;

{a0,a1.a2,33.a5,a6,a7,a8,n}
{0.233378 .5.-2. 1077, 1 .14696, 1 .3 1999.0.238202,-2.6 l 378.0,0.7846}

116

{h81.h83,11831.g812.g813}
{3.04766.0.90301 1.1183 1.1.35.0.13 }

(*****************IFIII***#II|**********************IIUIUI‘***************************)

(*" t82 is the Cauchy stress in direction 3 for the equilibrium case **)

“((182 =
((4\"(l + 81+ 2 82) 80
E"(((I/(2 ((82 + 4 83))»
<(<-al\"2> + 36 a3\"2 + 8 a3 85 + 12 a3 a6 - a6\"2 +
8 83 a7 + 8 a3 88 - 24 a3\"2 s3\"2 -
16 83 85 s3\"2 - 16 a3 a6 s3\"2 + 2 a6\"2 s3\"2 -
16 83 a7 s3\"2 + 4 a3\"2 s3\"4 + 8 a3 85 s3\"4 +
4 83 a6 s3\"4 - a6\"2 s3\"4 + 8 83 a7 s3\"4 -
16 a3 88 s3\"4 + 8 83 88 s3\"8 -
a2\"2 ((2 + s3\"4)> -
6 83 “<4 <(81 + 2 82)) ((82 + 4 a3)> +
<(81 - a6 + 82 33\"2 + a6 s3\"2 +
2 83 (((-3) + S3\"2))))\"2) -
a6 ‘<4 <(81 + 2 82)) <(82 + 4 83)) +
<(81 - a6 + 82 s3\"2 + a6 s3\"2 +
2 83 (((-3) + S3\"2))))\"2) +
2 83 s3\"2
“(4 ((81 + 2 82)) ((82 + 4 83)) +
<(81 - a6 + 82 s3\"2 + a6 s3\"2 +
2 83 (((-3) + S3\"2))))\"2) +
a6 s3\"2
“<4 <(81 + 2 82)) <(82 + 4 83)) +
<(81 - a6 + 82 s3\"2 + a6 53\"2 +
2 83 (((-3) + S3\"2))))\"2) +
81 <(<-4> 82-483+286+483s3\"2-
2 86 s3\"2 +
.(4 ((81 + 2 82)) ((82 + 4 83)) +
<(81 - a6 + 82 s3\"2 + a6 33\"2 +
2 83 (((-3) + s3\"2)>)>\"2>)> +
82 <(1083+285 +6a6+287+288-
4 85 s3\"2 - 6 a6 s3\"2 - 4 a7 s3\"2 -
2 a3 s3\"4 + 2 85 s3\"4 + 2 a7 s3\"4 -
4 a8 s3\"4 + 2 a8 s3\"8 +
s3\"2 “<4 <(81 + 2 82)) <(82 + 4 83)) +
<(81 - a6 + 82 s3\"2 + a6 s3\"2 +
2 83 <((-3) + s3\"2)>)>\"2>)>)
>)) ((
(1/<(82 + 4 83))\"2>
((S3\"2
<(81 - a6 + 82 s3\"2 + 86 s3\"2 +
2 83 (((-3) + S3\"2)> -
‘(4 ((81 + 2 82)) ((82 + 4 83)) +
<(81 - 86 + 82 s3\"2 + a6 s3\"2 +
2 83 (((-3) + S3\"2))))\"2)))"2)
)))"(((-8l) - 2 82))
<(81 <(82 + 4 a3 - 2 a3 s3\"2 + a6 s3\"2)> +
82\"2 ((2 + s3\"4)> -
S3V‘2 ((
<-12> 83\"2 - 8 83 85 - 8 83 a6 + 86\"2 -
8 a3 a7 + 4 a3\"2 s3\"2 + 8 a3 85 s3\"2 +
4 83 a6 s3\"2 - a6\"2 s3\"2 + 8 a3 a7 s3\"2 -

117

16 a3 88 s3\"2 + 16 83 88 s3\"6 +
2 83 ‘(4 ((81+ 2 82)) ((82 + 4 83)) +
<(81 - a6 + 82 s3\"2 + a6 s3\"2 +
2 83 (((-3) + S3\"2))))\"2) +
86 “<4 <(81 + 2 82)) ((82 + 4 83)) +
((81 - a6 + 82 s3\"2 + a6 s3\"2 +
2 83 ((<-3) + S3\"2))))\"2))) +
a2 ((8 83 + 2 85 s3\"2 + 3 a6 s3\"2 + 2 87 s3\"2 +
2 a3 s3\"4 - 2 85 s3\"4 - 2 87 s3\"4 +
4 88 s3\"4 - 4 88 s3\"8 -
s3\"2 “<4 ((81 + 2 82)) <(82 + 4 a3)> +
<(81 - a6 + 82 s3\"2 + a6 s3\"2 +
2 83 (((-3) + S3\"2))))\"2)))))))/
((S3 ((
81 - 86 + 82 s3\"2 + 86 s3\"2 + 2 83 (((-3) + s3\"2)> -
“(4 ((81 + 2 82)) ((82 + 4 83)) +
<(81 - 86 + a2 s3\"2 + a6 s3\"2 +
2 83 (((-3) + S3\A2))))\"2))))); ))

(******#*****$$*****************$*******#******$******************************)

(*"' tal is the Cauchy stress in direction 3 for rapid loading *"‘)
181 = (l/SBV‘Z)
((2 80 BA
((82 (((-3) + l/S3\"2 + 2 83)) -
2 87 (((-l> + s3\"2)> + 85 (((-l» + s3\"2)>\"2 +
81 (((-3> + %3 + s3\"2)> +
86 (((-l) + s3\"2)> (((-3) + 2183 + s3\"2)> +
83 (((-3) + 2/83 + S3\"2))\"2 +
87 (((-l) + s3\"4)> + 88 (((-l) + S3\"4)>\"2))
(((-l> + 83))
<(82 ((1 + $3 + s3\"2)> +
2 83 (((-l) + S3))\"2
((2 + 3 s3 + 3 s3\"2 + s3\"3)> +
$3 ((81 - 86 + 8183 - 86 $3 + 8133\“2 - 2 a6 s3\"2 +
2 85 s3\"3 + 2 86 s3\"3 + 2 87 s3\"3 +
2 85 s3\"4 + 2 86 s3\"4 + 2 a7 s3\"4 +
4 88 s3\"5 + 4 88 33V‘6 + 4 88 s3\"7 + 4 a8 s3\"8)
))))); \n
(""" p is the pressure during rapid loading *"')
p = (l/S3\"2)
((2 80 EA
((82 (((-3) + l/S3\"2 + 2 s3)> -
2 87 (((-l) + S3\"2)) + 85 (((-l) + S3\"2))\"2 +
81 (((-3) + 2/83 + S3\"2)) +
86 (((-l) + S3\"2)) <((-3) + %3 + S3\"2)) +
83 ((<-3> + %3 + s3\"2)>\"2 +
87 (((-l) + S3\"4)) + 88 (((-1) + s3\"4)>\"2)>
(((-l> + 83))
((S3 (((-81) + 86 + 86 83)) +
82 (((-l) - 83 + s3\"2)> +
2 83 (((-2) + $3 + S3\"2)))))); ))

(*************$$***********************$$*************************************)
(""' 03 is the Cauchy shear stress in the planes 2-3 **)

\!(<123 =
2 80 E\"(k\"2

118

<(81 + 82 + a7 + 83 k\"2 + 85 k\"2 + 86 k\"2 + a7 k\"2 +
9 88 k\"2 + 6 88 k\"4 + 88 k\"6)>> k
<(81+82+a7+283k\"2+285k\"2+286k\"2+
2 87 k\"2 + 18 88 k\"2 +18 88 k\"4 + 4 88 k\"6)> ; >>

(*"‘ Verify that stress is zero in the undeformed configuration **)
Simplify[taZ /. $3->l]
0

(** This> is the lateral stretch s1 in terms of the axial
stretch s3 *"')
“(lat =
‘<<-<81/<82 + 4 83>>> + <6 83/(82 + 4 83> +
86482 + 4 83> - <82 33\"2>/<82 + 4 83> -
<2 83 s3\"2>/<82 + 4 83> - <86 s3\"2>/<82 + 4 83> +
‘<<-4) <(<-8l> - 2 82)) <(82 + 4 83)) +
<(81 - 6 83 - a6 + 82 s3\"2 + 2 83 s3\"2 + 86 s3\"2)
>\"2»/<82 + 4 83>)” 2)
(** der2 and derl are the Young’s moduli for the equilibrium and rapid loading cases, respectively,
calculated according to equation (4.13) as a function of stretch, der3 is the out-of-plane Poisson's
ratio. cauculated according to equation (4.9) given the lateral deformation. contained in lat **)

derl = 33 D[t81.s3];

der’Z = 83 D[t82,83];

der3 = -s3 D[lat,s3] flat;

(** In the undeformed configuration are equal to the infinitesimal values **)

{derl.de{2.der3} /. s3->1
{ 163445080 13}

{der1,der2,der3} I. s3->0.65
{22.7181,1.27907 00527634}

{derl.det’2.dcr3} /. s3->0.75
{6.5401.0.964808,0.070948 }

{derl.der2.der3} I. s3->0.85
{2.90996.0.813855,0.0921555}

{derl.der2,det3} /. s3->O.95
{ 1.83184.0.778448,0.116551 }

{der1,de12.der3} /. s3->1.05
{ l.56418.0.84913,0.1443l4}

{der1.der2.der3} /. s3->l.15
{1.70437,l.04477,0. 175634}

{derl,der2.der3} /. s3->1.25
{2.19l36.l.41624.0.210699}

{der1,der2.der3} /. s3->1.35
{3.1 1671.2.06253.0.249689 }

{derl.der’2.der3} /. s3->1.45
{4.7 135.3. 16569.0.292754 }

(""'I der3 is equal to G13. calculated in terms of shear deformation *"‘)

119

der3 = D[t23,k];

(*"' In the undeformed configuration is equal to the infinitesimal Gl3 **)
der3 /. k->0

0.13

der3 /. k->0.15
0.138305

derB /. k->0.25
0.153601

der3 /. k->0.35
0.1779

der3 /. k->0.45
0.213069

der3 /. k->0.55
0.26204

der3 /. k->0.65
0.329303

der3 /. k->O.75
0.421692

(etwtemwmtemeweem**##m*******m****$*********#*************t******#*****************)

Notebook tiqua.nb was used to calculate the variations of the engineering param-

eter 13‘]? with stretch as required in the iterative procedure discussed in 4.2.2.

(*****************Ihillll‘fi'flllll'lfllllfllfl'*********************IF’IUII*******************************)

(*"' XX. yy and p82 come from the solutions of the equations t[[2.2]]=t[[3,3]]=0
which were used to enforce stress-free lateral forces in this uniaxial problem in the
plane of isotropy H)

(*"' xx and yy are operated manually in order to be able to solve this problem *"‘)

\!<<pa2 =
‘<<81 - 81 33\"2 + 82 <(2 - sl\"2 s3\"2)> -
S3\A2 <(
<-6> 83-285-486-287+2a3 sl\"2+
a6 s1\"2 + 2 83 s3\"2 + 2 85 s3\"2 + 2 a6 s3\"2 +
2 a7 s3\"2 - 4 88 s3\"2 + 4 a8 s3\"6)>>/<<(
82 + 2 83 + 86)) S3\"2)); >)
\!<< =
<81 - 81 s3\"2 + 82 <(2 - sl\"2 s3\"2)> —
S3V‘2 <(

<-6>83-285-4a6-287+28381\“2+a6sl\"2+
2 a3 s3\"2 + 2 85 s3\"2 + 2 a6 s3\"2 + 2 a7 s3\"2 -
4 a8 s3\"2 + 4 a8 s3\"6)>>/<<(82 + 2 83 + a6)>
S3V‘2); ))

120

\!«yy =
<( yy)) ((4 83)) -
<( <-8l> + 86 - 82 s1\"2 - 82 s3\"2 - a6 s3\"2 -
2 83 <(<-3> + SIV‘Z + S3\A2)))); ))

W = ”“2:

\!<<xx =
8 <(81 + 2 82)> 83 +
<(81 - a6 + 82 sl\"2 + a2 s3\"2 + 86 s3\"2 +
2 83 <(<-3) + sl\"2 + S3\"2)> ))\"2; ))

(*"‘ Now write both xx and yy in termos of numeric values of parameters **)

{xx,yy} /.
{80->0.0708959,81->7.3 ,82->—2.92739,83->1.56851,85->l.64,86->0.394l45,

a7->-3.45577.88->0.0} ;
{xl.yl} = %;

(*"' Verify that they are equal in the underformed configuration **)
{XLyl} I. {sl->1.s3->l}

{209969209969 }

x1 = Simplify[xl];
yl = Simplify[yl];

(n Now solve x1 == yl in order to obtain s3 in termos of $1 **)

Solve[xl == yl.s3];
sol = 33 /. %;
sol[[3]] /. s1->l

(*"‘ Verify that $3 is also equal to 1 when 31 is equal to l *"‘)
1.00000000000000004515194661153

(*¢****m***#**********$**#****************************$******************m**)

(*"‘ Bring now the axial force t[[1,1]] in terms of 31, $2 and 33 M)

\!<<t8xial =
2 80 E\"<<-2) 87 <(<-1) + S3\"2)> +

85 <(<-1> + S3\"2))\"2 +
81 <(<-3> + sl\"2 + 82\"2 + s3\"2)> +
86 <(<-l) + S3\"2)) <(<-3) + Sl\"2 + 32\"2 + s3\"2)> +
83 <(<-3> + sl\"2 + s2\"2 + s3\"2)>\"2 +
87 <(<-l> + S3\"4)> + 88 <(<-l) + S3\"4))\"2 +
82 <(<-3> + S2\"2 s3\"2 + Sl\"2 <(s2\"2 + S3\"2))))> s1

82 s3 <(sl\"2 $2\"2 s3\"2)>\"<<-1> - 81 - 2 82>

<(81 <(<-1) + sl\"2)> +
sl\"2 ((

<-6> 83 - a6 + 2 83 sl\"2 + 2 a3 s2\"2 + 2 83 s3\"2 +
a6 s3\"2)> + a2 <(<-2> + sl\"2 52\"2 + sl\"2 s3\"2)>)

>; ))

121

(*"' First replace s2 in terms of 51 and s3 **)
tax =taxial /. 82->p82;

(*"' Verify that it is zero in the undeformed configuration **)
tax /. {sl->l,s3->l}
0

tax =tax /.
{80->0.0708959,81->7.3,82->-2.92739.83->1.56851,85->1.64.a6->O.394145.
87->-3.45577.a8->0.0};

tax /. {sl->l,s3->l}
\!<<-1.5936416850436216“-16»

(** Now, replace s3 in terms of 51 **)

tax = tax I. s3->sol[[3]];

tax = Simplify[tax];

(** Verify that it is zero in the undeformed configuration **)
tax I. s1->l

\!<<-2.51872833700873943‘*"-l6>>

(** Now calculate the engineering parameter El l(s) according to 8 formula similar to (4.13) for
different values of the stretch **)

1.0 ((tax /. sl->l.0005) - (tax /. sl->0.9995))/0.001
1.24

O.65((t8x /. sl->0.6505) - (tax /. sl->0.6495))/0.001
3.29033

O.75((lax /. sl->0.7505) - (tax /. sl->0.7495))/0.001
1.88555

0.85((t8x /. sl->0.8505) - (tax /. sl->0.8495))/0.001
1.25618

095((t8x /. sl->0.9505) - (tax /. s1->0.9495))/0.001
1.1 1377

l.05((t8x /. sl->l.0505) - (tax /. sl->l.0495))/0.001
1.5682

l.15((tax I. sl->1.1505) - (18x /. sl->1.l495))/0.001
3.49511

(#31010!It!!!***********************It'll*****************************************Ilulnlnlnlfllulnhluk)
Notebook tilpeakrmb was used to calculate the variation with stretch of the engi-

neering parameter E11 as required in the iterative procedure discussed in 4.2.2.

(**********#*********$***********************¥***$******************$*************)

((*"' t8 is the term of the lateral stress which has to be zero to guarantee 8 stress free lateral
surface *"‘)

122

\!<(Ufl ==
<-2> 83 -
52\"2 <(81 - 6 a3 - 86 + <(82 + 2 a3)> 52\"2)> s3\"2 +
s2\"4 <(
81-2 <(3 83+85+286+a7)>+
<(82 + 2 83 + a6)> s2\"2)> s3\"6 +

2 <(83 + 85 + a6 + a7 - 2 88)> sZ\"4 s3\"8 +

4 88 S2\"4 S3\"12; ))
(** x is made equal to ta for numeric values for the coefficients **)
x = t8 I. {80->0.0708959,al->7.3.82->-2.92739,83->1.56851.85->1.64.

a6->0.394145.87->-3.45577,88->0.0} ;
x=FullSimplify[x];

(** With this equation we obtain 53 in terms of s2 in order to guarantee zero stress in direction 3 **)
Solve[x==0.s3];

res=s3 /. %;

(** Check that the res[[2]] is equal to l for 52 ->1 **)

res[[2]] I. $2->l

1.

(** The term t22 is the total stress in the plane of isotropy. in terms of s3 and $2 **)

\!<<t22 =
<1/<s2\"4 s3\"4>>
<(2 80 EA
<(<-2> 87 <(<-l) + S3\"2)) +
85 <(<-1) + S3\"2)>\"2 +
81 <(<-3) + S2\"2 + 1/<s2\"2 s3\"2> + S3\"2)) +
86 <(<-l) + S3\"2))
<(<-3> + 32\"2 + 1/<82\"2 s3\"2> + s3\"2)> +
83 <(<-3> + s2\"2 + 1/<52\"2 s3\"2> + s3\"2)>\"2 +
82 <(<-3> + 1/52\"2 + 1/S3V‘2 + 52\"2 s3\"2)> +
87 <(<-l> + S3\"4)) + 88 <(<-1) + s3\"4)>\"2)>
<(<-l) + S2\"4 S3\"2)>
<(32\"2 s3\"2 <(81 - 86 + 82 s3\"2 + 86 s3\"2)> +
2 a3 <(l + s2\"4 s3\"2 + $2\"2 s3\"2 <(<-3> + s3\"2)>)
>)))); ))
(** t2 if the total stress in the plane of isotropy in terms of numerical values for the coefficients **)
t2 = Q2 I.
{80—>0.0708959,a1->7.3,82->-2.92739,83->1.56851,85->l.64.
86->0.394l45,a7->-3.45577,a8->0.0} ;

(** Now, we obtain the uniaxial response by replacing s3 in terms of 52 according to res[[2]] **)
t2 = t2 /. s3->res[[2]];

t2 = Simplify[t2];

(** Verify that stress is zero in the undeformed configuration **)

t2 /. 52->1

0.f

(** These are the variations of the engineering parameter E11 in terms of the stretch **)

({t2 l. $2->l.0. l.0((t2/. $2->l.0005) - (t2/. 82->0.9995))/0.001}
{0.,1.4197}

{t2 /. 82->0.65.0.65((t2 I. $2->0.6505) - (t2 /. 82->0.6495))/0.001}
{-l.5541,12.3897}

{t2 /. $2->0.75.0.75((t2 I. 82->0.7505) - (t2 /. 82->0.7495))/0.001}

123

{-0.544092.3.73941 }

{t2 /. $2->O.85,0.85((t2 I. sZ->O.8505) - (t2 /. $2->0.8495))/0.001}
{-0.22725.l.71695}

{t2 /. $2->0.95.095((12 I. 52->0.9505) - ([2 /. 82->O.9495))/0.001}
{-0.0685001.1.29136}

{t2 I. 82->l.05.1.05((t2 I. sZ->l.0505) - (t2 I. sZ->l.0495))/0.001}
{0.0777403.1 .82614}

{t2 /. 82->l.15,1.15((t2/. sZ->1.1505)-(12/.32->l.l495))/0.001}
{0.331 159.4.3172}

B.5 Subroutine hypela.f used to define the components of the hypo-elastic tensor in
the program MARC

This is the subroutine hypela (contained in the file hypela.f) in which the variation
of the engineering parameters is specified as a function of deformation. The variation of
the engineering parameters can be obtained by employing the notebooks tiapp.nb,

tiqua.nb and tineaknb listed in B.4.

SUBROUTINE HYPEI.A(D.G,E,DE,S.TEMP.
1 DTEMP.NGENS.N,NN,KC,MATS,NDI.NSHEAR)
IMPLICIT REAL *8 (A-H, O-Z)
DIMENSION D(NGENS.NGENS).G(NGENS).E(NGENS),
l DE(NGENS).S(NGENS).TEMP(l),DTEMP(l),N (2)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

Definition of the engineering parameters

e3 : Young's modulus in the axial direction

e1 : Young‘s modulus in the plane of cartilage
v31 : Out-of-plane Poisson‘s ratio

v12 : In-plane Poisson’s ratio

gl3 : Shear modulus

000000000

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Definition of the out-of-plane Poisson‘s ratio for the rapid loading case
v31 = 0.49
C Definition of the Young’s modulus e1 in the plane of isotropy in terms of deformations,
C contained in the array E(I). These are logarithmic deformation which need to be converted
C to stretch, that’s why the exp function is used here
es = E(2)
s1 = (exp(es) + exp(E(3)))/2.0
if (sl.gt.1.3) then
el = 41.79
elseif (s 1 . gt. 1 .2) then
el = 16.32

124

elseif (sl.gt.1.l) then
el = 7.32
elseif (sl.gt.1.0) then
el = 3.91
elseif (sl.gt.0.9) then
el = 2.74
elseif (s 1 .gL0.8) then
el = 2.86
elseif (sl.gt.0.7) then
el = 4.56
elseif (s l .gt.0.6) then
el = 10.63
elseif (sl.gt0.5) then
el = 38.57
else
el = 45.5
endif
C Definition of the Young’s modulus e3 in the axis of loading
es = E(l)
$1 = exp(es)
if (sl.gt.1.3) then
e3 = 7.23
elseif (sl.gt.1.2) then
e3 = 4.13
elseif (sl.gt.1.l) then
e3 = 2.60
elseif (s 1.gt. 1.0) then
e3 = 1.94
elseif (sl.gt.0.9) then
e3 = 1.91
elseif (sl.gt.0.8) then
e3 = 2.62
elseif (s l .gt.0.7) then
e3 = 4.96
elseif (s 1 .gt.0.6) then
e3 = 13.3
elseif (sl.gt.0.5) then
e3 = 60.0
else
e3 = 122.2
endif
C Definition of the modulus g13 according to shear strain
es = 5(4)
if (es.lt.(-1.1).or.es.gt.1.1) then
g13 = 8.47
elseif (es.lt.(-1.0).or.es.gt.1.0) then
g13 = 5.08
elseif (es.lt.(-0.9).or.es.gt.0.9) then
g13 = 3.18
elseif (es.lt.(-0.8).or.es.gt.0.8) then
g13 = 2.07
elseif (es.lt.(-0.7).or.es.gt.0.7) then
g13 = 1.39
elseif (es.lt.(-0.6).or.es.gt.0.6) then
g13 = 0.96
elseif (es.lt.(-0.5).or.es.gt.0.5) then
g13 = 0.67
elseif (es.lt.(-0.4).or.es.gt.0.4) then
g13 = 0.47

125

elseif (es.lt.(-0.3).or.es.gt.0.3) then
g13 = 0.33
elseif (es.1t.(-0.2).or.es.gt.0.2) then
g13 = 0.24
elseif (es.lt.(-0. l ).or.es.gt.0. 1) then
g13 = 0.18
else
g13 = 0.15
endif
C Define the in-plane Poisson’s ratio either for rapid loading or for equlibrium
v12 = 1.0 - 0.5"'(el/e3)
C v12 = 0.1
C
C Determine now the components of the hypo-elastic tensor given the values of the engineering
C parameters. The hypo-elastic tensor is contained in the array D(IJ)
C
cl = e3*(l. - v12) - 2."‘el*v31*v31
C
D(l,l) = e3*e3*(1. - v12)/cl
D(l.2) = e1"‘e3*v3l/cl
D(l,3) = D(l,2)
D(1.4) = 0.0
D(2,l) = D(1.2)
D(2.2) = el*(e3 - el*v31*v3l)/(cl*(1. + v12))
D(2,3) = el*(e3*v12 + el*v3l"'v31)/(cl*(l. + v12))
D(2,4) = 0.0
D(3.1) = D(l,3)
D(3.2) = D(2,3)
D(3.3) = D(2.2)
D(3.4) = 0.0
D(4.l) = 0.0
D(4.2) = 0.0
D(4,3) = 0.0
D(4,4) = g13

(3(1) = 0.0

0(2) = 0.0

0(3) = 0.0

6(4) = 0.0
C Determine now the new value of the Cauchy stress S(I). DE(I) contains the increment of deformation
C during the current step

DO 10 I=l,4

DO 10 K=l,4

8(1) = 5(1) + D(I.K)*DE(K)

10 CONTINUE

RETURN

END

126

REFERENCES

Almeida, E. S., Spilker, R. L. and Holmes, M.H., 1995, “A transversely isotropic consititu-
tive law for the solid matrix of articular cartilage”, Proceedings of the 1995
Bioengineering Conference BED, Vol. 29, pp. 161-162.

Almeida, E. S. and Spilker, R.L., 1997, “Mixed and penalty finite element models for the
nonlinear behavior of biphasic soft tissues in finite deformation: Part H - Nonlinear
examples”, Computer Methods in Biomechanics and Biomedical Engineering, Vol.
1, pp. 151-170.

Anderson, DD. and Radin, R.L., 1991, “Stress wave effects in a finite element analysis of
an impulsively loaded articular joint”, Pmc. lnstn. Mech. Engrs., Vol. 205, pp. 27-
34.

Armstrong, C. G., Mow, V.C., and Wirth, CR, 1985, “Biomechanics of impact induced
microdamage to articular cartilage - a possible genesis for chondromalacia
patella”, In AAOS Symposium on Sports Medicine: The Knee (Edited by Finerman,
G.A.M.), pp. 70-84, C.V.Mosby Co., St Louis.

Armstrong, C. G., 1986, “An analysis of the stresses in a thin layer of articular cartilage in
a synovial joint”, Engng. Med, Vol. 15, pp. 55-61.

Armstrong, C. G., Lai, W. M. and Mow, V. C., 1984, “An analysis of the unconfined com-
pression of articular cartilage”, J. Biomech. Engng., Vol. 106, pp. 165-173.

Askew, M.J., Mow, V. C., 1978, “The biomechanical function of the collagen fibril ultra-
structure of articular cartilage”, J. Biomech. Engng., Vol. 100, pp. 105-115.

Ateshian, G. A., Lai, W. M., Zhu, W. B. and Mow, V. C., 1994 “An asymptotic solution for

127

the contact of two biphasic cartilage layers”, J. Biomechanics. Vol. 27, pp. 1347-
1360.

Ateshian, G.A., Warden, W.H., Kim, M., Grelsamer, RP. and Mow, V.C., 1997, “Finite
deformation biphasic material properties of bovine articular cartilage from con-

fined compression experiments”, J. Biomechanics, Vol. 30, pp. 1157-1164.

Athanasiou, K.A., Agarwal, A., and Constantinides, G., (1992) “Creep indentation of
human hip articular cartilage: comparison of biphasic finite element/optimization
and semi-analytical solutions”, Advances in Bioengineering, BED-Vol. 22, pp.
195-198.

Athanasiou, K.A., Rosenwasser, M.P., Buckwalter, J.A., Molinin, T.I. and Mow, V.C.,
1991, “Interspecies comparison of in situ intrinsic mechanical properties of distal

femoral cartilage”, Journal of Orthopaedic Research, Vol. 9, pp. 330-340.

Atkinson, T. S., Haut, R. C. and Altiero. N. J., 1997, “An investigation of failure criteria
for impact induced fissuring of articular cartilage”, Symposium Proceedings,

Injury Prevention Through Biomechanics, Wayne State University.

Atkinson, PJ. and Haut, RC, 1995, “Subfracture insult to the human cadaver patellofem-
oral joint produces occult injury”, Journal of Orthopaedic Research, Vol. 13, pp.
936-944.

Atkinson, P.J., Garcia, J.J., Altiero, NJ. and Haut, RC, 1997, “The influence of impact
interface on human knee injury: implications for instrument panel design and the

lower extremity injury criterion”, 41st Stapp Car Crash Conference Proceedings,
Orlando, Florida.

Bernstein, B., 1960, “ Hypo-elasticity and elasticity”, Archives of Rational Mechanics and
Analysis, Vol. 6, pp. 90-104.

128

Biot, M. A., 1941, “General theory of three-dimensional consolidation”, Journal of
Applied Physics, Vol. 12. pp. 155-164.

Brittberg, M., Lindahl, A., Nilsson, A., Ohlsson C., Isaksson 0., Peterson, L., 1994,
“Treatment of deep cartilage defects in the knee with autologous chondrocyte
transplantation”, New EnglandJ. Med, Vol. 331, pp. 889-895.

Brown, TD. and Vrahas, M.S., 1984, “The apparent elastic modulus of the juxtarticular
subchondral bone of the femoral head, J. Orthop. Res, Vol. 2, pp. 32-38.

Christensen, R. M., 1979, Mechanics of Composite Materials, Krieger Publishing Com-
pany, Malabar, Florida.

Clark, J ., 1990, “The organization of collagen fibrils in the superficial zones of articular
cartilage”, J. Anat., Vol. 171, pp.1 17-130.

Cohen, B., Gardner, T. R. and Ateshian, G. A., 1993, “The influence of transverse isotropy

on cartilage indentation behavior. A study on the human humeral head”, Trans.
Orthop. Res. Soc., Vol. 18, pp. 185.

Davis, M., Ettinger, W., Neuhaus, J ., Cho, S., Hauck, W., 1989, “The association of knee
injury and obesity with unilateral and bilateral osteoarthritis of the knee”, Ameri-
can J. Epidem.,Vol. 30(2), pp. 278-288.

Donahue, J .M., Buss, D., Oegema, TR, and Thompson, RC, 1983, “The effects of indi-
rect blunt trauma on adult canine articular cartilage”, The Journal of Bone and
Joint Surgery, Vol. 65(A), pp. 948-957.

Donzelli, P. S. and Spilker, R. L., 1996, “A finite element investigation of solid phase

transverse isotropy in contacting biphasic cartilage layers”, Advances in Bioengi-
neering, ASME, Vol. 33, pp. 349-350.

129

Eberhardt, A. W., Keer, L.M., Lewis, J .L. and Vithoontien, V., 1990, “An analytical model
of joint contact”, J. Biomech. Engng., Vol. 112, pp. 407-413.

Eberhardt, A. W., Lewis J. L. and Keer L. M., 1991, “Normal contact of elastic spheres
with two elastic layers as a model of joint articulation”, Journal of Biomechanical
Engineering, Vol. 113, pp. 410-417.

Garcia, J. J., Newberry, W. N., Atkinson, P. J ., Haut, R. C. and Altiero, N. J ., 1997, “Com-
parison of stresses in human and rabbit kness as they pertain to injuries observed
during impact experiments”, Proceedings of the 1997 Bioengineering Conference,
ASME, Vol. 35, pp. 315-316.

Garcia, J .J ., Altiero, NJ. and Haut, R.C. (In Press) “An approach for the stress analysis of
transversely isotropic biphasic cartilage under impact load”, Journal of Biome-

chanical Engineering.

Hamerrnan, D., 1989, “The biology of osteoarthrosis”, New England Journal of Medicine,
Vol. 320(20). PP. 1323-1330.

Haut, R.C., Ide, T.M., DeCamp, CE, 1995, “Mechanical responses of the rabbit patello-
femoral joint to blunt impact”, Journal of Biomechanical Engineering, Vol. 117,
pp. 402-408.

Hayes, W. C., Keer, L.M., Herrmann, G. and Mockros, LE (1972) “A mathematical anal-
ysis for indentation tests of articular cartilage”, Journal of Biomechanics, Vol. 5,
pp. 541-555.

Holmes, M.H. and Mow, V.C., 1990, “The nonlinear characteristics of soft gels and

hydrated connective tissues in ultrafiltration”, J. Biomechanics, Vol. 23, pp. 1145-
1156.

130

Holmes, M.H., 1986, “Finite deformation of soft tissue: analysis of a mixture model in

uni-axial compression”, Journal of Biomechanical Engineering, Vol. 108, pp. 372-
38 1 .

Jurvelin, J., Saamanen, A., Arokoski, J., Helminen, H.J., Kirivanta, I. and Tammi, M.,
1988, “Biomechanical properties of the canine knee articular cartilage as related to
matrix proteoglycans and collagen”, Engineering in Medicine, Vol. 17, pp. 157-
162.

Jurvelin, J.S., Buschmann, MD. and Hunziker, EB. (1997) “Optical and mechanical
determination of Poisson’s ratio of adult bovine humeral cartilage”, Journal of
Biomechanics, Vol. 30, pp. 235-241.

Khalsa, RS. and Eisenberg, SR, 1997, “Compressive behavior of articular cartilage is
not completely explained by proteoglycan osmotic pressure”, Journal of Biome-
chanics, Vol. 30, pp. 589-594.

Kelgren, J ., Lawrence, J ., 1958, “Osteoarthritis and disk degeneration in the urban popula-
tion”, Annals Rheum. Diseases, Vol. 17, pp. 386-396.

Kelly, P. A. and O’Connor, J. J ., 1996, “Transmission of rapidly applied loads through
articular cartilage. Part 1: uncracked cartilage”, Proc. Instn. Mech. Engrs., Vol.
210, pp- 27'37.

Kelkar, R. and Ateshian, G.A., 1995, “Contact creep response between a rigid imperme-
able cylinder and a biphasic cartilage layer using integral transforms”, Proceedings
of the 1995 Bioengineering Conference, ASME, Vol. 29, pp. 313-314.

Kwan, M.K, Lai, WM. and Mow, V.C., 1990,”A finite deformation theory for cartilage

and other soft hydrated connective tissues-I. Equilibrium results”, J. Biomechan-
ics, Vol. 23. PP. 145-155.

131

Li, X., 1994, A criterion to predict damage in articular cartilage due to blunt impact,
Ph.D. Dissertation, Department of Material Science and Mechanics, Michigan
State Universily.

Li, X., Haut, RC. and Altiero, N. J., 1995, “An analytical model to study blunt impact
response of the rabbit P—F joint”, Journal of Biomechanical Engineering, Vol. 117,
pp. 485-491.

Maroudas, A., 1975, “Fluid transport in cartilage”, Ann. Rheum. Dis, Vol. 34, pp. 77.

McCutchen, C. W., 1978, “Lubrication of joints”, The Joints and Synovial Fluid (Edited
by Sokoloff, L.), Vol. 1. pp. 437-483.

McCutchen, CW, 1982, “Cartilage is poroelastic, not viscoelastic (including an exact the-
orem about strain energy and viscous loss, and an order of magnitude relation for

equilibration time)”, Journal of Biomechanics, Vol. 15, pp. 325-327.

Mow, V. C. 8nd Lai, W. M., 1980, “Recent developments in synovial joint biomechanics”,
SIAM Rev.. Vol. 22, pp. 275-317.

Mow, V. C., Kuei, S. C. Lai, W. M. and Armstrong, C. G., 1980, “Biphasic creep and
stress relaxation of articular cartilage: Theory and experiment”, J. Biomech.
Engng., Vol. 102. pp. 73-84.

Mow, V.C., Holmes, M.H., Lai, WM ., 1984, “Fluid transport and mechanical properties of
articular cartilage: a review”, J. Biomechanics, Vol. 17, pp. 377-394.

Mow, V.C., Gibbs, M.C., Lai, W.M., Zhu W.B. and Athanasiou K.A., 1989, “Biphasic

indentation of articular cartilage-II. A numerical algorithm and an experimental

study”, Journal of Biomechanics, Vol. 22, pp. 853-861.

132

Mow, V.C., Zhu, W. and Ratcliffe, A., 1991, “Structure and function of articular cartilage

and meniscus”. In Basic Orthopaedic Biomechanics (Edited by Mow V.C. and
Hayes W.C.), Raven Press, New York.

Newberry, W.N. and Haut RC, 1996, “The effects of subfracture impact loading on the
patellofemoral joint in a rabbit model”, 40th Stapp Car Crash Conference Proc-
eedings, pp. 149-159.

Newberry, W.N., Zukosky, D. and Haut, RC, 1997, “Subfracture insult to a knee joint
causes alterations in bone and in the functional stiffness of overlying cartilage”,
Journal of Orthopaedic Research, Vol. 15, pp. 450-455.

Newberry, W.N., Mackenzie, CD. and Haut, RC, 1998, “Blunt impact causes changes in
bone and cartilage in regularly exercised animal model”, Journal of Orthopaedic
Research, Vol. 16, pp. 348-354.

Noll, W., 1955, “On the continuity of the solid and fluid states”, Journal of Rational
Mechanics and Analysis, Vol. 4, pp. 3-81.

Sheidegger, A. E., 1972, T he physics of flow through porous media, Third edition. Univer-

sity of Toronto Press, Toronto.

Silyn-Roberts, H. and Broom, N ., 1990, “Fracture behavior of cartilage-on-bone in
response to repeated impact loading”, Connective Tissue Research, Vol. 24, pp.
143-156.

Soltz, MA. and Ateshian, G.A., 1997, “Experimental measurement of cartilage interstitial
fluid pressurization under confined compression stress-relaxation”, Advances in

Bioengineering, BED, Vol. 36, pp. 159-160.

Spilker R.L., Suh J. and Mow V.C., 1990, “Effects of friction on the confined compressive

133

response of articular cartilage: a finite element analysis”, Journal of Biomechani-

cal Engineering, Vol. 112, pp. 138-146.

Suh J. and Spilker R.L., 1994, “Indentation analysis of biphasic articular cartilage: nonlin-

ear phenomena under finite deformation”, Journal of Biomechanical Engineering,
Vol. 116. PP. 1-9.

Terzaghi, K., 1943, Theoretical Soil Mechanics, John Wiley and Sons, Inc., New York.

Thompson, R.C., Oegema, T.R., Lewis, J.L., and Wallace, L., 1991, “Osteoarthritic
changes after acute transarticular load”, The Journal of Bone and Joint Surgery,
Vol. 73-A(7), pp. 990-1001.

Truesdell C., 1955, “Hypo-elasticity”, Journal of Rational Mechanics and Analysis, Vol.
4. pp. 83-133. 1019-1020.

Wu J. 2., Herzog W. and Ronsky J ., 1996, “Modeling axi-symmetrical joint contact with
biphasic cartilage layers-an asymptotic solution”, J. Biomechanics, Vol. 29, pp.
1263-1281.

Whipple R.R., Gibbs M.C., Lai W.M., Mow V.C., Mak AF. and Wirth CR. (1985)
“Biphasic properties of repaired cartilage at the articular surface”, Transactions of

the Orthopaedic Research Society, Vol. 10.

Xiao H., Bruhns O. T. and Meyers A., 1997, “Hypo-elasticity model based upon the loga-
rithmic stress rate”, Journal of Elasticity, Vol. 47, pp. 51-68.

Zhu W.B., Lai W.M., Mow V.C., 1986, “Intrinsic quasilinear viscoelastic behavior of the

extracellular matrix of cartilage”, Transactions of the Orthopaedic Research Soci-
ety, Vol. 11, pp. 407.

134

"illustrator