Nonlinear Micro Finite Element Analysis of Human Trabecular Bone

A study

—by Harun H. Bayraktar, Technical Support Engineer, ABAQUS Inc.

Trabecular bone must withstand the loads that arise during daily activities as well as those from trauma. Investigation of the mechanical properties of trabecular bone presents a challenge because of its high porosity and complex architecture, both of which vary substantially between anatomic sites and across individuals. While Micro Finite Element (µFE) analysis of this type of bone is the most commonly used method to analyze the bone’s mechanical behavior, the large size of these models has forced researchers to use custom codes and linear analysis.


Background
Found at the end of long bones (e.g., femur) and in cuboidal bones (e.g., spine), trabecular bone is a major load bearing biological tissue in the human skeleton. Its mechanical properties are of great clinical and research interest. Improved understanding of trabecular bone mechanical properties will provide insight into fracture mechanisms in bones as well as allow the assessment of the effects of aging, disease, and drug treatments. With over 85% porosity in the spine, the highly porous bone has a complex architecture that varies substantially between anatomic sites and across individuals. Therefore, mechanical property data from multiple specimens are required to determine statistically the mechanical properties of trabecular bone.

Micro finite element models (µFE) are used extensively to study the mechanical properties of trabecular bone, both at the continuum level and at the microstructural level. These models are obtained through high-resolution imaging of the bone specimens that can be automatically converted into finite element meshes with hexahedral elements. All elements in these meshes are identical and are typically 50 µm in size. At this level of discretion, a 5 mm cubic specimen µFE model will typically have over half a million degrees of freedom. µFE models of bone specimens that are similar to experimentally used specimens (8 mm dia. and 15 mm length) contain millions of degrees of freedom.

In the past the large scale of these problems forced many researchers to use custom codes that utilize element-by-element iterative solvers. Because of the complexity of nonlinear finite element modeling, these custom codes have been limited to linear elastic analysis. Although linear elastic finite element models cannot simulate failure behavior of bone, they are widely used by researchers to determine bone tissue elastic properties by calibrating them against experimental data (Ref. 1).

However, many questions remain unanswered regarding nonlinear mechanical behavior of trabecular bone.

FEA methodology
ABAQUS/Standard was well suited to these types of analyses because it solves large problems using parallel execution and includes sophisticated material models. In this technology brief we investigated the role of geometrical nonlinearities in trabecular bone mechanical behavior using ABAQUS/Standard. We also demonstrated the feasibility of solving large problems by examining the parallel performance (i.e., scalability) of a single linear elastic analysis for a model containing over 4 million degrees of freedom.

One human vertebral trabecular bone specimen with 9% volume fraction was imaged using microcomputed tomography (µCT 20, Scanco Medical AG, Bassersdorf, Switzerland) at 22 µm resolution. Two µFE models were created. First, the entire cylindrical specimen was meshed with hexahedral elements of size 44 µm. Second, a 5 mm cubic sub-region from the core of the cylinder was used to create a model with the same element size. Mesh quantities for both models are listed in Table 1.

The cylinder model was used to assess the parallel performance of the direct sparse solver. Frictionless displacement boundary conditions were applied at the top and bottom surfaces to simulate a compressive 1% strain. Linear elastic analyses were conducted using 1, 2, and 4 CPUs of an HP rx8620 computer.

For nonlinear analysis, the 5 mm cube model was used. A cube of this dimension is large enough to determine continuum level properties but small enough to make nonlinear analysis feasible. Bone tissue was modeled using the cast iron plasticity material model. Cast iron plasticity provides elastic-plastic behavior with different yield strengths and hardening in tension and compression and results in an unsymmetrical element stiffness matrix.

Therefore, the parallel sparse direct solver with unsymmetric storage was used. A tissue elastic modulus of 13.4 GPa was used with a Poisson’s ratio of 0.3 (Ref. 2). For the cast iron plasticity model, tissue yield stresses were 55.2 MPa in tension and 110.6 MPa in compression, based on the yield strains reported for human femoral trabecular tissue (Ref. 3). In both tension and compression, a hardening slope equal to 5% of the elastic modulus was used. Frictionless displacement boundary conditions were used to apply a 2% nominal strain in tension and compression. At such a low level of nominal strain, self-contact of the bone microstructure need not be considered. In addition, each model was run with and without the effects of geometrically nonlinear deformations. In total, four nonlinear analyses were performed and continuum level yield strains were calculated for comparison. All the analyses of the cube were performed using 2 CPUs of an IBM Power4 computer.

Results and conclusions
Linear analysis of the cylinder model took under 16 minutes wall clock time on 4 CPUs and used less than 11 GB of memory (Table 2).

The parallel direct solver scaling results are also shown in Table 2; the speed-up factor is based on the solve time. Nonlinear analyses of the cube µFE model with geometrical nonlinearities took less than 7.4 hours wall clock time and required 4.1 GB memory. Each nonlinear analyses required approximately 100 linear equation solves emphasizing the importance of the solver scalability. The localization of initial yielding within the bone structure makes convergence of nonlinear analyses particularly challenging (Figure 1).

The apparent stress (applied force/cross-sectional area = 25 mm²) is plotted against apparent strain (change in specimen length/original specimen length) in Figure 4. Initial yield is defined as the point at which 0.2% offset is reached. Similar to reported experimental data (Ref. 4), yield strains were higher in compression than tension.

Although the tissue material was hardening, softening was observed at the apparent level when geometric nonlinearities were included (Figure 2). In addition, the yield strains were similar to experimentally measured values, particularly in compression (Ref. 4). These results show that the different yield behaviors of the trabecular tissue in tension and compression and geometric nonlinearities need to be incorporated in µFE models to model the continuum level yield behavior of trabecular bone accurately.

Acknowledgements
ABAQUS, Inc., would like to acknowledge Professor Tony M. Keaveny of the University of California, the bone specimen image data and finite element mesh.

Harun Bayraktar joined ABAQUS from UC Berkeley, where his research concentrated on the mechanics of bone tissue and the use of finite element methods to determine the elastic and plastic properties necessary to study the non-linear behavior of human bone.  At ABAQUS, Harun provides advanced technical support for nonlinear applications of ABAQUS/Standard and ABAQUS/Explicit.

References
1. van Rietbergen, B.; H. Weinans; R. Huiskes; A. Odgaard, “A New Method to Determine the Trabecular Bone Elastic Properties and Loading Using Micromechanical Finite Element Models,” Journal of Biomechanics, vol. 28, pp. 69–81, 1995.
2. Rho, J. Y.; T. Y. Tsui; G. M. Pharr, “Elastic Properties of Human Cortical and Trabecular Lamellar Bone Measured by Nanoindentation,” Biomaterials, vol. 18, pp. 1325–1330, 1997.
3. Bayraktar, H. H.; E. F. Morgan; G. L. Niebur; G. E. Morris; E. K. Wong; T. M. Keaveny, “Comparison of the Elastic and Yield Properties of Human Femoral Trabecular and Cortical Bone Tissue,” Journal of Biomechanics, vol. 37, pp. 27–35, 2004.
4. Morgan, E. F.; and T. M. Keaveny, “Dependence of Yield Strain of Human Trabecular Bone on Anatomic Site,” Journal of Biomechanics, vol. 34, pp. 569–577, 2001.
5. ABAQUS User’s Manual, Version 6.4, ABAQUS, Inc., 2003.

Abaqus Inc.,
www.rsleads.com/409df-141