Applied Bionics and Biomechanics

Applied Bionics and Biomechanics / 2017 / Article

Research Article | Open Access

Volume 2017 |Article ID 4539178 | 19 pages | https://doi.org/10.1155/2017/4539178

Macrodamage Accumulation Model for a Human Femur

Academic Editor: Estefanía Peña
Received22 Apr 2017
Accepted19 Jun 2017
Published29 Aug 2017

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 () [1416]. 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 [1416]. 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.


Modulus of elasticityPoisson ratioShear modulus

Cortical boneν12 = 0.4
ν23 = 0.25
ν31 = 0.25
Trabecular boneν12 = 0.4
ν23 = 0.25
ν31 = 0.25

represents the maximum density, G12 max = 5.71 MPa, G23 max = 7.11 MPa, and G31 max = 6.58 MPa. The superscript numbers denote the following: 1 for the radial direction, 2 for the circumferential direction, and 3 for the longitudinal direction [1416].

Material group numberDensity (g/cm3)E1 (GPa)E2 (GPa)E3 (GPa)ν12ν23ν31G12G23G31

10.9971.1511.1511.8940.400.250.250.0530.0660.061
20.9981.1531.1531.8970.400.250.250.0530.0660.061
30.9991.1541.1541.9000.400.250.250.0530.0660.061
41.0001.1561.1561.9030.400.250.250.0530.0660.061
51.0011.1581.1581.9060.400.250.250.0530.0660.061
61.0021.1601.1601.9090.400.250.250.0530.0660.061
71.0021.1621.1621.9120.400.250.250.0530.0660.061
81.0031.1641.1641.9150.400.250.250.0530.0670.062
91.0041.1661.1661.9180.400.250.250.0540.0670.062
101.0051.1681.1681.9210.400.250.250.0540.0670.062
111.0061.1701.1701.9230.400.250.250.0540.0670.062
121.0071.1721.1721.9260.400.250.250.0540.0670.062
131.0081.1741.1741.9290.400.250.250.0540.0670.062
141.0091.1761.1761.9320.400.250.250.0540.0670.062
151.0101.1781.1781.9350.400.250.250.0540.0670.062
161.0111.1801.1801.9380.400.250.250.0540.0680.063
171.0121.1811.1811.9410.400.250.250.0540.0680.063
181.0131.1831.1831.9440.400.250.250.0540.0680.063
191.0141.1851.1851.9470.400.250.250.0550.0680.063
201.0151.1871.1871.9500.400.250.250.0550.0680.063
211.0161.1891.1891.9530.400.250.250.0550.0680.063
221.0161.1911.1911.9560.400.250.250.0550.0680.063
231.0171.1931.1931.9590.400.250.250.0550.0680.063
241.0181.1951.1951.9620.400.250.250.0550.0690.063
251.0191.1971.1971.9650.400.250.250.0550.0690.064
261.0201.1991.1991.9680.400.250.250.0550.0690.064
271.0211.2011.2011.9710.400.250.250.0550.0690.064
281.0221.2031.2031.9730.400.250.250.0550.0690.064
291.0231.2051.2051.9760.400.250.250.0560.0690.064
301.0241.2071.2071.9790.400.250.250.0560.0690.064
311.0251.2091.2091.9820.400.250.250.0560.0690.064
321.0261.2111.2111.9850.400.250.250.0560.0700.064
331.0271.2131.2131.9880.400.250.250.0560.0700.065
341.0281.2151.2151.9910.400.250.250.0560.0700.065
351.0291.2171.2171.9940.400.250.250.0560.0700.065
361.0301.2191.2191.9970.400.250.250.0560.0700.065
371.0301.2211.2212.0000.400.250.250.0560.0700.065
381.0311.2231.2232.0030.400.250.250.0560.0700.065
391.0321.2241.2242.0060.400.250.250.0570.0700.065
401.0331.2261.2262.0090.400.250.250.0570.0710.065
411.0341.2281.2282.0120.400.250.250.0570.0710.065
421.0351.2301.2302.0150.400.250.250.0570.0710.066
431.0361.2321.2322.0180.400.250.250.0570.0710.066
441.0371.2341.2342.0210.400.250.250.0570.0710.066
451.03812.36312.36320.2400.400.250.255.3676.6836.185
461.03912.38312.38320.2690.400.250.255.3776.6956.196
471.04012.40312.40320.2990.400.250.255.3866.7076.207
481.04112.42312.42320.3290.400.250.255.3966.7196.218
491.04212.44212.44220.3590.400.250.255.4066.7316.229
501.04312.46212.46220.3890.400.250.255.4156.7436.240
511.04412.48212.48220.4190.400.250.255.4256.7556.252
521.04412.50212.50220.4490.400.250.255.4356.7676.263
531.04512.52212.52220.4790.400.250.255.4446.7796.274
541.04612.54212.54220.5090.400.250.255.4546.7916.285
551.04712.56212.56220.5390.400.250.255.4646.8046.296
561.04812.58212.58220.5690.400.250.255.4746.8166.308
571.04912.60212.60220.5990.400.250.255.4836.8286.319
581.05012.62212.62220.6290.400.250.255.4936.8406.330
591.05112.64212.64220.6590.400.250.255.5036.8526.341
601.05212.66212.66220.6890.400.250.255.5136.8646.353
611.05312.68212.68220.7190.400.250.255.5226.8766.364
621.05412.70212.70220.7490.400.250.255.5326.8896.375
631.05512.72212.72220.7790.400.250.255.5426.9016.386
641.05612.74212.74220.8100.400.250.255.5526.9136.398
651.05712.76212.76220.8400.400.250.255.5626.9256.409
661.05812.78212.78220.8700.400.250.255.5726.9386.420
671.05812.80212.80220.9000.400.250.255.5816.9506.432
681.05912.82212.82220.9300.400.250.255.5916.9626.443
691.06012.84212.84220.9610.400.250.255.6016.9746.454
701.06112.86212.86220.9910.400.250.255.6116.9876.466
711.06212.88212.88221.0210.400.250.255.6216.9996.477
721.06312.90212.90221.0510.400.250.255.6317.0116.489
731.06412.92312.92321.0820.400.250.255.6417.0246.500
741.06512.94312.94321.1120.400.250.255.6507.0366.511
751.06612.96312.96321.1420.400.250.255.6607.0486.523
761.06712.98312.98321.1730.400.250.255.6707.0616.534
771.06813.00313.00321.2030.400.250.255.6807.0736.546
781.06913.02413.02421.2340.400.250.255.6907.0856.557
791.07013.04413.04421.2640.400.250.255.7007.0986.569
801.07113.06413.06421.2940.400.250.255.7107.1106.580

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 [1821]. 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:


Material group numberDensity (g/cm3)C11C22C33C12C13C23C44C55C66

10.9971.4231.4832.1870.5340.2970.2370.0660.0610.053
20.9981.4261.5412.4550.5350.2970.2380.0660.0610.053
30.9991.4281.5462.4670.5360.2980.2380.0660.0610.053
41.0001.4311.5512.4790.5360.2980.2380.0660.0610.053
51.0011.4331.5562.4910.5370.2990.2390.0660.0610.053
61.0021.4351.5612.5030.5380.2990.2390.0660.0610.053
71.0021.4381.5662.5160.5390.3000.2400.0660.0610.053
81.0031.4401.5722.5280.5400.3000.2400.0670.0620.053
91.0041.4421.5772.5400.5410.3010.2400.0670.0620.054
101.0051.4451.5822.5520.5420.3010.2410.0670.0620.054
111.0061.4471.5872.5650.5430.3020.2410.0670.0620.054
121.0071.4501.5932.5770.5440.3020.2420.0670.0620.054
131.0081.4521.5982.5890.5450.3030.2420.0670.0620.054
141.0091.4541.6032.6020.5450.3030.2420.0670.0620.054
151.0101.4571.6082.6140.5460.3040.2430.0670.0620.054
161.0111.4591.6142.6270.5470.3040.2430.0680.0630.054
171.0121.4621.6192.6400.5480.3050.2440.0680.0630.054
181.0131.4641.6242.6520.5490.3050.2440.0680.0630.054
191.0141.4661.6302.6650.5500.3060.2440.0680.0630.055
201.0151.4691.6352.6780.5510.3060.2450.0680.0630.055
211.0161.4711.6402.6910.5520.3070.2450.0680.0630.055
221.0161.4741.6462.7040.5530.3070.2460.0680.0630.055
231.0171.4761.6512.7170.5540.3080.2460.0680.0630.055
241.0181.4781.6562.7290.5540.3080.2460.0690.0630.055
251.0191.4811.6622.7430.5550.3090.2470.0690.0640.055
261.0201.4831.6672.7560.5560.3090.2470.0690.0640.055
271.0211.4861.6732.7690.5570.3100.2480.0690.0640.055
281.0221.4881.6782.7820.5580.3100.2480.0690.0640.055
291.0231.4911.6842.7950.5590.3110.2480.0690.0640.056
301.0241.4931.6892.8080.5600.3110.2490.0690.0640.056
311.0251.4951.6952.8220.5610.3120.2490.0690.0640.056
321.0261.4981.7002.8350.5620.3120.2500.0700.0640.056
331.0271.5001.7062.8490.5630.3130.2500.0700.0650.056
341.0281.5031.7112.8620.5630.3130.2500.0700.0650.056
351.0291.5051.7172.8760.5640.3140.2510.0700.0650.056
361.0301.5081.7222.8890.5650.3140.2510.0700.0650.056
371.0301.5101.7282.9030.5660.3150.2520.0700.0650.056
381.0311.5121.7332.9170.5670.3150.2520.0700.0650.056
391.0321.5151.7392.9300.5680.3160.2520.0700.0650.057
401.0331.5171.7452.9440.5690.3160.2530.0710.0650.057
411.0341.5201.7502.9580.5700.3170.2530.0710.0650.057
421.0351.5221.7562.9720.5710.3170.2540.0710.0660.057
431.0361.5251.7612.9860.5720.3180.2540.0710.0660.057
441.0371.5271.7673.0000.5730.3180.2550.0710.0660.057
451.03815.29515.93223.3695.7363.1862.5496.6836.1855.367
461.03915.31915.95823.4045.7453.1922.5536.6956.1965.377
471.04015.34415.98323.4385.7543.1972.5576.7076.2075.386
481.04115.36816.00923.4735.7633.2022.5616.7196.2185.396
491.04215.39316.03423.5075.7723.2072.5656.7316.2295.406
501.04315.41716.06023.5425.7813.2122.5706.7436.2405.415
511.04415.44216.08523.5765.7913.2172.5746.7556.2525.425
521.04415.46616.11123.6115.8003.2222.5786.7676.2635.435
531.04515.49116.13723.6465.8093.2272.5826.7796.2745.444
541.04615.51616.16223.6805.8183.2322.5866.7916.2855.454
551.04715.54016.18823.7155.8283.2382.5906.8046.2965.464
561.04815.56516.21423.7505.8373.2432.5946.8166.3085.474
571.04915.59016.23923.7845.8463.2482.5986.8286.3195.483
581.05015.61416.26523.8195.8553.2532.6026.8406.3305.493
591.05115.63916.29123.8545.8653.2582.6076.8526.3415.503
601.05215.66416.31623.8885.8743.2632.6116.8646.3535.513
611.05315.68916.34223.9235.8833.2682.6156.8766.3645.522
621.05415.71316.36823.9585.8923.2742.6196.8896.3755.532
631.05515.73816.39423.9935.9023.2792.6236.9016.3865.542
641.05615.76316.42024.0285.9113.2842.6276.9136.3985.552
651.05715.78816.44524.0625.9203.2892.6316.9256.4095.562
661.05815.81216.47124.0975.9303.2942.6356.9386.4205.572
671.05815.83716.49724.1325.9393.2992.6406.9506.4325.581
681.05915.86216.52324.1675.9483.3052.6446.9626.4435.591
691.06015.88716.54924.2025.9583.3102.6486.9746.4545.601
701.06115.91216.57524.2375.9673.3152.6526.9876.4665.611
711.06215.93716.60124.2725.9763.3202.6566.9996.4775.621
721.06315.96216.62724.3075.9863.3252.6607.0116.4895.631
731.06415.98716.65324.3425.9953.3312.6647.0246.5005.641
741.06516.01216.67924.3776.0043.3362.6697.0366.5115.650
751.06616.03716.70524.4126.0143.3412.6737.0486.5235.660
761.06716.06216.73124.4476.0233.3462.6777.0616.5345.670
771.06816.08716.75724.4826.0333.3512.6817.0736.5465.680
781.06916.11216.78324.5176.0423.3572.6857.0856.5575.690
791.07016.13716.80924.5526.0513.3622.6897.0986.5695.700
801.07116.16216.83524.5876.0613.3672.6947.1106.5805.710

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.


Macrodamage modelMeanSDVarianceProbability of failure

The proposed femoral model using Miner’s rule0.3332990.3931340.15455413.26%
Zioupos and Casinos [26]0.8066670.9447870.89262237.90%
Pattin et al. [25]0.557340.2468680.06094442.20%
Hambli [13]0.1719420.3931340.0112916.63%

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.

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).

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.


The modelsThe equationsR2

The cortical bone model0.973
The composite bone model0.982
The trabecular bone model0.992

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.

References

  1. R. B. Martin, D. B. Burr, and N. A. Sharkey, Skeletal Tissue Mechanics, Springer, New York, 1998.
  2. J. Y. Rho, L. Kuhn-Spearing, and P. Zioupos, “Mechanical properties and the hierarchical structure of bone,” Medical Engineering & Physics, vol. 20, no. 2, pp. 92–102, 1998. View at: Publisher Site | Google Scholar
  3. J. Keyak and S. Rossi, “Prediction of femoral fracture load using finite element models: an examination of stress- and strain-based failure theories,” Journal of Biomechanics, vol. 33, no. 2, pp. 209–214, 2000. View at: Publisher Site | Google Scholar
  4. D. Garcia, P. K. Zysset, M. Charlebois, and A. Curnier, “A 1D elastic plastic damage constitutive law for bone tissue,” Archive of Applied Mechanics, vol. 80, no. 5, pp. 543–555, 2010. View at: Publisher Site | Google Scholar
  5. M. Fondrk, E. Bahniuk, and D. Davy, “A damage model for nonlinear tensile behavior of cortical bone,” Journal of Biomechanical Engineering-Transactions of the Asme, vol. 121, no. 5, pp. 533–541, 1999. View at: Publisher Site | Google Scholar
  6. P. Zlámal, T. Doktor, O. Jiroušek, and I. Jandejsek, “Verification of numerical model for trabecular tissue using compression test and time-lapse X-ray radiography based on material model determined from three-point bending test of single trabecula,” Key Engineering Materials, vol. 586, pp. 265–269, 2014. View at: Publisher Site | Google Scholar
  7. A. N. Natali, E. L. Carniel, and P. G. Pavan, “Constitutive modelling of inelastic behaviour of cortical bone,” Medical Engineering and Physics, vol. 30, pp. 905–912, 2008. View at: Publisher Site | Google Scholar
  8. D. Garcia, P. Zysset, M. Charlebois, and A. Curnier, “A three-dimensional elastic plastic damage constitutive law for bone tissue,” Biomechanics and Modeling in Mechanobiology, vol. 8, no. 2, pp. 149–165, 2009. View at: Publisher Site | Google Scholar
  9. H. Ridha and P. Thurner, “Finite element prediction with experimental validation of damage distribution in single trabeculae during three-point bending tests,” Journal of the Mechanical Behavior of Biomedical Materials, vol. 27, pp. 94–106, 2013. View at: Publisher Site | Google Scholar
  10. M. Charlebois, M. Jirasek, and P. Zysset, “A nonlocal constitutive model for trabecular bone softening in compression,” Biomechanics and Modeling in Mechanobiology, vol. 9, no. 5, pp. 597–611, 2010. View at: Publisher Site | Google Scholar
  11. R. Hambli, “Micro-CT finite element model and experimental validation of trabecular bone damage and fracture,” Bone, vol. 56, no. 2, pp. 363–374, 2013. View at: Publisher Site | Google Scholar
  12. H. S. Hosseini, M. Horák, P. K. Zysset, and M. Jirásek, “An over-nonlocal implicit gradient-enhanced damage-plastic model for trabecular bone under large compressive strains,” International Journal for Numerical Methods in Biomedical Engineering, vol. 31, no. 11, 2015. View at: Publisher Site | Google Scholar
  13. R. Hambli, “Apparent damage accumulation in cancellous bone using neural networks,” Journal of the Mechanical Behavior of Biomedical Materials, vol. 4, no. 6, pp. 868–878, 2011. View at: Publisher Site | Google Scholar
  14. W. Taylor, E. Roland, H. Ploeg et al., “Determination of orthotropic bone elastic constants using FEA and modal analysis,” Journal of Biomechanics, vol. 35, pp. 767–773, 2002. View at: Publisher Site | Google Scholar
  15. L. Peng, J. Bai, X. Zeng, and Y. Zhou, “Comparison of isotropic and orthotropic material property assignments on femoral finite element models under two loading conditions,” Medical Engineering & Physics, vol. 28, no. 3, pp. 227–233, 2006. View at: Publisher Site | Google Scholar
  16. D. C. Wirtz, N. Schiffers, T. Pandorf, K. Radermacher, D. Weichert, and R. Forst, “Critical evaluation of known bone material properties to realize anisotropic FE-simulation of the proximal femur,” Journal of Biomechanics, vol. 33, no. 10, pp. 1325–1330, 2000. View at: Publisher Site | Google Scholar
  17. G. Bergmann, HIP98, Loading of the Hip Joint. Compact Disc, Freie Universität, Berlin, 2001.
  18. T. Goswami, “Creep-fatigue life prediction—methods and trends,” High Temperature Materials and Processes, vol. 14, no. 2, pp. 21–34, 1995. View at: Google Scholar
  19. T. Goswami, “Dwell sensitivity, part I behavior and modeling,” Mechanics of Materials, vol. 22, no. 2, pp. 105–130, 1996. View at: Google Scholar
  20. T. Goswami, “Low cycle fatigue life prediction—a new model,” International Journal of Fatigue, vol. 19, no. 2, pp. 109–115, 1997. View at: Google Scholar
  21. T. Goswami, “Dwell sensitivity modeling—a new damage model,” Material and Design, vol. 25, no. 3, pp. 191–197, 2003. View at: Google Scholar
  22. M. A. Meyers and K. K. Chawla, Mechanical Behavior of Materials, Cambridge University Press, Cambridge, NY, USA, 2009.
  23. S. Manson, “Fatigue: a complex subject-some simple approximations,” Experimental Mechanics, vol. 5, no. 4, p. 193, 1965. View at: Publisher Site | Google Scholar
  24. S. Fatihhi, M. Harun, M. Kadir et al., “Uniaxial and multiaxial fatigue life prediction of the trabecular bone based on physiological loading: a comparative study,” Annals of Biomedical Engineering, vol. 43, no. 10, pp. 2487–2502, 2015. View at: Publisher Site | Google Scholar
  25. C. Pattin, W. Caler, and D. Carter, “Cyclic mechanical property degradation during fatigue loading of cortical bone,” Journal of Biomechanics, vol. 29, no. 1, pp. 69–79, 1996. View at: Publisher Site | Google Scholar
  26. P. Zioupos and A. Casinos, “Cumulative damage and the response of human bone in two-step loading fatigue,” Journal of Biomechanics, vol. 31, no. 9, pp. 825–833, 1998. View at: Publisher Site | Google Scholar
  27. B. Elliott and T. Goswami, “Implant material properties and their role in micromotion and failure in total hip arthroplasty,” International Journal of Mechanics and Materials in Design, vol. 8, no. 1, pp. 1–7, 2012. View at: Google Scholar
  28. M. T. Makola and T. Goswami, “Hip implant stem interfacial motion, a finite element analysis,” International Journal of Experimental and Computational Biomechanics, vol. 1, no. 4, p. 343, 2011. View at: Google Scholar

Copyright © 2017 Farah Hamandi and Tarun Goswami. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

821 Views | 443 Downloads | 1 Citation
 PDF  Download Citation  Citation
 Download other formatsMore
 Order printed copiesOrder