Predicting the yield of the proximal femur using high-order finite-element analysis with inhomogeneous orthotropic material properties

Zohar Yosibash , David Tal , Nir Trabelsi

Abstract

High-order finite-element (FE) analyses with inhomogeneous isotropic material properties have been shown to predict the strains and displacements on the surface of the proximal femur with high accuracy when compared with in vitro experiments. The same FE models with inhomogeneous orthotropic material properties produce results similar to those obtained with isotropic material properties. Herein, we investigate the yield prediction capabilities of these models using four different yield criteria, and the spread in the predicted load between the isotropic and orthotropic material models. Subject-specific high-order FE models of two human femurs were generated from CT scans with inhomogeneous orthotropic or isotropic material properties, and loaded by a simple compression force at the head. Computed strains and stresses by both the orthotropic and isotropic FE models were used to determine the load that predicts ‘yielding’ by four different ‘yield criteria’: von Mises, Drucker–Prager, maximum principal stress and maximum principal strain. One of the femurs was loaded by a simple load until fracture, and the force resulting in yielding was compared with the FE predicted force. The surface average of the ‘maximum principal strain’ criterion in conjunction with the orthotropic FE model best predicts both the yield force and fracture location compared with other criteria. There is a non-negligible influence on the predictions if orthotropic or isotropic material properties are applied to the FE model. All stress-based investigated ‘yield criteria’ have a small spread in the predicted failure. Because only one experiment was performed with a rather simplified loading configuration, the conclusions of this work cannot be claimed to be either reliable or sufficient, and future experiments should be performed to further substantiate the conclusions.

1. Introduction

Many patient-specific finite-element (FE) analyses of the femur that mimic its mechanical response have been reported in the past 10 years, and validated by experimental in vitro experiments. In all these analyses an inhomogeneous isotropic material property was assumed with a constant Poisson ratio (Keyak et al. 1990; Martelli et al. 2005; Schileo et al. 2007; Trabelsi et al. 2009) and the FE computed strains (and in some publications also the displacements) were compared with those measured in the in vitro experiments. Overall, a good correlation is reported between the FE results and the measurements. The next natural step was to attempt to predict the ‘yield’ or ‘fracture’ load also by these FE analyses. For such predictions a ‘yield criterion’ had to be adopted and, once again, isotropic (at times also inhomogeneous) yielding criteria were assumed, such as the von Mises, the Drucker–Prager, maximum principal strain and maximum principal stress (Keyak et al. 2001; Keyak & Falkinstein 2003; Bessho et al. 2007; Schileo et al. 2008). In these cases the correlation between the predicted load and the experimental observations was not very strong. In Keyak et al. (2001), 18 pairs of femurs were loaded to failure and FE models constructed and failure predicted (using the von Mises criterion) under two loading conditions, one similar to joint loading during single-limb stance and one simulating impact from a fall. For the stance condition, the predicted and actual fracture locations agreed in 13 of the 18 cases (72% agreement). For the fall condition, the predicted and actual fracture locations agreed in 10 of the 15 cases where the actual fractures could be identified (67% agreement). In Bessho et al. (2007) 11 right femurs were loaded to failure at a 20° angle, and FE models constructed to mimic the experiment with the Drucker–Prager criterion employed to predict the yield loads. The correlation between the measured and predicted yield load was estimated to be: load (measured)=1592+0.809 load (predicted) (N), with a coefficient of linear regression r2≈0.88. Schileo et al. (2008) demonstrated that the maximum principal strain failure criterion in conjunction with FEMs correlates fairly well with failure predictions in three femurs subjected to stance loading.

In these studies, no quantitative comparison of stress-based and strain-based yield criteria was performed, nor realistic orthotropic material properties taken into consideration.

Recently, high-order FEMs (p-version of the FEM: p-FEMs) were suggested (Yosibash et al. 2007a,b) for the simulation of the proximal femur mechanical response. p-FEMs have many advantages over conventional FEMs: accurate surface representation and faster numerical convergence rates achieved by increasing the polynomial degree p of the shape functions over the same mesh, thus controlling numerical errors easily. Also, the inhomogeneous Young modulus may vary widely within each element, the elements may have large aspect ratios (required in cortical regions that are thin and long) and the elements may be far more distorted than in conventional FEMs (Szabó & Babuška 1991). These p-FE analyses of the human femur with inhomogeneous isotropic and orthotropic material properties have been shown to predict the mechanical response with high accuracy (Yosibash et al. 2007b, 2008; Trabelsi et al. 2009). The inhomogeneous orthotropic material properties were determined from quantitative computed tomography (QCT) scans based on micromechanics considerations, whereas the inhomogeneous isotropic material properties are empirically based, having Young’s modulus distributions according to the QCT scans with a constant Poisson ratio. We noticed that both the isotropic and orthotropic FE models predict a similar mechanical response when the principal strains along the shaft or the displacements are compared with experimental observations (Trabelsi et al. submitted). Because of the several possible isotropic ‘yielding criteria’ used to predict the ‘yield load’ we herein investigate the following.

  • — How large is the spread in the FE-predicted yield load when using the four different yield criteria when isotropic inhomogeneous material properties are used in the FE model?

  • — How large is the spread in the FE-predicted yield load when orthotropic inhomogeneous material properties are used in the FE model?

  • — Which material property (isotropic or orthotropic) assignment in combination with which ‘yield criteria’ best predicts the yield load obtained by a simple experiment?

To try and provide some answers to the above items, we constructed p-FE models from QCT data obtained from two different fresh frozen femurs. In each FE model both isotropic and orthotropic inhomogeneous material properties were assigned as follows.

The strains and displacements predicted by the p-FE simulations were extracted and used to determine the location at which failure starts and the magnitude of the yield force using four different yield criteria. We also performed a simple test on one proximal femur in which we increased the load on the femur head until fracture, monitoring strains and load and displacement were seen. The FE-predicted yield force and location were compared with the yield force measured in the experiment. To the best of our knowledge, this is the first attempt to use FE models of femurs with orthotropic material properties for the prediction of yield load.

2. Methods

(a) p-FE analyses for FF1 and FF3 femurs

Two fresh frozen femurs, one from a 30-year-old male donor, denoted as FF1, and another from a 54-year-old female donor, denoted as FF3 (the femur numbers match our previous publications; Trabelsi et al. 2009), were QCT scanned and manipulated to generate high-order FE bone models. The three-dimensional reconstruction of femoral geometry and the generation of the FE mesh for the bone model are detailed in Trabelsi et al. (2009) and Yosibash et al. (2007a,b) and will only be summarized herein. Two different FE models were considered for each femur: a model with orthotropic inhomogeneous micro-mechanics-based material properties as described in Trabelsi et al. (submitted) and a model with empirically based inhomogeneous isotropic material properties detailed in Keyak & Falkinstein (2003). The isotropic material properties are a constant Poisson ratio ν=0.3 and a varying Young’s modulus computed as follows (Keyak & Falkinstein 2003): Embedded Image 2.1 Embedded Image 2.2 Embedded Image 2.3and Embedded Image 2.4

The geometry of the femur was extracted from QCT slices and divided into cortical and trabecular regions using an in-house automatic algorithm detailed in Yosibash et al. (2007b) (with HU>600 taken as the cortical region). Exterior, interface and interior boundaries were traced at each slice. A smoothing algorithm was applied to generate smooth closed splines used for the solid body generation. The solid body was meshed by an auto-mesher with tetrahedral elements using the p-FE code StressCheck.1 The surfaces of the bone are accurately represented in the FE model by using the blending mapping method. The entire algorithm is illustrated in fig. 1 in Trabelsi et al. (submitted). The two resulting p-FE models are shown in figure 1.

Figure 1.

p-FE models of (a) FF1 and (b) FF3 femurs.

Both FF1 and FF3 models were clamped at the distal face and compressed by a 1000 N force on the femur head in the shaft direction (0° tilt). An FE analysis was performed, increasing the polynomial degree until p=4, resulting in 158 967 d.f. and 4486 p-elements for the FF1 model and 3741 elements and 133 743 d.f. for FF3. Because at p=4 the error in energy norm had converged beneath 2 per cent relative error, and the pointwise strains converged, we did not extend the polynomial degree beyond it.

For both FF1 and FF3 26 elements were selected at the upper and 26 at the lower side of the femur neck—see figure 2. These 52 elements were those in which the stresses and strains were the highest, and thus are considered to be the suspected elements in which yielding is most likely to occur. The four different yield criteria (details in the sequel) were applied to these 52 elements and, among them, the 10 most susceptible for yielding were determined: six at the upper part of the neck (tension) and four at the lower part of the neck (compression)—these elements are the numbered elements shown in figure 2. All further analyses were applied to these 10 elements.

Figure 2.

(a) FF1 and (b) FF3—the 52 elements with highest strains and stresses, and the 10 elements most susceptible to yielding using the various yielding criteria.

(b) Yield criteria

Experiments performed on small pieces of human femurs taken from either the cortical or the trabecular regions have shown that a compression yield stress (which is a function of the bone density) can be defined as follows (Keller 1994; Keyak et al. 1994):Embedded Image 2.5 whereas the tensile yield stress isEmbedded Image 2.6 Other studies measured the yield strain instead, assumed to be independent of the density, which for the cortical region is (Bayraktar et al. 2004)Embedded Image 2.7 andEmbedded Image 2.8

To determine whether the bone tissue reaches the yield stress/strain, we considered the four isotropic yield criteria used in Keyak & Falkinstein (2003), Bessho et al. (2007) and Schileo et al. (2008), summarized herein.

Maximum principal strain. The principal strains are computed, and Embedded Image is determined by the maximum tensile or compressive principal strain. The risk factor (RF) for the yield is the ratio between this maximum principal strain and the appropriate (tensile or compression) yield strain in equations (2.7) and (2.8), and the yield load is the ratio between the applied load and RF,Embedded Image2.9

Maximum principal stress. The principal stresses are computed, and Embedded Image is determined by the maximum tensile or compressive principal stress. The RF for yield is the ratio between this maximum principal stress and the appropriate (tensile or compression) yield stress in (2.5) and (2.6) (which depends on the density), and the yield load is the ratio between the applied load and RF,Embedded Image2.10

Drucker–Prager criterion. This criterion is an extension of the von Mises yield criterion taking into account also the influence of a hydrostatic stress component by inclusion of an additional term which is proportional to the hydrostatic stress. In brittle materials these play a more important role in yielding prediction. The equivalent yield stress is computed byEmbedded Image2.11where I1 denotes the first invariant of the stress tensor, and J2 is the second invariant of the deviatoric stress tensor,Embedded Image2.12The parameter α is set as 0.07 in Bessho et al. (2007) according to a former analysis on concrete (Kupfer & Gerstle 1973). It reflects the dilative potential of the material, related to the proportions of the volumetric and deviatoric strains (taking into consideration the hydrostatic pressure in addition to the deviator stress in von Mises’s failure law).

The RF for fracture is the ratio between σDP and the appropriate (tensile or compression) yield stress in equations (2.5) and (2.6) (which depends on the density), and the yield load is the ratio between the applied load and RF.

von Mises stress. The equivalent von Mises stress is σDP(α=0), and the procedure of determining RF and force to failure is identical to that described above. This criterion is suitable for ductile materials such as metal; however, it is not as suitable for bones as they are classified as brittle materials. Nevertheless, it has been used in some past publications and therefore will be considered herein also.

(c) FE prediction of the yield location and yield force

For each of the 10 elements shown to be the most susceptible to failure in §2.1 we evaluated three different measures: (i) averaged data on the outer face of an element (see figure 3a,b), (ii) averaged data in the element’s volume (see figure 3c,d), and (iii) the extremum value in the element.

Figure 3.

(a,b) Typical grid of points for averaging surface data and (c,d) typical grid of points for averaging volume data in an element.

The equivalent stress and principal stresses and strains (σeq,εi and σi) were extracted at 28 grid points on the element’s outer surface and 84 grid points in the interior of the element. Denoting by ε1 the maximum principal strain, the averaged maximum principal strain for locations under tension is computed byEmbedded Image The averaged maximum principal strain for locations under compression is the same, except that |ε3| is used instead of ε1. The extremum value is determined by Embedded Image or Embedded Image.

To determine the strength in tension or compression, the average density at the same n points was computed as follows:Embedded Image and the yield tension or compression stress is evaluated using (ρash)Avg.

Having all the required information, we computed the RF by four different criteria and the three different measures. The yield location is defined in the element in which the highest RF is obtained, and the fracture load is the applied load over the RF.

(d) Experiment on FF1—load to fracture

To have some confidence in the FE results, we performed one experiment (which of course is too few for validation purposes) on the FF1 femur. It was loaded in vitro until fracture at a constant rate of 1 cm min−1 with uniaxial strains along principal directions bonded at four locations in the neck area and head displacements in the two directions measured. Strain-gauge locations and directions were determined based on a past FE analysis of the same bone (Trabelsi et al. 2009). These measurements are intended to identify the ‘yield load’, i.e. when a deviation from a linear load–strain and load–displacement relationship is observed. It is important to mention that this fresh frozen femur was defrosted for 2 days approximately 4 years before this study and was used for mechanical testing with loads up to 1500 N. Thereafter, it was frozen again, kept in a deep-freezer and defrosted again for the fracture test. In figure 4 the test setting is presented, with the load being applied at an angle of 0° to the shaft.

Figure 4.

FF1: (a) test fixture, (b) fractured bone.

3. Results

(a) Yield load and yield location predicted by FE models with isotropic and orthotropic inhomogeneous material properties

Results from the various p-FE models with the four different material properties are obtained and compared with the experimental observations. All FE models show good convergence in the energy norm above p=4 (error<2%) as well as good convergence in pointwise stress and strain values.

The reasoning for the interest in predicting the femur’s yield load (as opposed to the fracture load) is as follows: (i) we wish to be able to trigger the instance of damage being created in the bone, which starts to accumulate immediately after yield (which is believed to be caused by micro-failures at the osteon and lamella level) and (ii) the use of linear elastic analysis is sufficient, therefore there is no need to address non-linear complex and non-validated models that are required to describe the damage evolution (unknown at present).

A summary of the computed yield loads by the four different criteria for the isotropic and orthotropic FE models at the 10 elements of interest in FF1 is presented in figure 5 (we include also the FE model with the elements of interest numbered on it so that the reader can easily identify the elements). VM is the von Mises criterion, DP is the Drucker–Prager criterion, Embedded Image is the maximum principal stress criterion, and Embedded Image is the maximum principal strain criterion.

Figure 5.

Predicted yield load in the isotropic and orthotropic FE models of FF1 for the three different measures. (a) Isotropic face average; (b) orthotropic face average; (c) isotropic element average; (d) orthotropic element average; (e) isotropic maximum in element; (f) orthotropic maximum in element; (g) the 10 elements most susceptible to yielding. VM, von Mises criterion; DP, Drucker–Prager criterion; Embedded Image, maximum principal stress criterion; Embedded Image, maximum principal strain criterion.

A summary of computed yield load by the four different criteria for the isotropic and orthotropic FE models at the 10 elements of interest in FF3 is shown in figure 6. The fracture location in the femurs is anticipated to be in the element in which the lowest yielding load is predicted. To determine the variation in the prediction of the yield load between stress-based and strain-based criteria, and between the isotropic and orthotropic models, we display the large amount of data presented in figures 5 and 6 in a summarized form. For example, when considering the face average method, we generated plots as shown in figure 7, summarizing the lowest predicted yield loads in the three most ‘critical’ elements in the isotropic and orthotropic models of FF1 and FF3. The maximum principal strain criterion is shown as a single circle whereas the three stress criteria are represented as a range bar.

Figure 6.

Predicted yield load in the isotropic and orthotropic FE models of FF3 for the three different measures. (a) Isotropic face average; (b) orthotropic face average; (c) isotropic element average; (d) orthotropic element average; (e) isotropic maximum in element; (f) orthotropic maximum in element; (g) the 10 elements most susceptible to yielding. VM, von Mises criterion; DP, Drucker–Prager criterion; Embedded Image, maximum principal stress criterion; Embedded Image, maximum principal strain criterion.

Figure 7.

Face average method: yield load in the three most critical elements by the strain and stress criteria in FF1 and FF3 in the (a) isotropic and (b) orthotropic FE models (solid line, stress-based failure criteria; open circle, strain-based failure criteria).

Then the lowest predicted yield load by the maximum principal strain criterion (a single value) and the range of the predicted yield load by the various stress-based criteria in both femurs is shown in figure 8 (where the experimental load for FF1 is also provided).

Figure 8.

Predicted yield load by the strain and stress criteria in FF1 and FF3 in the most critical elements. (a) Face average; (b) element average; (c) element maximum (solid line, stress-based failure criteria; open circle, strain-based failure criteria).

(b) Yield load and yield location in the experiment on FF1

To determine the yield load, we plotted the load versus displacement of the head, and load versus the strain measured by strain-gauge number 3 (closest to the failure location) in figure 9. At a given load, the femur’s response in terms of both strains and displacement deviated from a linear response. To determine in a more precise manner this deviation, we computed the intersection of the response with a 5 per cent deviation of the linear slope. This is because a 0.2 per cent strain (2000 μstrain) ‘yield strain’ would be far beyond the mechanical response of the bone. According to this criterion the yield load is 4200–4800 N. Our hypothesis is that fracture started approximately 1–2 mm away from the predicted ‘yield element’ no. 3—in figure 10 we superimposed the FE mesh on the actual picture of the fractured bone and marked element 3 at which it is believed that yielding occurs first (we did not use a high-speed camera to monitor the fracture evolution). It is worthwhile to note that no crushing was observed in the head area during and after the loading.

Figure 9.

(a) Load (N) versus displacement (mm), (b) load (N) versus strain SG no. 3 (solid line, experiment; grey dashed line, linear (1000–3000 trendline); black dashed line, linear (0.95 offset)).

Figure 10.

The FE mesh with the element at which yielding is predicted superimposed on the picture of the broken femur.

4. Discussion and conclusions

Although we have reported predictions based on a pointwise value, we feel that this measure is overconservative and very sensitive to small changes in geometry and material properties; therefore, it is not a good measure for yield initiation (especially since yielding is probably spread over a small area). Hence, the results based on pointwise maximum values were discarded.

The surface averaging technique provides the smaller spread between the various yielding criteria, a more logical interpretation (the failure is predicted to start where stresses/strains are the highest) and at the same time predicts a yield load that is the closest to the experimental observation.

Inspecting the FE results for the models with isotropic material properties, one may note that the four different failure criteria predict the yield load with a relative spread. The orthotropic models, on the other hand, predict a yield load that has a small spread between strain and stress-based failure criteria. It is worthwhile to emphasize that the quality of the FE results of both isotropic and orthotropic models were verified to be converged by increasing the polynomial degree. Therefore the numerical error was small enough so further mesh refinement or increase of the polynomial degree would not have changed the results.

View this table:
Table 1.

Summary of yield load (N) and the element number where yield is predicted (in parenthesis) for FF1 and FF3 using the average surface measure. VM, von Mises criterion; DP, Drucker–Prager criterion; Embedded Image, maximum principal stress criterion; Embedded Image, maximum principal strain criterion.

Yield load predicted by the strain-based criterion is lower than those predicted by the various stress criteria; therefore, this is a conservative criterion.

All criteria (stress-based and strain-based) predict the yield load at the location of actual fracture in the tested bone FF1. Therefore, it is not surprising that previous publications successfully predicted the location of fracture by these stress-based failure criteria (Keyak et al. 2001; Bessho et al. 2007) using isotropic material properties. The magnitude of the yield load was predicted with a moderate accuracy in these studies, and shown to be inferior to the strain-based criterion in Schileo et al. (2008) using isotropic material properties.

As expected, for the loading condition prescribed, the regions under tension are under a very large tensile stress, which implies yielding in tension and not in compression. Under the specific loading, it is conceivable that the yield starts under tension at the outer surface of the bone so that the averaged element surface measure is the most logical measure to predict the location and yield load. In table 1 we summarize the yield load predicted by the orthotropic and isotropic FE models. Indeed the strain criterion is the most appropriate (a conclusion supported by Schileo et al. (2008) also) for yield load prediction. This criterion also predicts well the location at which the fracture initiates when compared with the single available experimental result.

The orthotropic FE model in conjunction with the maximum principal strain criterion averaged over the element’s face area results in a very good prediction of the yield load and yield location (with good agreement between the FF1 experiment and predicted values). This strain criterion was also shown to be the least sensitive whether isotropic or anisotropic material properties were assigned to the FE model.

Our aim in this paper was not to develop a new yield criterion, but to investigate those previously applied to the prediction of ‘failure loads’. In this respect, the low value α=0.07 in Bessho et al. (2007) in the Drucker–Prager criterion resulted in a yield load which was very close to the one predicted by the von Mises criterion in all cases. This α is based on a correlation with yielding in concrete and should be calibrated to bone tissue specimens. Furthermore, we did not collect enough evidence to establish whether the strain-based criterion is the preferred one for situations involving compressive or shear stresses/strains.

The stress-based failure laws have two severe limitations that may bias considerably the yield load prediction in bones.

  • — The stresses are computed in the FE analysis by multiplying the strains by the material matrix, which is proportional to the Young modulus. On the bone’s outer surface, the Young modulus may be highly inaccurate if no special procedures are employed. This is because bone densities in CT voxels on the boundaries contain an average value of bone material and surrounding tissue; therefore, the estimated Young modulus may be inaccurate, and as a result the computed stresses may be inaccurate.

  • — The tension and compression yield stress is also given as a function of the bone density so, on the boundary, this value may also be highly inaccurate.

The strain criterion, on the other hand, is much more reliable because it is free of both limitations mentioned above. In summary, the use of the principal maximum strain criterion obtained by an orthotropic FE model averaged over the element surface is the best in predicting the yielding load and location and matches closely the experimental observations.

The main limitation of the present study is that only one femoral specimen was available for experimentation and no meaningful statistical analysis can be performed to substantiate the conclusions. Other loading configurations that induce high compressive and shear stresses and strains must also be considered in experiments to further validate the failure criteria in compression/shear. The exact ‘yielding load’ in such experiments is not well determined. We herein introduced a 95 per cent elastic slope intercept with the load–displacement curve as the ‘yielding load’.

Another limitation is associated with the relatively slow loading rate of 0.17 mm s−1 (168 μstrain s−1), which is one order of magnitude slower than the former reported physiological rates (2 mm s−1 for example in Cristofolini et al. (2007)). The slow strain rate may have been attributed to a more pronounced viscoelastic effect, enlarging the region beyond ‘yielding’. On the other hand, higher strain rates also necessitate a different relationship between the bone density and material properties when FE analyses are considered, and probably also different yield criteria which have to be strain-rate dependent. These topics have to be further investigated in future research activities.

It is likely that the ‘yield strain’ under tension (εyT=0.0073) in a purely cortical and trabecular bone are different. Thus using the single value for both bone tissues (especially being independent of the bone density) may introduce a large error.

Acknowledgements

The authors thank Dr Arie Bussiba & Mr Shimon Gruntman for their help with the bone experiment, and Prof. Charles Milgrom for supplying the bones and CT scans.

Footnotes

References

View Abstract