Title Page
 Table of Contents
 Modelling of the resection...
 Optimization theory
 Conclusions and further work
 Flow chart for the Hooke and Jeeves...
 Fortran source list: Hooke and...
 Biographical sketch

Title: Biomechanical analysis of a resection arthrodesis
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00086023/00001
 Material Information
Title: Biomechanical analysis of a resection arthrodesis
Physical Description: viii, 126 leaves : ill. ; 28 cm.
Language: English
Creator: Williamson, Norman R., 1952-
Publication Date: 1985
Subject: Arthrodesis   ( lcsh )
Bone-grafting   ( lcsh )
Tibia -- Surgery   ( lcsh )
Mechanical Engineering thesis Ph. D
Dissertations, Academic -- Mechanical Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis (Ph. D.)--University of Florida, 1985.
Bibliography: Bibliography: leaves 120-124.
Statement of Responsibility: by Norman R. Williamson.
General Note: Typescript.
General Note: Vita.
 Record Information
Bibliographic ID: UF00086023
Volume ID: VID00001
Source Institution: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: aleph - 000872593
oclc - 14562424
notis - AEG9856

Table of Contents
    Title Page
        Page i
        Page ii
        Page iii
        Page iv
    Table of Contents
        Page v
        Page vi
        Page vii
        Page viii
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
    Modelling of the resection arthrodesis
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
    Optimization theory
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
    Conclusions and further work
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
    Flow chart for the Hooke and Jeeves optimization algorithm
        Page 112
        Page 113
    Fortran source list: Hooke and Jeeves nonlinear optimization routines for global optimization of a constrained function
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
    Biographical sketch
        Page 125
        Page 126
        Page 127
        Page 128
Full Text









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.



ACKNOWLEDGMENTS . . . . ... ....... iii

ABSTRACT . . . . ... .... . ... vii


1 INTRODUCTION . . . . . ...1

The Tibial Turnback . . . . ... ... 1
Problem Statement . . . . ... .. 5
Methodology . . . . ... . ... 7
Overview of the Present Work . . . .. 9


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




Further Work . . . . . .
Conclusion . . . . . .

REFERENCES . . . . . . .

* 120
. 125


: : :


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



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 cross-sectional geometry of

the structure. This information was used to develop a finite element


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


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 length-leg

rotation-posterior 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.



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 long-term 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.


Knee Level- ---- ----


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 graft-host junctions (A and B) and the
proximal posterior graft-host junction (D) have cortical bone-to-cortical
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.


\ Fracture

Ant. Graft Post. Graft

IM Rod---

S Callus



Figure 1.2. Detail of healed tibiall turnback" resection
arthrodesis structure. Callus forms at the graft-host junctions of
the posterior graft, extending into the graft over a length of 3-4 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 diamond-shaped 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


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 anterior-posterior 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


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



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 bone-implant systems has been in the modelling of

the bone-implant 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 three-dimensional interface element which allowed the interface

regions to be modelled in a realistic manner using physiological stiffness


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.


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


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 bone-implant systems




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 three-dimensional 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 GP-3 sonic digitizer, Science Accessories Corporation) and

stored as three-dimensional coordinates in the laboratory's PDP-11/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

three-dimensional 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

(0R6-3 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 quasi-static 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 vector-fused

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, anterior-posterior (AP), and medial-lateral (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.




2 :.













I \



Figure 2.1-continued. (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 *.

I \t I
/ \ /1 ;: A Y

5 .

i i. I ; ...

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


The type of loading imposed on the leg at 15 per cent of the gait

cycle (varus-flexing 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


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

cross-sectional 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 cortical-cancellous 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 load-bearing cortical thickness in each

section. Possible error resulting from X-ray 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 cross-section. (a) Indicated
(A-B) 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.3-continued. (b) Indicated (A-B) 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 cross-sectional geometry of cortical bone

thus provided good accuracy and precision, within the limits of uncer-

tainty imposed by the cortical-cancellous 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 X-ray 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 three-dimensional 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 bone-cortical 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

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


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 graft-host 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


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 straight-line 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


to keeping the element distribution as regular as possible from section

to section, and a limit of 2-3 degrees was placed on deviations from the

18-degree 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 18-degree 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 cross-sections.

The three-dimensional 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.

Three-dimensional 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



3 ..

F ]i :


12 1:: 8-4
20. -,34

::I 1

D = 1
: :F=4. ?C
'.'F =3-7.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= '.'.





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


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 scan-derived 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-


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 three-dimensional 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 3-4 cm, fairing in grad-

ually over the inner 1-2 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 three-dimensional 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


0.01 cm4


0.05 cm4

Remodelled Graft:

Section Area

0.52 cm2


0.03 cm4


0.13 cm4

I(A-P) and I(M-L) are the area moments of inertia
axes, respectively; W(M-L) is the graft width (ML
is the graft thickness (AP dimension).

about the AP and ML
dimension), and T(A-P)


1.38 cm


0.77 cm


1.40 cm


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







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 Bone-Implant 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 three-dimensional interface elements at locations of

observed bone-rod 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 rod-intramedullary 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 location-dependent 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
TDE C= 1

:.JU= -1

T:F= 3. 2.


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
E MA::= 3 40

YI 14


Figure 2.7--continued. (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


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,



1 1:: 4
TDE := 1
PDE _= 1


-' I- 1
TI I '=

,'F = i
::F=4. 1


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 well-developed

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.



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-


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 three-dimensional


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 four-dimensional hypersurface which could not be

visualized readily: formal optimization methods must be employed in this


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


(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.,





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 two-dimensional 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


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 counter-productive 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


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

n-dimensional 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
(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 y-space, 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 derivative-free

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 Fletcher-Reeves (Fletcher and Reeves, 1964), Davidon-Fletcher-

Powell (Fletcher and Powell, 1963), Davidon (Davidon, 1967), and Fletcher

(Fletcher, 1970) have been shown to be somewhat more efficient than the

derivative-free 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.

Derivative-Free 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 Z-100 microcomputer; constraints were handled in both cases

using the transformation of variables method described above.)

With the above justifications, the Hooke-Jeeves 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 Hooke-Jeeves 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






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

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 reduction-search

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

Hooke-Jeeves 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.


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 foot-tibia-graft 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 anterior-posterior (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

medial-lateral (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.


'OiST1 -INP= H: I'
2. 0 3:3'
F' n-T 1
'TEF= 1

D I'-T= 41 .3
'F=3.3 1

IIt ID '"A=2.


E IST=4 1.
:F =4. L-2
,'F=-::. 1

r_,,:A=2 C.01


Figure 4.1. Two modes of bending which occurred in the grafts under
physiological loading: (a, left) compound bending in the antero-posterior
plane (ground reaction force component directed from left to right), and
(b, right) simple cantilever-type bending in the medial-lateral 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
antero-posterior (AP) plane, while internal leg rotation increases
medio-lateral (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.





I cm] / -
IM /
6.0 -

-2.0 -1.0 0 1.0 2.0 3.0 4.0

Figure 4.3. AP plane bending moment distributions in 18 cm posterior
graft for three leg positions.


\ \
15.0 -

12.0 -
9.0 -GRAFT
I [cm]

6.0 -\

3.0 -

-2.0 -1.0 0 1.0

Figure 4.4. ML plane bending moment
graft for three leg positions.

2.0 3.0 4.0

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 graft-host 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.6-4.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







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




Figure 4.5--continued, (b) Schematic of resection arthrodesis graft
structure, illustrating the locations of proximal and distal stress peaks
shown in Figure 4.5a.



EE -A--0

---- I--- EXTERNAL--
I I i



10 20


Figure 4.6. Peak stresses in posterior graft vs. leg rotation:
9 cm graft length.




II i










--|--I--- EXTERNAL -




10 20


Figure 4.7. Peak stresses in posterior graft vs. leg rotation:
12 cm graft length.
















---->1--- EXTERNAL-




10 20


Figure 4.8. Peak stresses in posterior graft vs. leg rotation:
18 cm graft length.




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.6-4.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 stress-leg rotation relationships illustrated in

Figures 4.6-4.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]




Optimum (Internal)

Leg Rotation [deg.]




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


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.














r -~ .
a :



a .

a :a .

20 (IN

a ;
0a *

a *
o *

Qo : critical
S'" i
t I ",.

20 40 60 80 100

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




Figure 4.10 is a composite of Figures 4.6-4.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


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


18 cm






-_- -0





-- E--- EXTERNAL--

Figure 4.10. Composite of Figures 4.6-4.8: posterior graft peak
stresses vs. leg rotation for three resection lengths. (For acritica
see text.) critical

9 cm



I r I r 1

18 cm







Line of inflection at graft length = 1 1.7 cm

L --- i ---- I ---- I ---- I ---- i ____



10 20


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 -



I I 1 r 1

downward. The transition between these general forms must occur in the

region of a line of inflection on the three-dimensional graft stress-rod

size-leg 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 stress-leg rotation

curves follow a quadratic form (Figures 4.6-4.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.12-4.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.

Table 4.2

Area Moments of Inertia for Three Modelled
IM Rod Sizes

Intramedullary (IM)

Rod Size [mm]




Area Moment of Inertia

I [cm4]










critical, graft









350 [MPa ]

C'critical, rod




-20 -10 0 10 20

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

[ MPa]

critical, graft








ROD (o )
ROD v /


[ MPa]

critical, rod




--* |<--- EXTERNAL-->

-10 0

10 20

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



[ MPal

critical, graft








350 [ MPa

Critical, rod



( )

I ]


-10 0

10 20


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.


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.12-4.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 post-revision.

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 10-20 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.12-4.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 rod-bone

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

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs