- Email: [email protected]
A structural model for the mechanical behavior of trabecular bone
1. Biomechunicr.
1973. Vol. 6, pp. 657470.
Pergamon Press.
Printed m Great Britam
A STRUCTURAL MODEL FOR THE MECHANICAL BEHAVIOR OF TRABECULAR Department of Orthopaedics,
BONE*
J. W. PUGH Hospital forJoint Diseases and Medical Center, New York, U.S.A.
R. M. ROSE Department of Metallurgy and Materials Science, Massachusetts Cambridge, Massachusetts, U.S.A.
Institute of Technology.
and E. L. RADIN Department of Orthopedic Surgery, Harvard Medical School. Childrens Hospital Medical Center. Boston, Massachusetts. U.S.A. Abstract- A quantitative model is developed for trabecuiar bone by approximating the trabecular geometry with a hypothetical network of compact bone. For the region immediately beneath the articular cartilage in the distal end of the femur, finite element analyses were performed with a high speed computer, assuming a physiological static load. The results indicate that bending and buckling of trabeculae are considerable in any elastic deformation of the bone: that fatigue fracture in some fraction of suitably oriented trabeculae is inevitable in normal ambulation; and that the stiffness varies considerably with lateral position across the subchondral plate. The latter depends totally on trabecular arrangement and may play a role in joint function and degeneration. The adjustments necessary to bring the gross stiffness into agreement with experiment imply that the intertrabecular soft tissues are of no consequence to the mechanical properties and that the compact bone of which trabeculae are made is probably not as stiff as cortical bone.
INTRODUCTION
AND BASIC ASSUMPTIONS
(a) Organization A WORKABLE semi-empirical model based on porous block geometry has been developed to explain the variation of the elastic (Young’s) modulus and mechanical strength of cancellous bone as a function of density (McElhaney, Alem and Roberts, 1970). Despite the crudity of the geometrical assumption, the model worked well when applied to cancellous bone from the cranium and vertebral bodies. Another assumption used in the above model, one which will be used as well in this paper, is that the properties of the trabecular material are identical to those of compact (cortical) bone. Such an assumption. that mature hard tissue has uniform microscopic properties, seems reasonable but is in fact unproven. *Rc,ccGvrd 29 Novrmber
1972.
(b) Mechanical properties of bone The accepted value for Young’s modulus of compact bone is 2 x lo6 psi and its strength properties are approximately 20,000 psi in compression and 18,000 psi in tension (Yamada, 1970; Currey, 1970; Evans, 1961a; Evans, 1957), the yield strength and the fracture strength being essentially the same. Poisson’s ratio for cortical bone is very close to 0.3 (McLeish and Habboobi, 1972). Trabecular bone has a stiffness of about 10” psi (McElhaney, Alem and Roberts. 1970) and a compressive strength of 500-1000 psi depending on where in the femur the bone is taken (Evans, 1961 b). Experiments have shown the mechanical properties of trabecular bone to be virtuahy unrelated to the fluid and fatty tissue present in the intertrabecular
658
J. W. PUGH, R. M. ROSE and E. L. RADIN
spaces (Swanson and Freeman, 1966). It is therefore reasonable to neglect, in a quantitative mechanical model, the presence of the intertrabecular soft tissues. Cartilage
CC) Structure
and mechanical function of trabecular bone The general structure of trabecular bone consists of a system of sheets of compact hard tissue, with a distribution of widths, the lower limit corresponding to rods. Typical photomicrographs that one observes of trabecular bone (Elias, 1971) can correspond stereologically only to such a distribution of interconnected sheets. Rodlike or beamlike trabeculae would appear very differently in cross-section. The sheet- and rod morphology has been confirmed for vertebral bodies by scanning electron microscopy (Dyson, Jackson and Whitehouse, 1970). The SEM also reveals that the open spaces present in the trabecular bone are interconnected. Considering condylar structures, it would seem that the function of the trabecular bone is to support and direct the operational loads through the ends of the long bones to the diaphyses. Evidence for the central role of the trabeculae in this function can be found in photoelastic studies which show that they are aligned in the directions of the normal stress trajectories (Kummer, 1966). As an example, in the human knee (Fig. l), the femoral condyles rest on the tibia1 plateaus and serve to transmit the body weight to the lower leg. The condyles are roughly cylindrical structures, and consist of a thick layer of cartilage, a thin layer of cortical bone referred to as the subchondral plate, and the region of trabecular bone immediately underneath the plate. The amount of trabecular bone present is large compared to the amount of cartilage and the thickness of the subchondral plate. The region between the cartilage of the condyles and plateaus is filled with synovial fluid, which serves no mechanical function other than lubrication (Radin and Paul, 1969). A mathematical model for the elastic proper-
-
Synwial fluid -Cartilage
I
Joint space
1
Fig. 1. Schematic of human knee joint. The cartilage is very thick compared to a trabecular thickness. The subchondral plate is roughly the same thickness as a trabecula.
ties of condylar structures (Burstein, Shaffer and Frankel, 1970) has shown that, at least under static conditions, the trabecular bone in the condyle is important in load-carrying; 80 per cent of the surface loading is transmitted through the trabecular region. The subchondral plate, furthermore, was found to be loaded with a combination of bending and membrane stresses. Finally, from the difference in properties between cortical and trabecular bone, it was concluded that, as a percentage of the mechanical strength, the stresses in the trabecular region are large compared to those in the cortical bone of the shaft of the femur. The plate-like trabecular geometry means that bending and buckling of individual elements can not be neglected even for small stresses, and therefore we must include bending and buckling in our computations. (d) Technique of computation The use of the computer to solve biomechanical problems is not new. Probably the most recent use is a program for the analysis of stresses in long bones (Piotrowski and Wilcox, 1971). This allows calculation of stresses and strains through the use of standard stress-strain relations and the Airy stress function for a bone consisting of any cross-
MECHANICAL
BEHAVIOR
sectional shape and variation of cross-section with length. An alternative which offers greater possibilities, as long as a high speed computer is available, is the use of finite-element techniques (Zienkiewicz, 197 1) which incidentally has led to a greater understanding of the mechanical behavior of composites (Hackett. 1971; Chen, 197 1). It is also applicable to a number of other systems and lends itself especially well to the analysis of very complex structures. The structural system is considered as an assemblage of idealized elastic elements joined together at discrete nodes. Using the direct stiffness concept, which allows adding together at each node the stiffness coefficients of adjacent elements, a stiffness matrix for the structural system is obtained. This matrix relates the external forces acting on the nodes to the displacements of the nodes: IF) =
[Sl{Sl.
(1)
An inversion of the stiffness matrix [S] allows calculation of the influence matrix, which gives the nodal displacements (6) as functions of the external loads acting on the structural system: IS> = [S]-‘(F).
(2)
This is solved for the nodal displacement, {S}, and once the nodal displacements are known, the solution is specified. Using the same strain pattern that we assumed in derivation of the stiffness matrix, one can derive the matrix of stress coefficients giving the stress in the nth element as a function of the nodal displacements: {a?,} = r ~?Il~hJ.
(3)
A continuum can be represented as an assemblage of elastic elements joined at nodes as above, and the equations for internal displacements and stresses can be derived using the matrix techniques outlined. The finite-element
OF TRABECULAR
BONE
659
method is in fact a solution that minimizes the total potential energy of the system (Zienkiewicz, 197 1). The system used in the present study is ICES STRUDL II (Logcher et al., 1968, 1971) developed at the M.I.T. Civil Engineering Systems Laboratory. This system allows the calculation of a stiffness analysis on a structure that is specified geometrically in conventional three-space. The stiffness analysis is a linear, elastic, static, small displacement analysis applying the finite-element method in an efficient algorithm. It requires the specification of element properties in acceptable form and treats the nodal joint displacements as unknowns. The ‘PBST’ plate bending and stretching element used in this analysis (McDonnell) is a six degree of freedom (six nodal unknowns), flat, triangular element which considers both the in-plane (membrane) and plate-bending forces. The element formulation uses the following displacement expansions: two for orthogonal stretching in-plane, two for inplane rotation, and two for bending displacements out-of-plane. (e) Small-dejection theory offlat plates The solution of the structure involves smalldeflection theory (Norris and Wilbur, 1960). This assumes that the thickness, h, of the plate is small compared with its length and breadth. In the generalized theory of elasticity, there are 6 independent components of stress necessary to define the state of stress of a body and three displacement components to define its position in space. Therefore a total of nine functions of x, y and z are necessary for solution of the elasticity problem. The smalldeflection theory (plane stress) allows the neglect of three stress components as being insignificant - (T,, r,, and 7zu- since the plates are thin compared to their length and breadth. The remaining stresses vary over the thickness of the plate. In order to avoid dealing with these distributed stresses, one may represent them by their resultant forces and moment
660
J. W. PUGH,
couples. The stress resultants N,,, NIV, and the moment couples M,,, M,, are given as follows:
R. M. ROSE and E. L. RADIN
N,, and Mvy and
The stress resultants and moment couples are given in force per unit length (lb/in.) and moment per unit length (in lb/in.) respectively. The displacements in the x-, y- and zdirections are given by u,, v, and W, respectively. The results of the STRUDL II analyses are given in the form of stress resultants and moment couples at the centroid of each element and displacements of each joint. Strains are presented in a calculation from the displacements. Reactions at the supports are given as forces and moments. In addition, the angular rotations in radians for each joint are given as three components as, CX,,(Y,for rotations about the X, y and z axes, respectively. The problem for thin plates thus separates into two solutions: the membrane problem involving N,,, N,,, N,, and the membrane displacements uOr v,,, w,; and the transverse bending problem involving the moment couples M,,, Myyr M,, (Fig. 2). The general solution to the problem is then represented by superposition (Fig. 3) of these two solutions (Kummer, 1966). From equations (4, 5 and 6), the membrane stresses are approximately given as follows: u,: = N&h
(10)
a!4 = N,,lh
(11)
7z-u = Nmlh.
(12)
Fig. 2. Stress resultants N and moment couples h4 of the type calculated in a STRUDL II analysis. This is an infinitesimal segment (dx, dy, h) of the thin sheet in plane stress. P
P
Fig. 3. The state of stress across the h-dimension of the thin sheet can be resolved into two components, the membrane stress cm and the bending stress CT&.
The loads for each member are just the stress resultants multiplied by the breadth of the plane: P, = N,, . b
(13)
P,=
(14)
N,,.b.
The stresses associated with the transverse bending of the plate are given by the following expressions (Savin, 196 1): (+
I
=$,f h3
X.X
(13
MECHANICAL u
?J
7x.v
BEHAVIOR
=122M h3
YY
=_122 Mmh3
(16)
(17)
The maximum values of the bending stresses are at values of z = +-h/2, which is at the surface of the plate. GEOMETRY OF THE MODEL
A structural model for trabecular bone was developed that, it is hoped, best incorporates the salient structural features of the real material and exhibits the same mechanical behavior, yet is soluble through the use of finite-element techniques. In developing this model, effort was made to have a structure very similar to the trabecular bone in the human knee. In particular, it was decided to model a small section of trabecular bone just adjacent to the cartilage in the distal end of the femur. As shown in Fig. 4, the model simulates the structure of the subchondral plate and the trabeculae immediately beneath the plate. From studies of longitudinal thin sections through human knees, it is evident that most trabeculae oriented are generally perpendicular to the subchondral plate with only a few degrees deviation, and that the thickness of the subchondral plate is roughly the same as a trabecula. However, some of the
OF TRABECULAR
BONE
661
trabeculae are inclined at a large angle to the vertical. The model therefore contains 5 vertically-oriented trabeculae and 4 inclined trabeculae for investigation of the behavior of these elements. The vertical trabeculae act as stiff supports because of lower internal bending moments. To reduce the average computed stiffness of the model to that of real trabecular bone, it was found necessary to have a slightly smaller volume fraction of trabeculae (larger trabecular spacing and thinner trabeculae) in the model as compared to the real bone (see Table 1). [It is quite possible that this is a reflection of the fact that the elastic moduli of the trabeculae are really less than those of cortical bone, as suggested by recent work (Pugh er al., 1973)]. The plane in Fig. 4 labelled 10 is the subchondral plate. Planes l-9 are the supporting trabeculae, and plane 11 simulates the lateral stabilizing trabeculae observed in longitudinal sections of trabecular bone. In scanning many thin sections of trabecular bone, it was decided to use a length of 0.1 in. for the distance between the subchondral plate and the lateral trabeculae, since this was about the scaled spacing between lateral trabeculae in the human knee. The area of the subchondral plate, plane 10, is O-4 X 0.06 in. To approximate the effect of and to allow investigation of the response of inclined trabeculae, planes 2
Fig. 4. Schematic of model. Plate 10 is the subchondral plate. Plates 1-9 are the supporting trabeculae. Plate 10 is loaded uniformly with 100 psi.
662
J. W. PUGH.
R. M. ROSE and E. L. RADIN
and 8 are tilted at an angle of 24” from the vertical, and planes 4 and 6 at an angle of 37”. Angles of this size have been observed in thin longitudinal sections. The portion of the model below plane 11 is a mirror image of the structure above plane 11, but is only one half of that structure. The model is thus symmetrical about plane 5. In the computation the model was considered to be loaded uniformly on the subchondral plate with a stress of 100 psi. This is based on an estimation of a weight-bearing area in the human knee of 2 in.2 and a body
(b)
Fig. 5 a, b. Plates 10 and 1 respectively as broken up into PBST finite elements. The cross-hatched regions represent the intersections of out-of-plane sheets of bone.
Table 1. Comparisons between model and human trabecular bone Model E, psi Thickness of trabeculae; in. Trabeculae/in. Vv bone CT Loading Properties of constituent material
Human (Pugh, 1972)
105
105
0.003 22.5 0.11 0.5 100 psi Same as compact bone
0~003-0~0079 316-63.5 0.2-0.47 o-9-0.55 100 psi (estimated) Unknown. But probably close to that of cortical bone, by assumption.
E=~x 106psi UCS = 20,000 psi UTS = 18,000 psi Y= 0.3
weight of 200 lb. It is assumed all the body weight is being carried on one leg, for example, during walking. The thick layer of cartilage in the knee is capable of evenly distributing the load to the subchondral plate (Frost, 1960) so the assumption of even loading of the subchondral plate is justified. To render the model soluble by STRUDL II, each of the planes numbered l-1 1 had to be broken up into triangular elements as in Fig. 5. The accuracy of the solution improves as the number of elements increases, but the practical limitation on number of elements is economy, as the cost to analyze a structure is roughly proportional to the square of the number of elements. This model has 244 nodal joints and 4 14 elements.
All joints were free to move in three-space except the supports, (the bottom ends of planes l-9). The joints at these positions were fixed in space. This was necessary to render the structure statically determinate. Each joint was rigid; that is bending moments were transferred across them. The joints at the support points also transmitted bending moments. Such joints are mechanically identical to those between real trabeculae. A program for a very simple structure is given in the Appendix as an example of a typical run. RESULTS
a. Stresses
AND DISCUSSION
on trabeculae
Figure 6 shows some of the points in the model at which computations of the stresses
MECHANICAL
BEHAVIOR
OF TRABECIJLAR
663
BONE
Uniformly applied stress
III
H
I
11
J
K
_+ _+ -o+
_-
I
LMN
_+ _+ _-
iI
0
_+
PO
i
_+ _-
-
Side view
8
Plane 23
il
4
67
5
0
9
Fig. 6. Key to Table 2. This is a side view of the model with stress positions labeled. The tensile (+) and compressive (-) sides of the sheet at each position are labeled.
and moments were made; Tables 2 and 3 show the values which resulted. It is apparent that, in the supporting trabeculae immediately under the subchondral plate, the highest bending moments exist where the trabeculae contact the plate. Indeed, a free-body diagram on a trabecula inclined at an angle CYfrom the perpendicular to the subchondral plate and loaded vertically has bending moments as a function of vertical distance x along the trabecula as follows: Mh= Ptan
(
+--.X
(1X)
i
where P is the vertical load and L the length of
the trabecula. This equation shows that a zero bending moment exists at x = L/2, and that M is a maximum for x = 0 and L. From the computer output (Fig. 7) the maximum bending moment occurs at the subchondral plate. The plot is not symmetrical along the length of plate 6 because of effects of the transverse plate present at the bottom and the inclined support beneath this. These data generally agree with the prediction of equation ( 18). Using the computer-generated data for point E and equations ( 10 and 15), the average normal compressive stress (membrane stress) in the plate is calculated as 2600 psi (see Table 2). the maximum bending stresses are H230
Table 2. Results of normal stress computation
Point A B
C D E F
G H
I J
: M N 0
Element number
Normal stress resultant (lb/in.)
N load (lb)
Normal moment couple M (in. lb/in.)
-7.8 -8.5 -6.1 2.2 2.2 2.1 - 2.2 I .I
-0-366 O-51 -0.47 +0.084 -0.47 -0.51 -0.366 0.132 0.132 0.126 -0.102 0.132
-4-7 x 10-Z’ -3.6 x IO-” 7.8 x IO-:’ I.3 x IO-” -7.8 x 10-a 3.6 x IO-” 4.7 x IO-3 -6.0 X 10-a I.3 x IO-’ I.7 x IO-” -4.8 I.8 Xx IO-” 10-y
-1.7 2.2 2.1 2.2 2.2
-0.102 O-132 0.126 O-132 O-132
1.8 x -4.8 x 7.7 x I .3 x -6.0 X
33 63 93 123 I53 183 213 219 286 291 298 303
-6.1 -8.5 -1.8
310 315 322 327 334
1.4
IO-’ IO--’ lo-:’ 10-z 10md
Normal membrane stress (psi) - 2040 -2830 -- 2600 1461 -2600 -2830 -2040 +733 +733 +700 +733 --567 -567 i733 +700 +733 +733
Max. normal bending stress (psi) 73140 i2400 25200 kO.867 35200 t2400 23140 T400 28670 -t5130 2 T320 12000 t12000 7320 t5130 k8670 7400
Normal outer-fiber stress (psi) -5180,
ill00
-430 -5230, +2600, -7800 +468, +466 -7800. +2600 -430, -5230 +1100. -5180 +333. +1133 +9403. -1937 +5830. -4430 +413, + 1053 +ll433, -12567 +11433, -12567 +413, + 1053 +5830, -4430 +9403, -1937 +333. +1133
664
J. W. PUGH,
R. M. ROSE
and E. L. RADIN
Table 3. Results of shear stress computation
Element number
Point A B C
D E F G
H I J K L M N 0
Shear stress resultant (lb/in.)
4Y,
Max. shear bending stress (psi)
-43 -80 -63 77
5133 267 a73 266
Shear outer-fiber stress (psi)
-0.13 -0.24 -0.19 0.23
153
-0.19
-1.1 x 10-4
-63
773
-136,
+10
183 213
-0.24 -0.13
-1.0x 10-a -2.0 x IO-4
-80 -43
767 7133
-147, -176,
-13 +90
279 286 291 298 303 310 315 322 327 334
0.068 -0.020 -0.0036 -0.092 0.019 -0*019 0.092 0.0036 0.020 -0.068
-2.5 x 3.6 X -3.4 x 2.9 x -3.6 x 3.6 x -2.9 x 3.4 x -3.6 x 2.5 x
23 -7 -1 -31
71670 k2400 72270 ?I940 T2400 k2400 i1940 +2270 72400 a1670
-1647, +2393, -2271, +1909, -2394, t2394, -1909, +2271, -2393, +1647,
+ 1693 -2407 +2269 -1971 +2406 -2406 +1971 -2269 +2407 -1693
2.0 1.0 1.1 9.9
x x x x
IO-4 10-d 10-a 10-S
10-S IO+
10-Z 1O-3 lo+ 1O-3 lo-$ 10-S 1O-3 10-S
r; in.
0 in/lb per in. x IO3
5
Fig. 7. Bending moment couple M,, as a function of distance from the horizontal supporting plate 11 for the inclined trabecula labeled plate 6. Note that the maximum occurs at the point of contact with the subchondral plate (y = 0.1). psi. Therefore, a compressive a tensile stress surface. These fracture stress
Shear membrane stress (psi)
33 63 93 123
I
-5
Shear moment couple M (in. lb/in.)
on one surface of the trabecula, stress of 7830 psi exists, while of 2630 psi exists on the other stresses are well below the of compact bone, and corre-
-: 31 1 7 -23
+90, -13, +10, + 143,
-176 -147 -136 +11
spond to the static load (but not the muscle forces) on this structure in uivo. For the same point, a similar calculation reveals that the shear stresses are -136 and +lO psi in the outer fibers of the trabecula (Table 3). Notice that the shear stresses are insignificant in the vertical trabeculae (points A-G) compared to the normal stresses. Thus, their neglection is justified. For the subchondral plate (points H-Q) the shear stresses are in some cases comparable to the normal stresses and should be considered. It is interesting to note the relatively large contribution of the bending component to the state of stress in the plates. This tends to verify the statement that plate-bending is an important consideration in a treatment micro-elasticity of trabecular bone.
of the
Uniformly applied stress
Fig. 8. Key to Table 3. Side view of the model with joints labeled.
MECHANICAL
b. Sti$ness
and macroscopic
BEHAVIOR
OF TRABECULAR
elastic modulus
Gross macroscopic measurements of the Young’s modulus of trabecular bone taken from subchondral regions generally result in values of 105psi and sometimes lower. In the framework of this paper, such ‘modulus’ measurements are really total stiffness measurements made on the multiply-connected structure called trabecular bone. Besides establishing that the results of the calculation agree with experiment, stiffness on a local scale may be important to the joint lubrication mechanism; both local and macroscopic stiffness are important to the function of trabecular bone as a shock absorber. The latter appears to be a key to the etiology of esteoarthritis (Radin, Paul and Rose, 1972) and is directly related to the stiffness since shock absorption in a purely elastic material has been shown to depend essentially on stiffness and maximum deflection capability (Crede, 1965). To compute the effective local modulus in our model the displacements are determined via equation (2) and the vertical displacements of the points in Fig. 6 are tabulated (see Table 4); these are divided by the total height of the model (0.15 in.) to obtain the strain; the stress (100 psi) is then divided by the strain to obtain the local modulus which is also tabulated for points A-M in Table 3. The local modulus is seen (in Table 4) to range from 6 X 103psi to 8.3 X 105psi, with an effective overall modulus (as predetermined above, see section on Geometry) of IO5psi*. The moduli at the edges of the model (8.3 X IO5psi) are probably not representative of the in viva situation due to geometry; however the large spatial variation of the modulus and the small positive deflection of the central region G may have more significance. It is known that in articular joints some areas of the bone are stiffer than others and that the cartiiage tends to compress into the less stiff areas and away
BONE
665
Table 4. Deflections
Point A B
c D E F
G H
I J K
L M
Joint number 2 30 54 58 82 86 114 142 166 170 194 198 226
Deflection X IO4{in.)
Stiffness x IO-5 psi
-0.18 8.33 -1-S 0.83 -12-O 0.12 -1.9 0.79 -27.0 0.06 -3.0 0.5 upheaval (small positive deflection) -3.0 0.5 -27.0 0.06 -1.9 0.79 -12.0 0.12 -1.8 0.83 -0.18 8.33
from the stiffer ones (Frost, 1960). Our model implies that this behavior is mechanically imposed on the softer cartilage by the behavior of the subchondral bone itself. It is important to note that the constant-stress boundary condition was chosen for our model to conform to the in vivo situation, where the load is transmitted to the subchondral plate by cartilage and synovial fluid; by way of contrast, testing machines are made of much harder and stiffer materials (e.g. steel) and impose a constant strain boundary condition on the bone surface. The latter is unphysiological (unless a surgical implant is present). Another consequence of Table 4 is that the absorption of shock will not be uniform. but will be concentrated in the low-modulus regions. c. Localfracture and collapse It must be remembered that the model was only loaded with a stress roughly equivalent to static body weight, with no consideration of muscle forces or the dynamic forces developed by walking or running. Therefore, the results presented in this section must be interpreted as minimum calculations. If the stress on the
_ *If our model were sectioned and statically tested (e.g. in an Instron machine) the overall measured stiffness would simply be determined by points A, B, D, F, G, H, J, L and M. the points of greatest stiffness. Neglecting points A and M because of edge effects. the ‘Instron’ stiffness is very close to IO5psi.
666
J. W. PUGH,
R. M. ROSE and E. L. RADIN
subchondral plate is increased by a factor of four, for example, fracture of some trabeculae would be expected, since the ultimate strength would be exceeded. To consider in more detail what is likely to happen, the discussion will be subdivided, and the three main issues, buckling, stress concentration and fatigue failure, discussed separately. Buckling. If a sheet is loaded uniaxially in compression, it can fail by buckling. The criterion for buckling can be phrased best by the following question: When a system is slightly disturbed from its equilibrium position, does it return to equilibrium or tend to depart even further? The critical uniaxial compressive load (Roark, 1943) for buckling of a plate similar in configuration and support to the vertical trabeculae (planes 1-9) in the model is given by p = 4 n2Ebh3 1212
(19)
where b is the breadth of the plate, h the thickness, and 1 the length-O-06, O-003 and 0.10 in. respectively -for the plates in the model. Using equation (19) one calculates a buckling load of 0.8 lb for the idealized plate. The stress of 100 psi on the subchondral plate is distributed evenly over an area of 0.06 X 0.4 in. The total load on the structure is 2.4 lb. Since there are 9 supporting trabeculae, each must support roughly 0.267 lb. This is below the critical buckling load for an individual trabecula. However, the calculations (see Table 2) show that each trabecula does nor support the same load. The highest load supported is that at point F. The effective load at this point is 0.51 lb. The. minimum load in any element in this trabecula is 0.45 lb. Under dynamic conditions, buckling of this trabecula is to be considered a possibility. The bending moments (Fig. 7) in the trabecula are such that the trabecula is deformed into an ‘S’. This would tend to alter the buckling behavior somewhat from that predicted by equation (19).
The fixity factor (coefficient of (,rr2Ebh3)/122) in equation ( 19) varies typically from l/4 to a maximum of 4 for different configurations of the plates (Drucker, 1967). The calculations in this section for a fixity factor of 4 can thus be regarded as an upper limit for the critical load for any geometry encountered. Of course in the real material the sheets are stabilized laterally at irregular intervals, so buckling in vim under dynamic loading is probably not as likely as indicated by the results of these calculations; however local buckling of certain trabeculae appears to be unavoidable. Effect ofstress concentrators. The structure of trabecular bone is known to consist of sheets with holes in them (Dyson et al., 1970). The holes provide interconnecting passageways for marrow and blood vessels. Therefore, the presence of stress concentrators as above is to be the rule rather than the exception. Additional types of stress concentrators present in trabecular bone could be of importance. The canaliculi, lacunae, and the lamellar structures are all present in the sheets of bone. The canaliculi and lacunae are so small and localized in that they do not extend all the way through the sheet of bone, that they would probably be of less importance, Lamellae that end inside the structure could cause obvious problems at their points of termination, depending on the mechanical properties of the interfaces between the lamellae. These intercould provide effective faces, however, barriers to crack propagation (Piekarski, 1970). The effect of such stress concentrators on the membrane stresses and bending stresses can be best summarized by the following two equations which follow directly from the standard formulae, (Savin, 1961): a:,=
a,(l+2K)
u; = fqo ( 1+ 0*79K)
(20) (21)
where K = c/a, the ratio of semi-axes of an elliptical hole and u, and ob are the membrane and bending stresses respectively. AS an example of stress increases due to typical
MECHANICAL
BEHAVIOR
OF TRABECULAR
BONE
667
concentrators, for a membrane stress of -26OOpsi and K = 2, uI, = -13,000psi. For bending stresses of 5230 psi and K = 2, (r; = + 13,500 psi. Therefore on one surface of the trabecula (point E in Fig. 5) there exists a tensile stress of 500 psi, while on the other surface, a compressive stress of 26500 exists, which exceeds the compressive strength of the bone. Local fracture would therefore be expected to occur at the end of the major axis of the elliptical hole on the compressive side of the trabecula (the positive-Z side of the plane in Fig. 9). This fracture itself causes stress concentration, and, under cyclic loading, crack propagation would be expected. Accepted yield criteria (McClintock and Argon, 1966) can be used to investigate the question of generalized yielding in the model. One such criterion is the following: * =
[H((T,- ‘T,J2+ (u,-~J2+
(~u-~d21
+ 3(.r;.u+Tsz+?y~]1’2.
(22)
When v equals the yield strength of the material, yielding occurs. As typical example, point H has the following maximum stress components: (T,= 1133psi v,~ = 3473 psi CX.U= 1693 psi and @ = 4250 psi. which is much less than the yield strength. In view of the fact that the shear stresses are negligible in the vertical trabeculae, and the contribution from the above calculation only minor in the subchondral plate, the use of such yield criteria in an analysis of this model is probably not justified. A good indication of the yield behavior is thus gained through reference only to the normal stresses. Fatigue fracture. The fatigue life of compact bone subjected to zero minimum stress as in viva (as opposed to zero mean stress) has been been shown to vary from lo3 cycles at a maximum stress of 11,700 psi to lo7 cycles at a
Fig. 9. Position of local fracture around macroscopic stress concentrator in an inclined trabecula.
maximum stress of 8,750 psi (Freeman, Swanson and Day, 1971). Reference to the outerfiber normal stresses in Table 2 shows that regions of vertical trabeculae would have a dynamic fatigue life of less than lo3 cycles if the dynamic stresses are only twice as great as the static stresses. (lo3 cycles represents a brisk 20-minute walk at roughly 2 steps per sec.) During this relatively short fatigue life, the natural turnover rate of the bone probably has no effect in healing incipient fatigue cracks. If, in addition, mild (K = 1) stress concentrators such as macroscopic holes are present, the occurrence of fatigue fractures would be more probable. Notice that even higher normal stresses exist in the subchondral plate (points H-Q) than in the vertical trabeculae. Even under static conditions, certain regions of the subchondral plate will have a fatigue life of less
668
J. W. PUGH.
R. M. ROSE and E. L. RADIN
than lo3 cycles. Considering the relative absence of macroscopic stress concentrators in the subchondral plate as compared to the trabeculae, the incidence of fatigue fracture is probably roughly the same for all regions of subchondral bone. IMPLICATIONS FOR THE ROLE OF CANCELLOUS BONE IN JOINT DEGENERATION It has already been pointed out above that the stiffness at the surface of the subchondral plate varies considerably with position according to the arrangement of the underlying trabeculae, and that this variation leads to nonuniform compression of the cartilage under load. It is possibly significant in this connection that scanning electron microscopy studies (Walker et al., 1969) show that, with age, cartilage develops large wavelength undulations; although no examination of the underlying cancellous bone was apparently made in the above study, the spacing of the cartilage undulations happens to correspond approximately to the average intertrabecular spacing. It is important to note that the least compliant regions of the cancellous bone also absorb less shock. Thus the cartilage above these regions is compressed more and is exposed to more shock loading; it is in these regions, then, that joint degeneration may begin. CONCLUSIONS
A quantitative model having the elastic modulus of trabecular bone can be constructed by approximating the trabecular geometry with a hypothetical network of compact bone, and then using finite-element analysis with a high-speed computer. Because the number and thickness of trabeculae in the model must be somewhat less than they are in normal trabecular bone in order to have equivalent stiffness, there are two implications: first, that the intertrabecular soft tissue is probably of no importance to the mechanical properties, and second, that the compact bone which constitutes the trabeculae is probably not as stiff as cortical bone, al-
though we previously assumed that it was (see Introduction). The computations themselves lead to several conclusions: (a) The in vivo elastic deformation of trabecular bone must involve considerable bending and buckling of trabeculae. (b) The bending (or buckling) stresses on trabeculae can be comparable to or greater than the membrane stresses. (c) For the loads due to normal activity, fatigue fracture seems inevitable in some trabeculae, even in normal, healthy bone; for some trabecular orientations the fatigue life may be lo3 cycles or less. (d) For any geometry similar to that of the model, the stiffness is spatially dependent. The nature of the spatial dependence depends totally on the trabecular arrangement, and large variations in stiffness may occur over a single intertrabecular spacing. Such behavior is totally suppressed by accepted mechanical testing methods, but should be very manifest in a normal articular joint, in at least two ways: the cartilage will be more severely compressed in regions where the subchondral bone is stiffer, and the transmission of shock loads to the cartilage will be much greater in the same region. Thus the spatial inhomogeneity of trabecular bone may in fact play some role in the functioning (and/or degeneration) of articularjoints. Acknowledgement-This
work was sponsored
by the
Office of Naval Research. REFERENCES Burstein, A. H., Shaffer, B. W. and Frankel, V. H. (1970) Elastic analvsis of condvlar structures. ASME Publication 70-WAjBHF-I. Chen, P. E. (1971) Strength properties of discontinuous fiber composites. Polym. Engng Sci. 11,5 l-56. Crede, C. E. (1965) in Shock and Vibration Concepts in Engineering Design, p. 128. Prentice-Hall, New Jersey. Currey, J. D. (1970) The mechanical properties of bone. Clin. Orthop. 73,2 1O-22 1. In Dorland’s Illustrated Medical Dictionary, (1965) p. 1598, Saunders, Phila. Drucker, D. C. (1967) In Introduction to Mechanics of DeformableSolids. Chap. 16. McGraw-Hill, New York.
MECHANICAL
BEHAVIOR
Dyson, E. D., Jackson, C. K. and Whitehouse, W. J. (1970) Scanning electron microscope studies of human trabecular bone. Nature 225.957-959. Elias, H. (1971) Three-dimensional structure identified from single sections. Scient 174,993-1000. Evans, F. G. (1957) In Stress and Strain in Bones, p. 176-202. Thomas, Springfield, Illinois. Evans, F. G. (1961a) Relation of the physical properties of bone to fractures. ASOS Instruc. Coursefor Lecfurers l&110-121. Evans, F. G. (1961b) Biomechanical Studies of the Musculo-Skeletal System, Chap. 3. Thomas, Springfield, Illinois. Freeman, M. A. R., Day, W. H. and Swanson. S. A. V. ( 197 1) Fatigue fracture in the subchondral bone of the human cadaver femoral head. Med. Biol. Engng 9, to joint biomechanics.
Henry Ford Hosp. Med. Bull. 8,415-432.
Hackett, R. M. (1971) Viscoelastic stresses in a composite system. Polym. Engng Sci. 11,220-225. Kummer, B. (1966) Photoelastic studies on the functional structure of bone. Folio Biotheoretica 6.3 1. Logcher. R. D. et al. (1968) ICES STRUDL I1 Vol. IFrume Analysis. Research Report R68-91, Department of Civil Engineering, M.I.T. Cambridge, Mass. 02 139; (1971) ICES STRUDL II Vol. II-Additional Design and Analysis Facilities.
F. A. and Argon, A. S. (1966) Mechanical pp. 276-279. Addison-Wesley. Reading, Mass. Research Report R70-77, Department of Civil Engineering, M.I.T..Cambridge,Mass. 02139. McDonnell, EC1 ICES STRUDL User’s Manual for PBST Finite Element. Engineering Computer International, Inc., 1073 Massachusetts Ave., Cambridge, Mass. 02139. McElhaney, J.. Alem, N. and Roberts, V. (1970) A Porous Block Model for Cancellous Bone. ASME Publication 70-WA/BHF-2. McLeish. R. D. and Habboobi (1972) Strain gauge techniques for cadaveric bone. Engng Med. 1,37. Norris, C. H. and Wilbur, J. B. (1960) Elementary Structural Analysis, pp. 558-566. McGraw-Hill, New York. Piekarski, K. (1970) Fracture of bone. J. uppl. Phys. 41. McClintock,
of Materials,
215-223.
Piotrowski. G. and Wilcox, G. A.,Jr.( 1971)TheSTRESS program: A computer program for the analysis of stress in long bones. J. Biomechanics 4.497-506. Pugh, J. W.. Runkle. J. C., Rose, R. M., Paul, I. L. and Trabecular Bone. Ph.D. Thesis, M.I.T., Cambridge, Mass. 02139. Pugh, J. W., Runkle. J. C.. Rose, R. M., Paul, I.L. and Radin, E. L.. Determination of the elastic modulus of individual trabeculae. To be presented at 1973 Meeting of Orthopnedic
Research
Society.
Radin, E. L. and Paul, I. L. (1969) Failure of synovial fluid to cushion. Nuture 222,999-1000. Radin. E. L.. Paul, I. 1,. and Rose, R. M. (1972) Hypothesis: Role of mechanical factors in the pathogenesis of primary osteo-arthritis. Lancet 1.5 19-522. Roark. R. J. ( 1943) In Formulas for Stress and Strain, p. 793. McGraw-Hill. New York.
BONE
669
Savin, G. N. (196 1) Stress Concentration Around Holes. pp. 73.305,32 1. Pergamon Press, Oxford. Swanson, S. A. V. and Freeman, M. A. R. (1966) Is Bone Hydraulically Strengthened? Med. Biol. Engng 4, 433-438.
Walker, P. et al. (1969) Lubrication Mechanism in Human Joints: A Study by Scanning Electron Microscopy. In Lubrication and Wear in Joints (Edited by Wright, Verna), pp. 49-56. Lippincott, London. Yamada. H. (19701 Strenath of Bioloeical Materials. pp. 19-34. Williams &nd Wilkins: Bait&ore. Zienkiewicz, 0. C. ( 197 1) The Finite-Element Method in Engineering Science. McGraw-Hill, London.
;
619-629.
Frost, H. M. (1960) Introduction
Behavior
OF TRABECULAR
NOMENCLATURE semi-minor axis of ellipse (in.) breadth (in.) semi-major axis of ellipse (in.) contiguity of trabeculae elastic modulus or stiffness (psi) force (lb) thickness (in.) ratio of semi-axis of an ellipse length (in.) length (in.) moment couple (in. lb/in.) bending moment (in. lb) stress resultant (lb/in.) load (lb) stiffness coefficients (lb/in.) stress coefficients (Ib/in.3) displacements in x-, y-. z-direction (in.) ultimate compressive strength (psi) ultimate tensile strength (psi) volume fraction (in.3/in.3) Cartesian coordinates (in.) angle (degree or radians) nodal displacement (in.) Poisson’s ratio normal stress (psi) yield strength under combined stresses (psi) stress due to presence of stress concentrator
C’; E F h K
1 L M Mb N P S T uo. vo. w0 ucs UTS V, x, VTz a 6 lJ rr ;r (7’
(psi) bending stress (psi) membrane stress (psi) shear stress (psi) column matrix square matrix. APPENDlX STRUDL ‘SIMPLE PROBLEM’ ‘BONE ANALYSIS’ TYPE SPACE FRAME DEBUG ALL DUMP TIME JOINT COORDINATES 1 1. 0. 0. s 2 0. 0. 0. s 3 1. 0. I. 4 0. 0. 1. 5 0. 1. 1. 6 1. I. 1. 7 1. I. 0. s
670
J. W. PUGH,
R. M. ROSE and E. L. RADIN
8 0. 1. 0. S ELEMENT INCIDENCES 1321 2234 3354 4536 5687 6865 UNITS LBS INCHES ELEMENT PROPERTIES 1 TO 6 TYPE ‘PBST’ THICKNESS 0.003 CONSTANTS E 2000000.0 ALL POISSON 0.3 ALL LOADING 1 ‘TOP’ ELEMENT LOADS 3 TO 4 SURFACE FORCES GLOBAL PZ -100.0 GLOBAL SURFACE FORCES STIFFNESS ANALYSIS WITH REACTIONS NJP 3
LIST DISPLACEMENTS, FORCES, REACTIONS, FINISH
STRESSES, ALL
z
100 psi
t
I
STRAINS.
1973. Vol. 6, pp. 657470.
Pergamon Press.
Printed m Great Britam
A STRUCTURAL MODEL FOR THE MECHANICAL BEHAVIOR OF TRABECULAR Department of Orthopaedics,
BONE*
J. W. PUGH Hospital forJoint Diseases and Medical Center, New York, U.S.A.
R. M. ROSE Department of Metallurgy and Materials Science, Massachusetts Cambridge, Massachusetts, U.S.A.
Institute of Technology.
and E. L. RADIN Department of Orthopedic Surgery, Harvard Medical School. Childrens Hospital Medical Center. Boston, Massachusetts. U.S.A. Abstract- A quantitative model is developed for trabecuiar bone by approximating the trabecular geometry with a hypothetical network of compact bone. For the region immediately beneath the articular cartilage in the distal end of the femur, finite element analyses were performed with a high speed computer, assuming a physiological static load. The results indicate that bending and buckling of trabeculae are considerable in any elastic deformation of the bone: that fatigue fracture in some fraction of suitably oriented trabeculae is inevitable in normal ambulation; and that the stiffness varies considerably with lateral position across the subchondral plate. The latter depends totally on trabecular arrangement and may play a role in joint function and degeneration. The adjustments necessary to bring the gross stiffness into agreement with experiment imply that the intertrabecular soft tissues are of no consequence to the mechanical properties and that the compact bone of which trabeculae are made is probably not as stiff as cortical bone.
INTRODUCTION
AND BASIC ASSUMPTIONS
(a) Organization A WORKABLE semi-empirical model based on porous block geometry has been developed to explain the variation of the elastic (Young’s) modulus and mechanical strength of cancellous bone as a function of density (McElhaney, Alem and Roberts, 1970). Despite the crudity of the geometrical assumption, the model worked well when applied to cancellous bone from the cranium and vertebral bodies. Another assumption used in the above model, one which will be used as well in this paper, is that the properties of the trabecular material are identical to those of compact (cortical) bone. Such an assumption. that mature hard tissue has uniform microscopic properties, seems reasonable but is in fact unproven. *Rc,ccGvrd 29 Novrmber
1972.
(b) Mechanical properties of bone The accepted value for Young’s modulus of compact bone is 2 x lo6 psi and its strength properties are approximately 20,000 psi in compression and 18,000 psi in tension (Yamada, 1970; Currey, 1970; Evans, 1961a; Evans, 1957), the yield strength and the fracture strength being essentially the same. Poisson’s ratio for cortical bone is very close to 0.3 (McLeish and Habboobi, 1972). Trabecular bone has a stiffness of about 10” psi (McElhaney, Alem and Roberts. 1970) and a compressive strength of 500-1000 psi depending on where in the femur the bone is taken (Evans, 1961 b). Experiments have shown the mechanical properties of trabecular bone to be virtuahy unrelated to the fluid and fatty tissue present in the intertrabecular
658
J. W. PUGH, R. M. ROSE and E. L. RADIN
spaces (Swanson and Freeman, 1966). It is therefore reasonable to neglect, in a quantitative mechanical model, the presence of the intertrabecular soft tissues. Cartilage
CC) Structure
and mechanical function of trabecular bone The general structure of trabecular bone consists of a system of sheets of compact hard tissue, with a distribution of widths, the lower limit corresponding to rods. Typical photomicrographs that one observes of trabecular bone (Elias, 1971) can correspond stereologically only to such a distribution of interconnected sheets. Rodlike or beamlike trabeculae would appear very differently in cross-section. The sheet- and rod morphology has been confirmed for vertebral bodies by scanning electron microscopy (Dyson, Jackson and Whitehouse, 1970). The SEM also reveals that the open spaces present in the trabecular bone are interconnected. Considering condylar structures, it would seem that the function of the trabecular bone is to support and direct the operational loads through the ends of the long bones to the diaphyses. Evidence for the central role of the trabeculae in this function can be found in photoelastic studies which show that they are aligned in the directions of the normal stress trajectories (Kummer, 1966). As an example, in the human knee (Fig. l), the femoral condyles rest on the tibia1 plateaus and serve to transmit the body weight to the lower leg. The condyles are roughly cylindrical structures, and consist of a thick layer of cartilage, a thin layer of cortical bone referred to as the subchondral plate, and the region of trabecular bone immediately underneath the plate. The amount of trabecular bone present is large compared to the amount of cartilage and the thickness of the subchondral plate. The region between the cartilage of the condyles and plateaus is filled with synovial fluid, which serves no mechanical function other than lubrication (Radin and Paul, 1969). A mathematical model for the elastic proper-
-
Synwial fluid -Cartilage
I
Joint space
1
Fig. 1. Schematic of human knee joint. The cartilage is very thick compared to a trabecular thickness. The subchondral plate is roughly the same thickness as a trabecula.
ties of condylar structures (Burstein, Shaffer and Frankel, 1970) has shown that, at least under static conditions, the trabecular bone in the condyle is important in load-carrying; 80 per cent of the surface loading is transmitted through the trabecular region. The subchondral plate, furthermore, was found to be loaded with a combination of bending and membrane stresses. Finally, from the difference in properties between cortical and trabecular bone, it was concluded that, as a percentage of the mechanical strength, the stresses in the trabecular region are large compared to those in the cortical bone of the shaft of the femur. The plate-like trabecular geometry means that bending and buckling of individual elements can not be neglected even for small stresses, and therefore we must include bending and buckling in our computations. (d) Technique of computation The use of the computer to solve biomechanical problems is not new. Probably the most recent use is a program for the analysis of stresses in long bones (Piotrowski and Wilcox, 1971). This allows calculation of stresses and strains through the use of standard stress-strain relations and the Airy stress function for a bone consisting of any cross-
MECHANICAL
BEHAVIOR
sectional shape and variation of cross-section with length. An alternative which offers greater possibilities, as long as a high speed computer is available, is the use of finite-element techniques (Zienkiewicz, 197 1) which incidentally has led to a greater understanding of the mechanical behavior of composites (Hackett. 1971; Chen, 197 1). It is also applicable to a number of other systems and lends itself especially well to the analysis of very complex structures. The structural system is considered as an assemblage of idealized elastic elements joined together at discrete nodes. Using the direct stiffness concept, which allows adding together at each node the stiffness coefficients of adjacent elements, a stiffness matrix for the structural system is obtained. This matrix relates the external forces acting on the nodes to the displacements of the nodes: IF) =
[Sl{Sl.
(1)
An inversion of the stiffness matrix [S] allows calculation of the influence matrix, which gives the nodal displacements (6) as functions of the external loads acting on the structural system: IS> = [S]-‘(F).
(2)
This is solved for the nodal displacement, {S}, and once the nodal displacements are known, the solution is specified. Using the same strain pattern that we assumed in derivation of the stiffness matrix, one can derive the matrix of stress coefficients giving the stress in the nth element as a function of the nodal displacements: {a?,} = r ~?Il~hJ.
(3)
A continuum can be represented as an assemblage of elastic elements joined at nodes as above, and the equations for internal displacements and stresses can be derived using the matrix techniques outlined. The finite-element
OF TRABECULAR
BONE
659
method is in fact a solution that minimizes the total potential energy of the system (Zienkiewicz, 197 1). The system used in the present study is ICES STRUDL II (Logcher et al., 1968, 1971) developed at the M.I.T. Civil Engineering Systems Laboratory. This system allows the calculation of a stiffness analysis on a structure that is specified geometrically in conventional three-space. The stiffness analysis is a linear, elastic, static, small displacement analysis applying the finite-element method in an efficient algorithm. It requires the specification of element properties in acceptable form and treats the nodal joint displacements as unknowns. The ‘PBST’ plate bending and stretching element used in this analysis (McDonnell) is a six degree of freedom (six nodal unknowns), flat, triangular element which considers both the in-plane (membrane) and plate-bending forces. The element formulation uses the following displacement expansions: two for orthogonal stretching in-plane, two for inplane rotation, and two for bending displacements out-of-plane. (e) Small-dejection theory offlat plates The solution of the structure involves smalldeflection theory (Norris and Wilbur, 1960). This assumes that the thickness, h, of the plate is small compared with its length and breadth. In the generalized theory of elasticity, there are 6 independent components of stress necessary to define the state of stress of a body and three displacement components to define its position in space. Therefore a total of nine functions of x, y and z are necessary for solution of the elasticity problem. The smalldeflection theory (plane stress) allows the neglect of three stress components as being insignificant - (T,, r,, and 7zu- since the plates are thin compared to their length and breadth. The remaining stresses vary over the thickness of the plate. In order to avoid dealing with these distributed stresses, one may represent them by their resultant forces and moment
660
J. W. PUGH,
couples. The stress resultants N,,, NIV, and the moment couples M,,, M,, are given as follows:
R. M. ROSE and E. L. RADIN
N,, and Mvy and
The stress resultants and moment couples are given in force per unit length (lb/in.) and moment per unit length (in lb/in.) respectively. The displacements in the x-, y- and zdirections are given by u,, v, and W, respectively. The results of the STRUDL II analyses are given in the form of stress resultants and moment couples at the centroid of each element and displacements of each joint. Strains are presented in a calculation from the displacements. Reactions at the supports are given as forces and moments. In addition, the angular rotations in radians for each joint are given as three components as, CX,,(Y,for rotations about the X, y and z axes, respectively. The problem for thin plates thus separates into two solutions: the membrane problem involving N,,, N,,, N,, and the membrane displacements uOr v,,, w,; and the transverse bending problem involving the moment couples M,,, Myyr M,, (Fig. 2). The general solution to the problem is then represented by superposition (Fig. 3) of these two solutions (Kummer, 1966). From equations (4, 5 and 6), the membrane stresses are approximately given as follows: u,: = N&h
(10)
a!4 = N,,lh
(11)
7z-u = Nmlh.
(12)
Fig. 2. Stress resultants N and moment couples h4 of the type calculated in a STRUDL II analysis. This is an infinitesimal segment (dx, dy, h) of the thin sheet in plane stress. P
P
Fig. 3. The state of stress across the h-dimension of the thin sheet can be resolved into two components, the membrane stress cm and the bending stress CT&.
The loads for each member are just the stress resultants multiplied by the breadth of the plane: P, = N,, . b
(13)
P,=
(14)
N,,.b.
The stresses associated with the transverse bending of the plate are given by the following expressions (Savin, 196 1): (+
I
=$,f h3
X.X
(13
MECHANICAL u
?J
7x.v
BEHAVIOR
=122M h3
YY
=_122 Mmh3
(16)
(17)
The maximum values of the bending stresses are at values of z = +-h/2, which is at the surface of the plate. GEOMETRY OF THE MODEL
A structural model for trabecular bone was developed that, it is hoped, best incorporates the salient structural features of the real material and exhibits the same mechanical behavior, yet is soluble through the use of finite-element techniques. In developing this model, effort was made to have a structure very similar to the trabecular bone in the human knee. In particular, it was decided to model a small section of trabecular bone just adjacent to the cartilage in the distal end of the femur. As shown in Fig. 4, the model simulates the structure of the subchondral plate and the trabeculae immediately beneath the plate. From studies of longitudinal thin sections through human knees, it is evident that most trabeculae oriented are generally perpendicular to the subchondral plate with only a few degrees deviation, and that the thickness of the subchondral plate is roughly the same as a trabecula. However, some of the
OF TRABECULAR
BONE
661
trabeculae are inclined at a large angle to the vertical. The model therefore contains 5 vertically-oriented trabeculae and 4 inclined trabeculae for investigation of the behavior of these elements. The vertical trabeculae act as stiff supports because of lower internal bending moments. To reduce the average computed stiffness of the model to that of real trabecular bone, it was found necessary to have a slightly smaller volume fraction of trabeculae (larger trabecular spacing and thinner trabeculae) in the model as compared to the real bone (see Table 1). [It is quite possible that this is a reflection of the fact that the elastic moduli of the trabeculae are really less than those of cortical bone, as suggested by recent work (Pugh er al., 1973)]. The plane in Fig. 4 labelled 10 is the subchondral plate. Planes l-9 are the supporting trabeculae, and plane 11 simulates the lateral stabilizing trabeculae observed in longitudinal sections of trabecular bone. In scanning many thin sections of trabecular bone, it was decided to use a length of 0.1 in. for the distance between the subchondral plate and the lateral trabeculae, since this was about the scaled spacing between lateral trabeculae in the human knee. The area of the subchondral plate, plane 10, is O-4 X 0.06 in. To approximate the effect of and to allow investigation of the response of inclined trabeculae, planes 2
Fig. 4. Schematic of model. Plate 10 is the subchondral plate. Plates 1-9 are the supporting trabeculae. Plate 10 is loaded uniformly with 100 psi.
662
J. W. PUGH.
R. M. ROSE and E. L. RADIN
and 8 are tilted at an angle of 24” from the vertical, and planes 4 and 6 at an angle of 37”. Angles of this size have been observed in thin longitudinal sections. The portion of the model below plane 11 is a mirror image of the structure above plane 11, but is only one half of that structure. The model is thus symmetrical about plane 5. In the computation the model was considered to be loaded uniformly on the subchondral plate with a stress of 100 psi. This is based on an estimation of a weight-bearing area in the human knee of 2 in.2 and a body
(b)
Fig. 5 a, b. Plates 10 and 1 respectively as broken up into PBST finite elements. The cross-hatched regions represent the intersections of out-of-plane sheets of bone.
Table 1. Comparisons between model and human trabecular bone Model E, psi Thickness of trabeculae; in. Trabeculae/in. Vv bone CT Loading Properties of constituent material
Human (Pugh, 1972)
105
105
0.003 22.5 0.11 0.5 100 psi Same as compact bone
0~003-0~0079 316-63.5 0.2-0.47 o-9-0.55 100 psi (estimated) Unknown. But probably close to that of cortical bone, by assumption.
E=~x 106psi UCS = 20,000 psi UTS = 18,000 psi Y= 0.3
weight of 200 lb. It is assumed all the body weight is being carried on one leg, for example, during walking. The thick layer of cartilage in the knee is capable of evenly distributing the load to the subchondral plate (Frost, 1960) so the assumption of even loading of the subchondral plate is justified. To render the model soluble by STRUDL II, each of the planes numbered l-1 1 had to be broken up into triangular elements as in Fig. 5. The accuracy of the solution improves as the number of elements increases, but the practical limitation on number of elements is economy, as the cost to analyze a structure is roughly proportional to the square of the number of elements. This model has 244 nodal joints and 4 14 elements.
All joints were free to move in three-space except the supports, (the bottom ends of planes l-9). The joints at these positions were fixed in space. This was necessary to render the structure statically determinate. Each joint was rigid; that is bending moments were transferred across them. The joints at the support points also transmitted bending moments. Such joints are mechanically identical to those between real trabeculae. A program for a very simple structure is given in the Appendix as an example of a typical run. RESULTS
a. Stresses
AND DISCUSSION
on trabeculae
Figure 6 shows some of the points in the model at which computations of the stresses
MECHANICAL
BEHAVIOR
OF TRABECIJLAR
663
BONE
Uniformly applied stress
III
H
I
11
J
K
_+ _+ -o+
_-
I
LMN
_+ _+ _-
iI
0
_+
PO
i
_+ _-
-
Side view
8
Plane 23
il
4
67
5
0
9
Fig. 6. Key to Table 2. This is a side view of the model with stress positions labeled. The tensile (+) and compressive (-) sides of the sheet at each position are labeled.
and moments were made; Tables 2 and 3 show the values which resulted. It is apparent that, in the supporting trabeculae immediately under the subchondral plate, the highest bending moments exist where the trabeculae contact the plate. Indeed, a free-body diagram on a trabecula inclined at an angle CYfrom the perpendicular to the subchondral plate and loaded vertically has bending moments as a function of vertical distance x along the trabecula as follows: Mh= Ptan
(
+--.X
(1X)
i
where P is the vertical load and L the length of
the trabecula. This equation shows that a zero bending moment exists at x = L/2, and that M is a maximum for x = 0 and L. From the computer output (Fig. 7) the maximum bending moment occurs at the subchondral plate. The plot is not symmetrical along the length of plate 6 because of effects of the transverse plate present at the bottom and the inclined support beneath this. These data generally agree with the prediction of equation ( 18). Using the computer-generated data for point E and equations ( 10 and 15), the average normal compressive stress (membrane stress) in the plate is calculated as 2600 psi (see Table 2). the maximum bending stresses are H230
Table 2. Results of normal stress computation
Point A B
C D E F
G H
I J
: M N 0
Element number
Normal stress resultant (lb/in.)
N load (lb)
Normal moment couple M (in. lb/in.)
-7.8 -8.5 -6.1 2.2 2.2 2.1 - 2.2 I .I
-0-366 O-51 -0.47 +0.084 -0.47 -0.51 -0.366 0.132 0.132 0.126 -0.102 0.132
-4-7 x 10-Z’ -3.6 x IO-” 7.8 x IO-:’ I.3 x IO-” -7.8 x 10-a 3.6 x IO-” 4.7 x IO-3 -6.0 X 10-a I.3 x IO-’ I.7 x IO-” -4.8 I.8 Xx IO-” 10-y
-1.7 2.2 2.1 2.2 2.2
-0.102 O-132 0.126 O-132 O-132
1.8 x -4.8 x 7.7 x I .3 x -6.0 X
33 63 93 123 I53 183 213 219 286 291 298 303
-6.1 -8.5 -1.8
310 315 322 327 334
1.4
IO-’ IO--’ lo-:’ 10-z 10md
Normal membrane stress (psi) - 2040 -2830 -- 2600 1461 -2600 -2830 -2040 +733 +733 +700 +733 --567 -567 i733 +700 +733 +733
Max. normal bending stress (psi) 73140 i2400 25200 kO.867 35200 t2400 23140 T400 28670 -t5130 2 T320 12000 t12000 7320 t5130 k8670 7400
Normal outer-fiber stress (psi) -5180,
ill00
-430 -5230, +2600, -7800 +468, +466 -7800. +2600 -430, -5230 +1100. -5180 +333. +1133 +9403. -1937 +5830. -4430 +413, + 1053 +ll433, -12567 +11433, -12567 +413, + 1053 +5830, -4430 +9403, -1937 +333. +1133
664
J. W. PUGH,
R. M. ROSE
and E. L. RADIN
Table 3. Results of shear stress computation
Element number
Point A B C
D E F G
H I J K L M N 0
Shear stress resultant (lb/in.)
4Y,
Max. shear bending stress (psi)
-43 -80 -63 77
5133 267 a73 266
Shear outer-fiber stress (psi)
-0.13 -0.24 -0.19 0.23
153
-0.19
-1.1 x 10-4
-63
773
-136,
+10
183 213
-0.24 -0.13
-1.0x 10-a -2.0 x IO-4
-80 -43
767 7133
-147, -176,
-13 +90
279 286 291 298 303 310 315 322 327 334
0.068 -0.020 -0.0036 -0.092 0.019 -0*019 0.092 0.0036 0.020 -0.068
-2.5 x 3.6 X -3.4 x 2.9 x -3.6 x 3.6 x -2.9 x 3.4 x -3.6 x 2.5 x
23 -7 -1 -31
71670 k2400 72270 ?I940 T2400 k2400 i1940 +2270 72400 a1670
-1647, +2393, -2271, +1909, -2394, t2394, -1909, +2271, -2393, +1647,
+ 1693 -2407 +2269 -1971 +2406 -2406 +1971 -2269 +2407 -1693
2.0 1.0 1.1 9.9
x x x x
IO-4 10-d 10-a 10-S
10-S IO+
10-Z 1O-3 lo+ 1O-3 lo-$ 10-S 1O-3 10-S
r; in.
0 in/lb per in. x IO3
5
Fig. 7. Bending moment couple M,, as a function of distance from the horizontal supporting plate 11 for the inclined trabecula labeled plate 6. Note that the maximum occurs at the point of contact with the subchondral plate (y = 0.1). psi. Therefore, a compressive a tensile stress surface. These fracture stress
Shear membrane stress (psi)
33 63 93 123
I
-5
Shear moment couple M (in. lb/in.)
on one surface of the trabecula, stress of 7830 psi exists, while of 2630 psi exists on the other stresses are well below the of compact bone, and corre-
-: 31 1 7 -23
+90, -13, +10, + 143,
-176 -147 -136 +11
spond to the static load (but not the muscle forces) on this structure in uivo. For the same point, a similar calculation reveals that the shear stresses are -136 and +lO psi in the outer fibers of the trabecula (Table 3). Notice that the shear stresses are insignificant in the vertical trabeculae (points A-G) compared to the normal stresses. Thus, their neglection is justified. For the subchondral plate (points H-Q) the shear stresses are in some cases comparable to the normal stresses and should be considered. It is interesting to note the relatively large contribution of the bending component to the state of stress in the plates. This tends to verify the statement that plate-bending is an important consideration in a treatment micro-elasticity of trabecular bone.
of the
Uniformly applied stress
Fig. 8. Key to Table 3. Side view of the model with joints labeled.
MECHANICAL
b. Sti$ness
and macroscopic
BEHAVIOR
OF TRABECULAR
elastic modulus
Gross macroscopic measurements of the Young’s modulus of trabecular bone taken from subchondral regions generally result in values of 105psi and sometimes lower. In the framework of this paper, such ‘modulus’ measurements are really total stiffness measurements made on the multiply-connected structure called trabecular bone. Besides establishing that the results of the calculation agree with experiment, stiffness on a local scale may be important to the joint lubrication mechanism; both local and macroscopic stiffness are important to the function of trabecular bone as a shock absorber. The latter appears to be a key to the etiology of esteoarthritis (Radin, Paul and Rose, 1972) and is directly related to the stiffness since shock absorption in a purely elastic material has been shown to depend essentially on stiffness and maximum deflection capability (Crede, 1965). To compute the effective local modulus in our model the displacements are determined via equation (2) and the vertical displacements of the points in Fig. 6 are tabulated (see Table 4); these are divided by the total height of the model (0.15 in.) to obtain the strain; the stress (100 psi) is then divided by the strain to obtain the local modulus which is also tabulated for points A-M in Table 3. The local modulus is seen (in Table 4) to range from 6 X 103psi to 8.3 X 105psi, with an effective overall modulus (as predetermined above, see section on Geometry) of IO5psi*. The moduli at the edges of the model (8.3 X IO5psi) are probably not representative of the in viva situation due to geometry; however the large spatial variation of the modulus and the small positive deflection of the central region G may have more significance. It is known that in articular joints some areas of the bone are stiffer than others and that the cartiiage tends to compress into the less stiff areas and away
BONE
665
Table 4. Deflections
Point A B
c D E F
G H
I J K
L M
Joint number 2 30 54 58 82 86 114 142 166 170 194 198 226
Deflection X IO4{in.)
Stiffness x IO-5 psi
-0.18 8.33 -1-S 0.83 -12-O 0.12 -1.9 0.79 -27.0 0.06 -3.0 0.5 upheaval (small positive deflection) -3.0 0.5 -27.0 0.06 -1.9 0.79 -12.0 0.12 -1.8 0.83 -0.18 8.33
from the stiffer ones (Frost, 1960). Our model implies that this behavior is mechanically imposed on the softer cartilage by the behavior of the subchondral bone itself. It is important to note that the constant-stress boundary condition was chosen for our model to conform to the in vivo situation, where the load is transmitted to the subchondral plate by cartilage and synovial fluid; by way of contrast, testing machines are made of much harder and stiffer materials (e.g. steel) and impose a constant strain boundary condition on the bone surface. The latter is unphysiological (unless a surgical implant is present). Another consequence of Table 4 is that the absorption of shock will not be uniform. but will be concentrated in the low-modulus regions. c. Localfracture and collapse It must be remembered that the model was only loaded with a stress roughly equivalent to static body weight, with no consideration of muscle forces or the dynamic forces developed by walking or running. Therefore, the results presented in this section must be interpreted as minimum calculations. If the stress on the
_ *If our model were sectioned and statically tested (e.g. in an Instron machine) the overall measured stiffness would simply be determined by points A, B, D, F, G, H, J, L and M. the points of greatest stiffness. Neglecting points A and M because of edge effects. the ‘Instron’ stiffness is very close to IO5psi.
666
J. W. PUGH,
R. M. ROSE and E. L. RADIN
subchondral plate is increased by a factor of four, for example, fracture of some trabeculae would be expected, since the ultimate strength would be exceeded. To consider in more detail what is likely to happen, the discussion will be subdivided, and the three main issues, buckling, stress concentration and fatigue failure, discussed separately. Buckling. If a sheet is loaded uniaxially in compression, it can fail by buckling. The criterion for buckling can be phrased best by the following question: When a system is slightly disturbed from its equilibrium position, does it return to equilibrium or tend to depart even further? The critical uniaxial compressive load (Roark, 1943) for buckling of a plate similar in configuration and support to the vertical trabeculae (planes 1-9) in the model is given by p = 4 n2Ebh3 1212
(19)
where b is the breadth of the plate, h the thickness, and 1 the length-O-06, O-003 and 0.10 in. respectively -for the plates in the model. Using equation (19) one calculates a buckling load of 0.8 lb for the idealized plate. The stress of 100 psi on the subchondral plate is distributed evenly over an area of 0.06 X 0.4 in. The total load on the structure is 2.4 lb. Since there are 9 supporting trabeculae, each must support roughly 0.267 lb. This is below the critical buckling load for an individual trabecula. However, the calculations (see Table 2) show that each trabecula does nor support the same load. The highest load supported is that at point F. The effective load at this point is 0.51 lb. The. minimum load in any element in this trabecula is 0.45 lb. Under dynamic conditions, buckling of this trabecula is to be considered a possibility. The bending moments (Fig. 7) in the trabecula are such that the trabecula is deformed into an ‘S’. This would tend to alter the buckling behavior somewhat from that predicted by equation (19).
The fixity factor (coefficient of (,rr2Ebh3)/122) in equation ( 19) varies typically from l/4 to a maximum of 4 for different configurations of the plates (Drucker, 1967). The calculations in this section for a fixity factor of 4 can thus be regarded as an upper limit for the critical load for any geometry encountered. Of course in the real material the sheets are stabilized laterally at irregular intervals, so buckling in vim under dynamic loading is probably not as likely as indicated by the results of these calculations; however local buckling of certain trabeculae appears to be unavoidable. Effect ofstress concentrators. The structure of trabecular bone is known to consist of sheets with holes in them (Dyson et al., 1970). The holes provide interconnecting passageways for marrow and blood vessels. Therefore, the presence of stress concentrators as above is to be the rule rather than the exception. Additional types of stress concentrators present in trabecular bone could be of importance. The canaliculi, lacunae, and the lamellar structures are all present in the sheets of bone. The canaliculi and lacunae are so small and localized in that they do not extend all the way through the sheet of bone, that they would probably be of less importance, Lamellae that end inside the structure could cause obvious problems at their points of termination, depending on the mechanical properties of the interfaces between the lamellae. These intercould provide effective faces, however, barriers to crack propagation (Piekarski, 1970). The effect of such stress concentrators on the membrane stresses and bending stresses can be best summarized by the following two equations which follow directly from the standard formulae, (Savin, 1961): a:,=
a,(l+2K)
u; = fqo ( 1+ 0*79K)
(20) (21)
where K = c/a, the ratio of semi-axes of an elliptical hole and u, and ob are the membrane and bending stresses respectively. AS an example of stress increases due to typical
MECHANICAL
BEHAVIOR
OF TRABECULAR
BONE
667
concentrators, for a membrane stress of -26OOpsi and K = 2, uI, = -13,000psi. For bending stresses of 5230 psi and K = 2, (r; = + 13,500 psi. Therefore on one surface of the trabecula (point E in Fig. 5) there exists a tensile stress of 500 psi, while on the other surface, a compressive stress of 26500 exists, which exceeds the compressive strength of the bone. Local fracture would therefore be expected to occur at the end of the major axis of the elliptical hole on the compressive side of the trabecula (the positive-Z side of the plane in Fig. 9). This fracture itself causes stress concentration, and, under cyclic loading, crack propagation would be expected. Accepted yield criteria (McClintock and Argon, 1966) can be used to investigate the question of generalized yielding in the model. One such criterion is the following: * =
[H((T,- ‘T,J2+ (u,-~J2+
(~u-~d21
+ 3(.r;.u+Tsz+?y~]1’2.
(22)
When v equals the yield strength of the material, yielding occurs. As typical example, point H has the following maximum stress components: (T,= 1133psi v,~ = 3473 psi CX.U= 1693 psi and @ = 4250 psi. which is much less than the yield strength. In view of the fact that the shear stresses are negligible in the vertical trabeculae, and the contribution from the above calculation only minor in the subchondral plate, the use of such yield criteria in an analysis of this model is probably not justified. A good indication of the yield behavior is thus gained through reference only to the normal stresses. Fatigue fracture. The fatigue life of compact bone subjected to zero minimum stress as in viva (as opposed to zero mean stress) has been been shown to vary from lo3 cycles at a maximum stress of 11,700 psi to lo7 cycles at a
Fig. 9. Position of local fracture around macroscopic stress concentrator in an inclined trabecula.
maximum stress of 8,750 psi (Freeman, Swanson and Day, 1971). Reference to the outerfiber normal stresses in Table 2 shows that regions of vertical trabeculae would have a dynamic fatigue life of less than lo3 cycles if the dynamic stresses are only twice as great as the static stresses. (lo3 cycles represents a brisk 20-minute walk at roughly 2 steps per sec.) During this relatively short fatigue life, the natural turnover rate of the bone probably has no effect in healing incipient fatigue cracks. If, in addition, mild (K = 1) stress concentrators such as macroscopic holes are present, the occurrence of fatigue fractures would be more probable. Notice that even higher normal stresses exist in the subchondral plate (points H-Q) than in the vertical trabeculae. Even under static conditions, certain regions of the subchondral plate will have a fatigue life of less
668
J. W. PUGH.
R. M. ROSE and E. L. RADIN
than lo3 cycles. Considering the relative absence of macroscopic stress concentrators in the subchondral plate as compared to the trabeculae, the incidence of fatigue fracture is probably roughly the same for all regions of subchondral bone. IMPLICATIONS FOR THE ROLE OF CANCELLOUS BONE IN JOINT DEGENERATION It has already been pointed out above that the stiffness at the surface of the subchondral plate varies considerably with position according to the arrangement of the underlying trabeculae, and that this variation leads to nonuniform compression of the cartilage under load. It is possibly significant in this connection that scanning electron microscopy studies (Walker et al., 1969) show that, with age, cartilage develops large wavelength undulations; although no examination of the underlying cancellous bone was apparently made in the above study, the spacing of the cartilage undulations happens to correspond approximately to the average intertrabecular spacing. It is important to note that the least compliant regions of the cancellous bone also absorb less shock. Thus the cartilage above these regions is compressed more and is exposed to more shock loading; it is in these regions, then, that joint degeneration may begin. CONCLUSIONS
A quantitative model having the elastic modulus of trabecular bone can be constructed by approximating the trabecular geometry with a hypothetical network of compact bone, and then using finite-element analysis with a high-speed computer. Because the number and thickness of trabeculae in the model must be somewhat less than they are in normal trabecular bone in order to have equivalent stiffness, there are two implications: first, that the intertrabecular soft tissue is probably of no importance to the mechanical properties, and second, that the compact bone which constitutes the trabeculae is probably not as stiff as cortical bone, al-
though we previously assumed that it was (see Introduction). The computations themselves lead to several conclusions: (a) The in vivo elastic deformation of trabecular bone must involve considerable bending and buckling of trabeculae. (b) The bending (or buckling) stresses on trabeculae can be comparable to or greater than the membrane stresses. (c) For the loads due to normal activity, fatigue fracture seems inevitable in some trabeculae, even in normal, healthy bone; for some trabecular orientations the fatigue life may be lo3 cycles or less. (d) For any geometry similar to that of the model, the stiffness is spatially dependent. The nature of the spatial dependence depends totally on the trabecular arrangement, and large variations in stiffness may occur over a single intertrabecular spacing. Such behavior is totally suppressed by accepted mechanical testing methods, but should be very manifest in a normal articular joint, in at least two ways: the cartilage will be more severely compressed in regions where the subchondral bone is stiffer, and the transmission of shock loads to the cartilage will be much greater in the same region. Thus the spatial inhomogeneity of trabecular bone may in fact play some role in the functioning (and/or degeneration) of articularjoints. Acknowledgement-This
work was sponsored
by the
Office of Naval Research. REFERENCES Burstein, A. H., Shaffer, B. W. and Frankel, V. H. (1970) Elastic analvsis of condvlar structures. ASME Publication 70-WAjBHF-I. Chen, P. E. (1971) Strength properties of discontinuous fiber composites. Polym. Engng Sci. 11,5 l-56. Crede, C. E. (1965) in Shock and Vibration Concepts in Engineering Design, p. 128. Prentice-Hall, New Jersey. Currey, J. D. (1970) The mechanical properties of bone. Clin. Orthop. 73,2 1O-22 1. In Dorland’s Illustrated Medical Dictionary, (1965) p. 1598, Saunders, Phila. Drucker, D. C. (1967) In Introduction to Mechanics of DeformableSolids. Chap. 16. McGraw-Hill, New York.
MECHANICAL
BEHAVIOR
Dyson, E. D., Jackson, C. K. and Whitehouse, W. J. (1970) Scanning electron microscope studies of human trabecular bone. Nature 225.957-959. Elias, H. (1971) Three-dimensional structure identified from single sections. Scient 174,993-1000. Evans, F. G. (1957) In Stress and Strain in Bones, p. 176-202. Thomas, Springfield, Illinois. Evans, F. G. (1961a) Relation of the physical properties of bone to fractures. ASOS Instruc. Coursefor Lecfurers l&110-121. Evans, F. G. (1961b) Biomechanical Studies of the Musculo-Skeletal System, Chap. 3. Thomas, Springfield, Illinois. Freeman, M. A. R., Day, W. H. and Swanson. S. A. V. ( 197 1) Fatigue fracture in the subchondral bone of the human cadaver femoral head. Med. Biol. Engng 9, to joint biomechanics.
Henry Ford Hosp. Med. Bull. 8,415-432.
Hackett, R. M. (1971) Viscoelastic stresses in a composite system. Polym. Engng Sci. 11,220-225. Kummer, B. (1966) Photoelastic studies on the functional structure of bone. Folio Biotheoretica 6.3 1. Logcher. R. D. et al. (1968) ICES STRUDL I1 Vol. IFrume Analysis. Research Report R68-91, Department of Civil Engineering, M.I.T. Cambridge, Mass. 02 139; (1971) ICES STRUDL II Vol. II-Additional Design and Analysis Facilities.
F. A. and Argon, A. S. (1966) Mechanical pp. 276-279. Addison-Wesley. Reading, Mass. Research Report R70-77, Department of Civil Engineering, M.I.T..Cambridge,Mass. 02139. McDonnell, EC1 ICES STRUDL User’s Manual for PBST Finite Element. Engineering Computer International, Inc., 1073 Massachusetts Ave., Cambridge, Mass. 02139. McElhaney, J.. Alem, N. and Roberts, V. (1970) A Porous Block Model for Cancellous Bone. ASME Publication 70-WA/BHF-2. McLeish. R. D. and Habboobi (1972) Strain gauge techniques for cadaveric bone. Engng Med. 1,37. Norris, C. H. and Wilbur, J. B. (1960) Elementary Structural Analysis, pp. 558-566. McGraw-Hill, New York. Piekarski, K. (1970) Fracture of bone. J. uppl. Phys. 41. McClintock,
of Materials,
215-223.
Piotrowski. G. and Wilcox, G. A.,Jr.( 1971)TheSTRESS program: A computer program for the analysis of stress in long bones. J. Biomechanics 4.497-506. Pugh, J. W.. Runkle. J. C., Rose, R. M., Paul, I. L. and Trabecular Bone. Ph.D. Thesis, M.I.T., Cambridge, Mass. 02139. Pugh, J. W., Runkle. J. C.. Rose, R. M., Paul, I.L. and Radin, E. L.. Determination of the elastic modulus of individual trabeculae. To be presented at 1973 Meeting of Orthopnedic
Research
Society.
Radin, E. L. and Paul, I. L. (1969) Failure of synovial fluid to cushion. Nuture 222,999-1000. Radin. E. L.. Paul, I. 1,. and Rose, R. M. (1972) Hypothesis: Role of mechanical factors in the pathogenesis of primary osteo-arthritis. Lancet 1.5 19-522. Roark. R. J. ( 1943) In Formulas for Stress and Strain, p. 793. McGraw-Hill. New York.
BONE
669
Savin, G. N. (196 1) Stress Concentration Around Holes. pp. 73.305,32 1. Pergamon Press, Oxford. Swanson, S. A. V. and Freeman, M. A. R. (1966) Is Bone Hydraulically Strengthened? Med. Biol. Engng 4, 433-438.
Walker, P. et al. (1969) Lubrication Mechanism in Human Joints: A Study by Scanning Electron Microscopy. In Lubrication and Wear in Joints (Edited by Wright, Verna), pp. 49-56. Lippincott, London. Yamada. H. (19701 Strenath of Bioloeical Materials. pp. 19-34. Williams &nd Wilkins: Bait&ore. Zienkiewicz, 0. C. ( 197 1) The Finite-Element Method in Engineering Science. McGraw-Hill, London.
;
619-629.
Frost, H. M. (1960) Introduction
Behavior
OF TRABECULAR
NOMENCLATURE semi-minor axis of ellipse (in.) breadth (in.) semi-major axis of ellipse (in.) contiguity of trabeculae elastic modulus or stiffness (psi) force (lb) thickness (in.) ratio of semi-axis of an ellipse length (in.) length (in.) moment couple (in. lb/in.) bending moment (in. lb) stress resultant (lb/in.) load (lb) stiffness coefficients (lb/in.) stress coefficients (Ib/in.3) displacements in x-, y-. z-direction (in.) ultimate compressive strength (psi) ultimate tensile strength (psi) volume fraction (in.3/in.3) Cartesian coordinates (in.) angle (degree or radians) nodal displacement (in.) Poisson’s ratio normal stress (psi) yield strength under combined stresses (psi) stress due to presence of stress concentrator
C’; E F h K
1 L M Mb N P S T uo. vo. w0 ucs UTS V, x, VTz a 6 lJ rr ;r (7’
(psi) bending stress (psi) membrane stress (psi) shear stress (psi) column matrix square matrix. APPENDlX STRUDL ‘SIMPLE PROBLEM’ ‘BONE ANALYSIS’ TYPE SPACE FRAME DEBUG ALL DUMP TIME JOINT COORDINATES 1 1. 0. 0. s 2 0. 0. 0. s 3 1. 0. I. 4 0. 0. 1. 5 0. 1. 1. 6 1. I. 1. 7 1. I. 0. s
670
J. W. PUGH,
R. M. ROSE and E. L. RADIN
8 0. 1. 0. S ELEMENT INCIDENCES 1321 2234 3354 4536 5687 6865 UNITS LBS INCHES ELEMENT PROPERTIES 1 TO 6 TYPE ‘PBST’ THICKNESS 0.003 CONSTANTS E 2000000.0 ALL POISSON 0.3 ALL LOADING 1 ‘TOP’ ELEMENT LOADS 3 TO 4 SURFACE FORCES GLOBAL PZ -100.0 GLOBAL SURFACE FORCES STIFFNESS ANALYSIS WITH REACTIONS NJP 3
LIST DISPLACEMENTS, FORCES, REACTIONS, FINISH
STRESSES, ALL
z
100 psi
t
I
STRAINS.