UFDC Home  myUFDC Home  Help 



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 YiChung 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 21. CPU tim e tests of distance evaluation m ethods...................................................... 17 22. Input parameters and output measures in the API between the contact model and dynam ic contact m odel .......................................................................... 21 31. Comparison of linear and nonlinear material model ......................... ...............47 41. Linear regressions to wear predictions. ....................................... ............... 62 A1. Maximum and average pressures during a cycle of gait simulation.....................71 B1. Wear, creep and damage depth in the TKR insert. .............................................72 C1. Damage area and volume lost in the TKR insert. .................................................73 LIST OF FIGURES Figure page 21. Fram work of contact model .............. ......... ............................ ..................8 22. Distance correction for two distancefinding methods ......................... .........13 23. Sensitivity of contact results to grid density.................................. ..................... 23 24. Visualization of example dynamic simulation for settling the femur onto the tibia.26 25. Contact area of natural knee at the static position .............................................. 27 26. Dynamic simulation of contact during a cycle of gait for Osteonics7000 implant design ..............................................................................2 9 27. D ynam ic sim ulation results. ............................................. ............................ 29 28. Contact areas during the contact phase of a gait cycle..................... ..............30 31. Overview of materials and methods. .................. .. ....................... ............ .. 38 32. Visualization of example dynamic simulation used to settle the femoral component onto the tibial insert .................. ...................................... .. ............ 42 33. Experimental setup and pressure measurement......................................................44 34. Comparison between experimental and predictions..........................................48 41. Overview of the computational framework for wear predictions.............................57 42. Comparison of retrieval and predictions. ............. ............ ................................59 43. The damage plots for five combinations of insert thickness and body mass. ..........61 44. Threedimensional 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 ultrahighmolecularweightpolyethylene (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 surfacesurface 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, RieggerKrugh 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 fulllimb or fullbody 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 subjectspecific models. The key problem that remains is how to develop subjectspecific 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 welldefined 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 ultrahigh materialweightpolyethylene (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. Subjectspecific 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, arthritisrelated 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 xray 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 noninvasive experimental approach does not exist for measuring in vivo joint loading. Consequently, computer simulations are needed to develop predictions. Since muscle cocontractions 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 twodimensional 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 threedimensional joint contact modeling. Wismans et al. (1980) were the first to develop a quasistatic threedimensional natural knee model using rigid body contact. Pandy et al. (1997) presented a quasistatic natural knee model using deformable contact with idealized (i.e., planes and polynomials) surfaces. AbdelRahman 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 fullbody dynamic simulation. Blankevoort et al. (1991), Cohen et al. (2001 and 2003), Kwak et al. (2003), and Elias et al. (2003) reported quasistatic natural knee models with deformable contact, where subjectspecific 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 quasistatic 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, threedimensional 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 fourstep 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. 21). 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 21. 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 surfacesurface 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 Bspline (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 threedimensional NURBS patch is described mathematically by a rectangular, twodimensional (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 uv 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 uv 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 uv space and reduce the number of uv 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 uv 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 preprocessing phase prior to beginning the dynamic simulation. In this way, many CPUintensive 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 twodimensional 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 uv 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 uv 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 22. Distance correction for two distancefinding 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 halfspace 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. 22A). The tangent plane is perpendicular to the unique line defined by where the normals on the two contact surfaces are colinear and antiparallel. 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. 22B) or conformal (Figure. 22C). This is achieved by using trigonometry to make the minimum distance vector (m in Figure. 22) slightly larger or the ray firing distance vector (r in Figure. 22) slightly smaller (Figure. 22C). The midsurface normal can be well approximated by a line (ADF in Figure. 22C) bisecting the minimum distance (AB in Figure. 22C) and ray firing (AC in Figure. 22C) 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 uv 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 userdefined 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 BreadthFirst 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 nonactive 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 BreadthFirst 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 21). 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 uv 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 uv lines to 96 patches with sparse uv 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 uv lines whenever possible. Table 21. CPU time tests of distance evaluation methods. Tolerance statistics (mm) CPU time (microsec) 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) 12patch surface 12patch surface 0.0486 0.00148 0.0019 2,175,888 2,461,545 (8x7) 24patch surface 24patch surface 0.0434 0.00202 0.0028 2,031,058 2,791,974 (5 x 6) 48patch surface 48patchsurface 0.0397 0.00163 0.0026 1,929,121 3,442,914 (4 x 5) 96patch surface 96patch 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 rayfiring or minimum distance from a point to a single surface. For multipatch surfaces, the rayfiring 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 halfspace theory, where the contacting bodies are semiinfinite 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 halfspace 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. 21 is used to calculate the pressure of any spring: (1 v)E p = d (21) (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 elementbyelement 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. 21 becomes: (1 v)E ( (22) 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 threeparameter nonlinear material model (Fregly et al. 2003): 1pl p E= E + (23) 2 po 2 po where E is strain, p is contact pressure, and so, Po, and n are material parameters. This is a nonlinear powerlaw 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. 23: E= 1/1 1+n (24) 2 po, Po Given the interpenetration dfor any spring, Eq. 24 can be substituted into Eq. 21 or Eq. 22 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 22), 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 bodythree 123 sequence (Kane et al. 1983). Table 22. 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 22). 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 userdefined 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 23. 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. 23). 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 degreeof 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.0T 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 semiautomatically 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 uv 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 timeframe response. Minimum distance was used for the distance calculations, the contact surface uv 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 24. 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 userdefined 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 flexionextension, internalexternal rotation, and anteriorposterior 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. 24). 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. 25). 3.0 MPa 2.5 0 2.0 0 1'0 0,5 0.0 Figure 25. 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 cruciateretaining 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 uv control points for contact surfaces was minimized in Rhinoceros, with the tolerance between the original and resurfaced 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 uv 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, flexionextension, internalexternal rotation, and anteriorposterior 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 offcenter to produce a 70% medial30% 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). Superiorinferior translation, mediallateral translation, and variousvalgus 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. 26), the contact forces in the medial and lateral compartments, and the net contact force and torque applied to both bodies (Figure. 27). 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. 28). 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 26. 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 27. 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 28. 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 halfspace 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, forcegenerating 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 elementbyelement 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 Kscan sensor. An Exatech Optetrak B posteriorstablized 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 postmortem 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 quasistatic, 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, nonhomogeneity 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 contactrejection 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 largescale 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 fullbody 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 postmortem, 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 postmortem 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 threedimensional 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 threedimensional 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 halfspace contact, where both bodies are elastic and semiinfinite. 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 (31) p d (31) (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 superiorinferior 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 31A). 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 31. 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 modulusstress 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. 31 to produce one nonlinear equation for p (Nufio and Ahmed 2001). Since each spring is independent of its neighbors, any standard nonlinear rootfinding algorithm can solve the resulting system of nonlinear pressure equations independently. For a linear material model, Eq. 31 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 nonzero 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 + (32) 2 ao 2 Uo where c is strain, a is stress, and co, uo, and n are material parameters. This is the nonlinear powerlaw 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 Cn1 E=11/ 1+nl o (33) Equation 33 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 31B; 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 mediallateral symmetry and a 9 mm minimum insert thickness (Figure 31C; 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 setup 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 degreeoffreedom (DOF) joint. Sagittal plane translations and varusvalgus 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 mediallateral 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 mediallateral 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 31D). Thus, 100 active elements per side were more than accurate enough for comparison with experimental data. > Time Figure 32. 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 32). 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 (le5 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 IIIM 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 builtin 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 33A). 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 varusvalgus 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 mediallateral 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 varusvalgus 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 33. Experimental setup and pressure measurement. A) Experimental contact pressure measurements made on the same implant using a Tekscan KScan pressure measuring system and a servohydraulic test machine. B) Sensitivity of experimental contact area to the pressure cutoff 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 KScan sensor has been shown to exhibit significant calibration drift over time (Otto et al. 1999), each experimental trial was selfcalibrated. Ramp loads from 20% to 100% of desired load were applied over 5 seconds, and these two points were used to perform the twopoint calibration procedure recommended by the manufacturer. The two loads required for postcalibration were measured during each trial by an inline 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 totrial calibration drift was eliminated. The second issue was determination of an appropriate pressure cutoff 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 cutoff values ranging from 0 to 2 MPa (Figure 33B). While changing the cutoff 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 cutoff 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 cutoff. 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 builtin 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 KScan 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 nondimensional 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 36% 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 31. Comparison of linear and nonlinear material model. Mean, standard deviation (SD), and rootmeansquare (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 31). Both models tracked the average contact pressure well, though the linear model tracked it better (Figure 34A). 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 34B). 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 34. 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 rigiddeformable 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 halfspace, 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 bestfit 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 elasticplastic 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 nonplanar contact patches. Both of these features are limitations in elastic halfspace 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 halfspace 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 posteriorstabilized 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 builtin elastic contact model uses either a Hertz formulation for approximately quadratic surfaces or an elastic halfspace 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 halfspace 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 quasistatic 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 subsurface stresses. While the model would be useful for wear predictions made from contact pressure and kinematic inputs, predictions involving subsurface 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 byelement 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 34B). 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 elasticplastic 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 ultrahighmolecularweightpolyethylene (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 welldefined 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 threestep process. Singleplane fluoroscopy was used to measure in vivo knee implant kinematics during treadmill gait and stair activities (Figure 41A). 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 cruciateretaining 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 41. Overview of the computational framework for wear predictions. A) In vivo fluoroscopic data used to provide patientspecific 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 elementbyelement based on interpenetration of the undeformed contact surfaces. The contact model was integrated into a multibody dynamic simulation of the implant components (Figure 41B) 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. Flexionextension, internalexternal rotation, and anteriorposterior translation were prescribed from the fluoroscopic kinematics, while superiorinferior translation, varusvalgus rotation, and mediallateral 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 offcenter 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 41C) was used to predict wear, creep, and damage (= wear + creep) on an elementbyelement 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 gaitstair activity were generated using a linear rulesof mixture. The entire modeling and simulation process was evaluated by comparing damage predictions from the model (Figure 42A) with damage observed on the tibial insert retrieved postmortem from the same patient who provided the in vivo kinematics (Figure 42B). Further details on the damage predictions can be found in (Fregly et al. 2004). A B C B 0.90 Figure 42. 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 42). 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 (41) 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 43. 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 44. Figure 44. Threedimensional 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 41. 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 41), indicating a strong bilinear relationship between insert thickness, body mass, and all predicted quantities in the region investigated (Figure 44). 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 damagerelated 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 tradeoff 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 damagerelated 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 tradeoff 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 camfollower 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 ligamentretaining 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 damagerelated 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 44) 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 68 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 41, 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 12 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 cocontraction 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 tradeoffs 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 nonlinear 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 A1. 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 B1. 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 Cl. 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 AbdelRahman, E. M. and Hefzy, M.S. (1998) Threedimensional dynamic behavior of the human knee joint under impact loading. Medical Engineering andPhysics 20, 276290. Ahmadi, N., Keer, L. M., and Mura, T. (1983) NonHertzian contact stress analysis for an elastic halfspace normal and sliding contact. International Journal of Solids and Structures 19, 357373. 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, 10131020. 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, 201231. Andriacchi, T. P. (1994) Dynamics of knee malalignment. Orthopaedic Clinics of North America 25, 395403. Archard, J. F. and Hirst, W. (1956) The wear of metals under unlubricated conditions. Proceedings of the Royal Society A236, 397410. Ateshian, G. A. (1993) A Bspline leastsquares surfacefitting method for articular surfaces of diarthrodial joints. Journal ofBiomechanical Engineering 115, 366 373. Banks, S. A. and Hodge, W. A. (1996) Accurate measurement of threedimensional knee replacement kinematics using singleplane fluoroscopy. IEEE Transactions on Biomedical Engineering 43, 638649. Bartel, D. L., Bicknell, V. L., and Wright, T. M (1986) The effect of conformity, thickness, and material on stresses in ultrahigh molecular weight components for total joint replacement. Bone and Joint Surgery 68A, 10411051. Bartel, D. L., Burstein, A. H. and Edwards, D. L. (1985) The effect of conformity and plastic thickness on contact stress in metalbacked plastic implants. Journal of biomechanical Engineering 107, 193199. 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, 7682. Bendjaballah, M. Z., ShiraziAdl, A., and Zukor, D. J. (1997) Finite element analysis of human knee joint in varusvalgus. Clinical Biomechanics 12, 139148. Benjamin, J., Szivek, J., Dersam, G., Persselin,S. and Johnson, R. (2001) Linear and volumetric wear of tibial inserts in posterior cruciateretaining knee arthroplasties. Clinical 0i Ithp q,,, li 1 and Related Research 392, 131138. Blankevoort, L., Kuiper, J. H., Huiskes, R., and Grrotenbeor, H. J. (1991) Articular contact in a threedimensional model of the knee. Journal ofBiomechanics 24, 10191031. 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, 253260. Brenan, K. E., Campbell, S. L., and Petzold, L. R. (1996) Numerical Solution of Initial Value Problems in DifferentialAlgebraic Equations. SIAM Classics in Applied Mechanics, 2nd edition, Philadelphia. Chao, E.Y. and Sim, F.H. (1995) Computeraided preoperative planning in knee osteotomy. Iowa Orthopedic Journal 15, 418. Cheng, R. C.K., Brown, T. D., and Andrews, J. G. (1990) Nonuniqueness of the bicompartmental contact force solution in a lumped parameter mathematical model of the knee. Journal ofBiomechanics 23, 353355. 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, 261263. Cohen, Z. A, Henry, J. H, McCarthy, D. M, Mow, V. C, and Ateshian, G. A. (2003) Computer simulations of patellofemoral joint surgery. Patientspecific models for tuberosity transfer. American Journal of Sports Medicine, 31, 8798 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, 480487. Conry, T. F. and Seirig, A. (1971) A mathematical programming method for design of elastic bodies in contact. Journal ofApplied Mechanics 38, 387392. 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, 95106. Cripton, P.A. (1993) Compressive characterization of ultrahigh 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, 39. Davoodi, R. and Loeb, G. E. (2001) Conversion of SIMM to SIMULINK for faster development of musculoskeletal models. Proceedings oflFESS, Cleveland, Ohio, 282284. Delp, S.L., Arnold, A. S., and Piazza, S. J. (1998) Graphicsbased modeling and analysis of gait abnormalities. Biomedical Materials and Engineering 8, 227240. 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, 11551166. 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, 232238. 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 tibiofemoral contact. Journal of Biomechanical Engineering 124, 273280. 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 threedimensional compliant contact model for dynamic simulation of total knee replacements. In Proceedings of the VIIth International Symposium on Computer Simulation in Biomechanics, 1013. August 57, 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, 16591668. 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 singleleg standingeffect of gender and anatomic parameters. Journal ofBiomechanics 34, 895905. 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, 842847. Glitsch, U. and Baumann, W. (1997) The threedimensional determination of internal loads in the lower extreme. Journal of Biomechanics, 30, 11231131. 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, 267275. 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, 435442. 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, S295S305. 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 2529, Key Blascayne, Florida, 551552. 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, 383393. 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, 951958. Hefzy, M. S. and Cooke, T. D. V. (1996) Review of knee models: 1996 update. Applied Mechanics Reviews 49, S 187193. Hoff, W. A., Komistek, R.D., Dennis, D.A., Gabriel, S. M., and Walker, S. A. (1998) Threedimensional determination of femoraltibial contact positions under in vivo conditions using fluoroscopy. Clinical Biomechanics 13, 455472. 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, 423430. 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), 18. 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 nonHertzian halfspace contact problems. Journal of Engineering M1A aheiuh tii \ 6, 193206. 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, 1622. 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, 6981. Kurtz, S. M., Bartel, D. L., and Rimnac, C. M. (1998) Postirradiation aging affects the stresses and strains in UHMWPE components for total joint replacement, Clinical O 1 /i1'(t hik % and Related Research 350, 209220. 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, 19071919. Kwak, S. D., Blankevoort, L., and Ateshian, G. A. (2000). A mathematical formulation for 3D quasistatic multibody models of diarthrodial joints. Computational Methods in Biomechanical and Biomedical Engineering 3, 4164. Landy, M., and Walker, P. (1998) Wear of ultrahighmolecularweight 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 molecularweight polyethylene. Journal of Biomedical Material Research 39, 261265. Lewis, G. (1998) Contact stress at articular surfaces in total joint replacements. Part I: Experimental methods. Biomedical Materials and Engineering 8, 91110. 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, 395400. 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, 635638. 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, 717721. 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, 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, 416421. 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, 3539. McGloughlin, T. M. and Kavanagh, A. G. (2000) Wear of ultrahigh molecular weight polyethylene (UHMWPE) in total knee protheses: a review of key influences. Proceedings of the Institute of Mechanical Engineers 214, 349359. 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 heeltoe running. Computer Methods in Biomechanics and Biomedical Engineering 3, 321334. Nuio, N. and Ahmed, A.M. (2001) Sagittal profile of the femoral condyles and its application to femorotibial contact analysis. Journal ofBiomechanical Engineering 123, 1826. Otto, J. K., Brown, T. D., and Callaghan, J. J. (1999) Static and dynamic response of a multiplexedarray 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, 2437. Pandy, M. G., Sasaki, K., and Kim, S. (1997) A Threedimensional musculoskeletal model of the human knee joint. Part 1: Theoretical Construction. Computer Methods in Biomechanics and Biomedical Engineering 1, 87108. 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, 543548. Perie, D. and Hobatho, M. C. (1998) In vivo determination of contact areas and pressure of the femorotibial joint using nonlinear finite element analysis. Clinical Biomechanics 13, 394402. Piazza, S. J. and Delp, S. L. (1999) Threedimensional simulation of total knee replacement motion during stepup tasks. In Proceedings of the 1999 Bioengineering Conference, BEDVol. 42, 423424. The American Society of Mechanical Engineers, New York. Piazza, S. J. and Delp, S. L. (2001) Threedimensional dynamic simulation of total knee replacement motion during a stepup task. Journal ofBiomechanical Engineering 123, 599606. 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 mediallateral conformity in total knee replacements does not minimize contact stress. Journal ofBiomechanics 35, 2734. Rhinoceros Version 1.0 user's guide (2001). Chapter 28: About NURBS. Robert McNeel & Associates, Seattle, WA. RieggerKrugh, 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, 233245. 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, 564571. Sathasivam, S. and Walker, P.S. (1999) The conflicting requirements of laxity and conformity in total knee replacements. Journal ofBiomechanics 32, 239247. 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 ScotchYoke mechanism. Journal of Tribology (Accepted). Schiehlen, W. (1997) Multibody system dynamics: roots and perspectives. Multibody System Dynamics 1, 149188. Schipplein, O. D. and Andriacchi, T. P. (1991) Interaction between active and passive knee stabilizers during level walking. Journal of Orthopaedic Research 9, 113119. 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, 3646. 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 80A (1), 5459. 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, 214. 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, 703708. Singh, K. P. and Paul, B. (1974) Numerical solution of nonHertzian elastic contact problems. Journal ofApplied Mechanics, 484490. 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, 225231. Tashman, S. and Anderst, W. (2003) Invivo measurement of dynamic joint motion using high speed biplane radiography and CTL: application to canine ACL deficiency. Journal ofBiomechanical Engineering 125, 238245. 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, 428437. Tetsworth, K. and Paley, D. (1994) Accuracy of correction of complex lowerextremity deformities by the Ilizarov method. Clinical 01 thIiope i~ % and Related Research 301, 102110. Waldman, S. D. and Bryant, J. T. (1997) Dynamic contact stress and rolling resistance model for total knee arthroplasties. Journal of Biomechanical Engineering 119, 254260. 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, 8389. 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, 677685. 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, 655659. 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), 109124. 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. 