|UFDC Home||myUFDC Home | Help|
This item has the following downloads:
EXPERIMENTAL EVALUATION OF A NATURAL KNEE CONTACT MODEL
USINTG RESPONSE SURFACE OPTIMIZATION
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
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
ACKNOWLEDGMENT S .........__.. ..... .__. .............._ iii..
LIST OF TABLES ........._.___..... .__. ...............v....
LIST OF FIGURES .............. ....................vi
AB STRAC T ................ .............. vii
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....
LIST OF TABLES
2-1. 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....
3-1. Comparison of response surface predictions to data points sampled from large and
small strain contact model s. ............. ...............20.....
LIST OF FIGURES
2-1. A human cadaver knee static experiment setup .........._._ ....... ..._ .............._.11
2-2. Original and segmented medical images ..........._..__......_. ......._._..........1
2-3. 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
3-1. Percent errors between response surfaces and predicted results from both contact
m odels .............. ...............21....
3-2. 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
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 computationally-efficient 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.
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
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 model-based 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 multi-dimensional 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.
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.
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 multi-dimensional 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 low-degree
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 low-degree 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 least-squares 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 D-optimality
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
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 pre-selected cut-off value of 10% (a typical value
for engineering analyses) is omitted and the RS re-generated 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 root-mean-square error ( RMSEa4 ).
Given the SSE, the root-mean 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 non-sampled 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):
RM~SE,, =,1 (4)
To provide a relative measure of fit quality, we also compute the percent adjusted RM\~SE
%RM~SEad =-1 (5)
where y~ represents the magnitude of the fitted quantity:
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
PRESS-based SSE called the PRESS statistic. From there, a PRESS-based 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 medial-lateral
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
self-align with the femur once an axial load was applied.
Using this set-up (Fig. 2-1), 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 K-scan sensor (Tekscan, South Boston, MA) inserted anteriorly into the medial
and lateral j oint space. The medial-lateral position of the femur was adjusted to produce
an approximately 70% medial-30% 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 2-1. A human cadaver knee static experiment setup. A) The knee was potted in
neutral alignment with six screws. B) The Tekscan K-scan sensor. C) The
knee was mounted with fixed 300 flexion angle in a servohydraulic test
machine with a sensor to measure intra-articular contact quantities. D) The
close-up view of the contact area.
Drift in the Tekscan sensor (Otto et al., 1999) was eliminated by post-calibrating
each trial with the manufacturer-suggested two-point 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 cut-off 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 cut-off 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 K-scan
sensor (Table 2-1) and the digitizer were averaged over the three trials to facilitate
contact model evaluation under two loading conditions.
Table 2-1. The averaged experimental measurements were collected from three
compression trials for both loads processing on the same cadaver knee via
servohydraulic test machine
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.0-T
GE Signa Horizon LX scanner with a quadrature knee coil. A T2-weighted 3D FastGRE
sequence was used with a 1 mm slice thickness, 256 x 256 image matrix (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. 2-2) using commercial image processing
software (SliceOmatic, Tomovision, Montreal, CA).
Figure 2-2. 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 semi-automatically 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
(Non-Uniform Rational B-Spline) 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. 2-3).
Figure 2-3. 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 degree-of-freedom (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
element-by-element 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
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
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 trial-and-error 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 least-squares
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.
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 3-1). 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. 3-1). 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. 3-2). 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 over-predicted by the contact models
on the lateral side and under-predicted 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 3-1. Comparison of response surface predictions to data points sampled from large and small strain contact models.
Side Large Small
Medial 4.67 4.24
Lateral 2.12 2.12
Max Pressure (MPa)
Ave Pressure (MPa)
Max Pressure (MPa)
Ave Pressure (MPa)
S500N Large strain
S500N Small strain
~i1000N Large strain
O 1000N Small strain
Figure 3-1. 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
S500N Large strain
2 20 500N Small strain
1000N Large strain
a- 1000N Small strain
Figure 3-2. 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
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 computationally-costly 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 gradient-based 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 off-the-shelf 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 gradient-based 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
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 non-conformal 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 under-prediction 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 time-dependent 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 best-fit 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 specimen-specific 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
post-calibrate 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
SUMMARY AND FUTURE STUDY
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
full-body dynamic musculoskeletal model. The resulting full-body 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
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 re-weighted 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 t-statistics 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 pre-selected 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
LIST OF REFERENCES
Ahn, C.-S. & Kim, K.-Y. (2003). Aerodynamic design optimization of a compressor rotor
with Navier-Stokes analysis. Proceedings of the hIstitution of M~echanical
Engineers, PartPPPPP~~~~~~~PPPPPP A: Journal of Power and Energy, 217, 179-1 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, 1013-1020.
Bei, Y. & Fregly, B. J. (2004). Multibody dynamic simulation of knee contact mechanics.
Medical Engineering & Physics (submitted).
Bendjaballah, M.Z., Shirazi-Adl, A., & Zukor, D.J. (1995). Biomechanics of the human
knee joint in compression: Reconstruction, mesh generation and Finite element
analysis. The Knee, 2, 69-79.
Bendjaballah, M.Z., Shirazi-Adl, A., & Zukor, D.J. (1997). Finite element analysis of
human knee j oint in varus-valgus. Clinical Bionzechanics, 12, 139-148.
Bendjaballah, M.Z., Shirazi-Adl, A., & Zukor, D.J. (1998). Biomechanical response of
the passive human knee joint under anterior-posterior forces. Clinical Bionzechanics,
Blankevoort, L., Kuiper, J.H., Huiskes, R., & Grootenboer, H.J. (1991). Articular contact
in a three-dimensional model of the knee. Journal ofBionzechanics, 24, 1019-1031i.
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, 6-19.
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 04-3 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, 87-98.
Dhaher, Y.Y. & Kahn, L.E. (2002). The effect of vastus medialis forces on patello-
femoralcontact: a model-based study. Journal of Bionzechanical Engineering, 124,
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 tibio-femoral contact. Journal of
Bionzechanical Engineering, 124, 273-280.
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, 63-79.
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, 1037-1047.
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, 1658-1668.
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, 236-244.
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
K-scan sensor and Fuji film. Journal ofBionzechanics, 32, 951-958.
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, 415-488
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 tibio-femoral contact. Journal of
Bionzechanical Engineering, 124, 273-280.
Herzog, W., Longino, D., & Clark, A. (2003). The role of muscles in joint adaptation
anddegeneration. Langenbeck 's Archives ofSurgely, 388, 305-315.
Holland, P.W., & Welsch, R.E. (1977). Robust regression using iteratively reweighted
least-squares. Conanunications in Statistics: Theory and Methods, 6, 813-827
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, 795-801.
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,
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, 129-140.
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, 173-179.
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, 669-75.
Kaufman, K.R., Kovacevic, N., Irby, S.E., & Colwell, C.W. (1996). Instrumented implant
for measuring tibiofemoral forces. Journal ofBiomechanics, 29, 667-671.
Keane, A.J. (2003). Wing optimization using design of experiment, response surface, and
data fusion methods. Journal of Aircraft, 40, 741-750.
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, 69-8 1.
Kuo, C.H. & Keer, L.M. (1993). Contact stress and fracture analysis of articular cartilage.
Biomedical Engineering Applications Basis Communications, 5, 5 15-521.
Kwak, S.D., Blankevoort, L., & Ateshian, G.A. (1999). A mathematical formulation for
3D quasi-static multibody models of diarthrodial joints. Computer M~ethods in
Biomedical Engineering, 3, 41-64.
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, 133-140.
Leary, S.J., Bhaskar, A., & Keane A.J. (2004). Global approximation and optimization
using adj oint computational fluid dynamics codes. A1AA Journal, 42, 63 1-641.
Li, G., Sakamoto, M., & Chao, E.Y.S. (1997). A comparison of different methods in
predicting static pressure distribution in articulating joints. Journal ofBiomechanics,
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, 657-662.
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, 341-346.
Liu, B., Haftka, R.T., & Akgun, M.A. (2000). Two-level composite wing structural
optimization using response surfaces. Structural and' M\ultidisciplinaly
Optimization, 20, 87-96.
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, 947-950.
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, 73-84.
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,
Mukherjee, N. & Wayne, J.S. (1998). Load sharing between solid and fluid phases in
articular cartilage: I--Experimental determination of in situ mechanical conditions in
a porcine knee. Journal ofBionzechanical Engineering, 120, 614-619.
Otto, J.K., Brown, T.D., & Callaghan, J.J. (1999). Static and dynamic response of a
multiplexed-array piezoresistive contact sensor. Experimental M~echanics, 39, 317-
Pandy, M.G. & Sasaki, K. (1998). A three-dimensional musculoskeletal model of the
human knee joint. Part 2: Analysis of ligament function. Computer M~ethods in
Bionzechanics and Biomed'ical Engineering, 1, 265-283.
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, 509-518.
Perie, D. & Hobatho, M.C. (1998). In vivo determination of contact areas and pressure of
the femorotibial joint using non-linear finite element analysis. Clinical
Bionzechanics, 13, 394-402.
Piazza, S.J. & Delp, S.L. (2001). Three-dimensional dynamic simulation of total knee
replacement motion during a step-up task. Journal of Bionzechanical Engineering,
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, 243-251.
Roux, W.J., Stander, N., & Haftka, R.T. (1998). Response surface approximations for
structural optimization. International Journal for Numerical M~ethods in
Engineering, 42, 517-534.
Sevant, N.E., Bloor, M.I.G., & Wilson, M.J. (2000). Aerodynamic design of a flying
wing using response surface methodology. Journal ofAircraft, 37, 562-569.
Schipplein, O.D. & Andriacchi, T.P. (1991). Interaction between active and passive knee
stabilizers during level walking. Journal of Orthopaedic Research, 9, 113-119.
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, 2-14.
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, 155-165.
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 pre-clinical tests. Journal ofBiomechanics, 35, 499-510.
Tashman, S. & Anderst, W. (2003). In-vivo measurement of dynamic joint motion using
high speed biplane radiography and CTL: application to canine ACL deficiency.
Journal ofBiomechanical Engineering, 125, 23 8-245.
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, 655-659.
Yi-Chung 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.