Simulation of Damage in Human Cortical Bone with Nonlinear Finite Element Analysis in MATLAB

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 cancel­ lous 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.


Introduction 2. Procedure
Accurate bone damage simulation can lead to better predic-2. 1 Making the Meshes tion of bone life, which can be used in prosthetic design and may First, the two-dimension finite element meshes were created aid clinical studies of disease such as osteoporosis. A number in the finite element analysis package Marc (MSC Software, of in vitro fatigue studies have been performed, establishing Santa Ana, CA, USA). There were four different meshes created: phenomenological models. One such model relates damage rate a tensile specimen, a bending specimen, a cancellous bone strut, per cycle, to the "normalized stress" (stress divided by initial and a slice of trabecular bone ( Figure 1). The tensile and bendmodulus) applied in cyclic loading. [1] ing specimens are modeled after fatigue tests specimens used by Finite element analysis is widespread in stress analysis of Zioupos. [3] The cancellous strut assumes symmetry and homo 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 FOR TRAN. [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 Figure 1. Finite element meshes of all four structures. From bone structures. Modeling damage with MATLAB can bring left to right, the tensile specimen (divided into a quarter of the study of bone damage to a larger population. the structure for symmetry), the bend test specimen (divided in half for symmetry), the cancellous strut, and the cancel lous slice. Green bars represent vertical axis constraints, orange bars represent horizontal constraints. Blue arrows represent the applied load.
geneity 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.

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.

Modeling the Damage Laws
The damage laws were coded in MATLAB functions and were applied at each cycle and repeated until the sample "broke", de fined 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 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 rou tines. The process was automated to repeat cyclically until the specimen failed.

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.

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 stud ies (strut and slice).

Modeling the Damage Laws
The damage function was first tested with the tensile mesh. The mesh was cyclically loaded with a constant maximum load.

D(n) =1− E
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 ele ments were adjusted whose major principle stress was negative, by multiplying the stress factor inserted into Equation (1) Figure 3. The error in damage levels are also plotted in Figure 4. As the number of cycles to failure decrease, the error increases. This is explained by the peak stresses in the shoulder of the sample seen in early cycles of Figure 2. With lower stresses and higher cycles to failure, these peak stresses and high damage areas smooth out and failure occurs in the gage length as the test is designed to do.

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.
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.
The strut was run with an apparent stress (net load divided by cross sectional area) of 7.39 MPa, which resulted in a maxi mum 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 compres sion. 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 (< 0.1%). This is consistent with earlier FE studies.

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.