Three-dimensional finite element modelling of bone: effects of element size J.H. Keyak* and H.B. Skinner+* *Rehabilitation Research and Development, Department of Veterans Affairs Medical Center, San Francisco and “Department of Orthopaedic Surgery, University of California, San Francisco, California, USA Received November 1991, accepted February 1992 ABSTRACT 7%s study quantifies the effectsof ebment s& on the stress/strain results ofjnite elnnent (FE) models of bone that are generated with a previously describedautomated method. Thti method uses cube-shapedhzxahedral elements, which enabledeltment shape and aspect ratio to be held constantwhile the effectiof elementse were studied. Threemodelsof a human proximal femur, each with a di$erent elementsize (3.7 mm, 3.8 mm and 4.8mm), were anulysed. Convqence in strain energyof the models had been vert$ed in previous work. lXe stressesand strains predictid by the models were compared on a pointwise basis using linear regressionanalyti. 7Iwre was a general decreasein the hovelof stress and strain when element s& was increased,even though convetgnce in strain ener~ had been achieved. An increasein element width from 3.1 mm to 3.8 mm &creased the predicted stressesby 13% to 29% overall; the predicted strains decreasedby 4% to 20% for the same increasein elementsize. Theseresults indicate thut linear cube-shapedhcxahedral elements must be very small (3 mm on a si& or smaller) to representth.eshnrp variations in mechanicalproperties thnt exist in bone,and that use of larger elements&creasesth predictcd stressesand strains. The elcmmts used in this study are similar to those typical@ used to representtrabeculurbone in conventional(non-automated)FE moaWing metho&. Therefore,the sensitivity of the stress/strain results to element se that was found fir trabeculurbone also applies to conventionalmodellingof such bone. This sensitivity to elements& implies that quantitative comparisonsof the stresses/ strains predicted in trabeculurboneby d@rent FE models may not be meaningfulif the elementsin those moa’elsare not the same s&. T?uzqualitative results of all three mo&h were in agreement, however, indicating that qualitative comparisonsof FE moo?&may be made. Keywords:
Finite element method, bone, computed tomography, biomechanics
INTRODUCTION Stress/strain anal sis using the finite element (FE) method requires x at the elements in the model be of an appropriate size. For typical engineering applications involving homogeneous materials, the adequacy of the FE mesh may be verified by confirming that convergence in strain energy has been achieved’. For applications involving inhomogeneous materials, the FE mesh must, in addition, adequately represent the variations in mechanical properties of the material. If the material is very inhomogeneous, as is bone24, satisfaction of the strain energy criterion alone may not ensure that the mechanical properties are modelled well enough to obtain meaningful stress/ strain data. Presently, FE models of bone are constructed without knowing the extent to which element size influences the re resentation roof mechanical perties and, there Pore, the stress/strain results. lX us, the reliability of these models is unknown. In preCorrespondence and reprint requests to: Dr Harry B. Skinner, Department of Orthopaedic Surgery, Room U-471, University of California, San Francisco, CA 94143-0728, USA 0 1992 Butterworth-Heinemann ol41-5425/92/06483-07
vious work, we developed an automated method of enerating FE models of bone using eight-noded elements of a user-specified Ylinear) cube-shaped size5. However, as is the case for conventional models, the effects of element size on the results of these models are unknown. The purpose of this study is to quantify the effects of element size on the stress/strain results of FE models of bone generated by this automated method. The cube-shaped hexahedral elements used by this method enable element sha e and as ect ratio to be held constant while the e fpects of e Pement size are studied. Because linear hexahedral elements are frequently used in three-dimensional FE models of bone&” the results of the present investigation have implicatibns for models generated by conventional methods.
METHODS Finite! element modelling The proximal portion of the left femur of a 4 l-yearold woman was modelled using the method described by Keyak et aL5. This method automatically generates
for BES J. Biomed. Eng. 1992, Vol. 14, November
483
Finite ebment modclling of bcme:J.H. Ktyak and H.B. Skinner
FE models of bone using eight-noded-cube-shaped elements of a constant user-specified size arranged in a cubic lattice. The elastic moduli of the elements are computed from a CT scan of the bone, as follows: The apparent bone density for each CT scan ixel within an element is estimated from a linear ca Pibration of the CT scan data. The elastic modulus for each pixel is then calculated from the estimated bone density using the density-modulus relation developed by Carter and Hayes’ . Finally, the eIastic modulus for each element is calculated by averaging the moduli for the pixels within the element. Use of this automated modellin technique enabled rapid generation of the FE mo f els. All elements in each model were the same size and shape, and roperties were inde envariations in mechanical dentl computed from tl! e same CT scan x ata, there isy permitting comparisons to be made without bias. Twenty-four contiguous transverse CT scan slices were obtained in vivo using a GE9800 Research Scanner operating at the following settings: 140 kVp; 70mA, 3 s; slice thickness, 3.Omm; pixel size, 1.38 mm; 320 X 320 matrix. Three models were generated. Model 1 had 1665 nodes and 1152 4.8 mm elements, Model 2 had 3034 nodes and 2232 3.8 mm elements, Model 3 had 4920 nodes and 3788 3.1 mm elements (f;igure 7). Convergence in strain energy for these models was verified in previous work5. The element moduli computed from the CT scan data ranged from 0 MPa to 15 308 MPa for Model 1, from 0 MPa to 16 739 MPa for Model 2, and from 0 MPa to 17 7 12 MPa for Model 3. There is an increase in maximum element modulus from Model 1 to Model 3 because smaller elements better represent peaks in the bone modulus. To make the model corn atible with the FE analysis software, the number of ePement moduli computed for each model was reduced to 50
1730 N 1670 N 310 N
1270 N
Figure 1 Model 2 consisted of 3034 nodes and 2232 3.8mm elements. The forces were applied at the nodes that are marked
484
J. Biomed. Eng. 1992, Vol. 14, November
values defined at equal intervals: Model 1, 153 MPa to 15 155 MPa; Model 2, 167MPa to 16571 MPa; Model 3, 177 MPa to 17 535 MPa. The Poisson’s ratio was assumed to be 0.4. The three-dimensional forces of gait during single limb su port were applied quasi-statically to each model tKrough a joint-reaction force at the femoral head and through the force of the abductors on the greater trochanter (Figure 7).The distal end of each model was restrained. The joint-reaction force, 1730 N (3.0 times body weight) directed through the centre of the femoral head, was derived from Case 2 of RydelP I’. This force was a plied to the surface of the models at nodes locate x within +_28” of the resultant force, measured from the centre of the femoral head. The nodal force components applied to a particular model were of equal magnitude, and each was directed through the centre of the femoral head. The resultant force imposed by the abductor muscles was estimated by requiring a quasi-static equilibrium force balance on the supporting limb. The derived force, 1270 N (2.2 times body weight), was applied in each of the models to the same portion of the superolateral surface of the greater trochanter. The force was equally divided over 17, 23, and 21 nodes in Models 1,2 and 3, respective1 . FE analyses of the models were performed on an I BM 3090 computer with SUPERB software (Structural Dynamics Research Corporation, Milford, OH), which uses isoparametric eight-noded hexahedral elements.
Evaluation of the finite element models To assess the effects of element size, the results of Model 3 were assumed to be ‘correct’ and were used as reference values: for each stress and strain component at each node in Model 3, the corresponding stress and strain corn onents predicted by Models 1 and 2 were found an B compared with those predicted by Model 3. Linear interpolation was used to evaluate the data of Models 1 and 2 between nodes. Because of differences in model geometries, some nodes of Model 3 were located outside the boundaries of Models 1 and/or 2 and could not be compared; 3843 locations in Model 1 and 4194 locations in Model 2 remained for examination. Linear regression analysis was used to evaluate the stress/strain comparisons for each component direction, with the data of Model 3 servin as the independent variable and the data of Mode Bs 1 and 2 servin as dependent variables. These regressions provi Bed measures of the effect of element size on the stresses/strains predicted over the entire model. In each analysis, the correlation coefficient indicates the extent of qualitative agreement between models, i.e., whether the models predict high stresses/strains in the same regions; deviation of the slope from 1 would indicate diminishing quantitative agreement. The standard error of the estimate (SEE) is a measure of the scatter of the data about the regression line and therefore reflects a component of uncertain in the stress/strain data. This uncertainty in the mo 2 elling is a direct consequence of defining elements with a particular size and specific boundaries. Initial plots of the results revealed that the regres-
Finite &runt modclling of bonc:JH. 6000
Table
4000
I
Slopes of regression
Stress/strain Component
Model l/Model Comparison
0.x UY
2000
*z TXY
TX< 7v
0
e, 4 6, YXY Yxz
-2000
YYT
-4000
+ -6000 -6000
I+
-4000
I -2000
I 0
1 2000
I
uaao
I 6000
sion lines were not representative of the data (Figwe 2); regression diagnostics indicated that snme data points exerted a disproportionate influence on the regression anaIysis’3. For most data points, the Cook’s distance (a measure of the influence of a particular point on a regression analysis) was below 0.0001; for the influential points, the Cook’s distance was much greater. To obtain regression equations that were more representative of the data, a second set of regression analyses that excluded the influential data points was performed. In each of these analyses, points for which the Cook’s distance was greater than 0.001 were omitted. This cut-off was selected because it was the highest value that would result in regression lines representative of the data. Fewer than 7.5% of the data points in each regression were omitted, with the actual number de ending on the stress/strain component and mode Ps that were examined. The locations of these points in the FE models were identified for future reference. The regression result%&iy indicate the extent of overall quantitative and qualitative agreement of the models. To indicate the quantitative differences in stress/strain data as afinction of location, the stresses/ strains in Model 3 were vectorial1 subtracted from those in Models 1 and 2 at each o r the locations that were compared previously. These differences in stress/strain components were then used to compute maximum normal differences in stress/strain within coronal sections, and these data were plotted.
RESULTS All regression equations were significant (Model 11 Model 3 comparison, R > 0.75, P< 0.001; Model 2/ Model 3 comparison, R > 0.88, P< 0.001). Tabh 7 and 2 summarize the results. Plots of two regressions are shown in Fi$mx 3a and b. Bonferrani t-tests
Ijnes* 3
Model 2/Model3
Comparison
Change
0.560 0.547 0.754 0.508 0.665 0.660
0.728 0.738 0.866 0.707 0.809 0.774
+23.1 +25.9 + 12.9 +28.1 + 17.8 + 14.7
0.824 0.901 0.906 0.611 0.802 0.780
0.914 0.963 0.951 0.795 0.919 0.850
+9.8 +6.4 +4.7 +23.1 + 12.7 +8.2
(%)’
*All slopes are significantly different from 1. +Percent change in slope is computed from:
(slope Model Z/Model
Table
2
SEE for regression
Stress/strain Component
Model 1/Model Comparison
Q, UY u
s. \
,
modclling of b0ncJ.H. Kq& at& H.B. skinm
- * ,
I
.
*
*
.,,.\.,,a..* . . ..\..I,..’ . . . . . . . . ..*.a . . . . . ..\..a . . . . . . . . , \. : . . . . . . \
.
.
.
.
.
.
.
.
,
i.:,
.
.
.
.
.
‘.
,
.
.
\
‘,,
:\ ::
...,
. . _ __,.i:
.
.
.
.
.
-
.
.
.
.
.
.
_.:;:..y.,
..
.
_ + ‘.::.’ . .
r-h.
.
.
.
.
-
. .
.
4
1
t i
\
t
.
.
:
:
1 ;
*
f
.
q
-
:
*
*
.
’
.
.
\
-
-
.
*
.
.
.
I
.
.
.
.
-
.
.:’
.
*
c
.
t
.
.
‘y;?_
Finite element method, bone, computed tomography, biomechanics
INTRODUCTION Stress/strain anal sis using the finite element (FE) method requires x at the elements in the model be of an appropriate size. For typical engineering applications involving homogeneous materials, the adequacy of the FE mesh may be verified by confirming that convergence in strain energy has been achieved’. For applications involving inhomogeneous materials, the FE mesh must, in addition, adequately represent the variations in mechanical properties of the material. If the material is very inhomogeneous, as is bone24, satisfaction of the strain energy criterion alone may not ensure that the mechanical properties are modelled well enough to obtain meaningful stress/ strain data. Presently, FE models of bone are constructed without knowing the extent to which element size influences the re resentation roof mechanical perties and, there Pore, the stress/strain results. lX us, the reliability of these models is unknown. In preCorrespondence and reprint requests to: Dr Harry B. Skinner, Department of Orthopaedic Surgery, Room U-471, University of California, San Francisco, CA 94143-0728, USA 0 1992 Butterworth-Heinemann ol41-5425/92/06483-07
vious work, we developed an automated method of enerating FE models of bone using eight-noded elements of a user-specified Ylinear) cube-shaped size5. However, as is the case for conventional models, the effects of element size on the results of these models are unknown. The purpose of this study is to quantify the effects of element size on the stress/strain results of FE models of bone generated by this automated method. The cube-shaped hexahedral elements used by this method enable element sha e and as ect ratio to be held constant while the e fpects of e Pement size are studied. Because linear hexahedral elements are frequently used in three-dimensional FE models of bone&” the results of the present investigation have implicatibns for models generated by conventional methods.
METHODS Finite! element modelling The proximal portion of the left femur of a 4 l-yearold woman was modelled using the method described by Keyak et aL5. This method automatically generates
for BES J. Biomed. Eng. 1992, Vol. 14, November
483
Finite ebment modclling of bcme:J.H. Ktyak and H.B. Skinner
FE models of bone using eight-noded-cube-shaped elements of a constant user-specified size arranged in a cubic lattice. The elastic moduli of the elements are computed from a CT scan of the bone, as follows: The apparent bone density for each CT scan ixel within an element is estimated from a linear ca Pibration of the CT scan data. The elastic modulus for each pixel is then calculated from the estimated bone density using the density-modulus relation developed by Carter and Hayes’ . Finally, the eIastic modulus for each element is calculated by averaging the moduli for the pixels within the element. Use of this automated modellin technique enabled rapid generation of the FE mo f els. All elements in each model were the same size and shape, and roperties were inde envariations in mechanical dentl computed from tl! e same CT scan x ata, there isy permitting comparisons to be made without bias. Twenty-four contiguous transverse CT scan slices were obtained in vivo using a GE9800 Research Scanner operating at the following settings: 140 kVp; 70mA, 3 s; slice thickness, 3.Omm; pixel size, 1.38 mm; 320 X 320 matrix. Three models were generated. Model 1 had 1665 nodes and 1152 4.8 mm elements, Model 2 had 3034 nodes and 2232 3.8 mm elements, Model 3 had 4920 nodes and 3788 3.1 mm elements (f;igure 7). Convergence in strain energy for these models was verified in previous work5. The element moduli computed from the CT scan data ranged from 0 MPa to 15 308 MPa for Model 1, from 0 MPa to 16 739 MPa for Model 2, and from 0 MPa to 17 7 12 MPa for Model 3. There is an increase in maximum element modulus from Model 1 to Model 3 because smaller elements better represent peaks in the bone modulus. To make the model corn atible with the FE analysis software, the number of ePement moduli computed for each model was reduced to 50
1730 N 1670 N 310 N
1270 N
Figure 1 Model 2 consisted of 3034 nodes and 2232 3.8mm elements. The forces were applied at the nodes that are marked
484
J. Biomed. Eng. 1992, Vol. 14, November
values defined at equal intervals: Model 1, 153 MPa to 15 155 MPa; Model 2, 167MPa to 16571 MPa; Model 3, 177 MPa to 17 535 MPa. The Poisson’s ratio was assumed to be 0.4. The three-dimensional forces of gait during single limb su port were applied quasi-statically to each model tKrough a joint-reaction force at the femoral head and through the force of the abductors on the greater trochanter (Figure 7).The distal end of each model was restrained. The joint-reaction force, 1730 N (3.0 times body weight) directed through the centre of the femoral head, was derived from Case 2 of RydelP I’. This force was a plied to the surface of the models at nodes locate x within +_28” of the resultant force, measured from the centre of the femoral head. The nodal force components applied to a particular model were of equal magnitude, and each was directed through the centre of the femoral head. The resultant force imposed by the abductor muscles was estimated by requiring a quasi-static equilibrium force balance on the supporting limb. The derived force, 1270 N (2.2 times body weight), was applied in each of the models to the same portion of the superolateral surface of the greater trochanter. The force was equally divided over 17, 23, and 21 nodes in Models 1,2 and 3, respective1 . FE analyses of the models were performed on an I BM 3090 computer with SUPERB software (Structural Dynamics Research Corporation, Milford, OH), which uses isoparametric eight-noded hexahedral elements.
Evaluation of the finite element models To assess the effects of element size, the results of Model 3 were assumed to be ‘correct’ and were used as reference values: for each stress and strain component at each node in Model 3, the corresponding stress and strain corn onents predicted by Models 1 and 2 were found an B compared with those predicted by Model 3. Linear interpolation was used to evaluate the data of Models 1 and 2 between nodes. Because of differences in model geometries, some nodes of Model 3 were located outside the boundaries of Models 1 and/or 2 and could not be compared; 3843 locations in Model 1 and 4194 locations in Model 2 remained for examination. Linear regression analysis was used to evaluate the stress/strain comparisons for each component direction, with the data of Model 3 servin as the independent variable and the data of Mode Bs 1 and 2 servin as dependent variables. These regressions provi Bed measures of the effect of element size on the stresses/strains predicted over the entire model. In each analysis, the correlation coefficient indicates the extent of qualitative agreement between models, i.e., whether the models predict high stresses/strains in the same regions; deviation of the slope from 1 would indicate diminishing quantitative agreement. The standard error of the estimate (SEE) is a measure of the scatter of the data about the regression line and therefore reflects a component of uncertain in the stress/strain data. This uncertainty in the mo 2 elling is a direct consequence of defining elements with a particular size and specific boundaries. Initial plots of the results revealed that the regres-
Finite &runt modclling of bonc:JH. 6000
Table
4000
I
Slopes of regression
Stress/strain Component
Model l/Model Comparison
0.x UY
2000
*z TXY
TX< 7v
0
e, 4 6, YXY Yxz
-2000
YYT
-4000
+ -6000 -6000
I+
-4000
I -2000
I 0
1 2000
I
uaao
I 6000
sion lines were not representative of the data (Figwe 2); regression diagnostics indicated that snme data points exerted a disproportionate influence on the regression anaIysis’3. For most data points, the Cook’s distance (a measure of the influence of a particular point on a regression analysis) was below 0.0001; for the influential points, the Cook’s distance was much greater. To obtain regression equations that were more representative of the data, a second set of regression analyses that excluded the influential data points was performed. In each of these analyses, points for which the Cook’s distance was greater than 0.001 were omitted. This cut-off was selected because it was the highest value that would result in regression lines representative of the data. Fewer than 7.5% of the data points in each regression were omitted, with the actual number de ending on the stress/strain component and mode Ps that were examined. The locations of these points in the FE models were identified for future reference. The regression result%&iy indicate the extent of overall quantitative and qualitative agreement of the models. To indicate the quantitative differences in stress/strain data as afinction of location, the stresses/ strains in Model 3 were vectorial1 subtracted from those in Models 1 and 2 at each o r the locations that were compared previously. These differences in stress/strain components were then used to compute maximum normal differences in stress/strain within coronal sections, and these data were plotted.
RESULTS All regression equations were significant (Model 11 Model 3 comparison, R > 0.75, P< 0.001; Model 2/ Model 3 comparison, R > 0.88, P< 0.001). Tabh 7 and 2 summarize the results. Plots of two regressions are shown in Fi$mx 3a and b. Bonferrani t-tests
Ijnes* 3
Model 2/Model3
Comparison
Change
0.560 0.547 0.754 0.508 0.665 0.660
0.728 0.738 0.866 0.707 0.809 0.774
+23.1 +25.9 + 12.9 +28.1 + 17.8 + 14.7
0.824 0.901 0.906 0.611 0.802 0.780
0.914 0.963 0.951 0.795 0.919 0.850
+9.8 +6.4 +4.7 +23.1 + 12.7 +8.2
(%)’
*All slopes are significantly different from 1. +Percent change in slope is computed from:
(slope Model Z/Model
Table
2
SEE for regression
Stress/strain Component
Model 1/Model Comparison
Q, UY u
s. \
,
modclling of b0ncJ.H. Kq& at& H.B. skinm
- * ,
I
.
*
*
.,,.\.,,a..* . . ..\..I,..’ . . . . . . . . ..*.a . . . . . ..\..a . . . . . . . . , \. : . . . . . . \
.
.
.
.
.
.
.
.
,
i.:,
.
.
.
.
.
‘.
,
.
.
\
‘,,
:\ ::
...,
. . _ __,.i:
.
.
.
.
.
-
.
.
.
.
.
.
_.:;:..y.,
..
.
_ + ‘.::.’ . .
r-h.
.
.
.
.
-
. .
.
4
1
t i
\
t
.
.
:
:
1 ;
*
f
.
q
-
:
*
*
.
’
.
.
\
-
-
.
*
.
.
.
I
.
.
.
.
-
.
.:’
.
*
c
.
t
.
.
‘y;?_
