Virginia Polytechnic Institute and State University, Department of Engineering Science and Mechanics
Blacksburg, VA 24061


This paper examines bone fatigue damage using Finite Element Analysis (FEA) to predict changes in modulus, damage and cycles to failure. The MATLAB-coded material model has been previously implemented in a full FEA package. It is based on published models of damage verified with experimental data. The MATLAB program inputs a finite element mesh and simulates the damage iteratively. Four meshes are examined: a machined tensile specimen of the bone, a machined bending specimen, a single strut of cancellous bone, and a 2-D slice of cancellous bone. This model is a fast testing method for ideas, allows easy comparison of different test geometries and material damage models, and may help to bring the study of bone damage to a broader audience.

Keywords: Bone mechanics, fatigue damage, finite element analysis

1. Introduction

Accurate bone damage simulation can lead to better prediction of bone life, which can be used in prosthetic design and may aid clinical studies of disease such as osteoporosis. A number of in vitro fatigue studies have been performed, establishing phenomenological models. One such model relates damage rate per cycle, to the "normalized stress" (stress divided by initial modulus) applied in cyclic loading.[1]

Finite element analysis is widespread in stress analysis of bone-implant systems. It allows complex shapes to be input to determine stresses, strains, and displacements. Earlier studies have incorporated damage models in fatigue FE studies of bone using MSC.Marc and user defined subroutines in FORTRAN.[2] This technique was powerful, but requires in-depth knowledge of expensive specialized software and limits the user to the corporate software designer. MATLAB is widely used at universities and elsewhere in the engineering community. The goal of this work is to implement the Cotton damage law into MATLAB and demonstrate its utility to simulate various 2-D bone structures. Modeling damage with MATLAB can bring the study of bone damage to a larger population.

2. Procedure

2.1 Making the Meshes

First, the two-dimension finite element meshes were created in the finite element analysis package Marc (MSC Software, Santa Ana, CA, USA). There were four different meshes created: a tensile specimen, a bending specimen, a cancellous bone strut, and a slice of trabecular bone (Figure 1). The tensile and bending specimens are modeled after fatigue tests specimens used by Zioupos.[3] The cancellous strut assumes symmetry and homogeneity of a regular hexagonal cellular solid, and was presented earlier by Cotton.[4] The cancellous slice has been previously modeled by Taylor.[2] The finished geometric mesh of nodes and elements were imported into MATLAB and converted into an input structure designed for the model. Material properties, nodal constraints, and loads were added to the input structure. For all models, Young's modulus of the bone tissue was taken as 15 GPa, Poisson's ratio of 0.35. Loads were applied such that peak stresses were on the order of 30 to 150 MPa.

2.2 Static Testing of the Meshes

All four meshes were simulated statically with a finite element analysis function modified from Zaicenco.[5] Our modifications included redesigning the data into a hierarchical structure. The FE function took the input structure and solved for the stresses and displacements in the mesh. The original graphical output was modified to better examine the stress tensors and nodal displacements. Results were tested and compared to calculated static results for verification for the tensile and beam meshes.

2.3 Modeling the Damage Laws

The damage laws were coded in MATLAB functions and were applied at each cycle and repeated until the sample "broke", defined as a stiffness reduction of 90%. The damage law[1] defined the damage rate, ΔD/ΔN, as a function of "normalized stress," σ/E*, by using a power law expression

where A and B where constants fit to experimental data, E* was the initial modulus and σ was the major (maximum magnitude) principle stress at the node. The use of normalized stress is widespread in bone mechanics, and allows for accurate extrapolation of material laws to bone of greatly different types. The damage, D(n), at cycle n, is defined as the fractional loss of stiffness

where E(n) is the modulus. The constants A and B were 1035.5 and 17, respectively. As these tests were in tension, values needed to be adjusted to represent compression. Compressive tests of bone show fatigue failure strengths and damage rates of roughly 50% of the respectively tensile results.[6] The elements were adjusted whose major principle stress was negative, by multiplying the stress factor inserted into Equation (1) by 0.69.

This model was run in MATLAB through a number of loading cycles. For each cycle, the FEA was used to determine stresses then the values were inserted into the damage law routines. The process was automated to repeat cyclically until the specimen failed.

2.4 Simulations Run

The tensile specimen was run through a simulation from zero to maximum tension (0-T) fatigue tests of different load magnitudes to verify accuracy. Stress versus cycles to failure (S-N) curves were generated and compared to published data. Bending, strut and cancellous samples were run to observe the effects of fatigue on these structures.

3. Results and Discussion

3.1 Static Testing of the Meshes

The four meshes created are presented, along with loads and constraints, in Figure 1. In a static analysis they all showed stresses and displacements predicted by analytic calculation (for tensile and beam), or that were consistent with earlier FE studies (strut and slice).

3.3 Other Models

Other models were also used to simulate fatigue damage in the beam, the bending specimen, strut and cancellous slice models. All models showed localized damage conforming to the location of maximum stress.

Error of the cycles to failure and damage Figure 3. Error of the cycles to failure and damage Error of the apparent strain for the tensile model Figure 4. Error of the apparent strain for the tensile model showing agreement with calculated results with lower applied normalized stress.

The beam was run with an applied load of 7.8 N, which resulted in a maximum initial strain of 5800 με in compression. The sample failed on the lower tensile surface at 2.25 million cycles. This is a much higher life than predicted by the tensile laws, due to the sample's outer surface damaging, dropping peak stresses there and raising stresses into the depth. These results are consistent with bending fatigue tests of bone, in that they fail in tension despite the additional compressive load from the load on top surface. It is worth noting that choosing a different definition of failure results in grossly different definitions of life, for example the 20% loss of stiffness used elsewhere would lead to a simulated life of 180,000 cycles. Damage distribution in the beam at 150,000 cycles is presented in Figure 5.

Damage in the beam specimen during the fatigue test Figure 5. Damage in the beam specimen during the fatigue test. The bottom tensile surface shows greater damage than the top compressive surface. Damage is highest midspan (right side of FE model) where the bending moment is highest. Damage in the cancellous strut during the fatigue test Figure 6. Damage in the cancellous strut during the fatigue test. Damage is highest near the arrows where the compression from the vertical section adds to the compression from bending in the horizontal trabeculae. Tension on the bottom surface is small enough not to damage noticeably

The strut was run with an apparent stress (net load divided by cross sectional area) of 7.39 MPa, which resulted in a maximum initial strain of 7441 με in compression. The sample failed in compression at 29,000 cycles. These results differed from the work of the previous simulations where tensile failure in bending occurred.[3] This is because this simulation constrained the strut model on both horizontal limits. Damage distribution in the strut during the test is presented in Figure 6.

The cancellous slice was run with an apparent stress (net force divided by whole cross sectional area) of 6.3 MPa, which resulted in a maximum initial strain of 11800 με in compression. Damage distribution in the slice at 379 cycles, is presented in Figure 7. It can be seen that only a few select areas experience any damage. Although these elements may have damage of 45% at this point, the entire structure has shown negligible loss of stiffness (

4. Conclusions

Bone damage simulation was performed in MATLAB, a widely used software in the science community. This simulation is more accessible and adaptable than typical FEA packages. The model is accurate with its predictions of the damage, modulus change, and cycles to failure. In the future, other damage laws[7, 8] could be added to the model and many different specimens configurations could be tested. This tool will readily allow both cross-comparison of different damage laws and interpretation of differences in test geometries and configurations, helping to broaden our understanding of bone damage.


John Cotton, Assistant Professor in Engineering Science and Mechanics at Virginia Tech, advised Mr. Kemeny with this project as an independent study and provided various modifications to the downloaded FE code. Frances Davis, a senior undergraduate in ESM at Virginia Tech, also worked on the FE code for a separate project and provided useful modifications of several of the graphics routines used.


[1] Cotton, J. R.; Winwood, K.; Zioupos, P.; Taylor M. Damage Rate is a Predictor of Fatigue Life and Creep Strain Rate in Tensile Fatigue of Human Cortical Bone Samples, Journal of Biomechanics, 2005, 127(2), 213-9.

[2] Taylor, M.; Cotton, J.; Zioupos, P. Finite element simulation of the fatigue behaviour of cortical and cancellous bone, Meccanica, 2002, 37, 419-429.

[3] Zioupos, P.; Wang, X.T.; Currey, J.D. The Accumulation of Fatigue Microdamage in Human Cortical Bone of Two Different Ages in Vitro, Clinical Biomechanics, 1999, 11(7), 365-375.

[4] Cotton, J. R.; Winwood, K.; Zioupos, P.; Taylor, M. Simulated material property degradation during fatigue of human cancellous bone, as modelled by a regular honeycomb structure, ASME Summer Bioengineering Conference, Snowbird, Utah, June 27-July 1, 2001.

[5] Zaicenco, A. "FEA for solid mechanics with MATLAB", MATLAB code downloaded from Matlab Central,, Posted 12/15/2004.

[6] Winwood K.; Zioupos, P.; Currey, J. D.; Cotton, J.R.; Taylor, M. Development rates of plastic and elastic components of strain in human cortical bone during tensile, compressive and shear fatigue loading and implications to bone biomechanics, Journal of Biomedical Materials Research, in press.

[7] Pattin, C.A.; Caler, W.E.; Carter, D.R. "Cyclic Mechanical Property Degradation During Fatigue Loading of Cortical Bone. Journal of Biomechanics, 1996, 29(1), 69-79.

[8] Griffin, L.V., et al., Model of Flexural Fatigue Damage Accumulation for Cortical Bone. Journal of Orthopedic Research., 1997, 15, 607-614.

About the Author

Steven Kemeny
 Steven Kemeny is from West Chester, Pennsylvania and went to Unionville High School. He is now a senior in Engineering Science and Mechanics at Virginia Tech. He is focusing his studies in biomechanics and mathematics. Steve wants to have a career in biomedical engineering. He is in the Marching Virginians, one of the marching bands at Virginia Tech. He plays the trumpet and has been a rank leader for two years.