Abstract
The objective of this study was to more fully understand the mechanical behavior of bone tissue that is important to find an alternative material to be used as an implant and to develop an accurate model to predict the fracture of the bone. Predicting and preventing bone failure is an important area in orthopaedics. In this paper, the macrodamage accumulation models in the bone tissue have been investigated. Phenomenological models for bone damage have been discussed in detail. In addition, 3D finite element model of the femur prepared from imaging data with both cortical and trabecular structures is delineated using MIMICS and ANSYS® and simulated as a composite structure. The damage accumulation occurring during cyclic loading was analyzed for fatigue scenario. We found that the damage accumulates sooner in the multiaxial than in the uniaxial loading condition for the same number of cycles, and the failure starts in the cortical bone. The damage accumulation behavior seems to follow a three-stage growth: a primary phase, a secondary phase of damage growth marked by linear damage growth, and a tertiary phase that leads to failure. Finally, the stiffness of the composite bone comprising the cortical and trabecular bone was significantly different as expected.
1. Introduction
In order to understand the bone fracture, it is very important to study the macrodamage of the bone with respect to mechanical and physiological loads. Bone tissue is a complex, multiphasic, heterogeneous, anisotropic, and highly hierarchized material structure. Predicting and preventing bone fracture is a very important area in orthopaedics given the volume of fractures that occurs annually. From a macroscopic point of view, bone tissue is divided into two types: the trabecular bone with 50–95% porosity [1] and the cortical bone with 5–10% porosity [1]. Bone tissue can be divided into five levels [2], which are macro, meso, micros, submicro, and nanostructure. The macrostructure level is the whole bone, which ranges from several millimeters to several centimeters, as shown in Figure 1. In this paper, an attempt has been made to establish a detailed understanding of the bone tissue mechanical behavior as it is important in the device design and to derive implant life. Correspondingly, an accurate damage prediction model for a bone tissue is needed in order to predict the fracture of the bone or the reliability of a bone-implant structure.
Numerous damage models were proposed using the macrostructure of the bone. However, each model has made an assumption regarding the mechanical properties, loading conditions, or the structure of the bone. These assumptions have not given realistic predictions for the damage accumulation in a bone. Depending on the mechanical properties of bone tissue, bone damage models can be divided into elastic-viscoplastic, elastoplastic, and plastic damage models. In addition, depending on the damage type, bone damage models can be divided into electromagnetic, fracture, bending, and fatigue damage models. The elastic-viscoplastic damage models take into consideration that the bone has elastic, plastic, and viscus material properties.
Recently, several models have been proposed that describe the damage model of the bone as an elastic viscoplastic model such as Keyak and Rossi [3]. They proposed fracture load by using finite element models and several failure theories [3]. However, they used isotropic material properties for bone tissue. Some studies proposed elastoplastic damage modes as well. These models take into account elastic and plastic material properties such as in the Garcia et al. study [4] and the Fondrk et al. study [5]. They proposed elastic plastic damage models for bone tissue and developed a model for cortical bone tissue only. Other studies proposed plastic damage models, which take into consideration that the bone has plastic material properties only. In addition to the mechanical properties, the loading conditions have a significant effect on the macrodamage accumulation of the bone. Some studies analyzed only tension, compression, or three-point bending [6].
Zlámal et al. proposed a numerical model for trabecular tissue using compression test and time-lapse X-ray radiography and three-point bending test of single trabecula [6]. Besides all of that, the main challenge that has been faced was to design a model for the bone that contains together the cortical and trabecular components of the bone. Some studies have worked only on the cortical component, such as Natali et al. [7]. Figures 2 and 3 show the use of a small sample from the femur to perform the finite element simulations. Other studies assumed that the damage starts at the trabecular components, so they created the damage models for the trabecular bone only, such as Charlebois et al. [10], Hambli [11], and Hosseini et al. [12]. Figures 4 and 5 show the use of a micro-CT to create small samples to perform the finite element simulations. In this paper, an attempt has been made to create a 3D model of the femoral bone that considers the anisotropic material properties of bone tissue and loads from realistic gait cycle to understand how damage accumulates in human bone tissue.
2. Material and Methods
2.1. Finite Element Modeling
Because of the difficulty in studying the macrodamage accumulation of the bone in vivo, mathematical and phenomenological models were used to simulate physiological conditions. A three-dimensional model of the femoral bone was created. Hip fractures are currently treated by trauma instrumentation. The choice of the biomaterial constituting the prosthesis determines the reliability. Hence, failure predictions in bone and bone-implant stability must be thoroughly investigated on computational models.
2.1.1. Creating the Model
A femur bone model was developed in three steps. Firstly, CT images for the femur were taken from a normal healthy femoral bone. Secondly, the CT images have been imported into the MIMICS 13.0 program to create a 3D model of the femoral bone, as shown in Figure 6. The cortical and trabecular components were created depending on the difference in density between them. Thirdly, the final model has been imported into SolidWorks (Dassault Systèmes SolidWorks Corp., Concord, MA, USA) to make the final improvements. The final model of the femoral bone that has both the cortical and the trabecular parts is shown in Figure 7.
2.1.2. Material Definition
The proposed material properties of the bone consider the anisotropic and nonhomogeneity of the bone with its two types, the cortical and the trabecular. The trabecular bone is a spongy region; its density is lower than that of the cortical region, which is the hard and dense part of the bone. There are various procedures that have been performed to approximate the modulus of elasticity (E) of the bone depending on Hounsfield units (HU) and density () [14–16]. To give a realistic approximation for the bone tissue material properties, nine elastic constants must be provided depending on the orientation of the principal axes of orthotropy. While it is straightforward to assign the principal axes to the cortical zone, it is very challenging for the trabecular zone. In this study, both the cortical and trabecular zones have been divided into eight smaller segments. Then, each segment has been divided into ten material groups.
Within the MIMICS program, the Hounsfield units (HU) were used to calculate the density () across each segment, and then the young’s modulus (E) has been calculated in the radial, axial, and circumferential directions. Figure 8 shows the HU distribution across the femoral bone CT images.
The mathematical relationship between Hounsfield units (HU) and effective density () that has been applied in MIMICS is as follows: where the unit for the effective density () is g/cm3. The CT slices were used to align the orientation of each segment material. Also, the orthotropic relationships between the elastic constants and the density are different for the cortical and trabecular parts as described earlier, as shown in Table 1. Also, Figure 9 shows the procedure that has been used to find the material properties for each part of the femoral bone that has been used in this study. Table 2 shows the 80 material groups with their densities and the nine elastic constants. The colors have been modified to be green-blue colors for the trabecular material groups and yellow-red colors for the cortical material groups, so one can differentiate between them, as shown in the last step. Finally, the material properties have been imported into ANSYS Workbench 16.2 (ANSYS Inc., Canonsburg, PA, USA). This procedure has been discussed in detail in the literature [14–16]. Additionally, to validate the importance of studying the bone as a composite material, the finite element simulation for each part of the bone has been done separately also. These simulations are extremely important in order to understand the effect of each zone on the whole bone.
2.1.3. Finite Element Mesh
Tetrahedral element type was used in this study with a minimum element size of 0.02 mm, as shown in Figure 10. The final model, which has been modified by SolidWorks, was imported into ANSYS Workbench 16.2. The finite element mesh was adapted automatically through the program. Following mesh convergence checks, the total number of elements was 26,898 for the whole femur, 22,328 elements in the cortical bone, and 4570 in the trabecular bone. In order to achieve repetitions of results within five percent, the meshing was refined with small increment size.
2.1.4. Loads and Boundary Conditions
To mimic physiological loading during normal walking, the reconstructed gait loads in the model were applied as a time-dependent analysis along its longitudinal axis. The gait cycle for walking was imported from the HIP98® program. In this program, total hip replacement joints on different patients were studied, and their movements were compared with the normal movements during different activities [17]. For the uniaxial loading, the equivalent maximum stress from the gait cycle was converted into a single load cycle. For the multiaxial loading condition from the gait pattern during walking, the initial applied triaxial load (Fz = vertical direction force, Fy and Fx = anterior–posterior and medial–lateral forces, resp.). A model fully fixed at the distal end was used in this study. The body weight acted on the femoral head and muscular force acted on the proximal femur. The hip contact, which transfers load from the upper body to the lower limbs, was investigated under static and dynamic conditions. The dynamic loads of the hip contact are shown in Figure 11 for walking condition. In this study, 106 numbers of cycles have been used assuming that the average number of human walking cycles in one year is 1,000,000/year.
2.2. Phenomenological Bone Macrodamage Model
Goswami investigated phenomenological life prediction methods in great detail [18–21]. The macroscopic deformable bodies can be described via continuum mechanics. The main assumption made considers the nonhomogeneous anisotropic material properties of the bone tissue for both the cortical and trabecular bones. Since we assumed our model to be a composite material with different Young’s moduli in the cortical and trabecular zones, an assumption was made that the strain in the cortical and trabecular zones is the same. Thus, we invoked a strain-based concept in damage modeling. The main material properties that are considered in creating the macrodamage accumulation model of the bone tissue are modulus of elasticity, fatigue strength coefficient, fatigue ductility coefficient, fatigue strength exponent, and fatigue ductility exponent. The first assumption in creating the model is to consider the elastic and plastic components as follows: where represents the total strain, the elastic strain, and the plastic strain. According to Hook’s law, where E denotes the modulus of elasticity and the total stress. As it is important to take the nonhomogeneity of the bone tissue, the total strain can be expressed as follows: where represents the strain into different directions.
According to the Coffin-Manson relation for the strain-life curve that is shown in Figure 12, the elastic and plastic parts can be expressed as follows: where is the fatigue strength coefficient, is the fatigue ductility coefficient, is the fatigue strength exponent, is the fatigue ductility exponent, and Cijkl are the elasticity tensor components. There are three types of fluctuating stresses. These are fully reversed, repeated, and fluctuating stresses. To create the model, the mean stress is taken into consideration. The mean stress exists when the loading is of a repeating or fluctuating type.
When considering the mean stress, the equation of the total strain is written as follows:
To find the mean stress and the ultimate stress , where and are the maximum and minimum von Mises stresses, respectively, during the loading cycle. As finding the empirical constants b and c needs experimental work, the universal slops method, shown in Figure 13, was used instead of the Coffin-Manson relation [23]. With the universal slope method, the fatigue strength exponent (b) is related to the ultimate tensile strength and ductility exponent (c) which is related to the true strain at the fracture of the material are replaced by average slope values of −0.12 and −0.6, respectively. The total strain relation is written as follows:
And by simplifying (8),
The amount of damage experienced by the body is quantified by a single damage variable D. The damage variable D = 0 when the material is undamaged, while the damage variable D = 1 when the material totally failed. According to Miner’s rule, the damage equation is as follows: where ni is the number of cycles of the occurred stress range and Nfi is the number of cycles to failure. In the case of anisotropic damage, the relation among the damage variable , stress, and strain is as follows: where D denotes the damage variable, σij are the stress components, εkl are the strains, and Cijkl are the elasticity tensor components (stiffness matrix), where where E denotes the Young’s modulus, G denotes the shear modulus, and denotes Poisson’s ratio. The superscript numbers denote the following: 1 for radial direction, 2 for circumferential direction, and 3 for longitudinal direction.
Also,
Table 3 shows the elasticity tensor components for each material group calculated by using (13) and (14). The proposed model of macrodamage accumulation of the bone tissue can be written as follows:
And by applying (15), the final equation for damage is as follows:
A single scalar damage variable is often insufficient to describe the variation in mechanical properties of damaged materials.
2.2.1. Applicability of Damage Models to the Femur
The gait cycle of the hip is used to predict the macrodamage accumulation for the femoral bone. Because the femoral bone is subjected to a complex loading, the rainflow method is used to simplify the counting of load cycles. This method is very accommodating with the use of Miner’s rule. The values of strength and ductility coefficients were used from the literature. The value of fatigue strength coefficient that was used is 6, and the fatigue ductility coefficient value that was used is 0.352 [24]. The procedure of using the rainflow method is shown in Figure 14.
A comparison between the proposed model and the three different macrodamage accumulations models was performed. The first model was for the cortical bone only, the second model for the trabecular bone only, and the third model for both cortical and trabecular composite bones. The first model is for the damage of the cortical bone from Pattin et al. [25]. In their study, thirty-two specimens of the cortical bone were used; the stress range (σ) = 83 MPa, number of cycles to failure (Nf) = 417, and the modulus (Ef) = 9.02 GPa. The other model is for the trabecular bone from Hambli [13]. In his study, five specimens were taken from the trabecular part of the head of the femoral bone; the stress range (σ) = 85 MPa, number of cycles to failure (Nf) = 107, and the modulus (Ef) = 0.17 GPa. The third model is for the damage of both the cortical and trabecular bone components from Zioupos and Casinos [26]. On the other hand, Miner’s rule and the finite element analysis data were used for the proposed model of the femoral bone that has both the cortical and trabecular components.
Figure 15 shows the relation between the damages of the bones in terms of cycle fraction (n/Nf) for the models, where n is the number of cycles at a specific stress range and Nf is the number of cycles to failure at the same stress range. The convex curve shows the damage of the cortical bone, while the concave curve shows the damage accumulation of the trabecular bone as the cycles increase. Monte Carlo simulation was performed using the results from deterministic analysis that shows damage accumulation with number of cycles, probabilistically. Monte Carlo simulation generates a set of random variables normally distributed about a mean and standard deviation. Monte Carlo simulation was carried out for the proposed model and the other three models. The mean and standard deviation for each macrodamage accumulation model have been measured by using the JMP program, as shown in Figure 16. The simulation for each model consisted of 200 random generated variables normally distributed. The probability of failure was calculated for each model. Table 4 shows the mean, standard deviation (SD), variance, and probability of failure for the four models.
3. Results
First, stress-strain analyses in uniaxial and multiaxial loading conditions are considered, then fatigue life prediction of the bone is carried out. The maximum von Mises stresses were obtained from both uniaxial and multiaxial loading conditions for static simulations, as shown in Figures 17 and 18, where the stresses are 78.7 and 99.4 MPa for the uniaxial and multiaxial loadings, respectively. Figure 19 shows the von Mises stresses for the dynamic simulation of both loading conditions, where the stresses are 105.8 and 124.2 MPa for the uniaxial and multiaxial loadings, respectively. In addition, the total life was obtained from both uniaxial and multiaxial loading conditions for the dynamic simulation assuming that the bone is not a synthetic material with regeneration/remodeling capabilities, as shown in Figure 20. The relation between the maximum stress and the number of cycles to failure is shown in Figure 21 for both the uniaxial and multiaxial loading conditions. The polynomial curve fitting (, for the multiaxial loading condition, and , for the uniaxial loading condition) proves that the stress decreases linearly with the increase in life or number of cycles to failure (Nf). In addition, Figure 22 shows that for the given life, the trabecular bone accumulated approximately 25% more plastic strain than the cortical bone. Also, the same trend was observed with elastic strain accumulation in the trabecular bone where it was approximately 6% higher than the cortical bone.
(a)
(b)
The finite element modeling of damage considers that the damage equals to zero when the element in the region of interest is undamaged. While, the damage is equal to one when the element failed. Figure 23 shows that the damage starts at the femoral neck after 106 cycles. To make a comparison between the cortical and trabecular components of the bone, each part has been evaluated individually. Figure 24 shows the relation between the damage and the fraction of fatigue lifetime (n/Nf).
(a)
(b)
(c)
Force versus displacement curves were presented from the finite element analysis, as shown in Figure 25. The polynomial curves fitting for the whole bone data (, ), for the cortical bone (, ), and for the trabecular bone (, ) suggest linear relation between the force and the displacement.
To measure the stiffness, data was generated from 26,898 elements and analyzed. It appears that the mean stiffness of the cortical bone was 7890; trabecular bone, 2860; and that of the whole bone, 4864 N/mm, as expected.
4. Discussion
Macrodamage accumulation of bone tissue was estimated for two different loading conditions. For the mathematical part, the elastic and plastic behaviors of the bone were taken into consideration. Also, the anisotropic and the nonhomogeneous material properties of the cortical and trabecular zones were included. The MIMICS program was used to create the material properties depending on the Hounsfield unit and the relation among the density of the bone and the modulus of elasticity and Poisson’s ratio assigned based on grayscale distribution across the 3D model of the femur.
To validate the importance of studying the bone as a composite material, a study on each part of the bone has been done separately also. This is very important in order to understand the effect of each zone on the whole bone. A comparison between the proposed model and the three different macrodamage accumulations models was performed. The first model [25] shows cortical bone behavior, and the study was done on a small sample of the femoral bone. The second model [13] was for the trabecular bone only, on a small sample of the bone. The third model [26] was for a portion of the bone that contained both cortical and trabecular parts. However, the material properties were simple, isotropic, and homogeneous for all the three models. Moreover, the macrodamage models were nonlinear in the first two models and linear in the last model. In the current study, Miner’s rule was used with the proposed femoral model that contained both the cortical and the trabecular components, and a linear relationship was assumed. Also, the rainflow method was used to simplify the gait cycle of normal walking. Figure 15 shows that the damage in the cortical bone is higher than that in the trabecular bone for the same fraction of fatigue cycles (n/Nf) at a particular stress range. For the cortical bone, the damage starts to decrease when (n/Nf) reaches 0.9, while for the trabecular bone, the damage keeps increasing till (n/Nf) reaches 1.
The probability of failure was calculated from the distribution of the random variables for each model by using Monte Carlo simulation. The probability of failure for the proposed model was 13.26%, while the probability of failure was 37.90% for the whole bone model and 42.20% for the cortical bone model. The reason for this large difference between the probability of failure of the proposed model and the other models is likely due to the entire femoral bone was studied in our study. On the other hand, the other models were on a small sample of the bone. The data clearly shows that the composite bone as considered in the present study has lower von Mises stresses and thus lower failure probability than elastic/plastic materials.
Furthermore, the finite element analysis allowed a deeper understanding for the macrodamage accumulation of bone tissue. A comparison between different loading conditions was evaluated. The first loading condition was a multiaxial loading, where the cycle for normal walking was used including Fx, Fy, and Fz; the other loading condition was the uniaxial loading, where the equivalent maximum stress from the gait cycle was converted into a single load cycle. The results showed a significant difference between the two loading conditions. In static finite element simulation, the maximum von Mises stress was 78 MPa for the uniaxial loading condition and 99 MPa for the multiaxial loading condition, respectively. These results were expected as the loads are higher in the multiaxial loading condition, which led to a greater amount of stress than those in the uniaxial loading condition. The advantage of the static simulation in this study is to confirm the validation of the 3D model of the femoral bone with the literature. In the dynamic finite element simulation, the maximum von Mises stresses were 105.8 MPa for the uniaxial loading condition and 124.2 MPa for the multiaxial loading condition, respectively.
The study showed that the failure starts faster in the multiaxial loading condition than that in the uniaxial loading condition for the same number of cycles. Furthermore, the finite element simulation showed that the relation between the stress and the strain stays the same till the stress reaches 65 MPa. Then, the stress starts to be higher for the multiaxial loading condition than that for the uniaxial loading condition for the same amount of displacement. In addition, the finite element simulation for the damage of the bone showed that the damage starts at the femoral neck. This result was expected, as the femoral neck is the weakest point in the femoral bone, and the study was done on a healthy bone that does not have any injury.
The anisotropic material properties were used in the finite element simulation of the proposed model. The damage accumulation process in a long bone may be described by a three-stage process, as shown in Figure 26.
Since stage II shows a linear behavior, stage I is reflective of the primary phase, where the damage developed in the cortical bone decreases as the cycle fractions increase. However, as stage I transitions to stage II, the damage accumulated in the cortical bone increases linearly until about a cycle fraction of 0.8; upon attaining this level of fatigue life, the damage mode transitions to a more rapid damage accumulation that cannot be described by a linear equation. This state, stage III, is known as the tertiary damage accumulation stage and must lead to the bone fracture. We proposed this behavior for the cortical bone, and the lower ranges of damage hold good for the trabecular bone as well, assuming that the bone is anisotropic and nonhomogeneous. However, the damage accumulated on the composite bone was derived from the material properties of both the cortical and trabecular bones. Three damage prediction equations were developed, as shown in Table 5, where B represents the fatigue cycle fractions (n/Nf) at a particular stress range. These equations can be used in deriving the bone fracture at a given stress range and fatigue life. The charts in Figure 26 show that R2 decreases as the damage increases. The failure starts in the cortical bone before the trabecular bone.
By comparing between the behaviors of the damage of the bone that were reported in Figure 15 versus Figure 27, our effort shows a very clear three-stage process. Therefore, the mathematical significance of our analysis is applicable in the engineering design.
The results from the FE analysis was used to determine mean stiffness. It appears that the mean stiffness of the cortical bone was 7890; trabecular bone, 2860; and that of whole bone, 4864 N/mm. Data generated from 26,898 elements was analyzed, and we observed a significant difference in the stiffness of each element. The stiffness is observed in Figures 25 and 27 for the whole bone and the cortical and trabecular bone components, respectively. The micromotions or displacements in the hip with implants were investigated [27, 28] and found to be 2.5 to 6 times higher in the composite bone than with the implants. This difference was a result of the mismatch between the E values of the bone and implant materials.
Our results are consistent with stress concentration on the bone surface via the body and surface stress. These stresses are concentrated on the first layers of the cortical bone which is several millimeters thick. Since we are assuming repeated cyclic loads in this study, damage likely concentrated on the surface comprised of the cortical bone. Since mechanism in the cancellous bone is displacement driven, the composite bone assumes that stress on both the zones will be same whereas the displacement will be different. Also, our results are consistent with femoral fractures observed clinically resulting from high stress.
The results of FEM analysis is presented in terms of both max von Mises stress and strain values (Figures 19 and 20), respectively, showing the composite laws and material properties as expected, that is, the displacement in the trabecular zone is higher, resulting in a higher strain than that in the cortical zone. Also, a higher total strain obtained life for the D equivalent of 1 is lower than at low total strains. A similar trend was found for the von Mises stress plot as well.
The limitations of this study can be the inability to validate the in vivo conditions in the absence of a biological self-healing environment. The second limitation is the computer, which makes it harder to apply more than 106 cycles during the damage simulation of the bone.
5. Conclusion
(1)Based on the nonlinear relationship of the macrodamage mathematical models of bone tissue, a conceptual model has been proposed and tested on a human femur. Monte Carlo simulation showed that the probability of failure for the proposed model was lower than that for the other models. The reason for this difference is that in this study the entire femoral bone was separated in terms of cortical and trabecular components.(2)The results have been validated using anisotropic material properties that showed the bone tissue damage cannot be expressed by only the cortical or the trabecular bone and both of them should be taken into consideration to develop a more realistic simulation.(3)Three damage prediction equations were developed (cortical, trabecular, and together cortical and trabecular). These equations can be used in deriving the bone fracture equations at a given stress range and fatigue life.(4)The study showed that the failure starts faster in the multiaxial loading condition than the uniaxial loading condition for the same number of cycles in the finite element simulation. Also, the damage starts at the femoral neck, as the femoral neck is the weakest part of the femoral bone.(5)The failure starts in the cortical bone before the trabecular bone. This means that the trabecular bone is more ductile while the cortical bone is more brittle.(6)The damage behavior seems to follow a three-stage regression; stage one was described by the primary phase of damage growth, stage two was described by the secondary phase of damage growth, and stage three was described by the tertiary phase of damage growth.(7)There is a significant difference in the stiffness of each element. Also, the stiffness of the cortical bone and the trabecular bone are significantly different as expected.Disclosure
An earlier version of this work was presented as an abstract at the 12th Annual Dayton Engineering Sciences Symposium (DESS 2016).
Conflicts of Interest
The authors declare that there is no conflict of interest of any nature regarding the publication of this paper.