UFDC Home  myUFDC Home  Help 



Full Text  
EXPERIMENTAL EVALUATION OF A NATURAL KNEE CONTACT MODEL USINTG RESPONSE SURFACE OPTIMIZATION By YICHIUNG LIN A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2004 Copyright 2004 by YiChung Lin ACKNOWLEDGMENTS First, I would like to sincerely thank Dr. Benj amin J. Fregly, my advisor, for the direction he has given me throughout my research. His encouragement and dedication have greatly contributed to the accomplishment of this thesis. I would also like to thank Dr. Raphael T. Haftka and Dr. Andrew Kurdila for being my committee members. It was a great honor to have them both. I would also like to thank all the Computational Biomechanics Lab members for their support in my past and present work, with special thanks to Dong, Yanhong, Jaco and Jeff for their assistance and suggestion. TABLE OF CONTENTS Page ACKNOWLEDGMENT S .........__.. ..... .__. .............._ iii.. LIST OF TABLES ........._.___..... .__. ...............v.... LIST OF FIGURES .............. ....................vi AB STRAC T ................ .............. vii CHAPTER 1 INTRODUCTION ................. ...............1.......... ...... 1.1 Need for Accurate Contact Model ................ ...............1.............. 1.2 Need for Efficient Contact Model Evaluation ................. .......... ................2 1.3 Approach............... ...............3 2 M ETHODS .............. ...............5..... 2.1 Response Surface Optimization............... .............. 2.2 Contact Pressure Experiments ................ ...............10................ 2.3 Knee Model Creation............... ...............12 2.4 Knee Model Evaluation ................ ...............16................ 3 RE SULT S ................. ...............18.......... ..... 4 DI SCUS SSION ................. ...............23................ 5 SUMMARY AND FUTURE STUDY ................. ...............28........... ... 5.1 Summary ................. ...............28................ 5.2 Future Study ................. ...............28........... ... LIST OF REFERENCES ................. ...............30........... .... BIOGRAPHICAL SKETCH .............. ...............35.... 1V LIST OF TABLES Table pg 21. The averaged experimental measurements were collected from three compression trials for both loads processing on the same cadaver knee via servohydraulic test m machine .............. ...............12.... 31. Comparison of response surface predictions to data points sampled from large and small strain contact model s. ............. ...............20..... LIST OF FIGURES Figure pg 21. A human cadaver knee static experiment setup .........._._ ....... ..._ .............._.11 22. Original and segmented medical images ..........._..__......_. ......._._..........1 23. The anterior and posterior views of the A) 3D point clouds and B) 3D NURBS model with six screws created from CT and MRI images .............. ....................14 31. Percent errors between response surfaces and predicted results from both contact m odels .............. ...............21.... 32. Percent errors between experimental and predicted results from both contact m odels.. ............ ...............22..... Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science EXPERIMENTAL EVALUATION OF THE NATURAL KNEE CONTACT MODEL USINTG RESPONSE SURFACE OPTIMIZATION By YiChung Lin August 2004 Chair: Benjamin Jon Fregly Major Department: Mechanical and Aerospace Engineering Finite element, boundary element, and discrete element models have been employed to predict contact conditions in human j points. When optimization is used to evaluate the ability of such models to reproduce experimental measurements, the high computational cost of repeated contact analysis can be a limiting factor. This thesis presents a computationallyefficient response surface optimization methodology to address this limitation. Quadratic response surfaces are fit to contact quantities (peak pressure, average pressure, contact area, and contact force) predicted by a joint contact model for various combinations of material modulus and relative pose (i.e., position and orientation) of the contacting bodies. The response surfaces are used as surrogates for costly contact analyses in an optimization that minimizes differences between measured and predicted contact quantities. The approach is demonstrated by evaluating a linear elastic discrete element contact model of the tibiofemoral j oint, where the model was created using CT and MRI data from the same cadaveric specimen used in static pressure experiments. For variations in material modulus and relative bone pose within the envelope of experimental uncertainty (+ 1 mm and lo), quadratic response surfaces accurately predicted contact quantities computed by the discrete element model. Using these response surfaces, 500 optimizations with different initial guesses were performed in less than 90 seconds. For a flexion angle of 300 and axial loads of 500 and 1000 N, the optimizations demonstrated that small and large strain versions of the contact model could match all experimentally measured contact quantities to within 10% error with the exception of peak contact pressure, which was in error by as much as 85%. Thus, discrete element models of natural j points may be best suited for predicting contact quantities that involve averaging across the surface but not quantities associated with specific locations on the surface. CHAPTER 1 INTTRODUCTION 1.1 Need for Accurate Contact Model According to the Arthritis Foundation, there were nearly 43 million Americans with arthritis or chronic j oint symptoms in 1998. The number went up to 70 million in 2001 and most likely will keep climbing due to the rising number of aging Baby boomers. This foundation also reported that arthritis limits daily activities such as walking, running and dressing for more than seven million Americans. There are many different types of arthritis including osteoarthritis, rheumatoid arthritis, gout, and juvenile arthritis. Among them, osteoarthritis is the most prevalent form of arthritis affecting more than 20.7 million Americans. Osteoarthritis, or degenerative joint disease, is multifactorial with genetic, biologic, and mechanical factors all playing a role. Of the mechanical factors involved, contact pressure within the joint has been shown to have an interactive effect on developing this disease (Hasler et al., 1998; Herzog et al., 2003). Thus, knowledge of in vivo contact forces and pressures in human joints would be valuable for improving the prevention and treatment of joint arthritis. Although dynamic imaging advances now is capable collecting accurate measurement of in vivo j oint kinematics (Komistek et al., 2003; Tashman and Anderst, 2003), joint contact forces and pressures are difficult to measure in vivo (Kaufman et al., 1996), therefore, necessitating modelbased analyses to develop predictions as well as static in vitro testing to evaluate these predictions. A variety of joint contact modeling methods have been used for this purpose, including finite element (Bendjaballah et al., 1995, 1997, 1998; Donzelli and Spilker, 1998; Donahue et al., 2002; Stolk et al., 2002), boundary element (Kuo and Keer, 1993; Haider and Guilak, 2000), and discrete element methods (Li et al., 1999, 2001; Pandy and Sasaki, 1998; Piazza and Delp, 2001; Dhaher and Kahn, 2002). 1.2 Need for Efficient Contact Model Evaluation Once a joint contact model has been created that represents an in vitro testing situation, its ability to reproduce experimentally measured contact quantities (e.g., peak pressure, average pressure, contact area, and contact force) must be evaluated before any further application. At least two factors complicate the process of evaluation. The first is uncertainties in the experimental measurements. These uncertainties can often be estimated and involve quantities such as the position and orientation (i.e., pose) of cadaveric bones measured by the test apparatus, contact pressures and areas recorded by a pressure sensor, and the articular surface geometry determined from medical imaging data. Unknown model parameters, such as material parameters in the contact model, present additional sources of uncertainty. A second complicating factor is the high computational cost of repeated contact analysis. Given an estimated envelope of uncertainty, optimization methods can be used to determine if a feasible combination of model parameters could be used to reproduce all experimental measurements simultaneously (Fregly et al., 2003). For example, an optimization can vary model material parameters and the relative pose of the contacting bodies within the envelope of uncertainty until the model produces the best match to the experimental contact data. The problem with this approach is that the high computational cost of repeated contact analysis can make such optimizations extremely time consuming and in some cases even impractical. Response surface (RS) methods have been utilized successfully in other situations to eliminate computational bottlenecks in optimization studies. Response surfaces are simply multidimensional linear regression curve fits to quantities of interest predicted by an engineering model. Once the mathematical form of the RS is specified, linear least squares is typically used to determine the coefficients that provide the best fit to each predicted quantity of interest (i.e., each response) as a function of the specified design variables. These surface approximations are then used as surrogates for costly engineering analyses when the optimization is performed. Outside of the biomechanics community, RS optimization methods have been used for structural design applications (Jansson et al., 2003; Liu et al., 2000; Rikards et al., 2004; Roux et al., 1998), aerodynamic designs (Ahn and Kim, 2003; Papila et al., 2002; Sevant et al., 2000), and fluid dynamics (Burman and Gebart, 2001; Keane, 2003; Leary et al., 2004). Within the biomechanics community, little work has been performed using RS optimization methods with the exception of recent studies by Jung and Choe (1996), Chang et al. (1999), and Hong et al. (2001). To our knowledge, no studies in the literature have used RS methods to perform optimizations of contact problems. 1.3 Approach The two goals of this thesis are first, to develop a computationally efficient RS optimization approach for evaluating a joint contact model's ability to reproduce static experimental contact measurements, and second, to apply this approach to the evaluation of a discrete element contact model of the tibiofemoral joint. Our specific hypotheses were that 1) quadratic RSs can accurately predict contact quantities (peak pressure, average pressure, contact area, and contact forces) computed by a discrete element contact model for small variations material modulus and relative pose of the contacting bodies, and 2) a discrete element contact model of the tibiofemoral j oint can reproduce experimentally measured contact quantities that involve averaging across the surface (i.e., average pressure, contact area, contact force) but not quantities associated with specific locations on the surface (i.e., peak pressure). Our study provides a new computationally efficient approach for contact model evaluation as well as a better understanding of the capabilities and limitations of discrete element contact models of natural human j points. CHAPTER 2 IVETHOD S 2.1 Response Surface Optimization This section provides a general overview of RS approximation methods as well as specific modifications required to apply these methods to contact analyses. The RS method can be defined as a collection of statistical and mathematical techniques useful for constructing smooth approximations to functions in a multidimensional design space. Once a mathematical form has been selected, the coefficients of the approximate function (responses surface) are determined using data from either physical experiments or numerical simulations. The most common mathematical form for a RS is a lowdegree polynomial. For example, a quadratic response surface with design variable inputs x, and x, and output y is formulated as y = Po + P x1 + P x2 + P x12 + Pqx2 + P x, x2 (1) where the p, (i = 0,...,5) are the unknown coefficients to be fitted. A low degree polynomial minimizes the number of unknown coefficients and tends to smooth out noise in the function. Response surface approximations work best when the number of design variable inputs is small (< 10), since a large number of design variables results in a complicated design space that is difficult to fit with lowdegree polynomials To develop RS approximations for contact problems, one must identify the design variable (DV) inputs, the outputs to be predicted, and the mathematical form of the RS relating them. For contact model evaluation with static experimental data, the design variables are the six relative pose parameters (i.e., three translations and three rotations  6 DVs) and material modulus (1 DV) of the contacting bodies. The experimental uncertainty of these quantities can be estimated and their values can be changed in the contact model. The RS outputs are peak pressure, average pressure, contact area, and contact force. These quantities can be calculated by the contact model and measured experimentally for comparison. The hypothesized mathematical form is a quadratic RS with one modification. For linearly elastic Hertzian point contact, the peak pressure, average pressure, contact force, and contact area are all functions of interpenetration (vertical translation) to a power less than two, while the material modulus (assumed to be the same for both bodies) linearly scales each quantity except for area (Johnson, 1985). Thus, data for the RSs are generated using a material modulus of one, only the six pose parameters are used as RS inputs, and the RS outputs (except area) are scaled by the desired modulus value. With the RS formulation specified, the next step is to determine a sampling scheme within the design space to provide data for fitting the RS. Since this sampling process is only preformed once to generate the RSs, the computational cost of the contact analyses is paid only once up front. A quadratic response surface using k design variables will possess p = (k +1l)(k + 2) /2 unknown coefficients, where k = 6 since only pose parameters are used as RS inputs. Consequently, a minimum of p = 28 data points must be sampled to perform the linear leastsquares fit. However, to cover the design space in a systematic manner, we select a larger number of sample points using design of experiments (DOE) theory. Several DOE sample criteria are available, including the factorial design, face centered central composite design (FCCCD), and the Doptimality design. We choose the FCCCD criteria for its ability to sample all regions of the design space. For a quadratic RS, this approach utilizes 2k + 2k +1 = 77 sample points, where the samples are taken at the center, the corners, and the face centers of a k dimensional hypercube. For contact analyses, we make two modifications to the FCCCD sampling scheme to improve the quality of the fit. The first modification accounts for infeasible points. A sampled point is deemed to be infeasible and is therefore omitted if the contact force and area predicted by the contact model are zero. This modification avoids fitting regions of the design space where no contact is occurring. The second modification accounts for outlier points. Once a RS is generated from feasible points, the RS output is compared to the computed value from the contact model for every sample point. The point with the largest absolute percent error above a preselected cutoff value of 10% (a typical value for engineering analyses) is omitted and the RS regenerated from the remaining sample points. The procedure is iterated until all sample points are below 10% error. This modification avoids fitting regions of the design space where only light contact is occurring, thereby providing a better fit in the regions of interest where the contact force is large. Omission of several sample points does not pose a problem to the fitting process since the FCCCD sampling scheme is highly over determined. After a RS is generated, the quality of the resulting fit must be assessed, since a poor quality fit indicates that a different mathematical form for the RS should be considered. We use three common error measures for this purpose. All of these measures make use of the sum of the squares of the errors SSE between predicted responses A from the RS and actual responses y, computed by the contact model, where n (n >> 28 and n < 77) is the number of sample points used to generate the RS: The first measure of fit quality is the adjusted rootmeansquare error ( RMSEa4 ). Given the SSE, the rootmean square error (RM~SE) can be calculated from RM~SE = (3) However, this measure will be zero if n = p (i.e., no redundant points), even though the errors would not necessarily be zero at nonsampled points. To address this limitation, we choose a more conservative adjusted RM~SE that uses n p (i.e., the number of degrees of freedom remaining in the fitting process) rather than n in the denominator of Eq. (3): ISSE RM~SE,, =,1 (4) Vnp To provide a relative measure of fit quality, we also compute the percent adjusted RM\~SE usmng 100 rSSE %RM~SEad =1 (5) where y~ represents the magnitude of the fitted quantity: 1 "6 The second measure of fit quality is the adjusted coefficient of determination (R2 aq ). The coefficient of determination R2 SUfferS from a, similar problem to RMSE in that a perfect fit will be indicated if n = p Consequently, we used the adjusted R2 value to account for the degrees of freedom n p remaining in the fit: SSE / (n p) RZ = 1 (7) where ]7 is the mean of the actual responses. The Einal measure of fit quality is the RM~SE calculated from the prediction error sum of squares (PRESS) statistic. To evaluate the predictive capability of a RS, one should ideally sample additional points distinct from those used to generate the RS. However, this approach would require a significant number of additional costly contact analyses. To circumvent this issue, the PRESS analysis excludes one sample point at a time from the set used to generate each RS. The RS is regenerated using the remaining n 1 sample points and the prediction error at the omitted sample point calculated. This process is repeated for all n sample points, and the resulting errors are used to compute a PRESSbased SSE called the PRESS statistic. From there, a PRESSbased RM~SE can be calculated from Eq. (3), where n rather than n p is used in the denominator since each error is calculated from a RS that omits that point. Once accurate RSs are generated for the output quantities of interest, they are used in an optimization to evaluate the contact model's ability to reproduce experimental measurements. Each time the optimization requires a peak pressure, average pressure, contact force, or contact area from the contact model, a response surface is used in place of a contact analysis to provide the value. By fitting quantities computed by the contact model, one can create any cost function that can be built up from the basic quantities. If the cost function was fitted directly using its own response surface, then additional contact analyses would be required to generate a new response surface each time the cost function was modified. With our approach, a wide variety of cost functions can be constructed without the need for any additional contact analyses. 2.2 Contact Pressure Experiments The response surface methodology described above was used to evaluate a natural knee contact model's ability to reproduce experimental contact measurements. The experiments were performed on a single cadaveric knee specimen cut approximately 15 cm above and below the joint line and showing no visible signs of degenerative joint disease. Institutional review board approval was obtained for the testing and subsequent modeling efforts. The menisci, fibula, and patella were removed from the specimen, and three titanium bone screws were inserted into the tibia and femur as landmarks for contact model alignment. The tibia and femur were potted in neutral alignment (MacWillams et al., 1998) and mounted in a MTS MiniBionix 858 servohydraulic test machine. The position and orientation of the femur were constrained using custom fixturing that allowed adjustment of the sagittal plane rotation and mediallateral translation relative to the ram of the test machine. The axial plane position and orientation of the tibia were unconstrained using a ball plate, thereby allowing the tibia to selfalign with the femur once an axial load was applied. Using this setup (Fig. 21), we collected four experimental quantities of interest from the medial and lateral compartments of the knee: contact force, peak pressure, average pressure, and contact area. The knee was fixed at a flexion angle of 30o and a Tekscan Kscan sensor (Tekscan, South Boston, MA) inserted anteriorly into the medial and lateral j oint space. The mediallateral position of the femur was adjusted to produce an approximately 70% medial30% lateral load split between the two sides (Hurwitz et al., 1998; Schipplein and Andriacchi, 1991). The specimen was subjected to three trials of a 4 second ramp load from 200 to 1000 N. At the end of each ramp, the locations of the six screw heads were digitized using a Microscribe 3DX digitizer (Immersion Corp., San Jose, CA) possessing an accuracy of 0.23 mm. Figure 21. A human cadaver knee static experiment setup. A) The knee was potted in neutral alignment with six screws. B) The Tekscan Kscan sensor. C) The knee was mounted with fixed 300 flexion angle in a servohydraulic test machine with a sensor to measure intraarticular contact quantities. D) The closeup view of the contact area. Drift in the Tekscan sensor (Otto et al., 1999) was eliminated by postcalibrating each trial with the manufacturersuggested twopoint calibration procedure using the known loads at the start and end of the ramp. Crinkling of the sensor (Harris et al., 1999), which introduces erroneous pressures on sensels outside the true contact area, was accounted for by determining the pressure cutoff value (0.05 MPa) above which little additional drop in contact area occurs when pressures below this value are set to zero (Fregly et al., 2003). Contact quantities measured with the Tekscan sensor were therefore calculated by ignoring all sensels with a pressure below 0.05 MPa. Following Tekscan sensor calibration and pressure cutoff determination, the four experimental quantities of interest were calculated on each side for applied loads of 500 and 1000 N. Peak pressure was calculated using the averaging function in the Tekscan software, thereby reducing the effect of local sensor "hot spots" on this quantity. The resulting data from the Kscan sensor (Table 21) and the digitizer were averaged over the three trials to facilitate contact model evaluation under two loading conditions. Table 21. The averaged experimental measurements were collected from three compression trials for both loads processing on the same cadaver knee via servohydraulic test machine Experimental Load Experimental Quantity Side 500 N 1000 N Force (N) Medial 317 & 4 658 & 5 Max Pressure (MPa) 4.10 & 0.05 7.94 & 0.14 Ave Pressure (MPa) 1.11 & 0.03 2.07 & 0.02 Area (mm2) 287 & 10 318 & 4 Force (N) Lateral 183 & 4 337 & 5 Max Pressure (MPa) 1.51 & 0.03 2.63 & 0.04 Ave Pressure (MPa) 0.79 & 0.01 1.33 & 0.03 Area (mm2) 229 & 4 252 & 5 2.3 Knee Model Creation Prior to experimental contact testing, MRI (magnetic resonance imaging) and CT (computed tomography) data were collected from the same cadaveric specimen for purposes of contact model creation. Sagittal plane MRI data were collected using a 3.0T GE Signa Horizon LX scanner with a quadrature knee coil. A T2weighted 3D FastGRE sequence was used with a 1 mm slice thickness, 256 x 256 image matrix (0.625 x 0.625 mm pixel size), and 160 x 160 mm field of view. Axial CT data were collected from the same specimen using a GE LightSpeed QX/i scanner in helical mode. The scanning parameters were a 1.25 mm overlapping slice thickness, 512 x 512 image matrix (0.313 x 0.313 mm pixel size), and 160 mm x 160 mm field of view. The tibia, femur, and bone screws in both data sets were segmented (Fig. 22) using commercial image processing software (SliceOmatic, Tomovision, Montreal, CA). Figure 22. Original and segmented medical images. A) Original CT slice. B) Segmented CT slice. C) Original MRI slice. D) Segmented MRI slice. The menisci were not segmented and were omitted from the model. Articular cartilage and subchondral bone surfaces were segmented manually from the MRI data, while cortical bone and bone screw surfaces were segmented semiautomatically from the CT data using a watershed algorithm. The 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 MRI and CT point cloud data into a combined geometric model for contact analysis. Point clouds from each imaging modality were imported separately and converted to polygonal surface models. The subchondral bone surfaces from MRI were registered automatically to the corresponding cortical bone surfaces from CT, creating a composite geometric model with articular cartilage surfaces from MRI and cortical bone and bone screw surfaces from CT. NURBS (NonUniform Rational BSpline) surfaces were fitted to the polygonal models, with the tolerance (mean + standard deviation) between the original point clouds from MRI and the final NURBS surfaces being 0. 18 f 0. 18 mm for the femur and 0.20 f 0.29 mm for the tibia (Fig. 23). Figure 23. The anterior and posterior views of the A) 3D point clouds and B) 3D NURBS model with six screws created from CT and MRI images The NURBS surfaces for the tibia and femur articularr cartilage, cortical bone, and bone screws) were imported into Pro/MECHANICA MOTION (Parametric Technology Corporation, Waltham, MA) to construct a multibody contact model. The mean digitized bone screw locations were also imported to determine the nominal alignment of the tibia and femur. For both bones, a stiff linear spring was placed between each screw head and its mean experimental location and a static analysis performed to determine the pose that best matched the experiments. Differences between the digitized and nominal bone screw locations were on the order of 1 mm. Starting from these nominal poses, the tibia was fixed to ground and the femur connected to it via a 6 degreeoffreedom (DOF) joint. Custom contact code was incorporated into the multibody model and used to solve for the medial and lateral contact conditions as a function of the 6 DOFs between the two bones (Bei and Fregly, 2004). The contact code implemented two versions of a linear elastic discrete element contact model. The first was a small strain version, where the contact pressure p for each contact element on the tibial articular surfaces was calculated from (An et al., 1990; Blankevoort et al., 1991; Li et al., 1997) (1 v)E d p = (8) (1+ v)(1 2v) h where E is Young' s modulus of the articular cartilage, v is Poisson' s ratio, h is the combined thickness of the femoral and tibial articular cartilage, and d is the interpenetration of the undeformed contact surfaces. Both h and d were calculated on an elementbyelement basis using the ACIS 3D Toolkit (Spatial Corporation, Wesminster, CO). For large strains, a second version of the model was implemented that accounted for geometric nonlinear behavior (Blankevoort et al., 1991): (1 v)E I'dl p = In1 (9) (1 + v)(1 2v) hi For both versions, a dense contact element grid of 50 x 50 was used for the medial and material articular surfaces of the tibia. The femoral and tibial articular cartilage were assumed to be linear elastic and isotropic with Poisson's ratio = 0.45 (Blankevoort et al., 1991). Young's modulus was set to 1 MPa to facilitate its use as a design variable in the response surface optimizations. 2.4 Knee Model Evaluation Seventy seven contact analyses were performed with the model to provide the sample points necessary to generate response surfaces using the FCCCD. Each sample point represented a different pose of the femur relative to the tibia within the neighborhood of the nominal pose. This neighborhood was defined to be & 1 mm and a 10 based on the estimated envelope of experimental pose uncertainty. Though this envelope appears small, it corresponds to large changes in contact conditions. Within this envelope, response surfaces were generated as described above for the medial and lateral contact force, peak pressure, average pressure, and contact area computed by the contact model. The optimization cost function g(x) constructed from these response surfaces sought to match experimental average pressures pospin both compartments simultaneously with a penalty term to ensure the experimental contact forces F were reproduced: g(x, E) = ( pm, E pove (x)) tedral (are E pove (x)) ar,,w; + w ~F EF(x)) medral + (F EF(x)) ateral1 (10) where ~(x) and pore, (x) are the force and average pressure predicted by response surfaces. In this equation, x represents the 6 pose parameter design variables, E is the material modulus seventh design variable, and w = 1000 is the weight of the penalty term. This cost function mimics the results of a static analysis, since contact force is matched closely while minimizing errors in the most reliable experimental measure of interest. The form of Eq. (10) was specified by following a trialanderror approach that included different quantities in the cost function. A larger value of w was not used since it resulted in a poor match to the other contact quantities. To seek the global minimum, we performed 500 nonlinear leastsquares optimizations using the response surfaces from each contact model and the Matlab Optimization Toolbox (The Mathworks, Natick, MA). Uniformly distributed random initial guesses were selected within the bounds + 1 for the first six design variables and 1 to 10 for the seventh. The best set of design variables was selected based on the smallest cost function value from the 500 optimizations. A final Pro/MECHANICA contact analysis was performed for both contact models using the optimized design variables to verify the accuracy of the response surface approximations. CHAPTER 3 RESULTS All %RMiS~ad and %RMSFvress results were less than 5% while all R2, vadV8UeS were greater than or equal to 0.995. Error measures from the small and large strain versions of the contact model were generally close. %RM~SEvress was always larger than %RM~SEad but not dramatically so. Furthermore, %RM~SEady and %RM~SEvress were of comparable magnitude to the variability in the corresponding experimental measurements (Table 31). For each RS, a minimum of 7 and maximum of 26 infeasible/outlier points were eliminated from the original FCCCD set of 77 sample points, and all of these points involved a superior translation of +1 mm (i.e., no contact or light contact). Based on the best results from the 500 RS optimizations, the contact quantities computed from response surfaces were compared to both versions of the contact model (Fig. 31). The relative errors from all contact quantities were below 10%. Both versions of the contact model could reproduce all experimentally measured contact quantities to within 10% error at both loads (Fig. 32). The one exception was peak pressure, which was only matched to within 50% error. Contact force was matched to within 1% error, consistent with the use of a penalty term on this quantity in the cost function. Average pressure and contact area errors were always in opposite directions, consistent with matching contact force closely. Peak pressure was overpredicted by the contact models on the lateral side and underpredicted on the medial side (worst error). The optimal value of Young' s modulus was 2.6 MPa for the small strain model and 2.3 MPa for the large strain one. The best result from 500 optimizations never hit the upper or lower bounds on the pose parameter design variables. Each set of 500 optimizations required approximately 90 seconds of CPU time on a 2.8 GHz Pentium 4 PC. Table 31. Comparison of response surface predictions to data points sampled from large and small strain contact models. RMSE, Side Large Small Medial 4.67 4.24 0.065 0.035 0.059 0.048 10.8 10.8 Lateral 2.12 2.12 0.030 0.018 0.026 0.024 4.54 4.54 %RMSEag RMSE,, %RMSE . Predicted Quantity Force (N) Max Pressure (MPa) Ave Pressure (MPa) Area (mm2) Force (N) Max Pressure (MPa) Ave Pressure (MPa) Area (mm2) Large 0.601 1.53 3.31 3.50 0.788 1.15 2.24 2.58 Small 0.674 1.16 3.19 3.50 0.886 0.815 2.27 2.58 Large 10.0 0.090 0.078 14.3 2.92 0.039 0.033 6.00 Small 6.13 0.048 0.064 14.3 2.85 0.023 0.030 6.00 Large 1.29 2.13 4.40 4.64 1.09 1.49 2.91 3.41 Small 0.976 1.58 4.25 4.64 1.19 1.06 2.97 3.41 Large 0.999 0.999 0.999 0.995 0.999 0.999 0.998 0.996 Small 0.999 0.999 0.999 0.995 0.999 0.999 0.998 0.996 Pmax B S500N Large strain S500N Small strain ~i1000N Large strain O 1000N Small strain Medial Lateral Pave Force 20 E 20 20 20 40 Medial Lateral Medial Lateral Figure 31. Percent errors between response surfaces and predicted results from both contact models. They were calculated for A) contact force, B) contact peak pressure, C) contact averaged pressure and D) contact area. Pose parameters and Young's modulus were defined based on the best result from 500 optimization runs. Force 20 A S500N Large strain 2 20 500N Small strain 1000N Large strain a 1000N Small strain 40 Medial Lateral Pave 20 40 Pmax c.B o o 0 Medial Lateral Area 01 D 1 Medial Lateral Medial Lateral Figure 32. Percent errors between experimental and predicted results from both contact models. They were calculated for A) contact force, B) contact peak pressure, C) contact averaged pressure and D) contact area. Pose parameters and Young's modulus were defined based on the best result from 500 optimization runs. CHAPTER 4 DISCUSSION This thesis has presented a novel response surface optimization methodology for evaluating a joint contact model's ability to reproduce static experimental contact measurements. The approach modifies traditional RS approximation methods for application to contact analyses. By replacing computationallycostly contact analyses with quadratic RS approximations, optimizations that vary the relative pose of the contacting bodies can be performed rapidly to minimize differences between model predictions and experimental measurements. Evaluation of a discrete element contact model of the tibiofemoral j oint constructed from CT and MRI data revealed that quadratic RSs can accurately approximate model outputs within a small envelope of relative pose uncertainty. Furthermore, optimization studies utilizing these RSs demonstrated that the model could reproduce the contact force, average pressure, and contact area. However, the peak pressure measured experimentally on the medial and lateral sides could not be reproduced. This finding suggests that discrete element contact models of natural human j points that utilize homogeneous, isotropic material properties may be best suited for analyses such as multibody dynamic simulations where obtaining correct contact forces by integrating over the surfaces is more critical than obtaining correct contact pressures at specific locations on the surfaces. Use of RSs to replace repeated contact analyses in optimization studies is worthwhile for several reasons. First, RS tends to smooth out noise in the design space. During optimization, the risk of entrapment in a local minimum is therefore reduced. Second, RS approximations are computationally efficient. Rather then repeating costly contact analyses during an optimization, a single set of contact analyses is performed once up front to generate the necessary RS approximations for output quantities of interest. Extremely fast function evaluations allow one to search for the global optimum using either repeated gradientbased optimization starting from multiple initial guesses, as performed in our study, or global optimization employing a large population size. Third, RS optimizations are convenient to implement. Optimized contact solutions can be founded utilizing any offtheshelf optimization algorithm once the RS approximations are constructed. Fourth, a variety of optimization problem formulations can be evaluated quickly. Each contact quantity that could potentially appear in the cost function or constraints can be fitted with its own RS. A wide variety of cost functions can then be constructed by weighting contributions from the different RSs. Fifth, RS approximations facilitate the calculation of analytical derivatives for gradientbased optimization. This benefit is the direct result of having analytical representations for the contact quantities of interest as a function of the design variables. The caveat to these benefits is that an appropriate mathematical form (polynomial or otherwise) must be identified to represent the responses of interest as a function of the design variable inputs, which is not always possible. The primary deficiency of the discrete element contact model was its inability to match the peak contact pressures measured experimentally. Peak pressures were not included in the cost function since they are sensitive to local inhomogeneities in cartilage material properties and Tekscan sensor response. When we replaced average pressures with peak pressures in the cost function for curiosity, we found that peak pressure errors decreased on the medial side by about 25% but increased on the lateral side by roughly 50%, while average pressure and contact area errors increased by approximately 30% and 20%, respectively. Thus, a variety of other factors likely contributed to the large peak pressure errors, including insertion of a relatively stiff sensor into the joint space (Wu et al., 1998), small inaccuracies in the surface geometry, and local variations in material properties. Underprediction of peak pressure by the model on the moderately conformal medial side and overprediction on the nonconformal lateral side is consistent with observations made by Wu et al. (1998) on how insertion of a relatively stiff sensor into the joint space affects peak pressure measurements as a function of contact conformity. Contact model limitations may have also contributed to the poor match of the peak pressure data. Though a strength of our model is that it accounts for local variations in cartilage thickness, a weakness is that it does not account for local variations in cartilage material properties. Though homogeneous, isotropic cartilage material models have been used in the literature (Bendj aballah et al., 1995; Blankevoort et al., 1991; Haut et al., 2002; Perie and Hobatho, 1998), recent studies have reported significant local variations in Young's modulus, Poisson's ratio and thickness for articular cartilage (Jurvelin et al., 2000; Laasanen et al., 2003). Mukherj ee and Wayne (1998) suggested that regions with the highest Young's modulus correspond to regions with the highest contact pressure and cartilage thickness. This observation fits our underprediction of peak pressure on the medial side, where an increase in Young's modulus would produce improved agreement with the large measured peak pressures. In addition, our elastic contact model does not account for the effect of timedependent fluid flow on contact pressures or areas as captured by biphasic models of cartilage (Han et al., 2004; Mow et al., 1980; Mow et al., 1982). However, for loading over a short time period, an elastic model may still provide a reasonable approximation of the in vivo situation depending on the intended application of the model (Donzelli et al., 1999; Mow et al., 1982; Shepherd and Seedhom, 1997;). The value of Young' s modulus predicted by our optimizations is consistent with studies found in the literature. Experimental studies of the compressive modulus of articular cartilage have reported values as low as 2.0 MPa for short time frame loading (Setton et al., 1999; Shepherd and Seedhom, 1997). Other studies have reported the modulus of the solid phase alone to be in the neighborhood of 0.3 MPa (Hasler et al., 1999). Thus, the bestfit modulus values of2.3 and 2.6 MPa found in our study are consistent with these numbers. Other discrete element knee studies have utilized an elastic modulus of 4 MPa (Cohen et al., 2003; Kwak et al., 2000) which again is close to our optimized modulus values. Furthermore, if our estimate of Poisson' s ratio is decreased slightly from 0.45 to 0.4, our optimized value of Young' s modulus will increase from 2.6 to 4.6 MPa, closer to the middle of the range commonly reported (Setton et al., 1999; Shephard and Seedhom, 1997). A weakness of our study is that only a single knee specimen was evaluated. Significant modeling effort is required to create a knee model with articular cartilage and subchondral bone geometry that matches a particular cadaver specimen. This may explain why other studies that utilized specimenspecific geometric models only analyzed a single specimen as well (Bendjaballah et al., 1995; Haut et al., 2002; Li et al., 1999). The knee used in our study was part of a larger experimental study involving 20 knees. However, we were able to collect CT and MRI data from only two of those knees. Since the contact area for the second specimen was not completely contained within the boundary of the Tekscan sensor at the maximum load of 1000 N, we were not able to postcalibrate those experimental contact data and use them for a second model evaluation. Comparison of model predictions with experimental measurements from additional knee specimens using a range of flexion angles and loads would provide further verification of the capabilities and limitations of discrete element models of natural joints. CHAPTER 5 SUMMARY AND FUTURE STUDY 5.1 Summary In summary, this thesis has presented a computationally efficient RS optimization for predicting knee contact measurements using a 3D simulation model. The more computational costly is each contact analysis; the more beneficial is the use of RS optimization method. The present implementation works for the tibiofemoral j oint of the natural knee with either large or small strain contact model. The RS optimization can predict contact forces, areas and average pressures well with the exception of peak pressures. Our ultimate goal is to incorporate the elastic foundation contact model into a fullbody dynamic musculoskeletal model. The resulting fullbody model can then be utilized to study joint contact pressures during motion. Future research is required to evaluate whether the peak pressure errors are due to experimental or contact model inaccuracies. 5.2 Future Study While our results showed great potential to predict contact quantities that involve averaging across the surface, there are some extra steps that could be taken to improve the accuracy of RS. Instead of using two modifications in this thesis to remove the outliers ,there is another more robust procedure called iteratively reweighted least square (IRLS) fitting (Holland and Welsch, 1977). It can be utilized to remove or weight down outliers. Additionally, varying numbers of insignificant coefficients for each RS was found based on tstatistics analysis. While it was not performed in this study, there are several statistical methods, (e.g., a stepwise regression procedure) that can be used to discard those insignificant coefficients to improve the prediction accuracy. In addition to 500 optimizations, another 50 groups of optimizations were run. The best result was selected from each group containing 100 optimizations with different initial guesses. Within 50 sets of optimized design variables, the maximum variations in translation and rotation were 1.90 mm and 1.950 respectively. Compared to the range of design space, these variations suggested that the final result from 500 optimizations may not be unique. Thus an alternative cost function or type of different RS method (e.g., neural networks and kriging) could be another candidate for future improvement. Given how well RSs could approximate contact model outputs for a small range of relative pose variations, our next step will be to investigate whether RSs can be used to generate extremely fast forward dynamic contact simulations that utilize a wider range of relative pose variations. The computational cost of repeated geometry evaluations is the current limiting factor to the incorporation of deformable contact models into multibody dynamic simulations of human j points (Bei and Fregly, 2004). However, if accurate response surfaces could be generated to represent the net force and torque calculated about a preselected point on the tibia (e.g., the origin of the tibial coordinate system), then extremely fast contact solutions could be produced for any given set of relative pose parameters. Though additional design variables would need to be included to account for friction or damping effects, these effects are expected to be small for most modeling situations involving human joints. Higher degree polynomials or some other mathematical formulation would likely be required to represent a wider range of relative pose variations. LIST OF REFERENCES Ahn, C.S. & Kim, K.Y. (2003). Aerodynamic design optimization of a compressor rotor with NavierStokes analysis. Proceedings of the hIstitution of M~echanical Engineers, PartPPPPP~~~~~~~PPPPPP A: Journal of Power and Energy, 217, 1791 8. An, K.N., Himeno, H., Tsumura, H., Kawai, T., & Chao, E.Y.S. (1990). Pressure distribution on articular surfaces: application to joint stability evaluation. Journal of Bionzechanics, 23, 10131020. Bei, Y. & Fregly, B. J. (2004). Multibody dynamic simulation of knee contact mechanics. Medical Engineering & Physics (submitted). Bendjaballah, M.Z., ShiraziAdl, A., & Zukor, D.J. (1995). Biomechanics of the human knee joint in compression: Reconstruction, mesh generation and Finite element analysis. The Knee, 2, 6979. Bendjaballah, M.Z., ShiraziAdl, A., & Zukor, D.J. (1997). Finite element analysis of human knee j oint in varusvalgus. Clinical Bionzechanics, 12, 139148. Bendjaballah, M.Z., ShiraziAdl, A., & Zukor, D.J. (1998). Biomechanical response of the passive human knee joint under anteriorposterior forces. Clinical Bionzechanics, 13, 625633. Blankevoort, L., Kuiper, J.H., Huiskes, R., & Grootenboer, H.J. (1991). Articular contact in a threedimensional model of the knee. Journal ofBionzechanics, 24, 10191031i. Burman, J. & Gebart, B.R. (2001). Influence from numerical noise in the objective function for flow design optimization. hIternational Journal of Nunterical M~ethods for Heat and Fluid Flow, 11, 619. Chang, P.B., Williams, B.J., Santer, T.J., Notz, W.I., & Bartel, D.L. (1999). Robust optimization of total joint replacements incorporating environmental variables. Journal ofBionzechanical Engineering, 121, 3 043 10. Cohen, Z.A., Henry, J.H., McCarthy, D.M., Mow, V.C., & Ateshian, G.A. (2003). Computer simulations of patellofemoral joint surgery. The American Journal of Sports M~edicine, 31, 8798. Dhaher, Y.Y. & Kahn, L.E. (2002). The effect of vastus medialis forces on patello femoralcontact: a modelbased study. Journal of Bionzechanical Engineering, 124, 758767. Donahue, T.L.H., Rashid, M.M., Jacobs, C.R., & Hull, M.L. (2002). A finite element model of the human knee joint for the study of tibiofemoral contact. Journal of Bionzechanical Engineering, 124, 273280. Donzelli, P. S., & Spilker, R. L. (1998). A contact finite element formulation for biological soft hydrated tissues. Computer M~ethods in Applied M~echanics and Engineerin, 153, 6379. Donzelli, P.S., Spilker, R.L., Ateshian, G.A., & Mow, V.C. (1999). Contact analysis of biphasic transversely isotropic cartilage layers and correlations with tissue failure. Journal ofBionzechanics, 32, 10371047. Fregly, B.J., Bei, Y., & Sylvester, M.E. (2003). Experimental evaluation of a multibody dynamic model to predict contact pressures in knee replacements. Journal of Bionzechanics, 36, 16581668. Haider, M.A., & Guilak, F. (2000). Axisymmetric boundary integral model for incompressible linear viscoelasticity: Application to the micropipette aspiration contact problem. Journal ofBionzechanical Engineering, 122, 236244. Han, S.K., Federico, S., Epstein, M., & Herzon, W. (2004). An articular cartilage contact model based on real surface geometry. Journal ofBionzechanics (in press). Harris, M.L., Morberg, P., Bruce, W.J.M., & Walsh, W.R. (1999). An improved method for measuring tibiofemoral contact areas in total knee arthroplasty: A comparison of Kscan sensor and Fuji film. Journal ofBionzechanics, 32, 951958. Hasler, E.M., Herzog, W., Wu, J.Z., Muller, W., & Wyss, U. (1998). Articular cartilage biomechanics: theoretical models, material properties, and biosynthetic response. Critical Reviews in Biomedical Engineering, 27, 415488 Haut, T.L., Hull, M.L., Rashid, M.M., & Jacobs, C.R. (2002). A finite element model of the human knee joint for the study of tibiofemoral contact. Journal of Bionzechanical Engineering, 124, 273280. Herzog, W., Longino, D., & Clark, A. (2003). The role of muscles in joint adaptation anddegeneration. Langenbeck 's Archives ofSurgely, 388, 305315. Holland, P.W., & Welsch, R.E. (1977). Robust regression using iteratively reweighted leastsquares. Conanunications in Statistics: Theory and Methods, 6, 813827 Hong, J.H., Mun, M.S., & Song, S.H. (2001). An optimum design methodology development using a statistical technique for vehicle occupant safety. Proceedings of the Institution of M~echanical Engineers, Part D: Journal of Automobile Engineering, 215, 795801. Hurwitz, D.E., Sumer, D.R., Andriacchi, T.P., & Sugar, D.A. (1998). Dynamic knee loads during gait predict proximal tibial bone distribution. Journal ofBiomechanics, 31, 423430. Jansson, T., Nilsson, L., & Redhe, M. (2003). Using surrogate models and response surfaces in structural optimization With application to crashworthiness design and sheet metal forming. Structural and M~ultidisciplinary Optimization, 25, 129140. Johnson, K.L. (1985). Contact Mechanics. Cambridge University Press, New York. Jung, E.S. & Choe, J. (1996). Human reach posture prediction based on psychophysical discomfort. International Journal oflndustrial Ergonomics, 18, 173179. Jurvelin, J.S., Arokoski, J.P., Hunziker, E.B., & Helminen, H.J. (2000). Topographical variation of the elastic properties of articular cartilage in canine knee. Journal of Biomechanics, 33, 66975. Kaufman, K.R., Kovacevic, N., Irby, S.E., & Colwell, C.W. (1996). Instrumented implant for measuring tibiofemoral forces. Journal ofBiomechanics, 29, 667671. Keane, A.J. (2003). Wing optimization using design of experiment, response surface, and data fusion methods. Journal of Aircraft, 40, 741750. Komistek, R. D., Dennis, D. A., & Mahfouz, M. (2003). In vivo fluoroscopic analysis of the normal human knee. Clinical Oi Iluspe~l~ iL a ~nd Rela.ted Resea~rch, 410, 698 1. Kuo, C.H. & Keer, L.M. (1993). Contact stress and fracture analysis of articular cartilage. Biomedical Engineering Applications Basis Communications, 5, 5 15521. Kwak, S.D., Blankevoort, L., & Ateshian, G.A. (1999). A mathematical formulation for 3D quasistatic multibody models of diarthrodial joints. Computer M~ethods in Biomedical Engineering, 3, 4164. Laasanen, M.S., Toyras, J., Korhonen, R.K., Rieppo, J., Saarakkala, S., Nieminen, M.T., Hirvonen, J., & Jurvelin, J.S. (2003). Biomechanical properties of knee articular cartilage, Biorheology, 40, 133140. Leary, S.J., Bhaskar, A., & Keane A.J. (2004). Global approximation and optimization using adj oint computational fluid dynamics codes. A1AA Journal, 42, 63 1641. Li, G., Sakamoto, M., & Chao, E.Y.S. (1997). A comparison of different methods in predicting static pressure distribution in articulating joints. Journal ofBiomechanics, 30, 635638. Li, G., Gil, J., Kanamori, A., & Woo, S.L.Y. (1999). A validated 3D computational model of a human knee joint. Journal ofBiomechanical Engineering, 121, 657662. Li, G., Lopez, O., & Rubash, H. (2001). Variability of a 3D finite element model constructed using MR Images of a knee for joint contact stress analysis. Journal of Bionzechanical Engineering, 123, 341346. Liu, B., Haftka, R.T., & Akgun, M.A. (2000). Twolevel composite wing structural optimization using response surfaces. Structural and' M\ultidisciplinaly Optimization, 20, 8796. MacWilliams, B.A., DesJardins, J.D., Wilson, D.R., Romero, J., & Chao, E.Y.S. (1998). A repeatable alignment method and local coordinate description for knee joint testing and kinematic measurement. Journal ofBionzechanics, 31, 947950. Mow, V.C., Kuei, S.C., Lai, W.M., & Armstrong, C.G. (1980). Biphasic creep and stress relaxation of articular cartilage in compression: Theory and experiments. Journal of Bionzechanical Engineering, 102, 7384. Mow, V.C., Lai, W.M., & Holmes, M.H. (1982). Advanced' theoretical and' experimental techniques in cartilage research. In Bionzechanics: Principles and' Applications (edited by R. Huiskes, D. van Campen, and J. DeVijn), Martinus Nijhoff Publishers, Boston, 6369 Mukherjee, N. & Wayne, J.S. (1998). Load sharing between solid and fluid phases in articular cartilage: IExperimental determination of in situ mechanical conditions in a porcine knee. Journal ofBionzechanical Engineering, 120, 614619. Otto, J.K., Brown, T.D., & Callaghan, J.J. (1999). Static and dynamic response of a multiplexedarray piezoresistive contact sensor. Experimental M~echanics, 39, 317 323. Pandy, M.G. & Sasaki, K. (1998). A threedimensional musculoskeletal model of the human knee joint. Part 2: Analysis of ligament function. Computer M~ethods in Bionzechanics and Biomed'ical Engineering, 1, 265283. Papila, N., Shyy, W., Griffin, L., & Dorney, D.J. (2002). Shape optimization of supersonic turbines using global approximation methods. Journal of Propulsion and' Power, 18, 509518. Perie, D. & Hobatho, M.C. (1998). In vivo determination of contact areas and pressure of the femorotibial joint using nonlinear finite element analysis. Clinical Bionzechanics, 13, 394402. Piazza, S.J. & Delp, S.L. (2001). Threedimensional dynamic simulation of total knee replacement motion during a stepup task. Journal of Bionzechanical Engineering, 123, 599606. Rikards, R., Abramovich, H., Auzins, J., Korjakins, A., Ozolinsh, O., Kalnins, K., & Green, T. (2004). Surrogate models for optimum design of stiffened composite shells. Composite Structures, 63, 243251. Roux, W.J., Stander, N., & Haftka, R.T. (1998). Response surface approximations for structural optimization. International Journal for Numerical M~ethods in Engineering, 42, 517534. Sevant, N.E., Bloor, M.I.G., & Wilson, M.J. (2000). Aerodynamic design of a flying wing using response surface methodology. Journal ofAircraft, 37, 562569. Schipplein, O.D. & Andriacchi, T.P. (1991). Interaction between active and passive knee stabilizers during level walking. Journal of Orthopaedic Research, 9, 113119. Setton, L.A., Elliot, D.M., & Mow, V.C. (1999). Altered mechanics of cartilage with osteoarthritis: human osteoarthritis and an experimental model of joint degeneration. Osicoat~l th; iti\ and' Cartilage, 7, 214. Shepherd, D.E.T. & Seedhom B.B. (1997). A technique for measuring the compressive modulus of articular cartilage under physiological loading rates with preliminary results. Journal of Engineering in M~edicine, 211, 155165. Stolk, J., Verdonschot, N., Cristofolini, L., Toni, A., & Huiskes, R. (2002). Finite element and experimental models of cemented hip joint reconstructions can produce similar bone and cement strains in preclinical tests. Journal ofBiomechanics, 35, 499510. Tashman, S. & 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, 23 8245. Wu, J.Z., Herzog, W., & Epstein, M. (1998). Effects of inserting a pressensor film into articular joints on the actual contact mechanics. Journal of Biomechanical Engineering, 120, 655659. BIOGRAPHICAL SKETCH YiChung Lin was born on August 14, 1976, in Kaohsiung, Taiwan. In 1998, he received the bachelor' s degree in the Department of Mechanical Engineering of the National Cheng Kung University, Taiwan. He j oined the Computational and Biomechanics Laboratory of Assistant Professor B. J. Fregly in the fall of 2002. He is planning to enter the doctoral program in mechanical and aerospace engineering after finishing the master's program. 