Introduction

Bone fracture often requires surgical treatments, especially in elderly population due to the aging process. Bone mechanical properties are substantively altered with age and the bone capacity to withstand the loading forces is decreased resulting in a worse bone resistance. One of the most interesting areas to focus on is the proximal femur due to the significant percentage of permanent disability and mortality after a hip fracture.5 Through numerical models based on the finite element method, it is possible to analyse the femur mechanical behaviour, accounting for the elastic range and fracture process. By this way, numerical models can help in the biomechanical and clinical fields to provide a deeper insight into the mechanical behaviour of bone.

A thorough understanding of bone fracture patterns in hip fracture is still a challenge and an active field of research. In this context, finite element models are a valuable tool for fracture behaviour prediction, although must be accompanied by experimental validation. The combination of experimental tests and numerical models is an important point of interest, because of their complementary character.8

In this work, we will use an artificial simulant of a human femur, made of composite material. The composite femur provided by Sawbones (Sawbones, Pacific Research Laboratories Inc., Vashon, USA) has been used by different researchers in order to analyse its behaviour and its performance. It is important to emphasise that this kind of specimens is designed to simulate the biomechanical properties of young healthy femurs2,9,13 providing advantages for model validation by avoiding the variability of properties inherent to biological tissues. These composite bones are useful to develop controlled analysis, thanks to their homogeneous properties, smoothed surface and non-variability between specimens.

In the literature, the composite femur has been analysed by different authors. In Ref. 25 a composite femur was used to measure the implant cup-bone relative displacements and stem-bone micro motion, measuring surface strains at different locations and orientations of the proximal femur. These measurements were compared with results from finite element (FE) models. In Ref. 26 the authors used the mechanical properties provided by the composite manufacturer.

A composite femur was used in Ref. 11 to predict the overall stiffness and peak bone stress in the same femur after injury, repair, and healing with respect to its intact condition. These authors concluded that the stiffness of a synthesized femur after using a plate and screws for repairing a 5 mm fracture gap was half that of the intact femur.

A femur analogue composite was also analysed in Ref. 17. In this case the finite element models were validated against strain data from digital image correlation (DIC). The predicted principal strains proved to be in good agreement with the experimental recordings from the DIC system, which provided a large number of measurements easily. The composite femur has also been used in other types of research, such as new positions of screws in displaced neck femur fractures27 or new infrared thermography technique to measure stresses in bone surface.30

In Ref. 7 the proximal femur zone was studied to identify its behaviour under different loading conditions and the corresponding fracture paths. This method demonstrated its ability to replicate spontaneous fractures in vitro. The work established that it is not essential to include the muscles when investigating head-neck fractures.

In Ref. 32 a femur micro-finite element model was developed distinguishing between cortical and trabecular bone. The stresses and strains were quantified for the loading case corresponding to a side fall onto the greater trochanter. The authors concluded that the cortical bone in the femoral neck is highly stressed in the osteoporotic case.

In Ref. 19 the human proximal femur behaviour was studied under several strain-rate loading conditions found in practice. Authors tested 12 cadaveric femurs up to failure and checked that in all specimens fracture started on the lateral side of the neck. They concluded that the proximal femur presents a linear elastic behaviour up to fracture.

A sideways loading configuration was studied in Ref. 15 to evaluate a subject-specific finite element modelling technique from CT data. Numerical and experimental studies were performed to compare finite element predicted strains with strain-gauge measurements from three cadaveric proximal femurs. Twelve loading configurations were carried out, spanning a wide cone.

In Ref. 1 five specimen-specific femur models were developed obtaining the geometry and materials from a CT scan. The specimens were analysed in an elastic range and the fracture load in experimental tests was compared vs. numerical models. In this case small fracture paths were obtained, due to the convergence problems that appear in certain implementation of the XFEM. Other authors have employed the material property degradation technique in their numerical models, leading to longer fracture paths.18

Biomechanical properties of the bone obtained from other works in the literature are commonly used when modelling proximal femur fracture. However it is worth noting that accurate predictions strongly depend on a truthful bone behaviour characterization. There is a wide dispersion of numerical values of mechanical properties of bone in the literature due to changes in terms of the ageing, disease, nutrition and other factors.10,14,24 On the other hand, when weakening effects of osteoporosis exist, the importance of characterizing the bone behaviour increases, because of the enlarged scattering of the properties. These facts indicate the importance of accounting for the actual mechanical behaviour of bone, in order to obtain accurate numerical predictions.

The main goal of this research is to describe a combined methodology (experimental and numerical) for the bone fracture analysis through the use of a composite bone. The work provides the procedures followed to develop consistent FE models that can be used to ascertain different fracture conditions. This work includes two types of experimental tests: at local level for characterization of the mechanical properties of the bone and at global level for fracture testing of the proximal femur. A numerical model is developed that implements the mechanical properties previously obtained from the characterization tests as input properties.

In the global numerical model, the crack onset in the fracture process is modelled using an XFEM approach and the crack initiation is predicted using a maximum principal strain criterion. Although there are studies that use other fracture criteria, this criterion is one of the most used in the literature. Many works use FE models to evaluate fracture patterns and loads making use of a fracture criterion. A complete review of these approaches is given in Doblaré et al.10 There is a lack of agreement about the use of strain-based or stress-based failure theories, although many authors prefer the first one, probably due to the quasi-brittle behaviour of bone.10

In real human femurs, fracture behaviour is essentially a tensile failure mode for a stance loading and compressive failure mode for a fall configuration.28 Therefore, either a maximum or a minimum stress/strain criterion is used to localise the fracture point in the numerical models, with different fracture limits for tension and compression.3

In artificial femurs similar to the one studied in this work, Ebrahimi et al. analyse the fracture load by means of FE.11 They use a von Mises equivalent stress criterion to find the location of the peak stress.

Usually, a maximum principal strain criterion is proposed in the literature to estimate the failure risk in human femur models. Schileo et al.29 concluded that the maximum principal strain criterion can be defined as a suitable method to evaluate fracture risk in human femur. Other authors follow this approach, see e.g.,1,4,29 In some of these numerical models, tension or compression failures are differentiated.16,31 Other failure criteria include an equivalent strain criterion to model fracture paths18 and even more complex criteria, such as a simplified Tsai-Wu criterion to predict failure of trabecular bovine femur.12

In this work, only crack initiation and the early stages of crack propagation are considered due to the difficulties associated with the numerical modelling of long propagating cracks in bone fracture. The actual path during the last stages of crack propagation depends on the presence of local defects and heterogeneities which cannot be represented in the models here presented. Also, the load-controlled nature of the problem leads to unstable crack propagation. Our main interest focuses on the correct prediction of crack location during the first stages of growth. This may provide insight about the origin of fracture when dealing with real bones.

Materials and Methods

Experimental Work

The experimental work was carried out using a bone analogue made of composite material presenting a mechanical behaviour similar to healthy human femur.2,26 A fourth generation composite femur (Model no. 3406, Sawbones) was used. This model is made of two distinct materials, see Fig. 1a. It has an exterior zone which simulates cortical bone (with variable thickness) and interior foam with similar properties to trabecular bone.

Figure 1
figure 1

(a) Composite femur provided by Sawbones; (b) CT scanner from the femur proximal zone; (c) Bone geometry obtained from the CT scanner.

Full size image

As commented above, experiments include the mechanical characterization of the composite bone and fracture testing. In addition, a bone density scanner (CT) of the composite femur and image processing was carried out prior to model development and mechanical testing. The experimental work is described below.

Scanner and Geometry Generation

The development of the numerical model requires the introduction of the real bone geometry in the FE code together with its density distribution. This fact is especially important when human bones are studied due to the great dispersion in geometry and mechanical properties due to age, gender, size and other factors. In this work a composite bone is used and a CT scan is proposed in order to illustrate the methodology directly applicable to human bone.

The composite femur was introduced in a CT-scanner (SIEMENS Somaton) with a pixel size of 0.44 mm and a slice thickness of 1 mm. The scanning parameters allowed the differentiation of the materials that mimic the cortical and trabecular zones (see Fig. 1b). These areas required a special image treatment due to the small cortical thickness. The small cortical thickness (especially in the upper side of the femoral neck) made the model meshing operation difficult, as explained below.

The image treatment, the bone modelling and the numerical model meshing were carried out through the software ScanIP (Simpleware, Exeter, UK) from a CT scanner (Fig. 1b). A series of filters were carefully applied to generate the correct geometry, paying special attention to the femoral neck, due to small thickness of the cortical bone. The two different volumes (cortical and trabecular bone) were clearly differentiated by means of the application of several filters (having into account the differences in Hounsfield Units values of the two materials due to their different densities). The zones were separated from the rest depending on their grey scale value for each voxel, which is related to Hounsfield Units. Finally Gauss filters were applied to smooth the bone surface, obtaining the final geometry. Material properties were input into Abaqus FE code, with constant Young’s modulus for each material that simulate mimics the cortical and trabecular zones, since they are considered homogeneous in this work.

Finally the proximal 3D geometry was obtained from the CT scanner (Fig. 1c), with a surface mesh consisting of about 650,000 triangles tiling the surface. Through this software a FE mesh with refined zones was generated, as explained in the numerical model section of this paper.

Mechanical Characterization

As commented above, the mechanical characterization is essential since the material modelling clearly influences the numerical results obtained from the FE model. It is important to carry out mechanical tests on the same bone specimen for validation, because there is a great scatter of the properties found in the scientific literature.

A three-point bending test was used to obtain the Young’s modulus and the cortical bone critical strain, as detailed in what follows.

Three-Point Bending Test for E Calculation

Although the manufacturer of the composite bone provided material properties that are typical of bone, several experimental tests were performed to obtain the Young’s modulus of the composite material that is used to replicate the cortical bone.

Four specimens were prepared from the cortical zone in the femur diaphysis to carry out three-point bending tests (Fig. 2a). These were obtained after performing the femur fracture test, since the diaphysis was not damaged during breakage. All samples were rectangular prism shaped; although, due to the irregular geometry of the femur they presented different dimensions as detailed in Table 1.

Figure 2
figure 2

(a) Specimens used in the three point bending tests to estimate E; (b) three point bending test; (c) notched specimen used in the tests for indirect measurement of critical strain; (d) three point bending test with notched specimen; (e) dimensions of the notched specimen; (f) specimen 6 after failure.

Full size image
Table 1 Dimensions of the intact specimens.
Full size table

Tests have been carried out using a testing machine (MTS Criterion C42.503, load cell 100 N). An increasing load was applied on the specimen in the mid-section with an actuator speed of 0.05 mm/s, what can be assumed equivalent to quasi-static conditions (Fig. 2b). Assuming an isotropic and elastic behaviour, the Young’s modulus E can be obtained from the beam deflection d and other specimen properties through the expression:

$$E = \frac{{F \cdot Z^{3} }}{48 \cdot d \cdot I}$$
(1)

where F is the load, Z is the length between supports and I is the section moment of inertia.

Progressive load was applied up to 90 N and during this load application five measurements were obtained at the load application point. By averaging all measurements with the four specimens (see Table 1) the estimated Young’s modulus is 10404 MPa, with a standard deviation of 1071 MPa.

Three Point Bending Test to Obtain Critical Strain

The critical strain governing the onset of failure of the bone is an important parameter for the modelling of crack pattern and fracture using the XFEM method. Moreover in the literature there is a great dispersion of the values for critical strain of cortical bone. Up to the authors’ knowledge, there is no reported value for the composite bone used in the experiments of this work.

Here the maximum principal strain has been considered as the critical strain to initiate the fracture onset. This criterion has been commonly employed by other authors1 based on the hypothesis that the fracture starts in the upper region of the femoral neck, which is loaded in a tensile state. The critical strain was obtained from specimens including a notch in the middle section (see Figs. 2c and 2d) subjected to three point bending tests on the same machine. Six specimens were extracted from the diaphysis after the fracture test, all from the cortical zone since the numerical model has only into account the fracture at its first stages. The small semi-cylindrical notch was machined in the specimen central zone in order to generate a stress concentration leading to failure in this area (Fig. 2c). Therefore the crack always starts at the central region, thus having a controlled failure location.

The specimen dimensions are shown in Table 2, while a sketch of the specimen is given in Fig. 2e.

Table 2 Dimensions and critical strain (ε c) of the specimens with notch.
Full size table

The load was progressively applied up to failure (with an actuator speed of 0.1 mm/s, i.e., quasi-static conditions). The crack onset occurred at the notch as expected (Fig. 2f). Once the maximum load leading to failure was obtained, a finite element analysis for each of the six specimens was modelled in order to inversely obtain the critical failure strain by inverse correlation. This analysis will be explained in detail in the section dealing with numerical modelling. The values calculated for the critical failure strain are included in Table 2.

Femur Mechanical Tests

Two different tests were carried out with the complete composite femur. First, the global stiffness of the femur was initially checked. Next, a fracture test that involved a load increase up to bone failure. The aim of these tests is to obtain experimental results for the numerical model validation. The stiffness test enables the measurement of strains at different points of the femur during loading, whereas the fracture tests provide the crack location and maximum failure load.

Stiffness Test

The composite femur was loaded in elastic regime in a 100 kN universal hydraulic testing machine (INSTRON 8801, load cell 100 kN). Strain gages were adhered to the femur in order to measure strain during the test. Eight strain gages measurements were used to study the strain in a local direction. The gages were stuck on the top and middle of the diaphysis (4 gages) and on the femoral neck (4 gages, see detailed location of the strain gages in Fig. 3a).

Figure 3
figure 3

(a) Location of strain gages in the femur; (b) stance loading configuration; (c) experimental rig for stiffness evaluation with strain gages; (d) detail of the experimental rig; (e) broken specimen after fracture test.

Full size image

The loads were applied to the femur in the stance loading direction8 (Fig. 3b). This position is the common configuration analysed in the literature because it corresponds to the physiological scheme of a standing-up human position.

The bone was aligned by rotating the long axis of the femur to 8° adduction in the frontal plane by means of an appropriate rig support (see Fig. 3d). The load was applied on the femoral head, on a spherical region of diameter 35 mm. The distal diaphysis was potted in surgical cement (Surgical Simplex P, STRYKER, Mahwah, NJ, USA), with a distance between the load and the cement surface of 250 mm. Three different loads were considered to analyse the femur deformation in an elastic range. These loads were 250, 500, and 750 N. The data acquisition was carried out at once when the load was reached in order to avoid relaxation of the material.

Load was increased up to the different values considered (250 N, 500 N and 750 N with an actuator speed of 0.3 mm/s in order to keep quasi-static conditions). The femur linear elastic behaviour was verified through this test, as reported by other authors.8,19 The strains obtained from the gages experimentally and strains from the numerical model are compared and discussed in the results section. Figures 3c and 3d show the complete experimental rig, including measurement channels.

Fracture Test

The femur did not undergo any damage in the previous test since it was loaded in the elastic range, thus the fracture test was carried out in the same sample. The load was applied in the same configuration (8° from vertical axis) with an actuator speed of 0.3 mm/s.

From this test, the maximum load and fracture pattern were obtained. The failure occurred in the upper zone of the femoral neck. The maximum load registered at fracture was 6330 N. The results from the experimental test and the comparison with the numerical model predictions are shown in the results section. The fracture path caused in this test is shown in Fig. 3e.

Numerical Modelling

Two different models are used in the present work. The first one simulates the three-point bending tests of the notched specimens, and it was developed in order to estimate the critical strain. The second model simulates the behaviour of the whole femur, both in elastic regime and in the fracture test. These models are described in the following sub-sections.

Numerical Modelling of Notched Three-Point Bending Specimens

Model Description

Six models corresponding to the dimensions of the notched specimens given in Table 2 were developed in Abaqus 6.12 in order to infer the cortical critical strain. A two dimensional approach was used since the thickness is enough to assume a plane strain condition. The Young’s modulus obtained from the unnotched three-point bending specimens was input into this model in order to obtain the most realistic value for critical strain.

The numerical model was loaded to the registered experimental failure load. It was developed using the commercial FE code Abaqus/Standard, using two dimensional quadrilateral 8-node elements with 2nd order interpolation under plane strain (CPE8R in Abaqus). These 2nd order elements combined with a reduced integration provide enough accuracy for the FE analysis. A previous mesh sensitivity analysis was carried out to guarantee convergence, and therefore an accurate estimation of the principal strain at the notch. This analysis was performed for the numerical model corresponding to specimen number 1, varying the element size from 1 to 0.1 mm. Finally a size element of 0.25 mm was chosen, since the variation in the estimated critical strain was negligible for smaller elements sizes next to the notch zone.

The model was composed of about 4000 elements and 12,000 nodes. The boundary conditions reproduced the three-point bending test as shown in Fig. 4.

Figure 4
figure 4

Model of the notched specimen No. 1 and detail of the notch area showing maximum level of strain.

Full size image
Numerical Results: Critical Strain

The 2D model allowed the calculation of the strain field for each notched specimen. The load applied in the simulations was the maximum load producing failure for each case. The maximum principal strains in the vicinity of the notch were estimated from the numerical model (see Fig. 4 for the notched-specimen No. 1).

The analysis of the five cases yielded a mean value of critical failure strain equal to 0.0165 (values for the six samples are shown in Table 2) with a standard deviation equal to 0.0014. The value obtained in the experiments for the cortical material of the composite bone is higher than the critical strain given in Ref. 25, where a critical strain of 0.006 for human cortical bone is reported. The critical strain ε c is an essential quantity of interest that is required as an input in the following macro model to predict the onset of crack initiation under an XFEM approach.

Specimen deflection at fracture onset has been measured in the six notched three-point bending specimens (both from experimental tests and numerical models). The good agreement between numerical and experimental results validates the procedure for the Young’s modulus estimation (E = 10,400 MPa) and the numerical models used for critical strain calculation.

Numerical Femur Model

Model Description

The development of the femur model requires the previous scanner of the specimen, the treatment of the medical images and the generation of the mesh. The software Scan IP was used to create the 3D femur geometry. The software allowed the differentiation between cortical and trabecular zones and the assignment of different mechanical properties depending on local density.

A mesh sensitivity study was developed prior to define the suitable element size in the femur. This step avoids an excessive computational cost and, at the same time, guarantees an adequate accuracy in the predicted parameters.

The element size was set equal to 4.3 mm in the cortical diaphysis zone and a refined element size of 2 mm in the cortical proximal zone (a similar element size was reported in Ref. 1). The neck zone undergoes elevated stresses and usually experiences the onset of fracture. Therefore, a refined mesh at this zone with an element size equal to 1 mm is required to achieve an accurate solution in the fracture area. The trabecular zone was meshed with an element size of 3 mm since it has less influence on the global stiffness of the femur.

Altogether the femur was meshed with 184,400 quadratic tetrahedral elements (type C3D10 in Abaqus) and 295922 nodes. Figure 5a shows the mesh with the refined zones.

Figure 5
figure 5

(a) Femur meshing showing the refined zone at the neck; (b) principal strains in femur, under stance loading (750 N); (c) fracture onset predicted with the numerical model for stance loading.

Full size image

The boundary conditions and the loading state are set as close as possible to the conditions imposed during the experiments. The distal part of the diaphysis is constrained to the testing rig simulating the embedment in surgical cement and the load was applied in the proper direction similar to the configuration used in the experiments. The testing rig was also included in the numerical model as it was checked that its stiffness has an influence about 10% on the strains on the bone surface (this is in agreement with the numerical models developed by other authors15).

Mechanical properties of the bone used in the simulations are summarized in Table 3. As commented above, the Young’s modulus and the failure strain were obtained from the experiments developed in the present work and inferred from the three-point bending tests. The properties related to fracture energy that must be input to the XFEM numerical model were obtained from Refs. 6 and 33. Critical energy values for different fracture modes have been obtained from the cortical fracture toughness K C, which is related to bone density through Eq. (2), given in Ref. 6. Geometric nonlinearity was considered in the FE analysis when fracture process was modelled. The following expressions determine these relationships, and their values are shown in Table 3.

Table 3 Composite femur mechanical properties.
Full size table
$$K\left( {N m^{ - 1.5} } \right) = 0.7413 \cdot 10^{6} \cdot \rho^{1.49}$$
(2)
$$G\left( {J m^{ - 2} } \right) = \frac{{K^{2} \left( {1 - \nu^{2} } \right)}}{E}$$
(3)
$$G_{\text{IIC}} /G_{\text{IC}} = G_{\text{IIIC}} /G_{\text{IC}} = 0.33$$
(4)

The crack onset was predicted through the initiation criterion available in the code Abaqus based on the maximum principal strain. The propagation was simulated using the XFEM capability available in Abaqus, using the Virtual Crack Closure Technique (VCCT) with mixed mode behaviour based on the Benzeggagh-Kenane expression. The XFEM method enables the introduction of crack surfaces that are independent of the mesh geometry (no need to conform to element sides) which is a great advantage for crack modelling using the finite element method.

Results and Discussion

Numerical Results: Elastic Regime and Failure

The loading state described in the experimental work was reproduced by the numerical macro model of the femur, as shown in Fig. 5b, showing the principal strain contour under stance loading of 750 N. The FE strains were obtained at the same points where the strain gages were positioned.

Numerical model predictions were correlated (R 2 = 0.99) to the experimental values with a slope of the regression line of 0.81 for 24 strains (eight strain gages for each load case) at different loads (250, 500 and 750 N), see Fig. 6. Average relative errors between model and experimental strains were about 9%.

Figure 6
figure 6

Comparison between the numerical model and experiment (microstrains). The numbers next to the markers indicate the gage position number given in Fig. 3a.

Full size image

The mean error for the predicted strains (about 9%), is deemed reasonable when compared with other results in the literature (see for instance1). Several factors lead to differences between measured and predicted values. One source of error is the position of the strain gage that can slightly differ from the zone where the numerical strain is registered. Also the boundary conditions imposed in the numerical model induce some errors due to the compliance of the attachment at the base of the bone.

The review of works developed by other authors with composite femurs evidences a dispersion of the results in the elastic regime. Only Ebrahimi et al.11 get a slope close to 1.0. Grassi et al.16 obtain similar results to the ones shown in this work, with a slope about 0.8, underestimating predicted strains, whereas Pal et al.26 overestimated their numerical strains, with percentage errors about 15%.

Concerning the fracture loading, the experimental test yielded a maximum load of 6330 N. The fracture test occurred at the neck, as it is also found in the works by other authors.1,7 Our numerical simulations also predict the failure at this zone; see the fracture path calculated for the stance loading case in Fig. 7a.

Figure 7
figure 7

Boundary conditions, fracture paths obtained with the numerical model and details with principal strains in different configurations: (a) stance loading7; (b) trochanter fall32; (c) sideways fall.22

Full size image

The prediction of the fracture load requires the statement of a failure criterion. In this case, fracture load has been calculated when a node reaches critical strain, obtaining a predicted failure load of 6069 N in the numerical model, which is close to the load of 6330 N found in the experimental test.

Regarding the predicted fracture location and early stages of crack growth shown in Fig. 5c, it can be seen that it is similar to the failure region observed in Fig. 3e. Moreover, the numerical results showed the same trends obtained by other authors also with composite femurs.17

To these authors’ knowledge, there are no works in the literature dealing with failure in composite femur using finite element models. Some researchers have analysed experimentally the fracture load of composite bones in order to compare it to real femurs,13 concluding that the failure mode is similar to published results for human femur. Moreover, fracture load for a stance loading test is similar for composite and real femurs. In this configuration, fracture load in a human femur is about 6000 N in experimental results obtained by Refs. 1 and 28, although obviously real femurs may show a high dispersion.

The use of XFEM to model fracture of a human femur has been accomplished by Ref. 1 with results similar to the ones presented in this work. Like in this study, Ali et al. only simulated the initial steps of femur fracture.

Other Load Configurations

In addition to the stance loading scenario used for model validation, other interesting loading conditions have been studied in this work in order to analyse different fracture paths. Falling to the side has been analysed by several authors with different loading conditions. In this work loading configurations denoted “sideways fall” (introduced in Ref. 22) and “fall to the side or trochanter fall” (analysed in Ref. 32) have been studied using the numerical model previously developed. The different configurations proposed in Refs. 32 and 22 have been applied in our composite femur and the fracture load has been calculated with the numerical model. Boundary conditions for different fall configurations are shown in Fig. 7. The contact against the impacting surface is modelled by constraining the normal displacement to this surface in a certain number of nodes of the trochanter region. In the trochanter fall loading the angle between the femoral shaft and the horizontal is 10° and the femur is internally rotated by 15°. The sideways fall loading condition consists in a load applied to the femoral head at 20° anteversion and 30° rotation. Values for the predicted fracture load are detailed in Table 4 using the value of critical strain at one node as failure criterion. The predicted fracture paths using XFEM in our model are shown in Fig. 7.

Table 4 Fracture loads for different loading conditions.
Full size table

As concluded also by other authors,19 the numerical model predicts the fracture initiation on the lateral side of the femoral neck. Concerning other loading configurations, the sideways fall causes a fracture closer to the trochanteric region. Results for the fracture load show that the sideways fall is the most critical loading configuration. Other authors have analysed these fall configurations in human femurs, but not in composite femurs. Experimentally, and restricting to human femurs, most authors have obtained fracture loads that are lower for fall conditions22 than for the stance fracture loads. This is also the case in the numerical simulations by Ref. 1.

Despite the numerical model has been validated through the stance configuration, further experimental tests involving other configurations should be carried out in the future. These tests will provide experimental confirmation in order to assess which loading configurations are the most critical, leading to significant conclusions for clinical applications.

Influence of Material Properties on Fracture Load

As explained above one of the contributions of this work is to highlight the importance of the estimation of proper mechanical properties to be input in the FE numerical model. In this section the influence of the Young’s modulus (both at the cortical and trabecular zone) and the critical failure strain was analysed through a parametric study. Different values of these parameters obtained from real bones in the literature have been input in the numerical model, leading to strong differences in the estimated load.

Table 5 shows the fracture load obtained for the nominal values of mechanical properties obtained in the present work, and the variations in the numerically predicted fracture load when these properties are varied in the range given in the literature. The fracture load has been obtained for the configuration corresponding to stance loading. The predicted fracture load changes dramatically, up to about 100%, simply by changing each property separately in the range provided in the literature (see Fig. 8). This result highlights the importance of a proper mechanical characterization of the bone used in the tests for the accurate simulation of fracture initiation.

Table 5 Predicted fracture load obtained considering the variation of mechanical properties.
Full size table
Figure 8
figure 8

Variations of fracture load (stance loading configuration) with mechanical properties of real bones found in the literature.

Full size image

Figure 8 shows the fracture load dependency when the most important parameters are varied in the numerical model. A critical principal strain variation produces the highest fracture load variation (its slope is the greatest). On the other hand, elastic modulus variation of the trabecular bone has little effect on the predicted fracture load. Obviously when the Young’s modulus or the critical principal strain increases, the fracture load increases too. Similar studies analysing the sensitivity of femur fracture load to these parameters have not been found in the literature. By means of this study, the influence of some important material properties has been assessed, being one of the objectives of this work.

Of course, the modelling of real bones should include the non-isotropic behaviour and heterogeneity found in real bone tissues, in contrast to the isotropic and homogeneous character of bone simulant, at least at each of its layers. However, this is not a shortcoming of the proposed procedure. The FE models can include non-isotropic material properties and a local definition of material properties at the element level, provided the local properties are characterized, which is probably the most difficult task. These aspects can have an important influence on the computation of the actual crack location and propagation.

Conclusions

The main objective of this work is to present a methodology for the analysis of femur fracture, including XFEM modelling and experiments involving both bone characterization and validation of fracture simulation.

The experimental tests enabled the mechanical characterization of the proximal femur. Three-point bending tests were performed in order to obtain the cortical bone Young’s modulus. The critical principal strain corresponding to failure was obtained from notched specimens under a three-point bending test and using numerical modelling to infer the property through inverse analysis. The methodology was applied to a composite bone, establishing a procedure applicable to a human bone.

The main conclusions derived from this work are summarized below:

  • The experimental tests proposed to measure bone properties demonstrated their ability to obtain accurate values of Young’s modulus and critical strain, at least for the simulant specimen considered in this work. The calculation of critical strain required also the development of a numerical model simulating the tests using notched specimens and estimating ε c through inverse analysis.

  • The numerical model of the femur based on a FE approach using XFEM provided good predictions both for elastic behaviour and for the fracture crack path under stance loading. It was highlighted the importance of an adequate calibration of the fracture criterion.

  • The influence of mechanical properties (Young’s modulus and critical principal strain), was evidenced by a sensitivity analysis simulating fracture with the XFEM model for a wide range of real bone values commonly used in the literature. The importance of a proper characterization of bone is proven by the great differences found in the predicted fracture load when using a range of values from the literature.

  • Once calibrated, the model was used to predict fracture paths in other loading conditions. The most unfavourable case is the load configuration corresponding to a sideways fall. This leads to the lowest fracture load.

Our study presents the limitation that only one artificial composite femur has been used and this sample fits better young bones. However, the principal goal of this work is to set a methodology that can be applied to human femurs, being the use of composite bone appropriate for this objective, as validated by the results.

The analyses presented in this work are focused on predicting the onset of fracture and the early stages of crack propagation. The complete modelling of bone fracture is a difficult task, since the actual crack path for the last stages of propagation may depend on many local factors (local defects and heterogeneities). Further improvement of the present work should involve the application of the proposed methodology to human femurs, including the material micro-characterization of the bone tissues in different regions. When human femurs are studied the non-isotropic behaviour and heterogeneity should be taken into account, both in experimental characterization and numerical modelling. In this latter case, a local definition of material properties at element level is necessary.