Journal of Biomedical Engineering and Biosciences (JBEB)

ISSN: 2564-4998

1. Introduction

The cervical spine is a complex section of the spinal column. It has several components that work together to provide humans with a large range of motion [1-8] in the neck and protects the spinal cord. Therefore, the recent literature shows how to ensure the health, strength and mechanical properties of the vertebrae [9-28]. The objective of this research is to expand fracture toughness considerations in cortical bone. Inherent challenges with these types of studies involve varying structures of the cortical bone and the interfacial regions between cortical and the trabecular structure. Several features create voids in cortical bone including Haversian canals, Volkmann's canals canaliculai and randomly distributed pores.

Several studies have investigated the limit of the fracture toughness of cortical bone [9-16]. The stress intensity factor of bone can be described in terms of longitudinal and transverse fracture toughness. Longitudinal stress intensity is associated with the long axis of long bones, i.e. the femur. Irregular bones specifically the vertebrae have a less well defined longitudinal and transverse direction. Physiologically the load sharing between the trabecular core and the endplate is complex phenomena. When a healthy disc is present, the endplate is loaded under tensile stresses. However, when the disc is replaced the forces transferred cause a compressive stress in the endplate. However, some regions of the endplate may experience tension [29]. The cortical endplate of the vertebral body was modeled assuming a simply supported beam, to determine the strength of the endplate and its mechanical behavior.

A probabilistic assessment of the toughness was developed in the presence of pores [30]. The pores were assumed from femur studies compiled from literature [31-34] and assumed normally distributed. The probabilistic methodology developed here considers fracture toughness limits. A beam model was constructed with predefined cracks present in the tension region of a beam. The tension region of a simply supported beam is the area opposite side of the neutral axis from the load. The beam models were also considered plain stress models. Four fracture models were considered that defined the stress intensity at crack tips. Three of these models Vashishth [12], Dowling [35] and Feng [15] considered compact tension fracture, in plane strain conditions. The fourth model was a three point beam bending model constructed by Yan [13]. The crack size was examined to determine the probability when it would exceed the fracture toughness of the bone.

2. Methods

A unit width was assumed with an initial beam length to simplify calculation steps. A mid-sagittal section was assumed with a thickness function that was developed from studies of the endplate measurements. As a result the cross sectional area varied in depth only with position along the beam length. A second order polynomial fit was developed applied to the data presented in Table 1. The model is consistent with findings of the cortical endplate distributions. The thickest regions are the peripheries of the endplates or the ends of the beam by the support and the thinnest regions of the endplate are in the center or the middle of the beam. Once the geometry was established, an assumption was made for the elastic modulus of the cortical bone. Several values were investigated and the typical range is between 16-20 GPa. This model uses 18.6 GPa as the modulus of elasticity (E). The beam is assumed to behave linear elastically. Modulus data is presented in Table 2.

The following assumptions were made for this paper; cortical bone acts as an isotropic, homogenous, linear elastic material, the endplate behaves like a beam, the radial cortical shell provides pinned supports for the beam, single cracks are normally distributed, influence of cracks do not interact with each other, shape factors for compact sections apply to individual beam segments, the shape function for the beam applies to cracks at any position x, and finally uniaxial stress conditions prevail.

Two load cases were considered as shown in figure 1:

Load Case I

• - Simply supported beam
• - Pin/roller supports
• - Two symmetrically placed static point loads
• - Load simulates disc with poor placement mismatch of implant and vertebrae radius

Load Case II

• - Simply supported beam
• - Pin/roller supports
• - Uniformly distributed static load
• - Load placed symmetrically on middle 3/4 of beam
• - Load simulates disc with good implant and vertebrae radius match

The load applied to the beam is of equal force for both cases. It is distributed in two different conditions as previously described. The first case models an implanted disc in contact with the vertebral body in only two points of contact (simulating drastic curvature differences the design aid of two equal concentrated loads symmetrically placed was used) and the second with complete contact between an artificial disc and vertebral body (where contact is continuous along the length of the implant the uniform load partially distributed design aid was used.). Physiologically the vertebral body shares load between the cortical shell and the trabecular core [38]. The magnitude of the load at a minimum comes from the head and the contribution of the length of the neck to the level of the vertebrae under consideration. Conservatively the entire length of the neck and the weight of the head can be considered to act on the superior endplate of the vertebrae [42].

The actual load applied was a static load of 200 N. The load is also considered purely static without cycling. According to anthropometric measurements the weight of the head and neck is 8 percent of the total body weight of the person under consideration [42]. For a 200 pound person that head and neck segment weighs 16 pounds. In SI units the head and neck weighs 71.2 N. This is just the load from the head and neck which the spine sees constantly in an upright position. A larger load was considered to investigate more extreme loads. 200 N was used as a benchmark. An endplate width and depth was estimated from measurements taken of the cervical spine [43]. A depth of anterior to posterior was estimated at 20 mm. If it is assumed that the entire endplate distributes load evenly then on a per unit width the 200 N load is divided by 20 mm giving 10 N/mm through the depth. Since this theoretical beam is 1 mm thick the load applied to the unit width under investigation is 10 N. This ensures K_Ic from reported literature can be used as opposed to a lower K_Ic due to cyclic loading.

The stress was calculated through the length of the beam assuming basic beam bending equations in which planar sections remain plane. This was done to limit the crack opening stresses in the transverse plane and eliminating shear that would cause the opening of cracks in directions inducing Mode II or Mode III behavior. The random pore distribution from cortical femur bone was used to develop a distribution of pores within the vertebrae. Pore sizes in the vertebrae were assumed consistent with femur data [31-34]. This information was used as a basis to form a crack. The pores from the data were assumed to be circular. This allowed for the calculation of the diameter of the pore that was used as the basis of the length of the crack. The diameter was then considered the nonunion of adjacent differential slices on the beam in the tension section creating a sharp, through edge crack condition. This crack length was important in Mode I crack opening and respective stress intensity factors. The crack was also considered a small edge crack in a homogenous material. The crack direction was perpendicular to the beam section at every point consistent with transversely oriented pores in vertebral cortical bone [19]. A Monte Carlo extension module, RiskAMP [44] (Structured Data, New York), was used to create the crack simulation. RiskAMP is a random number generator. A truncated normal distribution was used to generate the crack lengths. It consisted of a mean, standard deviation, maximum and minimum crack length. This model assumed a crack from a sample size of 500. The Monte Carlo simulator establishes the crack distribution parameters for every crack position along the length of the beam. The large number means the distributions are consistent in terms of mean and standard deviation. Secondly, it increases the accuracy of the results. Porosity data is presented in Table 3.

The initial crack length is established the stress intensity was calculated for each of the initial crack lengths. Four models were developed following ASTM standards [13-16, 45]. Each of the models introduced a shape function: Norman [14] and Feng [15] used a shape function for the stress intensity factor developed by Behiri and Bonfield. Others [16,35] used a compact tension shape function while Yan used a three point beam bending function per ASTM E-399 [45].The shape functions were, with respect to Mode I crack opening, in terms of 'a', the initial crack length, and W the width of the specimen at that point. In this case the specimen width corresponds to the height of the beam. While the previously mentioned models are a function of a/W the Dowling model is a function of x position along the beam. Therefore the beam was divided into 100 differential slices with which to determine the stress intensity. The same spots were used with each respective a/W ratio for comparison.

Each of these shape functions and test methods had specific specimen size requirements on them according to ASTM E399 testing standards. Typically the requirement for a/W is in the range of 0.45 to 0.55. Feng [15] expanded the a/W initial condition to 0.7. Other studies have shown that an increase in the ratio of a/W (longer initial crack) increases the resulting Kc value with that associated initial condition. ASTM E399 [45] recommends test specimen configuration. The point of load application is where W is considered to start when measuring the ratio of a/W. This theoretical model assumes that the load application to each compact specimen is supplied by a bending moment perpendicular to the crack orientation. The beam model has restrictions on the size of the span to width of the specimen. This ratio is on the order of 4 to 5. It was assumed that the equations would hold for the decreased specimen width and plain strain conditions prevail.

An initial crack length 'a' was estimated and corresponding stress intensity factor calculated. A limit crack length was then determined based on the fracture toughness of transverse cortical bone. Once the critical crack length was determined the probability of a crack existing with a shorter length was determined. The analysis shows Haversian canals long axis oriented in the axial plane [46]. Because of the orientation of the osteons cracks would grow perpendicular to the long axis of the osteon. However when a crack grows into an osteon, it would wrap around the osteon thereby increasing the stress intensity value as compared to a crack growing parallel to the osteons long axis, crack would open the osteons like a zipper. Once the size of the transition crack length was determined a comparison made to the distribution of cracks at the same spot along the beam developed by the Monte Carlo simulations.

3. Results

The probability of a crack length shorter than the critical length was determined for both, Load Case 1 and Load Case 2. These results are shown in Figures 2 and 3. It can be seen that the Yan beam model [16] is the most conservative model. For case 1 and case 2, 90 percent of the beam was too thin for the presence of any crack. Stated another way the minimum length of cracks sizes present created an a/W ratio was too large for the depth of the beam at the corresponding point. A three point beam model does not accurately reflect the support of the trabecular core, which reduces the span length. While the support is not the same as in a bend test a modification to the span coefficient may yield results closer to that of the other model.

Corresponding to the position of the beam allowable crack length values were presented for the two load cases (Figures 4 and 5). A critical crack length was determined by solving for the crack length at the transverse critical stress intensity. Yan's model [13] was the least likely to have subcritical crack lengths, with the exception of small distances from the supports. Additionally Dowling's model is the most likely to have a subcritical crack present with Vashishth's model in the middle.

For both load crack present would be critical thereby exceeding the critical transverse stress intensity. Dowling's mode scenarios it is clear that all three models show similar behavior. Vashishth's model [16] indicates that a large portion of the beam model would contain subcritical cracks. However towards the center of the beam negative crack lengths arose. This indicates any l [35] created a ceiling for the amount of permissible cracks. The crack limit threshold was also greater for every point along the beam than the Vashishth and Yan models.

Figures 6 through 7 show the probability of the presence of a supercritical crack. For both load case 1 and 2 data is accumulated around 6 MNm^3/2. This indicates that at each stress intensity value there is a probability of exceeding or not exceeding bone's transverse critical stress intensity. When a prediction equation is fit to the data the trends for all four fracture models using the probabilistic methods were to the left hand side of a majority of the original models. This indicates a conservative prediction for the existence of a supercritical crack. This would mean that device design based on this method would be more conservative in predicting subsidence and the onset of yielding and/or fracture. For both load cases 1 and 2 the most conservative prediction model was Dowling's model followed by a close tie between Feng and Vashishth's model with Yan's model being the least conservative.

4. Discussion and Conclusion

This research develops a probabilistic assessment method for determining the presence of cracks of sufficient length to exceed the critical fracture toughness of bone. This technique is needed because the porous structure of bone causes stress concentrations that initiate failure before bone reaches its yield stress. This is important because bone failure cannot be characterized by yield mechanics of gross material properties alone. The work expands bone fracture research to a non ASTM standard model. This expanded work increases the applicability of the fracture work to areas of biological crack mechanics not well described by ASTM specimen size requirements. It also establishes the presence of pores that change the type of failure that is considered mechanically

This model assumes that plain strain stress intensity applies to the beam model because of the size constraints of the vertebral cortical bone. It was assumed that there would be no interaction between plain strain and plain stress. Mode I is assumed to act independently of Mode's II and III. The entire applied load was considered to activate Mode I only. The Vashishth '04 [16] model predicts 52% of the Case 1 and 65% of the Case 2 beams cannot meet a minimum required crack limit. The Dowling, which is also based on a compact tension model, has the highest probability of a crack below the fracture limit. In Case 1 the entire beam had a chance of the strength being governed by yielding. This is most likely due to the maximum moment. The max moment from Case 1 is 12.5 Nmm while the max moment for Case 2 is 31.3 Nmm. While the load was the same the moment distribution tends to increase. Dowling's model also consists of a discrete moment at each section of the beam. This moment can be found directly from the bending equations.

Differences in the probability are most likely attributed to the shape factor. Feng's models have an entirely random distribution that seems to be more dependent on crack length than it does on the amount of stress at each point. The Feng model differs from the previous three models in this respect. The other three models tend to be more dependent on the stress applied to the differential sections as opposed to being controlled by the crack length as with the Feng model. The Yan beam model [13] is the least conservative limit in both case 1 and case 2 for the stress intensity factor. This seems to be due in strong part to the span multiplier present in the equation for KI. For the case of vertebral endplate tissue separated from the cancellous core, a value of 20 mm would be correct if compared to three-point beam testing. This is an overestimate of the stress intensity value at points along the length because the trabecular core supports the vertebral endplate at intervals closer than the 20 mm span. The stress intensity values calculated by the Vahishth and Feng models are the most similar with the Dowling model being the least conservative model for both load case 1 and 2. The shape function of Dowling's model undervalues applied stress compared to other models. That does not necessarily mean that a crack above the crack limit threshold exists there. The probability of exceeding the threshold crack limit was calculated verse the KI value. A Weibull cumulative distribution plot (CDF) was used to describe the probability as a function of KI. The Weibull equation and the probability domain were graphically cut off after the point at which the crack threshold maintained 100% probability. This theoretical model produced the probability of a threshold crack based on a limiting Kc of 6 MNm^3/2. The randomized values are grouped around the transverse Kc limit of 6. The probability of a sufficient crack was fit with a Weibull continuous distribution function (CDF) [47] in the form of the equation:

(1)

This model is also limited to Mode I, tensile, failure. Bone is subjected to more than tensile forces. Expanding the study to consider Mode II and III fracture and the associated probabilities of fracture would improve bone strength assessments. This research does not consider the overall reduction in strength; however, it could be included in the analysis as a reduction to bending stress or deflection limits. It is also recommended to further examine the length multiplier for the Yan beam model. A shorter unsupported length would give more a representative stress intensity values per crack length.

References