Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine http://pih.sagepub.com/

Development and validation of a distal radius finite element model to simulate impact loading indicative of a forward fall Timothy A Burkhart, Cheryl E Quenneville, Cynthia E Dunning and David M Andrews Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine published online 10 February 2014 DOI: 10.1177/0954411914522781 The online version of this article can be found at: http://pih.sagepub.com/content/early/2014/02/07/0954411914522781 A more recent version of this article was published on - Mar 13, 2014

Published by: http://www.sagepublications.com

On behalf of:

Institution of Mechanical Engineers

Additional services and information for Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine can be found at: Email Alerts: http://pih.sagepub.com/cgi/alerts Subscriptions: http://pih.sagepub.com/subscriptions Reprints: http://www.sagepub.com/journalsReprints.nav Permissions: http://www.sagepub.com/journalsPermissions.nav

Version of Record - Mar 13, 2014 >> OnlineFirst Version of Record - Feb 10, 2014 What is This?

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

Original Article

Development and validation of a distal radius finite element model to simulate impact loading indicative of a forward fall

Proc IMechE Part H: J Engineering in Medicine 1–14 Ó IMechE 2014 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/0954411914522781 pih.sagepub.com

Timothy A Burkhart1, Cheryl E Quenneville2, Cynthia E Dunning1 and David M Andrews3

Abstract The purpose of this work was to develop and validate a finite element model of the distal radius to simulate impact loading. Eight-node hexahedral meshes of the bone and impactor components were created. Three separate impact events were simulated by altering the impact velocity assigned to the model projectile (pre-fracture, crack and fracture). Impact forces and maximum and minimum principal strains were calculated and used in the validation process by comparing with previously collected experimental data. Three measures of mesh quality (Jacobians, aspect ratios and orthogonality) and four validation methods (validation metric, error assessment, fracture comparisons and ensemble averages) assessed the model. The element Jacobians, aspect ratios and orthogonality measures ranged from 0.08 to 12, 1.1 to 26 and 270° to 80°, respectively. The force and strain validation metric ranged from 0.10 to 0.54 and 0.35 to 0.67, respectively. The estimated peak axial force was found to be a maximum of 28.5% greater than the experimental (crack) force, and all forces fell within 62 standard deviation of the mean experimental fracture forces. The predicted strains were found to differ by a mean of 33% across all impact events, and the model was found to accurately predict the location and severity of bone damage. Overall, the model presented here is a valid representation of the distal radius subjected to impact.

Keywords Radius, impact, validation metric, forward fall

Date received: 7 May 2013; accepted: 15 January 2014

Introduction While experimental testing on human participants and cadaveric specimens provides the most realistic response of the human body and the most accurate injury patterns, these types of testing are not always the most feasible.1–4 For example, the safety of participants limits in vivo testing to submaximal loads, and the destructive nature of damage tests on cadaveric specimens can become costly.1 In comparison, finite element models (FEMs) provide a feasible and practical alternative for calculating the stress and strain response of bone under a variety of loading conditions.3 Obtaining accurate outputs from a FEM is dependent on the development of a high-quality mesh.5 Given the complex geometry of long bones such as the radius, elements with large distortions often occur and are potential sources of low accuracy and solution instability.6 It has also been suggested that, in general,

a FEM composed entirely of hexahedral elements will lead to more accurate results, is more computationally efficient, and mimics bone geometry better than a model meshed with shell, tetrahedral or voxel elements.7,8 While several researchers have developed FEMs of the distal radius,1,9–19 they have not included

1

Department of Mechanical and Materials Engineering, Western University, London, ON, Canada 2 Department of Mechanical Engineering, McMaster University, Hamilton, ON, Canada 3 Department of Kinesiology, University of Windsor, Windsor, ON, Canada Corresponding author: Timothy A Burkhart, Department of Mechanical and Materials Engineering, Western University, 1151 Richmond Street, London, ON N6A N9B, Canada. Email: [email protected]

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

2

Proc IMechE Part H: J Engineering in Medicine

measurements of the mesh quality, have provided little to no information about how they were validated (a limitation not exclusive to radius bone models), and have generally been composed of a combination of element shapes. Furthermore, the static loading conditions used in past studies (those listed above) do not make the findings applicable to situations where dynamic impact loading occurs (i.e. arresting a forward fall with the distal upper extremity). These models have also not considered the kinematics of a forward fall, such as the angle the forearm makes with the ground or the angle of the wrist at impact.20 Therefore, the purpose of this study was to describe the development of a human distal radius FEM used to simulate dynamic impact loading indicative of a forward fall. In addition, a thorough assessment of the model’s mesh quality is provided, and its validity is established using several metrics.

Methods Experimental protocol The experimental protocol used to validate the model has previously been described in detail.20,21 Briefly, eight cadaveric radius specimens (Table 1) were potted in poly(vinyl chloride) (PVC) pipe using DenstoneÒ cement (Modern MaterialsÒ, Heraeus Holding GmbH, Hanau, Germany) at an angle of 75° (Figure 1(b)), the most commonly reported angle that the forearm makes with the impact surface during a forward fall–initiated impact.22–25 The potted specimens were placed within a pneumatically controlled impact system such that a 6.8-kg projectile impacted a distal bracket, transferring the force through a load cell and model lunate/scaphoid onto the articular surface of the radius (Figure 1(b)).26 A model lunate/scaphoid was used experimentally to avoid fractures to the carpal bones. Prior to potting, the contours of the lunate and scaphoid were aligned with the respective fossae on the articular surface of the radius using a custom-designed potting jig. Also, to further ensure accurate radius–carpal alignment, the position of the model lunate/scaphoid on the

load cell was fully adjustable. The velocity of the projectile was controlled by adjusting the pressure in the acceleration tube. Pilot testing revealed that an impact occurring at 20 J (2.4 m/s) was low enough not to cause any damage to the radius specimens. Therefore, each specimen experienced an initial 20-J impact (pre-fracture impact event) followed by subsequent impacts increasing in 10-J increments until the specimen was fractured into at least two distinct fragments (fracture impact event). Following each impact, the specimens were carefully inspected to ensure that no damage had occurred before moving onto the next impact. Finally, prior to fracture, if an impact resulted in external nonpropagating damage, the trial was categorized as a crack impact event. Each specimen was subjected to a mean (standard deviation (SD)) of 3.5 (1.2) and 4.2 (1.7) impacts to experience a crack and fracture event, respectively.20 Prior to experimental testing, the locations of three strain gauge rosettes (Vishay MicroMeasurements, Vishay Precision Group, Malvern, PA, USA; grid resistance = 350 O; gauge factor = +1.3%) and two accelerometers (MMA220KEG (x- and y-axes) and MMA1210 (z-axis); Freescale Semiconductor, Ottawa, ON, Canada) were identified on each of the specimens by digitizing the center of each transducer (Microscribe G2X; Immersion Corporation, San Jose, CA, USA). Three screw holes placed in the PVC pipe potting device were also digitized to create a local coordinate system. Transformations were conducted to locate the x–y–z coordinates of the center of each strain gauge and accelerometer on the solid model. This allowed elements in these regions of the mesh to be selected so that strains and accelerations could be directly compared to those measured experimentally.

Mesh development The eight potted specimens were scanned using computed tomography (CT; GE LightSpeed VCT; General Electric Healthcare, Chalfont, St Giles, UK) at 120 kVp, 100 mA, with a slice thickness of 0.625 mm. A section of hydroxyapatite was also included in the scan volume to assist with material registration when

Table 1. A summary of the ranking system that was used to select the most representative radius specimen for FEM purposes. Specimen

Age

Rank

BMD

Rank

Frykman score

Rank

Fracture force

Rank

Sum

Overall rank

07017L 07030R 08007L 07043L 08010L 07007R 04011R 07004L Mean

76 54 60 68 46 57 54 73 61

7 3 1 3 7 2 3 6

0.39 0.58 0.48 0.42 0.47 0.47 0.53 0.54 0.48

8 7 1 6 2 3 4 5

7 7 7 7 7 7 3 3 6

1 1 1 1 1 1 7 7

2021 4340 1493 2519 1941 1980 3848 1627 2471

2 8 6 1 4 3 7 5

26 23 9 15 22 11 25 29

7 5 1 3 4 2 6 8

The bone mineral density (BMD) is taken from the distal third of the radius determined from the Dual-energy X-ray absorptiometry (DXA) scans.

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

Burkhart et al.

3

Figure 1. (a) Finite element meshes and (b) experimental setup of the radius bone and impactor components. The impactor exits the acceleration tube from the left contacts the impact plate, and the force is transferred through the load cell and onto the articular surface of the radius. Also shown are the locations of the three accelerometer axes (a) proximally and distally (modified with permission from Burkhart et al.20,21).

extracting the surface geometry. From the eight scanned radius specimens, one was chosen as a representative specimen from which a single FEM was built. Selection of the most representative sample was based on ranking the specimens according to donor age, bone mineral density of the distal third of the radius, the fracture force, and the Frykman fracture score. Specimens were ranked on each of these variables according to how well they agreed with the mean values of each variable (Table 1).4 The CT scan files of the representative specimen were imported into MimicsÒ (Materialise, Leuven, Belgium) medical imaging software to segment and build a solid model of the cortical, cancellous and marrow regions of the bone. The threshold of the densest region of bone (i.e. the cortical bone) was determined in comparison to the section of hydroxyapatite that was included in the scan. This resulted in Hounsfield units from 621 to 2693 for the cortical bone and \ 621 for the cancellous bone. These values compare well with previously reported data.15 The transition from cancellous bone to bone marrow was not clearly defined (i.e. fragments of cancellous bone found throughout the intermedullary canal); therefore, as an assumption, this transition was identified as the axial location where the cavity was not entirely occupied by cancellous bone.

Surface geometries were extracted from the threedimensional solid bone models in the form of stereolithography (STL) files and were imported into TrueGridÒ meshing software (XYZ Scientific Inc. Livermore, CA, USA). Here, eight-node solid hexahedral elements were used to generate volumetric meshes of the three bone components (compact bone, cancellous bone and bone marrow) and all components of the experimental impact system (Figure 1(a)). The geometrical accuracy of the models was retained through the projection method that TrueGrid relies on, in which the nodes of blocks of mesh are directly attached to the surface of the STL models. Hexahedral elements were selected over tetrahedral elements, as they are considered to be more accurate in dynamic simulations.5,27–30 This is a result of tetrahedrals being inherently stiff, incompressible and predisposed to mesh locking.

Material models and properties A nonlinear elastic–plastic material model was selected for the cortical and cancellous bone, while a linear elastic model was chosen to represent the bone marrow. Bone material properties were selected from the literature and, when possible, were radius, age and sex specific (Table 2). All bone components were assigned a

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

4

Proc IMechE Part H: J Engineering in Medicine

Table 2. Material properties assigned to each of the model components. Part

Number of elements

Mass density (kg/m3)

Rigid part volume (mm3)

Elastic modulus (GPa)

Plastic modulus (GPa)

Poisson’s ratio

Yield stress (GPa)

C (s21)

p

Cortical Cancellous Marrow Lunate/scaphoid Projectile Pot/cement Proximal bracket Distal bracket Foam Rear bearing Front bearing

55,800 180,000 697,500 40,000 5000 37,800 6400 36,000 8000 2688 1000

1150a 970d 1070g 1060 930 1.0E+4 205 220 52 5.0E+6 5.0E+6

– – – – 7302.0 253.9 20945.1 19,489 – – –

25.1b 1.8d 0.02g 3.0 210.0 210.0 210.0 210.0 0.008 210.0 210.0

1.25 0.09e – – – – – – – – –

0.3c 0.3c 0.49g 0.38 0.3 0.3 0.3 0.3 0.3 0.3 0.3

0.124 b 0.001f – – – – – – – – –

64.7 64.7 – – – – – – – – –

6.4 6.4 – – – – – – – – –

a

Schonenau et al.31 Burstein et al.32 c Boutroy et al.16 d Troy and Grabiner.15 e Imai et al.33 f Kim et al.34 g Peng et al.35 b

constant stress element formulation, and while this element formulation is computationally efficient, it is prone to hourglassing (i.e. deformations that occur in a state of zero stress). Therefore, hourglass control was implemented in the model such that small internal forces are applied to avoid negative volume elements from occurring. To control for strain-rate effects, the Cowper– Symonds strain-rate formulation was selected for the three types of bone. The Cowper–Symonds model scales the static yield stress (Figure 2(a)) according to the factor36  1p e_ 1+ C

ð1Þ

where e_ is the strain rate, and C and p are materialspecific parameters that are determined experimentally. For the current model, bone-specific C and p parameters were defined based on the experimental data set presented by McElhaney (1966).37 The Cowper– Symonds model (equation (1)) can be rewritten as  0 p so 1 ð2Þ e_ = C so where s0o is the dynamic strain rate–specific yield stress and so is the static yield stress. Equation (2) can be written as  0  s ln e_ = p ln o  1 + ln C ð3Þ so which represents the equation of a straight line. The above forms of the yield stress and strain rates are then plotted against each other (Figure 2(b)), and the C and p parameters can be determined, such that ln C is the

intercept and p is the slope of the regression line describing the above relationship. Based on this data set, ln C = 4.176 s21, C = 65.1 s21 and p = 6.4. All of the above parameters were selected based on the results of a parametric investigation reported previously.38 To measure model bone surface strains in the model, sets of elements representing the centers of the three experimental strain gauges were created (Figure 3). Model accelerations were also computed at the locations that correspond to the experimental accelerometers on the bone (Figure 1(a)). The locations of the strain gauges and accelerometers were identified on the bones prior to experimental testing by digitizing (Microscribe G2X; Immersion Corporation, San Jose, CA, USA) the center of each strain gauge and accelerometer. Three screw holes placed in the PVC potting were also digitized (and were visible in the CT scans) to create a local coordinate system. Transformations were conducted, and the x–y–z coordinates of the center of each strain gauge (Figure 3) and accelerometer (Figure 1(a)) were located on the solid model and exported to a text file. This allowed elements in these regions to be selected so that model strains and accelerations could be compared to those measured experimentally. Also, given the angle of the bone within the impactor and the orientation of the accelerometers experimentally, it was necessary to define local coordinate systems within the model at the locations of the accelerometers (Figure 1(a)). This ensured that the accelerations from the model were calculated along axes that corresponded to the experimental accelerometers. The foam that was placed between the projectile and distal bracket (Figure 1(a)) (used to ensure an adequate impulse duration experimentally) was modeled as having a low density, as it has been shown to have a good correlation with experimental results and is appropriate

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

Burkhart et al.

5

Figure 2. (a) A series of strain-rate (s21) dependent stress– strain curves calculated using the Cowper–Symonds strain-rate effect model whose coefficients (C and p) are calculated from (b) the transformed McElhaney37 (1966) data.

Figure 3. Locations of the experimental strain gauges in comparison to the location of the numerical strain calculations.

for simulating highly compressible foams (i.e. a foam that is typically used to simulate seat cushion foams).39,40 To avoid unrealistically large deformations, compared

to what was seen experimentally, an internal contact force was included in the definition of the foam material. The modulus of the foam (Table 2) was determined from experimental compression tests (InstronÒ materials testing machine, model 8872, Canton, MA, USA). The carpal bones, specifically the lunate and scaphoid, were composed of high-density polyethylene (SawBonesÒ; Pacific Research Laboratories Inc., Vashon, WA, USA). While polyethylene undergoes relatively small deformations when compressed, it was modeled as an elastic material so that three force components (Fx: medial–lateral, Fy: inferior–superior and Fz: axial (Figure 1(a))) could be calculated for the node set representing the load cell. The impact forces were calculated at a cross section that was defined through the second layer of elements of the carpal bones. This location was chosen given its close proximity to the experimental load cell. A node set was formed that defined the nodes through the cross section, and an element set was constructed that included all of the elements to one side of the cross-section node set. The forces from the elements in the element set were then summed to compute the cross-section forces. Finally, the remaining five impactor components (projectile, proximal and distal brackets and front and rear bearings; Figure 1) were modeled as rigid materials and were assigned the material properties of stainless steel. The masses of the model impactor components were matched to the corresponding experimental components by calculating the volumes of each rigid material (measured in LS-PrePostÒ, Livermore Software Technology Corp., Livermore, CA) and adjusting the densities assigned to each (Table 2). A number of boundary conditions were required to ensure that the motion of the model impactor components matched those seen experimentally. To constrain the motion of the distal bracket within the x–y plane, a planar joint was defined between the distal bracket and the distal bearing. This required creating a set of nodes from both parts that was used to define the direction of the normal force. Two separate three-node sets were then used to form the coordinate system within each of the parts. A force was applied to the distal bracket, simulating the frictional force between the distal bracket and linear rail that results from the large bending moment resisted by the bearings. A coefficient of friction of 0.5 was used, based on the data from a previous experimental investigation using the impactor described here.26 The motion of the projectile and the proximal bracket were also constrained to the z-direction. An initial velocity was applied to the projectile that was varied to mimic the pre-fracture (2.1 m/s), crack (2.7 m/s) and the fracture (3.3 m/s) impact events, resulting in three impact event simulations.

Mesh quality assessment Several mesh quality assessment techniques were used to evaluate the quality of the mesh prior to model simulations.5 First, the element Jacobian was examined,5,41

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

6

Proc IMechE Part H: J Engineering in Medicine

utilizing an assessment criterion of being a positive value, preferably greater than 0.2 and with less than 5% falling below 0.7.5,42,43 The orthogonality of the mesh was also assessed by measuring each of the three angles formed by the vertices at each of the eight nodes for all elements.5 The aim was to avoid absolute deviations greater than 70°.4,5 In some cases, this is not always feasible and a secondary orthogonality criterion was used, such that less than 5% of the orthogonality measures were allowed to deviate from the above value.5,43 Finally, the last method of mesh quality assessment was the calculation of element aspect ratios.5 The most numerically accurate solution is found when the edges of the elements are equal in length.5 Therefore, the goal was to achieve aspect ratios as close to one as possible while avoiding ratios greater than 10 and with less than 5% exceeding 3.4,5,42,43

Energy balance assessment The global energy balance of the model was analyzed to make certain that there were no major inconsistencies in the energy of the system. This was achieved by ensuring that the total energy equals the sum of internal, kinetic, sliding, hourglass, system damping and rigid wall energies.5,44 Since hourglass control was implemented, a specific analysis of the hourglass energy was conducted to confirm that it did not contribute more than 10% to the total energy.5,45,46 Similarly, the sliding energy was also investigated, as negative sliding energies are undesirable.5,44 Finally, a mass scaling technique was employed, and, therefore, an analysis of the total added mass to the system was conducted.5 The mass added specifically to the bone components was also monitored to ensure that excessive mass was not added that could affect the structure and stiffness of these components.5

Model and simulation validation To assess the ability of the model to accurately predict the response of the distal radius to impact loading, four different validation methods (validation metric, percentage errors, fracture comparisons and ensemble averages) were used to compare the numerical results to those measured experimentally. The validation metrics, percentage errors and fracture assessments were determined by comparing the model results to the experimental findings of the specific specimen that was used to develop the model, while the time domain curves were compared to the ensemble averages calculated from all eight experimental specimens. The agreement between the experimental and model force, acceleration and strain data were assessed over the entire impact time interval (;0.2 s) using a validation metric proposed by Oberkampf and Trucano.5,47 The validation metric is computed using

  N y(tn )  Y(tn ) 1X   V=1  tanh  N n=0 Y(tn )

ð4Þ

where V is the validation metric, N represents the total number of samples, tanh is the hyperbolic tangent trigonometric function, y(tn) is the numerical measurement of the dependent variable at time t and Y(tn) is the experimental measurement of the dependent variable at time t. The major advantage of this metric is that it measures the agreement between experimental and numerical results in a way that positive and negative errors cannot cancel each other out. The validation metric produces a value of 1 when there is perfect agreement between experimental and numerical results and approaches 0 as the differences increase.5,48 When calculating the validation metric, the model and experimental signals were aligned at approximately 2 ms prior to the impulse onset.5 The percentage difference between the model and experimental results was calculated for the peak force, load rate, impulse, and impulse duration for all three force components, and the peak acceleration and acceleration rate in the axial (parallel with the long axis of the radius) and off-axis (perpendicular to the long axis of the radius) directions (Figure 1(a)), for both the distal and proximal accelerometer, and the maximum and minimum principal strains at all three strain gauge locations were also calculated. The percentage differences were calculated for all three impact events (pre-fracture, crack and fracture). The third validation method used the von Mises stresses to visually compare the location and severity (an increase in volume of failed elements would indicate a more severe fracture) of bone failure during the three impact events.5 Numerically, an element was assumed to have failed when the calculated stresses exceeded a pre-determined critical stress level (134 and 5.3 MPa for the cortical and cancellous bone, respectively).42 Although these ultimate stress values were initially determined from the lower extremity, these data do not currently exist specifically for the distal radius. Finally, to provide an assessment of the generalizability of the model, data corridors (i.e. ensemble averages) were developed based on the entire sample of experimental data.5,49,50 As force is the most commonly presented validation variable, data corridors of the Fx, Fy and Fz forces were created by aligning all of the experimental signals (within each force axis) at the time of impulse onset and the mean 6 2 SD were calculated and represented the upper and lower boundaries of the data corridors. The simulation results are considered to be a good representation of the sample if they fall within the upper and lower boundaries.5,49,50

Results Overall, the entire model was composed of approximately 1.04 million elements and 1.13 million nodes

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

Burkhart et al.

7

Table 3. Energy balances and validation metrics of the radius model across the three impact simulations.

Energy balance (J) Total energy Hourglass energy (% total energy) Kinetic energy Internal energy Damping Sliding energy Total mass added (% of total mass) Validation metrics Force Fx Fy Fz Acceleration Distal axial Distal off-axis Proximal axial Proximal off-axis Strain Gauge 1 Maximum principal Minimum principal Gauge 2 Maximum principal Minimum principal Gauge 3 Maximum principal Minimum principal

Pre-fracture

Crack

Fracture

15.2 1.1 (7.2) 12.5 1.2 0 0.2 2.3 (1.2E203)

25.8 2.1 (8.1) 22.1 2.0 0 0.4 2.7 (1.3E203)

38.5 3.2 (8.3) 29.9 3.3 0 2.2 3.0 (1.4E203)

0.10 0.34 0.43

0.21 0.22 0.54

0.22 0.24 0.23

0.40 0.42 0.37 0.36

0.37 0.35 0.30 0.36

0.40 0.32 0.37 0.39

0.56 0.34

0.58 0.48

0.34 0.39

0.63 0.67

0.42 0.41

0.35 0.39

0.63 0.67

0.42 0.41

0.35 0.39

(Table 2). The Jacobians of the radius mesh ranged from 0.08 to 12 and 0.18 to 8 in the interior (cancellous and marrow regions) and the exterior (cortical region), respectively. The deviation in the cancellous marrow mesh ranged from 270° to 85° (corresponding to absolute interior angles of 20°–175°), while the cortical mesh had deviations ranging from 270° to 80° (corresponding to absolute interior angles of 20° to 170°). While some angles were found to deviate more than 70°, the mean deviation was zero in all of the meshes and less than 50 (out of more than 1.04 million) elements were found to violate this rule. The aspect ratios for the cancellous and marrow mesh section ranged from 1.2 to 12.4, while aspect ratios between 1.1 and 26 were found for the cortical mesh. While the maximum aspect ratio for the cortical bone was far greater than the cut-off value of 10, only one element was shown to have this value, and it was located away from the area of interest. The remaining elements had aspect ratios that ranged from 1.1 to 15 (;90% of elements with an aspect ratio \ 10). Finally, for each mesh assessment criterion, the violating elements accounted for less than 5% of all bone elements (cortical, cancellous and marrow). Both the total global and hourglass energies increased across the three impact simulations (pre-fracture, crack and fracture). However, the hourglass energy only accounted for a maximum of 8.3% of the total energy (fracture simulation; Table 3). Across all three impact events, the sum of the total energy was balanced by the individual energy contributions

(Table 3). Similarly, the maximum added mass (3.0 kg) accounted for a very small proportion of the total mass (1.4E203%) (Table 3), and only 3% of this was added to the bone components. The validation metrics ranged from 0.10 (Fx prefracture) to 0.67 (strain gauges 2 and 3 minimum principal strain; Table 3). Overall, the strain gauge and force data resulted in the highest (0.46 (0.12)) and lowest (0.28 (0.13)) mean (SD) validation metrics, respectively. The greatest differences between the model and experimental peak forces were consistently found along the x-axis (medial–lateral), such that the model tended to overestimate the peak force by a maximum of approximately 197%. While the model also overestimated peak forces along the y- and z-axes, they tended to agree relatively well with the experimental peak forces during the pre-fracture and crack impact events. The estimated peak forces along the y-axis were approximately 4.8%, 3.5% and 58.9% greater than the experimental y-axis forces, and the peak model forces along the z-axis were found to be 19.1%, 28.5% and 64.2% greater than the peak forces determined from the experimental testing for pre-fracture, crack and fracture events, respectively (Figure 4(a), (c) and (e)). While the model Fx impulse overestimated the experimental impulse by a maximum of 170%, relatively good agreement was found for the Fy (maximum: 63.6%, pre-fracture) and Fz (maximum: 24.8%, fracture) impulses. Finally, small differences were found for the duration of the impulse across all three

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

8

Proc IMechE Part H: J Engineering in Medicine

Figure 4. Comparison of (a, c and e) the model and experimental force and (b, d and f) strain data for the (a and b) pre-fracture, (c and d) crack and (e and f) fracture impact events.

force axes and impact events (maximum: 38%, Fx, fracture). A comparison of the model forces to the mean experimental forces (an ensemble average calculated across all eight specimens) suggests a relatively good relationship, with the model forces falling

consistently within 62 SD in both the Fy and Fz directions across the pre-fracture and crack impact events (Figure 5). Across all the force axes, the model fracture event forces fell outside of 2 SD of the experimental ensemble averages, and are therefore not shown here.

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

Burkhart et al.

9

Figure 5. Comparison of the model (a and b) Fx, (c and d) Fy and (e and f) Fz forces (solid black line) to the mean 6 2 SD (dotted black lines) for the (a, c and e) pre-fracture and (b, d and f) crack experimental impact events. Data from the fracture impact events are not shown here.

Overall, the model was capable of predicting the peak accelerations relatively well. Distally, the peak acceleration differences ranged from 1.8% (offaxis, crack) to 231.7% (off-axis, pre-fracture), while the differences in the peak proximal accelerations ranged from 22.4% (off-axis, fracture) to 262.9% (axial, crack). When the model and experimental

acceleration rates were compared, a range of 3.9% (distal off-axis, pre-fracture) to 102.8% (distal axial, crack) was found. The maximum and minimum principal strains measured from the site of gauge 1 (radial styloid) and gauge 2 (ulnar side of radius) in the model were found to differ by a mean of approximately 33% across all

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

10

Proc IMechE Part H: J Engineering in Medicine

impact events, when compared to the strains recorded experimentally (Figure 4(b), (d) and (f)). The smallest strain differences occurred at the site of gauge 1 during the crack event (3.0%), while the largest differences were recorded at the site of gauge 2 during the fracture event (75%). Relatively large differences were found at the site of gauge 3 (proximal), ranging from 29.6% (maximum principal strain, pre-fracture) to 430% (minimum principal strain, fracture). Finally, during the pre-fracture impact event, there were no elements in the cortical bone that exceeded the critical limit, while a cluster of cancellous elements located in the volar region on the ulnar side of the radius indicated that ‘‘damage’’ had occurred. Similarly, no cortical ‘‘damage’’ was noted during the crack event, while there was an increase in the volume of elements in the cancellous bone that exceeded the critical limit (Figure 6). However, during the fracture event, the cortical bone showed areas of ‘‘damage’’ on the volar aspect of the intra-articular surface on the ulnar side of the radius as well as into the sigmoid fossa (Figure 6). The cancellous bone also showed significant signs of ‘‘damage’’ during the fracture event (Figure 6). The pattern of ‘‘damage’’ seen during the fracture event for the cortical bone and during all events for the cancellous bone corresponded very well with the damage that was observed experimentally; this was evident not only in the specimen that was used to generate the model (Figure 6), but across all tested specimens as well.20

Discussion To date, this is one of the first FEMs of the distal radius composed entirely of hexahedral elements that has simulated dynamic impact loading. Taking into consideration all of the mesh quality assessments, and the entirety of the validation results (including the validation metrics, the comparison to specimen-specific experimental results, mean experimental results and comparisons of fracture patterns), the model presented here is a good representation of the distal radius, and it responds accurately to the application of loads consistent with a forward fall arrest. While previous studies have used FEMs of the radius to predict the response to loading, none have included an assessment of the mesh quality. The mesh diagnostics presented here, when viewed collectively across all three bone components, suggest that the current mesh is of relatively high quality.5 While there were some elements that did not meet all of the diagnostic criteria, these elements were not located in areas of high interest (i.e. the articular surface of the radius and the distal radial diaphysis) and accounted for less than 0.5% of all the mesh elements.5 It is not expected that these violating elements would have a significant effect on the stresses and strains that were calculated in and around the articular region. Although the final bone meshes presented here were highly discretized, consisting of over 900,000 elements,

a mesh density sensitivity analysis was not conducted (and therefore, it is possible that a coarser mesh would have performed as well). The difficulties in altering the mesh density in a manually generated mesh for both tetrahedral and hexahedral elements have often been noted.5,51 Given the high areas of curvature, the complex geometry and the thin cortical shell associated with the radius, it becomes even less feasible to perform a mesh sensitivity analysis on a manually generated hexahedral mesh.5 Furthermore, the multipart butterfly method used in TrueGrid creates nodal compatibility issues between merged faces within and between parts that are not conducive to efficiently changing the number of elements in the mesh once established. Nodal compatibility issues also make it difficult to manually adjust the shape of individual elements. However, the mesh assessment methods used here dictated the coarseness of the mesh to some degree, by systematically increasing/decreasing the number of elements to optimize the mesh quality metrics.5 Specifically, it is important that the Jacobians are, at a minimum, positive, as negative Jacobians are representative of negative volume elements, and this will prevent the simulations from running.5 This protocol agrees with previously published dynamic hexahedral meshes of the lower extremity4,42 (specifically the tibia), and a recent review suggests that this is an acceptable alternative to performing a mesh sensitivity analysis.5 In finite element simulations, a phenomenon known as hourglassing can occur where an element undergoes a deformation in the absence of strain.5 Hourglassing can lead to inaccurate results and, in severe cases, can result in negative volume elements. To overcome hourglassing, a stiffness form hourglass control was implemented in which a small elastic stiffness was added, allowing the elements to resist hourglassing. However, inaccuracies can also result from the addition of too much hourglass energy. In the simulations reported here, the hourglass energy increased across the three impact events, but remained below the recommended 10% of total global energy.5,39,45,46 This, combined with the energy balance findings and the small percentage of added mass (including only a small increase in the mass of the bone components), suggests that the results of the model were not affected by the inclusion of hourglass control and provides strong verification of the model.5 Overall, the model presented here is a good representation of the distal radius, and the impact loads agreed relatively well with the experimental findings, primarily during the pre-fracture and crack impact events, when no initial damage is present in the bone. The Fy and Fz loads predicted by the model during the pre-fracture and crack events agreed well with the specimen-specific experimental loads, and the differences between the model and experimental findings may be partially attributed to the selection of material properties. The elastic modulus of cortical and cancellous bone used

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

Burkhart et al.

11

Figure 6. Maximum von-Mises stresses that exceeded the critical level during the (a and b) pre-fracture, (c and d) crack and (e and f) fracture events in the (a, c and e) cortical and (b, d and f) cancellous bone regions and a comparison of the damage to the (g) sigmoid fossa and (h) intra-articular reported from the experimental testing. Elements that were predicted to have failed (cortical = 134 MPa; cancellous = 5.3 MPa) are highlighted black.

here were taken from previous works,15,31–33 and while care was taken in selecting these values (i.e. values were selected to match the age and sex of the specimen donor), they are still estimates of the material properties of the specific human bone used here. In addition, although homogenous material mapping was used for the current investigation, it has been shown by Taddei et al.52 that applying CT-based inhomogenous material properties has a minimal effect on the accuracy of strain predictions for simulating slow loading rates. It is the opinion of the current authors that differences between inhomogenous and homogenous material mapping would be minimal when simulating for impact loading to predict fracture.

The peak experimental forces along the x-axis were negligible when compared to the model, resulting in the greatest error in the model outputs. This error can be attributed to the procedure that was used to create the mesh of the carpal bones. This specific mesh was created based on a highly smoothed version of the scanned lunate/scaphoid and resulted in a continuous mesh (i.e. void of any holes that were present from the scanning). However, although great care was taken to ensure that the carpals were properly aligned with the radius, the surface geometry may have been altered enough to lose some of the detailed morphology that is helpful in aligning the lunate/scaphoid with the radius. Although the irregular anatomical alignment of the carpals with the

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

12

Proc IMechE Part H: J Engineering in Medicine

radius bones was present in both the experimental and modeling work, it was much more pronounced in the model, resulting in a greater medial–lateral shift (along the Fx-axis) of the carpals on the intra-articular surface and increased forces along this axis. While this modelexperimental disparity in the x-axis forces may have affected the calculation of the stresses and strains of the distal radius, the effect is predicted to be minimal. This is evident (described in detail below) by the similarity in the fracture patterns between the experimental specimens and those predicted by the model (Figure 6). During the experimental impact testing, the forces (along all axes) tended to decrease between the crack and fracture impact events, which may be attributed to a change in the material properties of the bone (i.e. decreased bone stiffness) as a result of the repeated impacts and the damage incurred during the crack event.20,53–55 However, in the current model, the changes in material properties were not simulated, thus the stiffness of the bone was maintained and the impact forces continued to increase; these differences were most notable along the Fz-axis. When the fracture forces from the current investigation are compared to those from previous finite element modeling studies, they tend to agree well. For example, Buchanan and Ural18 found a fracture load of approximately 4100 N for a simulation of a 50-year-old radius loaded at 18°, which represent conditions similar to those tested herein, aside from the underlying damage. This further suggests that the differences seen between the model and experimental forces during fracture are partially a result of the underlying changes to the bone seen experimentally that cannot, at this time, be incorporated into the modeling process. While the model and experimental acceleration data agree relatively well, the variation between the model and experimental accelerations is most likely a result of where the accelerations were calculated. The model accelerations were calculated at a single node directly on the surface of the bone, while the experimental accelerations were calculated from external sensors, whose centers were located approximately 0.5 cm away from the bone surface. These differences are analogous to the differences that exist between bone-mounted and skin-mounted accelerometers.56 Furthermore, although care was taken to clear the bone of all surrounding tissues, the experimental accelerations could be affected by underlying residual tissue. Differences between the experimental and model strains might also be explained by differences in the locations of where the strains were measured. While the numerical strains were measured from an element that corresponds to a location within the surface layer of the bone, the experimental strains by contrast were offset from the bone by the thickness of the strain gauge and strain gauge glue (;0.5 mm). Furthermore, the large differences between the model and experimental strains seen at the site of gauge 3 (proximal radius) are most likely a result of the complex interaction that occurs between the bone and cement just

proximal to this location, which was not considered in the model (i.e. this was modeled as a single rigid part). A number of different methods have been proposed to simulate bone fracture in FEMs. However, given the difficulty in simulating actual bone fracture and the subjectivity of the majority of fracture prediction methods, no generally accepted validated theory is currently available.57 The Coulomb–Mohr criteria,15,57 cohesive elements with a fracture plane,58 von-Mises stresses and maximum principal strains4 have been used previously. Quenneville and Dunning4 found that the vonMises stress criteria accurately predicted the location of fracture, and as such, this method was used in the current study. Despite the differences that were found between the model and experimental results, the model accurately predicted the location and magnitude of the damage inflicted by the impacts. The elements in the cortical bone that exceeded the critical von-Mises stresses, indicating fracture, were primarily located on the intra-articular surface and into the sigmoid fossa. The failed elements in the cancellous bone exhibited the same visual pattern as those in the cortical bone, but over a larger area. While there were no failed cortical elements in the space between the volar, ulnar aspect of the intra-articular surface and the sigmoid fossa, the volume of failed elements in these regions (;200– 400 mm3) suggests that propagation of the failure would have occurred.15

Conclusion This study presents a FEM of the distal radius and its response to dynamic impact loading. Three mesh quality assessment methods, and four validation techniques were used, and when reviewed collectively, they suggest that this model is a good representation of the distal radius. Furthermore, this model was found to simulate the pre-fracture and crack impact response well,5 and most importantly, the model was successful in predicting the location and severity of fractures that were obtained experimentally across all three impact events. Acknowledgement We thank Dr Bill Altenhof for providing much appreciated insight into this project. Declaration of conflicting interests The authors declare that there is no conflict of interest. Funding This study was funded by NSERC.

References 1. Rogge RD, Adams BD and Goel VK. An analysis of bone stresses and fixation stability using a finite element

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

Burkhart et al.

2.

3.

4.

5.

6.

7.

8.

9.

10.

11.

12.

13.

14.

15.

16.

13

model of simulated distal radius fractures. J Hand Surg Am 2002; 27: 86–92. Morgan EF and Bouxsein MK. Use of finite element analysis to assess bone strength. BoneKEy Osteovision 2005; 2: 8–19. Cristolfolini L, Schileo E, Juszczyk M, et al. Mechanical testing of bones: the positive synergy of finite-element models and in vitro experiments. Philos Trans A Math Phys Eng Sci 2010; 368: 2725–2763. Quenneville CE and Dunning CE. Development of a finite element model of the tibia for short-duration highforce axial impact loading. Comput Methods Biomech Biomed Engin 2011; 14: 205–212. Burkhart TA, Andrews DM and Dunning CE. Finite element modeling mesh quality, energy balance, and validation methods: a review with recommendations associated with the modeling of bone tissue. J Biomech 2013; 46: 1477–1488. Valle L and Ray MH. Development and validation of a 50th percentile male human femur: attachment A. Worchester, MA: National Highway Traffic Safety Administration, Worchester Polytechnique Institute, 2005. Tautges TJ. The generation of hexahedral meshes for assembly geometry: survey and progress. Int J Numer Meth Eng 2001; 50: 2617–2642. Fyllingen O, Hopperstad OS, Hanssen AG, et al. Brick elements versus shell elements in simulations of aluminum extrusions subjected to axial crushing. In: 7th European LS-DYNA conference, Salzburg, 14–15 May 2009, http://www.dynalook.com/european-conf-2009/GI-01.pdf Ulrich D, Van Rietbergen B, Laib A, et al. Load transfer analysis of the distal radius from in-vivo high-resolution CT-imaging. J Biomech 1999; 32: 821–828. Pistoia W, Van Rietbergen B, Lochmuller E-M, et al. Estimation of distal radius failure load with micro-finite element analysis models based on the three-dimensional peripheral quantitative computed tomography images. Bone 2002; 30: 842–848. Pistoia W, Van Rietbergen B and Ruegsegger P. Mechanical consequences of different scenarios for simulated bone atrophy and recovery in the distal radius. Bone 2003; 33: 937–945. Pistoia W, Van Rietbergen B, Lochmuller E-M, et al. Image-based micro-finite-element modeling for improved distal radius strength diagnosis: moving from ‘‘bench’’ to ‘‘bedside.’’ J Clin Densitom 2004; 7: 153–160. Carrigan SD, Whiteside RA, Pichora DR, et al. Development of a three-dimensional finite element model for carpal load transmission in a static neutral posture. Ann Biomed Eng 2003; 31: 718–725. Anderson DD, Deshpande BR, Daniel TE, et al. A threedimensional finite element model of the radiocarpal joint: distal radius fracture step-off and stress transfer. Iowa Orthop J 2005; 25: 108–117. Troy KL and Grabiner MD. Off-axis loads cause failure of the distal radius at lower magnitudes than axial loads: a finite element analysis. J Biomech 2007; 40: 1670–1675. Boutroy S, Van Rietbergen B, Sornay-Rendu E, et al. Finite element analysis based on in vivo HR-pQCT images of the distal radius is associated with wrist fracture in postmenopausal women. J Bone Miner Res 2008; 23: 392–399.

17. MacNeil JA and Boyd SK. Bone strength at the distal radius can be estimated from high resolution peripheral quantitative computed tomography and the finite element method. Bone 2008; 42: 1203–1213. 18. Buchanan D and Ural A. Finite element modeling of the influence of hand position and bone properties on the Colles’ fracture load during a fall. J Biomech Eng 2010; 132: 081007-1–081007-8. 19. Edwards WB and Troy KL. Finite element predictions of surface strain and failure load at the distal radius using simplified boundary conditions. Med Eng Phys 2012; 34: 290–298. 20. Burkhart TA, Andrews DM and Dunning CE. Failure characteristics of the intact distal radius in response to impact loading. J Orthop Res 2012; 30: 885–892. 21. Burkhart TA, Dunning CE and Andrews DM. Determining the optimal system specific cut-off frequencies for filtering in-vitro upper extremity impact force and acceleration data by residual analysis. J Biomech 2011; 44: 2728–2731. 22. Greenwald RM, Janes PC, Swanson SC, et al. Dynamic impact response of human cadaveric forearms using a wrist brace. Am J Sports Med 1998; 26: 825–830. 23. Staebler MP, Moore DC, Akelman E, et al. The effect of wrist guards on bone strain in the distal forearm. Am J Sports Med 1999; 27: 500–506. 24. Meyers ER, Hecker AT, Rooks DS, et al. Geometric variables from DXA of the radius predict forearm fracture load in vitro. Calcif Tissue Int 1993; 52: 199–204. 25. Meyers ER, Sebeny EA, Hecker AT, et al. Correlation between photon absorption properties and failure load of the distal radius in vitro. Calcif Tissue Int 1991; 49: 292–297. 26. Quenneville CE, Fraser GS and Dunning CE. Development of an apparatus to produce fractures from shortduration high-impulse loading with an application to the lower leg. J Biomech Eng 2012; 132: 014502-1–014502-4. 27. Tadepalli SC, Gandhi AA, Fredericks DC, et al. Cervical laminoplasty construct stability: an experimental and finite element investigation. Iowa Orthop J 2011; 31: 207–214. 28. Benzley SE, Perry E, Merkley K, et al. A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis. In: 4th international meshing roundtable, Albuquerque, NM, 16–17 October 1995, http://www.imr.sandia.gov/papers/Roundtable.agenda .html 29. Raut P. Impact of mesh quality parameters on elements such as beam, shell and 3D solid in structural analysis. Int J Eng Res App 2012; 2: 99–103. 30. Shim VB, Boshme J, Vaitl P, et al. An efficient and accurate prediction of the stability of percutaneous fixation of acetabular fractures with finite element simulation. J Biomech Eng 2011; 133: 094501-1–09501-4. 31. Schonenau E, Neu CM, Rauch F, et al. Gender-specific pubertal changes in volumetric cortical bone mineral density at the proximal radius. Bone 2002; 31: 110–113. 32. Burstein A, Reilly DT and Martens M. Aging of bone tissue: mechanical properties. J Bone Joint Surg Am 1976; 58: 82–86. 33. Imai K, Ohnishi I, Bessho M, et al. Nonlinear finite element model predicts vertebral bone strength and fracture site. Spine 2006; 31: 1789–1794. 34. Kim YS, Choi HH, Cho YN, et al. Numerical investigations of interactions between the knee-thigh-hip complex

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014

14

35.

36.

37. 38.

39. 40.

41.

42.

43.

44.

45.

Proc IMechE Part H: J Engineering in Medicine with vehicle interior structures. Stapp Car Crash J 2005; 49: 85–115. Peng L, Bai J, Zeng X, et al. Comparison of isotropic and orthotropic material property assignments on femoral finite element models under two loading conditions. Med Eng Phys 2006; 28: 227–233. Jones N. Strain rate sensitive behaviour of materials 1990. In: Jones N (ed.) Structural Impact. 1st ed. New York: Cambridge University Press, 1990, pp.333–384. McElhaney JH. Dynamic response of bone and muscle tissue. J Appl Physiol 1966; 21: 1231–1236. Quenneville CE. Experimental and numerical assessments of injury criteria for short-duration, high-force impact loading of the tibia. Doctoral Dissertation, Western University, London, ON, Canada, 2009. LSTC. LS-DYNA user’s manual, V.971. Livermore, CA: LSTC, 2007. Sambamoorthy B and Halder T. Characterization and component level correlation of energy absorbing (EA) polyurethane foams (PU) using LS-DYNA material models. In: 3rd European LS-DYNA user’s conference, Paris, 18–19 June, 2001. Zhang Y, Hughes TJR and Bajaj CL. Automatic 3D mesh generation for a domain with multiple materials. In: Proceedings of the 16th international meshing roundtable, Seattle, WA, 14–17 October 2007, http://www.imr .sandia.gov/papers/imr16.html Untaroiu C, Darvish K, Crandall J, et al. A finite element model of the lower limb for simulating pedestrian impacts. Stapp Car Crash J 2005; 49: 157–181. Ray MH, Mongiardini M, Atahan AO, et al. Recommended procedures for verification and validation of computer simulations used for roadside safety applications. Report number: 22–24, October 2008, Worchester MA: National Cooperative Highway Program, Worchester Polytechnique Institute. Schinkel-Ivy A, Altenhof W and Andrews DM. Validation of a full body finite element model (THUMS) for running-type impacts to the lower extremity. Comput Methods Biomech Biomed Engin 2014; 17: 137–148. Brewer JC. Effects of angles and offsets in crash simulation of automobiles with light trucks. SAE technical paper no. 2001-06-0077, 2001.

46. Cheng ZQ, Thaker JG, Pilkey WD, et al. Experiences in reverse-engineering of a finite element automobile crash model. Finite Elem Anal Des 2001; 37: 843–860. 47. Oberkampf WL and Trucano TG. Verification and validation in computational fluid dynamics. Prog Aerosp Sci 2002; 38: 209–272. 48. Jin SY, Majumder A, Altenhof W, et al. Axial cutting of AA6061-T6 circular extrusion under impact using singleand dual-cutter configurations. Int J Impact Eng 2010; 37: 735–753. 49. Bir C, Viano D and King A. Development of biomechanical response corridors of the thorax to blunt ballistic impacts. J Biomech 2004; 37: 73–79. 50. Craig M, Bir C, Viano D, et al. Biomechanical response of the human mandible to impacts of the chin. J Biomech 2008; 41: 2972–2980. 51. Gatti CJ, Marrat JD, Palmer ML, et al. Development and validation of a finite element model of the superior glenoid labrum. Ann Biomed Eng 2010; 38: 3766–3776. 52. Taddei F, Cristofolini L, Martelli S, et al. Subject-specific finite element models of long bones: an in vitro evaluation of the overall accuracy. J Biomech 39: 2457–2467. 53. Martin RB, Gibson VA, Stover SM, et al. Residual strength of equine bone is not reduced by intense fatigue loading: implications for stress fracture. J Biomech 1997; 30: 109–114. 54. Reilly GC and Currey JD. The effects of damage and micro-cracking on the impact strength of bone. J Biomech 2000; 33: 337–343. 55. Gupta HS and Zioupos P. Fracture of bone tissue: the ‘‘hows’’ and ‘‘whys.’’ Med Eng Phys 2008; 30: 1209–1226. 56. Lafortune MA, Henning E and Valiant GA. Tibial shock measured with bone and skin mounted transducers. J Biomech 1995; 28: 989–993. 57. Keyak JH and Rossi SA. Prediction of femoral fracture load using finite element models: an examination of stress- and strain-based failure theories. J Biomech 2000; 33: 209–214. 58. Ural A. Prediction of Colles’ fracture load in human radius using cohesive finite element modeling. J Biomech 2009; 42: 22–28.

Downloaded from pih.sagepub.com at Uppsala Universitetsbibliotek on October 3, 2014