Numerical description and experimental validation of a rheology model for non-Newtonian fluid flow in cancellous bone

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

Available online at www.sciencedirect.com

www.elsevier.com/locate/jmbbm

Research Paper

Numerical description and experimental validation of a rheology model for non-Newtonian fluid flow in cancellous bone Rene´ P. Widmer Soykaa,n, Alejandro Lo´pezb, Cecilia Perssonb, Luca Cristofolinic,d, Stephen J. Fergusona a

Institute for Biomechanics, ETH Zurich, Schafmattstrasse 30, HPP O 14, 8093 Zurich, Switzerland Division of Applied Materials Science, Department of Engineering Sciences, The Ångström Laboratory, Uppsala University, Uppsala, Sweden c Laboratory for Medical Technology, Rizzoli Orthopaedic Institute, Bologna, Italy d Department of Industrial Engineering, University of Bologna, Bologna, Italy b

art i cle i nfo

ab st rac t

Article history:

Fluids present or used in biology, medicine and (biomedical) engineering are often

Received 3 April 2013

significantly non-Newtonian. Furthermore, they are chemically complex and can interact

Received in revised form

with the porous matrix through which they flow. The porous structures themselves display

6 June 2013

complex morphological inhomogeneities on a wide range of length scales. In vertebro-

Accepted 13 June 2013

plasty, a shear-thinning fluid, e.g. poly(methyl methacrylate) (PMMA), is injected into the

Available online 28 June 2013

cavities of vertebral trabecular bone for the stabilization of fractures and metastatic lesions. The main objective of this study was therefore to provide a protocol for numerically investigating the rheological properties of PMMA-based bone cements to predict its spreading behavior while flowing through vertebral trabecular bone. A numerical upscaling scheme based on a dimensionless formulation of the Navier–Stokes equation is proposed in order to relate the pore-scale rheological properties of the PMMA that were experimentally estimated using a plate rheometer, to the continuum-scale. On the pore length scale, a viscosity change on the order of one magnitude was observed whilst the shear-thinning properties caused a viscosity change on the order of only 10% on the continuum length scale and in a flow regime that is relevant for vertebroplasty. An experimental validation, performed on human cadaveric vertebrae (n¼9), showed a significant improvement of the cement spreading prediction accuracy with a non-Newtonian formulation. A root mean square cement surface prediction error of 1.53 mm (assuming a Newtonian fluid) and 1.37 mm (assuming a shearthinning fluid) was found. Our findings highlight the importance of incorporating the nonNewtonian fluids properties in computational models of porous media at the appropriate length scale. & 2013 Elsevier Ltd. All rights reserved.

n

Corresponding author. Tel.: +41 44 633 4060. E-mail addresses: [email protected], [email protected] (R.P. Widmer Soyka).

1751-6161/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jmbbm.2013.06.007

44

1.

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

Introduction

Many of the fluids present or used in biology, medicine and (biomedical) engineering are often significantly nonNewtonian. They demonstrate in differing degrees a varying time- and shear-dependence of the viscosity. Frequently, they are multi-component and chemically complex, and may react with the porous domain through which they are flowing and with the other fluids with which they come into contact within the porous structure (Pearson and Tardy, 2002). These structures are themselves geometrically highly complex and display architectural inhomogeneities over a wide range of length scales. For example, poly(methyl methacrylate) (PMMA) bone cements are used in orthopaedics for vertebral augmentation procedures such as vertebroplasty, one of the most important techniques for the stabilization of osteoporotic compression fractures and the fixation of other weakening lesions such as tumours (Heini et al., 2000). During the procedure, the cement – typically a mixture of powders (poly(methyl methacrylate) and an X-ray contrast agent, e.g. barium sulphate) and liquid (methyl methacrylate monomer) – is injected through a cannula. Cement leakage is a common occurrence and represents the main complication risk for vertebroplasty (Laredo and Hamze, 2004). Many factors have an influence on the injection biomechanics and, ultimately, on whether or not cement leakage occurs (Baroud et al., 2004). Cement viscosity has been identified as an important parameter for controlling the spreading pattern and therefore the risk of cement leakage (Baroud et al., 2004, 2006; Cotten et al., 1996; Bohner et al., 2003). Furthermore, the spreading pattern is believed to influence the mechanical outcome (Liebschner et al., 2001; Tschirhart et al., 2005). Even though the importance of the bone cement rheology has been acknowledged as a key parameter of the injection biomechanics of vertebroplasty, only few studies have addressed it. The uniformity of the cement filling of the vertebral body correlates highly to the cement viscosity (Baroud et al., 2006; Loeffel et al., 2008), but the local tissue morphology, that varies among the different patients being treated, represents a high variability in cement infiltration behavior in general. Additionally, pore-scale and continuumscale viscosity are unique, due to the nonlinear dependence of the viscosity on local tissue morphology and fluid deformation rate. Consequently, as a general conclusion, determining the infiltration behavior for bone augmentation a priori based on clinical or cadaveric studies is difficult (Bohner et al., 2003). Measuring the local tissue morphology absolutely is de facto an impossible task to handle. X-ray dosage limitations and geometrical constraints do not allow to scan the patients using micro-computed tomography ðμCTÞ prior to the treatment. This is a certain limitation that inhibits the utilization of computational models describing the pore-scale rheological behavior of the injected cements (Beaudoin et al., 1991; Baroud and Yahia, 2004; Chui et al., 2006; Lian et al., 2008). The rheology models presented in Baroud and Yahia (2004) are moreover limited to the characterization of the bone cement flow through the injection cannula. The geometry of the cannula is considered to be the most important factor influencing the injection force needed

to drive the bone cement into the vertebral body (Baroud et al., 2004; Baroud and Steffen, 2005; Gisep and Boger, 2009). On the other hand, the design of the cannula shaft does not influence the filling pattern. Therefore, these models are only of minor importance when it comes to the description of the cement infiltration process in vertebroplasty. Bohner et al. (2003) presented an analytical and experimental model to characterize the injection of bone cements into a porous structure with a special focus on modeling the cement extravasation. This model overcomes the aforementioned limitations, but assumes a spherical spreading pattern and Newtonian fluid properties, i.e. an independence of the fluid viscosity on the fluid deformation rate. This implies that the risk of extravasation does not depend on the injection speed. Studies have shown that PMMA-based bone cements are, in addition to the time-dependency of the viscosity, shearthinning (Farrar and Rose, 2001), thus cement viscosity decreases with an increase of the fluid deformation rate, and the risk of cement leakage is increased. Clinically highly relevant models of the infiltration process were presented by Teo et al. (2007a) and the authors themselves (Widmer and Ferguson, 2011). Specimen-specific and continuum-scale computational fluid dynamics (CFD) models were used to predict the bone infiltration process and filling pattern combined with methods to assess the mechanical outcome of the vertebroplasty. The authors’ model focused on the simulation of the displacement of the bone marrow by the injected biomaterial and employed a simplified fluid rheology model. In the model of Teo et al. (2007a), the viscosity of the infiltrating fluid was related to its shearing rate according to the power law equation μ ¼ C  γ_ r−1 ;

ð1Þ

where μ refers to the pore-scale or true viscosity, C to the consistency index, γ_ to the fluid deformation or shearing rate and r to the power law index. The shearing rate was calculated on the basis of the strain rate tensor γ_ , which was obtained from the kinematic equations from the CFD code according to rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 γ_ : γ_ T ¼ γ_ : γ_ ; ð2Þ γ_ ¼ 2 2 since γ_ is a symmetric tensor. The strain rate tensor is given by γ_ ≈∇q þ ∇q T

ð3Þ

where q is the apparent or Darcy velocity vector that represents the average fluid velocity in the porous medium. Hence, Eq. (3) is truly only an approximation. The effective pore-scale velocity and deformation rates might differ by orders of magnitude for fluid flowing along the tortuous paths of the porous medium and fluid flow observed at the continuum length scale. The challenge is that Eq. (1) is defined at the pore length scale while Eq. (3) estimates the deformation rates at the continuum length scale and consequently is not able to express the fine details of the effective and highly tortuous fluid flow paths. This misunderstanding of the varying effects, occurring at different length scales, is potentially the main reason why in the past researchers failed to successfully describe and characterize the rheology of non-Newtonian

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

bone cements that are injected into the cancellous bone of vertebral bodies during the augmentation process. The main motivation for this study is to improve the understanding of the viscosity as the key parameter of the vertebroplasty procedure and to overcome the technical challenges and limitations of the current state of the art. The objective is therefore to provide a protocol to numerically investigate the rheological properties of bone cements and predict its spreading behavior while flowing through vertebral trabecular bone. The protocol consists of the following steps: (a) simulate the flow of PMMA-based bone cement through representative bone structures at the pore-scale length; (b) estimate a scale-length invariant number in order to perform an upscaling of the pore-scale solution; (c) repeat the steps (a) and (b) while systematically varying the shearthinning effect; (d) curve fitting for deriving an approximate continuum-scale model that describes the cement viscosity at the continuum length scale. The continuum-scale model is validated in an experimental setup, i.e. by injecting PMMAbased bone cements into cadaveric human vertebral bodies. We hypothesize that the proposed upscaling scheme has a positive effect on the prediction accuracy of the computational flow model that has been proposed previously by the authors (Widmer and Ferguson, 2011).

approach is based on the assumption – confirmed by all detailed pore-scale models – that the definitions and intrinsic properties of the fluid velocity, fluid pressure and the permeability of the porous matrix will not be altered by the rheology of the fluid. Consequently, all non-Newtonian effects will be considered in a suitable definition of μ at the continuum length scale which still has the dimensions of Newtonian viscosity (Pearson and Tardy, 2002). This apparent or Darcy viscosity μ is therefore determined from the pressure drop due to steady flow through a sample of the porous matrix and is a function of the rheological properties of the fluid (including the dependency on the shearing rate and the temperature), the apparent fluid velocity q and the morphology of the porous matrix. In a clinical setting, only a weak correlation between viscosity during hardening and the ambient temperature was found (Boger et al., 2009). Hence, temperature effects are neglected. Summarizing, from Eq. (4) we obtain for the apparent viscosity μ ¼−

q ∥q ∥2

Methods

2.1.

Theoretical

Most often, mathematical and numerical models that quantify the transport process within a porous medium are based on continuum descriptions that hold on average. A first attempt to characterize the flow of Newtonian fluids, i.e. fluids where the shear stress exerted by the fluid is proportional to its deformation rate and the fluid viscosity is the constant of proportionality, is Darcy's law (Darcy, 1856) that relates the gradient of the pore pressure to the apparent fluid velocity within a rigid porous medium q ¼−

ks μ

∇p;

ð4Þ

where q refers to the apparent fluid velocity vector, k s to the intrinsic permeability tensor, μ to the fluid viscosity and ∇p to the pressure gradient. However, the validity of this model is given only if (a) the fluid is Newtonian and (b) the flow is dominated by viscous forces, i.e. the ratio of inertial to viscous forces (continuum-scale Reynolds number Re) is considerably smaller than one. Most often, q , p and k s are functions of the position in space x and time t (Pinder and Gray, 2008). An approximation that has been adopted for simplicity throughout this study, is that both the solid and fluid phases are incompressible (Pearson and Tardy, 2002). Consequently, the principle of mass conservation reduces to the principle of volume conservation ∇  q ¼ 0:

ð5Þ

A major limitation of Darcy's law is that it is valid only for Newtonian fluids. Hence, it is necessary to extend Eq. (4) to cover fluids with non-Newtonian rheological behavior. Our

k s ∇p:

ð6Þ

For trabecular bone, the continuum length scale Reynolds number Re can be estimated (Cowin and Telega, 2003) according to Re ¼ 4

2.

45

ρ∥q ∥ μSv

;

ð7Þ

where ρ is the fluid density and Sv is the specific surface of the bone volume. All variables of Eq. (7) can be estimated at the continuum length scale (Cowin and Telega, 2003). However, Re is traditionally estimated at the pore length scale. This correspondence provides the basis of the proposed upscaling scheme to relate pore-scale viscosity μ to Darcy viscosity μ. Generally, the dimensionless and steady-state Navier–Stokes equations are given by Reðu ′  ∇Þu ′ ¼ −Re∇p′ þ ∇2 u ′ ;

ð8aÞ

∇  u ′ ¼ 0;

ð8bÞ

u ′ ¼ u =u0 ;

ð8cÞ

p′ ¼ p=ðρu20 Þ:

ð8dÞ

u ′ is the dimensionless velocity, p′ the normalized pressure and u0 a characteristic velocity. For non-Newtonian fluids, Re is defined (Murdock, 1976) as 2−r

Re ¼

ρu0 Dr ; C

ð9Þ

where D is a characteristic length that can be determined from the porosity β of the bone and the specific surface according to Cowin and Telega (2003) D¼6

1−β Sv

Sv ¼ 323  β−939  β2 þ 1340  β3−1010  β4 þ 288  β5

ð10aÞ ð10bÞ

The upscaling is achieved through homogenization, i.e. averaging of the pore-scale Reynolds number over a reference volume and relating it to a Darcy flow regime of equal

46

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

characteristics Re≡

∮u′ ^ Reðx; y; zÞ ds : ∮u′ ^ ds

ð11Þ

^ of the The line integral is evaluated along the ridges u′ three-dimensional fluid velocity field and u0 is defined as the fluid velocity at the spatial location ðx; y; zÞ. Eq. (11) must not be treated in terms of an exact mathematical statement, thus should be understood in a more relaxed sense.

2.2.

Experimental

2.2.1.

Pore-scale rheology

Rheological data was acquired using an AR 2000 rheometer (TA Instruments, New Castle, DE, United States) with a 19 mm diameter parallel-plate titanium geometry. Commercial Simplex P was prepared with a modified formulation of 10 g PMMA powder and 10 g monomer liquid, to resemble an injectable cement for vertebroplasty, since Simplex P is intended for prosthesis fixation. 1 g BaSO4 was added on top of the powder (gross liquid-to-powder ratio of 0.909 mL/g) to increase the radio-opacity. The powder and the liquid phases were mixed for 30 s in 50 mL falcon tubes, using a capvibrator (Ivoclar Vivadent AB, Solna, Sweden). Specimens consisting of approximately 2 mL of the mixture were scooped onto the lower plate to avoid pre-shearing, and the gap between the plates was set at 1000 μm. The equilibration time was 30 s after which the sampling started and continued every 5 s for up to 15 min from the start of mixing. The viscosity over time for the modified Simplex P was obtained at six steady shear rates ðf0:1; 0:5; 1:0; 10:0; 25:0; 100g s−1 Þ and a 25 1C lower plate temperature using a flow test in peak hold mode. Three specimens, each specimen coming from a separate mix, were tested per shear rate. The chosen shear rates and temperatures were intended to encompass a range of flow and temperature conditions to which the bone cement might be subjected during injection. A power law in the form of μð_γ ; tÞ ¼ CðtÞ_γ rðtÞ−1 ;

ð12aÞ

þ φrheo  t; CðtÞ ¼ φrheo 0 1

ð12bÞ

þ φrheo t rðtÞ ¼ φrheo 2 3

ð12cÞ

was fitted to the experimental data, where t is the relative …φrheo are the fitting time to the start of the mixing and φrheo 0 3 coefficients.

2.2.2.

Continuum-scale injection trials

Nine cadaveric vertebral bodies (four thoracic vertebrae (T11) and five lumbar vertebrae (4  L1 and 1  L5)) from four male donors (aged 70–88) were provided by the International Institute for the Advancement of Medicine (IIAM, 175 May Street, Edison, NJ 08837). Computed tomography (CT) scans (BrightSpeed, GE Healthcare, 3000 N. Grandview Blvd. Waukesha, WI 53188) of the specimens were taken prior to and after the cement injections at the Istituto Ortopedico Rizzoli. The CT scans were performed in helical mode with 0.195 mm pixel spacing in the transverse plane and 0.625 mm slice thickness. The same Simplex P bone cement formulation as used for the rheology measurements was utilized. The solid and liquid phases were mixed for 60 s. 6 mL of the mixture were filled

into the syringe barrel. The cement was delivered with the support of a computer-assisted injection device (Loeffel, 2007) through 100 mm long biopsy needles with an inner diameter of 3 mm. For four specimens, a bipedicular access was chosen (alternating pedicles approx. after every 1 mL cement delivery). For the remaining five specimens, the cement was delivered through one pedicle (2  left, 3  right). The purpose of alternating the access strategy was to accommodate the requirements of a parallel study and the effect was not examined in this study. The injection was manually governed through the force feedback-controller of the injection device. The injection pressure and volume were continuously recorded. The injection was stopped either after the complete delivery of the mixture portion, or after the occurrence of any cement leakage found through visual inspection of the cortical shell.

2.3.

Computational

2.3.1.

Pore-scale flow simulations

Two cadaveric human vertebral trabecular bone specimens (T12 and L1 level, two individual donors) were harvested and scanned using a high resolution computed tomography scanner (μCT 40, Scanco Medical AG, Bassersdorf, Switzerland) at the Department of Biomedical Engineering of the Eindhoven University of Technology. The scans, performed in air, were taken at a nominal isotropic resolution of 41 μm and segmented using a global threshold. Each specimen scan was divided into totally 28 sub-volumes of 6  6  6 mm size. The 3D SLICER software package was used to convert the sub-volumes into a smooth surface mesh. The smoothing procedure consisted of a pre-filtering step (Gaussian filter kernel, s ¼ 0:8). For each sub-volume, the iso-value of the marching-cubes algorithm was adjusted to restrict the difference between the volume enclosed by the surface mesh and the original bone volume to 2%. The surface meshes were then converted into volume meshes (μCT finite element models) using the tetrahedral mesh generator NETGEN (NETGEN V4.4, http://www.hpfem.jku.at/netgen). Further represen tative bone specimens were generated using a cellular model of trabecular bone (Kim and Al-Hassani, 2002). The horizontal and vertical trabecular spacing was varied from 0:80 mm ≤ Tb:Spðh;vÞ ≤1:60 mm whereas the trabecular thickness was varied from 0:10 mm ≤Tb:Thðh;vÞ ≤0:28 mm (three equally spaced values per variable) to encompass a physiological set of trabecular bone morphology (Hulme et al., 2007). Representa tive examples of μCT and cellular models are depicted in Fig. 1. The residual form of Eq. (8) was implemented into the open source Finite Element (FE) library libMesh (Kirk et al., 2006). The fixed-point algorithm (Kerkhoven and Saad, 1991) was applied to solve the non-linear problem that arises from Eqs. (1), (8) and (9). The fluid flow was simulated for each principal axis by applying a volumetric flow rate boundary condition of magnitude Q at the inlet and a reference pressure p ¼ p0 at the outlet, respectively, thus enforcing a fluid flow across the principal direction under investigation. The flow rate Q at the inlet was iteratively adjusted by using the bisection algorithm for matching the predetermined continuum-scale Reynolds number Re. On the bone-void interface, a no-slip boundary condition ðu ¼ 0Þ was enforced.

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

47

Fig. 1 – The volume rendering (grey) of the trabecular bone and the corresponding tetrahedral mesh (green) of a representative μCT finite element model is shown in (a). (b) shows a close-up view on the unit cell model. The blue box indicates the boundaries of the unit cell for the subsequent analyzes. The unit vectors e i , defining the principal axis directions, are also indicated. The transversal direction is a synonym for the x; x- and z; z-directions and the longitudinal direction for the y; y-direction. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

On the model walls, a zero-shear slip boundary condition was applied (∂u=∂n ¼ 0, where n is the coordinate normal to the wall) to mimic the flow through a semi-infinite plane. Postprocessing of the CFD results was carried out using ParaView 3.10.0 (Kitware, Inc., New York, USA). To establish the correspondence between the continuum-scale Reynolds number and the (normalized) Darcy viscosity, a regression model in the form of up up up μ up up ¼ φ0 ðφ1  ReÞφ2 Reþφ3 Cþφ4 C

ð13Þ up

up

was fitted to the predicted viscosity data. φ0 …φ4 fitting coefficients of the regression model.

2.3.2.

are the

Continuum-scale flow simulations

The surface models of the vertebral bodies were converted into tetrahedral volume meshes using NETGEN, and each model consisted of approximate 50,000 elements. For each element, the bone volume ratio BV=TV was determined from the Hounsfield numbers of the CT stacks (Teo et al., 2007b) and used to determine the anisotropic permeability tensor (Widmer and Ferguson, 2012). A reference pressure of 0 Pa was applied to the whole cortical shell, but not to the endplates, of the vertebrae. The fluid flow through the endplates was inhibited, as the specimens were embedded into PMMA caps at the superior and inferior aspect. Flow rate boundary conditions were applied at the cannula tip to match the volume profile that had been recorded during the experiment. The flow was governed by Eq. (4) and implemented into the libMesh FE library. A mixed boundary formulation was used to track the cement/bone marrow fluid interface (Widmer and Ferguson, 2011). The simulations were performed with and without consideration of the shear-thinning effect.

2.4.

Error quantification and statistical evaluation

The vertebral bodies and, in case of the post-augmentation scans, the injection cannulae and cement patterns were manually segmented using the 3D SLICER software package (3D SLICER V.3.6.21.0, http://www.slicer.org). A rigid surface

registration was performed to realign the post-augmentation models to the pre-augmentation models. After the re-align ment, the error distribution of the predicted and experimen tal cement distribution was quantified using the Hausdorff distance (Huttenlocher et al., 1993). First, the errors were assigned to groups according to the study factors pseudoplasticity (considered; PP+/neglected; PP−) and spine region (thoracic/lumbar), for detecting differences between the prin cipal study groups. Second, the errors were reallocated to the groups thoracic/lumbar (pseudo-plasticity ignored; Thoracic PP7) or pseudo-plasticity considered/neglected while ignoring the spine level, for detecting differences between the anatomi cal site or the Darcy viscosity model. The Kruskal–Wallis test with post-hoc multiple comparison (Tukey–Kramer signifi cance criterion; per family error rate of PFE ¼ 0:05) was used to detect the differences between the groups. The Spearman correlation coefficients were computed for the normalized viscosity, the rheological parameters C and r and the morpho logical parameters bone volume fraction BV=TV, degree of anisotropy DA, the trabecular spacing Tb:Spðh;vÞ and the trabecular thickness Tb:Thðh;vÞ . p-values for Spearman's ρ were determined using the exact permutation distributions, a value of po0:05 was considered as significant. The fitting error indicators mean squared error (MSE), root mean squared error (RMSE) and mean percentage error (MPE) were computed to quantify the goodness of the upscaling model fit.

3.

Results

The measured native viscosity and the fitted power law are shown in Fig. 2 and illustrate the shear-thinning behavior of the PMMA-based bone cement. Due to self-expulsion of the specimens, the measurements became unstable after 280 s at the highest deformation rate (100 1/s) and the data were excluded for the nonlinear fitting process on the pore-scale viscosity regression model (Eq. (12)). The fitting process indicated …φrheo ¼ f−115:07; 1:1798; −1:8051; 0:0005g. the coefficients φrheo 0 3 All coefficients differ statistically significant from zero

48

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

Fig. 2 – (a) Measured and fitted viscosity vs. the time measured at 25 1C lower plate temperature. For better readability, only every 3rd measured data point is drawn. Red markers refer to measurements that were excluded from the fitting process. (b) The predicted viscosity μðγ_ ; tÞ vs. the measured viscosity μ. R2 ¼ 0:987, MPE¼25.5%. The viscosity is over-predicted at the lower end whereas it is under-predicted at the upper end of the scale on average. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

Table 1 – Spearman correlation coefficients of different parameters and the viscosity ratio μ=C. The individual lines correspond to the different bone models ðμCT=cellular modelÞ and flow directions (transverse/longitudinal directions). The isotropic models neglect the influence of the flow direction. All values differ statistically significant from zero. Model

Flow dir.

r

Re

BV=TV

DA

Tb:SpðhÞ

Tb:SpðvÞ

Tb:ThðhÞ

Tb:ThðvÞ

Unit cell

Transv. Longitud. Isotropic Transv. Longitud. Isotropic

−0.30 −0.35 −0.31 −0.34 −0.48 −0.38

−0.88 −0.85 −0.86 −0.86 −0.79 −0.82

−0.09 −0.15 −0.10 −0.04 −0.02 −0.03

0.02 0.08 0.04 −0.03 0.03b −0.01c

0.05 0.03 0.04 0.09 0.05 0.07

0.08 0.14 0.10 0.08 0.05 0.08

−0.03 −0.05 −0.04 −0.04a −0.02a −0.03a

−0.03 −0.05 −0.03

μCT

a b c

A classification according to the structural orientation was not made available by the image processing software. po0:0001; po0:00001 otherwise. po0:01.

ðpo0:00001Þ. The regression model explained 98.7% of the measured variance at the MPE level of 25.5%. The Spearman correlation coefficients of different morphological and rheological parameters and the normalized Darcy viscosity μ=C are presented in Table 1. The individual lines correspond to the different bone models (μCT/cellular model) and flow directions (transverse/longitudinal directions). μCT models were derived from μCT scans of representative bone samples whereas the cellular models were derived from the analytical description of a hexagonal-like structure. The isotropic models neglect the dependence on the flow direction. The μCT scans were analyzed by using the BoneJ plugin (BoneJ V1.3.8; Doube et al., 2010). This software did not provide algorithms to calculate Tb:Thðh;vÞ classified according to the structural orientation. Therefore, only Tb:Th was calculated for the μCT models. All numbers are statistically significant (po0:01). A strong correlation was found for the continuum Reynolds number Re. The power law index r was moderately correlated to the Darcy viscosity μ. Only a

weak correlation was found for the morphological parameters BV=TV, DA, Tb:Spðh;vÞ and Tb:Thðh;vÞ . Fig. 3 shows the normalized viscosity μ=C as a function of the continuum-scale Reynolds number Re for r ¼ f0:0; 0:4; 0:8g and the fitted upscaling model. The error bars indicate min/ max values and the markers correspond to the median value predicted for the different cellular and μCT models of a particular fluid type and flow regime. To improve the readability of the data, the data of the transverse and longitudinal directions were aggregated into one graph. For a fluid with the strongest pseudo-plastic behavior ðr ¼ 0:0Þ, the shearthinning effect causes an alteration of −6.38% (Re ¼ 100 , cellular model) to 17.2% (Re ¼ 10−8 , μCT model). For 10−8 ≤ Re ≤10−1 , a log-linear trend is observable. Inertial forces start to dominate the flow for Re410−1 , which is manifested in the deviation of the predicted values from the almost perfect loglinear trend. The fitting coefficients of the Darcy viscosity model (Eq. (13)) for a pseudo-plastic fluid are shown in Table 2. All coefficients differ statistically significant from

49

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

Fig. 3 – The normalized viscosity fraction μ=C vs. the Darcy Reynolds number Re. The results referring to the individual principal directions were pooled into the same graph. With increasing pseudo-plasticity effect ðr-0Þ, the curve deviates from a horizontal line, i.e. the response of a pure Newtonian fluid. With the Reynolds number increasing, inertial forces start to play a role. This is expressed by the increasing steepness of the curve for large Reynolds numbers. (a)–(c) indicate median 7 min/max values. (d) shows the relationship if all data were aggregated into one dataset. up

up

Table 2 – The fitting coefficients of the Darcy viscosity models for 0 ≤r ≤1. φ0 …φ4 are the fitting coefficients. All coefficients up up differ statistically significant from zero (φ2…4 ) or one ðφ0…1 Þ, respectively ðpo0:00001Þ. The individual lines correspond to the different bone models ðμCT=cellular modelÞ and flow directions (transverse/longitudinal directions). The isotropic models neglect the influence of the flow direction. The pooled model accumulates all collected data into one regression fit. Notice that the fitting coefficients depend on the morphological properties of the porous matrix but are independent of the specific rheological properties of the shear-thinning fluid. up

up

up

up

up

Model

Flow dir.

φ0

φ1

φ2

φ3

φ4

Unit cell

Transverse Longitudinal Isotropic

1.0068 1.0082 1.0072

121.83 51.750 90.909

0.0111 0.0113 0.0112

−0.0109 −0.0110 −0.0109

−0.0025 −0.0038 −0.0029

μCT

Transverse Longitudinal Isotropic

1.0061 1.0083 1.0066

53.776 4.6998 22.972

0.0103 0.0099 0.0101

−0.0101 −0.0096 −0.0099

−0.0024 −0.0061 −0.0032

Pooled

Isotropic

1.0069

47.085

0.0106

−0.0104

−0.0030

up

up

zero ðφ2…4 Þ or one ðφ0…1 Þ, respectively ðpo0:00001Þ. The coefficients depend on the morphological properties of the porous solid but are valid for all pseudo-plastic fluids within

the boundaries of the investigated flow regime. The goodness of the fit is quantified by the error indicators shown in Table 3. The fitting coefficients and goodness of fit were

50

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

Table 3 – Different error indicators (MSE, RMSE, MPE) that summarize the quality of fit of the power law that relates the continuum-scale Reynolds number Re to the predicted Darcy viscosity μ. The individual lines correspond to the different bone models ðμCT=cellular modelÞ and flow directions (transverse/longitudinal directions). The isotropic models neglect the influence of the flow direction. The pooled model accumulates all collected data into one regression fit and deviates the most from the predicted data on the average. Model

Flow dir.

MSE ½10−5 s2ð1−r Þ 

RMSE ½10−3 s1−r 

Unit cell

Transverse Longitudinal Isotropic

4.52 10.9 8.21

6.72 10.5 9.06

5.25 11.6 8.22

μCT

Transverse Longitudinal Isotropic

4.10 4.80 10.4

6.38 6.93 10.2

4.49 5.54 10.2

Pooled

Isotropic

9.57

9.78

MPE ½10−3 %

215

Fig. 4 – (a) demonstrates typical pressure and injection volume curves of a PMMA-based bone cement injection using the motor-driven delivery device. The recorded cement volume profile was subsequently used to drive the continuum-scale model. During the progress of the injection, the cement hardens and a higher injection pressure is needed to force the cement into the vertebral body. (b) Segmented experimental cement cloud colored by the prediction error of the continuum-scale model. The principal source of error is uncertainty in the location of the cannula during the in-vitro experiments.

Fig. 5 – Box plots of the analyzed error groups. (a) shows the error characteristics of the injection trials grouped according to the level and the underlying continuum-scale viscosity model (principal study groups). (b) shows the group characteristics if the errors were grouped purely based on the spine level or the continuum-scale viscosity model (intermixed study groups). PP+ indicates the consideration of the shear-thinning/pseudo-plasticity effect whereas PP− indicates the neglect of the pseudo-plasticity. PP7 indicates that the group is intermixed in relation to pseudo-plasticity. The statistical tests indicate for both the principal and intermixed study groups differences between the subgroups ðpo0:00001Þ. nPair-wise differences between subgroups suggested by the post-hoc analysis.

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

calculated for the individual models and a pooled set of data. A maximum MPE of 0.215% was estimated for the pooled model. A representative profile of injection pressure and volume that has been recorded during one of the vertebral injection experiments, is depicted in Fig. 4(a). During the progress of the injection, the pressure required to force the cement into the vertebral body increases. This characteristic is typical for the injection mechanics of vertebroplasty. Fig. 4(b) furthermore shows the cement distribution that has been reconstructed from the post-augmentation CT scan, for a bipedicular injection strategy. The superimposed colors refer to the cement surface prediction error of the continuum-scale model. The maximum error was estimated as 4.60 mm and 5.20 mm (with and without consideration of the shearthinning effect, respectively). The group RMSE were estimated as 1.93 (Thoracic PP−), 1.33 (Lumbar PP−), 1.71 (Thoracic PP+), 1.29 (Lumbar PP+), 1.83 (Thoracic PP7), 1.31 (Lumbar PP7), 1.37 (PP+) and 1.53 (PP−). All numbers are statistically significant. Fig. 5 shows the box-plots of the cement/bone marrow interface position prediction error. In Fig. 5(a), the errors are grouped according to the anatomical type of the specimen and whether the pseudo-plastic behavior was taken into account. In Fig. 5(b), the error data were grouped only according to the anatomical type or the consideration of the shear-thinning effect. The group medians were estimated as 1.58 (Thoracic PP−), 1.12 (Lumbar PP−), 1.40 (Thoracic PP+), 1.00 (Lumbar PP+), 1.48 (Thoracic PP7), 1.07 (Lumbar PP7), 1.12 (PP+) and 1.26 (PP−). The statistical tests indicate for both the principal and intermixed study groups differences between the subgroups ðpo0:00001Þ. The post-hoc analysis indicated mean rank differences between all tested subgroup-pairs.

4.

Discussion

The combination of a theoretical and a computational approach has been presented to describe the rheological behavior of PMMA-based bone cement flow through vertebral trabecular bone. An improvement, compared to the initial models (Baroud and Yahia, 2004; Teo et al., 2007a; Widmer and Ferguson, 2011), was achieved by the introduction of a multi-scale approach, i.e. the rigorous distinction and treatment of the physical phenomena at the pore and continuum length scale. The Reynolds number has been proposed to relate the pore-scale rheological and flow properties to the continuum length scale. Representative rheology data were obtained using a plate rheometer and applied to pre-compute the corresponding continuum-scale/Darcy viscosity by simulating the cement flow through representative specimens derived either from μCT scans or a cellular model of trabecular bone (Kim and Al-Hassani, 2002). This procedure facilitated the simulation of cement infiltration of cancellous bone. Cadaveric experiments of continuum-scale cement infiltration were conducted and compared to the computational results. The microscale rheology measurements revealed a high initial slope of the viscosity vs. time dependence that grades into a log-linear relationship. A pure linear time-dependence

51

as found by Baroud et al. (2003) could not be confirmed. However, this can be partially explained by the different formulations used during the experiments. The comparison of the cement distributions derived from the in-vitro and the computational simulations indicate an RMSE of the cement surface position on the order of 1.5 mm. In relation to a representative cement sphere diameter of 3.5 cm, this results in a relative prediction error of 5%. In an earlier study (Loeffel et al., 2008), the cement was injected into vertebral bone surrogates; the comparison to a computational model (Widmer and Ferguson, 2011), that had neglected the shearthinning behavior of the PMMA-based bone cement, indicated a typical overall error on the order of 15%. Although this level of error is acceptable, more attention should be drawn to the maximum prediction error which is 3.5 times higher than the RSME. It was found that in most cases the main error mode was in the direction of the cement delivery cannula(e), see also Fig. 4(b). The authors conclude therefore, that the localization and recovery of the needle tip position in the post-augmentation scan is the main contributor to the error. Teo et al. (2007a) did not perform any validation of their computational model, hence a comparison and a qualification of the herein presented model is not possible. The statistical tests and post-hoc analysis indicate that all study subgroups differ under the assumption of a per family error rate of PFE ¼ 0:05 (see also Fig. 5). The medians of the groups, where the shear-thinning effect was considered, were in all cases smaller compared to their counterparts. The improvement was moderate (7.87% overall RSME and 11.5% maximum reduction). The introduction of the proposed upscaling scheme transforms Darcy's law (Eq. (4)) into a nonlinear problem. However, the matrix-free approach (Kerkhoven and Saad, 1991) that was used demonstrated robustness and a good convergence rate. Therefore, the trade-off between computational complexity and accuracy improvement is in favor of the accuracy. The spinal level of the cadaver specimens had a significant effect on the prediction error, regardless of the effect of pseudo-plasticity. Although no significant systematic differences exist in the bone architecture along the T10 to L2 spine levels (Amling et al., 1996; Lochmüller et al., 2008), local variation might explain the dependence of the prediction accuracy on the anatomical site. The pore-scale flow analyzes were based on data derived from μCT data (T12 and L1 level) and a cellular rod model that is representative for lumbar cancellous bone. The experiments and continuum-scale simulations, however, were performed using T11 specimens. This mismatch of the anatomical site might imply that the proposed upscaling scheme is very sensitive to the micro-structure of the trabecular bone, or more generally, the solid phase of the porous structure. The coefficients in Table 1 however indicate that the viscosity is strongly correlated to the flow regime parameter Re and moderately correlated to the fluid type r. Only a weak correlation to the morphology of trabecular bone architecture (BV=TV, DA, Tb:Spðh;vÞ , Tb:Thðh;vÞ ), contrary to the finding above, was found. This must be considered as the most important finding of this study. The morphological parameters that were herein investigated are based on a rodlike architecture of the cancellous bone that is typical for the

52

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

lumbar levels. However, there is a gradual change to a more plate-like micro-architecture towards the cervical spine (Grote et al., 1995). The parameters that were used in this study might therefore not be adequate to reflect the differences in micro-architecture at all spine levels. The effect of the morphology on the viscosity is also reflected in Fig. 3. The statistically significant viscosity differences estimated for the cellular and the μCTmodel emphasize the importance of this aspect. The literature reports that Darcy's law is valid up to the transition from a creeping to a laminar flow regime, where inertial forces start to dominate the flow (Bear, 1972). This transition typically happens at Re ¼ 1. Fig. 3, however, suggests that this transitional regime starts one decade earlier (deviation from the pure log-linear trend). Two explanations, potentially amongst others, are proposed for this observation. First, the continuum-scale Reynolds number Re is an average estimate for a porous matrix, where the velocities and the pore-scale Reynolds number might locally differ by magnitudes. Consequently, a deviation from the log-linear trend can be expected for Re⪡1. Furthermore, the relationship to relate the pore to the continuum-scale Reynolds number in Eq. (11) might bias the prediction because of its informal foundation. However, the proposed upscaling scheme provides the appropriate degrees of freedom to compensate for this bias. In summary, we found that the pseudo-plastic properties associated with the PMMA-based bone cement should not be neglected in a computational model to simulate the infiltration of cancellous bone. The proposed upscaling scheme extends the Darcy law, which is adequate for Newtonian fluids, to non-Newtonian fluids. This finding, and our proposed solution, could be applied potentially to any general problem of a non-Newtonian fluid flowing through a porous media. The shear-rate dependency of the viscosity transforms the linear Darcy equation into a non-linear problem. The moderate gain of accuracy justifies the increased computational demand. It has been found that the pseudo-plastic rheology induces a viscosity change at the pore-scale on the order of one magnitude, whereas at the continuum-scale the viscosity change is on the order of 10%, which highlights the necessity for a multi-scale treatment of the problem.

Acknowledgments Funding for this research project was provided by the European Union (project FP7-ICT-223865-VPHOP), which is gratefully acknowledged. The authors thank Bert van Rietbergen and Javad H. Marangalou, Eindhoven University of Technology, for providing the bone scans within the VPHOP framework.

references

Amling, M., Pösl, M., Ritzel, H., Hahn, M., Vogel, M., Wening, V.J., Delling, G., 1996. Architecture and distribution of cancellous bone yield vertebral fracture clues. Archives of Orthopaedic

and Trauma Surgery 115, 262–269, http://dx.doi.org/10.1007/ BF00439050 ISSN 0936-8051. Baroud, G., Steffen, T., 2005. A new cannula to ease cement injection during vertebroplasty. European Spine Journal 14 (June (5)), 474–479, http://dx.doi.org/10.1007/s00586-0040822-1. Baroud, G., Yahia, F.B., 2004. A finite element rheological model for polymethylmethacrylate flow: analysis of the cement delivery in vertebroplasty. Proceedings of the Institution of Mechanical Engineers 218 (5), 331–338. Baroud, G., Wu, J.Z., Bohner, M., Sponagel, S., Steffen, T., 2003. How to determine the permeability for cement infiltration of osteoporotic cancellous bone. Medical Engineering & Physics 25 (May (4)), 283–288 ISSN 1350-4533. Baroud, G., Bohner, M., Heini, P., Steffen, T., 2004. Injection biomechanics of bone cements used in vertebroplasty. BioMedical Materials and Engineering 14 (4), 487–504. Baroud, G., Crookshank, M., Bohner, M., 2006. High-viscosity cement significantly enhances uniformity of cement filling in vertebroplasty: an experimental model and study on cement leakage. Spine (Phila Pa 1976) 31 (October (22)), 2562–2568, http://dx.doi.org/10.1097/01.brs.0000240695.58651.62. Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, New York ISBN 0-444-00114-X. Beaudoin, A.J., Mihalko, W.M., Krause, W.R., 1991. Finite element modelling of polymethylmethacrylate flow through cancellous bone. Journal of Biomechanics 24 (2), 127–136. Boger, A., Wheeler, K.D., Schenk, B., Heini, P.F., 2009. Clinical investigations of polymethylmethacrylate cement viscosity during vertebroplasty and related in vitro measurements. European Spine Journal 18 (September (9)), 1272–1278, http: //dx.doi.org/10.1007/s00586-009-1037-2. Bohner, M., Gasser, B., Baroud, G., Heini, P., 2003. Theoretical and experimental model to describe the injection of a polymethylmethacrylate cement into a porous structure. Biomaterials 24 (July (16)), 2721–2730 ISSN 0142-9612. Chui, C.-K., Teo, J., Wang, Z., Ong, J., Zhang, J., Si-Hoe, K.-M., Ong, S.-H., Yan, C.-H., Wang, S.-C., Wong, H.-K., Anderson, J.H., Teoh, S.-H., 2006. Integrative haptic and visual interaction for simulation of PMMA injection during vertebroplasty. Studies in Health Technology and Informatics 119, 96–98. Cotten, A., Dewatre, F., Cortet, B., Assaker, R., Leblond, D., Duquesnoy, B., Chastanet, P., Clarisse, J., 1996. Percutaneous vertebroplasty for osteolytic metastases and myeloma: effects of the percentage of lesion filling and the leakage of methyl methacrylate at clinical follow-up. Radiology 200 (2), 525–530. Cowin, S., Telega, J., 2003. Bone mechanics handbook. Applied Mechanics Reviews 56, 672–689. Darcy, H., 1856. Les fontaines publiques de la ville de Dijon. Dalmont, Paris. Doube, M., Klosowski, M.M., Arganda-Carreras, I., Cordelières, F.P., Dougherty, R.P., Jackson, J.S., Schmid, B., Hutchinson, J.R., Shefelbine, S.J., 2010. Bonej: free and extensible bone image analysis in imageJ. Bone 47 (December (6)), 1076–1079, http://dx. doi.org/10.1016/j.bone.2010.08.023. Farrar, D.F., Rose, J., 2001. Rheological properties of PMMA bone cements during curing. Biomaterials 22 (November (22)), 3005–3013. Gisep, A., Boger, A., 2009. Injection biomechanics of in vitro simulated vertebroplasty - correlation of injection force and cement viscosity. Bio-Medical Materials and Engineering 19 (January (6)), 415–420, http://dx.doi.org/10.3233/BME-20090607. Grote, H., Amling, M., Vogel, M., Hahn, M., Pösl, M., Delling, G., 1995. Intervertebral variation in trabecular microarchitecture throughout the normal spine in relation to age. Bone 16 (3), 301–308, http://dx.doi.org/10.1016/8756-3282(94)00042-5 ISSN 8756-3282.

journal of the mechanical behavior of biomedical materials 27 (2013) 43 –53

Heini, P.F., Wälchli, B., Berlemann, U., 2000. Percutaneous transpedicular vertebroplasty with PMMA: operative technique and early results. European Spine Journal 9, 445–450, http://dx.doi.org/10.1007/s005860000182 ISSN 09406719. Hulme, P.A., Boyd, S.K., Ferguson, S.J., 2007. Regional variation in vertebral bone morphology and its contribution to vertebral fracture strength. Bone 41 (December (6)), 946–957, http://dx. doi.org/10.1016/j.bone.2007.08.019. Huttenlocher, D., Klanderman, G., Rucklidge, W., 1993. Comparing images using the Hausdorff distance. IEEE Transactions on Pattern Analysis and Machine Intelligence 15 (9), 850–863. Kerkhoven, T., Saad, Y., 1991. On acceleration methods for coupled nonlinear elliptic systems. Numerische Mathematik 60, 525–548, http://dx.doi.org/10.1007/BF01385735 ISSN 0029599X. Kim, H.S., Al-Hassani, S.T.S., 2002. A morphological model of vertebral trabecular bone. Journal of Biomechanics 35 (August (8)), 1101–1114. Kirk, B., Peterson, J., Stogner, R., Carey, G., 2006. libMesh: a C++ library for parallel adaptive mesh refinement/coarsening simulations. Engineering with Computers 22 (December (3)), 237–254. Laredo, J.D., Hamze, B., 2004. Complications of percutaneous vertebroplasty and their prevention. Skeletal Radiology 33, 493–505, http://dx.doi.org/10.1007/s00256-004-0776-8 ISSN 0364-2348. Lian, Z., Chui, C.-K., Teoh, S.-H., 2008. A biomechanical model for real-time simulation of PMMA injection with haptics. Computers in Biology and Medicine 38 (March (3)), 304–312, http://dx.doi.org/10.1016/j.compbiomed.2007.10.009. Liebschner, M.A., Rosenberg, W.S., Keaveny, T.M., 2001. Effects of bone cement volume and distribution on vertebral stiffness after vertebroplasty. Spine (Phila Pa 1976) 26 (July (14)), 1547–1554. Lochmüller, E.-M., Pöschl, K., Würstlin, L., Matsuura, M., Müller, R., Link, T.M., Eckstein, F., 2008. Does thoracic or lumbar spine bone architecture predict vertebral failure strength more accurately than density. Osteoporosis International 19 (April (4)), 537–545, http://dx.doi.org/10.1007/s00198-007-0478-x. M. Loeffel, 2007. Computer Assisted High Pressure Cement Injection in Spinal Interventions. PhD Thesis, University of Bern.

53

Loeffel, M., Ferguson, S.J., Nolte, L.-P., Kowal, J.H., 2008. Vertebroplasty: experimental characterization of polymethylmethacrylate bone cement spreading as a function of viscosity, bone porosity, and flow rate. Spine 33 (May (12)), 1352–1359, http://dx.doi.org/10.1097/BRS.0b013e3181732aa9. Murdock, J.W., 1976. Fluid Mechanics and its Applications. Houghton Mifflin Company. Pearson, J.R.A., Tardy, P.M.J., 2002. Models for flow of nonNewtonian and complex fluids through porous media. Journal of Non-Newtonian Fluid Mechanics 102 (2), 447–473, http://dx. doi.org/10.1016/S0377-0257(01)00191-4 ISSN 0377-0257. URL 〈http://www.sciencedirect.com/science/article/pii/ S0377025701001914〉. Pinder, G.F., Gray, W.G., 2008. Essentials of Multiphase Flow and Transport in Porous Media. John Wiley & Sons, Inc.. Teo, J., Wang, S.C., Teoh, S.H., 2007a. Preliminary study on biomechanics of vertebroplasty: a computational fluid dynamics and solid mechanics combined approach. Spine 32 (May (12)), 1320–1328, http://dx.doi.org/10.1097/ BRS.0b013e318059af56. Teo, J., Si-Hoe, K., Keh, J., Teoh, S., 2007b. Correlation of cancellous bone microarchitectural parameters from microCT to CT number and bone mechanical properties. Materials Science and Engineering: C 27 (2), 333–339 ISSN 0928-4931. Tschirhart, C.E., Roth, S.E., Whyne, C.M., 2005. Biomechanical assessment of stability in the metastatic spine following percutaneous vertebroplasty: effects of cement distribution patterns and volume. Journal of Biomechanics 38 (August (8)), 1582–1590, http://dx.doi.org/10.1016/j.jbiomech.2004.07.023. Widmer, R.P., Ferguson, S.J., 2011. A mixed boundary representation to simulate the displacement of a biofluid by a biomaterial in porous media. Journal of Biomechanical Engineering 133 (5), 051007, http://dx.doi.org/10.1115/ 1.4003735 URL 〈http://link.aip.org/link/?JBY/133/051007/1〉. Widmer, R.P., Ferguson, S.J., 2012. On the interrelationship of permeability and structural parameters of vertebral trabecular bone: a parametric computational study. Computer Methods in Biomechanics and Biomedical Engineeringhttp://dxdoi.org/ 10.1080/10255842.2011.643787 URL 〈http://www.tandfonline. com/doi/abs/10.1080/10255842.2011.643787〉.