Experimental determination of whole long bone sectional properties

J. Llionwhwics Vol. IS. No. 4. pp. 297-303. Printed in Great Britatn.

0021-92901821040297-07 SO3.WO 0 1982 Pergamon Press Ltd.

1982.

EXPERIMENTAL DETERMINATION OF WHOLE LONG BONE SECTIONAL PROPERTIES*? ALICE A. GIE~ and DENNIS R. CARTER Harvard/MIT Division of Health Sciences and Technology, and Orthopaedic Research Labs, Massachusetts General Hospital, Boston, MA 02114, U.S.A. Abstract-A strain gage based experimental method is presented which provides an alternative to use of area section41 properties for obtaining the cross sectional antroids, flexural rigidities and axial stiffness of whole bone spkcimens. While area sectional property computations require detailed records on the geometry of cross se$tions, the experimental method presented here requires only gage locations plus the applied loads and recbrded strains. Both techniques are used to analyze the midshaft sectional properties of canine radii. The resjllts demonstrate the advantage of the experimental method for the analysis of heterogeneous cross sections with an unknown distribution of bone elastic modulus. By applying the experimental and area methods together, effective elastic modulus values for the sections were obtained.

properties in structural analyses of long bones. Such an application, however, will not account for cross section heterogeneity. Relaxation of the homogeneity assumption in area sectional property computation can be accomplished by multiplying each area segment by a weighting factor related to its elastic modulus. The density ratio relative to compact bone was used for this purpose by Piziali et al. (1980) on segments containing significant proportions of cancellous bone. The approach taken here, in contast, is to extract sectional properties directly from the forms taken by beam theory equations when the elastic modulus is explicitly assumed to vary over a cross section. The spatial coordinates of three or more strain gages and the strain response of those gages to known bending and eccentric axial loads are the required inputs to the analysis. The results of the analysis may be taken by themselves or combined with area analysis to yield further insights into the elastic modulus distribution in the cross section.

INTRODUCTION

Knowkdge of structural properties of bone cross sections is fundamental to attempts to analyze the mechanical reslponse of long bones. Improved techniques to determine the details of geometric and material prope?ty distributions for whole bone cross sections will be of value in research on the mechanical behavior of musculoskeletal systems. The method to be presented here provides a convenient way to quantify many pf the structural properties of long bone cross sections by simple experiment with strain gages in vitro. Studies of bone cross sectional properties to date have focused primarily on models of the bone as a homogeneous bterial. Sectional properties have then been determine@ on the basis of the distribution of area represented by bone geometry (Nagurka and Hayes, 1980). Values fbr elastic moduli, when required, have generally been assumed. This simplification has allowed much progress on characterization of bone systems in vitro and identification of forces and moments in viva. Piziali et al. (1980) for example, used area methods ain matched human leg bones to obtain data which cohld be used as input in a beam-type finiteekment nlodel. Rybicki et al. (1977)and Carter et al. (1981) ma& use of area sectional properties in analyzing in vivo strain gage information. Huiskes et al. (1981) used arta sectional properties in theoretical stress analyses of a femoral shaft to compare in vim strain gage data to theoretical predictions. Carter et al. (1981) described a strain gage based experimental technique which can be used to establish an effective elastic modulus for a whole bone cross section. The tecihnique provides a reasonable estimate for elastic modlhlus to be combined with area sectional

ANALYSIS APPROACH

l Received 23 June 1981. t Supported by NIH Grant AM 27117 and NASA Grant NAG 2-6.

The overall strategy of the analysis involves the expression of quantities of interest in terms of simultaneous systems of linear equations. The equations used emerge from the theory of straight beams of arbitrary cross section. Since the geometric curvature of whole long bone specimens is small, theory for a straight beam can appropriately be applied. In practice, more equations than unknowns are generally available so the systems of interest are overdetermined. Application of the least squares criterion provides estimated solutions of the systems. The excess of data over the requirements for determinancy is advantageous in opening the door for statistical inference. Confidence levels associated with the solution estimates rise as more data is collected. The cross sectional plane under investigation is

297 BM15:.- E

ALICE A. GIESand DENNISR.

298

CARTER

(0,b) 1:x = j 4x, Y) y2 dA I,, = JE(x,y)x’dA 1:; = j 15(x,y) xy dA A* = jE(x,y)dA e+ I’X’X I*. Y Y’ Fig. 1. Schematic representation of bone cross section to be analyzed. established by placing three or more strain gages around the periphery of the bone, all at the same distance along a proximal-distal bone axis. An example .of such a cross section for the case when four gages are used is shown in Fig. 1. An arbitrary global coordinate system with axes’ Cand r) has been imposed to specify the gage locations. Distributing the gages around the bone periphery as shown rather than grouping them close together enhances the accuracy of the results. Also shown in Fig. 1 is an effective centroid (a, b) whose coordinates are calculated in the analysis. Effective centroidal axes x and y parallel [ and q. The effective centroid possesses the desired symmetry prop erties for a nonhomogeneous cross section, that is F(x,y)ydA

E(x, y) x dA = 0,

=

effective centroid flexural rigidity about x fective centroidal axis flexural rigidity about y fective centroidal axis generalized product of inertia axial stiffness principal angle flexural rigidity about principal axis flexural rigidity about principal axis

efef-

x’ y’

The analysis involves the processing of recorded strain data from bending and axial loading Obtaining the imposed bending curvatures experienced by the beam model from the gage positions and strains is the key to the bending portion of the analysis. Computation based on these curvatures then provides the centroid and flexural rigidities. Processing of the axial loadings involves a procedure of isolating the pure axial component of strain and bending curvatures due to load eccentricity. The axial stiffness is then computed from the pure axial strain component. Further processing of the sectional properties calculated in bending provides the principal angle and flexural rigidities with respect to the principal axes. 1. Bending curvatures Let the bending curvatures ~~ and K= of a beam subjected to biaxial bending be defined by the beam theory equation relating the (x, y) position and strain at the position,

I

J

E = (Ky)x + (K,)Y,

in contrast to the area centroid where xdA = 0.

ydA = s

I

Applying a series of bending loads and then a series of axial loads of unknown eccentricity and recording the associated strains at each gage establishes values for the following variables. number of gages N, number of axial loadings N AL number of bending loadings N axial loads !‘;I;. = 1, 2, 3,..., N,, (M&, j = 1,2, 3, . . .. N,, bending moment component about x axis for jth bending load (M,),, j = 1.2, 3, . . ., N, bending moment component about y axis for jth bending load Strain recorded at the ith &lj gage under the jth loading The values of these variabks together with those of intermediate variables in the computation (to be defined as they are used) will be shown to determine the following sectional properties.

(1)

where K,,and K* are the bending curvatures about the y and x axes respectively. The coordinates of the strain gages (Fig 1) are, xi = Ci - a

(24

Yi = it - b

WI

where a and b are the coordinates of the effective centroid and subscript i denotes the strain gage under consideration. Comparing the other gages to one gage, say the first, equations (1) and (2) may be combined as rl - cI = (K&0(, - cl)+

(K&i

- ttl).i=

1,23,...,N, (3)

where the subscript i denotes the gage number. The unknowns are K,. and K,. For N, r 3 (i.e. three or more gages) the system can be solved for the least squares best estimates of K~ and K= for any applied bending load. Specifically, solving such systems for each ap plied bending load j = 1, 2, 3, . .., N,, yields the curvatures (K~), and (K=)j experienced under the jth loading : Gj -

Elj=

(Ky)j(Ci-CI)

+

(Kz)i('ti-$~)

(4)

Fig. 2. Canine radius with strain gages on the midshaft. Cantikver bending moment is applied by a hanging weight at kft. (From Carter et a/.. 1981).

299

Fig 3. Midshaft canine radius sections for which both experimental and area cross sectional properties we: obtained. (a) Compact specimen. (b) Relatively porous, highly remodeled specimen. White axes are princip. axes (x’, y’) determined from experimental method, the white star is the area centroid, and arrows indica strain gage location. Top and kft directions are cranial and medial aspects, respectively.

301

Determination of whok Ion~gbone sectional properties i = 1,2,3,.. ., N,

j=1,2,3

,..., N,,.

It is these curvatures (K?), and (K,)~ that are required in the computation of the centroid and flexural rigidities.

Since the strain h at each gage and the coordinates (x,, yi) of each gage are known, the system may be solved (using appropriate methods for overdeterminancy, if present) for the unknowns &,,, K,. and K,. Equilibrium then requires that the axial stiffness is A* = PJc,

2. Centroid With the values of the curvature coefficients (K,.), and known, the general beam bending equation (1) may be rewritten for a single gage, again say the first, in the following form : (K&

Clj=tKy)jXl

+(K,)jYl,i=1,2,3,...,~~~

(5)

For N,, 2 2 the system may be solved for The system is overdetermined for N,, > unknown x1 and y, are estimated with accuracy as NgL + X. The coordinLtes of the effective centroid other gages are then obtained as: a=(,

x1 and y,. 2 and the increasing and of the

-x,

(64

b = VI - Yl

3. Flexural

5. Principal angle The principal angle is given by

e* = 0.5 tan - ’ [21&/(1,:. - I&)].

i = 2,3,....N,

yi=qi-b

i=2,3

,...,

Va) N,.

(7b)

This may be verified by considering Mohr’s circle for the generalized moments and product of inertia fzX,I;,, and I$. Alternatively it may be obtained by differentiating the expressions for principal flexural rigidities below. 6. Flexural

rigidities with respect to principal axes

The principal flexural rigidities are

I*X’X’= I:;, cosz 0* + I& sin2 0* - I&, sin2 0*

rigidities

(13a)

The flexural egidities I:, and Z&and the generalized product of inertia I& are determined from the relevant beam theory equations using the curvatures determined in Section 1 above. For a beam in biaxial bending with bending moment components M, and M, M, =

KxIzx

M, =

-K&

+

K&

(W

-

(8b)

K&

Again considering the N BLbending loadings j = 1,2, 3,..., N,, the following system may be written. (M.z)j = (K.Ji 1.L + (KY),I$

(9a)

(My), = - (K,)j I.& - (Ky),I&

(9b)

1.2.3 ,..., N,,.

I*. YYs= f,:.cos2 8* + I:, sin’@ + I:,. sin2 6*.

7. Eflective elastic modulus values When area sectional properties are available, effective ellstic modulus values may be found by dividing the experimentally determined sectional properties by the corresponding area sectional properties :

EL = We,, E,:, = Ip*ylI,, E:y = W,,

E:.,. = I:y/I,~,.

4. Axial stijhess

APPLICATIONS

The axial stiq‘ness can be obtained by analysis of the

results of axial iloading by a procedure similar to that which yielded the curvatures in subsection 1. Experimental applicaticin of a nominally axial load to a whok long bone gen+ally results in application of a somewhat eccentric axial force. This produces biaxial bending e5ect$ plus pure axial deformation. The equation relating strain and (x, y) position becomes +

KVXi+Kxy,,

i= 1,2,3 ,..., N,.

(13b)

These expressions can be derived from first principles by using the integral definitions for flexural rigidities together with rotation of axes formulae for x’ and y’ in terms of x and y.

Using the applied bending loads (IM*)~and (M,), and the curvatures (KJ, and (K,)~, the system may be solved for the three upknowns I:=, I&, and I:,,. Again least squares techniques provide means for dealing with the overdeterminaticy of the system.

Ci=Cg

(12)

(6b)

Xi = (Ii -a

j=

(11)

where P is the applied load. If N,, axial loadings are used, N,,. estimates of the axial stiffness are obtained. Comparison of thesevalues provides an indication of their accuracy level.

(10)

E+ = E+/A

E;.y, = I;y~lIyy

The analytical technique described was applied to ten canine radius’cross sections. To demonstrate its usefulness, the result of its application to two cross sections will be presented. The whok bones were prepared by bonding three or more strain gages at the mid-shaft as reported by Caler er al. (1981). One end of each specimen was potted in a square prism of polyester resin to assure reproducibility of load orientations with respect to specimen

ALICE A. GIESand DENNISR. CARTER

302

and ground (Fig. 2). Bending loads were applied by hanging weights from the free end of the cantilevered bone. By clamping the plastic prism on its various faces and on its comers in a prepared jig, eight discrete orientations of bending for a given applied weight were obtained for each specimen. Axial loading was accomplished by hanging the bone vertically from the prism and suspending weights from the free end. The strain at each gage was recorded for all loading conditions (Garter et al.. 1981). At the completion of the testing, the entire diaphysis of each bone was embedded in an extension of the prism. This provided a fixed frame of reference to which geometrical measurements could be related. Thin sections from the strain gaged diaphysis region were mounted as slides and photographically enlarged. Global coordinate systems were imposed and the gage locations measured. From the global gage coordinates, the applied loads, and the strains recorded under each loading, the sectional properties were computed as described above. Area sectional properties for each cross section were also computed using the method described in Nagurka and Hayes (1980).A PDP 1l/O3 microcomputer interfaced with a digitizing tablet sufficed for all computations. The effectivecentroid, flexural rigidities and generalized product of inertia, axial stiffness,principal angle, and flexural rigidities with respect to principal axes were obtained by the experimental method. The area method provided the area centroid, area moments of inertia and product of inerti& area principal angle, and principal moments of inertia. Differences of location between the effective centroids and area centroids reflected the departure of specimens from cross sectional homogeneity. Division of each experimentally determined sectional property by the corresponding area sectional property (e.g. IV, = E&) provided information on the spatial directional dependence of the effectiveelastic modulus. RESULTS The results of the analysis on ten canine radii were reasonable based upon literature elastic modulus values and histological examination of the specimens. The general shape of the radius sections was elliptical with major and minor axes of approximately 18 and 12mm respectively. The effective centroids of all sections were within 1 mm of the area centroid and were consistently in a more cranial (positive y) dtrection. The effective elastic modulus for bending in the ~ranial-caudal direction (about the x-axis) was consistently lower than that for bending in the lateralmedial direction (about the y-axis). The effective axial elastic modulus for the sections was in the range of 15 to 27 GPa. Two representative cross sections are shown in Fig. .1. Figure 3(a) shows a specimen consisting of compact bone with low porosity. Figure 3(b) shows a specimen

Table 1. Exparimmtal and area sectional properties Section A (0, b) (cm) Area antroid A* (MN) A (cm’) EL (GPa) IL (Nm’) 1, (cm’) E:, (GPa) I$ (Nm’) 1, (cm’) E& (GPa) If (Nm’) I,, (cm4) E& (GPa) B* (degrees) tJ (degrees) I:.,, (Nm’) I,.,. (cm4) E:.,. (GPa) I;,. (Nml)

I,.,. (cm4) E:.,. (GPa)

Section B

0 (cm)

(-L-)62, -0.014) 2.557 1.006 25.4 28.89 0.1137 25.4 56.79 0.2100 27.0 - 14.35 -0.0434 33.1 22.9 21.0 22.83 0.0970 23.5

62.85 0.2267 27.7

0

(-ZO18. - 0.043) 1.502 0.866 17.4 13.90 0.883 15.7 29.61 0.1742 17.0 - 1.346 0.0012 112.17 4.9 -5.1 13.78 0.0883 15.6

29.72 0.1742 17.1

with extensive remodeling and relatively high porosity. The computed sectional properties for these two cross sections are presented in Table 1. With the exception of E!& the effectiveelastic modulus values for the section shown in Fig. 3(a) were between 25.4 and 27.7GPa. Similar effective elastic modulus values for the section in Fig. 3(b) were between 15.6 and 17.4GPa. The unreasonably high values found for E!&are due to the smallness of l$ and I,, relative to I&. I&, and I,,, I respectively. Small errors introduced in the experimdf tal and area methods used will thus cause large errors in the calculation of E&, The principal angle of about 22 degrees for the section in Fig. 3(a) is a result of the rotated orientation of the specimen in its plastic enclosure visible in Fig. 3(a). In contrast, the section shown in Fig. 3(b) with approximately cranial uppermost orientation produced small principal angles.

DISCUSSION It has been shown that sectional properties can be determined by experimental means. Experimentally determined sectional properties can be compared with each other for a given specimen or used to compare one section to another. When both experimental and area sectional properties have been computed, estimates of the effective elastic modulus under various loading conditions can be obtained. The experimental method has been presented here without error estimates. The theory of these is a topic in itself. We are currently working on a method to provide error information through nonlinear curve fitting The agreement between sectional properties obtained by the experimental method and by area

Determination

303

of whole long bone sectional properties

methods substantiates the validity of the experimental method. No homogeneity assumptions are built into the experimenta! method and one need not assume a value for elastic modulus. Sectional properties obtained by the experimental method should, therefore, be more revealing than property estimates based on cross section geometry. We have consistently found that the area centroid of the canine radius midshaft is caudal to the effective area centroid. In addition, the effective moduli E& E!& are less than E;,. and E& respectively. These findings are eonsistent with the histological appearance of the cross sections. In the caudal area of the sections there is invariably more extensive bone remodeling and often greater intracortical porosity. Numerous researchers have found that more extensive bone remodeling is accompanied by a decrease in elastic modulus (Carter and Spengler, 1978). This view is substantiated by the results of our mechanical analyses. The effective axial elastic modulus values (E*) of 25.4 GPa and 17.4 GPa for the two sections presented are in line with bone elastic modulus values in the literature. The difference of these values is reflected in the greater porosity of the section in Fig. 3(b) as seen histologically. Mcrduius values of the other radius sections that we have analyzed also tend to correlate well with histological appearance. Insufficient information has been gathered, however, to indicate whether remodeling and porosity effects are the dominant factor in determining in situ bone modulus or whether aspects of bone chemical composition are of overriding importance. This and a continuum description of the distribution of elastic modulus remain goals for future research. The experimental method of sectional property determination requires application of gages and acquisition of strain data under loading but relatively little geometric information. In contrast, area methods require extensive documentation of bone geometry and, where inhomogeneity is important, histology. Thus the experimental method offers advantages in studies where departures from homogeneity are important. Further refinements will be required, however, before it will be possible to determine using the

information as the spatial distribution dulus for whole bone cross sections.

of elastic mo-

Acknowledgemenrs-We thank William E. Caler and Eva Shimaoka for technical assistance.

REFERENCES

Caler, W. E., Carter, D. R. and Harris, W. H. (1981) Techniques for implementing an in uioo bone strain gage system. J. Biomechanics 14, 503-507. Carter, D. R. and Spengkr. D. M. (1978) Mechanical properties and chemical composition of cortical bone. C/in. &r/lop. 135, 192-217. Carter, D. R., Vasu, R., Spengler, D. M. and Dueland, R. T. (1981) Stress fields in the unplated and plated canine femur calculated from in oiuo strain measurements. J. Biomech. 14.63-70.

Carter,D. R.,Cakr, W. E. and Harris, W. H. (I981)Resultant loads and elastic modulus calibration of long bone cross sections. J. Biomechak 14, 739-745. Huiskes, R., Janssen, J. D. and Slooff, T. J. (1981) A detailed comparison ofexperimental and theoretical stress-analyses of a human femur. Mechanical Properries 01 Bone (Edited by Cowin, S. E.), ASME publication AMD-Vol. 45, pp. 211-234. American Society of Mechanical Engineering, New York. Nagurka, M. L. and Hayes, W. C. (1980) An interactive graphics package for calculating cross-sectional properties of complex shapes. J. Biomechanics 13,59-64. Pixiali, R. L, Hight. T. K. and Nagel, D. A. (1980) Geometric properties of human leg bones. J. Biomech. 13.881-885. Rybicki, E. F., Mills, E. J., Turner, A. S. and Simonen, F. A. (1977) In uiuo and analytical studies of forces and moments in equine long bona. J. Biomechanics 10, 701-705.

NOMENCLATURE (C, q)

(a, b) ‘A”,y) E

I M P

c ;

global coordinates (measured) location of atroid in global coordinates centroid axes coordinates (computed) area elastic modulus moment of inertia bending moment axial load (force) strain curvature principal angle

Subscripts i in&x of gage number

experimental method, all the detail which could be found currently by a painstaking area study incorporating modulus weighting factors for tiny area

i xv Y

in&x for loading number axis about which there is bending or rotation

elements. Perhaps ingenious combination of these methods will be required for the study of such elusive

*

denotes principal axes generalized property defined by R* where R is any sectional property A, I,,, I,, I,, I,.,., I,.,

Superscripts