BIOMECHANICAL ANALYSIS OF A
RESECTION ARTHRODESIS
By
NORMAN R. WILLIAMSON
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1985
TO MY PARENTS,
JOHN AND INGER
ACKNOWLEDGMENTS
It has been a privilege to work under the insightful eye and
encouraging wing of my committee chairman, Dr. George Piotrowski,
whose assistance in matters technical and otherwise I gratefully
acknowledge. Thanks are extended to Dr. Gary Miller for suggesting this
problem, and for providing guidance and assistance throughout this work.
I am indebted to Dr. Chester Tylkowski for his patience and perseverence
in introducing me to the worlds of orthopedic surgery and gait analysis,
and for the encouragement and enthusiasm shown throughout this effort.
Thanks also go to Dr. Dempsey Springfield, who performed the sham
resection arthrodesis procedure and provided guidance on matters
pertaining to surgical procedure and the healing process. Dr. George
Sandor, for many years an unfailing source of wisdom and inspiration,
introduced me to the methods of optimization; for both I am very
grateful. Above all, I would like to thank my entire committee for
their friendship throughout the years of our association.
Dr. R. William Petty, Chairman, provided full access to the
facilities and resources of the University of Florida Department of
Orthopedics, for which I am grateful and without which this work
would not have been possible. The computed tomography imaging was
made possible through the generous assistance of Dr. Terry Hudson,
formerly of the University of Florida Department of Radiology and
currently at Massachusetts General Hospital in Boston.
The generous support for this effort by Dr. John Swanson and his
staff at Swanson Analysis Systems, Inc., is gratefully acknowledged.
I was fortunate to have the assistance of Mr. Art Smith of the
statistics section of IFAS, and of Mr. Allen Gregalot of the Department
of Engineering Sciences, at some very critical times during the course
of this work. Many thanks also go to Mr. Lance Jackson, Manager of the
Harvard University Science Center Computing Facility and a true Renais
sance figure, for his unending support and assistance.
Finally, I would like to thank Mrs. Edna Larrick for her assistance
in preparing the manuscript, and Mrs. Katarzyna Piercy and Mr. Dave Peace
for their help with many of the illustrations.
This research was aided by a grant from the Orthopedic Research
and Education Foundation.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS . . . . ... ....... iii
ABSTRACT . . . . ... .... . ... vii
CHAPTERS
1 INTRODUCTION . . . . . ...1
The Tibial Turnback . . . . ... ... 1
Problem Statement . . . . ... .. 5
Methodology . . . . ... . ... 7
Overview of the Present Work . . . .. 9
2 MODELLING OF THE RESECTION ARTHRODESIS . . .. .11
Introduction. . . . .... . 11
Structural Loading . . . .. 12
System Geometry Data Acquisition . . .... .21
Finite Element Modelling . . . .... .25
Some Final Considerations . . . .... .41
3 OPTIMIZATION THEORY . . . . ... 43
Introduction . ..... . . 43
The Optimization Problem . . . .... 44
Objective Function Formulation . . ... 47
Optimizing a Constrained Nonlinear Function . . 50
Constrained Optimization. . . . .... .51
Unconstrained Search Techniques . . .... .56
Global Optimization . . ..... ...63
Conclusions: Constrained Nonlinear Optimization Methods .. 63
4 RESULTS . . . . ... ....... .64
Part One: The Effects of Graft Length and Leg Position 65
Summary: Part One . . . . ... ... 84
Part Two: The Effect of IM Rod Structural Stiffness
and Leg Position in Long Resections . . ... 85
Summary Part Two . . . . 91
Part Three: Optimization of Rod Size and Leg Position
in Long Resections . . . ..... 92
Application of the Methods . . . ... 95
Summary: Part Three . . . .... 97
Summary of Results . . . . ... .. 100
TABLE OF CONTENTS (Continued)
CHAPTERS
5 CONCLUSIONS AND FURTHER WORK . . . .
Further Work . . . . . .
Conclusion . . . . . .
APPENDIX AFLOW CHART FOR THE HOOKE AND JEEVES
OPTIMIZATION ALGORITHM . . .
APPENDIX BFORTRAN SOURCE LIST: HOOKE AND JEEVES NONLINEAR
OPTIMIZATION ROUTINES FOR GLOBAL OPTIMIZATION
OF A CONSTRAINED FUNCTION . . .
REFERENCES . . . . . . .
BIOGRAPHICAL SKETCH . . . . . .
* 120
. 125
Page
: : :

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
BIOMECHANICAL ANALYSIS OF A
RESECTION ARTHRODESIS
By
Norman R. Williamson
May 1985
Chairman: George Piotrowski
Major Department: Mechanical Engineering
The tibiall turnback" resection arthrodesis of the knee is a method
of surgical intervention for patients with tumors located in the distal
femur. This procedure uses autogenous cortical bone grafts from the
anterior tibial shaft and the ipsilateral fibula to bridge the defect
created by the resection. A fluted intramedullary (IM) rod is used to
provide stable fixation.
Fatigue fractures of the posterior fibularr) graft have occurred
primarily in resections longer than 12 cm. An analysis of this struc
ture was undertaken to investigate the effects of graft length on the
stresses in the posterior graft. In addition, the effects of leg (tibia
plus grafts) rotation and IM rod size on the stresses in the graft and
rod were analyzed.
The physiological loads on the fused limb during normal walking were
obtained via gait analysis. Computed tomography (CT) scanning of a sham
resection arthrodesis was used to obtain the crosssectional geometry of
the structure. This information was used to develop a finite element
vii
model of the structure, consisting of 847 elements with 4700 degrees of
freedom. Parametric simulations were performed with this model, varying
graft length, leg position, and IM rod size. A computer program for the
constrained nonlinear optimization of these variables was developed and
implemented.
A stress peak was found to occur proximally in the posterior graft,
which correlated well with the predominant in vivo fatigue fracture loca
tion. This stress peak was minimized with an internal leg rotation (with
respect to the femur) of 20 degrees.
An analytical result possibly related to the clinically observed
"12 cm length effect" was an inflection on the cubic graft lengthleg
rotationposterior graft stress response surface at a graft length of
12 cm. The stresses were greater and more sensitive to leg position
in the longer grafts than in grafts less than 12 cm in length.
The largest feasible IM rod, combined with internal leg rotation,
minimized the likelihood that fatigue fracture of the IM rod or the
posterior graft would occur in the modelled structure.
viii
CHAPTER 1
INTRODUCTION
The Tibial Turnback
The tibiall turnback" resection arthrodesis is a method of surgical
intervention for patients with malignant tumors located in the distal
femur (thigh). This procedure uses autogenous cortical bone grafts to
bridge the defect created by the resection, and an intramedullary (IM)
rod to provide stable fixation (Enneking, 1982). Figure 1.1 illustrates
the resulting knee fusion structure and its principal components: the
fluted IM rod (E), which extends over the entire length of the structure
and which is supported by cortical bone proximally in the femur and dis
tally in the tibia; the anterior cortical graft (F), taken from the ante
rior cortex of the tibia; and the posterior graft (G), taken from the
ipsilateral fibula. (The reader is referred to Enneking (1982) for a
detailed description of the surgical procedure.)
This procedure has been found to provide excellent results in most
patients, from a functional as well as an oncological point of view.
Nevertheless, a survey of resection arthrodesis patients (Enneking
et al., 1980) revealed certain shortcomings which need to be addressed.
Fatigue fractures in the posterior fibularr) graft occurred in a number
of cases, seriously compromising the longterm reliability of the proce
dure. These fractures consistently occurred in the proximal portion of
the graft at or near the distal extreme of the proximal callus, as illus
trated in Figure 1.2.
Femur
Knee Level  
A N
Tibia
Figure 1.1. Schematic illustration of the tibiall turnback"
resection arthrodesis structure. Anterior graft (F) is taken from the
tibial cortex, as shown; posterior graft (G) is taken from the ipsi
lateral fibula. Both anterior grafthost junctions (A and B) and the
proximal posterior grafthost junction (D) have cortical bonetocortical
bone contact; distal posterior junction (C) consists of cortical bone of
the fibular graft resting in the cancellous bone of the tibial plateau.
Fluted intramedullary rod (E) extends over the entire length of the
structure, and is supported by cortical bone in the proximal femur and
distal tibia.
Callus
\ Fracture
Ant. Graft Post. Graft
IM Rod
S Callus
iI
STibia
Figure 1.2. Detail of healed tibiall turnback" resection
arthrodesis structure. Callus forms at the grafthost junctions of
the posterior graft, extending into the graft over a length of 34 cm.
Fatigue fractures of the posterior graft have been observed to occur
near the distal end of the proximal callus, as shown.
Eleven cases were surveyed in which the specific procedure and IM
rod shown in Figure 1.1 were used; four patients sustained graft frac
tures. In another three cases involving a similar procedure (femoral
turnback, for resection of tumors in the proximal tibia), one graft frac
ture occurred. Combining the above cases with fifteen earlier procedures
in which a weaker, more flexible diamondshaped Street or cloverleaf
shaped Kuntscher IM rod was used, a total of 29 cases yielded 13 graft
fatigue fractures for an overall failure rate of 45 per cent.
Resection (and hence, graft) length was found to have a major
influence on the failure rate of the graft. A graft length of 12 cm
or less demonstrated a fracture rate of 14 per cent (1 out of 7), whereas
grafts longer than 12 cm suffered a failure rate of 54 per cent (12 out
of 22). This "12 cm length effect" could be seen even in the limited
sample of procedures using the stronger and more rigid fluted IM rod.
No graft fractures occurred in the three cases with a graft length less
than or equal to 12 cm, whereas a failure rate of 45 per cent (5 out of
11) obtained in the longer resections.
Another factor affecting the reliability of the resection arthro
desis procedure was the type and size of intramedullary rod used to
provide fixation. The use of the fluted IM rod was prompted by the high
rate of failure (31 per cent; 4 out of 13) of the Kuntscher and Street
rods used in the initial cases. Enneking reported no failures of the
stronger fluted rod; a recent case, however, saw the fatigue fracture
of a 9 mm fluted rod in a 15 cm femoral resection.
Problem Statement
A mechanical analysis of the tibial turnback resection arthrodesis
structure was undertaken to study the above effects. Specifically, the
aims of this study were
(1) to establish an analytical basis for the clinically observed
"12 cm length effect" described above,
(2) to find the relationship between graft length and the stresses
which occur in the posterior graft,
(3) to determine the effect of fluted IM rod size on the stresses
generated in the rod and posterior graft, and
(4) to establish criteria by which the incidence of fatigue frac
tures might be reduced for this procedure.
The primary goal of this study therefore was to provide some
guidance to the surgeon, based on engineering principles, which would
necessarily take the form of specifications for parameters over which
the surgeon had some control.
The only specifications which would be
advanced would be those which could be practically implemented, or which
represented physiological conditions the surgeon might be expected to
encounter.
Nevertheless, a number of parameters were open to study. Changes
in each parameter would not only have an effect on the stresses in the
graft but, as in many other systems, would in all likelihood have an
influence on the effects of other parameters as well. Thus the effects
of the parameters would be coupled, and a thorough analysis would have to
take this into account. In addition, the restriction that proposed
specifications be limited to practically achievable values implied that
constraints would have to be placed on the parametric variables. The
problem, formulated in this manner, was seen to be a constrained optimi
zation problem, and was therefore approached on that basis.
The specific parameters to be studied here were chosen based on the
concerns expressed by the surgeons, and included those most directly
under their control in all applications of this procedure, i.e., leg
position and fluted IM rod size. "Leg position" in the present context
refers to surgically imposed axial rotation of the tibia with respect to
the femur. The anterior and posterior grafts, illustrated in Figure 1.1,
are oriented in the anteriorposterior plane of the tibia during assembly
of the arthrodesis structure. Axial rotation of the tibia thus places
the grafts in a similarly rotated orientation with respect to the femur.
The primary significance of this is that the loads applied to the grafts
vary according to the amount and direction of rotation of the tibia
graft structure. In the present work the influence of leg position on
minimizing the stresses in the posterior graft for various resection
lengths and the significance of leg position and IM rod size for avoiding
fatigue fractures of the graft and rod in longer resections were
addressed.
A second goal of this work was to advance a specific philosophy
regarding the study of biomechanical problems of this type, founded on
the principles of optimization and guided by a desire for practical
results. It is felt that the methods and routines developed here would
be directly applicable to the analysis of many important orthopedic
problems.
Methodology
Finite Element Analysis
The complexity of the structure to be studied led immediately
to the use of the finite element method. This method has been used
extensively to study stress distributions in various biomechanical
structures, particularly the hip and knee joints. An excellent review
of the literature concerning these applications may be found in Huiskes
and Chao (1983); see also, McNeice (1980). The use of this method
facilitates the modelling of complex geometries and material distribu
tions, and has been found to be the only method capable of realisti
cally modelling orthopedic systems (Rohlmann et al., 1980a).
A recurrent limitation in previous applications of the finite element
method to the study of boneimplant systems has been in the modelling of
the boneimplant interface. Previous models have assumed perfect bonding
between the bone and implant, thus transmitting tensile as well as
compressive forces across the interface. In the physiological system,
only compressive forces are transmitted across the boundary; furthermore,
Miller (1977) has shown that a finite interfacial stiffness exists in
this region, and has demonstrated its crucial role in the load transfer
which takes place between the implant and the supporting bone. The ANSYS
finite element program (Swanson Analysis Systems, Inc.) used in this work
provided a threedimensional interface element which allowed the interface
regions to be modelled in a realistic manner using physiological stiffness
values.
The geometry of the resection arthrodesis structure was obtained
via computerized tomography (CT). This method provided sequential cross
sectional views of the structure, upon which both the structural shape
and the locations of IM rod support were based. The physiological loads
imposed on the structure during normal walking were obtained using the
techniques of gait analysis.
The use of the above methods led to the development of a finite
element model of the resection arthrodesis which was as realistic as
possible within the constraints of the available computer resources.
The verification of a model usually depends on a comparison of the
model's results with those obtained experimentally using, e.g., strain
gage techniques. This type of explicit experimental stress analysis was
not possible here, since a physiologically healed structure was prereq
uisite. The model was therefore verified implicitly on the basis of its
ability to predict correctly the likelihood and location of the clini
cally observed fatigue fractures.
Optimization
The overall approach taken in this work was founded on the principles
of optimization: a study of the coupled effects of a number of parameters
on the resulting stresses in the posterior graft would be undertaken, with
the aim of minimizing these stresses. This took the form of a parametric
study wherein the effects of parametric variations on the graft stresses
were ascertained via sequential simulations using the finite element
model. A number of authors (see, e.g., Huiskes (1979), and Rohlmann et al.
(1980a,b)) have stressed that the fundamental advantage of the finite
element method in biomechanical research lies in its ability to model
accurately the qualitative behavior of the physiological structure. This
is based on the fact that the trends observed in a parametric study are
often experimentally verifiable, even though the numerical values may be
unrealistic. The use of the methods of optimization increases the utility
of the finite element method in this context as they operate on the form
of the observed response, providing verifiable information on the optimum
properties of any number of structural parameters.
The methods of optimization have been applied in several previous
biomechanical studies. Ghista and Elanjovan (1977) combined optimization
with the finite element method to design a femoral fracture fixation plate
based on minimizing the implant's weight. Lin et al. (1978) used
optimization and the finite element method in an effort to determine the
material properties of the intervertebral disc, based on minimizing the
error between the parametric results obtained with the model and exper
imental data. The studies of Crowninshield (1978) and Hardt (1978) used
optimization to investigate the indeterminate problem of specifying the
force contributions of various muscles during normal activity. However,
no study could be found in the literature which used the formal methods
of optimization to determine the optimum configuration of an orthopedic
structure based on minimizing stress levels in some region.
Overview of the Present Work
The modelling of the resection arthrodesis structure is discussed
in Chapter 2. Particular emphasis is given to the methods of data acqui
sition, and to the assumptions made during the modelling process.
The principles of constrained optimization have been only cursorily
addressed in the biomechanical literature. The choice of appropriate
optimization methods can be crucial to the success of this approach, espe
cially in the constrained nonlinear case where a number of options are
available. For this reason, a rather rigorous introduction to the prin
ciples and methods is given in Chapter 3, followed by a discussion of the
methods selected for use in this work and the routines developed to imple
ment them.
The results presented in Chapter 4 include specific recommendations
to the surgeon regarding optimum leg position and IM rod size which were
shown in the model to reduce the likelihood of posterior graft and IM
rod fatigue fractures, particularly in longer resections. In addition,
an analytical basis for the aforementioned "12 cm length effect" is
proposed.
The work concludes with a consideration of other salient parameters
and effects which might be investigated in an expanded study using the
model and routines developed here. This, it is felt, would further the
understanding of the resection arthrodesis structure, and might well
lead to greater insight into the behavior of boneimplant systems
generally.
CHAPTER 2
MODELLING OF THE RESECTION ARTHRODESIS
Introduction
A realistic model of the tibial turnback resection arthrodesis was
developed in three stages: (1) a determination of the physiological loads
on the structure, (2) the obtaining of sufficient geometrical data to
describe fully the threedimensional structure, and (3) the actual imple
mentation of these data using the finite element method. (Material prop
erty specifications for the various structural components were taken from
the literature.)
The ties between the model and the system it represents must be those
of plausible association; a tractable model must include simplifications
and idealizations. In the present case, for instance, certain simplifying
assumptions were made regarding material properties, due to a lack of firm
data to support a more complex specification as well as to the success,
demonstrated in the literature, of these assumptions in yielding experi
mentally verifiable results.
Due to the material redistribution which takes place during healing,
the graft geometries obtained initially (for a normal fibula and anterior
tibial cortex) were unrealistic insofar as the modelling of a healed
structure was concerned. Assumptions regarding the geometries of fully
incorporated grafts were therefore made on the basis of clinical and
radiographic observations.
The above assumptions, as well as others deemed necessary or
desirable for the goals of this study, will be addressed in this chapter
in the course of describing the model's formulation. The chapter con
cludes with a discussion on the choice of a primary variable of interest
for the description of graft behavior (von Mises stress). The final
result of this effort was a model which proved to be quite efficacious
in yielding valuable information regarding posterior graft and intra
medullary rod behavior, via the proposed optimization strategies.
Structural Loading
The information required to describe the loads on the resection
arthrodesis structure during normal walking was obtained at the
University of Florida Gait Analysis Laboratory. This facility pro
vided for the simultaneous acquisition of motion and ground reaction
force data which was used to describe the kinematics and kinetics of
the subject's walking pattern. The methods used to obtain these data
are well developed and have been shown to yield reliable quantitative
information (Simon et al., 1977; Tylkowski et al., 1982).
Three high speed 16 mm film cameras (Locam 2, Redlake Corporation),
run at 50 frames per second, provided the initial kinematic information.
Standard anatomical landmarks in each of the three views were digitized
(Graf Pen GP3 sonic digitizer, Science Accessories Corporation) and
stored as threedimensional coordinates in the laboratory's PDP11/34
minicomputer (Digital Equipment Corporation). Subsequent processing
of these data provided a visual representation of each 1/50 second of
gait cycle in the form of stick figures which were used to determine the
threedimensional position of the leg at each interval of the gait cycle.
The dynamic ground reaction forces which acted on the supporting
limb during the stance phase of gait were obtained using a force plate
(0R63 dynamometer, Advanced Mechanical Technology Incorporated) set in
the floor of the laboratory. The three orthogonal components of the
ground reaction force vector were measured directly at a sampling rate
of 500 Hz. The resultant vector's magnitude and orientation were com
puted from these components; its origin was found via "apparent moment"
data which provided additional information sufficient to compute the
point of application of the vector on the force plate. In addition to
plots and printouts of the component force values, a vectorial repre
sentation was overlaid on the kinematic stick figure plots from which
the force vector's relationship to the supporting limb throughout the
gait cycle was ascertained. This allowed the quasistatic axial and
bending loads on the arthrodesis structure to be calculated.
The above methods were used to evaluate three resection arthrodesis
patients. The stick figures were used to obtain the force vectorfused
limb (tibia) relationships in the sagittal and frontal planes. The
relationships between the leg and the applied load obtained for Subject 1
at 15 per cent of the gait cycle (see Figures 2.la and 2.1b) were used
in the subsequent finite element study; the three force vector component
values were obtained directly from the electrical data file (rather than
measuring the illustrated vector) to improve accuracy. Based on these
values and the above relationships, the ground reaction force vector
components were transformed into the three tibial component directions,
viz., axial, anteriorposterior (AP), and mediallateral (ML), for
application to the model.
Figure 2.1. Stick figure representation of a resection
arthrodesis subject, with overlay of ground reaction force
vector. (a) Sagittal Plane. Note that the applied flexing
moment about the fused knee is general throughout the gait
cycle except at foot strike (0 percent of gait cycle).
Forces applied to the FEM model were based on data obtained
for this subject at 15 per cent of gait cycle.
4ILOHWS SODW4 AR PECENIFSE3
:94R
/
2 :.
4.
THETA W6.
PHI SI.
DREGA 1.8
2I
II!
I
'94'~
I
29
55<
30
594
$23.T?
I \
I'
1
Figure 2.1continued. (b) Frontal Plane. Note that the
applied varus moment about the fused knee in this plane is gen
eral throughout the gait cycle except at foot strike (0 per cent
of gait cycle). Forces applied to the FEM model were based on
data obtained for this subject at 15 per cent of gait cycle.
&NHS!FS 31OWI RE. Ptt tNIFS1 5 SZ.7 7 TR *.
1HCIA fH.
I \t I
/ \ /1 ;: A Y
5 .
i i. I ; ...
sis
so. ss
The resulting axial force and ML force in the coronal plane were
near their maximum values for this subject at 15 per cent of the gait
cycle. At this point in the gait cycle, the AP force in the sagittal
plane was approximately 70 per cent of the maximum value which occurred
at 45 per cent of the gait cycle. Thus a rather conservative (i.e., low)
AP load was applied to the model. This led to conservative estimates,
quantitatively, of the effects of bending loads in the saggital plane
on the magnitudes of the stresses developed in the posterior graft.
It was assumed that the positions of the grafts were determined
by specific locations on the tibial plateau; thus the orientation of
the grafts with respect to the tibia and foot was assumed to be fixed
and invariant with axial rotation of the combined foot, tibia, and graft
portions of the arthrodesis structure. The effect of leg position
(i.e., axial rotation of the tibia and grafts with respect to the
neutrally positioned femur) on graft stresses was therefore modelled
by performing an equal but opposite rotation of the AP and ML force
vector components about the axis of the leg (see Figure 2.2). Internal
leg rotation was modelled by external force vector rotation, and vice
versa.
The type of loading imposed on the leg at 15 per cent of the gait
cycle (varusflexing bending moment about the fused knee) was general,
save at footstrike (see Figure 2.1). In addition, this relationship
held for the other two subjects which were evaluated. The loads
applied to the model were therefore qualitatively representative of
those which would be applied in vivo for any resection arthrodesis
patient with this general gait pattern. The computed force vector
Ext. Rotation Int. Rotation
jI
Ground Reaction
Force Vector
(Orthogonal Component)
Figure 2.2. Relationship between orientation of the applied load
and leg position. External leg rotation increases the loading in the AP
plane, while internal leg rotation increases ML plane loading. Varia
tions in leg position were modelled by an equivalent (opposite) rotation
of the force vector, i.e., internal leg rotation was modelled by external
rotation of the applied load and vice versa.
Table 2.1
Load Vector Component Values (Transformed) of Subject
of Figures 2.la and 2.1b at 15 percent of Gait Cycle,
Neutral Leg Rotation
Axial Force ML Force AP Force
510. N 60. N 70. N
(Body Weight = 520. N)
component values which were derived for the subject of Figure 2.1 and
subsequently applied to the model are listed in Table 2.1.
System Geometry Data Acquisition
The tibial turnback resection arthrodesis procedure was performed
on a fresh cadaver of a stature similar to the subject of Figure 2.1,
using a 12 mm fluted rod and a 12 cm resection length. A Philips com
puted axial tomography (CT) scanner was employed to obtain the bone
crosssectional geometry. The fused limb was CT scanned at 6.2 mm
intervals over its entire length, for a total of 132 transverse sections.
The CT scan images were optimized for corticalcancellous bone
contrast (Baxter and Sorenson, 1981; Koehler et al., 1979). A full scale
window width (CT number = 3200), centered at the average radiographic
density of the cortical and cancellous bone (CT number = 900), minimized
the error in reproducing the loadbearing cortical thickness in each
section. Possible error resulting from Xray beam hardening due to the
centrally located IM rod was found not to be a factor with these settings.
Verification of the image geometry was accomplished by cutting
several windows in the femoral cortex, measuring with a ruler the cor
tical thickness at those points, and comparing the measured values to
the computed image values at the same locations. Figure 2.3 illustrates
the computed image values for cortical bone thickness at two locations
(4.0 and 7.7 mm); the thicknesses as measured with a ruler were 4 and
8 mm, respectively. Due to the gradual transition from cortical to
cancellous bone, the best precision to which cortical thickness could
be measured was found to be approximately 1 mm. This agrees with the
Figure 2.3. CT scan image of femoral crosssection. (a) Indicated
(AB) dimension of anterior cortical thickness at this location in the
scanned femur was measured to be 4.0 mm; calculated value (above) based on
scan image equals 4.0 mm (1.3x magnification in source image).
Figure 2.3continued. (b) Indicated (AB) dimension of lateral
cortical thickness at this location in the scanned femur was measured
to be 8.0 mm; calculated value (above) based on scan image equals 7.7 mm
(1.3x magnification in source image).
uncertainty reported by Rohlmann (1980a). The scan image settings used
here for obtaining the crosssectional geometry of cortical bone
thus provided good accuracy and precision, within the limits of uncer
tainty imposed by the corticalcancellous transition region, and
appeared to avoid successfully any deleterious effects due to beam
hardening. The results obtained were undoubtedly enhanced by the high
resolution (6.2 mm section thickness) used for the scans.
CT scan images at an actual magnification of 1.3x were obtained
for all sections using the settings described above. Observations of
the salient bony features of the tibia led to a longitudinal geomet
rical development based on every other scanned section, i.e., every
12.4 mm, an interval which still retained sufficient detail of all
prominent tibial features. Note, however, that although the longitu
dinal tibial geometry consisted of linear extrapolations over 12.4 mm
intervals, the individual sections were based on CT scan images which
represented an averaging of Xray data over just 6.2 mm. This strat
egy preserved the salient structural details while keeping the model's
size to within tolerable limits.
The appropriate CT scan images were traced and discretized into
radially disposed elements according to the considerations discussed
in the next section. The discretized section tracings were then digi
tized to obtain the threedimensional nodal coordinate values for the
structure. Longitudinal development was effected by keeping all sec
tions in proper registration, using the IM rod's centerline to provide
translational (AP and ML) alignment. Rotation of the structure was
of course held constant throughout the imaging process so that
sectional image rotational alignment only required that the edges
of successive images be kept parallel.
Finite Element Modelling
The portion of the ground reaction load carried by the grafts is
transmitted to them by the tibia's cortical shell. Anterior graft load
ing is transmitted along the tibial cortex directly to the anterior cor
tical graft. No such cortical bonecortical bone contact exists, how
ever, at the distal posterior graft junction (location "C" in Figure 1.1).
Here the cortical graft is implanted in the cancellous bone of the tibial
plateau, so that the loading path traverses a region of cancellous bone en
route from the tibial cortex to the posterior graft. Proximally, full
corticalcortical contact exists at both graft junctions with the femur.
Further, the cortical shell thickness of the femur is much greater than
that of the tibia, which translates into much greater structural
stiffness.
Some of the applied bending load is carried directly from the tibia
to the femur via the IM rod; the transfer of this load to and from the
rod is effected at points of cortical intramedullary support, distally
in the tibia and proximally in the femur.
The observations discussed above led to two assumptions regarding
modelling of the proximal portion of the structure. The first assump
tion was that the proximal grafthost junctions could be realistically
modelled by "fixing" the proximal ends of the grafts, using displace
ment constraints to prevent both translation and rotation of these
modes. The second assumption was that the proximal IM rod support
could be modelled with appropriate displacement constraints on the rod
(allowing only axial translation) where cortical support was seen to
occur. Additional justification for the latter assumption was that
its effects would be limited in scope to the proximal femur (according
to St. Venant's principle), and therefore of no consequence in this
study.
It may be seen, however, that no such assumptions could be made
in the distal portion of the structure. Lying "upstream" in the ground
reaction force path, the significance of the complicated geometry of the
tibia (including the anterior defect at the host site of the anterior
graft bone) as well as of the possible changes in distal cortical rod
support with structural deformation could not be predicted. An effort
was therefore made to model the tibia as accurately as possible, con
sistent with computer memory, run time, and storage limitations.
Each transverse section of the tibia was discretized by drawing
a series of radial lines on the CT scan image tracings, intersecting
at the intramedullary rod's centerline. The maximum angle of sweep
between these radial lines which would preserve the salient geometri
cal features of the tibia's bony structure was observed to be approx
imately 20 degrees. It has been suggested, however, that an arc of
about 15 degrees should be used when modelling cylindrical structures
with straightline elements (DeSalvo and Swanson, 1983). The above
stated limitations led to a compromise and the use of an average angle
of 18 degrees per element, for a total of 20 elements per transverse
section for the tibial cortex.
Rohlmann et al. (1980b) have demonstrated the deleterious effects
which irregular element shapes can have on calculated stress distribu
tions. To the extent that these effects changed the force transmission
characteristics of the structure, the resulting loads on the grafts and
the IM rod could be artifactually altered. Attention was therefore paid
27
to keeping the element distribution as regular as possible from section
to section, and a limit of 23 degrees was placed on deviations from the
18degree average elemental sweep angle to accommodate bony features.
It has been shown (Hayes et al., 1978) that the function of cancellous
bone is to transfer and distribute the applied load to and from the corti
cal shell. In an axisymmetric finite element model of the proximal tibia,
Hayes et al. showed that the transfer of stresses from the articulating
surfaces to the tibial cortex is essentially complete proximal to the
tibial shaft per se. In the present work, the cancellous bone of the
proximal tibia served to transfer the applied load from the tibial cortex
to the posterior graft. The proximal 37.2 mm of cancellous bone (wherein
the load transfer would be expected to occur) was therefore modelled in
four concentric layers with the same 18degree average elemental sweep
described above, for a total of 80 elements per transverse section. Once
again attention was paid to maintaining a regular element arrangement from
section to section, especially in the proximal inner layers where direct
load transfer to the posterior graft took place. In the remainder of the
model, the structural contributions of cancellous bone would be insignif
icant and were therefore neglected. Graft and IM rod element nodes were
located at the centroids of their respective crosssections.
The threedimensional model geometry thus obtained was manually
entered into the ANSYS preprocessor in the form of nodal coordinate
values, resulting in the nodal configuration shown in Figure 2.4.
Threedimensional isoparametric solid elements (ANSYS STIF 45; see
DeSalvo and Swanson (1983) and Kohnke (1983) for a complete descrip
tion of this and other cited element types) were used to model the
PEP7 INP=
,.'
3 ..
]
F ]i :
PESECTIOI HFPTHFPIE'S IS
HI IS','S
12 1:: 84
20. ,34
FFEP7 HOrlES'
LITII S ALIri
::I 1
D = 1
['IT=34.6
: :F=4. ?C
'.'F =37.4
7F =4. "'
Figure 2.4. Computer plot of nodes used to define the three
dimensional geometry of the resection arthrodesis structure during
the finite element modelling process. Geometry was obtained from
sequential computed tomography (CT) scanning of the structure (see
text). (Z axis is directed from medial to lateral side, X axis
from anterior to posterior).
'PEP? I IP= '.'.
F'F'EP? ELEIHE 1HT
TSET
lUTO SCHLIiF
;.U=1
DIST=17.5
T:F=4.31
I'F=13.7
ZF=4.41
HIDDEN
I
I PESECTIOInI HPTHPF'ODESI'. : 1
Figure 2.5. Computer plot of elements used to model the tibial
portion of the resection arthrodesis structure. Note anterior defect
created by removal of graft bone.
cortical and cancellous bone of the tibia. Figure 2.5 illustrates
the tibial portion of the finite element model.
The 12 mm fluted IM rod was modelled as a right circular cylinder,
using ANSYS STIF 9 cylindrical elements and geometrical information
derived from Miller (1977). The effect of the flutes on the rod's
bending stiffness was included by specifying a cylindrical geometry
(i.e., outer diameter and wall thickness) which provided an area moment
of inertia equal to that given by Miller for a fluted rod of similar
nominal diameter. The proximal curvature of the rod was neglected as
it was felt to be of no structural consequence (its sole function is
to mimic the curvature of the femoral shaft, thus facilitating rod
insertion).
It is known that the shapes of the grafts are altered during the
course of the healing process (H. Burchardt, Personal Communication,
Aug. 1984). A direct transfer of the CT scanderived geometry of the
unincorporated fibular and tibial grafts would therefore not be real
istic. Radiographic observations of fully incorporated grafts indicated
that they take on a much more regular shape along their length, the ante
rior graft becoming more prismatic and the posterior graft, more cylin
drical.
Wolff's law (see, e.g., Fung (1981)) asserts that changes in the shape
of bone may occur as a response to the stress (and/or strain) acting in
it. Based on the assumption that the remodelling in this case might, at
least in part, effect a change in graft shape to minimize the stresses at
a given location, a series of graft cross sections were analyzed using
the SCADS computer program (Piotrowski, 1971) to determine which partic
ular cross sections exhibited the least bending stress in response to
a given load. The geometrical properties of these cross sections were
then used as a baseline for further development.
Anteriorly, it was observed that while the width (ML dimension)
of the graft undergoes little or no change, a slight thickening (AP dimen
sion) seems to occur along the full length of the graft. Therefore, the
base line anterior graft geometry was changed to reflect this by adding
10 per cent to the graft thickness. This graft was then modelled using
ANSYS STIF 4 threedimensional beam elements, using the geometrical
information described above and listed in Table 2.2.
The healed posterior graft exhibits callus formation near its
proximal and distal junctions, appearing not unlike that which forms
during fracture healing. This callus has material properties similar
to cortical bone (H. Burchardt, Personal Communication, Aug. 1984),
and extends into the graft over a length of 34 cm, fairing in grad
ually over the inner 12 cm. It has been observed that a diametral
increase of about 25 percent (over the baseline graft) obtains along
the thicker portion of the callus, with no appreciable change in the
graft's wall thickness (Springfield, 1984).
Guided by the above considerations, the posterior graft was
modeled using ANSYS STIF 9 threedimensional cylindrical elements.
The baseline geometry described above was used for the midsection.
The proximal and distal 3.6 cm were modelled as callus, with a 25 per
cent increase in diameter over the outer 2.4 cm, tapering down lin
early over 1.2 cm to the baseline graft diameter as illustrated in
Figure 2.6. The geometrical information used to describe this graft
is listed in Table 2.3.
Table 2.2
Anterior Graft Geometrical Properties
Base Line Graft:
Section Area
0.31 cm2
I(ML)
0.01 cm4
I(AP)
0.05 cm4
Remodelled Graft:
Section Area
0.52 cm2
I(ML)
0.03 cm4
I(AP)
0.13 cm4
I(AP) and I(ML) are the area moments of inertia
axes, respectively; W(ML) is the graft width (ML
is the graft thickness (AP dimension).
about the AP and ML
dimension), and T(AP)
W(ML)
1.38 cm
T(AP)
0.77 cm
W(ML)
1.40 cm
T(AP)
0.85 cm
.'_
Figure 2.6. Schematic illustration of posterior graft element
geometry. Callus (see Figure 1.2) was modelled in a stepwise fashion
(see Table 2.3 for defining specifications).
CALLUS 3
CALLUS 2
CALLUS 1
GRAFT MIDSECTION
(BASE LINE)
CALLUS 1
CALLUS 2
CALLUS 3
6*
Table 2.3
Posterior Graft Geometrical Specifications
Outer Diameter Wall Thickness
Base Line Graft: 0.92 cm 0.15 cm
Callus 1 (+ 8.3 per cent): 1.00 cm 0.15 cm
Callus 2 (+ 16.7 per cent): 1.08 cm 0.15 cm
Callus 3 (+ 25.0 per cent): 1.16 cm 0.15 cm
Base line graft dimensions apply to graft midsection; Callus 3 refers
to the callus at the outermost 2.4 cm of the graft; Callus 1 and 2 are
the linear extrapolations for tapering of the callus between Callus 3
and the base line graft.
A nodal spacing, and hence, element length, of 0.6 cm was used for
the grafts to provide a reasonably fine stress distribution resolution,
giving a total of 30 elements for each 18 cm graft. Graft length varia
tions during the course of parametric studies were made in the midsection,
preserving the callus elements at the extremities of the posterior graft.
Finally, it should be noted that the posterior graft was extended 1.2 cm
into the cancellous bone of the tibial plateau, corresponding to the pre
scribed surgical procedure and thus modelling as realistically as possible
the actual load transfer path between the graft and the tibia. To avoid
confusion, however, "graft length" and "resection length" will here be
used interchangeably (as they are by the clinicians), and the imbedded
portion of the posterior graft therefore omitted from the graft length
specifications given in this work.
The BoneImplant Interface
The CT scan images of the sham resection arthrodesis described
previously were used to determine the regions of cortical bone support
for the IM rod. Proximally these regions were modelled by placing
appropriate displacement constraints on the rod, as discussed earlier.
Axial translation of the rod was permitted; AP and ML displacements
were restricted according to the support seen to occur in the femur.
Cortical support of the IM rod in the tibia was modelled using
ANSYS STIF 52 threedimensional interface elements at locations of
observed bonerod contact. These elements allowed the rod to slide
on, and separate from, the supporting cortical bone while transmit
ting only compressive forces. The defining specifications for these
elements included a 0.5 mm interference fit (before load application,
corresponding to the actual fluted rodintramedullary canal reaming
size specifications), and a compressive interface stiffness value of
6.0 x 106 N/m (Miller, 1977). A total of 32 interface elements in
seven adjacent cross sections were used to provide the observed distal
IM rod support. Figure 2.7 illustrates the IM rod with its associated
supporting elements and constraints.
Material Properties
Bones are themselves structures, consisting of both cortical and
cancellous bone. Although bone exhibits nonhomogeneous, anisotropic
behavior, a complete and consistent characterization of these effects has
not yet been made. It has been shown, however, that assumptions of
homogeneity and isotropy can yield experimentally verifiable results
in the analysis of whole bones (Huiskes, 1982). This, in addition to
the need to restrict the complexity of the model to a level consistent
with the available computer resources, led to the invoking of these sim
plifying assumptions.
Values for the material properties of cortical and cancellous bone
were taken from the literature and are listed in Table 2.4. The stiff
ness value used for cancellous bone deserves comment: reports in the
literature have demonstrated consistently a locationdependent range of
values for this parameter. The value chosen is consistent with the
reported values, and in particular corresponds to the modulus given by
Goldstein et al. (1983) for the region of the tibial plateau which in
the present case supports the posterior graft.
The material properties for the (steel) IM rod were taken from
Miller (1977) and are listed in Table 2.4.
'REPF INP= 12 1 : :4
21.3031
PPEP? ELEMEI TS
TSET
T5ET
TDE C= 1
PDEC=1
AUTO SCALING
:.JU= 1
I'3T=23.3
T:F= 3. 2.
YF'=42.6
ZF=3.8S
RE'ECTIOII HPTHPODESIS 3.1
Figure 2.7. Computer plot of intramedullary (IM) rod elements and
associated support. (a) Displacement constraints were imposed proximally
in regions of cortical support of the IM rod as observed on CT scan images.
'PEP? IHP= 12 1 :4
2 1 .7347
PPEPT ELEIHEIT'S
T'SET
E MA::= 3 40
ITlO SC' LItG3
YI 14
PESE.,TIOH HPTHPOIES= 1.1
Figure 2.7continued. (b) Detail of Figure 2.7a, illustrating the
interface elements used to model the distal tibial cortical support for
the IM rod. Support was provided based on CT scan image observations.
Table 2.4
Material Property Specifications
E(Young's
Cortical Bone: 18.
(After Yamada and Evans, 1970;
Cancellous Bone: 0.2
(After Carter and Hayes, 1977;
1982; Goldstein et al., 1983)
Fluted IM Rod: 200.
(After Miller, 1977)
Modulus) v(Poisson's Ratio)
GPa 0.3
Fung, 1981)
GPa 0.2
Hayes et al., 1978; Williams and Lewis,
GPa
0.3
F'EPF IHIF'=
1 1:: 4
c0.'?0Uh1
PPEP7 ELEHEHITS
TDE := 1
PDE _= 1
FE'C=1
HUTO ',: HLIIHI'
' I 1
TI I '=
,'F = i
::F=4. 1
PESE':TIOI AHPTHPCIDESIS .'
Figure 2.8. Computer plot of finite element model of the tibiall
turnback" resection arthrodesis structure. Complete model consisted of
847 elements, with a total of approximately 4700 degrees of freedom.
The Finite Element Model
The complete finite element model of the tibial turnback resection
arthrodesis is illustrated in Figure 2.8 and consisted of 847 elements
with approximately 4700 degrees of freedom. A total of 26 runs in two
series were made to obtain the results described in Chapter 4. The model
was a nonlinear one due to the bilinear interface elements, and there
fore required an iterative approach; each run took from 4 to 12 hours
of CPU time on a VAX 11/750 minicomputer, depending on the number of
iterations required to achieve interface element stabilization.
Some Final Considerations
A finite element analysis may invariably be depended on in one
respect: to produce a surfeit of data. Discrimination based on knowl
edge of the actual physical system's behavior as well as on the goals
of the study must be used to choose a primary variable of interest.
In the present study, the primary regions of interest were the
posterior graft and the intramedullary rod; the overall model was
designed to represent the loads on these components in a realistic man
ner. The goal of this modelling effort was in effect to build a "black
box" by which the influence of changes in certain parameters (e.g., leg
position and IM rod geometry) on the likelihood of graft or rod fracture
could be elucidated via the optimization methods described in Chapter 3.
This approach was predicated on the assumption that the "state of stress"
at a point in the graft or rod could be plausibly represented by a single
value for each combination of the parameters (design variables). The von
Mises stress meets this requirement, taking into account the effect of all
three principal stresses to arrive at a single value which may be said
to represent the state of stress at a point. The von Mises stress has
been used in this context in several previous biomechanical finite
element analyses (see, e.g., Oonishi et al. (1983), Huiskes (1980),
Crowninshield et al. (1980), and Harris et al. (1978)).
A second consideration was that of the known mode of failure of the
graft: fatigue fracture. It has been posited (see, e.g., McClintock
(1966)) that the nucleation of a fatigue crack is due primarily to shear
stress. Furthermore, Reilly (1974) has shown that fractures in bone are
initiated by shear stresses. Thus, the criterion by which failure of the
graft might be predicted must take into account the shear stresses at
that point. Although calculated on the basis of the three principal
stresses at a point, it may be shown that the von Mises stress is actu
ally a function of the three major shear stresses at that point (Burr,
1981). Finally, the von Mises stress is the basis for a welldeveloped
theory of failure based on a consideration of the distortion energy in an
isotropic material. It is therefore reasonable, under the prevailing
assumption of isotropic behavior, to expect that the von Mises stress
could serve well as a predictor of the likelihood of fatigue fracture
in the posterior graft and the IM rod.
On the basis of the above considerations, the von Mises stress was
chosen to be the primary variable of interest in this study for analyz
ing the mechanical behavior of the posterior graft and IM rod, and is
the basis upon which the results in Chapter 4 will be presented.
CHAPTER 3
OPTIMIZATION THEORY
Introduction
The overall objective of this work was to find ways by which locally
high stresses in the graft portion of the resection arthrodesis structure
could be minimized. The finite element model discussed in the preceding
chapter was designed to allow for variations in a number of structural
parameters such as leg position (rotation) and IM rod geometry (size),
with the goal of finding the specific values for these parameters which
minimized the stresses in the posterior graft.
The optimization methods and routines to be discussed in this chapter
enhanced the utility of the parametric finite element analysis in two
ways. First, they facilitated the planning of future experiments (com
puter runs) based on the trends observed in previous results. This
allowed for the most efficient use of resources while maximizing the
amount of useful information obtained from the entire series of exper
iments.
These methods also expanded the capabilities of the parametric
finite element analysis by allowing the coupled effects of any number
of parameters on posterior graft stresses to be studied. In the pres
ent case two parametric or design variables, i.e., leg position and
IM rod geometry, were studied to determine their relationship to the
maximum von Mises stress in the posterior graft. The graphical
representation of this relationship was therefore a threedimensional
43
surface with the two design variables on the horizontal axes and the
graft stress represented by the surface topology. A number of other
design variables, discussed in the final chapter of this work, remain to
be investigated. It is clear, however, that the addition of just one more
variable would yield a fourdimensional hypersurface which could not be
visualized readily: formal optimization methods must be employed in this
case.
The routines necessary to accomplish this expanded analysis were
therefore written and implemented at this stage in order to (1) demon
strate their usefulness, and (2) provide a complete foundation upon which
future research might be built.
This chapter will begin with the formulation of the present study
as an optimization problem. The methods required to solve this problem
will then be selected and discussed.
The Optimization Problem
The formulation of the optimization problem may be summarized in
the following seven steps:
(1) Define a suitable objective for the problem of interest.
The stated objective of this study was to find ways by which locally
high stresses in the graft portions of the resection arthrodesis struc
ture could be minimized.
(2) Examine the system to be evaluated to determine the inter
relationship of the various system elements. In the present case the
complexity of the apparent interrelationships led to the use of the
finite element method. This allowed for the conceptualization of
a model that would not require the more drastic and overt simplifications
of a purely analytical model.
(3) Construct a model for the system which will allow the objective
to be formulated in terms of the design variables. We have seen in
Chapter 2 how the finite element method could be used to obtain a model of
this very complex system. The result was a model with which the effects
of changes in various system parameters (design variables) on the graft
stresses could be experimentally evaluated.
(4) Define the restrictions which must be placed on the design
variables. Restrictions were placed on the design variables based on
feasibility both from a physiological and a technical point of view:
leg rotation was restricted to 20 degrees from neutral, both internally
and externally, while IM rod size was constrained to the available range
of nominal diameters, i.e., 9 mm to 14 mm.
(5) Carry out a simulation by expressing the objective in terms of
the design variables, using the model developed in (3) above. This
expression is the objective function. A series of parametric finite
element runs was performed, varying the values of the design variables
within the feasible region (defined in (4), above) and observing the
effect on stress peaks in the graft. In this manner a number of exper
imental data points relating the design variables to the stresses in
the graft were obtained. A least squares fit of a quadratic surface
to these data points resulted in a formal mathematical expression, the
objective function, with which numerical optimization could proceed.
(6) Verify that the proposed model does in fact represent the
system being studied. It is unrealistic to suppose that a model could
be developed which would simulate precisely every nuance of the physio
logical structure. In the present study, verification of the results
was made on the basis of the model's ability to predict accurately
failure of the grafts in the manner, and under the conditions, observed
clinically.
(7) Determine the optimum solution for the system, subject to
the restrictions defined in (4), above. That is, find those values of
the design variables which will minimize the value of the objective
function, restricting the variables to values within the feasible region
defined in (4). The restrictions on the design variables and the qua
dratic form of the objective function meant that the methods of con
strained nonlinear optimization would be required. An examination of
the available methods led to the use of the Hooke and Jeeves procedure
for the unconstrained optimization, combined with the transformation of
variables method for handling the constraints and the random restart
method for global optimization.
This process yielded the specific values of the design variables,
within the feasible region, which minimized the peak posterior graft
von Mises stresses in the model, thus minimizing the risk of graft failure
under physiological loading conditions. The methods implemented to accom
plish the above will be described in the following sections.
Objective Function Formulation
We are seeking here a method of obtaining a relationship between the
design variables X (a vector) and the objective criterion f(X), here
designated to be the von Mises stress peak in the graft. The design
variables may be n in number, i.e.,
x(1)
x(2)
X*
x(n)
We will restrict ourselves in this discussion to n = 2, for two reasons.
First, it will be possible to visualize the relationships if we restrict
the overall dimensionality to three (i.e., a twodimensional design space
yields a functional surface in three dimensions). Second, we will be con
cerned in this work with the effect of two design variables (load orien
tation and intramedullary rod geometry) on the graft response (von Mises
stress). The algorithms developed here are nevertheless general, in that
an extension to n dimensions will only require that more experiments be
carried out. This is important because we are seeking here not only to
solve a specific problem, but also to develop a general philosophy
which may be applied to other similar biomechanical investigations.
The process we will use to obtain the above relationships is called
experimental response surface fitting. This procedure seeks to represent
a set of "black box" outputs in terms of the inputs by a mathematical
relationship. In our case, the "black box" is the finite element model;
the inputs are the various values for the design variables, and the
outputs are the predicted responses at the graft, i.e.,
load orientation [x(1)] FEM maximum
model posterior graft
IM rod geometry [x(2)] + "Black Box" stress f(X)
obtained by repeated finite element runs using different design variable
sets. It can be seen that this will result in a set of experimental data
points in three dimensions. The mathematical representation of this sur
face will be considered valid if the approximation is of sufficient qual
ity to preserve the essential features of the experimentally observed
response surface.
A least squares fit of a general quadratic surface was used to
approximate the experimental response surface in the region of interest,
which was restricted to physically and physiologically feasible values for
the two design variables. Two limitations of this mathematical model
should be noted. First, the approximating surface was based on data
acquired at discrete points. An approximating function produced in
this manner can be depended on to interpolate accurately functional
values between data points in the region of interest. Experiments have
shown, however, that the derivatives of the approximating function can
be unreliable if used explicitly to approximate the local form of the
actual surface. Optimization methods which require explicit evaluations
of objective function derivatives should therefore be avoided in this
instance.
Secondly, the experimental data were gathered over a limited portion
of the design space. There is no basis for assuming a more general appli
cability for the approximating function outside this region. Indeed, the
evaluation of this function outside these defined limits should be
avoided. This is not necessarily a limitation, however, as the limits
were chosen to represent the range of feasibility for the design vari
ables and it would be counterproductive to seek an optimum which couldn't
be achieved in the actual system. We will therefore strictly constrain
our search for an optimum design to this bounded region of feasibility.
We are now in a position to state formally the constrained nonlinear
optimization problem in terms of the objective function and its associated
constraints:
Find X = {x(1), x(2)} [transpose]
2 2
To minimize f(X) = al + a2x1 + ax2 + a4x1 + a5x2 + a6x1x2
Subject to gl(X) = A < x(1) s B
g2(X) = C < x(2) < D
Where x(1) = load vector orientation
x(2) = IM (intramedullary) rod geometry (~EI)
A, B = physiological limits on foot rotation
C, D = technical feasibility limits for IM rod.
The question of how to find the minimum(s) of this objective function
subject to the above constraints will now be addressed.
Optimizing a Constrained Nonlinear Function
Many numerical optimization methods are available for the general
ndimensional constrained nonlinear problem, and all have their advocates.
No one method can be singled out as best for all cases, and each must be
evaluated in relation to the problem under consideration. All utilize
an iterative sequential search of the analytical surface involving the
successive calculation of new values of the objective function and the
comparison of these values with the "best" value obtained up to that
point. The basic differences in the methods are in the philosophy dic
tating the direction of search, and in the determination of if and when
that direction should be changed in the pursuit of the minimum objective
function value. It is in the choosing of directions of search that a
method is efficient, converging rapidly to a solution, or inefficient.
This numerical approach is analogous to starting off in a mountain
valley, blindfolded. The search will involve a climb to the highest peak
which must be achieved by repeated test moves, following each with an
altimeter reading. Now it will be seen that random moves would get you
there eventually, but a strategy which allowed you to use previous expe
rience to modify your search direction, i.e., an extrapolation approach
of some sort, would be helpful. (This analogy is based on finding a max
imum because the reverse doesn't work so well, gravity being what it is.
On the other hand, one can readily see that finding a maximum is the
same as finding a minimum, turned upside down. This would be formally
accomplished by multiplying the objective function by (1), i.e.,
minimum [f(X)] = maximum [f(X)].)
Another problem which will need to be considered is the following.
Suppose you've reached a peak. How do you know it is the highest one?
(Remember, you're still blindfolded). Your only recourse is to (1) walk
around a bit and see if you go any higher, then (2) start again from a
different location and see if you end up in the same "highest" place.
Do this a number of times, and when you always end up calling the same
peak "highest," call that the "global optimum." Now remember we have
fences (constraints) to deal with, too. Our search method must allow us
to follow a fence if that is the best direction, and must be able to
recognize a high point along the fence even if the actual peak is on the
other side of it. Finally, we have no knowledge whatsoever regarding the
terrain on the other side of the fence (it may be a sheer drop for all
we know). To avoid unpleasant surprises we will therefore wish to
confine our activities to the area enclosed by the fence (i.e., the
feasible region).
Constrained Optimization
Techniques for numerically searching a constrained function may be
divided into two classes: direct and indirect. Direct methods attempt to
deal with active constraints in a direct manner, either by rebounding and
continuing in the feasible region or by following them. Indirect methods,
on the other hand, seek to redefine the objective function in such a
way that the new function has an unconstrained optimum at the same point
as the optimum of the given constrained problem; here the constraints are
used implicitly in an unconstrained search.
Direct Methods
Well developed methods for direct optimization of functions with
linear (Generalized Reduced Gradient or GRG (Wolfe, 1967; Beale, 1967),
Rosen (Rosen, 1960), Complex (Box, 1965)) and nonlinear (GRG (Abadie and
Carpentier, 1969), Complex) constraints are available. Generally they are
either very complicated to implement and suitable primarily for large
problems with few nonlinear terms (GRG (Beale, 1969)), or quite ineffi
cient near an optimum and prone to premature termination in deep narrow
valeys (Complex (Guin, 1967; Beveridge and Schechter, 1970)). Further
more, the Complex method can be prone to early termination along a
boundary (Reklaitis et al., 1983). In all cases, nonlinear constraints
can require many changes in search direction which can overwhelm com
pletely any efficiency in the line search and lead to very slow progress.
In addition, the more highly developed and very efficient unconstrained
search techniques cannot be used. Nevertheless, these methods do attempt
to find the minimum of the true function without recourse to distortion
(as in the penalty methods, below) which may be an advantage for certain
problems (Box, 1965; Davies and Swann, 1969).
Indirect Methods
Indirect methods for handling constraints include a number of
techniques, all of which modify the original objective function to
"internalize" the constraint conditions. There are two ways to accom
plish this: penalty functions and transformations of variables.
The driving motivation here is that the field of unconstrained numerical
search methods is highly developed, and many efficient algorithms are
available including several which are more powerful than the
limited number available for direct search of constrained nonlinear
functions (Box, 1966). Thus it is the prevailing view that if constraints
can be removed from the explicit formulation, this should be done (Davies
and Swann, 1969; Beveridge and Schechter, 1970).
The penalty function method combines the original objective function
and its associated constraints in a single modified "penalty function"
which may be handled by unconstrained search techniques. However, most
penalty methods require function evaluations outside the feasible region.
This is a decided disadvantage in the present case, as we wish to search
only within the feasible region due to limitations on our knowledge of the
experimentally derived response surface. Further discussion on this
method may be found in Reklaitis et al. (1983).
The transformation of variables method is fundamentally a means of
restating the original problem in such a way that the constraints are
implied. In this way a constrained optimization problem can sometimes
be reduced to a form in which no constraints appear explicitly, so that
the powerful methods of unconstrained optimization may be used. Thus the
motivation is the same as for the aforementioned penalty method, but the
approach is different.
The transformation method has a number of advantages over the penalty
methods, from the standpoint of overall accuracy as well as efficiency and
ease of implementation. These advantages are most pronounced in the class
of problems in which the constraints consist of bounds on the design var
iables, as in the subject problem of this work.
Davies and Swann (1969) endorse the use of transformations and
Powell (1969) strongly favors the method. Ibbett (1969) ran a sequence
of tests comparing a penalty method (constrained Rosenbrock) with the
transformation method (using the unconstrained Rosenbrock search tech
nique (Rosenbrock, 1960)), and found that the latter yielded a better
final solution. Box (1966) performed a more extensive comparison
between the transformation method and the penalty method of Rosenbrock
(Rosenbrock and Storey, 1966), the Complex method (Box, 1965), and a
modified Rosenbrock method (RAVE, described in Box, 1965). In all six
test problems the transformation method yielded the most accurate
result in the fewest number of function evaluations. Further, Box states,
The advantage of these transformations is that they have been
found to result in the correct solutions being obtained easily
for problems for which alternative methods made only slow
progress, or ceased to make any progress whatsoever once one
or two constraints were reached, even when the current point
was still a long way from the optimum. (p. 73)
Box has used the method with success for problems of up to twenty
variables, and concludes,
For many problems no transformations can be found, but the
scope for ingenuity is considerable, and the advantages of
the transformation equally so. (p. 75)
The Transformation of Variables Method
The object here is to transform the independent variables while
leaving the objective function unaltered (Box, 1965 and 1966). The
most important transformations for our purposes are
2
(1) xi t 0 2 x. = y
(2) xi > 0 xi= lil
(3) L X < U = Li + ( L) sin2
(3) Li < x < Ui + xi = Li + (Ui Li) sin yi"
A constrained f(X) problem may be transformed into an unconstrained f(Y)
problem by rewriting the original objective function in terms of the
transformation variables (Y). Transformations (1) and (2) restrict
x(i) to positive values, while (3) constrains x(i) between constant
lower (L) and upper (U) bounds. In addition, one may use (3) to enforce
an equality constraint by setting L U. (L may not be equal to U, but may
be set arbitrarily close). The constraints in the present work in all
cases consist of bounds on the design variables. This is therefore a
class of problems admirably suited to the transformation of variables
method. (Note that while transformation (3) leads to a periodic solution
in yspace, difficulties may easily be avoided provided the search
technique does not take steps so large as to jump from peak to peak;
careful application of the search will insure that this does not occur.)
An important property of these transformations is that although
some are not 1:1 mappings (e.g., the examples given above), they cannot
introduce additional local optima (Box, 1966). Many other transforma
tions are available, and the possibilities for invention are great.
Finally, it is worth noting that the transformation of variables method
does not require function evaluations outside the feasible region.
This is important here because, as discussed previously, the objective
function will be developed based exclusively on information obtained
inside this region. Further, the validity of the experimental model
cannot be guaranteed beyond the physiological limits imposed by the
actual structure. Thus in this application the method of transforma
tions holds a distinct advantage over the aforementioned methods which
rely on (in this case, unreliable) external function evaluations.
This method, combined with one of the powerful unconstrained search
techniques, will provide an effective tool for solving the biomechanical
optimization problem under consideration.
Unconstrained Search Techniques
These methods admit of division into two categories: gradient
methods, which require explicit (or numerical approximations to) function
derivatives with respect to the design variables, and derivativefree
methods. Iterative line searches are required in all of these methods,
so greater efficiency is primarily achieved by more accurate choices in
search direction. It is in the manner of choosing these directions that
the methods find their distinctive personalities.
Gradient Methods
Generally speaking, the explicit derivative information required in
these methods is used to direct the search along lines calculated to yield
the best local approximation to the minimum. Certain gradient methods
such as FletcherReeves (Fletcher and Reeves, 1964), DavidonFletcher
Powell (Fletcher and Powell, 1963), Davidon (Davidon, 1967), and Fletcher
(Fletcher, 1970) have been shown to be somewhat more efficient than the
derivativefree methods when explicit derivatives of the objective
function are available, especially on very complex problems (Fletcher,
1969; Himmelbau, 1972). They can, however, prove quite unreliable when
used with numerical approximations to the derivatives (Sargent and
Sebastian, 1972; Kuester and Mize, 1973), and are far more difficult and
time consuming to implement than the several very efficient derivative
free methods. Thus we will henceforth consider only the latter methods.
DerivativeFree Search Methods
Historically it was found that the simplest concepts, those of
tabulation, random search, or improving each variable in turn (univariate
search) were very inefficient and often unreliable. Improved methods were
soon devised such as the simplex search (Spendley et al., 1962), refined
by Nelder and Mead (1965), and the method of rotating coordinates due to
Rosenbrock (1960). Both methods are heuristic, which is to say they rely
on an ad hoc rather than a strictly theoretical approach to the problem.
A third method, due to Powell (1964), uses theoretical concepts based on
optimizing a quadratic function as a foundation for both choosing a likely
direction and performing the subsequent line search.
One salient feature of the present problem is that, as stated
previously, all constraints will consist of bounds on the design vari
ables. The implication of this is that all constraints will lie parallel
to the coordinate axes. A search method which could easily realign itself
with one or more coordinate axes as the search progressed would hold a
distinct advantage in this type of problem. The Hooke and Jeeves direct
search method (Hooke and Jeeves, 1961) has this property.
Powell (1970), in an excellent survey of the field, described the
Hooke and Jeeves method as an efficient one that is particularly adept at
following curved valleys in the objective function. This method has found
wide acceptance in the engineering community (see, e.g., Kramer and Sandor
(1975); Reinholz (1983)).
In practice, the method of Hooke and Jeeves proved to be very
efficient in the aforementioned class of constrained problems. A series
of evaluations were run using several (unconstrained) test functions
taken from the literature, including Powell's function (Powell, 1964),
Eason and Fenton's (gear inertia) Problem No. 10 (Reklaitis et al., 1983),
and Wood's function (Himmelbau, 1972), which represent a wide range of
nonlinear optimization problems. In addition, bounding constraints were
imposed on the Eason and Fenton and Wood problems to test the constraint
handling aspects of the optimization algorithms. The method of Hooke
and Jeeves consistently demonstrated an efficiency equal to or surpass
ing that of Powell's method. In particular, Hooke and Jeeves' method
proved to be very robust on the constrained problems, quickly converging
to the solution in cases which had brought Powell's method to a complete
standstill. (Both algorithms were coded in Microsoft Fortran 77 and run
on a Zenith Z100 microcomputer; constraints were handled in both cases
using the transformation of variables method described above.)
With the above justifications, the HookeJeeves method was selected
for use here, and will now be described briefly; the original paper
(Hooke and Jeeves, 1961) should be consulted for details regarding its
development and application to other types of numerical problems.
The HookeJeeves Method for Nonlinear Optimization
The "pattern search" method of Hooke and Jeeves is a particularly
elegant, conceptually simple heuristic procedure which derives much of
its flexibility from the fact that it makes no assumptions about the
objective function being evaluated (save, continuity within the feasible
region). It is founded on the observation that points along a narrow
valley are often nearly colinear, and that very often the most efficient
way to find the minimum is to follow such a valley. It therefore tries
to identify the direction of a valley and uses that as a search direction.
Now if a long straight step is taken along a curved narrow valley, the end
of the step will probably fall somewhere on the slope, but this doesn't
matter because it is easy to identify a direction back down to the val
ley's floor. The actual search proceeds via two stages: pattern moves,
in which long steps are taken along the valley, and exploratory moves
which climb back down to the floor of the valley. Successively longer
pattern moves are made when previous moves along the established pattern
meet with success; thus, the search for the minimum accelerates based on
information obtained in previous moves. Using this strategy, the method
is particularly adept at optimizing objective functions having deep,
narrow curved valleys (the "worst case" scenario for an unconstrained
optimization routine).
This method is also well suited to optimizing nonlinear objective
functions with bounding constraints (such as those employed in this
study). This propitious property arises from the method's ability to
"restart" itself in the presence of unexpected developments. Such a
restart is made by establishing an entirely new pattern based solely on
newly acquired information obtained in an exploratory search in the
coordinate directions.
In the present case where the true minimum of a function might well
lie outside the bounds of feasibility, the above procedure allows the
method to take the fastest route toward the optimum, continually accel
erating until a fence is encountered. At this point, the search will
turn and proceed directly along the fence to the constrained optimum.
Once the search nears the optimum, the method converges on the
solution by taking smaller and smaller steps until the step size matches
the requested accuracy for the design variables. (Note that in the
present study we are concerned primarily with obtaining accurate values
for the design variables based on trends in the objective function, not
vice versa. Thus it is a specified accuracy in the design variables that
we seek, and which is provided in the above procedure.)
The first or exploratory stage in the algorithm consists of small
steps away from the initial base point X in each of the coordinate
directions, i.e. (see Figure 3.1),
X' = X0 + ae
where X' = trial point,
Xo = new base point,
a = step size, and
e = unit coordinate directions.
Each direction is tested in turn, saving only those moves which result in
a decrease in the objective function. The sum of the successful explora
tory moves becomes the new base point, and the starting point for the pat
tern move stage. In this stage, a move is tried in the direction defined
by the initial and new base points and covering an equal distance, i.e.,
X '= X + (X XH )
where X' = trial point,
X = new base point, and
XH = home point (previous base point).
If this first pattern move is a failure, an exploratory search is carried
out at the new base point in order to define a new direction. If, how
ever, the first pattern move is successful, the trial point becomes the
base point for the next exploratory stage, while the terminus of the
previous exploration becomes the new home point. Following the next
(FAILURE)
SOLUTION
'"' BOUNDARY
S........... EXPLORATORY MOVES
 PATTERN MOVES
Xe
Figure 3.1. Hooke and Jeeves method for nonlinear optimization.
Functional search proceeds via exploratory moves in the coordinate
directions, followed by accelerating pattern moves toward the optimum.
Note ability to turn abruptly and follow a boundary to the constrained
optimum.
exploratory move, another pattern move is made in the direction defined
by the new base point (terminus of the last exploration) and the home
point (terminus of the penultimate exploration). This sequence of
exploratory moves followed by pattern moves is repeated, each sequence
imparting an acceleration to the pattern move, until a pattern search
fails. It is at this point that the algorithm "resets" itself by aban
doning the established pattern and starting again from the new base
point. Thus sharp changes in direction from one pattern to the next are
easily accomplished. This allows the search, for instance, to proceed
rapidly along any constraints it may encounter en route to the solution,
as illustrated in Figure 3.1.
When the search arrives at a point from which no further improvement
in the objective function can be found in either exploratory or pattern
moves, the exploratory step size is reduced and the sequence resumed.
Final convergence is obtained by repeating this step reductionsearch
process until the step size has been reduced beyond the requested accuracy.
The pattern moves are the primary venue for function minimization,
and it may be seen that these moves have the inherent potential of moving
directly down a valley or along a constraint, rapidly approaching the
minimum. This potential has been borne out in tests on both unconstrained
and constrained functions, as described above. A flow chart for the
HookeJeeves nonlinear unconstrained optimization algorithm as implemented
by the author may be found in Appendix A.
Global Optimization
Throughout the above discussion on constrained optimization an
implicit assumption was made that once a minimum was found, it could be
considered to be the overall solution. This is seldom a good general
assumption, as a nonlinear function may exhibit more than one local min
imum, i.e., it may be multimodal. There currently exist no straight
forward means to deal with this situation (Kunzi et al., 1968; Herson,
1975). The best approach (Kramer and Sandor, 1975) once an optimum has
been found is to generate several starting points randomly distributed
throughout the feasible region of the design space and to compare the
results of each search with previous solutions. If an improved optimum
is found it is then used as the reference for another search sequence.
If no improvement can be found, it may be assumed that the calculated
minimum is the global optimum.
Conclusions: Constrained Nonlinear Optimization Methods
The basic theory and methods of nonlinear optimization were reviewed.
Procedures for obtaining an approximating objective function based on
experimental data were developed, and the inherent limitations of this
objective function were discussed. Based on these limitations and the
desire for an efficient algorithm, the Hooke and Jeeves method for
optimizing an unconstrained nonlinear function was selected, combined
with the transformation of variables method for handling bounding con
straints and the random sequential search method for global optimization.
Listings of the Fortran routines written to implement these algorithms
may be found in Appendix B.
CHAPTER 4
RESULTS
Two parametric studies were conducted to investigate the resection
arthrodesis structure. The first study, the results of which are
described in Part One of this chapter, comprised a series of thirteen
finite element simulations performed to evaluate the effects of graft
(resection) length and leg rotation on the peak stresses that occurred
in the posterior graft. As discussed in Chapter 2, the orientation of
the grafts with respect to the tibia and the foot was assumed to be
fixed. Leg rotation, for the purpose of this study, was therefore
defined as axial rotation of the combined foottibiagraft structure
with respect to the neutrally positioned femur. It was further assumed
that leg rotation would not materially affect the gait and force pat
terns. The model was verified during this series by noting that a
prominent stress peak was found in the model at the location in which
the aforementioned fatigue fractures were known to occur.
An additional thirteen finite element simulations were performed
in the second study to assess the effects of changes in intramedullary
(IM) rod structural stiffness on the peak posterior graft stresses
in long resections with various leg rotations. The results of this
study are presented in Part Two. The above data were then used to
develop the analytical response surface for the formal optimization
which was performed to obtain the optimum values of rod structural
stiffness (geometry) and leg rotation. These results are described
in Part Three. As discussed in Chapter 2, "stress" will henceforth
refer to the von Mises stress at the point in question.
Part One: The Effects of Graft Length and Leg Position
Two distinct modes of bending were seen in the graft at all lengths
and leg positions, as illustrated in Figure 4.1. The compound bending
exhibited in the anteriorposterior (AP) plane (Figure 4.la) is due to
the AP graft locations, which confine the relative displacement of the
tibia (with respect to the femur) primarily to translation. The signif
icance of this mode of bending is increased by external leg rotation
which places a greater proportion of the force vector in this plane
(see Figure 4.2).
Figure 4.3 illustrates the effect of leg rotation on the linear
AP bending moment distribution along the graft. Note that the effects
of this mode of bending are greatest proximally, and that external leg
rotation most strongly influences the more proximal bending moments.
Figure 4.1b illustrates the second bending mode, occurring in the
mediallateral (ML) plane of the structure. In this case the two grafts
and the rod act together as a simple cantilever to which the tibia is
attached. The significance of this mode is increased with internal leg
rotation (see Figure 4.2), as more force is applied in this plane. The
effects of leg rotation on the linear ML bending moment distribution in
the graft are illustrated in Figure 4.4. Here the strongest influence is
seen distally, both in terms of the resulting ML bending moments and of
the effect of leg rotation on those moments.
66
'OiST1 INP= H: I'
2. 0 3:3'
F' nT 1
POST
'TEF= 1
ITEP=12
DI';'PLHi:EI'ENT
HlUTO SCAHLIHG
ZV= I
D I'T= 41 .3
::F=4.62
'F=3.3 1
IIt ID '"A=2.
HLITOI SfliLI[H';
E IST=4 1.
:F =4. L2
,'F=::. 1
r_,,:A=2 C.01
FESECTIO I IHTH :DESIS 3.1
Figure 4.1. Two modes of bending which occurred in the grafts under
physiological loading: (a, left) compound bending in the anteroposterior
plane (ground reaction force component directed from left to right), and
(b, right) simple cantilevertype bending in the mediallateral plane
(applied force directed from right to left).
xt. Rotation < I Int. Rotation
Y oel. p Ione
\ Ground Reaction
Force Vector
(Orthogonal Component)
Figure 4.2. Relationship between orientation of the applied load
and leg position. External leg rotation increases the loading in the
anteroposterior (AP) plane, while internal leg rotation increases
mediolateral (ML) plane loading. Variations in leg position were
modelled by an equivalent (opposite) rotation of the force vector, i.e.,
internal leg rotation was modelled by external rotation of the applied
load and vice versa.
18.0
15.0
PROXIMAL
12.0
9.0
I cm] / 
IM /
6.0 
20 (INTERNAL) ROTATION
DISTAL  00 (NEUTRAL) ROTATION
3.0 / +200 (EXTERNAL) ROTATION
2.0 1.0 0 1.0 2.0 3.0 4.0
BENDING MOMENT
[Nm]
Figure 4.3. AP plane bending moment distributions in 18 cm posterior
graft for three leg positions.
18.0
\ \
15.0 
PROXIMAL
12.0 
POSITION \ \
ALONG \ ,
POSTERIOR
9.0 GRAFT
I [cm]
6.0 \
D IST AL
3.0 
I I I
2.0 1.0 0 1.0
BENDING
[Nr
Figure 4.4. ML plane bending moment
graft for three leg positions.
2.0 3.0 4.0
MOMENT
n]
distributions in 18 cm posterior
We have seen that two distinct modes of bending occur in the graft.
Bending in the ML plane causes a maximum bending moment to occur dis
tally, which may be reduced by external leg rotation. It must be noted,
however, that a reduction of the distal bending moments by external leg
rotation must necessarily be at the expense of increased proximal (AP)
bending moments. Bending in the AP plane produces the largest moments
proximally. These proximal moments may be reduced by internal leg rota
tion, but only at the expense of increased distal bending moments.
This situation might be further clarified by turning now to the
consequent von Mises stresses, concentrating primarily on the stress
peaks which occur in the posterior graft. The moment distributions of
Figures 4.3 and 4.4 indicate that these peaks will occur at the extrem
ities of the graft, i.e., at the grafthost junctions. However, the
callus which develops in these regions as a result of the healing
process effectively presents a zone of increased area moment of inertia
(I) to the bending moments in these regions; this serves to shift the
stress peaks centrally, away from the junctions. Thus we find that the
peak stresses occur at the interior boundary of the graft callus (i.e.,
at the distal end of the proximal callus and vice versa), as illustrated
in Figure 4.5. It will be seen that the proximal stress peak is
invariably the greater; since this is where the fatigue fractures have
been found to occur, the behavior thus predicted by the model may be seen
to correspond well to that observed clinically.
Figures 4.64.8 illustrate the effect of leg rotation on the
proximal and distal stress peaks (as described above) in the graft.
Three resection (and hence, graft) lengths were investigated: 9 cm,
representing the shortest clinically useful resection length
PROXIMAL
POSITION
ALONG
POSTERIOR
GRAFT
DISTAL
P
d
VON MISES STRESS
Figure 4.5. Stress peaks in posterior graft under physiological
loading. (a) General form of von Mises stress distributions in posterior
graft. Note proximal (p) and distal (d) stress peaks located at the
interior borders of callus (see Figure 4.5b).
Proximal von Mises Stress Peak
Ant. Graft Post. Graft
IM Rod
i Distal von Mises Stress Peak
Callus
.171I
Tibia
I I
Figure 4.5continued, (b) Schematic of resection arthrodesis graft
structure, illustrating the locations of proximal and distal stress peaks
shown in Figure 4.5a.
> INTERNAL
I I
_DISTAL
DISTAL
EE A0
 I EXTERNAL
I I i
10
PROXIMALAL
10 20
[DEGREES]
Figure 4.6. Peak stresses in posterior graft vs. leg rotation:
9 cm graft length.
0"
(MPa]
30
20
II i
75
70
65
60
45
'0,.
( INTERNAL
PROXIMALL
.DISTAL
",o
I EXTERNAL 
30
20
10
10 20
DEGREES
[DEGREES]
Figure 4.7. Peak stresses in posterior graft vs. leg rotation:
12 cm graft length.
/
0
[MPa]
L
75
70
65
60
55
50
45
40
35
> INTERNAL
/ PROXIMAL
DISTAL
>1 EXTERNAL
I I J
30
20
10
10 20
(DEGREES]
Figure 4.8. Peak stresses in posterior graft vs. leg rotation:
18 cm graft length.
[MPa]
II
I
i I
(Figure 4.6); 12 cm, a median resection length (Figure 4.7); and 18 cm,
representing the longest viable graft length for this procedure
(Figure 4.8). Two salient features may be noted on these plots. First,
the proximal stress peak dominates for all graft lengths and for all
feasible leg positions. This peak is markedly higher for longer resec
tions, but progressively increased internal rotation is quickly rewarded
by reductions in the proximal stresses.
A second feature, discussed above in terms of the bending moments
but perhaps illustrated more graphically in Figures 4.64.8, is the com
promise one must make when attempting to reduce the stresses in the
graft. Here we see once again how internal leg rotation can be used to
effect an improvement in the proximal stress peak, albeit at the expense
of an increase in the distal peak. However, we may also note that in
finding the specific amount of leg rotation for which the proximal and
distal stress peaks are equal, we will have found also the situation
for which the overall stress distribution in the graft is the most uni
form. Thus by simultaneously minimizing the proximal and distal stress
peaks we will also minimize the probability of fatigue fractures.
The least squares method was used to obtain analytical expressions
describing the peak stressleg rotation relationships illustrated in
Figures 4.64.8 for each graft length. These expressions were solved to
find the optimum leg rotation (i.e., that for which the proximal and
distal von Mises stress peaks were equal). The results are given in
Table 4.1, and show that the optimum leg position for all graft lengths,
in terms of minimizing the stresses in the graft, is an internal rota
tion of approximately 20 degrees.
Table 4.1
Optimum Leg Positions for Three Graft Lengths
Graft Length [cm]
9.0
12.0
18.0
Optimum (Internal)
Leg Rotation [deg.]
21.7
21.2
19.4
Figure 4.9 illustrates the effect of leg rotation on the overall
stress distribution in the 18 cm graft. The proximal peak is substan
tially reduced with internal rotation, paid for by a concomitant
(slight) increase in the distal stresses. The overall resultwith
internal rotation is a more uniform, and therefore more desirable,
stress distribution.
The critical stress in the graft, i.e., that which must not be
exceeded if fatigue fractures are to be avoided, is the fatigue
strength. This may be taken in the present case to be the product
of the yield strength of cortical bone and the appropriate strength
reduction factorss. The latter must be assumed; however, a factor of
0.5 would be reasonable to account for geometry changes at the end of
the callus, surface porosity due to the healing process (including
revascularization), etc. Burstein et al. (1976) have given a value
of
a = 126. MPa
as an appropriate yield strength for cortical bone. Thus we obtain
as a fatigue strength estimate for the graft,
of = 63. MPa.
This approximate value is shown as a gray band in Figure 4.9 (and
subsequent figures), and serves to emphasize the danger of fatigue
fracture which exists proximally with external leg rotation. One can
also clearly see the dramatic overall improvement which is obtained
with internal rotation.
18.0
PROXIMAL
15.0
12.0
POSITION
ALONG
POSTERIOR
GRAFT
[cm]
DISTAL
9.0
6.0
3.0
a
0
a
M
r ~ .
a :
a
.
a
a .
a :a .
0
20 (IN
0
a ;
0a *
a *
o *
S::*
Qo : critical
S'" i
t I ",.
20 40 60 80 100
0
[MPa I
Figure 4.9. Von Mises stress distributions in posterior graft for
two leg positions: 18 cm graft length. Peak stresses occur at interior
borders of callus (see Figure 4.5), and are minimized with internal
leg rotation of 20 degrees. (For critical, see text.)
I
ETERNAL) ROTATION
EXTERNAL) ROTATION
Figure 4.10 is a composite of Figures 4.64.8, and several features
should be noted. By again superimposing the graft fatigue strength
(defined above) one may see how, for longer graft lengths, the proximal
stress peaks are reduced to levels below the critical stress with inter
nal leg rotation. Also note that the stress peak in the 9 cm graft
never reaches the critical stress; this is borne out clinically by the
fact that no fatigue fractures have been reported for grafts of this
length.
An even more propitious consequence of internal leg rotation for
longer graft lengths may be observed in Figure 4.10. At 20 degrees of
internal rotation the maximum stresses in the longer grafts are less
than those in the short (9 cm) graft. This lends further support to the
assertion that fatigue fractures might be avoided in the longer grafts
by placing the leg in internal rotation at approximately 20 degrees.
A further observation to be made in Figure 4.10 lies in the overall
behavior of the curves relating the proximal stress peaks for each
length to leg rotation within the range of feasible leg positions; the
proximal stress peaks are replotted in Figure 4.11. Note that the gen
eral slope of the curve for the 9 cm graft is small: the stresses rise
gradually from internal to external leg rotation. The general slope
of the curves for the longer graft lengths, on the other hand, are much
steeper: the stresses rise sharply from internal to external rotation,
exhibiting greater sensitivity to leg position. Note also that the
curve for the short graft is concave upward, while the curves for the
longer grafts are concave downward. We may therefore observe two forms
of behavior in these curves: shallow, concave upward and steep, concave
a PROXIMAL STRESS PEAK
oDISTAL STRESS PEAK
18 cm
75
70
65
critical
60
_ 0
$> ( INTERNAL
30
20
10
 E EXTERNAL
6
DEGREES]
Figure 4.10. Composite of Figures 4.64.8: posterior graft peak
stresses vs. leg rotation for three resection lengths. (For acritica
see text.) critical
9 cm
0
[MPa]
55
C
I r I r 1
18 cm
75
70
65
critical
60
~~~~0~
\
Line of inflection at graft length = 1 1.7 cm
C INTERNAL l EXTERNAL
L  i  I  I  I  i ____
20
10
10 20
[DEGREE
[DEGREES]
Figure 4.11.
rotation for three
length of 11.7 cm.
Proximal peak stresses in posterior graft vs. leg
graft lengths, and line of inflection at a graft
(For critical' see text.)
9 cm 
0M
[MPa]
30
I I 1 r 1
downward. The transition between these general forms must occur in the
region of a line of inflection on the threedimensional graft stressrod
sizeleg position surface.
The graft length at which the inflection occurs may be found by
setting the second partial derivative (with respect to rotation) of the
stress surface equation equal to zero and solving for length. A third
degree polynomial is the lowest which may exhibit such a line of inflec
tion. Having already established that the graft stressleg rotation
curves follow a quadratic form (Figures 4.64.8), a least squares fit of
the data to a (lengthwise) cubic surface, i.e.,
o= f(, 6)
met with excellent results (R2 = 0.9998). The second derivative of the
resulting expression with respect to rotation was taken, and solved for
the length at which the line of inflection occurs, giving
a= 11.7 cm.
This line is shown in Figure 4.11.
It was stated earlier that a "12 cm length effect" was seen to exist
in the posterior graft, in that grafts longer than this critical length
had a much greater propensity for proximal fatigue fracture than did
shorter grafts. Considering the evidence discussed above, it is possible
that the "12 cm length effect" might be related to the changes in the
form of the graft stress vs. leg position curves on the analytical stress
surface. The transition from a generally flat, concave upward curve for
short grafts to a steep, concave downward relationship for long grafts
takes place in the region of the inflection which was found to occur at
a graft length of about 12 cm. It may be observed in Figure 4.11 that
grafts longer than this critical length reach critical levels of stress
with little or no external rotation, while shorter grafts require a good
deal of external rotation before nearing the critical zone. In addition,
as described above, the stresses in grafts longer than 12 cm are much
more sensitive to leg position than are the stresses in shorter grafts.
This implies that the risk of fatigue fracture in longer grafts increases
sharply with external rotation, as does the advantage to be gained with
internal rotation.
Summary: Part One
The investigation of the effects of graft length and leg rotation
on peak graft stresses led, in addition to a verification of the model,
to:two primary conclusions. First, the peak graft stress, and thus the
risk of fatigue fracture, is minimized in the model for all graft lengths
with an internal leg rotation of approximately 20 degrees. In addition,
a "12 cm length effect" was observed in the form of an inflection on the
cubic graft stress surface at a graft length of 11.7 cm. Stresses in
grafts longer than this critical length were observed to be much more
sensitive to leg position than were the stresses in shorter grafts.
Longer grafts in neutral to external rotation showed proximal peaks
which exceeded a critical stress level. Thus the model predicts a high
probability of proximal fatigue fracture in these cases, which corre
sponds well to the high rate of proximal fatigue fracture observed
clinically in longer resections.
Part Two: The Effect of IM Rod Structural Stiffness
and Leg Position in Long Resections
The effect of IM rod structural stiffness on the proximal graft
stress peak in long resections was investigated by varying the cross
sectional size (and therefore the area moment of inertia, I) of the
implant. In addition to the 12 mm nominal diameter rod described in
Chapter 2 and used in Part One (above), two other sizes (9 mm and 14 mm)
were chosen which represent the range of IM rod sizes available for this
procedure. The range of area moments of inertia thus provided for study
are shown in Table 4.2, and are based on factory specifications (see also
Miller (1977)).
This series was primarily concerned with the effects of rod size on
posterior graft stresses in longer resections; therefore the graft length
was held constant at 18 cm. Figures 4.124.14 illustrate the results for
the three IM rod sizes. The proximal posterior graft von Mises stress
peak data are represented here by triangles. The circles represent the
maximum von Mises stresses in the IM rod, which occurred proximally
near the region of femoral cortical support. As previously described,
the critical graft stress was taken to be one half of the yield strength
of cortical bone. For the IM rod, the endurance limit is the critical
value. A conservative estimate of the endurance limit in bending for
wrought steel is usually taken to be one half of the ultimate strength
(Burr, 1981). Assuming an ultimate strength of 560. MPa for the
constituent material (316 LVM stainless steel), the critical IM rod
stress was therefore estimated to be
Gf,rod = 280. MPa.
f,rod
Table 4.2
Area Moments of Inertia for Three Modelled
IM Rod Sizes
Intramedullary (IM)
Rod Size [mm]
9.0
12.0
14.0
Area Moment of Inertia
I [cm4]
0.03
0.10
0.18
90
65
80
75
GRAFT\ )
[MPa]
65
critical, graft
60
55
50
550
5
500
450
400
ROD(O)
350 [MPa ]
300
C'critical, rod
250
200
150
20 10 0 10 20
6
[DEGREES]
Figure 4.12. Proximal peak stresses in posterior graft and
fluted IM rod vs. leg rotation: 18 cm graft length, 9 mm rod.
(For critical see text.)
0
GRAFT
[ MPa]
65
Cr
critical, graft
60
55
50
550
500
450
400
ROD (o )
ROD v /
350
[ MPa]
300
critical, rod
250
200
150
* < EXTERNAL>
10 0
10 20
6
[DEGREES]
Figure 4.13. Proximal peak stresses in posterior graft and
fluted IM rod vs. leg rotation: 18 cm graft length, 12 mm rod.
(For critical, see text.)
20
RAFT
GRAFT(
[ MPal
65
cr
critical, graft
60
55
50
550
500
450
400
'ROD
350 [ MPa
300
Critical, rod
250
200
( )
I ]
150
10 0
10 20
6
[DEGREES]
Figure 4.14. Proximal peak
fluted IM rod vs. leg rotation:
(For critical see text.)
stresses in posterior graft and
18 cm graft length, 14 mm rod.
20
As before, these critical stresses are shown as bands in the figures;
the graft and rod stresses are plotted as a function of leg rotation.
In Part One it was asserted that the graft stresses in the model
could be minimized by placing the leg in internal rotation at 20
degrees; Figures 4.124.14 graphically illustrate the consequences of
this maneuver. Again we see how internal rotation reduces the peak
stresses in the graft, for all IM rod sizes. Here we can also observe
that, with increased internal rotation, more of the load is carried by
the rod. Thus with internal rotation our concerns regarding structural
failure must shift from the graft to the rod.
Looking first at the results for the 9 mm IM rod (Figure 4.12), we
see that the proximal graft stress peak lies well above the critical
stress for all leg positions. In addition, the rod stress becomes
critical at about 5 degrees of external rotation, increasing with inter
nal rotation. In this case we may only ask which will fail first, the
graft or the rod; failure of one (followed in all likelihood by the
other) is certain. This prediction is supported by clinical experience
at the University of Florida. The only reported failure of a fluted IM
rod in a tibial turnback procedure occurred when a 9 mm rod was used
with internal leg rotation (20 degrees) in a 15 cm resection. Subsequent
to a proximal fatigue fracture of the rod, a revision was performed
replacing the original implant with a 10 mm rod of similar length and
placing the leg in external rotation. No IM rod failure has been
reported for this case 2 1/2 years postrevision.
Turning now to the results for the 12 mm rod (Figure 4.13), we can
see the now familiar graft stress pattern which rapidly improves with
internal leg rotation. Critical graft stress levels occur with neutral
to external leg rotation, becoming more favorable with increasing inter
nal rotation. Also to be observed are the maximum rod stresses which
increase linearly with internal leg rotation. In this case the maximum
stress in the rod falls just short of its critical value at the optimum
(for the graft) leg position of minus 20 degrees.
Finally, looking at the results for the 14 mm rod (Figure 4.14), one
may observe that safe stress levels are predicted almost throughout the
range of viable leg positions in both the rod and the graft, becoming
critical (in the graft) only with 1020 degrees of external rotation.
Summary: Part Two
The consequence of reducing the proximal stress peak in the
posterior graft via internal leg (tibia plus grafts) rotation was shown
to be an increase in the IM rod stresses, due to a shift of the bending
load from the graft to the rod. It was observed that a 9 mm nominal
diameter rod provided insufficient support in the long (18 cm) resec
tion arthrodesis model: graft peak stresses were considerably higher
than the critical level for all leg positions, and the maximum rod
stresses became critical with neutral to internal leg rotation. Thus
the model predicted failure of the graft and/or rod for this case.
The optimum leg position in the model for 12 mm and 14 mm nominal
diameter IM rods was 20 degrees of internal rotation with respect to
the neutral femur. This leg position minimized the stresses in the
posterior graft, a strategy made viable by the greater bending strength
of the larger rods.
Part Three: Optimization of Rod Size
and Leg Position in Long Resections
In this section we will use the algorithms and routines developed
in Chapter 3 to formally optimize the surface represented schematically
in Figure 4.15 (a composite of Figures 4.124.14). It might be suggested
that this is superfluous, since one may see by inspection that the opti
mum values for leg position and rod size, i.e., those values that mini
mize the proximal stress peak in the graft, are minus 20 degrees
(internal rotation) and 14 mm, respectively. In addition, the trends are
clear: for any value of leg rotation the largest rod size is always the
best choice, while for any rod size the maximum feasible amount of inter
nal rotation (20 degrees) should be used.
It is intended, however, that the methods developed here should
be directly applicable to the optimization of more than just these two
variables. As discussed in the final chapter of this work, there are
many other effects which should be investigated such as rodbone
interface stiffness and friction, alternative rod and graft geometries,
etc. Any or all of these might lead to greater understanding of the
operative biomechanical interactions in the structure, as well as to
further improvements in the reliability of this procedure. All of these
effects, however, influence each other as well as the structural behavior
generally, probably in a nonlinear fashion such as demonstrated in the
work already described. Therefore, future parametric studies must not
only take into account the isolated effects of a single design variable,
but rather must add to the information base already obtained regarding
