## 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 *r*^{2}≈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*. 2007*a*,*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*. 2007*b*, 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.

— Inhomogeneous isotropic using empirical Hounsfield Units (HU)–Young’s modulus relations with a constant Poisson ratio according to Keyak & Falkinstein (2003).

— Inhomogeneous orthotropic based on microstructure with material principal directions coinciding with strain principal directions according to Hellmich

*et al*. (2008) and Trabelsi*et al*. (submitted).

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*. (2007*a*,*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):
2.1
2.2
2.3and
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*. (2007*b*) (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.

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.

### (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):
2.5
whereas the tensile yield stress is
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)
2.7
and
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 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,2.9

*Maximum principal stress.* The principal stresses are computed, and 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,2.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 by2.11where *I*_{1} denotes the first invariant of the stress tensor, and *J*_{2} is the second invariant of the deviatoric stress tensor,2.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 3*a*,*b*), (ii) averaged data in the element’s volume (see figure 3*c*,*d*), and (iii) the extremum value in the 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 by
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 or .

To determine the strength in tension or compression, the average density at the same *n* points was computed as follows:
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.

## 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, is the maximum principal stress criterion, and is the 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.

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

### (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.

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

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

↵1 StressCheck is trademark of Engineering Software Research and Development, Inc., St Louis, MO, USA.

One contribution of 13 to a Theme Issue ‘The virtual physiological human: computer simulation for integrative biomedicine I’.

- © 2010 The Royal Society