<%BANNER%>

Dynamic Stimulation of Knee Joint Contact during Human Movement


PAGE 1

DYNAMIC SIMULATION OF KNEE JOINT CONTACT DURING HUMAN MOVEMENT By YANHONG BEI A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2003

PAGE 2

This document is dedicated to my husband, Liu Jian.

PAGE 3

ACKNOWLEDGMENTS 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 studies. 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. iii

PAGE 4

TABLE OF CONTENTS Page ACKNOWLEDGMENTS.................................................................................................iii LIST OF TABLES.............................................................................................................vi LIST OF FIGURES..........................................................................................................vii ABSTRACT.....................................................................................................................viii CHAPTER 1 INTRODUCTION......................................................................................................1 1.1 Background...........................................................................................................1 1.2 Objective...............................................................................................................2 1.3 Dissertation Organization.....................................................................................3 2 METHODOLOGY FOR DYNAMIC CONTACT SIMULATION...........................5 2.1 Introduction...........................................................................................................5 2.2 Multibody Knee Contact Methodology................................................................7 2.2.1 Contact Surface Preparation.....................................................................8 2.2.2 Efficient Distance Calculations...............................................................11 2.2.3 Contact Model Development..................................................................17 2.2.4 Dynamic Model Construction.................................................................20 2.3 Sample 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 Materials and Methods........................................................................................37 3.2.1 Elastic Contact Model.............................................................................37 3.2.2 Parametric Material Model.....................................................................39 iv

PAGE 5

3.2.3 Dynamic Implant Model.........................................................................40 3.2.4 Contact Pressure Experiments................................................................43 3.3 Results.................................................................................................................47 3.4 Discussion...........................................................................................................49 3.4.1 Contact Model Accuracy........................................................................49 3.4.2 Contact Model Advantages.....................................................................50 3.4.3 Contact Model Limitations.....................................................................52 4 PREDICTED SENSITIVITY OF KNEE IMPLANT WEAR TO INSERT THICKNESS AND BODY MASS...........................................................................54 4.1 Introduction.........................................................................................................54 4.2 Materials and Methods........................................................................................56 4.2.1 Nominal Wear Prediction.......................................................................56 4.2.2 Neighboring Wear Predictions................................................................59 4.3 Results.................................................................................................................60 4.4 Discussion...........................................................................................................63 5 CONCLUSIONS AND FUTURE RESEARCH......................................................69 APPENDIX A MAXIMUM PRESSURES DURING SIMULATIONS..........................................71 B WEAR, CREEP AND DAMAGE DEPTH..............................................................72 C DAMAGE AREA AND VOLUME LOST..............................................................73 LIST OF REFERENCES...................................................................................................74 BIOGRAPHICAL SKETCH.............................................................................................83 v

PAGE 6

LIST OF TABLES Table page 2-1. CPU time tests of distance evaluation methods........................................................17 2-2. Input parameters and output measures in the API between the contact model and dynamic contact model.............................................................................................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 vi

PAGE 7

LIST OF FIGURES Figure page 2-1. Framework 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.......................................................................................................................29 2-7. Dynamic simulation 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 mass..........................................................................................62 vii

PAGE 8

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 MOVEMENT By Yanhong Bei December 2003 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 viii

PAGE 9

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 patients body mass or by losing weight, given the insert thickness. ix

PAGE 10

CHAPTER 1 INTRODUCTION 1.1 Background 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 1

PAGE 11

2 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, DLima et al. 2001, Pri 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, Nuo and Ahmed 2001, Pandy et al. 1997) is a simpler and more efficient way to predict the knee joint contact pressures. 1.2 Objective Combining deformable joint models with rigid body system is one of the important future research directions in multibody dynamics (Schiehlen 1997). However, no such

PAGE 12

3 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 artificial knee. 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

PAGE 13

4 dynamics software program. The models 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 patients 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 wear performance. Chapter 5 summarizes the research with recommendations for future improvements to the contact model and simulation efficiency.

PAGE 14

CHAPTER 2 METHODOLOGY FOR DYNAMIC CONTACT SIMULATION 2.1 Introduction 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. 5

PAGE 15

6 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.

PAGE 16

7 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 manufacturers 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

PAGE 17

8 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. 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.

PAGE 18

9 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 possible. 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 possible. 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

PAGE 19

10 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

PAGE 20

11 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 y 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 the y 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).

PAGE 21

12 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 at y = 0 on which a rectangular grid of elements is created. The center points of these elements are projected in the y 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

PAGE 22

13 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. A B C D Fixed Surface Fixed Surface Moving Surface Fixed Surface Moving Surface Moving Surface Tangent Plane 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

PAGE 23

14 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

PAGE 24

15 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

PAGE 25

16 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

PAGE 26

17 surface should be merged into a single NURBS patch with a minimal number of u-v lines whenever possible. Table 2-1. CPU time tests of distance evaluation methods. Tolerance statistics (mm) CPU time (micro-sec) Geometry Maximum Average STD Ray firing Min.distance Single surface (50) 4,459,921 879,488 Single surface (18) 0.0412 0.00448 0.0053 3,789,363 437,182 12-patch surface (8) 0.0486 0.00148 0.0019 2,175,888 2,461,545 24-patch surface (5) 0.0434 0.00202 0.0028 2,031,058 2,791,974 48-patch surface (4) 0.0397 0.00163 0.0026 1,929,121 3,442,914 96-patch surface (3) 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 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.

PAGE 27

18 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 p of any spring: dhEp)21)(1()1( (2-1)

PAGE 28

19 where E is Youngs modulus of the elastic layer, is the Poissons 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: hdEp 1ln)21)(1()1( (2-2) For nonlinear materials, E can be a function of the pressure p (Nuo and Ahmed 2001), which we derive from a three-parameter nonlinear material model (Fregly et al. 2003): noooopppp2121 (2-3) where is strain, p is contact pressure, and o , and 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 op n 1 n This material model has been shown to fit experimental data for polyethylene extremely well (Cripton 1993,

PAGE 29

20 Fregly et al. 2003). For an estimated value for p, ddpE c an be calculated from Eq. 2-3: 1121/1noooppnpE (2-4) Given the interpenetration d for any spring, Eq. 2-4 can be substituted into Eq. 2-1 or Eq. 2-2 to produce a single nonlinear equation for p 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 concerned. 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.,

PAGE 30

21 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 elastic half space) Material model (small or large strain) Material parameters (material properties, thickness offset) Friction parameters (static and dynamic friction coefficients) Damping parameter (nonlinear) Distance calculation method (ray firing or minimum distance) Element grid density (number of divisions along each direction) Output flag (off or on, and output format) Kinematics Contact force (x, y, z components and magnitude) on medial and lateral sides Contact torque (x, y, z components and magnitude) on medial and lateral sides Maximum and average pressure on medial and lateral sides Contact area Maximum interpenetration Output loads (equal and opposite for the two contacting bodies) 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

PAGE 31

22 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 concerned. 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 Matlabs 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.

PAGE 32

23 A B D C 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 grids. 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.

PAGE 33

24 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,

PAGE 34

25 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

PAGE 35

26 analysis. The articular cartilage and subchondral bone surfaces needed by the DLL were read in separately. The meniscus was omitted in the model. Youngs modulus for articular cartilage was set to 4 MPa (Cohen et al. 2003) and Poissons 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 femur. 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

PAGE 36

27 (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 mm. 002.0 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

PAGE 37

28 dynamic simulation. The contact surfaces were read in separately by the DLL. Youngs modulus for the ultrahigh molecular weight polyethylene tibial insert was set to 463 MPa (Kurtz et al. 2002) and Poissons 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

PAGE 38

29 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. 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. A B Figure 2-7. Dynamic simulation results. Contact forces A) and torques B) during a gait cycle of Osteonics7000 implant design.

PAGE 39

30 Figure 2-8. Contact areas during the contact phase of a gait cycle. 2.4 Discussion 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.,

PAGE 40

31 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. 2004). 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

PAGE 41

32 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 significant.

PAGE 42

33 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

PAGE 43

34 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 resources. 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.

PAGE 44

CHAPTER 3 EXPERIMENTAL EVALUATION OF AN ELASTIC FOUNDATION MODEL TO PREDICT CONTACT PRESSURES IN KNEE REPLACEMENTS 3.1 Introduction 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, DLima et al. 2001, Otto et al. 2001, Pri 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 35

PAGE 45

36 al. 1991, Li et al. 1997, Nuo 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 ultimately wear.

PAGE 46

37 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 for any spring can be calculated from (Johnson 1985, An et al. 1990, Blankevoort et al. 1991) p dhpEp)21)(1()()1( (3-1) where is Youngs modulus of the elastic layer, which can be a nonlinear function of )(pE p is Poissons ratio of the elastic layer, is the layer thickness at the spring location, h

PAGE 47

38 and 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). d A B C D 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 Youngs modulus-stress data reported by Cripton (1993) at 23C. 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.

PAGE 48

39 Given the known value for the deflection d of any spring at each instant in time, the contact pressure p for the spring can be easily calculated. For a nonlinear material model, an equation for as a function of (Cripton 1993) can be substituted into Eq. 3-1 to produce one nonlinear equation for (Nuo 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 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 pressures. E pp p 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: noooo2121 (3-2) where is strain, is stress, and o o and 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 n

PAGE 49

40 studies, and that choosing n = 1 produces a standard linear material model. For a given value of the current value of ddE for any spring can be found from 1121/1nooonE (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 23C with 0 = 0.0257 and 0 = 15.8; R2 = 0.925 at 37C with 0 = 0.0597 and 0 = 18.4). For comparison purposes, both linear (n =1) and nonlinear (n = 3) UHMWPE material models were used in our study. For both models, 0 was fixed and 0 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 0 = 1, 0 = 400, corresponding to a constant Youngs 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 0 = 0.0257, 0 = 15.9, very close to material parameters determined from Criptons data at 23C. Poissons 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

PAGE 50

41 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

PAGE 51

42 0.35 mm and maximum axial rotation was 1.5, well within the expected experimental variations. 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-defined tolerance. 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

PAGE 52

43 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 (1e-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 MOTIONs 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 90 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.

PAGE 53

44 A B 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 90) 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 23C, 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

PAGE 54

45 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

PAGE 55

46 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 90

PAGE 56

47 where the contact area is smaller, these errors were deemed to be sufficiently small not to hinder comparison with model predictions. 3.3 Results 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 material models. Errors (MPa) Linear Nonlinear Mean -2.2 -6.7 SD 1.2 3.6 Peak pressure without averaging RMS 2.5 7.5 Mean 2.0 -2.5 SD 1.4 2.7 Peak pressure with averaging RMS 2.4 3.6 Mean 0.31 0.56 SD 0.49 1.1 Average pressure 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

PAGE 57

48 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 90, which was consistent with the experiments. A 0 30 60 90 B 0 30 60 90 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.

PAGE 58

49 3.4 Discussion 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 implant designs. 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, DLima 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 23C found that a linear material model with a Youngs 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 Youngs modulus of UHMWPE to be 463 MPa using a miniature

PAGE 59

50 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 Youngs modulus is needed to produce the same total deformation. This explains the best-fit Youngs 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

PAGE 60

51 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 Youngs 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 MOTIONs 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

PAGE 61

52 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).

PAGE 62

53 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 models 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.

PAGE 63

CHAPTER 4 PREDICTED SENSITIVITY OF KNEE IMPLANT WEAR TO INSERT THICKNESS AND BODY MASS 4.1 Introduction 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 54

PAGE 64

55 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.

PAGE 65

56 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 insert thickness. 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

PAGE 66

57 thickness of 6.8 mm. The patients 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 analysis. 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 Youngs modulus (463 MPa) and Poissons 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

PAGE 67

58 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 Archards 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).

PAGE 68

59 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. 2004). 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

PAGE 69

60 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 = a0 + a1 Insert Thickness + a2 Body Mass (4-1) where a0, a1 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. 4.3 Results 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

PAGE 70

61 (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 thickness.

PAGE 71

62 Figure 4-4. 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 = a0 + a1 Insert Thickness + a2 Body Mass. Medial side Prediction a0 a1 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 Lateral side Prediction a0 a1 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

PAGE 72

63 predicted quantities except for damage area (a1 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 (|a1 / a2| > 1). The a0 values were not close to zero, indicating that the fit was only reasonable in the region investigated. 4.4 Discussion 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

PAGE 73

64 (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 factors. 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

PAGE 74

65 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 reasonable. 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

PAGE 75

66 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

PAGE 76

67 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 load. By dividing the a2 column by the a1 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 predicted trends. 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

PAGE 77

68 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 level. 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.

PAGE 78

CHAPTER 5 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 69

PAGE 79

70 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.

PAGE 80

APPENDIX A MAXIMUM PRESSURES DURING SIMULATIONS Table A-1. Maximum and average pressures during a cycle of gait simulation. Maximum pressures (MPa) Medial thickness (mm) Lateral thickness (mm) Body mass (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) Medial thickness (mm) Lateral thickness (mm) Body mass (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 71

PAGE 81

APPENDIX B WEAR, CREEP AND DAMAGE DEPTH Table B-1. Wear, creep and damage depth in the TKR insert. Medial thickness (mm) Lateral thickness (mm) Body mass (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 75 0.50 0.46 0.43 0.41 0.39 0.57 0.51 0.46 0.43 0.40 100 0.61 0.56 0.53 0.50 0.48 0.66 0.59 0.54 0.50 0.47 125 0.70 0.65 0.61 0.58 0.56 0.75 0.67 0.61 0.57 0.54 Wear depth (mm) 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 75 0.58 0.52 0.45 0.40 0.38 0.33 0.29 0.26 0.25 0.23 100 0.67 0.56 0.51 0.48 0.46 0.38 0.34 0.31 0.29 0.27 125 0.74 0.63 0.59 0.56 0.54 0.43 0.38 0.35 0.33 0.31 Creep depth (mm) 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 75 1.07 0.98 0.88 0.81 0.77 0.89 0.80 0.73 0.67 0.63 100 1.28 1.12 1.04 0.98 0.94 1.04 0.93 0.85 0.79 0.74 125 1.44 1.29 1.20 1.15 1.10 1.17 1.05 0.96 0.90 0.84 Damage depth (mm) 150 1.57 1.45 1.37 1.30 1.24 1.30 1.16 1.07 0.99 0.94 72

PAGE 82

APPENDIX C DAMAGE AREA AND VOLUME LOST Table C-1. Damage area and volume lost in the TKR insert. Medial thickness (mm) Lateral thickness(mm) Body mass (kg) 6 8 10 12 14 6 8 10 12 14 50 345 359 370 383 391 296 310 322 333 341 75 368 384 397 406 417 321 334 346 354 359 100 386 403 415 425 436 336 348 359 370 379 125 401 417 429 442 455 349 361 373 385 393 Damage area (mm2) 150 413 431 446 460 472 360 372 384 398 405 50 94 90 87 85 84 42 40 39 38 37 75 132 127 124 120 118 59 56 55 53 52 100 170 163 158 155 152 75 72 70 68 67 125 205 197 192 188 185 90 87 84 83 82 Volume lost (mm3) 150 239 235 229 224 220 105 101 99 97 96 73

PAGE 83

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 and Physics 20, 276-290. 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 of Biomechanical Engineering 115, 366-373. 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. 74

PAGE 84

75 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 Orthopaedics 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 Orthopaedics 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 of Biomechanics 24, 1019-1031. 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 Orthopaedics and Related Research 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 of Biomechanics 23, 353-355. Chillag, K. J., and Barth, E. (1991) An analysis of polyethylene thickness in modular total knee components. Clinical Orthopaedics 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 of Applied Mechanics 38, 387-392.

PAGE 85

76 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. Queens 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 and Physics 25, 3-9. Davoodi, R. and Loeb, G. E. (2001) Conversion of SIMM to SIMULINK for faster development of musculoskeletal models. Proceedings of IFESS, Cleveland, Ohio, 282-284. 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 of Biomechanics 28, 1155-1166. DLima, D.D., Chen, P.C., and Colwell, C.W. (2001) Polyethylene contact stresses, articular congruity, and knee alignment. Clinical Orthopaedics 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, Calgary, Canada. 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 of Biomechanics 36, 1659-1668.

PAGE 86

77 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 of Biomechanics 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, 842-847. 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 of Biomechanics 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 of Biomechanics 33, 435-442. Guilak, F., Butler, D. L., and Goldstein, S. A. (2001) Tissue engineering, cells, scaffolds, and growth factors: functional tissue engineering. Clinical Orthopaedics and Related Research 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 Orthopaedics 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 of Biomechanics 32, 951-958. Hefzy, M. S. and Cooke, T. D. V. (1996) Review of knee models: 1996 update. Applied Mechanics Reviews 49, S187-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.

PAGE 87

78 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 Mathematics 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 Orthopaedics 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 Orthopaedics and Related Research 350, 209-220. Kurtz, S. M., Jewett, C. W., Bergstrm, 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-85. Lee, K. Y. and Pienkowski, D. (1998) Compressive creep characteristics of extruded ultrahigh molecular-weight polyethylene. Journal of Biomedical Material Research 39, 261-265.

PAGE 88

79 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 of Biomechanics 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 OConnor, J. J. (1999) Bone position estimation from skin marker co-ordinates using global optimization with joint constraints. Journal of Biomechanics 32, 129134. 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 Orleans, LA. 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 Francisco, CA. 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.

PAGE 89

80 Nuo, N. and Ahmed, A.M. (2001) Sagittal profile of the femoral condyles and its application to femorotibial contact analysis. Journal of Biomechanical Engineering 123, 18-26. 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-323. Otto, J. K., Callaghan, J. J., and Brown, T. D. (2001) Mobility and contact mechanics of a rotating platform total knee replacement. Clinical Orthopaedics 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-1008. Paul, B. and Hashemi, J. (1981) Contact pressures on closely conforming elastic bodies. Journal of Applied Mechanics 48, 543-548. Pri, 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 of Biomechanical Engineering 123, 599-606. 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-1194. Rawlinson, J. J. and Bartel, D. L. (2002) Flat medial-lateral conformity in total knee replacements does not minimize contact stress. Journal of Biomechanics 35, 27-34. Rhinoceros Version 1.0 users guide (2001). Chapter 28: About NURBS. Robert McNeel & Associates, Seattle, WA.

PAGE 90

81 Riegger-Krugh, C, Gerhart, T. N., Powers, W. R., and Hayes, W. C. (1998) Tibiofemoral Contact Pressures in Degenerative Joint Disease. Clinical Orthopaedics 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 of Biomechanics 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. O., dela Rosa, M., Favae, 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 Orthopaedics and Related 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), 54-59. 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. Osteoarthritis 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 of Applied Mechanics, 484-490.

PAGE 91

82 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, 225-231. 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 of Biomechanical 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-848. 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 of Arthroplasty 13, 428-437. Tetsworth, K. and Paley, D. (1994) Accuracy of correction of complex lower-extremity deformities by the Ilizarov method. Clinical Orthopaedics and Related Research 301, 102-110. Waldman, S. D. and Bryant, J. T. (1997) Dynamic contact stress and rolling resistance model for total knee arthroplasties. Journal of Biomechanical Engineering 119, 254-260. 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 of Biomechanics 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, 677-685. Wu., J. Z., Herzog, W., and Epstein, M. (1998) Effects of inserting a pressensor film into articular joints on the actual contact mechanics. Journal of Biomechanical 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-895. Zajac, F. E. (1993) Muscle coordination of movement: A perspective. Journal of Biomechanics 26 (suppl 1), 109-124.

PAGE 92

BIOGRAPHICAL SKETCH 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 masters 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. 83


Permanent Link: http://ufdc.ufl.edu/UFE0002364/00001

Material Information

Title: Dynamic Stimulation of Knee Joint Contact during Human Movement
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0002364:00001

Permanent Link: http://ufdc.ufl.edu/UFE0002364/00001

Material Information

Title: Dynamic Stimulation of Knee Joint Contact during Human Movement
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0002364:00001


This item has the following downloads:


Full Text












DYNAMIC SIMULATION OF KNEE JOINT CONTACT
DURING HUMAN MOVEMENT

















By

YANHONG BEI


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


2003

































This document is dedicated to my husband, Liu Jian.















ACKNOWLEDGMENTS

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

studies.

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

Page

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

CHAPTER

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

APPENDIX

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


Table page

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

Figure page

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
MOVEMENT

By

Yanhong Bei

December 2003

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.














CHAPTER 1
INTRODUCTION

1.1 Background

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.

1.2 Objective

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

artificial knee.

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

wear performance.

Chapter 5 summarizes the research with recommendations for future improvements

to the contact model and simulation efficiency.














CHAPTER 2
METHODOLOGY FOR DYNAMIC CONTACT SIMULATION

2.1 Introduction

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

possible.

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

possible.

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







13


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.

A
Fixed Surface



Tangent Plane Moving Surface
D A

B E
Fixed Surface

m r

Moving Surface
B
C E
Fixed Surface C



Moving Surface


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

whenever possible.

Table 2-1. CPU time tests of distance evaluation methods.


Tolerance statistics (mm) CPU time (micro-sec)
Geometry
Maximum Average STD Ray firing Min.distance


Single surface
Single surface 4,459,921 879,488
(50x50)
Single surface
nglesurface 0.0412 0.00448 0.0053 3,789,363 437,182
(18 x 14)
12-patch surface
12-patch surface 0.0486 0.00148 0.0019 2,175,888 2,461,545
(8x7)
24-patch surface
24-patch surface 0.0434 0.00202 0.0028 2,031,058 2,791,974
(5 x 6)
48-patch surface
48-patchsurface 0.0397 0.00163 0.0026 1,929,121 3,442,914
(4 x 5)
96-patch surface
96-patch surface 0.0276 0.00134 0.0020 1,908,546 3,577,557
(3 x5)

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:

(1- v)E
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)
S(l+v)(1 2v)

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.

2003):


1pl p
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.

2-3:


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

concerned.

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
format)
Kinematics

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

concerned.

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.












A B

1600 / j 38 -
1500 36
1400 34.
S1300 32
o 0
o 1200. 30 .
0 0
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


C
D

S23.5 17

S22.5 11 5
E 215. a- 13
S20,5
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
grids.

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

femur.











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).


3.0 MPa

2.5

0 2.0



0 1'0

0,5

0.0


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

15
0 _--- --- ---- 15 -
.10

-500
5

-1000 0 -/" -

-5 I

SFx -10. A TX
-115 \ /
-Fz Tz
-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.










30 MPa
25

15


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.

2.4 Discussion

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.

2004).

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

significant.









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

resources.

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.














CHAPTER 3
EXPERIMENTAL EVALUATION OF AN ELASTIC FOUNDATION MODEL TO
PREDICT CONTACT PRESSURES IN KNEE REPLACEMENTS

3.1 Introduction

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

ultimately wear.









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).

A B

1 ^1500-

2 1000 -- Experiment
Model
500-


0 10 20 30 40 50
m Stress (MPa)
C D
30-
0" Peak
2 2-
Average



0 20 40 60 80 100
Active Elements

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

pressures.

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:

1 1
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

Cn-1
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

variations.

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.









> Time

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-
defined tolerance.

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.










A












B
200-

E 150- 0.2 MPa



I 50-
o
0
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.

3.3 Results

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
material models.
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








48



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
Load (N)
B Experiments
A Linear material model
Nonlinear material model


1500 2250 3000
90






1500 2250 3000
Load (N)


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 Discussion

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

implant designs.

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.














CHAPTER 4
PREDICTED SENSITIVITY OF KNEE IMPLANT WEAR TO INSERT THICKNESS
AND BODY MASS

4.1 Introduction

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

insert thickness.

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

analysis.

caMa








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.

2004).

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.

4.3 Results

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.


0.6m


,0.4

0.2

0


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
thickness.














E


04


1Figure 4-4.



Figure 4-4.


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
Mass.
Medial side

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
Lateral side

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.

4.4 Discussion

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

factors.

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

reasonable.

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

load.

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

predicted trends.

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

level.

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.














CHAPTER 5
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.















APPENDIX A
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)
mass
(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)
mass
(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















APPENDIX B
WEAR, CREEP AND DAMAGE DEPTH

Table B-1. Wear, creep and damage depth in the TKR insert.
Body Medial thickness (mm) Lateral thickness (mm)
mass
(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
Damage
depth 100 1.28 1.12 1.04 0.98 0.94 1.04 0.93 0.85 0.79 0.74
(mm)
(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















APPENDIX C
DAMAGE AREA AND VOLUME LOST

Table C-l. Damage area and volume lost in the TKR insert.
Body Medial thickness (mm) Lateral thickness(mm)
mass
(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,
276-290.

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-
373.

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,
1019-1031.

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,
3-9.

Davoodi, R. and Loeb, G. E. (2001) Conversion of SIMM to SIMULINK for faster
development of musculoskeletal models. Proceedings oflFESS, Cleveland, Ohio,
282-284.

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
28, 1155-1166.

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,
Calgary, Canada.

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,
842-847.

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
33, 435-442.

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-
85.

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
32, 129-134.

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
Orleans, LA.

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
Francisco, CA.

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
123, 18-26.

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-
323.

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-
1008.

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
123, 599-606.

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-
1194.

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),
54-59.

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,
225-231.

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-
848.

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
301, 102-110.

Waldman, S. D. and Bryant, J. T. (1997) Dynamic contact stress and rolling resistance
model for total knee arthroplasties. Journal of Biomechanical Engineering 119,
254-260.

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,
677-685.

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-
895.

Zajac, F. E. (1993) Muscle coordination of movement: A perspective. Journal of
Biomechanics 26 supplyl 1), 109-124.















BIOGRAPHICAL SKETCH

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.