|UFDC Home||myUFDC Home | Help|
This item has the following downloads:
DYNAMIC SIMULATION OF KNEE JOINT CONTACT
DURING HUMAN MOVEMENT
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
This document is dedicated to my husband, Liu Jian.
I would like to express my gratitude to my advisor, Dr. Benjamin J. Fregly, for his
patient guidance, constant encouragement, and endless support throughout my doctoral
My thanks also go to Dr. Nam Ho Kim, Dr. Andrew J. Rapoff, Dr. Bhavani V.
Sankar and Dr. W. Gregory Sawyer, for serving on my committee. I thank them for
providing valuable comments and help in the progress of my research. I want to thank Dr.
Scott A. Banks (from The Biomotion Foundation of West Palm Beach, Florida) for
helpful fluoroscopy data and comments on the manuscripts.
I thank all of the students from the Computational Biomechanics Laboratory
(especially Jeffery A. Reinbolt and Yi-Chung Lin) for their help with contact
visualization and natural knee tests. Thanks also go to Matthew A. Hamilton from the
Triobology Laboratory for his help with wear analysis.
Last, I would like to thank my husband for his understanding and love during the
past few years. His skill and experience in programming was a great help with my model
implementation. My parents receive my deepest gratitude and love for their many years
of support and encouragement.
TABLE OF CONTENTS
A C K N O W L E D G M E N T S ................................................................................................. iii
LIST OF TABLES ....................................................... ............ ....... ....... vi
L IST O F F IG U R E S .... ...... ................................................ .. .. ..... .............. vii
A B S T R A C T .......................................... .................................................. v iii
1 IN T R O D U C T IO N ................. ................................ ........ ........ .. ............. .
1.1 B background ......... ...... .................................................................. ....................
1.2 O objective ............... .................................................................... ................... 2
1.3 D issertation O organization ............................................................................. 3
2 METHODOLOGY FOR DYNAMIC CONTACT SIMULATION.........................
2 .1 In tro d u ctio n .................... .................... .. .. ................ ................ .. 5
2.2 M ultibody K nee Contact M ethodology ........................................ .....................7
2.2.1 C contact Surface Preparation ........................................ .....................8
2.2.2 Efficient Distance Calculations.......................................... 11
2.2.3 Contact M odel Development................................ ...... ............... 17
2.2.4 Dynamic Model Construction..... ................................20
2.3 Sam ple Applications ......... .. .............. .... ..........................24
2.3.1 Static Analysis of Natural Knee Contact........................................... 24
2.3.2 Dynamic Simulation of Artificial Knee Contact .............. .....................27
2.4 Discussion .................................... ................................ .........30
3 EXPERIMENTAL EVALUATION OF AN ELASTIC FOUNDATION MODEL
TO PREDICT CONTACT PRESSURES IN KNEE REPLACEMENTS ..............35
3.1 Introduction..................................................................... .... .. .......... 35
3.2 M materials and M ethods............................................... ............................. 37
3.2.1 E plastic C contact M odel ...................................... .................. .... ........... 37
3.2.2 Param etric M material M odel .......................................... ............... 39
3.2.3 D ynam ic Im plant M odel .................................. ..................................... 40
3.2.4 Contact Pressure Experim ents ............................................ ............43
3 .3 R esu lts ...................................... ............................ .... ........ ...... 4 7
3 .4 D iscu ssio n .................. ...................................... ......... ................ 4 9
3.4.1 Contact M odel A accuracy ............................................. ............... 49
3.4.2 Contact Model Advantages .......................................... ...............50
3.4.3 Contact M odel Lim stations .......................................... ............... 52
4 PREDICTED SENSITIVITY OF KNEE IMPLANT WEAR TO INSERT
THICKNESS AND BODY M ASS.................................... ......................... 54
4.1 Introduction..................................................................... .... .. .......... 54
4.2 M materials and M ethods............................................... ............................. 56
4.2.1 N om inal W ear Prediction ............................................ ...............56
4.2.2 Neighboring W ear Predictions..... .......... .......................................59
4 .3 R results ......... ......... .............................................. ............................60
4.4 Discussion ............. ...... ..... ..... .......... ..... .............. .........63
5 CONCLUSIONS AND FUTURE RESEARCH ............................................... 69
A MAXIMUM PRESSURES DURING SIMULATIONS .......................................71
B WEAR, CREEP AND DAMAGE DEPTH .......................................................72
C DAMAGE AREA AND VOLUME LOST ................................... .................73
L IST O F R E FE R E N C E S ............. ......... ............ ................................................74
B IO G R A PH IC A L SK E TCH ..................................................................... ..................83
LIST OF TABLES
2-1. CPU tim e tests of distance evaluation m ethods...................................................... 17
2-2. Input parameters and output measures in the API between the contact model and
dynam ic contact m odel .......................................................................... 21
3-1. Comparison of linear and nonlinear material model ......................... ...............47
4-1. Linear regressions to wear predictions. ....................................... ............... 62
A-1. Maximum and average pressures during a cycle of gait simulation.....................71
B-1. Wear, creep and damage depth in the TKR insert. .............................................72
C-1. Damage area and volume lost in the TKR insert. .................................................73
LIST OF FIGURES
2-1. Fram work of contact model .............. ......... ............................ ..................8
2-2. Distance correction for two distance-finding methods ......................... .........13
2-3. Sensitivity of contact results to grid density.................................. ..................... 23
2-4. Visualization of example dynamic simulation for settling the femur onto the tibia.26
2-5. Contact area of natural knee at the static position .............................................. 27
2-6. Dynamic simulation of contact during a cycle of gait for Osteonics7000 implant
design ..............................................................................2 9
2-7. D ynam ic sim ulation results. ............................................. ............................ 29
2-8. Contact areas during the contact phase of a gait cycle..................... ..............30
3-1. Overview of materials and methods. .................. .. ....................... ............ .. 38
3-2. Visualization of example dynamic simulation used to settle the femoral component
onto the tibial insert .................. ...................................... .. ............ 42
3-3. Experimental setup and pressure measurement......................................................44
3-4. Comparison between experimental and predictions..........................................48
4-1. Overview of the computational framework for wear predictions.............................57
4-2. Comparison of retrieval and predictions. ............. ............ ................................59
4-3. The damage plots for five combinations of insert thickness and body mass. ..........61
4-4. Three-dimensional surface plots of the total damage depth as a function of insert
thickness and body m ass. ............................................... .............................. 62
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
DYNAMIC SIMULATION OF KNEE JOINT CONTACT DURING HUMAN
Chair: Benjamin J. Fregly
Major Department: Mechanical and Aerospace Engineering
The knee is one of the most important joints in the human body. Degenerative
damage in this joint sometimes results in function loss, and leads to total knee
replacements (TKR). However, mild wear of ultra-high-molecular-weight-polyethylene
(UHMWPE) tibial inserts greatly affects the longevity of TKRs. To understand the
dynamics of the natural knee and to improve TKR implant designs, it is essential to
develop proper tools to study the contact and wear mechanism of the knee. This
dissertation provides the conceptual and computational details of a methodology for
investigating contact and wear in the knee during human movements. It includes four
steps: articular geometry preparation, efficient surface-surface distance evaluation, three-
dimensional contact model development, and dynamic contact model construction. The
geometry of the articular surfaces is obtained from CT and MRI images for the natural
knee or from CAD models for the implant designs. The contact model is incorporated
into the dynamic simulation system. The dynamic simulation is driven by in vivo
fluoroscopy data of gait or stair. Wear is predicted by a computational wear model using
the dynamic contact solutions. Sample analyses compare well to experiment results and
TKR insert retrievals with reliable accuracy within reasonable CPU time.
This methodology is applied to the study of wear sensitivity of TKR polyethylene
to insert thickness and patient body mass. The simulations of twenty five combinations of
insert thickness (6, 8, 10, 12 and 14 mm) and body mass (50, 75, 100, 125 and 150 kg)
are performed in the neighborhood of a nominal simulation that predicts in vivo damage
well both quantitatively and qualitatively. Each simulation predicts maximum wear,
creep, and damage depth, along with damage area and volume lost. When the
polyethylene thickness increases, maximum wear depth, creep depth, damage depth, and
volume lost decrease while wear area increases. The regression equations are fit to the
results and can be used to estimate the wear benefit achieved by using a thicker insert
given the patient's body mass or by losing weight, given the insert thickness.
The knee is one of the most important and most studied joints in the human body.
Unfortunately, it is highly susceptible to osteoarthritis which causes pain and loss of
function and in some cases leads to total knee replacement (TKR). Contact mechanics
influences osteoarthritis progression in natural joints (Harman et al 1998, Riegger-Krugh
et al. 1998) and wear in artificial joints (Schmalzried et al. 2000). Deformable joint
models are desirable for studying knee contact mechanics under in vivo dynamic loading
conditions. Moreover, deformable knee joint models should be integrated into full-limb
or full-body human movement models for three main reasons. First, human movement is
dynamic, indicating that the loading experienced by the knee joint will be a function of
the motions occurring during the activity. Second, dynamic systems are coupled; and
movement at one joint cannot be studied apart from movement of the entire system
(Zajac 1993). Third, muscles and ligaments play an important part in determining joint
contact forces and motion (Hsieh and Draganich 1998, Li et al. 1999). These observations
suggest that models combining skeletal dynamics (rigid multibody dynamics), muscle
actuators, ligament bundles, and deformable joint contact are needed to predict knee
contact conditions during movement.
As shown above, motion and contact affect dynamics in natural knees and TKRs.
However, contact conditions cannot currently be measured in vivo. Thus, dynamic
simulations of knee mechanics would be valuable for predicting in vivo contact pressures
and analyzing wear. To develop such simulations, deformable knee models that can
accurately predict contact pressures need to be integrated into a framework that combines
multibody dynamics with musculoskeletal modeling.
Significant rigid body computational tools already exist for modeling and
simulating of human movement (Davoodi and Loeb 2001, Davoodi et al. 2003, Delp et
al. 1995). These tools can be used to build subject-specific models. The key problem that
remains is how to develop subject-specific deformable joint contact models and
incorporate them into larger human movement models.
Finite element analysis (FEA) has been used to study knee joint contact mechanics
under static loading conditions (Bartel et al. 1986 and 1995, Bendjaballah et al. 1997,
D'Lima et al. 2001, Perie and Hobatho 1998, Otto et al. 2001, Rawlinson and Bartel
2002, Sathasivam and Walker 1998 and 1999). Dynamic FEA has recently been applied
to simulations of knee implant components under well-defined loading conditions
(Giddings et al. 2001,Godest et al. 2002, Halloran et al. 2003). One drawback of these
analyses is their intensive use of CPU time, typically requiring hours or days to predict
motion and loads simultaneously. This makes them impractical for incorporation into a
larger dynamic musculoskeletal model to be used in design sensitivity or optimization
studies. Furthermore, detailed stress analyses carried out by FEA are unimportant in gross
movement simulations. Compared to FEA, an elastic foundation contact model
(Blankevoort et al. 1991, Li et al. 1997, Nufio and Ahmed 2001, Pandy et al. 1997) is a
simpler and more efficient way to predict the knee joint contact pressures.
Combining deformable joint models with rigid body system is one of the important
future research directions in multibody dynamics (Schiehlen 1997). However, no such
computational framework exists to date, primarily because of the lack of a fast and
accurate contact solver as well as the technical difficulties involved in implementing such
a solver within a multibody dynamics environment. Our study aimed to develop a
methodology to study contact in the knee joint during human movement, thus enabling us
to study osteoarthritis development in natural knees and wear development in ultra-high-
material-weight-polyethylene (UHMWPE) of artificial knees.
1.3 Dissertation Organization
Chapter 2 gives conceptual and computational details are presented for a
methodology to combine deformable contact modeling with multibody dynamic
simulation for specific application to the knee. Four steps are involved in the
methodology. Subject-specific geometry of articular knee surfaces is obtained from CT
and MRI images for natural knees or from CAD models for artificial knees. An efficient
distance evaluation method is developed. The deformable contact model is developed
based on elastic foundation theory. Dynamic contact simulations system are created by
incorporating the contact model into commercial multibody dynamic software. Two
sample applications of the methodology to natural and artificial knees are shown and
discussed. One dynamic simulation is used to demonstrate a static analysis of a natural
knee. A second dynamic simulation is used to demonstrate a dynamic analysis of an
In Chapter 3, the methodology is evaluated by performing static pressure
experiments with a commercial knee implant in neutral alignment using flexion angles of
0, 30, 60, and 90 and loads of 750, 1500, 2250, and 3000 N. Using manufacturer CAD
geometry for the same implant, an elastic foundation model with linear or nonlinear
polyethylene material properties was implemented within a commercial multibody
dynamics software program. The model's ability to predict experimental peak and
average contact pressures simultaneously is evaluated by performing dynamic
simulations to find the static configuration.
In Chapter 4, the methodology is applied to study the wear sensitivity of a TKR to
insert thickness and body mass. Although a thicker insert may minimize wear in heavier
patients, the remaining bone stock is also reduced. Thus, estimation of wear as a function
of insert thickness and patient body mass could be helpful for surgical planning. Such
estimates are developed by performing simulations in the neighborhood of a nominal
simulation that predicts in vivo damage well, both quantitatively and qualitatively.
Twenty five combinations of insert thickness (6, 8, 10, 12 and 14 mm) and body mass
(50, 75, 100, 125 and 150 kg) were studied in gait simulations. Output from each
simulation was used in a computational wear model to predict maximum wear, creep, and
damage depth, along with damage area and volume lost. All predicted quantities were fit
using multiple linear regression. The regression results can be used to estimate the wear
benefit achieved by using a thicker insert given the patient's body mass, or alternatively,
by losing weight given the insert thickness. Furthermore, we studied the correlation
between maximum pressures and wear quantities predicted by the simulations to
determine the extent to which pressure predictions alone can provide information about
Chapter 5 summarizes the research with recommendations for future improvements
to the contact model and simulation efficiency.
METHODOLOGY FOR DYNAMIC CONTACT SIMULATION
According to the Arthritis Foundation, arthritis-related problems are second only to
heart disease as the leading cause of work disability. Mechanical loading, and especially
dynamic loading, is believed to play a major role in the degenerative process (Setton et
al. 1998, Tetsworth and Paley 1994). Furthermore, motion (i.e., kinematics) and loading
(i.e., contact pressures) appear to have an interactive effect on disease progression
(Andriacchi 1994; Tashman and Anderst 2003). A similar situation exists in artificial
joints, where motion and loading both influence wear progression (Archard and Hirst
1956) which can lead to osteoylsis, and ultimately to implant failure. Even for tissue
engineering applications, the mechanical environment of the replacement or repair cells
must be understood to optimize the repair process (Guilak and Goldstein 2001).
Knowledge of in vivo joint motion and loading during functional activities would
therefore help address these clinically significant issues.
While dynamic x-ray imaging advances now permit accurate measurement of in
vivo joint motion (Banks and Hodge 1996, Hoff et al. 1998, Kanisawa et al. 2003,
Komistek et al. 2003, Tashman and Anderst 2003), a non-invasive experimental approach
does not exist for measuring in vivo joint loading. Consequently, computer simulations
are needed to develop predictions. Since muscle co-contractions increase joint
compressive forces, multibody dynamic musculoskeletal models are needed to provide
the estimated muscle forces (Anderson and Pandy 1999, Davoodi et al. 2003, Delp et al.
1998, Neptune et al. 2000). However, the types of joints used in these models have a
dramatic effect on predicted forces (Glitsch and Baumann 1997). Furthermore, three-
dimensional anatomical models produce muscle force estimates that are more consistent
with experimental EMG data (e.g., antagonistic muscle activity) than estimates obtained
from two-dimensional models (Li et al. 1999). These observations suggest that for some
joints, applied muscle forces and joint contact loads should be computed simultaneously
using joint models that include the articular surface geometry (Hefzy and Cooke 1996, Li
et al. 1999). For statically indeterminate joints such as the knee (Cheng et al. 1990), these
joint models require a deformable rather than a rigid body contact theory to produce
accurate contact force estimates (Piazza and Delp 1999).
For the specific case of the knee, numerous studies have highlighted the value of
combining multibody simulation methods with three-dimensional joint contact modeling.
Wismans et al. (1980) were the first to develop a quasi-static three-dimensional natural
knee model using rigid body contact. Pandy et al. (1997) presented a quasi-static natural
knee model using deformable contact with idealized (i.e., planes and polynomials)
surfaces. Abdel-Rahman and Hefzy (1998) used similar idealized surfaces with rigid
body contact to create a dynamic model of a natural knee. Piazza and Delp (1999)
extended this work by applying rigid body contact to an artificial knee within a full-body
dynamic simulation. Blankevoort et al. (1991), Cohen et al. (2001 and 2003), Kwak et al.
(2003), and Elias et al. (2003) reported quasi-static natural knee models with deformable
contact, where subject-specific contact surfaces were generated from medical imaging
data. Fregly et al. (2003 and 2004) developed dynamic models of artificial knees using
deformable contact, with static and dynamic analyses performed via dynamic simulation.
As indicated by the brief review above, multibody models incorporating three-
dimensional knee contact have progressed from quasi-static to dynamic, from rigid body
contact to deformable contact, and from simple surface geometry to complex surface
geometry. However, few studies have merged multibody dynamic simulation methods
with deformable knee contact models, largely because of the computational cost of
repeated geometry evaluations during numerical integration.
Our study presents a computationally efficient methodology for incorporating
deformable, three-dimensional knee contact models into a multibody dynamic simulation
environment. The methodology combines a variety of geometric modeling, theoretical,
and numerical approaches to create a deformable knee contact model that can be
incorporated into any multibody dynamic simulation software used for musculoskeletal
modeling. Though the current model is tailored to the natural and artificial tibiofemoral
joint, similar methodology can be applied to other joints as well. Sample applications are
presented for a natural knee contact model created from MRI and CT data; and for an
artificial knee contact model produced from manufacturer's CAD data. By giving the
details of the computational methodology, we hope to make this approach more widely
available for use in dynamic musculoskeletal models that seek to predict in vivo joint
motion and loading.
2.2 Multibody Knee Contact Methodology
We propose a computationally efficient four-step methodology for incorporating
deformable contact models of human joints into a multibody dynamics. The four steps
involve preparing of the articular surface geometry, developing of efficient methods for
repeated evaluation of the contact geometry, developing of an efficient contact solver that
accounts for the unique characteristics of human joints, and specifying of an application
programming interface (API) and appropriate numerical integration methods that will
work within any multibody dynamic simulation environment framework (Figure. 2-1).
The result is a deformable contact modeling methodology that can be easily incorporated
into larger multibody dynamic models of the musculoskeletal system containing elements
such as muscles, ligaments, neural controllers, and ground contact models. Although the
development here is for deformable contact models of the human tibiofemoral joint, a
similar methodology can be applied equally well to other joints.
Dynamic Simulation System Elastic Contact Model Articular Surfaces
t I Femur
\[ L' J Tibia
Geometry Evaluatin Numerical Routines Contact Solver
Figure 2-1. Framework of contact model. Geometry of processed articulate surfaces is
input to the contact model. The contact model has three functions: geometry
evaluation, numerical calculation and contact solver. Contact model is
integrated into the dynamic contact simulation system.
2.2.1 Contact Surface Preparation
Repeated distance evaluations between contacting surfaces are the primary
computational bottleneck when incorporating elastic contact into a multibody framework.
For any elastic contact problem solved via numerical methods, distances between the
undeformed contact surfaces must be sampled for each relative pose of the contacting
bodies specified by the numerical integrator, regardless of the model used for calculating
contact pressures. From our experience, this typically requires about 80% of the total
CPU time during the course of a dynamic simulation, and possibly more depending on
the efficiency of the contact solver. Thus, efficient methods for surface-surface distance
evaluation are essential.
One factor that can significantly affect distance calculation efficiency is the manner
in which the contacting surfaces are modeled. Contact surfaces are commonly
represented mathematically, using either polygonal approximations (McGuan et al. 1998,
Piazza and Delp 1999) or nonuniform rational B-spline (NURBS) patches (Ateshian
1993; Farin 1997). The advantage of polygonal (or tessellated, surface geometry is that
distance calculations can be performed very efficiently. However, discontinuities in the
surface geometry can lead to discontinuities in contact forces and torques (Puse and
Laursen 2002), which can dramatically slow numerical integration. Thus, we chose to
work with NURBS patches, to represent the contact surface geometry as accurately as
The efficiency with which distances between NURBS patches can be calculated
depends on at least four factors: the number of spline curves used to represent each spline
patch, the number of spline patches used, whether or not the spline patches are trimmed
(i.e., part of the surface is considered to be removed), and the method used for calculating
distances. In computer aided design (CAD) software, each three-dimensional NURBS
patch is described mathematically by a rectangular, two-dimensional (i.e., u, v)
parametric space. The speed with which a NURBS patch can be evaluated is proportional
to the number of spline curves used to represent it in the u and v directions. Thus, one
goal is to create smooth and accurate surface geometry using as few u-v curves as
Tradeoffs exist between use of a single NURBS patch or multiple patches to
represent a contact surface. If a single NURBS patch is used, a larger number of spline
curves will be required in the u and v directions than if multiple NURBS patches are
used. However, one disadvantage of multiple patches is that small "holes" can exist at the
seams between two patches due to geometry tolerance issues. This can produce infinite
distance calculations in some geometry modeling packages. In either case, trimmed
surfaces, where a portion of the surface is removed, take longer to evaluate than non-
trimmed surfaces due to the computation required determining if a point lies inside or
outside the trimming. The general goal therefore becomes to represent each contact
surface with as few untrimmed NURBS patches as possible, each with the minimal
number of u-v curves necessary to obtain the desired surface accuracy.
Creation of single or multiple untrimmed NURBS patches can be achieved using a
similar process for natural and artificial knee contact geometry. For both types of knees,
the general process involves finding points on the contact surfaces, converting these
points to a polygonal surface representation, and fitting NURBS patches to the polygonal
surfaces. For natural knees, points defining the articular cartilage and subchondral bone
surface geometry can be obtained by segmenting MRI data using commercial image
processing software. For artificial knees, points defining the tibial insert and femoral
component contact surfaces can be obtained by laser scanning. The point clouds can then
be converted to polygonal representations using commercial reverse engineering software
such as Geomagic Studio (Raindrop Geomagic, Research Triangle Park, NC). For
artificial knees, polygonal surface geometry can also be obtained directly from
manufacturer CAD models typically containing trimmed surfaces. Once polygonal
geometry is obtained, reverse engineering software can be used to create a rectangular
grid of multiple NURBS patches. These patches can be either used directly or, since they
are untrimmed, merged to create a single NURBS patch for each contact surface. Surface
modeling software such as Rhinoceros (Robert McNeel & Associates, Seattle, WA) can
resample the u-v space and reduce the number of u-v lines for each surface patch. Fewer
lines are typically required in one direction than the other. Deviations between the
original point cloud data (or CAD geometry) and final NURBS surfaces can be quantified
in Geomagic Studio, and the number of u-v lines modified until the desired surface
accuracy is achieved.
As a final surface preparation step, all geometry is transformed to facilitate
subsequent geometry calculations. For both knee types, the positive axis is made to be
superior, the positive z axis medial, and the positive x axis defined by y cross z (posterior
for a right knee, anterior for a left knee). For a natural knee, an anatomic coordinate
system following this convention can be constructed for each bone such using estimated
joint centers and anatomic landmarks (Chao and Sim 1995). For an artificial knee,
implant dimensions can be used to define the coordinate system origin in the x and z
directions. For they direction, the bottom surface of the tibial insert and the dimensions
of the bone box of the femoral component can be used to define the origin location.
2.2.2 Efficient Distance Calculations
The two primary goals for distance calculations are to minimize the number of
calculations required and to make each calculation as computationally efficient as
possible. An additional requirement is development of a distance sampling method that
works for conformal (i.e., one surface with positive curvature and the other with negative
curvature) and nonconformal (i.e., both surfaces with positive curvature) situations
equally well. Both situations exist in human joints and may even exist at different
locations within the same joint (e.g., a patellofemoral joint model).
Given the relative pose of the two contacting bodies as specified by the numerical
integrator at any instant, four strategies are used to minimize the number of distance
calculations. One strategy is to perform all geometry initialization in a pre-processing
phase prior to beginning the dynamic simulation. In this way, many CPU-intensive
operations are performed only once per simulation. These include reading in the contact
surface geometry for both bodies (e.g., tibia and femur), detecting connectivity between
the NURBS patches defining individual contact surfaces, and determining whether the
surfaces are from a natural knee or an artificial knee. A natural knee will have twice as
many contact surfaces as an artificial knee due to the subchondral bone surfaces needed
to calculate local cartilage thickness. Contact pairs and associated data structures are then
created, where each pair consists of a "fixed" body and "moving" body for calculating
relative kinematics. For the fixed body, an oriented bounding box placed around each
contact surface is used to define a plane aty = 0 on which a rectangular grid of elements
is created. The center points of these elements are projected in they direction onto the
corresponding contact surface to create an approximately rectangular grid of contact
elements (Paul and Hashemi 1981). Note that not all of the projected elements will have
corresponding contact elements on the fixed body contact surface. With this method, each
contact element can be assigned a two-dimensional index that facilities storage of the
element thickness, area, and local coordinate system for later use.
One additional step is required to account for cartilage thickness in a natural knee.
For each element on the fixed body, the distance between the contact surface and the
corresponding subchondral bone surface is calculated and stored as the element thickness.
For the moving body, each NURBS patch in the contact surface is subdivided based on
its u-v space, and the distance between the sampled points and the corresponding
subchondral bone surface calculated and stored in a separate data structure. During a
simulation, once a u-v location is found on a particular NURBS patch of the moving
body, bilinear interpolation is used to estimate the cartilage thickness at that location.
Tangent Plane Moving Surface
Fixed Surface C
Figure 2-2. Distance correction for two distance-finding methods.
Another strategy for minimizing the number of distance calculations is to calculate
distances from one surface directly to the other. In nonconformal numerical contact
models using elastic half-space theory, it is common to create a tangent plane (i.e., a
midsurface) between the two contact surfaces, with all contact calculations being
performed relative to this midsurface (Figure. 2-2A). The tangent plane is perpendicular
to the unique line defined by where the normals on the two contact surfaces are co-linear
and anti-parallel. There are several drawbacks to this approach. Finding the tangent plane
requires solving a nonlinear root finding problem that is sensitive to the initial guess and
fails frequently for complex geometry, halting the entire dynamic simulation. A tangent
plane does not work well for conformal joints such as the hip, since the contact
midsurface is highly curved and may even be "wavy" for some anatomic joints.
Furthermore, a single tangent plane does not work well in situations where multiple
contacts occur between the same pair of surfaces. Finally, a tangent plane necessitates
two sets of distance calculations one from the plane to each of the contact surfaces,
with the total distance between the surfaces being the sum of the two distances.
To eliminate these issues, we have developed an approach to calculate the
distances between two contact surfaces using implicit rather then explicit midsurface
creation. For any contact element on the fixed body, the distance to the corresponding
contact surface on the moving body is calculated directly and then corrected, cutting the
number of distance evaluations in half. This is done using either a minimum distance or
ray firing method starting from the fixed body. For any element on the fixed body,
minimum distance finds a distance vector that is perpendicular to the moving body, while
ray firing finds a larger magnitude distance vector that is perpendicular to the fixed body.
However, the distance vector used to calculate contact pressures should lie along the
normal of the midsurface between the two contact surfaces, regardless of whether the
contact is nonconformal (Figure. 2-2B) or conformal (Figure. 2-2C). This is achieved by
using trigonometry to make the minimum distance vector (m in Figure. 2-2) slightly
larger or the ray firing distance vector (r in Figure. 2-2) slightly smaller (Figure. 2-2C).
The midsurface normal can be well approximated by a line (ADF in Figure. 2-2C)
bisecting the minimum distance (AB in Figure. 2-2C) and ray firing (AC in Figure. 2-2C)
vectors. If minimum distance is used, the ray firing direction is found by evaluating the
surface normal on the fixed body. Once the bisecting direction is calculated, its
magnitude (AD) is defined such that its projection along the minimum distance direction
produces the minimum distance vector (AB). If ray firing is used, the minimum distance
direction (AB) can be approximated by evaluating the surface normal on the moving
body (CE) at the u-v location hit by the ray. If the distance vector from ray firing (AC) is
projected onto the minimum distance direction (AE), a corrected distance (AF) can be
found such that its projection along the minimum distance direction produces the
approximated minimum distance vector (AE). In practice, distance results from the two
approaches are nearly identical, with the corrected distances always lying between the
minimum distance and ray firing results.
A third strategy for minimizing the number of distance calculations is to copy
previous results whenever the relative kinematics change little. Numerical integrators
frequently put the contacting bodies in the same relative pose several times while
performing trial steps. If the change in kinematic variables (i.e., generalized coordinates
and their time derivatives) is below a small user-defined tolerance, pressure results stored
from the previous pass are copied to the current solution. If only the derivatives of the
generalized coordinates change, the elastic portion of the previous pressure results is
copied, and only velocity is recalculated.
A final strategy for minimizing distance calculations is to search for new contacting
elements starting from the results of the previous pass. Searching all elements for
interpenetration, as required during the first pass, is extremely expensive
computationally. Since numerical integration changes the relative pose of the contacting
bodies only incrementally as the simulation progresses, it is reasonable to expect that the
new contact area will contain some of elements from the previous solution. Thus, if the
change in relative kinematics is considered too large to copy the previous solution, the
new contact area is found from the previous one using an efficient Breadth-First search
algorithm. This algorithm starts by testing all of the previous active elements one at time.
If a tested element is still active, all of its neighbors are searched systematically, with the
search progressing outward until a ring of non-active elements is identified. Remaining
elements on the list that have not yet been evaluated are tested in a similar manner. In
addition to minimizing the number of distance evaluations, this method is able to identify
multiple contact regions that form by the splitting of a single contact region. If no active
elements are found by the Breadth-First search, all elements on the fixed body are tested
to ensure that no contact is occurring elsewhere.
To make each distance calculation as computationally efficient as possible, we
quantified the CPU time for the minimum distance and ray firing methods using single
and multiple patch NURBS surfaces from the femoral component of an artificial knee.
Distance calculations were performed using the ACIS 3D Toolkit (Spatial Corporation,
Westminster, CO), and CPU times were quantified using Rational Quantify (IBM
Corporation, White Plains, NY) (Table 2-1).
For single surface patches of comparable tolerance, minimum distance was 5.1 to
8.7 times faster than ray firing, with the fastest surface being the one with the fewest u-v
lines. For multiple surface patches of comparable tolerance, ray firing CPU times
decreased by a factor of 2.3 when going from a single patch with dense u-v lines to 96
patches with sparse u-v lines. The majority of the performance gain from multiple patches
was achieved with only 12 patches. Though a larger number of patches resulted in
slightly more accurate geometry, computational considerations indicate that each contact
surface should be merged into a single NURBS patch with a minimal number of u-v lines
Table 2-1. CPU time tests of distance evaluation methods.
Tolerance statistics (mm) CPU time (micro-sec)
Maximum Average STD Ray firing Min.distance
Single surface 4,459,921 879,488
nglesurface 0.0412 0.00448 0.0053 3,789,363 437,182
(18 x 14)
12-patch surface 0.0486 0.00148 0.0019 2,175,888 2,461,545
24-patch surface 0.0434 0.00202 0.0028 2,031,058 2,791,974
(5 x 6)
48-patchsurface 0.0397 0.00163 0.0026 1,929,121 3,442,914
(4 x 5)
96-patch surface 0.0276 0.00134 0.0020 1,908,546 3,577,557
The single surface was a merged femur articular surface from the implant CAD
model. Very dense grids (50 x 50) were used to create the single patch surface. The other
geometries were generated by rebuilding or splitting the single surface along its
isoparametric lines (lines that are parallel to U or V direction in the UV space). The
tolerance of the other geometries from the single patch was measured by three
parameters: the maximum deviation, the average deviation and the standard deviation
(STD). For single surfaces, deformations could be found by either ray-firing or minimum
distance from a point to a single surface. For multi-patch surfaces, the ray-firing and the
minimum distance from a point to a body could be used.
2.2.3 Contact Model Development
Once distances are sampled between the undeformed contact surfaces, a variety of
elastic contact theories can be used to generate contact pressure results numerically.
Elastic half-space theory, where the contacting bodies are semi-infinite and possess
homogenous linear elastic material properties, is extremely well developed in the
literature (Ahmadi et al. 1983, Conry and Seireg 1971, Johnson 1985, Kalker and Van
Randen 1972, Paul and Hashemi 1981). However, it violates many characteristics of
natural and artificial joints. Specifically, the contacting surfaces are not approximately
planar in the region of contact, an elastic layer of finite thickness covers one or both
contacting bodies, and the breadth of the contacting bodies is of comparable dimensions
to the contact area.
An alternative that satisfies these characteristics is elastic foundation (or "bed of
springs") theory (An et al. 1990, Blankevoort et al. 1991, Li et al. 1997, Pandy et al.
1997). Derived from plane strain elasticity theory for an elastic layer bonded to a rigid
substrate (Johnson 1985), this contact modeling approach scatters independent springs
over the contact surfaces, where the springs represent an elastic layer of known thickness
covering one or both bodies. Unlike elastic half-space theory, this simplified contact
model does not account for how pressure applied at one surface location produces
deformations at all locations, thereby eliminating the integral nature of contact problems.
However, the benefits of this simplification are faster pressure calculations and facilitated
analysis of conformal geometry, layered contact, and nonlinear materials.
We select different forms of the elastic foundation model depending on the
magnitude of surface deformations and the type of material model. For small
deformations such as in artificial knees (An et al. 1990, Blankevoort et al. 1991, Li et al.
1997, Pandy et al. 1997), Eq. 2-1 is used to calculate the pressure of any spring:
p = d (2-1)
(1+ v)( 2v)h
where E is Young's modulus of the elastic layer, v is the Poisson's ratio of the elastic
layer, h is the layer thickness, and d is the spring deformation. Both h and d are calculated
on an element-by-element basis, with d defined as the interpenetration of the undeformed
contact surfaces in the direction of the approximated midsurface normal. If both bodies
possess an elastic layer of the same material, then the two layers are treated as a single
layer of combined thickness. Thus, for natural knees, the local thickness of the tibial and
femoral articular cartilage are added to determine h for any element, while for an
artificial knee, h is due to the local thickness of the tibial insert assuming the femoral
component is rigid. For larger deformations as in natural knees, the stiffness of the elastic
layer increases with surface deformation because of geometric nonlinear behavior
(Blankevoort et al. 1991), so that Eq. 2-1 becomes:
(1- v)E ( (2-2)
For nonlinear materials, E can be a function of the pressure p (Nufio and Ahmed
2001), which we derive from a three-parameter nonlinear material model (Fregly et al.
E= -E +- (2-3)
2 po 2 po
where E is strain, p is contact pressure, and so, Po, and n are material parameters.
This is a nonlinear power-law material model (Johnson 1985) with the addition of a linear
term which results in a standard linear material model when n = 1. This material model
has been shown to fit experimental data for polyethylene extremely well (Cripton 1993,
Fregly et al. 2003). For an estimated value forp, E = dp/de c an be calculated from Eq.
E= 1/1 1+n (2-4)
2 po, Po
Given the interpenetration dfor any spring, Eq. 2-4 can be substituted into Eq. 2-1 or Eq.
2-2 to produce a single nonlinear equation forp that can be solved using standard
nonlinear root finding methods.
Once p is known for every element, the net contact force and torque acting on the two
bodies are computed using the principle of replacement (Kane and Levinson 1985). The
contact force on each element in the direction of the approximated midsurface normal is
found by multiplying the element pressure by the element area. Then the calculated
element forces are replaced with an equivalent contact force vector, defined as the vector
sum of the individual element forces applied at the origin of the fixed body coordinate
system, along with a contact torque vector, defined as the sum of the moments of the
individual element forces about the fixed body coordinate system origin. In this way, the
contact model looks like any other load in the model as far as the numerical integrator is
2.2.4 Dynamic Model Construction
We have developed an application programming interface (API) to allow our
contact modeling methodology to function within any multibody dynamic simulation
environment. The approach described above was coded in C++ and compiled as a
dynamic link library (DLL). All geometry evaluations are performed using the ACIS 3D
Toolkit (Spatial Corporation, Westminster, CO). In addition to model parameters (e.g.,
materials, friction, damping, number of elements) and model options (e.g., choice of
contact model, material model, distance sampling method, output quantities see Table
2-2), six generalized coordinates and their time derivatives must be passed to the DLL by
the multibody code. The generalized coordinates specify the relative pose of the moving
body with respect to the fixed body, with the rotations defined using a body-three 1-2-3
sequence (Kane et al. 1983).
Table 2-2. Input parameters and output measures in the API between the contact model
and dynamic contact model.
Input parameters Output measures
Contact model (elastic foundation or Contact force (x, y, z components and
elastic half space) magnitude) on medial and lateral sides
Material model (small or large Contact torque (x, y, z components
strain) and magnitude) on medial and lateral
Material parameters (material sides
properties, thickness offset) Maximum and average pressure on
Friction parameters (static and medial and lateral sides
dynamic friction coefficients) Contact area
Damping parameter (nonlinear) Maximum interpenetration
Distance calculation method (ray Output loads (equal and opposite for
firing or minimum distance) the two contacting bodies)
Element grid density (number of
divisions along each direction)
Output flag (off or on, and output
Outputs from the DLL to the multibody code include a net contact force and torque
applied to each body, along with measured quantities of interest (e.g., peak and average
contact pressures, contact forces, contact areas, and number of active elements see
Table 2-2). Output data files can also be generated depending the setting of an output
flag. The net contact force and torque are computed using the principle of replacement
(Kane and Levinson 1985). The calculated element forces (elastic, frictional, and
damping) are replaced with an equivalent contact force vector, defined as the vector sum
of the individual element forces applied at the origin of the fixed body coordinate system,
along with a contact torque vector, defined as the sum of the moments of the individual
element forces about the fixed body coordinate system origin. In this way, the contact
model looks like any other load in the model as far as the numerical integrator is
Three primary dynamic analyses can be performed with this framework: static
analyses, forward dynamics simulations, and inverse dynamics analyses. Static analyses
are listed in this category since they can be performed easily and efficiently via a
modified dynamic simulation method. Static analyses are useful for settling the
contacting bodies together prior to performing a dynamic simulation, and also for
predicting contact pressures, areas, and forces for comparison with experimental data
(Fregly et al. 2003). For multibody dynamic systems, static analyses are typically
performed by placing large numbers of dampers throughout the system. The equations of
motion are then numerically integrated until the maximum acceleration in the system falls
below some user-defined tolerance. The difficulty with this approach is determining the
placement of the dampers and the values of the damping coefficients. An alternate
approach, and the one employed here, is to use numerical rather than physical damping to
settle the bodies into a static configuration. This can be easily achieved by employing a
low order implicit integrator (e.g., an implicit Euler integrator or Matlab's ODE15s
integrator set to first order) with loose tolerances for the numerical integration. Though
the settling motion predicted by the integrator will be inaccurate, the final static
configuration of the system will be correct because both the acceleration and velocities
are within a small tolerance at the end of the static analysis.
1600 / j 38 -
o 1200. 30 .
20 20 20 -20
O 40 4z.40 4.4 40 40
,qo 60 80 ..... 60 0 -:"
"80so 8- 80 "- s --- 80 --- :..
S 100 100 I'D 100 100 0
S22.5 11 5
E 215. a- 13
2 0 0<00
20 20 0 20 20
0 40 40 0. 40 ---. 40
q.. 6. 60 6,* 60 -0 --....
8 80 oc80 --""-O80 P40.
S 100100 100 100 100 U0''"
Figure 2-3. Sensitivity of contact results to grid density. Sensitivity of A) contact force,
B) contact torque, C) maximum pressure and D) average pressure to fixed
Forward dynamic simulations using the contact DLL also require the use of
implicit integration, since the elastic contact model makes the system equations
numerically stiff. A higher order implicit integrator is needed here to maintain simulation
accuracy. In addition, to maximize computational performance, a coarse element grid is
used to reduce the number of distance calculations. This is a reasonable since contact
forces and torques are much less sensitive to grid density than are contact pressures and
areas (Figure. 2-3). Turning off generation of output data files also significantly improves
performance. Once a forward dynamics simulation has been generated, an inverse
dynamics analysis using a fine element grid can be performed at little computational cost.
This permits more accurate prediction of contact pressures and areas. Output data files
can also be generated with data reported only at the desired intervals.
2.3 Sample Applications
Two sample applications are provided to demonstrate that the proposed
computational methodology works well within a multibody dynamic simulation
environment. The first application involves static analysis of a natural knee model, and
the second dynamic simulation of an artificial knee model. A commercial multibody
dynamics simulation code (Pro/MECHANICA MOTION, Parametric Technology
Corporation, Waltham, MA) was used as the test bed for both knee models, where the
tibia was welded to ground and the femur connected to the tibia via a six degree-of-
freedom joint. This joint provided the relative kinematic inputs required by the contact
DLL, which was interfaced to the multibody code via a custom load. For the static
analysis, an implicit Euler numerical integrator was used to provide numerical damping,
while the implicit integrator DASSL (Brenan et al. 1996) was used for the dynamic
simulation. All simulations were performed on a 2.4 GHz Pentium IV workstation.
2.3.1 Static Analysis of Natural Knee Contact
The natural knee contact model was created from MRI and CT data collected from
a single cadaver knee. Institutional review board approval was obtained for the data
collection and subsequent modeling. Prior to scanning, the specimen was cut
approximately 15 cm above and below the joint line. Sagittal plane MRI data were
collected using a 3.0-T GE Signa Horizon LX scanner with a quadrature knee coil. A T2-
weighted 3D FastGRE sequence was used with a 1 mm slice thickness, 256 x 256 image
matrix, and 160 mm x 160 mm field of view. Axial plane CT data were collected from
the same specimen using a GE LightSpeed QX/i scanner (General Electric Corporation,
Milwaukee, WI) in helical mode. The scanning parameters were a 1.25 mm overlapping
slice thickness, 512 x 512 image matrix, and 160 mm x 160 mm field of view. The tibia
and femur in both data sets were segmented using commercial image processing software
(SliceOmatic, Tomovision, Montreal, CA). Articular cartilage and subchondral bone
surfaces were segmented manually from the MRI data, while cortical bone surfaces were
segmented semi-automatically from the CT data using a watershed algorithm. The
resulting point clouds from both scans were exported for subsequent surface creation.
Commercial reverse engineering software (Geomagic Studio, Raindrop Geomagic,
Research Triangle Park, NC) was used to convert the point cloud data into a combined
geometric model for contact analysis. Point clouds from each bone and imaging modality
were loaded separately and automatically converted to polygonal surface models. The
subchondral bone surfaces from MRI were then automatically aligned to the
corresponding cortical bone surfaces from CT, creating a composite geometric model
with artricular cartilage surfaces from MRI and cortical bone surfaces from CT. NURBS
surfaces were fitted to the polygonal models, and the articular cartilage and subchondral
bone surfaces from MRI merged into single NURBS patches for subsequent contact
analysis. The number of u-v control points for the articular cartilage and subchondral
bone surfaces was minimized in commercial surface modeling software (Rhinoceros,
Robert McNeel & Associates, Seattle, WA). The tolerance between the original point
cloud from MRI and the final NURBS contact surfaces was 0.19 and 0.14 mm on the
average for the tibia and femur respectively.
The articular cartilage and cortical bone surfaces for the tibia and femur were
imported into Pro/MECHANICA MOTION to construct a contact model for the static
analysis. The articular cartilage and subchondral bone surfaces needed by the DLL were
read in separately. The meniscus was omitted in the model. Young's modulus for
articular cartilage was set to 4 MPa (Cohen et al. 2003) and Poisson's ratio to 0.45
(Blankevoort et al. 1991) to represent a relatively short time-frame response. Minimum
distance was used for the distance calculations, the contact surface u-v spaces were 30 x
30 (medial side) and 20 x 31 (lateral side) for the tibia, and 11 x 13 (both sides) for the
Initial Position Static Position
Figure 2-4. Visualization of example dynamic simulation for settling the femur onto the
tibia. Time flows from left (initial configuration) to right (final static
configuration). The dynamic simulation is terminated when all translational
and rotational accelerations are less then a small user-defined tolerance.
For a static analysis, the mass properties do not affect the final solution and so were
set to large enough values to permit dynamic simulation. The femur was placed above the
tibia and subjected to a 1000 N vertical load to bring it into contact with the tibia. Drivers
were added to the flexion-extension, internal-external rotation, and anterior-posterior
translation. The remaining degrees of freedom were adjusted by the dynamic simulation
until the femur reached its final static configuration. The contact element grid was set to
20 x 15 in both compartments. The static analysis successfully settled the femur onto the
tibia, adjusting the three free degrees of freedom as necessary to achieve the final pose
(Figure. 2-4). The static analysis was slow initially when the femur was out of contact,
since all elements on the tibia are searched when no contact is occurring. Once both sides
were in contact, the numerical integrator settled the femur into its final static
configuration within 62 seconds, with the contact area on the medial side being more
anterior than on the lateral side (Figure. 2-5).
Figure 2-5. Contact area of natural knee at the static position
2.3.2 Dynamic Simulation of Artificial Knee Contact
The artificial knee contact model was created from manufacturer CAD geometry of
an Osteonics 7000 cruciate-retaining knee implant (Stryker Howmedica Osteonics, Inc.,
Allendale, NJ). The contact surfaces were extracted from the tibial insert and femoral
component and imported into Geomagic Studio for resurfacing. After conversion to
polygonal surface models, a new network of untrimmed NURBS patches was created and
merged for each contact surface. The number of u-v control points for contact surfaces
was minimized in Rhinoceros, with the tolerance between the original and re-surfaced
geometry being + 0.002 mm.
The original CAD geometry for the tibial insert, tibial tray, and femoral component
were imported into Pro/MECHANICA MOTION to construct a contact model for
dynamic simulation. The contact surfaces were read in separately by the DLL. Young's
modulus for the ultrahigh molecular weight polyethylene tibial insert was set to 463 MPa
(Kurtz et al. 2002) and Poisson's ratio to 0.46 (Bartel et al. 1995). Minimum distance was
used for the distance calculations, the tibial insert and femoral component contact surface
u-v spaces were 40 x 25 and 15 x 20, respectively, and the contact element grid was set to
30 x 30 in both compartments. Mass properties for all bodies were calculated from their
volumes assuming uniform mass density. To simulate a clinically realistic situation,
flexion-extension, internal-external rotation, and anterior-posterior translation of the
femoral component were prescribed based on fluoroscopically measured in vivo
kinematics from a patient with the identical implant (Fregly et al. 2004). The femoral
component was loaded axially off-center to produce a 70% medial-30% lateral load split
(Hurwitz et al. 1998, Schipplein and Andriacchi 1991) using a vertical ground reaction
force curve scaled to between 0.25 and 3.0 body weight (Lu et al. 1997, Taylor et al.
1998, Taylor and Walker 2001). Superior-inferior translation, medial-lateral translation,
and various-valgus rotation were predicted via numerical integration of the equations of
motion, with a preliminary static analysis performed to settle the femoral component onto
the tibial insert.
The dynamic simulation was able to predict the motion of the three free degrees of
freedom (Figure. 2-6), the contact forces in the medial and lateral compartments, and the
net contact force and torque applied to both bodies (Figure. 2-7). Since the load split was
fixed, the knee adduction torque predicted by the contact model mirrored the axial load
applied to the femoral component, similar to experimental adduction torque
measurements (Hurwitz et al. 1998, Schipplein and Andriacchi 1991). A subsequent
inverse dynamics analysis using a refined element grid predicted the contact area and
pressure distribution throughout the gait cycle (Figure. 2-8). The forward dynamic
simulation required approximately 10 minutes of CPU time to complete and the inverse
dynamics simulation required less than 30 seconds.
t = 0.0 sec t = 0.17 sec t = 0.34 sec t= 0.51 sec
t = 0.6 sec t = 0,85 sec t = 1 02 sec t = 122 sec
Figure 2-6. Dynamic simulation of contact during a cycle of gait for Osteonics7000
implant design. The time period of the gait cycle is 1.22 seconds.
500 A 20 B
0 _--- --- ---- 15 -
-1000 0 -/" -
SFx -10. A TX
-115 \ /
-2500 -20 .
0 02 04 06 08 1 12 14 0 02 04 06 08 1 12 14
Time (sec) Time (sec)
Figure 2-7. Dynamic simulation results. Contact forces A) and torques B) during a gait
cycle of Osteonics7000 implant design.
Heel Strike Foot Flat jHeel Lift 0 Toe Off
I Comact Period 2 J7 M.lsrTance Peric 7 Pou e PrOPiO ioe Per 1 00%
Figure 2-8. Contact areas during the contact phase of a gait cycle.
The methodology of dynamic simulation of contact in the human movement in both
natural and artificial knees has been presented in details in previous sections. Some of the
applications of the methodology shown that it is able to meet the requirement of studying
the contact pressure dynamically efficiently, which is important for the investigation of
the contact mechanism of the joint.
Once the element distances are available, a variety of contact models can be
implemented. We have focused on two versions of an elastic foundation model since this
model is in line with the characteristics of human joints. However, we have also
implemented an elastic half-space model within the same infrastructure. It is not
discussed here since we have found that it does not reproduce experimental results when
applied to artificial knees (Fregly et al. 2003).
The DLL of contact model can be easily incorporated into any multibody dynamic
simulation environment using the specified API. This includes dynamic musculoskeletal
models generated with SIMM and the Dynamics Pipeline (Musculographics Inc.,
Evanston, IL), which has been popular in the field of human body modeling. Within this
larger simulation environment, force-generating elements representing muscles and
ligaments can be included. The predicted time histories of contact pressures, slip
velocities, or other quantities of interest can be output on an element-by-element basis
and used for subsequent analyses, such as computational wear analyses (Fregly et al.
Though our study has focused on the computational methodology used to produce
working dynamic simulations incorporating elastic joint contact, comparisons with
experimental data indicate that the approach produces accurate contact pressure
predictions as well. Fregly et al. (2003) compared contact pressures calculated by the
model with static pressure measurements made with a Tekscan K-scan sensor. An
Exatech Optetrak B posterior-stablized knee implant (Exatech Corporation, Gainesville,
FL) was used in our study. Both the linear and nonlinear material models were evaluated,
and the linear model was able to reproduce the peak and average experimental contact
pressures simultaneously for 16 combinations of flexion angle and axial load. Fregly et
al. (2004) used dynamic simulations of the Osteonics 7000 knee implant to predict mild
wear on the tibial insert. Model predictions were compared to the tibial insert retrieved
post-mortem from the patient who provided the in vivo fluoroscopic kinematic inputs.
The maximum wear depths (medial and lateral) predicted by the model were within 0.1
mm of the retrieved insert, the locations of maximum wear matched the retrieval, and the
damage areas and volumes were qualitatively and quantitatively similar to the retrieval.
Ongoing work is evaluating the ability of the large deformation model to reproduce
experimentally measured static contact pressures in a natural knee model without the
meniscus. Preliminary results indicate that contact forces and average pressures will be
well predicted, while maximum pressures will be less well matched. Such evaluations are
valuable for determining the capabilities and limitations of the current contact models so
that they can be refined and improved as new understanding is gained.
On of the limitations of the methodology is that the contact model is elastic and
therefore quasi-static, meaning that the contact surfaces are assumed to equilibrate
infinitely fast as the dynamic simulation progresses. Inclusion of viscoelastic effects
would require the introduction of additional states to track the changes over time of the
individual element deformations. This is not a serious limitation if only contact forces or
pressures are of interest and an appropriate aggregate modulus is chosen for the cartilage.
Plasticity is also currently lacking in the model (Giddings et al. 2001, Kurtz et al. 1998),
though this feature could be easily added via a new material model.
No meniscus is present in the current formulation. The meniscus is difficult to
model even in finite element software (Donahue et al. 2002). An approach similar to that
of Li et al. (1999), where the menisci are modeled as spring elements, could be employed
as a starting point. Bones treated as rigid for the natural knee model. However, using a
detailed finite element model of a natural knee, Donahue et al. (2002) recently reported
that changing the bone material properties from deformable to rigid changed contact-
related prediction from the model by less than 2%. Thus, this does not appear to be a
serious limitation. For a natural knee, non-homogeneity and orthotropy of articular
cartilage may be important, depending on the intended application. Modifications to the
contact model could be made to attempt to address these issues if they are found to be
Multiple contact areas will not be detected if they are not present at the start of the
simulation when the contact elements are initialized or if they do not split off from other
contact regions so that they can be tracked. Merging of two known contact regions into is
handled properly. In most cases, only one contact region will exist between the two
contact surfaces. No algorithm has been included to quickly detect out of contact
situations. A fact contact-rejection algorithm based on bounding boxes could be
implemented to detect obvious cases, but since the contact surfaces remain in contact in
most situations, this is not a serious issue.
In the geometry preprocessing, an element surface be created on a simple shape
such as a plane, cylinder, or sphere. While the planar element surface currently
implemented works well for the tibiofemoral joint if the tibia is the fixed body and the
patellofemoral joint if the patella is the fixed body, it will not work well for the ankle or
hip. A cylindrical element surface could work well for the ankle if the talus is treated as
the fixed body, while a spherical element surface could be appropriate for the hip with
either the femoral head or acetabular cup treated as the fixed body. A fixed grid size is
currently used for the contact elements. Though the size can be specified by the user at
run time, this still limits the accuracy of the computed contact forces since fewer
elements will be in contact for shallow penetrations. Ideally, an adaptive element scheme
could be implemented, though this would increase the complexity of the solution process.
For now, the element grid must be specified such that a sufficient number are active in
both compartments for the lightest load.
Distance calculations are still the primary computational bottleneck in the
simulations. To further reduce CPU time, multithreading of the distance calculations on
multiprocessor machines has been investigated though not yet implemented. Since
contact pairs are independent during the contact calculations, each of them could be
analyzed with a separate thread to take full advantage of available computational
Though the current implementation is efficient computationally, additional
improvements would facilitate the use of the model in optimization studies involving full-
body musculoskeletal models. Dynamic simulations requiring 10 minutes of CPU time
currently necessitate the use of parallel processing to perform optimization studies
utilizing large-scale musculoskeletal models.
In summary, a detailed computational methodology for incorporating a deformable
contact model of the knee into multibody dynamic simulation software has been
presented. The current implementation works for the tibiofemoral joint of artificial knees
or natural knees without the meniscus. The methodology can predict contact forces,
pressures, and areas and is sufficiently fast computationally to perform dynamic
simulations of the tibiofemoral joint in as little as 10 minutes. The specified API allows
the contact model to be incorporated into any multibody dynamics simulation code.
Future work will incorporate this contact model into a full-body dynamic musculoskeletal
model generated by SIMM and the Dynamics Pipeline to study knee joint contact
mechanics during functional activity.
EXPERIMENTAL EVALUATION OF AN ELASTIC FOUNDATION MODEL TO
PREDICT CONTACT PRESSURES IN KNEE REPLACEMENTS
Wear remains a primary factor limiting the life span of TKRs. Liberated
polyethylene wear debris can initiate osteolysis (i.e., bone destruction) resulting in pain
and implant loosening. Researchers currently have three basic options for studying wear:
1. Analyze implants retrieved after failure, 2. Analyze implants retrieved post-mortem, or
3. Analyze implant wear test results. Ideally, implants prone to failure would be identified
before such designs are used in patients. While revision and post-mortem retrievals are
valuable for studying insert damage modes (Bartel et al. 1986), they can be difficult to
obtain and take years before becoming available (Harman et al. 2001). Physical wear
testing is essential, and recent knee simulator designs are becoming more successful at
reproducing the wear patterns observed in retrievals (Walker et al. 1997). However, a
single test can cost tens of thousands of dollars and take months to run.
A computational wear model is an attractive solution to these limitations
(Sathasivam and Walker 1998). Required inputs to such a model are in vivo tibial insert
surface kinematics and contact pressures. Deformable body contact analyses, such as
finite element analyses (FEA) (Bartel et al. 1986 and 1995, Bendjaballah et al. 1997,
D'Lima et al. 2001, Otto et al. 2001, Perie and Hobatho 1998, Rawlinson and Bartel
2002, Sathasivam and Walker 1998 and 1999), elasticity analyses (Bartel et al. 1986, Jin
et al. 1995, Rawlinson and Bartel 2002), and elastic foundation analyses (Blankevoort et
al. 1991, Li et al. 1997, Nufio and Ahmed 2001; Pandy et al. 1997), can predict contact
pressures but usually only under static conditions. In contrast, rigid body contact analyses
using multibody dynamic simulation methods can predict knee motion efficiently (Abdel-
Rahman and Hefzy 1998, Godest et al. 2000, Piazza, and Delp 2001) but cannot predict
contact pressures (Cheng et al. 1990). Dynamic FEA codes have begun to bridge the gap
between these two approaches but take hours or days of CPU time to predict motion and
contact pressures simultaneously. Thus, no fast simulation approach exists to provide the
required wear model inputs, making it difficult to perform design sensitivity or
optimization studies of TKR wear.
Our study experimentally evaluates a novel modeling approach for predicting knee
implant contact pressures during a dynamic task. The approach combines the best
features of the rigid and deformable body modeling methods mentioned above. It
integrates traditional multibody dynamics to predict large overall motions with
deformable body contact to predict contact pressures between general three-dimensional
surfaces. The resulting approach is sufficiently fast computationally to perform dynamic
knee simulations in minutes rather than hours or days. Static experimental contact
pressure data collected from a commercial knee implant were used to evaluate the
applicability and limitations of the current formulation, which uses an elastic foundation
contact model with linear and nonlinear polyethylene material properties. The ability to
predict static pressures is a necessary first step to predicting dynamic pressures and
3.2 Materials and Methods
3.2.1 Elastic Contact Model
An elastic foundation contact model (Johnson 1985; also called a rigid body spring
model An et al. 1990, Blankevoort et al. 1991, Li et al. 1997) was implemented within a
commercial multibody dynamics software program (Pro/MECHANICA MOTION,
Parametric Technology Corporation, Waltham, MA). The contact model is a dynamic
link library that can be integrated into any multibody dynamics code and uses the ACIS
3D Toolkit (Spatial Corporation, Westminster, CO) to perform geometry evaluations
between general three-dimensional surfaces. The model uses a "bed of springs" scattered
over the surfaces of the contacting bodies to push the surfaces apart (Johnson 1985). The
springs represent an elastic layer of known thickness covering a rigid substrate on one or
both bodies, where each spring is independent from its neighbors. If both bodies possess
an elastic layer of the same material, then the two layers may be treated as a single layer
of combined thickness (Li et al. 1997). Layered contact is in contrast to half-space
contact, where both bodies are elastic and semi-infinite. The assumption of independent
springs eliminates the integral nature of contact problems, thereby greatly simplifying the
analysis of conformal geometry (e.g., a sphere in a spherical cup) or nonlinear materials.
For a rigid femoral component contacting an ultrahigh molecular weight
polyethylene (UHMWPE) tibial insert of finite thickness, the contact pressure p for any
spring can be calculated from (Johnson 1985, An et al. 1990, Blankevoort et al. 1991)
p= (1- v)E(p) d (3-1)
p- d (3-1)
(1+ v)(1- 2v)h
where E(p) is Young's modulus of the elastic layer, which can be a nonlinear function of
p, v is Poisson's ratio of the elastic layer, h is the layer thickness at the spring location,
and d is the spring deflection, defined as the interpenetration of the undeformed surfaces
in the direction of the local surface normal. Note that d can be computed at each instant in
time given the current position and orientation of the tibial insert and femoral component.
In our study, the thickness h was calculated separately for each spring as the local insert
thickness in the superior-inferior direction. Springs were distributed approximately
uniformly over the tibial insert contact surfaces by projecting a planar rectangular
element grid onto the insert surfaces. The number of elements was adjusted manually for
each simulation such that approximately 100 springs were always "active" on each
contact surface at the final static position (Figure 3-1A).
2 1000 -- Experiment
0 10 20 30 40 50
m Stress (MPa)
0 20 40 60 80 100
Figure 3-1. Overview of materials and methods. A) 3D pressure plot demonstrating the
uniform distribution of springs over the tibial insert contact surfaces. B)
Comparison of the parametric nonlinear material model to experimental
UHMWPE Young's modulus-stress data reported by Cripton (1993) at
230C. C) Dynamic contact model using commercial knee implant CAD. D)
Sensitivity of predicted average and peak contact pressures to the number of
active elements in a contact region. Convergence to within 5% error by 50
elements. Assumed "true" values determined using 400 active elements per
side are indicated by gray lines.
Given the known value for the deflection d of any spring at each instant in time, the
contact pressure for the spring can be easily calculated. For a nonlinear material model,
an equation for E as a function of p (Cripton 1993) can be substituted into Eq. 3-1 to
produce one nonlinear equation for p (Nufio and Ahmed 2001). Since each spring is
independent of its neighbors, any standard nonlinear root-finding algorithm can solve the
resulting system of nonlinear pressure equations independently. For a linear material
model, Eq. 3-1 can be solved directly for p. The calculated pressures can then be
multiplied by their corresponding areas to produce a set of point forces. Finally, these
forces can be replaced with a single equivalent force and torque applied to the rigid
bodies for purposes of dynamic simulation (Kane and Levinson 1985). For comparison
with the experiments, average pressure was calculated by averaging all non-zero
3.2.2 Parametric Material Model
To facilitate the use of optimization in this solution process, a parametric material
model with a minimal number of parameters was developed to represent the nonlinear
properties of UHMWPE:
Se= o- -+-- (3-2)
2 ao 2 Uo
where c is strain, a is stress, and co, uo, and n are material parameters. This is the
nonlinear power-law material model discussed in Johnson (1985) but with the addition of
a linear term. This model has the advantages that only three parameters are required to
define a material, making the model easy to use for design sensitivity and optimization
studies, and that choosing n = 1 produces a standard linear material model. For a given
value of the current value of E = do/ds for any spring can be found from
E=11/ 1+nl o (3-3)
Equation 3-3 was fit to the experimental E versus -data reported by Cripton
(1993) at 23 and 37C. For both temperatures, n = 3 provided the best fit to the
experimental data (Figure 3-1B; R2 = 0.966 at 230C with so = 0.0257 and = 15.8; R2 =
0.925 at 37C with S0 = 0.0597 and a= 18.4).
For comparison purposes, both linear (n =1) and nonlinear (n = 3) UHMWPE
material models were used in our study. For both models, so was fixed and 6 tuned to
match the average contact pressure from the 0 flexion/3000 N experiment (see below),
since the largest load would subject the material to the greatest range of pressures. The
tuned parameters for the linear model were So = 1, M = 400, corresponding to a constant
Young's modulus of 400 MPa, which is close to the value of 463 MPa recently reported
by Kurtz et al. (2002). For the nonlinear model they were So = 0.0257, a = 15.9, very
close to material parameters determined from Cripton's data at 230C. Poisson's ratio for
the polyethylene was 0.46 for both material models (Bartel et al. 1995). These tuned
parameters were used to simulate the remaining 15 experimental cases (see below).
3.2.3 Dynamic Implant Model
Using this contact algorithm, a dynamic model with elastic contact was created
using CAD geometry from a commercial knee implant possessing medial-lateral
symmetry and a 9 mm minimum insert thickness (Figure 3-1C; Optetrak B, Exactech
Corporation, Gainesville, FL). To simply the geometry evaluations and eliminate
potential problems caused by seams between surface patches, medial and lateral femoral
contact surfaces (four per side) were replaced with single surface approximations using
Geomagic Studio (Raindrop Geomagic, Research Triangle Park, NC) and Rhinoceros
(Robert McNeel & Associates, Seattle, WA) software. The tolerance between the original
and approximate surfaces in the regions of contact was measured with Geomagic Studio
to be + 0.002 mm, well within manufacturing tolerance.
The dynamic model was constructed to emulate the experimental set-up described
below, which was designed to produce approximately equal contact forces, peak
pressures, average pressures, and contact areas in the medial and lateral compartments.
Deformable contact was permitted between the medial femoral condyle and insert
surface, the lateral femoral condyle and insert surface, and the femoral cam and tibial
post. A fixed force was applied vertically downward to the femoral component at the
approximate point of load application in the experiments. The femoral component was
connected to the tibial insert via a six degree-of-freedom (DOF) joint. Sagittal plane
translations and varus-valgus rotation were free while the remaining three DOFs were
fixed. The flexion angle was fixed at each of the four experimental values. To account for
small variations in experimental positioning (Liau et al. 1999), an optimization was
performed for each flexion angle to find the medial-lateral translation and internal-
external rotation that best matched the experimental peak and average contact pressures
simultaneously. This procedure was performed only for the 3000 N load, with the
resulting translation and rotation being used for the remaining three loads. The optimized
position was the same for both material models. Maximum medial-lateral translation was
0.35 mm and maximum axial rotation was 1.50, well within the expected experimental
With any numerical approach involving elements, an important issue is the number
of elements required to produce an accurate solution. To address this issue, we performed
a global sensitivity study to investigate the convergence of the peak and average contact
pressures as a function of the number of active springs. The sensitivity study was
performed using the 0 flexion/3000 N load case. Assuming the "true" solution
corresponded to 400 active elements per side, peak and average contact pressure were
within 10% error by only 25 active elements and 5% error by 50 active elements (Figure
3-1D). Thus, 100 active elements per side were more than accurate enough for
comparison with experimental data.
Figure 3-2. Visualization of example dynamic simulation used to settle the femoral
component onto the tibial insert. Time flows from left (initial configuration)
to right (final static configuration). The dynamic simulation is terminated
when all translational and rotational accelerations are less then a small user-
With the relative component positioning and accuracy requirements established, a
global sensitivity study was performed for each flexion angle to predict the variation in
peak and average contact pressure with load. For each load, the sensitivity study
performed a forward dynamic simulation using implicit integration with numerical
damping (rather than physical damping) to cause the components to settle together into a
static configuration (Figure 3-2). Mass and inertia of the femoral component were
computed from its geometry assuming uniform density. Each dynamic simulation was
terminated when all translational and rotational accelerations were less than a small user-
defined tolerance (le-5 mm/s2 and rad/s2 in our simulations). Typical computation time
for a single dynamic simulation with 100 active elements per side was less than 30
seconds on a 1.2 GHz Pentium III-M laptop computer. In all, 16 experimental cases (four
flexion angles with four loads per flexion angle; see below) were predicted by the
sensitivity studies. All optimization and sensitivity studies were performed with
Pro/MECHANICA MOTION's built-in design study capabilities.
3.2.4 Contact Pressure Experiments
A commercial knee implant identical to the CAD model described above was used
in the experiments (Figure 3-3A). Fixturing allowed the femoral component to be
positioned at flexion angles of 0, 30, 60, and 900 relative to the tibial component. The
femoral component was connected to the ram of an MTS servohydraulic test machine via
a pin joint that allowed varus-valgus rotation, thereby permitting equilibration of the
frontal plane contact moment (Harris et al. 1999). A lockable slider joint was located
above this pin joint to allow medial-lateral positioning of the vertical load axis above the
center of the femoral component (Harris et al. 1999). The tibial component was mounted
to the top of a tilt table with lockable varus-valgus rotation and horizontal plane
translations. To facilitate comparison with the model, the knee was tested in neutral
alignment with the degrees of freedom adjusted via trial and error until approximately
equal contact force, peak pressure, average pressure, and contact area were obtained in
the medial and lateral compartments.
E 150- 0.2 MPa
0 0.5 1 1.5 2
Pressure Cutoff (MPa)
Figure 3-3. Experimental setup and pressure measurement. A) Experimental contact
pressure measurements made on the same implant using a Tekscan K-Scan
pressure measuring system and a servohydraulic test machine. B) Sensitivity
of experimental contact area to the pressure cut-off value selected for the
Tekscan sensor. The change in contact area stabilizes by about 0.2 MPa, as
indicated by the gray line.
To provide a wide range of conditions for evaluating the model, 16 experimental
cases were used, composed of all possible combinations of four flexion angles (0, 30, 60,
and 900) and four loads (750, 1500, 2250, and 3000 N). The loads were representative of
approximately one to four times body weight. All tests were performed at 230C, and three
separate trials were performed for each case. Contact pressure and area data were
collected from the medial and lateral compartments of the implant using a Tekscan K-
Scan pressure measuring system with a fresh sensor (Harris et al. 1999). Though the
sensor will affect the pressure measurements (Wu et al. 1998), the effect is expected to be
much less than in natural joints where the sensor is stiffer than the contacting surfaces.
Results from both compartments were averaged (total of six measurements) to produce
mean and standard deviation data for evaluation of model predictions.
Four issues related to the Tekscan sensor required special attention. The first was
the development of an accurate calibration procedure. Since the K-Scan sensor has been
shown to exhibit significant calibration drift over time (Otto et al. 1999), each
experimental trial was self-calibrated. Ramp loads from 20% to 100% of desired load
were applied over 5 seconds, and these two points were used to perform the two-point
calibration procedure recommended by the manufacturer. The two loads required for
post-calibration were measured during each trial by an in-line load cell. Thus, the total
load measured by the sensor always matched the load applied by the MTS machine.
Because each trial was calibrated using load cell measurements from the same trial, trial-
to-trial calibration drift was eliminated.
The second issue was determination of an appropriate pressure cut-off value. Due
to the moderately conformal nature of the knee design tested, the sensor "crinkled"
slightly during the experiments (Lewis 1998). These crinkles introduce erroneous contact
pressures on sensels outside the true contact area, making it necessary to determine a
pressure below which all measured sensel pressures are set to zero. To determine this cut-
off, we plotted the experimentally measured contact area for the 0 flexion/1500 N load
case using pressure cut-off values ranging from 0 to 2 MPa (Figure 3-3B). While
changing the cut-off had little effect on the measured contact force or peak pressure, it
had a dramatic effect on the measured contact area and hence average pressure. In
particular, changing the pressure cut-off from 0 to 0.2 MPa resulted in a 50% drop in
contact area. Since little additional drop in area occurred beyond 0.2 MPa, we chose this
as our cut-off. Average experimental pressures were therefore calculated by ignoring all
sensels with pressures below 0.2 MPa.
The third issue was the methodology used to calculate peak pressures. Because
peak pressure can be sensitive to measurements from a single sensel, this quantity was
determined using two approaches. The first used the maximum pressure recorded by a
single sensel. The second applied the built-in averaging function available in the Tekscan
software prior to determining the sensel with maximum pressure. This function performs
weighted averaging of each sensel with its eight neighbors, which eliminates local "hot
spots" from a single sensel due to sensor crinkling, small surface imperfections, or non-
uniform response from the sensels. Since averaging reduces peak pressure measurements,
results without averaging can be viewed as an upper bound on the peak pressure and
results with averaging as a lower bound.
The final issue was estimation of pressure measurement errors due to sensor
discretization. The K-Scan sensor measures pressures on discrete sensels with an area of
1.61 mm2. Thus, both peak and average pressure measurements will be influenced by the
sensel size. To estimate the magnitude of discretization errors, a theoretical power law
relationship was used (see Fregly and Sawyer, 2003 for details). This relationship
estimates peak and average pressure errors as a function of a non-dimensional area
variable, defined as the ratio of the number of perimeter sensels in the contact patch to
the total number of sensels with pressure. At 0 flexion, calculation of this non-
dimensional variable from the active Tekscan sensels indicated errors on the order of 1-
2% for peak and 3-6% for average pressure. Though errors will be slightly larger at 900
where the contact area is smaller, these errors were deemed to be sufficiently small not to
hinder comparison with model predictions.
The linear and nonlinear contact models were evaluated by their ability to match
experimental peak and average contact pressures simultaneously. Matching both is
necessary for the predicted pressure distribution to match the experiments. Contact force
comparisons are not reported since the force is always matched exactly at the final static
configuration. Contact area comparisons are not reported since the area follows similar,
but inverted, trends to the average pressure (e.g., a slightly high predicted average
pressure means a slightly low predicted contact area).
Table 3-1. Comparison of linear and nonlinear material model. Mean, standard deviation
(SD), and root-mean-square (RMS) errors in predicted peak and average
contact pressures for the 16 experimental cases using the linear and nonlinear
Errors (MPa) Linear Nonlinear
Mean -2.2 -6.7
Peak pressure without averaging SD 1.2 3.6
RMS 2.5 7.5
Mean 2.0 -2.5
Peak pressure with averaging SD 1.4 2.7
RMS 2.4 3.6
Mean 0.31 0.56
Average pressure SD 0.49 1.1
RMS 0.57 1.2
The linear material model predicted the experimental data more closely than did the
nonlinear material model, with mean, standard deviation and RMS errors typically being
two to three times smaller for the linear model (Table 3-1). Both models tracked the
average contact pressure well, though the linear model tracked it better (Figure 3-4A).
Since the contact force was always matched exactly, the contact area was also predicted
well by both material models. In contrast, only the linear model tracked the peak pressure
well. The peak pressure predicted by the linear model was generally between the two
experimental measurements (with and without averaging), whereas that predicted by the
nonlinear model was consistently below the two experimental measurements except at
750 N (Figure 3-4B). Both models predicted posterior cam contact only at 900, which was
consistent with the experiments.
750 1500 2250 3000 750
201 60 20
750 1500 2250 3000 750
A Linear material model
Nonlinear material model
1500 2250 3000
1500 2250 3000
750 1500 2250 3000
750 1500 2250 3000 750 1500 2250 3000
Load (N) Load (N)
Experiments with averaging
Experiments without averaging
A Linear material model
Nonlinear material model
Figure 3-4. Comparison between experimental and predictions. A) average and B) peak
contact pressures for all possible combinations of four flexion angles and
four loads. Predictions were made with a single contact model using linear
and nonlinear material models. Experimental peak contact pressure
measurements were made with and without averaging of neighboring
sensels. Error bars on experiment data indicate + 1 standard deviation.
3.4.1 Contact Model Accuracy
Our study experimentally evaluates a hybrid rigid-deformable dynamic modeling
approach for predicting contact pressures in total knee replacements. Though the dynamic
simulations were used to predict final static configurations, the evaluation validates the
computational efficiency of the approach compared to dynamic finite element methods. It
also validates the accuracy of the approach for predicting static contact pressures with a
single set of material parameters, which is an essential first step toward predicting
dynamic contact pressures and eventually wear. Overall, these findings indicate that this
novel modeling approach is well suited to sensitivity and optimization studies of knee
An important issue affecting contact model accuracy is the choice of material
model. Some finite element studies of knee replacements have used nonlinear UHMWPE
material models (Bartel et al. 1995, D'Lima et al. 2001, Godest et al. 2002, Otto et al.
2001, Rawlinson and Bartel 2002) while others have used linear material models (Bartel
et al. 1986, Sathasivam and Walker 1998 and 1999). We are unaware of a previous study
that has performed a thorough contact pressure comparison between model and
experiments using actual implant geometry to validate either material model. At least one
other study performed at 230C found that a linear material model with a Young's
modulus of 500 MPa provided the best match to experimental contact area data measured
using loads between 500 and 2500 N (Stewart et al. 1995). In those experiments, a
spherical glass indenter contacted a thick, wide UHMWPE block that approximated a
half-space, with contact area predictions being made with a Hertz contact model. A more
recent study measured Young's modulus of UHMWPE to be 463 MPa using a miniature
specimen shear punch test (Kurtz et al. 2002). For an elastic foundation contact model,
since pressures on surrounding elements do not contribute to the deformation of an
element, a slightly lower value of Young's modulus is needed to produce the same total
deformation. This explains the best-fit Young's modulus of 400 MPa found in our study.
3.4.2 Contact Model Advantages
As with any engineering model, the proposed dynamic contact model has its
advantages and limitations. The primary advantages of the elastic foundation contact
model are its simplicity and versatility. Since each spring is treated as independent, more
computationally intensive coupled contact solutions involving either quadratic
programming (Conry and Seirig 1971, Kalker and Van Randen 1972) or repeated linear
system solution (Singh and Paul 1974) are not needed. This makes the model ideal for
incorporation into a multibody dynamic simulation framework. The independent nature
of the springs also makes it easy to incorporate nonlinear material models into the
formulation as described above. An elastic-plastic material model could be incorporated
as well by following a solution process similar to the nonlinear material model. In
addition, the formulation accounts for layered contact between bodies of finite and
variable thickness and breadth, and conformal contact involving non-planar contact
patches. Both of these features are limitations in elastic half-space contact models
(Johnson 1985). Because conformity is not an issue, the current contact model will be
applicable to more or less conformal knee implant designs than the one used here, as well
as to rotating platform designs.
Another advantage of this computational framework is that it is not restricted to the
elastic foundation contact model evaluated here. Within the same framework, we have
also implemented an elastic half-space contact model (Conry and Seirig 1971, Kalker and
Van Randen 1972) with modifications to account for conformal contact (Paul and
Hashemi 1981). Currently, this model cannot handle layered contact, contact patches with
dimensions comparable to the contacting bodies, or nonlinear materials. However, unlike
the elastic foundation model, it accounts for how pressure applied at one location
produces deformations at all locations. Data from this contact model are not presented
here since a single value of Young's modulus could not be found that matched the
experimental pressures for all loads and flexion angles. Nonetheless, this model
demonstrates that more complicated contact models can easily be implemented within the
same framework to address limitations in the elastic foundation model discovered
through future experimental evaluation.
Previous studies have also used commercial multibody dynamics codes to simulate
contact conditions in knee replacements. Similar to our study, the elastic contact model
developed for ADAMS (MSC Software, Santa Ana, CA) uses discrete compressive
springs on the surfaces of the contacting bodies. This model has been used to simulate
stability tests (McGuan et al. 1998), knee simulator machine motions (Rullkoetter et al.
1999), and cam engagement in posterior-stabilized knees (Metzger et al. 2001). However,
it predicts only contact forces, not contact pressures, and uses tessellated surface
approximations that can cause force discontinuities (Puse and Laursen 2002) that slow
numerical integration. Pro/MECHANICA MOTION's built-in elastic contact model uses
either a Hertz formulation for approximately quadratic surfaces or an elastic half-space
boundary element formulation for general surfaces. This model has been used to simulate
contact pressure as a function of flexion angle (Fregly 1999) and peak contact pressure
sensitivity to femoral component malrotation (Fregly 2001). While the use of the actual
surface geometry reduces contact force discontinuities, such half-space formulations have
significant limitations as discussed above. Neither of these contact models has been used
to predict experimental pressure data as done here.
3.4.3 Contact Model Limitations
Several important limitations are also present in the current formulation. First, the
contact model is quasi-static and does not account for viscoelasticity (Waldman and
Bryant 1997). A viscoelastic contact model would require adding states to keep track of
the deformation of each element, which would make the contact calculations more
complicated. Alternatively, if the elastic parameters are tuned to match dynamic contact
pressure data collected at a physiological loading rate and temperature, an equivalent set
of linear material parameters could be chosen to approximate the viscoelastic situation for
those loading conditions. Second, the model does not account for how pressure applied at
one location produces deformations at all locations. However, this was not a serious
limitation for matching our experiment data. Third, the model provides predictive
information on contact pressures but not surface tensile stresses or sub-surface stresses.
While the model would be useful for wear predictions made from contact pressure and
kinematic inputs, predictions involving sub-surface stresses would require finite element
analyses (Sathasivam and Walker 1998). Fourth, the current model does not account for
friction (Sathasivam and Walker 1997). However, since local slip velocity information is
available for each element, a Coulomb friction model can easily be added on an element-
by-element basis as a reasonable approximation (Johnson 1985). Fifth, the model has
only been evaluated for a knee implant, but given the versatility of the approach, it would
likely be applicable to other artificial joints as well (An et al. 1990, Li et al. 1997, Genda
et al. 2001).
The applicability of the linear material model to larger peak contact pressures
requires further investigation. The force applied to each side was a maximum of 1500 N.
During gait, it is possible that the entire load could be borne on one condyle (Stiehl et al.
1999). At 1500 N per condyle, the peak contact pressure in our experiments was between
25 and 35 MPa (Figure 3-4B). Assuming polyethylene has a yield stress of 14 MPa
(Bartel et al. 1995), yielding will be initiated below the surface when the peak contact
pressure reaches approximately 1.6 times the yield stress (Johnson 1985) or about 22.5
MPa. Thus, significant plastic deformation may not have occurred over the course of the
experiments. This hypothesis is supported by the model's prediction of only 4%
maximum strain, making it less surprising that the linear material model worked so well.
Implementation of an elastic-plastic contact model, as described above, may therefore be
required to model higher contact pressures.
PREDICTED SENSITIVITY OF KNEE IMPLANT WEAR TO INSERT THICKNESS
AND BODY MASS
Total knee replacement (TKR) longevity is limited in part by mild wear of the
ultra-high-molecular-weight-polyethylene (UHMWPE) tibial insert. Mild wear is affected
by factors such as the choice of polyethylene (Kurtz et al. 1998), femorotibial contact
stresses (Bartel et al. 1986), tibial insert thickness (Bartel et al. 1985), patient body mass
via its affect on joint loads (Schmalzried et al. 1998), patient activity level (Cornwall et
al. 2001, Harman et al. 2001, Schmalzried et al. 1998), and arthroplasty design (Mahoney
et al., 2002; McGloughlin and Kavanagh., 2000). In clinical practice, surgeons must take
all of these factors into consideration simultaneously, even though their interactive effect
on mild wear has not been quantified.
Previous studies have yielded valuable insights into how contact pressures and
surface kinematics influence mild wear. Static finite element and elasticity solutions
using simplified knee implant geometry have predicted increased surface damage with
increased contact pressures resulting from decreased insert thickness (Bartel et al. 1985,
Bartel et al. 1986, Bartel et al. 1995, Chillag and Barth 1991). However, patient activity
produces relative motion between the femoral component and tibial insert that plays an
important role in wear as well (Blunn et al. 1991, Schmalzried et al. 1998, Schmalzried et
al. 2000, Shepherd et al. 1999, Zahiri et al. 1998). This has been verified by in vitro wear
testing and fluoroscopic analysis of TKRs that were later retrieved, establishing an
important link between surface kinematics and surface damage (Cornwall et al. 2001,
Harman et al. 2001, Schmalzried et al. 1998). Nonetheless, surface kinematics alone are
not sufficient to predict damage depths.
More recently, actual implant geometry and dynamic modeling methods have been
used to study the combined effect of contact pressures and surface kinematics on mild
wear. Dynamic finite element analyses (FEA) have been used to simulate knee implant
components under well-defined loading conditions in a knee simulator machine
(Giddings et al. 2001, Godest et al. 2000, Godest et al. 2002). Though these studies have
predicted contact pressures and subsurface stresses, they have not been coupled with
wear analysis. An alternate approach combining multibody dynamic simulation with
elastic contact and wear modeling was able to predict maximum in vivo damage depths to
within 0.1 mm for one patient with good qualitative agreement for damage area and
volume lost as well. Simulation inputs were in vivo fluoroscopic measurements combined
with an assumed axial loading pattern. The tibial insert retrieved postmortem from the
same patient was used for validation (Fregly et al. 2004).
Among the factors that influence TKR mild wear, insert thickness and patient body
mass are two over which the surgeon and patient have control. Large body mass has a
tendency to increase wear, which is why surgeons routinely counsel TKR patients with
excessive body mass to lose weight prior to surgery (Parvizi et al. 2000). However, no
clinical data exits to estimate the wear reduction that can be achieved through a reduction
in body mass. These data would be helpful for assessing the relative benefit of increasing
insert thickness verses decreasing body mass, especially given the negative consequence
of increased bone resection to accommodate a thicker insert.
Our study predicted the combined effect of insert thickness and body mass on mild
wear of a TKR starting from a previously validated in vivo simulation. Computational
methods were used to evaluate 25 combinations of insert thickness (6, 8, 10, 12 and 14
mm) and patient body mass (50, 75, 100, 125, and 150 kg) during simulated gait activity.
The maximum wear, creep, and damage depth, wear area, and volume lost were predicted
for each combination. Multiple linear regression was used to determine simple equations
that can estimate mild wear under clinical conditions. The results are useful for
understanding the extent to which in clinical practice increased insert thickness can make
up for increased body mass, or alternatively, decreased body mass can permit decreased
4.2 Materials and Methods
The combined influence of insert thickness and body mass on mild wear was
predicted using computational methods. A nominal wear prediction was generated first
and evaluated against a retrieved tibial insert for which in vivo fluoroscopic data were
available. Neighboring wear predictions were then generated for additional combinations
of insert thickness and body mass. Details of the nominal and neighboring wear
predictions are provided below.
4.2.1 Nominal Wear Prediction
The development of a nominal wear prediction under in vivo conditions involved a
three-step process. Single-plane fluoroscopy was used to measure in vivo knee implant
kinematics during treadmill gait and stair activities (Figure 4-1A). The fluoroscopy data
were collected from a single patient (female, age 65 at surgery, height 170 cm, mass 70
kg) 21 months after surgery. The patient received an Osteonics 7000 cruciate-retaining
implant (Stryker Howmedica Osteonics, Inc., Allendale, NJ) with a minimum insert
thickness of 6.8 mm. The patient's tibial insert was retrieved postmortem after 51 months
of implantation when she was dead. The study was approved by the institutional review
board and the patient gave informed consent for both the fluoroscopic and retrieval
Figure 4-1. Overview of the computational framework for wear predictions. A) In vivo
fluoroscopic data used to provide patient-specific kinematic inputs to a
dynamic contact model of the same knee implant design. B) Multibody
dynamic model with elastic contact used to predict contact pressures and
slip velocities on the tibial insert surface during one gait cycle. C)
Computational wear model used to predict wear, creep, and total damage at
individual locations on the tibial insert surfaces based on outputs from the
dynamic contact model.
The fluoroscopic kinematics were then used to drive multibody dynamic
simulations incorporating elastic contact between the implant components. The contact
model utilized elastic foundation theory with linear material properties to calculate
contact pressures (Johnson 1985, An et al. 1990, Blankevoort et al. 1991, Li et al. 1997 ).
The values for the Young's modulus (463 MPa) and Poisson's ratio (0.46) were taken
from the literature (Bartel et al. 1995, Kurtz et al. 2002). The tibial insert contact surfaces
were discretized into an approximately rectangular grid of about 2500 elements, with
pressure computed on an element-by-element based on interpenetration of the
undeformed contact surfaces. The contact model was integrated into a multibody
dynamic simulation of the implant components (Figure 4-1B) constructed within
commercial software (Pro/MECHANICA MOTION, Parametric Technology
Corporation, Waltham, MA). Six degrees of freedom were permitted between the moving
femoral component and the fixed tibial insert, both of which were treated as rigid bodies
that deformed only in the regions of contact. Flexion-extension, internal-external rotation,
and anterior-posterior translation were prescribed from the fluoroscopic kinematics, while
superior-inferior translation, varus-valgus rotation, and medial-lateral translation were
predicted via forward dynamic simulation. One cycle was simulated for both gait and
stair, with the simulations predicting the time history of contact pressures and sliding
conditions on the tibial insert surface elements. An off-center vertical load was applied to
the femoral component to produce a 70% medial, 30% lateral load split (Hurwitz et al.
1998, Schipplein et al. 1991). Load magnitude for gait and stair activities was determined
by scaling vertical ground reaction forces curves from a patient of similar age, height,
weight, and knee flexion characteristics to be between 0.25 and 3 times body (Taylor et
al. 1998, Taylor et al. 2001).
Finally, a computational damage model (Figure 4-1C) was used to predict wear,
creep, and damage (= wear + creep) on an element-by-element basis from the time
history of contact pressures and sliding conditions predicted by the gait and stair
simulations. Predicted wear, creep and damage depths were defined as the material
removed along the local normal direction of each element on the tibial insert surface. The
wear model used Archard's wear law, while the creep model was derived from
experimental data in the literature (Archard and Hirst 1956, Lee and Pienkowski 1998).
Damage predictions for one cycle were extrapolated to the total time of implantation
assuming one million cycles per year of gait or stair activity (Schmalzried et al. 1998).
Predictions for combined gait-stair activity were generated using a linear rules-of-
mixture. The entire modeling and simulation process was evaluated by comparing
damage predictions from the model (Figure 4-2A) with damage observed on the tibial
insert retrieved postmortem from the same patient who provided the in vivo kinematics
(Figure 4-2B). Further details on the damage predictions can be found in (Fregly et al.
A B C B 0.90
Figure 4-2. Comparison of retrieval and predictions. A) insert damage observed upon
retrieval. B) insert damage predicted by the computational methodology
assuming an 85% gait, 15% stair partitioning of activities and C) insert
damage predicted by the computational methodology assuming 100% gait.
"X" indicates the locations of maximum damage.
4.2.2 Neighboring Wear Predictions
Staring from this nominal case, wear predictions were generated for 25 neighboring
cases comprised of all possible combinations of five insert thicknesses (6, 8, 10, 12 and
14 mm) and five patient body masses (50, 75, 100, 125, and 150 kg). These values were
chosen to cover a broad spectrum of clinical situations. All 25 simulations followed the
same steps as the nominal case for generating wear predictions. The simulations assumed
that the prescribed kinematics were identical for each test and that the vertical load was
scaled to patient body mass in the same way as in the nominal case.
Unlike the nominal case, the neighboring wear predictions were performed only for
gait activity. Although implanted components see a wide spectrum of activities, gait is
the primary physical activity undertaken by TKR patients and has the greatest effect on
wear (Schmalzried et al. 1998). Furthermore, for the nominal case, the gait simulation
alone produced damage predictions that were close to the retrieval in terms of location
and magnitude of maximum damage (Figure 4-2).
Once damage predictions were generated by the 25 gait simulations, multiple linear
regression was performed to quantify the sensitivity of predicted quantities to insert
thickness and body mass. The predicted quantities were maximum wear, creep, and
damage depth, wear area, and volume lost on both the medial and lateral side. The
regression functions took the form:
Prediction = ao + al x Insert Thickness + a2 x Body Mass (4-1)
where ao, al and a2 are the unknown regression coefficients. Insert thickness was in
millimeters and body mass in kilograms. Reporting the results in this format simplifies
the task of estimating wear under desired clinical conditions with body mass being able to
compare to insert thickness directly.
The gait simulations predicted the most damage for the thinnest insert with largest
body mass and the least damage for the thickest insert with smallest body mass (Figure 4-
3). Maximum wear depth was similar on the medial and lateral sides, ranging from 0.82
mm for the worst case to 0.30 mm for the best. Maximum creep depth on the medial side
was almost identical to maximum wear depth (maximum of 0.78 mm, minimum 0.32
mm), whereas on the lateral side it was much smaller (maximum 0.47 mm, minimum
0.18 mm). Consequently, the maximum damage depth was greater medially (maximum
1.57 mm, minimum 0.62 mm) than laterally (maximum 1.30 mm, minimum 0.51 mm).
Volume lost followed a similar trend between the two sides, with medial values
(maximum 239 mm3, minimum 83.5 mm3) being more than twice as large as lateral ones
(maximum 105 mm3, minimum 36.9 mm3). In contrast, wear area was largest for the
thickest insert with largest body mass, ranging from 472 to 345 mm2 on the medial side
and from 405 to 296 mm2 on the lateral one.
Figure 4-3. The damage plots for five combinations of insert thickness and body mass.
A) Best case scenario with thickest insert (T = 14 mm) and smallest body
mass (M = 50 kg). B) Corner scenario with thickest insert (T = 14 mm) and
largest body mass (M = 150 kg). C) Central scenario with median insert
thickness (T = 10 mm) and body mass (M = 100 kg). D) Corner scenario
with thinnest insert (T = 6 mm) and smallest body mass (M = 50 kg). E)
Worst case scenario with thinnest insert (T = 6 mm) and largest body mass
(M = 150 kg)."M" stands for patient body mass and "T" minimum insert
Three-dimensional surface plots of the total damage depth as a function of
insert thickness and body mass. A) damage depth on the medial side and B)
damage depth on the lateral side.
Table 4-1. Linear regressions to wear predictions. Linear regression results for
maximum wear, creep, and damage depth, along with damage area and
volume, for (a) the medial side and (b) the lateral side. The regression
equation is of the form Prediction = ao + al x Insert Thickness + a2 x Body
Prediction ao al a2 R2
Wear depth (mm) 0.307 -0.015 0.004 0.990
Creep depth (mm) 0.467 -0.022 0.003 0.975
Damage depth (mm) 0.767 -0.037 0.007 0.987
Damage area (mm2) 272.6 6.380 0.736 0.988
Damage volume (mm3) 38.52 -2.032 1.407 0.999
Prediction ao al a2 R2
Wear depth (mm) 0.458 -0.023 0.003 0.982
Creep depth (mm) 0.271 -0.014 0.002 0.983
Damage depth (mm) 0.730 -0.036 0.005 0.982
Damage area (mm2) 238.5 5.450 0.625 0.984
Damage volume (mm3) 18.70 -0.945 0.605 0.998
The predictions from all 25 simulations were extremely well fitted by multiple
linear regression. All R2 values were greater than or equal to 0.975 (Table 4-1), indicating
a strong bilinear relationship between insert thickness, body mass, and all predicted
quantities in the region investigated (Figure 4-4). Increasing insert thickness decreased all
predicted quantities except for damage area (ai columns all negative except for damage
area), while increasing body mass increased all predicted quantities (a2 columns all
positive). In general, a one millimeter change in insert thickness had a much larger
influence on damage-related quantities than did a kilogram decrease in body mass (lai /
a2l > 1). The ao values were not close to zero, indicating that the fit was only reasonable
in the region investigated.
Polyethylene insert thickness and patient body mass are important factors when
planning TKR surgeries. Their combined effect on polyethylene mild wear has not been
studied in previous investigations. Thus, no quantitative data exist to show how insert
thickness and body mass combine to affect wear and what the trade-off would be between
them. Starting from a validated nominal wear prediction, the current study used a
computational approach to investigate the mild wear sensitivity of polyethylene to these
two factors and predicted that an approximately bilinear relationship exists between them.
The modeling and simulation approach used in our study involved several
assumptions and limitations that could affect interpretation of the results. One assumption
was that the kinematic data and activity level were the same for all of the simulations. It
is possible that patients with different body mass would exhibit different knee kinematics
given the same implant design and surgical alignment. Furthermore, different patient
activity levels may affect the kinematics, loading, and resulting wear of the implant
(Schmalzried et al. 2000, Zahiri et al. 1998). The patient used in the present study could
be described as active, while patients with higher body mass tend to exhibit less activity
(Schmalzried et al. 1998). Since an increased applied load can be offset by a lower
activity level, increased body mass has not been correlated to increased wear clinically
(McClung et al. 2000). Thus, the present damage predictions can be viewed as a
conservative overestimate of the clinical situation. By holding kinematics and activity
level constant, a computational approach makes it possible to compare the effects of
insert thickness and body mass on damage-related quantities without other confounding
Although large deformation (Bartel et al. 1995) and plasticity (Giddings et al. 2001,
Kurtz et al. 2002) are important mechanical behaviors of polyethylene, a small strain
linear elastic polyethylene material model was used in the present study. There were three
reasons for this choice. First, peak and average pressures predicted by the linear elastic
material model under static conditions compared extremely well to experimentally
measured pressure data for a wide range of flexion angles and loads (Fregly et al. 2003).
A nonlinear elastic material model investigated in the same study did not perform as well.
Second, the nominal damage predictions generated using the linear elastic material model
were in close quantitative and qualitative agreement with the insert retrieved from the
same patient who provided the fluoroscopic kinematics (Fregly et al. 2004). Thus, there is
reason to believe that the neighboring predictions are reasonable as well. Third, fully
plastic deformation would not have occurred in any of the 25 simulations. Although
plastic deformation is initiated below the surface when the average contact pressure
reaches the yield stress, fully plastic deformation does not occur until the average
pressure reaches about three times the yield stress (Johnson 1985). The largest average
pressures in our simulations ranged from 5.76 to 27.9 MPa. Thus, assuming the yield
stress of polyethylene to be between 14 and 23 MPa (McGloughlin and Kavanagh 2000),
the highest average pressure at any point in any of the simulations was still less than
about twice the smaller yield stress value. In such situations, the contact pressures are still
given by elastic theory (Johnson 1985), suggesting that the simulation predictions were
Another assumption was that the tibial insert geometry did not need to be
progressively modified as the damage evolved. For some systems, accurate wear
prediction requires accounting for the coupled evolution of wear, kinematics, and load.
However, simultaneous coupling would require geometry updating after every cycle
along with simulation over millions of cycles, which is not currently practical
computationally. One trade-off would be to update the geometry a limited number of
times. For example, the first simulation could be used to extrapolate damage out to half
the total number of cycles. After updating, the geometry would reflect the new damage
state, and a second simulation could be performed for the remaining cycles. Validated
theoretical studies have been performed for cam-follower mechanisms to determine the
number of cycles that can be extrapolated from a single simulation before geometry
updating is necessary (Sawyer et al. 2002). Updating frequency was a function of
changes in the applied load and surface geometry. In knee implants, the applied loads are
not expected to change significantly as the polyethylene wears, and changes in the
surface geometry are relatively small compared to the principal radii of curvature of the
components. Thus, the system is expected to be only weakly coupled (Fregly et al. 2004)
so that omission of geometry updating should not affect the predictions significantly.
An important limitation of the study is that the results were generated for only a
single implant design, a sagittally nonconforming posterior cruciate ligament-retaining
implant. For more sagittally conforming designs, contact pressures will decrease, which
would tend to improve wear performance. However, larger contact areas and relative
motions would also occur, which would tend to degrade wear performance. The overall
result is not clear. A study on retrieved inserts reported a mean linear wear rate of 0.35
mm/year for curved inserts and 0.41 mm/year for flat ones, though this difference was not
statistically significant (Benjamin et al. 2001). Thus, it is not clear the extent to which the
present results would apply to more conformal knee designs.
As a valuable side benefit, the gait simulations can also be used to assess the ability
of contact pressure data alone to predict surface damage. A separate linear regression
analysis was performed to determine the correlation between the maximum pressure
predicted over each simulation and each of the damage-related quantities. While
maximum wear, creep, and damage depths were highly correlated with maximum
pressure (R2 > 0.976), wear volume (R2 > 0.732) and especially wear area (R2 > 0.084)
were less strongly correlated. Furthermore, for each simulation, the location of maximum
damage on each side tended to be very close to the location of maximum pressure during
the simulation. Linear regression of maximum average pressure calculated over each
simulation resulted in similar trends. These findings suggest that contact pressure
predictions generated via dynamic simulation may be useful surrogates for prediction of
maximum damage depth and location. However, a variety of activities would contribute
to overall wear characteristics under clinical conditions.
Though the linear regressions produced large R2 values for the range of parameters
used in the simulations, three dimensional surface plots (Figure 4-4) revealed a slightly
nonlinear trend for each predicted variable. Previous studies (Bartel et al. 1985, Bartel et
al. 1986) have shown that as insert thickness is decreased to small values (less than 4
mm), contact pressures increase rapidly. However, when the thickness is large (more than
6-8 mm), contact pressures became less sensitive to insert thickness. For static loads with
a fixed flexion angle of zero degrees, the same trend was produced by our model (not
shown). Direct comparison of the pressures predicted in our study with those reported in
the literature is difficult due to differences in implant design, kinematics, and applied
By dividing the a2 column by the al column in Table 4-1, one can predict how
losing body mass (which reduces joint loads) compares to increasing insert thickness. For
wear depth, each 10 kg of body mass lost is equivalent to gaining 1 to 2 mm of insert
thickness. Similar trends exist for other quantities with the exception of total damage
volume, where 6 to 7 mm of effective insert thickness is gained for each 10 kg lost. Since
the predicted volume loss is due to both wear and creep, the effective insert thickness
gained will be smaller for wear alone. Any prediction, however, is limited by a
recommended 6 to 8 mm minimum insert thickness to avoid high pressures on the
polyethylene (Bartel et al. 1985, Bartel et al. 1986). A 10 kg decrease in body mass
appears to be large compared to a 1-2 mm increase in insert thickness. Though it is easier
to reduce wear by increasing insert thickness, greater bone resection can complicate
potential revision surgery. While decreasing body mass poses no such future risk, it is
more difficult to achieve in clinical practice. Furthermore, factors besides body mass,
such as muscle co-contraction or joint stability, also affect joint loads and thus these
Comparison of the damage predictions with previously reported data is difficult due
to the dearth of quantitative results in the literature. To our knowledge, no retrieval data
exist to indicate the relationship between insert thickness, body mass, and wear when
patient activity level is taken into account. Some studies have reported no correlation
between wear rates and patient weight (Landy and Walker 1988), while others have
reported a negative correlation (Benjamin et al. 2001). However, applied loads in TKRs
increase and activity level decreases with increasing body mass. Therefore, the lack of a
clear body mass influence in retrieval studies may be largely due to decreased activity
In summary, the effect of insert thickness and body mass on mild wear of TKR
polyethylene inserts was investigated using computational methods. Wear was predicted
by dynamic simulation using in vivo fluoroscopic kinematic data in conjunction with
elastic contact and computational wear models. The results showed that increased body
mass (resulting in increased joint load) leads to higher wear while increased insert
thickness reduces it. This suggests that when the anatomy permits, there may be a benefit
to using a thicker insert for active patients and those with higher body mass. Multiple
linear regression of the predictions facilitates estimation of mild wear trade-offs for
particular clinical situations.
CONCLUSIONS AND FUTURE RESEARCH
This dissertation has presented a methodology to simulate deformable contact in
the knee joint during human movement. Details of the methodology are provided along
with sample applications for both natural and artificial knees. The methodology can
predict contact forces, pressures, and areas and perform dynamic simulations of the
tibiofemoral joint in a computationally efficient manner in as little as 10 minutes of CPU
time. The specified application programming interface allows the contact model to be
incorporated into any multibody dynamics simulation code. The contact simulations
developed based on the methodology are evaluated against experimental contact pressure
results and retrieved tibial insert from TKRs. The predictions compare well with
experimentally measured pressures and retrieval wear. The study of the sensitivity of
TKR insert wear to patient body mass and polyethylene thickness shows that the
simulations can be used to address clinicially important issues.
Although the methodology generates satisfactory results, some of its limitations are
addressed below, which will help direct future research to improve both model accuracy
and computational efficiency.
Despite its simplicity, the elastic foundation model is able to predict contact
efficiently with satisfactory accuracy, which makes it an excellent choice for the purpose
of human movement simulation. However, as pointed by some studies, plasticity
(Giddings et al. 2001, Kurtz et al. 1998) and large deformation (Bartel et al. 1995) play
an important role in the mechanical behavior of polyethylene in artificial knees. In natural
knees, the cartilage, meniscus and bone are all non-linear by nature. The inclusion of
nonlinear material properties makes the current contact model more realistic.
Fixed grids are currently used in the contact model. It is efficient and easy to use
the fixed grids, since once they are constructed in a preprocessing phase, no further
changes need to be made during later calculations. However, if the actual contact area is
very small, the contact area defined by the elements of the fixed grids will not be
accurate, making the contact pressures inaccurate as well. Implementation of adaptive
grids would be a good solution to this problem. Such grids would be constructed for
every contact area. If each grid had the same number of elements and the elements kept
the same aspect ratios, then all the elements would have a constant relationship to one
another and could be easily scaled. Adaptive grids can match small contact areas
accurately. However, they require a method for identifying potential contact areas
efficiently and distinguishing multiple contact areas if they exist.
Geometry preprocessing and deformation evaluation are the two aspects that take
most of the CPU time during a dynamic contact simulation. Faster distance finding
methods would bring a remarkable improvement in computational efficiency.
Multithreading of the distance calculations would be another option. Contact pairs are
independent during contact calculation and each of them could be processed and
evaluated with a separate thread.
MAXIMUM PRESSURES DURING SIMULATIONS
Table A-1. Maximum and average pressures during a cycle of gait simulation.
Maximum pressures (MPa)
Body Medial thickness (mm) Lateral thickness (mm)
(kg) 6 8 10 12 14 6 8 10 12 14
50 24.8 21.8 20.0 18.7 17.6 16.3 13.9 12.3 11.2 10.4
75 30.6 27.6 25.3 23.6 22.2 18.9 16.4 14.7 13.5 12.6
100 36.2 32.5 29.9 28.0 26.5 21.3 18.8 16.9 15.6 14.6
125 41.1 36.9 34.2 32.2 30.6 23.6 20.9 19.0 17.6 16.4
150 45.6 41.4 38.4 36.1 34.2 25.8 22.9 20.9 19.3 18.1
Average pressures (MPa)
Body Medial thickness (mm) Lateral thickness (mm)
(kg) 6 8 10 12 14 6 8 10 12 14
50 14.5 13.1 12.0 11.3 10.8 8.0 7.1 6.6 6.1 5.8
75 18.7 16.9 15.4 14.1 13.5 10.2 9.0 8.4 7.8 7.4
100 22.5 19.8 18.1 17.1 16.4 11.9 10.8 9.9 9.2 8.8
125 25.4 22.7 21.1 20.0 18.4 13.9 12.3 11.4 10.4 9.9
150 27.9 25.7 24.1 21.8 20.7 15.2 13.9 12.5 11.7 11.3
WEAR, CREEP AND DAMAGE DEPTH
Table B-1. Wear, creep and damage depth in the TKR insert.
Body Medial thickness (mm) Lateral thickness (mm)
(kg) 6 8 10 12 14 6 8 10 12 14
50 0.38 0.35 0.33 0.31 0.30 0.45 0.41 0.37 0.34 0.32
Wear 75 0.50 0.46 0.43 0.41 0.39 0.57 0.51 0.46 0.43 0.40
depth 100 0.61 0.56 0.53 0.50 0.48 0.66 0.59 0.54 0.50 0.47
(mm) 125 0.70 0.65 0.61 0.58 0.56 0.75 0.67 0.61 0.57 0.54
150 0.80 0.74 0.70 0.66 0.64 0.82 0.74 0.68 0.63 0.60
50 0.46 0.41 0.38 0.35 0.32 0.28 0.23 0.21 0.20 0.18
Creep 75 0.58 0.52 0.45 0.40 0.38 0.33 0.29 0.26 0.25 0.23
depth 100 0.67 0.56 0.51 0.48 0.46 0.38 0.34 0.31 0.29 0.27
(mm) 125 0.74 0.63 0.59 0.56 0.54 0.43 0.38 0.35 0.33 0.31
150 0.78 0.71 0.67 0.63 0.60 0.47 0.42 0.39 0.36 0.34
50 0.84 0.76 0.70 0.66 0.62 0.73 0.64 0.58 0.54 0.51
Damage 75 1.07 0.98 0.88 0.81 0.77 0.89 0.80 0.73 0.67 0.63
depth 100 1.28 1.12 1.04 0.98 0.94 1.04 0.93 0.85 0.79 0.74
(mm) 125 1.44 1.29 1.20 1.15 1.10 1.17 1.05 0.96 0.90 0.84
150 1.57 1.45 1.37 1.30 1.24 1.30 1.16 1.07 0.99 0.94
DAMAGE AREA AND VOLUME LOST
Table C-l. Damage area and volume lost in the TKR insert.
Body Medial thickness (mm) Lateral thickness(mm)
(kg) 6 8 10 12 14 6 8 10 12 14
50 345 359 370 383 391 296 310 322 333 341
Damage 75 368 384 397 406 417 321 334 346 354 359
area 100 386 403 415 425 436 336 348 359 370 379
(mm2) 125 401 417 429 442 455 349 361 373 385 393
150 413 431 446 460 472 360 372 384 398 405
50 94 90 87 85 84 42 40 39 38 37
Volume 75 132 127 124 120 118 59 56 55 53 52
lost 100 170 163 158 155 152 75 72 70 68 67
(mm3) 125 205 197 192 188 185 90 87 84 83 82
150 239 235 229 224 220 105 101 99 97 96
LIST OF REFERENCES
Abdel-Rahman, E. M. and Hefzy, M.S. (1998) Three-dimensional dynamic behavior of
the human knee joint under impact loading. Medical Engineering andPhysics 20,
Ahmadi, N., Keer, L. M., and Mura, T. (1983) Non-Hertzian contact stress analysis for an
elastic half-space normal and sliding contact. International Journal of Solids and
Structures 19, 357-373.
An, K. N., Himenso, S., Tsumura, H., Kawai, T., and Chao, E. Y. S. (1990) Pressure
distribution on articular surfaces: application to joint stability analysis. Journal of
Biomechanics 23, 1013-1020.
Anderson, F. C. and Pandy, M. G. (1999) A dynamic optimization solution for vertical
jumping in three dimensions. Computer Methods in Biomechanics and Biomedical
Engineering 2, 201-231.
Andriacchi, T. P. (1994) Dynamics of knee malalignment. Orthopaedic Clinics of North
America 25, 395-403.
Archard, J. F. and Hirst, W. (1956) The wear of metals under unlubricated conditions.
Proceedings of the Royal Society A236, 397-410.
Ateshian, G. A. (1993) A B-spline least-squares surface-fitting method for articular
surfaces of diarthrodial joints. Journal ofBiomechanical Engineering 115, 366-
Banks, S. A. and Hodge, W. A. (1996) Accurate measurement of three-dimensional knee
replacement kinematics using single-plane fluoroscopy. IEEE Transactions on
Biomedical Engineering 43, 638-649.
Bartel, D. L., Bicknell, V. L., and Wright, T. M (1986) The effect of conformity,
thickness, and material on stresses in ultra-high molecular weight components for
total joint replacement. Bone and Joint Surgery 68-A, 1041-1051.
Bartel, D. L., Burstein, A. H. and Edwards, D. L. (1985) The effect of conformity and
plastic thickness on contact stress in metal-backed plastic implants. Journal of
biomechanical Engineering 107, 193-199.
Bartel, D. L., Rawlinson, J. J., Burstein, A.H., Ranawat, C. S., and Flynn, W.F. (1995)
Stresses in polyethylene components of contemporary total knee replacements.
Clinical 0i Ith/lvoi, ,li \ and Related Research 317, 76-82.
Bendjaballah, M. Z., Shirazi-Adl, A., and Zukor, D. J. (1997) Finite element analysis of
human knee joint in varus-valgus. Clinical Biomechanics 12, 139-148.
Benjamin, J., Szivek, J., Dersam, G., Persselin,S. and Johnson, R. (2001) Linear and
volumetric wear of tibial inserts in posterior cruciate-retaining knee arthroplasties.
Clinical 0i Ithp q,,, li 1 and Related Research 392, 131-138.
Blankevoort, L., Kuiper, J. H., Huiskes, R., and Grrotenbeor, H. J. (1991) Articular
contact in a three-dimensional model of the knee. Journal ofBiomechanics 24,
Blunn, G. W., Walker, P. S., Joshi, A., and Hardinge, K. (1991) The dominance of cyclic
sliding in producing wear in the total knee replacement. Clinical 0i Ithlqipti \ and
RelatedResearch 273, 253-260.
Brenan, K. E., Campbell, S. L., and Petzold, L. R. (1996) Numerical Solution of Initial-
Value Problems in Differential-Algebraic Equations. SIAM Classics in Applied
Mechanics, 2nd edition, Philadelphia.
Chao, E.Y. and Sim, F.H. (1995) Computer-aided pre-operative planning in knee
osteotomy. Iowa Orthopedic Journal 15, 4-18.
Cheng, R. C.-K., Brown, T. D., and Andrews, J. G. (1990) Non-uniqueness of the
bicompartmental contact force solution in a lumped parameter mathematical model
of the knee. Journal ofBiomechanics 23, 353-355.
Chillag, K. J., and Barth, E. (1991) An analysis of polyethylene thickness in modular
total knee components. Clinical 0i ilthl vi i' \ and Related Research 273, 261-263.
Cohen, Z. A, Henry, J. H, McCarthy, D. M, Mow, V. C, and Ateshian, G. A. (2003)
Computer simulations of patellofemoral joint surgery. Patient-specific models for
tuberosity transfer. American Journal of Sports Medicine, 31, 87-98
Cohen, Z. A, Roglic, H, Grelsamer, R. P, Henry, J. H, Levine, W. N, Mow, V. C, and
Ateshian, G. A. (2001) Patellofemoral stresses during open and closed kinetic chain
exercises. An analysis using computer simulation. American Journal of Sports
Medicine 29, 480-487.
Conry, T. F. and Seirig, A. (1971) A mathematical programming method for design of
elastic bodies in contact. Journal ofApplied Mechanics 38, 387-392.
Cornwall, G. B., Bryant, J. T., and Hannsson, C. M. (2001) The effect of kinematic
conditions on the wear of ultrahigh molecular weight polyethylene (UHMWPE) in
orthopaedic bearing applications. Proceedings of the Institute of Mechanical
Engineers 215, 95-106.
Cripton, P.A. (1993) Compressive characterization of ultra-high molecular weight
polyethylene with applications to contact stress analysis of total knee replacements.
Master of Science Thesis. Queen's University, Kingston, Ontario.
Davoodi, R., Brown I. E., and Loeb G. E. (2003) Advanced modeling environment for
developing and testing FES control systems. Medical Engineering andPhysics 25,
Davoodi, R. and Loeb, G. E. (2001) Conversion of SIMM to SIMULINK for faster
development of musculoskeletal models. Proceedings oflFESS, Cleveland, Ohio,
Delp, S.L., Arnold, A. S., and Piazza, S. J. (1998) Graphics-based modeling and analysis
of gait abnormalities. Biomedical Materials and Engineering 8, 227-240.
Delp, S. L., Kocmond, J. H., and Stern, S. H. (1995) Tradeoffs between motion and
stability in posterior substituting knee arthroplasty design. Journal ofBiomechanics
D'Lima, D.D., Chen, P.C., and Colwell, C.W. (2001) Polyethylene contact stresses,
articular congruity, and knee alignment. Clinical 01 Ili, q,,pli \% and Related
Research 392, 232-238.
Donahue, T. L. H., Rashid, M. M., Jacobs, C. R., and Hull, M. L. (2002) A finite element
model of the human knee joint for the study of tibio-femoral contact. Journal of
Biomechanical Engineering 124, 273-280.
Elias, J. J., Wilson, D. R., Adamson, R., and Cosgarea, A. J. (2003) Evaluation of a
computational model used to predict the patellofemoral contact pressure
distribution. Journal of Biomechanics (in press).
Fregly, B. J. (1999) A three-dimensional compliant contact model for dynamic simulation
of total knee replacements. In Proceedings of the VIIth International Symposium on
Computer Simulation in Biomechanics, 10-13. August 5-7, University of Calgary,
Fregly, B.J. (2000) Effect of femoral component malrotation on contact stress in total
knee replacements. In Proceedings of the 24th Annual Meeting of the American
Society of Biomechanics. University of Illinois at Chicago, Chicago, IL.
Fregly, B. J., Bei, Y., and Sylvester, M. E. (2003) Experimental evaluation of a
multibody dynamic model to predict contact pressures in knee replacements.
Journal ofBiomechanics 36, 1659-1668.
Fregly, B. J., Sawyer, W. G., Harman, M. K., and Banks, S. A. (2004) Computational
wear prediction of a total knee replacement from in vivo kinematics. Journal of
Biomechanics, special issue on the knee (in press).
Genda, E., Iwasaki, N., Li, G., MacWilliams, B.A., Barrance, P.J., and Chao, E.Y.S.
(2001) Normal hip joint contact pressure distribution in single-leg standing-effect
of gender and anatomic parameters. Journal ofBiomechanics 34, 895-905.
Giddings, V.L., Kurtz, S.M., and Edidin, A.A. (2001) Total knee replacement
polyethylene stresses during loading in a knee simulator. Journal of Tribology 123,
Glitsch, U. and Baumann, W. (1997) The three-dimensional determination of internal
loads in the lower extreme. Journal of Biomechanics, 30, 1123-1131.
Godest, A. C., Beaugonin, M., Haug, E., Taylor, M., and Gregson, P. J. (2002)
Simulation of a knee joint replacement during a gait cycle using explicit finite
element analysis. Journal ofBiomechanics 35, 267-275.
Godest, A.-C., Simonis de Cloke, C., Taylor, M., Gregson, P. J., Keane, A. J.,
Sathasivam, S., and Walker, P.S. (2000) A computational model for the prediction
of total knee replacement kinematics in the sagittal plane Journal ofBiomechanics
Guilak, F., Butler, D. L., and Goldstein, S. A. (2001) Tissue engineering, cells, scaffolds,
and growth factors: functional tissue engineering. Clinical 0i Ithelve i 1i, and
RelatedResearch 391, S295-S305.
Halloran, J. P., Eas;ey, S. K., Penmetsa Janaki, Laz, P. J., Petrella, A. J., and Rullkoetter,
P. J. (2003) Efficient dynamic finite element rigid body analysis of TJR. 2003
Summer Bioengineering Conference, June 25-29, Key Blascayne, Florida, 551-552.
Harman M. K., Banks S. A., and Hodge W. A. (2001) Polyethylene damage and knee
kinematics after total knee arthroplasty. Clinical 0i Ithlvp, ,, li, and Related
Research, 392, 383-393.
Harris, Morberg, P., Bruce, W. J. M., and Walsh, W.R. (1999) An improved method for
measuring tibofemoral contact areas in total knee arthroplasty: a comparison of K-
scan sensor and Fuji film. Journal ofBiomechanics 32, 951-958.
Hefzy, M. S. and Cooke, T. D. V. (1996) Review of knee models: 1996 update. Applied
Mechanics Reviews 49, S 187-193.
Hoff, W. A., Komistek, R.D., Dennis, D.A., Gabriel, S. M., and Walker, S. A. (1998)
Three-dimensional determination of femoral-tibial contact positions under in vivo
conditions using fluoroscopy. Clinical Biomechanics 13, 455-472.
Hurwitz, D. E., Sumer, D. R., Andriacchi, T. P., and Sugar, D. A. (1998) Dynamic knee
loads during gait predict proximal tibial bone distribution. Journal of
Biomechanics 31, 423-430.
Jin, Z. M., Dowson, D. and Fisher, J. (1995) Contact pressure prediction in total knee
joint replacements. Part I: general elasticity solution for ellipticall layered contacts.
Jounral of Engineering in Medicine, 1995, 209(H1), 1-8.
Johnson, K. L. (1985) Contact Mechanics. Cambridge University Press, Cambridge.
Kalker, J. J. and Van Randen, Y. (1972) A minimum principle for frictionless elastic
contact with application to non-Hertzian half-space contact problems. Journal of
Engineering M1A aheiuh tii \ 6, 193-206.
Kane, T. R. and Levinson, D. A. (1985) Dynamics: Theory and Applications. McGraw
Hill, New York.
Kane, T. R., Likins, P. W. and Levinson, D. A. (1983) Spacecraft dynamics. McGraw-
Hill Book Co., New York.
Kanisawa, I., Banks, A. Z., Banks, S. A., Moriya, H., and Tsuchiya, A. (2003) Weight-
bearing knee kinematics in subjects with two types of anterior cruciate ligament
reconstructions. Knee Surgery, Sports Traumatology, Arthroscopy 11, 16-22.
Komistek, R. D., Dennis, D. A., and Mahfouz, M. (2003) In vivo fluoroscopic analysis of
the normal human knee. Clinical Oi th 1iop/itlii \ and Related Research 410, 69-81.
Kurtz, S. M., Bartel, D. L., and Rimnac, C. M. (1998) Post-irradiation aging affects the
stresses and strains in UHMWPE components for total joint replacement, Clinical
O 1 /i1'(t hik % and Related Research 350, 209-220.
Kurtz, S. M., Jewett, C. W., Bergstrom, J. S., Foulds, J. R., and Edidin, A.A. (2002)
Miniature specimen shear punch test for UHMWPE used in total joint
replacements. Biomaterials 23, 1907-1919.
Kwak, S. D., Blankevoort, L., and Ateshian, G. A. (2000). A mathematical formulation
for 3D quasi-static multibody models of diarthrodial joints. Computational
Methods in Biomechanical and Biomedical Engineering 3, 41-64.
Landy, M., and Walker, P. (1998) Wear of ultra-high-molecular-weight polyethylene
components of 90 retrieved knee prostheses. Journal of Arthroplasty Suppl. 3, 73-
Lee, K. Y. and Pienkowski, D. (1998) Compressive creep characteristics of extruded
ultrahigh molecular-weight polyethylene. Journal of Biomedical Material
Research 39, 261-265.
Lewis, G. (1998) Contact stress at articular surfaces in total joint replacements. Part I:
Experimental methods. Biomedical Materials and Engineering 8, 91-110.
Li, G., Rudy, T. W., Sakane, M., Kanamori, A., Ma, C. B. and Woo, S. L.-Y. (1999)
Importance of quadriceps and hamstring muscle loading on knee kinematics and in-
situ forces in the ACL. Journal ofBiomechanics Vol. 32, No. 4, 395-400.
Li, G., Sakamoto, M., and Chao, E.Y.S. (1997) A comparison of different methods in
predicting static pressure distribution in articulating joints. Journal of
Biomechanics 30, 635-638.
Liau, J.-J., Cheng, C.-K., Huang, C.-H., Lee, Y.-M., Chueh, S.-C., and Lo, W.-H. (1999)
The influence of contact alignment of the tibofemoral joint of the prosthesis in in
vitro biomechanical testing. Clinical Biomechanics 14, 717-721.
Lu, T. W. and O'Connor, J. J. (1999) Bone position estimation from skin marker co-
ordinates using global optimization with joint constraints. Journal of Biomechanics
Mahoney, O. M., McClung, C. D., Rosa, M.A. and Schmalzried, T. P. (2002) The effect
of total knee arthroplasty design on extensor mechanism function. Journal of
Arthroplasty 17, 416-421.
McClung, C. D., Zahiri, C. A., Higa, J. K., Amstutz, H. C. and Schmalzried, T. P. (2000)
Relationship between body mass index and activity in hip or knee arthroplasty
patients. Journal of Orthopaedic Research 18, 35-39.
McGloughlin, T. M. and Kavanagh, A. G. (2000) Wear of ultra-high molecular weight
polyethylene (UHMWPE) in total knee protheses: a review of key influences.
Proceedings of the Institute of Mechanical Engineers 214, 349-359.
McGuan, S., Jasty, M., and Kaufman, M. (1998) Total knee system performance
measurement through computerized intrinsic stability testing. Proceedings of the
65th Annual Meeting of the American Academy of Orthopaedic Surgeons. New
Metzger, R., Lombardi, A.V., Mallory, T.H., Fada, R.A., Hartman, J.F., Adamsn, J.B.,
and McGuan, S.P. (2001) Late versus early engagement of posterior stabilized
prostheses: Effect on extensor moment arm and resultant extensor loads. In
Proceedings of the 47th Annual Meeting of the Orthopaedic Research Society. San
Neptune, R. R., Wright, I. C., and Bogert, A. J. van den (2000) A method for numerical
simulation of single limb ground contact events: application to heel-toe running.
Computer Methods in Biomechanics and Biomedical Engineering 3, 321-334.
Nuio, N. and Ahmed, A.M. (2001) Sagittal profile of the femoral condyles and its
application to femorotibial contact analysis. Journal ofBiomechanical Engineering
Otto, J. K., Brown, T. D., and Callaghan, J. J. (1999) Static and dynamic response of a
multiplexed-array piezoresistive contact sensor. Experimental Mechanics 39, 317-
Otto, J. K., Callaghan, J. J., and Brown, T. D. (2001) Mobility and contact mechanics of a
rotating platform total knee replacement. Clinical 01 thIipa dil \ and Related
Research 392, 24-37.
Pandy, M. G., Sasaki, K., and Kim, S. (1997) A Three-dimensional musculoskeletal
model of the human knee joint. Part 1: Theoretical Construction. Computer
Methods in Biomechanics and Biomedical Engineering 1, 87-108.
Parvizi, J., F. R. C. S., Trousdale, R. T., and Sarr, M. G. (2000) Total joint arthroplasty in
patients surgically treated for morbid obesity. Journal of Arthroplasty 15, 1003-
Paul, B. and Hashemi, J. (1981) Contact pressures on closely conforming elastic bodies.
Journal ofApplied Mechanics 48, 543-548.
Perie, D. and Hobatho, M. C. (1998) In vivo determination of contact areas and pressure
of the femorotibial joint using non-linear finite element analysis. Clinical
Biomechanics 13, 394-402.
Piazza, S. J. and Delp, S. L. (1999) Three-dimensional simulation of total knee
replacement motion during stepup tasks. In Proceedings of the 1999
Bioengineering Conference, BED-Vol. 42, 423-424. The American Society of
Mechanical Engineers, New York.
Piazza, S. J. and Delp, S. L. (2001) Three-dimensional dynamic simulation of total knee
replacement motion during a step-up task. Journal ofBiomechanical Engineering
Praemer, A., Furner, S., and Rice, D. P. (1999) Musculoskeletal Conditions in the United
States. American Academy of Orthopaedic Surgeons, Rosemont, IL.
Puse, M. A. and Laursen, T. A. (2002) A 3D contact smoothing method using Gregory
patches. International Journal for Numerical Methods in Engineering 54, 1161-
Rawlinson, J. J. and Bartel, D. L. (2002) Flat medial-lateral conformity in total knee
replacements does not minimize contact stress. Journal ofBiomechanics 35, 27-34.
Rhinoceros Version 1.0 user's guide (2001). Chapter 28: About NURBS. Robert McNeel
& Associates, Seattle, WA.
Riegger-Krugh, C, Gerhart, T. N., Powers, W. R., and Hayes, W. C. (1998) Tibiofemoral
Contact Pressures in Degenerative Joint Disease. Clinical 0i Ith/lpi& 1i and
Related Research, 348, 233-245.
Rullkoetter, P. J., McGuan, S., and Maletsky, L. P. (1999) Development and verification
of a virtual knee simulator for TKR evaluation. In Proceedings of the 45th Annual
Meeting of the Orthopaedic Research Society, Anaheim, CA.
Sathasivam, S. and Walker, P. S. (1998) Computer model to predict subsurface damage
in tibial inserts of total knees. Journal of Orthopaedic Research 16, 564-571.
Sathasivam, S. and Walker, P.S. (1999) The conflicting requirements of laxity and
conformity in total knee replacements. Journal ofBiomechanics 32, 239-247.
Sawyer, W. G., Diaz, K. I., Hamilton, M. A., and Micklos, B.(2002) Evaluation of an
analytical model for the evolution of wear and load in a Scotch-Yoke mechanism.
Journal of Tribology (Accepted).
Schiehlen, W. (1997) Multibody system dynamics: roots and perspectives. Multibody
System Dynamics 1, 149-188.
Schipplein, O. D. and Andriacchi, T. P. (1991) Interaction between active and passive
knee stabilizers during level walking. Journal of Orthopaedic Research 9, 113-119.
Schmalzried, T. P., Shepherd, E. F., Dorey, F. J., Jackson, W. 0., dela Rosa, M., Fa'vae,
F., McKellop, H. A., McClung, C. D., Martell, J., Moreland, J. R., and Amstutz, H.
C. (2000) Wear is a function of use, not time. Clinical 0i I/thi qpl1li \ andRelated
Research 381, 36-46.
Schmalzried, T. P., Szuszczewicz, E. S., Northfield, M. R., Akizuki, K. H., Frankel, R.
E., Belcher, G.and Amstutz, H. C. (1998) Quantitative Assessment of Walking
Activity after Total Hip or Knee Replacement. Bone and Joint Surgery 80-A (1),
Setton, L. A., Elliott, D. M. and Mow, V. C. (1998) Altered mechanics of cartilage with
osteoarthritis: Human osteoarthritis and an experimental model of joint
degeneration. Oe,,ui thil itii and Cartilage 7, 2-14.
Shepherd, E. F., Toloza, E., McClung, C. D. and Schmalzried, T. P. (1999) Step Activity
monitor: increased accuracy in quantifying ambulation activity. Journal of
Orthopaedic Research 17, 703-708.
Singh, K. P. and Paul, B. (1974) Numerical solution of non-Hertzian elastic contact
problems. Journal ofApplied Mechanics, 484-490.
Stewart, T., Jin, Z. M., Shaw, D., Auger, D. D., Stone, M., and Fisher, J. (1995)
Experimental and theoretical study of the contact mechanics of five total knee joint
replacements. Proceedings of the Institution of Mechanical Engineers Part H 209,
Tashman, S. and Anderst, W. (2003) In-vivo measurement of dynamic joint motion using
high speed biplane radiography and CTL: application to canine ACL deficiency.
Journal ofBiomechanical Engineering 125, 238-245.
Taylor, S. J. G. and Walker, P. S. (2001) Force and moments telemetered from two distal
femoral replacements during various activities. Journal of Biomechanics 34, 839-
Taylor, S. J. G., Walker, P. S., Perry, J. S., Cannon, S. R. and Woledge, R. (1998) The
forces in the distal femur and the knee during walking and other activities measured
by telemetry. Journal ofArthroplasty 13, 428-437.
Tetsworth, K. and Paley, D. (1994) Accuracy of correction of complex lower-extremity
deformities by the Ilizarov method. Clinical 01 thIiope i~ % and Related Research
Waldman, S. D. and Bryant, J. T. (1997) Dynamic contact stress and rolling resistance
model for total knee arthroplasties. Journal of Biomechanical Engineering 119,
Walker, P. S., Blunn, G. W., Broome, D. R., Perry, J., Watkins, A., Sathasivam, S.,
Dewar, M, and Paul, J. P. (1997) A knee simulating machine for performance
evaluation of total knee replacements. Journal ofBiomechanics 30, 83-89.
Wismans, J., Veldpaus, F., Janssen, J., Huson, A. and Struben, P. (1980) A three-
dimensional mathematical model of the knee joint, Journal of Biomechanics, 13,
Wu., J. Z., Herzog, W., and Epstein, M. (1998) Effects of inserting a pressensor film into
articular joints on the actual contact mechanics. Journal ofBiomechanical
Engineering 120, 655-659.
Zahiri, C. A., Schmalzried, T. P., Szuszczewicz, E. S., and Amstutz, H.C. (1998)
Assessing activity in joint replacement patients. Journal of Arthroplasty 13, 890-
Zajac, F. E. (1993) Muscle coordination of movement: A perspective. Journal of
Biomechanics 26 supplyl 1), 109-124.
Yanhong Bei was born in Beijing, China. She received her Bachelor of Engineering
degree in spacecraft design from the Beijing University of Aeronautics and Astronautics
in 1996. Three years later, she finished her study at the Institute of Spacecraft System
Engineering of the Chinese Academy of Space Technology with a master's degree in
spacecraft system design. In 1999, she entered the Ph.D program at the University of
Florida, in the Department of Mechanical and Aerospace Engineering, working on
dynamic simulation of contact in the knee joint during human movement.