Finite element analysis of tibial fractures

Fractures of the tibial shaft are relatively common injuries, when e.g. front seat occupants are involved in frontal crashes [1]. Although the fracture mechanism is predominantly caused by axial compression of the leg, a substantial bending moment develops in the tibial shaft due to its natural curvature and to transverse loading when the tibial shaft collides with the dashboard [2]. This can lead to tibial shaft fractures [3]. In general, experimental crash tests have shown that tibial fractures may be expected when the sum of the normal stresses due to axial compression and the resultant bending moment exceed the strength of the tibial cortical bone (tibia index) [4]. In a review study of tibial shaft fractures, 88.9% were located in middle third or distal third part of the shaft [3]. These fractures were caused by a variety of occurrences including falls, various sports, pedestrian crashes, etc. All the studies referred to above gave an indication that tibial shaft fractures share a number of features in terms of site, type and local fracture mechanisms.

In this study, we aimed to set up an advanced, mathematical computer-based model using finite element (FE) analysis of the bones of the lower leg to examine whether an adequate computer model for predicting lower leg fracture initiation could be generated. We wanted to substantiate the validity of the model by predicting the fracture location to get an indication of the morphology of tibial shaft fractures for various common mechanical fracture loads (axial compression, torsion and bending).

MATERIAL AND METHODS

Finite element model

FE analyses were performed using the finite element code COSMOS/M. The three-dimensional "lower leg" FE model (LL FE model) used included the tibia and the fibula. The bony structures were generated by segmentation of a data set of computed tomographies from the Visible Human Project [5]. Segmentation is a data reduction tool used to simplify the geometrical shape of the bony structures of the structure without losing accuracy. Segmentation was necessary for the FE code to handle the large amounts of data involved. The cortical bone was simulated using tetrahedral elements (TETRA 4). The LL FE model comprised a total of 21,219 volumetric elements. The material property for trabecular bone was a uniform distribution of Young’s modulus of 100 MPa and a Poisson’s ratio of 0.2. Correspondingly, Young’s modulus for cortical bone was 12,000 MPa and a Poisson’s ratio of 0.3 [6]. The thicknesses of cortical and trabecular bone were taken from the Visible Human . The interosseos membrane and the tibiofibular joints were simulated using binding of the proximal and distal contact areas of the bones.

The fracture locations for this particular FE model were determined as the superficial areas of maximal Von Mises stress (σ Von Mises) in the cortical bone. Moreover, the σVon Mises were compared with a ratio of failure (RF) between σVon Mises and ultimate stress (σlim ) as in Scheilo et al [7]. The RF had to exceed 1 (RF > 1) for failure. The RF was determined as:

The σlim was determined in a linear relationship with the ash density in accordance with Keyak et al [8]:

Two different ash densities were determined from two different density-elasticity relationships between the apparent density by Carter and Hayes [9] and Scheilo et al [7]. Ash density was obtained from a ratio of 0.6 between ash and the apparent density, within the range 0.55-0.63 identified in the literature [8, 10]. The two density-elasticity relationships were specific for the hip region, but the results for both density-elasticity relationships with regard to the ultimate stress were similar (106 and 110 MPa, respectively). These results coincide with experimental results for ultimate stress for the tibial cortical bone, which is approximately 110 MPa for the tibia cortical bone with a yield stress amounting to 100 MPa [11].

Loading and boundary conditions

Three types of load were applied to the LL FE model:

1. Axial compression load. A total load of 800 N was distributed on the upper surface of the lower leg FE model. The compression load was applied as 60% to the medial side and 40% to the lateral side [12], as the knee contact force is not uniformly shared between the condyles of the tibia [13]. The lower end of the lower leg FE model was completely immobilized in all directions.

2. Torsional load. Torsional loads of 1,000 N in external rotation were applied to the upper and lower ends of the LL FE model. We were not able to obtain data with values for the torsional loads by which fractures to the tibia are known to occur, so the loads were assumed to be an adequate rotational moment (for the sake of convenience). The middle part of the tibial shaft was immobilized in the axial direction for analysis stability.

3. Transverse load.. Bending load was applied as nodal loads to the middle parts of the tibial and fibular shafts. A 5,000 N load was applied to the tibia and 2,500 N was applied to the fibula. The upper and lower ends of the LL FE model were completely immobilized in all directions.

Static analyses were performed for all three loading cases. Fracture locations and patterns were estimated by evaluating the maximal nodal Von Mises stresses as described previously. The LL FE model is illustrated in Figure 1.

Convergence test for the FE model was performed to ensure that the FE model had an appropriate number of elements. During this test, the number of elements in the FE model was increased until the point at which the calculated results converge to one exact solution, thus yielding the appropriate number of elements for in the FE model.

RESULTS

The distribution of Von Mises stress of the torsional load is illustrated in Figure 2.

Maximal Von Mises stress was located in the distal third of the lower leg, and Von Mises stresses were distributed in a spiral pattern. The RF was 0.46.

The distribution of Von Mises stress of the axial load is illustrated in Figure 3.

Maximal Von Mises stress was located in the middle third of the lower leg, and Von Mises stresses were distributed in an oblique pattern. The RF was 0.26.

The distribution of Von Mises stress of the lateral load is illustrated in Figure 4.

Maximal Von Mises stress was located in the middle third of the lower leg, and Von Mises stresses were distributed in a transverse pattern. The RF was 0.46.

Discussion

The purpose of the present study was to evaluate the use of the maximal Von Mises stress calculated by finite element analysis as a measure for examination of the fracture mechanisms of the lower leg. We included experimental material data and geometric data of the bones of the lower leg and were able to simulate adequate fracture localizations and to obtain an indication of fracture morphology of the most common lower leg fractures [3]. We applied adequate fracture loads with regard to direction and were able to generate fractures located in the middle and distal third of the lower leg, as expected. Moreover, the fracture patterns generated by the Von Mises stress distribution showed adequate spiral, transverse and oblique type patterns for torsional, lateral direct and axial loads, respectively. While static analysis may not be the best way to obtain true predictions of fracture patterns, the von Mises stress distribution method, on the other hand, yields clinically correct locations of fracture initiations prior to initial failure. After a fracture, stress redistribution at the fracture site would be expected, which would influence subsequent fracture patterns. Static FE analysis gives an estimate of the location of the fracture onset and the level of fracture risk. In this study, RF values for the various loads were somewhat lower than expected, which is presumably so because the estimated loads were too small. To the best of our knowledge, no in vivo data are available on this subject. However, certain in vitro mechanical test data on the whole tibia bone indicate that we have probably underestimated the loads – at least for axial compression (which had the smallest RF) – by a factor 4 [11]. Furthermore, studies have indicated that other mechanical signals, such as maximal principal strain or distorsion energy [7, 14], may be better suited for fracture analysis. This may also explain some of the discrepancy in RF values. However, the generated FE results were substantiated by known clinical and experimental studies with regard to locations and patterns [2, 3], thereby indicating that simulation of tibial shaft fractures is a valid path for further research, i.e. in fracture prevention, simulation of different types of osteosynthesis and the process of bony healing. We acknowledge the need for further improvement of the LL FE model as outlined above. To summarize, we have performed a general evaluation of the FE model and the computed results of a reaction force balance and the performed convergence tests, but have not correlated the computed results with an experimental cadaver test. The present LL FE model was based on geometrical data from one test subject (from The Visible Human Project) and we acknowledge that the computed stress results, and thus the fracture morphology and site, may be influenced by anatomical variation, i.e. the curvature of the tibia.

In this study, a linear analysis was chosen since the aim was not to investigate the whole fracture process, but to examine if we could reproduce mechanisms of fracture at onset. The use of static analysis is justifiable, because bone has been observed to have a quasi-brittle failure mode when tested experimentally [7]. Earlier studies have also indicated that even if bone is considered to be of anisotropic or orthotropic nature, it is justified to use isotropic material properties for maximal Von Mises stress analysis for long bones and also for uniaxial yield strain for the lower leg [15, 16]. However, the assumption of isotropic behaviour may have influenced the computed maximal Von Mises stress distribution, and thereby the fracture estimation. All in all, when considering the relative crudity of the LL FE model and the estimated loads, the use of isotropic material properties instead of anisotropy seems justifiable taking into account the lack of overall knowledge of anisotropic fracture properties of bone in FE models [14].

In this study, we excluded the effects of lower leg muscles. It has been demonstrated that the calf muscles influence the forward bowing (bending) of the tibia, thus affecting fracture development [17]. However, the muscles were not included, since – in the context of a fracture – the effects caused by muscles are currently unknown. The tibiofibular joints and the interosseos membrane were simulated by binding the tibiofibular joints. Again, appropriate mechanical properties have not been reported, hence justifying the crude simulation of these parts. In conclusion, we acknowledge that we would need to address the above mentioned considerations in future studies. However, in the present study we have been able to simulate fracture mechanisms, obtain adequate fracture locations and get an indication of fracture morphology. This indicates that the LL FE model may be used to examine e.g. fracture prevention and to simulate different types of osteosynthesis and the process of bony healing in future studies.

Correspondence: Christian Wong, Department of Orthopaedics, University Hospital of Hvidovre, 2650 Hvidovre, Denmark. E-mail: chwo123@gmail.com

Accepted: 24 February 2010

Conflicts of interest: None