An Idealised Biphasic Poroelastic Finite Element Model of a Tibial Fracture

The outcome of a bone fracture partly depends upon the mechanical environment experienced by the fracture callus (reparative tissue) during the healing. Therefore biomechanics of bone fracture healing has been examined in many clinical or biological, mathematical or finite element studies (Cheal et al. 1991, DiGioia et al. 1986, Claes et al. 1999, Doblare et al. 2004 and Oh et al. 2010). Most of the studies model the components of bone fractures as monophasic, homogenous materials, which may not be appropriate considering the large inter fragmentary displacements and high porosity of the reparative tissue. Therefore, this study describes an idealised mathematical model of a healing bone fracture with biphasic approach when the callus bone is modelled as mixture of solids and fluids. Markel et al. (1990) reported that the porosity of the callus in a healing canine osteotomy decreased from 99.6% at 2 weeks to 38% at 12 weeks. Therefore, the biphasic, poroelastic model for fracture callus and bone has been suggested in the literature (Carter et al. 1998, Simon et al. 1992, Prendergast et al1997, Spilker et al. 1990). Biphasic poroelastic models for soft tissues (Mow et al. 1980, Simon et al. 1985, Van Driel et al. 1998, Prendergast et al. 1997, Spilker et al. 1990) have been developed and applied to model cartilage (Mow et al. 1980) and intervertebral discs (Simon et al. 1985). Van Driel et al. (1998) and Prendergast et al. (1997) modelled tissue adjacent to prostheses using poroelastic material properties to investigate tissue differentiation. In the field of fracture healing however, only monophasic material properties of callus have been simulated (Carter 1988, Carter 1998, Blenman 1989, Cheal 1991, DiGioia 1986, Claes 1999, Gardner 1998 and 2000). This is probably because of the paucity of data in the literature on the values of parameters required to define the biphasic material properties of fracture callus. Simulation of a biphasic, compressible, anisotropic, linear poroelastic material model requires forty material constants (Simon 1992). Even the very simplified simulation of an isotropic material requires a minimum of five material constants. However, the number of material constants required to simulate a biphasic, poroelastic medium can be further reduced to three if the solid and fluid media are assumed to be incompressible (Simon 1992). These three independent material parameters are Lame's material stiffness parameters (λ and μ) and hydraulic permeability (k). Alternatively, Zienkiewicz and Taylor (1994b) suggested a method to model poroelastic behaviour under `undrained' condition using the modulus of elasticity, Poisson's ratio, the


Introduction
The outcome of a bone fracture partly depends upon the mechanical environment experienced by the fracture callus (reparative tissue) during the healing. Therefore biomechanics of bone fracture healing has been examined in many clinical or biological, mathematical or finite element studies (Cheal et al. 1991, DiGioia et al. 1986, Claes et al. 1999, Doblaré et al. 2004and Oh et al. 2010. Most of the studies model the components of bone fractures as monophasic, homogenous materials, which may not be appropriate considering the large inter fragmentary displacements and high porosity of the reparative tissue. Therefore, this study describes an idealised mathematical model of a healing bone fracture with biphasic approach when the callus bone is modelled as mixture of solids and fluids. Markel et al. (1990) reported that the porosity of the callus in a healing canine osteotomy decreased from 99.6% at 2 weeks to 38% at 12 weeks. Therefore, the biphasic, poroelastic model for fracture callus and bone has been suggested in the literature (Carter et al. 1998, Simon et al. 1992, Prendergast et al1997, Spilker et al. 1990). Biphasic poroelastic models for soft tissues (Mow et al. 1980, Simon et al. 1985, Van Driel et al. 1998, Prendergast et al. 1997, Spilker et al. 1990) have been developed and applied to model cartilage (Mow et al. 1980) and intervertebral discs (Simon et al. 1985). Van Driel et al. (1998) and Prendergast et al. (1997) modelled tissue adjacent to prostheses using poroelastic material properties to investigate tissue differentiation. In the field of fracture healing however, only monophasic material properties of callus have been simulated , Carter 1998, Blenman 1989, Cheal 1991, DiGioia 1986, Claes 1999and 2000. This is probably because of the paucity of data in the literature on the values of parameters required to define the biphasic material properties of fracture callus. Simulation of a biphasic, compressible, anisotropic, linear poroelastic material model requires forty material constants (Simon 1992). Even the very simplified simulation of an isotropic material requires a minimum of five material constants. However, the number of material constants required to simulate a biphasic, poroelastic medium can be further reduced to three if the solid and fluid media are assumed to be incompressible (Simon 1992). These three independent material parameters are Lame's material stiffness parameters (λ and µ) and hydraulic permeability (k).
Alternatively, Zienkiewicz and Taylor (1994b) suggested a method to model poroelastic behaviour under `undrained' condition using the modulus of elasticity, Poisson's ratio, the porosity of the matrix, and the bulk modulus of the fluid phase. In the present study, a finite element model (FEM) based on the poroelastic behaviour of the `undrained' callus at four temporal stages of healing is developed by modifying the theory proposed by Zienkiewicz and Taylor (1994). This model was developed to examine the influence of fluid pressure on the pattern of healing and to compare the distribution of stresses in the callus with the monophasic solutions developed for the same subject at the same temporal points reported earlier by Gardner et al. (2000).

Materials and methods
2D, monophasic, plane stress, FEM's of a mid-diaphyseal tibial fracture were developed at four stages during healing (4, 8, 12, and 16 weeks post operation) by Gardner et al. (2000). The geometry and the regionalisation of the callus are shown in Figure 1. The geometry, finite element meshing, boundary conditions ( Figure 2) and applied displacements (Table 1) used in the study of Gardner et al. (2000) are adopted in the biphasic poroelastic models of the present study. The tissue histology and calculated elastic moduli of regions of callus at four stages of healing are shown in Table 2.   However, the callus tissue in the present study was idealised as a homogenous, fully saturated, linear poroelastic medium consisting of a matrix of solid and incompressible fluid, as opposed to the monophasic material properties in the study of Gardner et al. (2000). The following section describes the theory employed for calculating equivalent poroelastic material properties of the callus from the modulus of elasticity, Poisson's ratio, porosity and bulk modulus of the fluid of the callus tissue. Unlike a monophasic medium, the normal stress acting across a plane within a biphasic poroelastic mass will have two components, an inter-granular pressure known as effective pressure or effective stress, and a fluid pressure called the pore pressure. The sum of these two will constitute the total normal stress. The volume change characteristics and the strength of poroelastic mediums are controlled by the effective stress not by the total stress. Thus the only difference between the present study and the previous study of Gardner et al. (2000) is that the constitutive equation of the callus is changed from monophasic medium to a biphasic poroelastic medium.

Mathematical description of the model
For 2D plane stress analysis, the constitutive equation for a monophasic material (Zienkiewicz and Taylor (1994a)) is Where σ is total stress. σ x , σ y are normal stresses in the x and y directions, τ xy is shear stress in the x-y plane, ε x , ε y are normal strains in the x and y directions and γ xy is the shear strain in the x-y plane. For plane stress analysis the constitutive matrix D (Zienkiewicz and Taylor 1994a) is defined as Thus a linear elastic plane strain D-matrix is obtained. The equation for static equilibrium in 2D (Dawe 1984) is: ( ∂σ x /∂x ) + (∂τ xy /∂y ) + R x = 0 (5a) ( ∂σ y /∂y ) + (∂τ xy /∂x ) + R y = 0 (5b) Using the definition of effective stress (Zienkiewicz and Taylor (1994b)) where σ X , σ Y, τ XY are total stresses; R x , R y are body forces; σ' X , σ' Y , τ' XY are effective stresses (positive value for tensile stresses and negative for compressive stresses), p is pore pressure of the fluid, which is conventionally described as positive for compressive pressure and negative for tensile pressure.
The combined seepage and conservation of fluid equation described by Zienkiewicz and Taylor (1994b) is: Where k is the permeability of the fluid, Q is the ratio of the bulk modulus of fluid to the porosity of the media and ε X, ε Y represent the volumetric strain rates of the solid skeleton. The superscript 'dot' denotes differentiation with respect to time. Since the gait cycle frequency of this clinical fracture was approximately 1 Hz, it is reasonable to assume that little or no seepage from the callus occurs during loading (Carter 1998. Therefore, because the time scale is short, if the local `undrained' condition of the callus is assumed, permeability 'k' can also be assumed to be effectively zero. Under these conditions, Equation 7 becomes: Integrating with respect to time and assuming homogenous initial conditions (p = 0, Substituting this value in equation (6a) and (6b) gives: www.intechopen.com Expanding for two-dimensional plane strain analysis (σ z ≠ 0) gives: Noting that the material behaviour is controlled by the effective stress in a poroelastic medium, then Combining equation (13) and (14): If the modified stress-strain relationship is defined as: and where K f is the bulk modulus of the fluid phase and n is the porosity of the porous media. Thus the D' matrix (Equation 17) for a poroelastic medium can be calculated by substituting Q values in the D matrix of corresponding plane strain analysis (Equation 4C).
The D matrix can be expanded to calculate the new set of material properties, E and ν, corresponding to the poroelastic material behaviour under undrained condition.

Development of the poroelastic FEM from the monophasic FEM
Using the above theory, a new set of material properties (E and ν) were calculated from the refined values of E and ν of the monophasic model of Gardner et al. (2000) and the values of Q described in this section. At first, the proportion of calcified tissue in the callus was extrapolated from data (at E = 800 MPa 20% calcification; at E = 2000 MPa 40% calcification; www.intechopen.com at E = 8000 MPa 70% calcification and at E = 18000 MPa 100% calcification occurs ) taken from Davy and Connoly (1982) for the new bone at intermediate densities corresponding to the elastic moduli of the callus adopted from the study of Gardner et al. (2000 ). The calculated calcification of the different regions of callus at 4, 8, 12 and 16 weeks is shown in Table 3. The porosity of the callus tissue was then assumed to be inversely related to the proportion of the calcified tissue present in the callus (Carter 1977), and therefore was found to vary from 0.9 for soft callus (<5% calcification) to 0.3 for woven bone (100% calcification).
As the porosity of 0.8 was used by Van Driel et al. (1998) for soft fibrous, cartilage and bone tissues, by comparison 0.9 appeared valid for the softer tissues of the present model. The porosity of 0.3 was suggested for woven bone (Carter 1977) which seems to be valid for the harder tissues of the present model. All intermediate values of porosities were linearly extrapolated from these two values at corresponding values of calcification shown in Table  3. The interstitial fluid in the callus was assumed to have the bulk modulus (K f ) of salt water (2.3 GPa) (Cowin 1999 Table 4. Longitudinal stress ranges in the monophasic model of Gardner et al. (2000) and effective stress ranges in the biphasic poroelastic model.

Discussion and conclusions
The magnitudes of fluid pressure and effective stress are very high and therefore appear unrealistic. In particular, the high tensile stress would produce cavitation or may lead to gas in the pore fluid, and the callus matrix may be ruptured. These conditions are typical of the undrained simulation behaviour under large deformations. In reality, no matter how small the permeability of the poroelastic medium or how rapid is the loading, fluid flow will occur under such high-pressure gradient. Therefore the absolute value of pressure and the presence of high tensile pressure in the callus are more an artifact of the modelling technique than are the patterns of the distribution of pressure and their trend in variation. Thus, for this technique of modelling the presence of spatial or temporal pressure gradients within the callus may be a valid indicator of the flow of fluid in the form of blood, nutrients or waste products. For example, the tensile fluid pressure regions of the present models at peak loading during walking are likely to show reduced magnitudes of tensile pressure or compressive pressure during unload phases. The hydraulic gradient will be reversed and the fluid flow will be in the opposite direction. This alternate inflow and outflow of fluid may be related to the transport mechanism of inflow of nutrients and oxygen, and outflow of waste products and carbon dioxide from the callus. Such inflow and outflow may enhance the growth of capillary blood vessels, thus accelerating the healing process. It can be envisaged from the above that if the movement is too small, the change in fluid pressure will be small and the beneficial effect will also be small. Therefore the results of the present study corroborate those of other studies (Kenwright 1998, Sarmiento and Latta 1995, Goodship and Kenwright 1985, Kenwright and Goodship 1989 suggesting that fracture site movement is necessary for efficient secondary healing. However, if the movement is too large the fabric of the callus matrix could be damaged because of the cyclical expansion and compression, and this damage could hinder healing. Furthermore, if the frequency of the movement were too high, the time interval between pressure gradient reversals would be too small for the fluid to flow. On the other hand, if the frequency is too low, fluid will penetrate the tissues before any substantial hydraulic gradient can be developed. Since the present study was limited to a single temporal point during the gait cycle it is not possible to define the optimum magnitude and frequency of movement beneficial for healing. Models that simulate the different temporal points of the gait cycle may provide more information about the optimum movements for fracture patients. At Week 4 (Figure 3a), there are high compressive fluid pressures in the inter fragmentary gap regions because the undrained model is under large compressive displacement and fluid is unable to flow out side the system boundary. In this condition the loads is predominantly taken by an incompressible fluid that controls the motion of the bone fragments and provides support for intact tissues. It does this by increasing the stiffness of the limb and it also protects the fracture from further damage (Sarmiento and Latta 1995). Since fluid pressure is a function of movement, the cortical gap locations are expected to undergo higher pressures than locations further away from the gap, as shown in Figure 3a.
It is worth noting that at this stage the callus is comprised of more than 90% fluid, therefore the incompressible fluid will resist high pore pressure. At Week 8, reduction in the porosity (Table 3) and compressive interfragmentary displacement (Table 1) (Table 1), resulting in elevated fluid pressures in the adjacent callus. At Week 16, both the applied compressive interfragmentary displacement and the fluid pressure in the callus reduce. Therefore fluid pressures appear to be more sensitive to longitudinal interfragmentary displacements than to the callus-porosity in the present study. Fluid pressure distribution patterns correlate with the general pattern of ossification, as reported by others (Blenman 1989, Sarmiento 1995, Yamagishi 1955. Blenman and Carter (1989) suggested that ossification progresses through the stages of `bone tuft', `bone wedge' and `bone bridge' and poroelastic models appear to corroborate this if it is believed that ossification may not take place in the regions of high fluid pressures. At Week 4, high pore pressure regions of the interfragmentary gap divide the proximal and distal callus. Therefore it appears that ossification is possible only in the low-pressure regions of the periosteal callus away from the interfragmentary gap, forming the `wedge' shaped ossified callus. At Weeks 8 and 12, fluid pressure in the callus is reduced almost to zero at the level of the gap, thus allowing the formation of a `bridge' of ossified callus between the `wedges'. Since the present study started at 4 weeks post operation, a pattern similar to that of the `tufts' theory may have also been present before this stage. The similarity between the magnitude and pattern of total stress diagrams of the biphasic poroelastic model and the corresponding longitudinal stress diagrams of the monophasic model of Gardner et al. (2000) are expected. This is because total stress is a function of the total force and total cross section area (Wood 1990) of the callus, which remain similar in both the models. The greatest disparities between the monophasic and biphasic solutions occur in initial healing at Week 4. This disparity exists because the soft callus has a high porosity initially and is subjected to high tensile stresses during the large applied initial displacements. As the callus calcifies, its porosity and the fluid pressures decrease so that total stress is closer to the effective stress. This effect of porosity is evident from Table 4, where the maximum difference between stresses from the two models is in the high-porosity central callus and the minimum difference is in the low-porosity peripheral callus. The patterns of effective stress in the biphasic models also differ from the corresponding monophasic models. In the biphasic models, substantial variations of stress are evident at the cortical gap, sub periosteal and endosteal callus. Whereas in general, the monophasic models reported by Gardner et al. (2000) predicts similar stress regimes throughout the central and adjacent callus. Therefore, if the differentiation and maturation of the callus are believed to be influenced by the preceding stress environments, then the biphasic models appear to predict more realistic patterns of tissue differentiation and maturation.

Limitations of the study
The results of the present study should be evaluated under modelling limitations. Firstly, solutions are valid only for the `undrained' condition that assumes that no fluid moves out of the system boundary during loading, but in reality a small amount of fluid may drain through the pores of the callus. However this drainage may not have invalidated the results of the present study as the gait cycle frequency is approximately 1 Hz, and physiological loading periods during the stance phase of gait are around 0.3 to 0.5 seconds, which is probably too rapid for significant drainage of fluid to take place (Gardner et al. , 2000. Secondly, the callus has been idealised as a linear, elastic, fully-saturated porous medium throughout healing. These idealised conditions may also be responsible for the high magnitudes of fluid pressure, whereas actual pressures are believed to be lower than predicted by the present model. However, the pattern and trend of temporal variations in fluid pressures during progressive ossification are unlikely to change significantly as a result of applying slightly different poroelastic material properties and constitutive equations. Despite the limitations of the present model, these results indicate that the biphasic material properties of the callus are more appropriate to the initial soft callus stage of healing and support the suggestion of Sarmiento and Latta (1995) "The incompressible fluid effect, or hydraulics is most important in early post injury period. We feel that hydraulics is responsible for the control of motion of fragments before callus has developed and that it provides the significant degree of stiffness observed in loaded limbs with fresh fractures fit with fracture braces."

Acknowledgement
The author wishes to acknowledge (a) the financial support obtained by commonwealth Commission UK during the study (b) the major input and supervision from Prof AHC Chan and Dr Trevor Gardner, University of Birmingham UK.