Computer aided stress analysis of long bones utilizing computed tomography

COMPUTER AIDED STRESS ANALYSIS OF LONG BONES UTILIZING COMPUTED TOMOGRAPHY SHLOMO A. MAROSI and MARTIN J. LINDEN New Jersey Institute of Technology. 323 Dr King Boulevard. Newark. NJ 07102. U.S.A. Abstract --A computer aided analysis method has been developed which utilizes computed tomography (CT) and a finite element (FE) computer program to determine the stress-displacement pattern in a long bone section. The CT data file provides the geometry. the apparent density and the elastic properties for the three-dimensional FE model. A developed pre-processor generates the FE model of a human diaphyseal tibia section which is then analyzed by the SAP IV finite element program. The results obtained are sorted and displayed by a developed post-processor and compared with stresses and deformations from the literature. The model generation method was verified by applying it to a model of simple geometry and boundary conditions. then comparing the results with the analytical solution of the same problem. The convcrgcncc behavior of nodal displacements was tested as a function of mesh refinement. This method provides an automatic. versatile. non-invasive and accurate tool of long bone modeling for finite element stress analysis.

lyzed by these methods were limited in gcncrality and accuracy and sometimes bones have an intricate geometrical structure

liuman

composed linear,

of non-homogeneous,

viscoclastic

materials.

anisotropic,

Numerical

arc often the best way to realistically and deformations

the major dilhcultics clcmcnt analysis gcomrtrical

non-

nificant

in modcling long hones for finite

density

propcrtics

of

the thcsc

have recently documented sig-

associations

used tmnsverse sections

of epoxy

to describe femur.

ci1stx

and an optical

tibia and fibula

McElhancy

c’r (I/.

with co-workers (1976) and alone (1980) bound-

bctwecn the apparent

of bone tissue and its strength

and I-laycs, 1970. 1977; Carter

Scvcral methods have been usod to dctcrminc long

digitizer

positive

(11 ul.. 19X1). A porous

bone geometry for finite elcmcnt models. tlight

portional

(81crl., 1979; Harrigan

block model. considcrcd

c’/ trl. (lY70).

ical propcrtics.

suggests

such as the elastic modulus,

arc pro-

to the apparent density raised to the third

power. Extending

this relationship,

Carter and )laycs

(1976) suggested that compact and canccllous bone can be viewed as a single material, property differences principally

(1984) used X-ray photographs

by

that local mechan-

Brown

DiGiogia

and stilfncss

cr N/., 1970; Carter

arics. Chand er cl/.( 1970). Murase cr ul. (19X2). and also md

of

by the wide range of stress

(Galantc rr al.. I97(k McHlhancy

structures.

(1975). Piziali

the destruction

in the litcraturc.

Many investigators

One of

mechanical

is emphasized

analysis results

analyze strcsscs

has been accurately measuring

and

method

proccdurcs

for such complex structures.

required

the bone sample. The need for an improved modeling

with mechanical

a function of variation

of live long bones for the same purpose. Lovejoy and

in the apparent density. Compression

Hurstcin

and bovine compact and cancellous bone specimens

(1977) used non-invasive

laminography

to

tests of human

obtain geometrical properties of bone sections. Minns

were conducted in these studies. The results yield an

er ul. (1975),

empirical

(1979). Askiew Littlc

Hayes et (11. (1978). Stein and Lewis

(1981).

cr (I/. (1986) used combinations

and Gmnik

Lewis

(1977) and

of physical bone

or casttng sections, X-ray photography

and digitizing

relationship

and strain modulus

between the apparent density

rate of the bone tissue in compression.

meaningful

predictions

techniques to determine the long bone geometry for

radiographic

IWO- or three-dimensional

computed tomography.

tinitc element models.

The bone elastic properties which arc needed for the finite

clcmcnt analysis

were often

mechanical and ultrasonic

evaluated

using

methods for the regions of

Computed

density

and its elastic

relationship

of bone stilfncss

mrasuremcnts

tomography

sively as a diagnostic

This

as obtained

is currently

results of such studies wcrc reviewed by Evans (1973).

et al. (1980) reported on a discrctization

Reilly

anatomical

Buskirk

and

Ashman

1975). Katz

(1981).

Gibson

(1980).

Van

(1985)

and

others. The finite element models generated and ana-

scans for bone modeling. Huang

structures

cf al. (1983) utilized

based on CT CT

procedure of scans. Mayott

scans of a tibia to calculate

arca and mass properties along the bone. Woolson

et

al. (1984) developed stereoscopic CT

to

display a three-dimensional sectional 3Y9 Em I,:,..

used exten-

mechanical research. A few attempts have been made recently to use CT

(1974.

by

tool and for medical and bio-

cortical and trabecular bone and of the cartilage. The and Rurstein

allows

based on

scans.

Nelson

techniques

image of bone using cross cr ul. (1985)

used CT

to

S. A. M~rnou

400

determine the geometry of individual

and bt. J. !_INDEN

femurs and used

the models as input for a computerized

implant design

In this study, a CT data file was used to provide the bone geometry

and elastic

properties

for a three-

dimensional

finite element model. The CT method is

an excellent

approach

for this purpose

for several

-There

is very little

distortion

between

the bone

section and its CT image. images

analog

are stored

to digital

in digital

conversion

evant data are extracted

first was a directory

radiographic

A pre-processor

files: the

was followed

was developed

finite element program

by a

for the SAP IV

with the main task of reading

format

and

no

is needed when rel-

from CT images.

density information

a three-dimensional

finite element model of this bone.

The

consists

pre-processor

-image

-input

generator display file generator

for the FE program.

The CT computer

facilities are generally available used the method

CT method

in hospitals and

is safe and routine. the in-

is feasible for determining

uiro cross sectional geometry

of human

structing graphic VAX/l

terminal

aided analysis method developed for

this research uses an auxiliary

software system. includ-

ing a pre- and post-processor,

for a selcctcd finite

integer

readable

computer

by the VAX/l

wcrc rcconstructcd array calculated dcr: the

generation),

model

numbering

and

display.

result

node

and

intcrprctation.

clcmcnt The

re-

pre-pro-

it on a

environment.

The

was used in this study

and floating

verted utilizing

as gcomctry definition,

(mesh

of recon-

to

the cross sectional images. The CT infor-

element program. The system performs such functions geometry discrctization

is capable

in the ECLIPSE

I mini-computer

reconstruct

bones.

procedure

the scanned image and displaying

mation. The computer

and display

needed for stress

analysis. as currently

main

detector

-mesh

properties of bone. This provides the finite element model with the elastic properties

reconstruction

-edge

(CT numbers

of the following

sections:

-model

-The

in consecutive

file which

number of display fires.

of image pixels) can be related to the mechanical

-CT

on the tape in 4096 word

records organized

the CT data file of a long bone section and generating

reasons:

-The

data were formatted

physical

program.

-CT

tibia section was archived onto a magnetic tape. The CT

pixel

rcconstructcd

point

data,

procedures

were con-

into a format

I. The cross sectional images

lint by line in pixels. using a MAP

with

V;I~IICS

siLc and

area. A VAX

display the rcconstructcd

avaihtblc in the hlc hca-

the s -_r diamctcrs lint printer

of the

was used to

images by a developed com-

ccssor reads the CT data file of a scanned long bone

puter routine. This routine utilized numeric character

section and gcncratcs a three-dimensional

codes to rcprcscnt the radiographic

finite cle-

density range of

mcnt model which is used as input for the sclcctcd

the image pixels, a higher numeral

tinitc clement program.

density range (Fig. I). The CT pixel data comprising the reconstructed

I’ROCEDURE

dimensional ccss data

indicates a higher

images wcrc rcorganizcd

matrix

structure.

to provide The

matrix

in a thrce-

an cflicicnt direct achas the form

NCT

A human male tibia was separated from a cadaver limb and cleaned of external scanned by a General the

University

School,

Hospital

Newark,

New

the

tissues. The tibia

Electric CT/T

produced

by

ECLIPSE

mini-computer.

of

New

Jersey. The

dedicated

was

8800 scanner at Jersey

Medical

scan data

Data

General

were S/200

The tibia was immersed in

water during scanning to reduce the beam hardening artifacts (Brooks

and DiChiro.

scanned simultaneously of dipotassium the density correction (Cann

1976). The bone was

with four standard

hydrogen

phosphate

range of human

cortical

for possible drift of the CT

and Ganant,

0.7 mm diameter

solutions

(KrHPO,)

with

bone, to allow number

scale

1980). A pair of steel wires of

was asscmblcd in parallel along the

bone as a geometric

reference. The tibia was scanned.

slice by slice, in designated spacings of 2.5 mm starting at the proximal

end using the CT automatic

spacing

option. The first CT data file produced contained radiographic

density information

images, each with approximately

64,000 pixels, 0.8 by

0.8 mm. After being processed by the ECLIPSE puter, the information

for the

the

of 44 cross sectional com-

1IO mm long diaphyseal

Fig. I. Numeric character codes rcprescnting radiographic density ranges in a reconstructed bone cross sectional image.

Computer

(K, L. M) where K indicates the cross section number, L the line and M the pixel. A first order edge detection routine was developed to determine the compact bone boundaries in each cross section. This routine searches the image pixels line by line. testing each pixel to determine whether it is in the range of compact bone density or not. The detected boundaries are stored and displayed by the VAX line printer. A mesh generation procedure was developed in the pre-processor to create the three-dimensional finite element model from the two-dimensional consecutive cross sectional images. Commercial mesh generation routines were not selected to discretize the model in this study. Using such routines will result in losing most of the pixel information which is used later to determine the material properties needed for the stress analysis. The mesh generator automatically creates a model of user-defined size and mesh refinement. Since the CT provides consecutive parallel cross sectional images, the mesh generation strategy is first to create a two-dimensional mesh in each individual image and then connect them to form a three-dimensional mesh. A constant nun&r of selected finite element nodes are defined in each cross section. The mesh generator consistsof a node selection routine, an clcmcnt dcfinition routine and an elastic definition routine. The node sclcction routine i~ulo~~:ttic;tlly divides the user sclcctcd number ofcross sectional images into the user-dcfincd grid. A global cartcsinn coordinalc system was selcctcd to dcline the node locations. A Cartesian approach was sclcctcd due to the pixel data organization in the three-dimensional matrix which reflects the pixel arranpmcnt in the images. The node selection algorithm applies a discrete division routine on the d&ctcd compact bone boundaries, to dcfinc

Fig. 2. A three-dimensional

-WI

aided stress analysis of long bones

YGRID number of ‘selected node lines’. Each node line is then discretized by the same method to XGRID number of nodes. The node index N is assignedto each selectednode as a function of its geometric location K. f.. M in the global system and calculated by the following expression: N=(K-l)+C+(Lfor:

I)*XGRID+M

(I)

K = I to KTOT cross sections (Z direction) L = I to YGRID number of lines (Y direction) AI = 1 to XGRID number of pixels (X direction)

where the constant C is selected by C>XGRlD* YGRlD to determine the node index range in each cross section. The selected nodes are filed in a two-dimensional matrix NODE(N. f) where the node index N can have a maximum value of N = KTOT+ YGRID*XGRID. The index I stores the R. L, M geometric information of each filed node. The real geometric location of each node is then calculated using the pixel size and the distance between the scanned cross sections. This method of node numbering allows an automatic and eflicient generation of three-dimensional elements for each two adjacent cross sections. The selected nodes are displayed graphically by the VAX line printer for each cross section, using a dcvcloped character code printing

routine.

An elcmcnt definition routine defining three-dimcnsional tight-node isoparamctric clcments follows the node sclcction. The strategy of automatically forming the elements for each two adjacent cross sections is coordinated with the node numbering method. Sclection of a constant number of nodes in each cross section allows no free nodes after the elements are formed. The elcmcnt dclinition algorithm takes into

FE mesh as generated for the tibia mid-shaft section.

402

S. A. MAR~Mand M. J. LINDEN

account the intrinsic bone geometry and yields a general procedure for modeling a human long bone section. The elastic properties definition routine calculates an average elastic modulus for each element in the model. This procedure consists of the following steps: I. An average radiographic density is calculated for each element based on information for each pixel within the element boundaries. Border pixels are tested as mostly within the element borders, and are thereby included in the averaging process. 2. The relationship between the radiographic density and the apparent density of the bone is determined using the information of the KrHPO* standard solutions, scanned together with the bone. Averaging the pixel information of the standard images yields four data points from which a linear relationship is obtained by applying a least squares method. 3. An empirical relationship between the compressive elastic modulus and the bone apparent density was obtained by Carter and Hayes (1976. 1977). Using this relationship. an average elastic modulus is calculated for each element in the model. This property is valid for a given strain rate. A relationship between the Poisson ratio and the apparent density in bone was not found in the literature. Thcrcforc. a constant value for the Poisson ratio was ass&cd to the cntirc model (Reilly and Burstcin, 1974). The shear modulus is computed in SAP IV for eachelrment. using the given elastic propcrtics. Future investigations of bone property relationships may provide improved elastic properties for the generated model. The three-dimensional model is now formatted to create an input file for the SAP IV finite element program. This program was selected due to availability of the source code and the capability of analyzing large size models. The model is automatically generated with a user-defined mesh size while the boundary conditions, loads and initial displacements, are specitied for each load case analysis. A thrcedimensional display of the analyzed model is obtained by the plotting module of ANSYS finite element

Fig. 3. A typical cross sectional mesh of the tibia.

software package (Fig. 2). The model is also displayed slice by slice on a Tektronix 4107 graphic terminal using PLOT IO graphic software routines (Fig. 3). These displays provide the means to visually examine both the slice series and the entire model. The output generation routine of SAP IV was modified to create an additional output file containing the model information and the computed stressesand displacements of the model elements and nodes. This file was used as input to a post-processor developed to sort the obtained analysis results to find maximum and minimum values of displacements and stresses.In addition, the stress pattern is displayed at any userselected bone slice using PLOT IO graphic routines, giving the option of displaying any combination of slice number and stress direction. The general procedure of model generation was verified by modeling and analyzing a three-dimensional object of simple geometry. A homogeneous isotropic beam with a uniform rectangular cross section was scanned and reconstructed by the developed pre-processor to form a model of I08 isoparametric brick elements. The beam was analyzed as a cantilever with a uniformly distributed load. The model generation procedure was verified by comparing the computer analysis results with the analytical solution of the same model and boundary conditions. The results obtained differed by 1.4% in the maximum deflection and 2.9% in the maximum bending stress. A convergence test was conducted to examine the behavior of nodal displacements as a function of mesh refinement. The model of the scanned diaphyseal tibia was generated, loaded and analyzed with eight ditferent mesh sizes. As occurs in a valid finite element model, the computed displacements of the four cxamined nodal points converged to asymptotic result as the mesh was relined. A demonstration problem was performed using the scanned tibia diaphysis. A section 45 mm long, starting 55 mm from the proximal tibia plateau was modeled, fixed at its distal end and loaded by a uniform compressive load of 3.5 times average body weight (Harrington. 1976) to simulate axial load on the tibia in a simple gait. The results of the computerized analysis were compared with stressanalysis results of a three-dimensional finite element study of the upper tibia, recently obtained by Little and co-workers (1986). In that study the model geometry was determined by casting a cadaver tibia in a polyester resin, slicing it in transverse planes and then photographing and digitizing the slices. Assuming a heterogeneous and isotropic model, the stiffness of the cancellous bone was determined by indentation tests and the stiffness of the compact bone by values from the literature. A three-dimensional finite element model was generated and axially loaded by a total compressive load of 3.5 times average human weight. The resultant stressesand displacements of both analyses were in good agrecmcnt in order of magnitude and in their distribution in the models.

Computer

around

The

computer

aided

analysis

method

yields an

efficient modeling tool for finite element stress analysis and includes the following -The

main features:

processor automatically

scanned by a GE CT/T

generates a three-di-

8800 scanner.

model size and the mesh size are user-deter-

mined. -The

generated model may be viewed

at the image

reconstruction,

model completion -The

and

stages.

stress analysis results are sorted and displayed

by the post-processor -The

and examined

node selection.

at user-selected locations.

processors are coded in FORTRAN

ing debugging information

77. includ-

routines to indicate improper

or ambiguous

modeling

input

and analysis

attempts. The computer reliable

aided

modeling

anrrlysis method

tool,

validated

model and by comparing the literature.

nodal displacements In conclusion, versatile,

a

by the verifying

its output with results from

Convergence

our of the generated

provides

tests validated the bchavi-

FE model

by testing selected

as a function of mesh rclincmcnt.

this method

non-invasive

and

provides

an automatic,

accurate

procedure

for

long bone linitc clcmont modcling. This may serve as a research

tool

towards

biomcchanical

bcttcr

undsrstanding

of the

of load bearing long bones.

bchaviour

The uscfulncss of this method can bc extcndsd by further

studies in the arca of modcling analysis and

applications.

The modeling method can &

cxtcndcd

to apply to the entire long bone. tibia or femur. Using a more advanced finite element software package may improve

the analysis

capabilities

efiiciency,

the model display

and the post-analysis

features of result

interpretation,

An extension

of the computer

analysis

may involve

the use of normal

study

pathologic results

aided and

bone samples for modeling. The analysis

would

provide

the data for determining

the

influence of selected factors on the bone mechanical response. The bone with

study of special cases such as a long

an implant

improved

implant

can provide

design

a method

and implant-bone

for

inter-

action. A long range goal of this study is to extend the analysis method to serve as a research and diagnostic tool for in cico bone modeling. At that stage, a long bone of a specific patient

will be scanned to provide

the finite element

while gait analysis will be

conducted

model.

to provide

the boundary

conditions.

may scrvc as a research and diagnostic standing

Iong

bone

biomechanical

This

tool in underfunction

a prostheses in the proximal

tibia. J. biomcch.

Enyny IO3.239-245. Brooks. R. A. and DiChiro, G. (1976) Beam hardening in X-ray reconstructive tomography. Phys. Med. Bid. Zl(3). 390-398. Brown, T. D. and DiGioia. A. M. (1984) A contact coupled finite element analysis of the natural adult hip. J. Bio-

mensional finite element model of a long bone shaft -The

403

aided stress analysis of long bones

and

remodeling.

mvchunizs 17161 437-448. Cann. C. E. and Genant. H. K. (1980) Precise measurement of vertebral mineral content using computed tomography. J. Compirr. Tomos(. 4(4l. 493-500. Carter, D. R. and Hayes, W. C. (19761 Bone compressive strength: The influence of density and strain rate. Science 194. 1174-t 176. Carter, D. R. and Hayes, W. C. (1977) The compressive behavior or bone as a two-phase porous structure. J. Bone Jt Surq. 59A. 953-962. Carter, D. R.. Schwab. G. H. and Spengler. D. M. (1979) The etlect of apparent density on the tensile and compressive properties of cancellous bone. Truns. Or~hop. Rrs. Sot. 4. 87. Chand. R., Haug, E. and Rim. K. (19761 Stresses in the human knee joint. J. Biomec%onics 9. 417422. Evans. F. G. (IY72) hfedwnictd Propcvties r!/. Bww. Charles Thomas. Springfield. Illinois. Galame. J.. Rostoker. P. S. and Ray, R. D. (1970) Physical properties of trabecular bone. C&f Ti‘sst~t* Ht*s. 5. 236-246. Gibson. L. J. (t9XS) The mechanical behavior of cancellous bone. J. Mio,?trc.hunic.s1%. 3 17-32X. Harrigan.T. P.. Carter, D. R.. Mann, R. W. and Harris, W. H. (19x1) The intlucnce of apparent density and trabccular orientation on the elastic modulus of canccllous bone. liu/~s. Whop. Res. Sot. 6, 277. tlarrington. I. J. (tY76) A bioengineering analysis of force actions at the knee in normal and pathological gait. Wiwncd. Eflgng I I, I67 ‘.172. llaycs. W. C.. Swcnstrn. L. W. Jr and Schurman, D. J. (IY7X) Axisymmetric linitc element analysis of the lateral tibiai plateau. J. Siomrchcrnics I I, 2 t -33. Hight, T. K., Nagcl. D. A. and Piriali. R. L. (IY75) A detailed structural analysis of the human tibia. ASME Appl. hfrch. Dia. IO, Y7-9Y. tluang. tl. K.. Suarez, F. R. and Toridis. T. G. (1980) Utilization ofcomputerized tomographic scans as input lo tinite elcmcnt analysis. Inrrrttufionul Cof~/rrcnce Procwdinqs: Finirc Ekmcnrs in liomuchunics 2. 797-1115. Katz, J. L. (1980) The structure and biomechanicx of bone. Symp. Sot. Exp. Bid. 34. 137-168. Lewis, J. L. (1977) Stress analysis or some features or knee prosthesis by finite elements. Truns. Orthop. Res. SW 2.55 Little. E. G.. Wcvers. H. W.. Sin. D. and Cooke, T. D. V. (1986) The three dimensional finite element analysis of the upper tibia. J. biomuch. Enyny 1011, I I I -I 19. Lovejoy, C. 0. and Burstein. A. H. (1977) Geometrical properties of bone sections determined by baminography and physical section. J. Biomechunics IO, 527-528. Mayott. C. W.. Langrana. N. A., Alexander. H. and Curtis. G. (1983) Geometric and mass properties of bone as measured using computer-aided tomography. Adrunrvs in Bioengine&q (Edited by Langrana. N. A.), pp. 28-31. The American Society of Mechanical Engine-ring. New York. McElhancy. J.. Alem. N. and Roberts V. (1970) A porous block model for canccllous bone. Am. Sot. Mrch. Ettyny ‘JO-WA/BHF-2, l-9. Minns, R. J., Bremble, G. R. and Campbell. J. (1975) The geometrical properties the human tibia. J. Biomechonics 8, 253-255. Murase, K.. Crowningshield. R. D.. Pederson. D. R. and Chang. T. S. (1982) An analysis of tibia1 component design in lotal knee arthoplasty. J. Binmrchunics 16. 13-22. Nelson. P. C.. Robertson. D. D. and Walker J. W. (1985) A computerized femoral intramedullary implant design

or

REFERENCFS Askicw, L. M. and Lewis. J. L. (1981) Analysis of model variables and fixation post length etTects on stresses

401

S. A. MAROMand M. J. LIUE~

package utilizing computed tomography data. IEEE I I I. &xa-201. Piiali. R. L. (1980) Geometric properties of human long bones. 1. Biomtvhanics 13. 8(11-885. Piziali. R. L.. Hight. T. K. and Nage), D. A. (1976) An extended structural analysis of long bones - applications to the human tibia. J. Biomrchunics 9. 695-701. Reilly. D. T. and Burstein. A. H. (1974) The properties OC cortical bone. 1. Bone Jr 10014022. Reilly. D. T. and Burstein. A. H. (1975) The ultimate properties ofcompact bone tissue. J. its 8, 393405.

mechanical Sury. 56% elastic and Biomechun-

Stein, I. D. and Granik, C. (1979) The human tibia: a simplified method of radiographic analysis of its cross section with anthropometric correlations. Ann. biomed. Engny 7, 103-I 16. Van Buskirk. W. C. and Ashman, R. B. (1981) The elastic moduli of bone. AS.\IE Appl. 61rc.h. Div. 45. 131-143. Woolson. S. T., Young. S. W.. White. D. N.. Dev. P. and Wood. S. L. (1984) Three dimensional bone imaging and modelling from computer assisted tomography. Trans. Orrhop. Res. Sec. 9. 59.