"Applied Science, Faculty of"@en . "Mechanical Engineering, Department of"@en . "DSpace"@en . "UBCV"@en . "Ariza, Oscar Rodrigo"@en . "2014-04-07T14:16:33Z"@en . "2014"@en . "Master of Applied Science - MASc"@en . "University of British Columbia"@en . "There is currently a disagreement between the hips that today\u00E2\u0080\u0099s screening techniques identify as likely to fracture and those that actually fracture. Specimen-specific finite element (FE) analysis based on computed tomography (CT) data has been presented as a more sensitive technique than standard screening based on bone mineral density (BMD) to screen for and identify hips predisposed to fracture. However, published studies using this technique have applied loading scenarios that do not sufficiently resemble the falls that typically lead to fracture. In particular, the effects of dynamics and strain rate stiffening have been neglected in all relevant FE studies.\n\nIt is my objective to explore here these previously neglected aspects. This thesis is divided into three sections focusing respectively on: 1) the effect of material stiffness and mapping techniques on a simplified quasi-static model; 2) the feasibility, accuracy, and model sensitivity of explicit dynamic FE analysis to simulate a fall impact on the femur; and 3) the feasibility of a novel multi-scale FE technique to analyse femoral bone behaviour at specific sites with a level of detail that resolves the trabecular structure. All models are compared against experimental data obtained from mechanical tests of femurs at appropriate loading rates.\n\nResults demonstrate that both dynamic and multi-scale modelling are feasible techniques that can be developed into powerful tools. Both quasi-static and dynamic FE models demonstrated sensitivity to material modulus-density relationship, with further work required to estimate the increased stiffness of bone during impact loading. For dynamic models, surface strain patterns matched fracture locations observed from high-speed video. The multi-scale technique demonstrated its value by identifying a potential stress concentration around a vascular hole that was otherwise hidden.\n\nThis study demonstrates that explicit dynamic FE analysis is a feasible technique for modelling the altered bone behaviour observed during femoral impacts, yet further research is required to design material models appropriate for impact rates. It is also demonstrated that nested multi-scale FE modelling is not only feasible but also potentially useful to identify microstructural features of bone that might make it prone to fracture even in the setting of relatively high BMD."@en . "https://circle.library.ubc.ca/rest/handle/2429/46324?expand=metadata"@en . " A NOVEL APPROACH TO FINITE ELEMENT ANALYSIS OF HIP FRACTURES DUE TO SIDEWAYS FALLS by Oscar Rodrigo Ariza B.A.Sc. University of Waterloo, 2010 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate and Postdoctoral Studies (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2014 ? Oscar Rodrigo Ariza, 2014 ii Abstract There is currently a disagreement between the hips that today?s\t ?screening techniques identify as likely to fracture and those that actually fracture. Specimen-specific finite element (FE) analysis based on computed tomography (CT) data has been presented as a more sensitive technique than standard screening based on bone mineral density (BMD) to screen for and identify hips predisposed to fracture. However, published studies using this technique have applied loading scenarios that do not sufficiently resemble the falls that typically lead to fracture. In particular, the effects of dynamics and strain rate stiffening have been neglected in all relevant FE studies. It is my objective to explore here these previously neglected aspects. This thesis is divided into three sections focusing respectively on: 1) the effect of material stiffness and mapping techniques on a simplified quasi-static model; 2) the feasibility, accuracy, and model sensitivity of explicit dynamic FE analysis to simulate a fall impact on the femur; and 3) the feasibility of a novel multi-scale FE technique to analyse femoral bone behaviour at specific sites with a level of detail that resolves the trabecular structure. All models are compared against experimental data obtained from mechanical tests of femurs at appropriate loading rates. Results demonstrate that both dynamic and multi-scale modelling are feasible techniques that can be developed into powerful tools. Both quasi-static and dynamic FE models demonstrated sensitivity to material modulus-density relationship, with further work required to estimate the increased stiffness of bone during impact loading. For dynamic models, surface strain patterns matched fracture locations observed from high-speed video. The multi-scale technique demonstrated its value by identifying a potential stress concentration around a vascular hole that was otherwise hidden. This study demonstrates that explicit dynamic FE analysis is a feasible technique for modelling the altered bone behaviour observed during femoral impacts, yet further research is required to design material models appropriate for impact rates. It is also demonstrated that nested multi-scale FE modelling is not only feasible but also potentially useful to identify microstructural features of bone that might make it prone to fracture even in the setting of relatively high BMD. iii Preface The investigation presented in this thesis is part of a collaborative project between the University of British Columbia (UBC) ? through the Orthopaedic and Injury Biomechanics Group (OIBG) and the Centre for Hip Health and Mobility (CHHM) ? and the Swiss Federal Institute of Technology (ETH) in Zurich, Switzerland ? through the Institute for Biomechanics (IfB). The work presented in chapter 4 was included as part of an accepted publication: B. Helgason, S. Gilchrist, O. Ariza, J. D. Chak, G. Zheng, R. Widmer, S. J. Ferguson, P. Guy, and P. A. Cripton, ?Development\t ?of\t ?a\t ?balanced\t ?experimental-computational approach to understanding the mechanics\t ?of\t ?proximal\t ?femur\t ?fractures,?\t ?Medical\t ?Engineering\t ?&\t ?Physics,\t ?no.\t ?In\t ?Review,\t ?2014.\t ?I created and processed all the simulations, updating work previously done for an unpublished paper, and then I re-evaluated and expanded the analysis previously done. My responsibilities in this project involved work in both UBC and ETH. I wrote the semi-automated code required to build all finite element models, as well as post-processing code used to create figures and data analysis. All writing in this thesis is my own, although a few paragraphs\t ?in\t ?the\t ?introduction\t ?sections\t ?have\t ?been\t ?based\t ?on\t ?Dr.\t ?Helgason?s\t ?writing\t ?for\t ?related\t ?publications. Unless otherwise stated in the caption, every figure in this thesis is my own creation. However, for figures describing the mechanical testing, I may have only added labels over photographs or video captures created by Seth Gilchrist, the PhD candidate at UBC whose experimental work forms the basis of the research here presented. Dr. Benedikt Helgason co-supervised my work throughout this project. He also created early versions of the FE modelling code. For chapter 2, he developed a cartilage model that improved my base modelling technique. Dr. Helgason guided me through all the decisions and analysis necessary for this research. Seth Gilchrist designed, executed, post-processed, and analysed all the experimental work used for validation, including all mechanical tests and QCT scans presented in chapters 2 and 3. For chapter 4, Seth post-processed and analysed experimental data collected by Jason Chak, a research engineer at OIBG, and also developed and executed the digital image correlation technique used for strain comparison. He also helped analyse all dynamic results. iv Dr. Peter Cripton co-supervised my work throughout this project and guided the focus and analysis directions of my work. He also provided much feedback on how to appropriately present\t ?the\t ?project?s\t ?results\t ?in\t ?this\t ?thesis. v Table ?of ?contents Abstract ............................................................................................................................................................... ii Preface ................................................................................................................................................................ iii Table of contents ............................................................................................................................................. v List of tables ..................................................................................................................................................... vii List of figures ................................................................................................................................................. viii List of abbreviations ....................................................................................................................................... x Acknowledgements ........................................................................................................................................ xi Dedication ........................................................................................................................................................ xii 1 Introduction .............................................................................................................................................. 1 1.1 Motivation ......................................................................................................................................................... 1 1.2 Anatomy of the proximal femur ............................................................................................................... 2 1.3 Mechanics of hip fracture ........................................................................................................................... 4 1.4 Mechanical testing of femurs .................................................................................................................... 7 1.5 Explicit and implicit finite element methods...................................................................................... 7 1.6 Finite element analysis of hip fracture ............................................................................................... 10 1.7 Multi-scale nested FE modelling ............................................................................................................ 11 1.8 Objectives ........................................................................................................................................................ 13 1.9 Scope ................................................................................................................................................................. 14 2 Quasi-static linear simulation ..........................................................................................................15 2.1 Introduction ................................................................................................................................................... 15 2.2 Methodology .................................................................................................................................................. 17 2.2.1 Experimental setup ........................................................................................................................... 18 2.2.2 FE model geometry ............................................................................................................................ 20 2.2.3 Heterogeneous material stiffness mapping ............................................................................. 21 2.2.4 Boundary conditions ......................................................................................................................... 23 2.2.5 Comparison criteria........................................................................................................................... 25 2.3 Results .............................................................................................................................................................. 26 2.4 Discussion ....................................................................................................................................................... 28 3 Dynamic nonlinear simulation .........................................................................................................33 3.1 Introduction ................................................................................................................................................... 33 3.2 Methodology .................................................................................................................................................. 34 3.2.1 Experimental setup ........................................................................................................................... 35 3.2.2 Geometry and material mapping ................................................................................................. 37 3.2.3 Material properties ............................................................................................................................ 38 3.2.4 Boundary conditions ......................................................................................................................... 41 vi 3.2.5 Comparison criteria........................................................................................................................... 44 3.2.6 Parameters tested for sensitivity ................................................................................................. 46 3.3 Results .............................................................................................................................................................. 47 3.3.1 Ex-vivo fall simulations .................................................................................................................... 47 3.3.2 Force, displacement, and time relationships .......................................................................... 49 3.3.3 Stiffness, peak force, and energy regressions ......................................................................... 54 3.3.4 High strains and fracture locations ............................................................................................. 56 3.3.5 Sensitivity analyses ........................................................................................................................... 61 3.4 Discussion ....................................................................................................................................................... 64 4 Multi-scale nested simulation ..........................................................................................................68 4.1 Introduction ................................................................................................................................................... 68 4.2 Methodology .................................................................................................................................................. 69 4.2.1 Ex-vivo and organ-level models ................................................................................................... 69 4.2.2 Micro geometry and material properties ................................................................................. 70 4.2.3 Boundary conditions ......................................................................................................................... 71 4.2.4 Comparison criteria........................................................................................................................... 72 4.3 Results .............................................................................................................................................................. 72 4.3.1 Ex-vivo fall simulation ...................................................................................................................... 72 4.3.2 Organ-level model .............................................................................................................................. 73 4.3.3 Micro-level models ............................................................................................................................ 73 4.4 Discussion ....................................................................................................................................................... 75 5 Conclusions .............................................................................................................................................77 References ........................................................................................................................................................79 Appendix A ? Energy balance during a drop tower test ..................................................................87 Appendix B ? Maximum principal strains from dynamic model ..................................................89 vii List ?of ?tables Table 2-1: General data for femoral specimens .................................................................................................. 18 Table 2-2: Mesh characteristics for each bone .................................................................................................... 20 Table 2-3: Selected modulus-density relationships .......................................................................................... 21 Table 3-1: General data for femoral specimens .................................................................................................. 35 Table 3-2: Mesh characteristics for each bone .................................................................................................... 38 Table 3-3: Comparison of initial experimental fracture locations ? as marked by Gilchrist ? against compressive principal FE strain patterns ................................................................................. 57 viii List ?of ?figures Figure 1-1: Proximal femur location and anatomy from anterior (middle) and posterior (right) views. Overview image on the left adapted from BodyParts3D, ? The Database Center for Life Science, licensed under CC Attribution-Share Alike 2.1 Japan .................................................. 3 Figure 1-2: CT-derived coronal section of proximal femur showing the internal structure of a bone with normal density (left) and of one classified as osteoporotic (right) ............................ 4 Figure 1-3: Relative incidence and approximate location of the different types of hip fracture. Adapted from local Baltimore data collected by Michelson [19], which is consistent with estimates found by other authors [18]. ....................................................................................................... 5 Figure 1-4: Stress reversal on the femoral neck occurring during sideways falls compared to regular upright loading. Adapted from Turner [27] with permission. ........................................... 6 Figure 2-1: Sample femur placed in positioning fixture ................................................................................. 19 Figure 2-2: Comparison of modulus-density relationships over bone density range ........................ 22 Figure 2-3: Example of digitized points used for registration (The registered CT-derived mesh and calculated geometry are shown for clarity) .................................................................................... 23 Figure 2-4: Sample FE model and boundary conditions ................................................................................. 25 Figure 2-5:\t ?Definition\t ?of\t ?femoral\t ?stiffness.\t ?F\t ?represents\t ?the\t ?applied\t ?force\t ?and\t ??\t ?the\t ?displacement\t ?of the greater trochanter ................................................................................................................................. 26 Figure 2-6: Comparison of experimental to FE stiffness results for each modulus-density relationship and material mapping method ............................................................................................ 27 Figure 2-7: Comparison of experimental to FE strain results for each modulus-density relationship and material mapping method ............................................................................................ 28 Figure 3-1: Overview (left) and detail (right) of fall simulator and test setup. The white straps shown in the overview picture, meant to prevent prolonged foam compression, were removed before testing. ................................................................................................................................... 36 Figure 3-2: Overview of nonlinear material properties .................................................................................. 40 Figure 3-3: FE boundary conditions ........................................................................................................................ 42 Figure 3-4: Sample frames from displacement tracking video showing the superior surface of the neck. The arrows mark the approximate tracked position of the impactor block (red, pointing left) and PMMA-femur interface (green, pointing right) ................................................. 44 Figure 3-5: Energy comparison from FE data using time and displacement matching. Experimental results are in black, FE results in blue ........................................................................... 46 Figure 3-6: Comparison between stiff, mean, and compliant foam behaviour ...................................... 47 Figure 3-7: Samples of three behaviours observed experimentally. The dashed lines represent the time at peak force ........................................................................................................................................ 48 Figure 3-8: Force-time comparison between experimental and FE results ........................................... 51 Figure 3-9: Force-displacement comparison between experimental and FE results ......................... 52 Figure 3-10: Displacement-time comparison between experimental and FE results ........................ 53 Figure 3-11: Stiffness results regression ............................................................................................................... 54 ix Figure 3-12: Peak force results regression based on time-matching (left) and displacement matching (right) .................................................................................................................................................. 55 Figure 3-13: Energy of deformation regression based on time-matching (left) and displacement matching (right) .................................................................................................................................................. 56 Figure 3-14: Comparison of sensitivity to foam properties .......................................................................... 62 Figure 3-15: Comparison of sensitivity to nonlinear material behaviour ............................................... 63 Figure 3-16: Comparison of results using two modulus-density relationships .................................... 64 Figure 4-1: Location of VOIs on organ-level model ........................................................................................... 71 Figure 4-2: Comparison of organ-level surface strains to DIC-based strains calculated by Gilchrist .................................................................................................................................................................................... 73 Figure 4-3: Comparison of minimum principal strains and deformations resulting from micro-level FE models to fracture locations from high-speed video .......................................................... 75 x List ?of ?abbreviations aBMD Areal bone mineral density BMD Bone mineral density BV/TV Bone volume over total volume (i.e. bone volume fraction) CT Computed tomography DIC Digital image correlation DXA Dual energy X-ray absorptiometry ETH Swiss Federal Institute of Technology (German: Eidgen?ssische Technische Hochschule) FE Finite element FEA Finite element analysis FH Femoral head GT Greater trochanter HA Hydroxyapatite NURBS Non-uniform rational b-spline PMMA Poly(methyl-methacrylate) PVC Polyvinyl chloride QCT Quantitative computed tomography RMSE Root mean squared error VOI Volume of interest xi Acknowledgements I would like to express my gratitude to my supervisor Dr. Peter Cripton at UBC for his input and guidance throughout the project and especially for his comments on this manuscript. I also extend thanks to Dr. Benedikt Helgason at ETH for his guidance and encouragement on each step of this investigation. Similarly, I am grateful to Seth Gilchrist for his monumental effort developing the experimental work on which I based the simulations here presented. This thesis would not exist without the contributions from all three researchers. I gratefully acknowledge the help of my lab mates at both UBC and ETH. I owe special thanks to Hannah Gustafson, Ren? Widmer, and Giacomo Marini for their assistance with various portions of my research. I am also grateful to Dr. Stephen Ferguson for accommodating me in his lab for a full year as an exchange student. I would like to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for its financial support during the first year of research. I thank my mom, my siblings, and my friends for their support and encouragement. I am especially grateful to my friend Shuo Yang for helping me stay motivated while writing this thesis. Finally, I must thank my girlfriend Ching Chi for her patience, love, and support throughout these years. xii Dedication To my parents, who taught me to be patient and resilient. 1 1 Introduction 1.1 Motivation Every year, there are more than 1.6 million hip fractures worldwide [1], more than 23 000 of them occurring in Canada [2]. This type of fracture not only has a very negative effect on a person?s\t ?mobility\t ?and\t ?quality\t ?of\t ?life,\t ?but\t ?it\t ?is\t ?also\t ?related\t ?to\t ?high\t ?levels\t ?of\t ?mortality, with 20% of people dying within one year after suffering a hip fracture [3]. The financial burden of care associated with hip fractures has also been estimated around $650 million per year in Canada [4]. All of these numbers are expected to increase over time, as baby boomers are reaching old age and worldwide life expectancy is rising, leading to higher proportions of elderly people prone to hip fracture [5]. Much of the current biomechanical research on the femur is geared toward improved fixation or bone replacement after a femoral fracture has occurred. While this research is necessary, it is always preferable to prevent a fracture in the first place. While there are some existing preventive\t ?methods,\t ?an\t ?accurate\t ?screening\t ?procedure\t ?is\t ?necessary\t ?to\t ?determine\t ?a\t ?person?s\t ?relative risk of fracture in order to efficiently administer the appropriate treatment. Currently, a person?s\t ?risk\t ?is\t ?most often estimated according to their bone mineral density (BMD), as measured through dual-energy X-ray absorptiometry (DXA) [6]. While it has been shown that low BMD is associated with elevated hip fracture risk [7], studies have also found that there is a large overlap between the BMD of people who suffer a fracture and the BMD of people with matching age and gender who do not [8]. While it has been shown that BMD is a good predictor of\t ?hip\t ?fracture\t ?risk\t ?at\t ?a\t ?population\t ?level,\t ?an\t ?individual?s\t ?BMD\t ?is\t ?not sufficient to determine his or her risk. The results of these previous studies hint at the existence of factors beyond BMD that significantly contribute to an\t ?individual?s\t ?hip fracture risk [9]. It has been suggested that architecture, i.e. the general orientation and position of load-bearing trabecular paths, plays a significant role in the ability of a femur to avoid fracture in a fall [10]. Given the variability of femoral architecture among the population, patient-specific finite element (FE) models have been proposed to provide more sensitive predictions of hip fracture locations and the loads at which fracture initiation occurs [11]?[14]. These models are generally developed from quantitative computed tomography (QCT) scans, allowing heterogeneous representation of the internal geometry and density of the bone. This technique has the advantage of allowing parametric representation of the material and geometric complexity of 1 ? Introduction 2 femurs ? a complexity that is difficult to model using other techniques. Similarly, unlike experimental techniques, the parametric representation gives researchers and clinicians the ability to subject a particular bone to fracture under multiple loading conditions or property distributions [11]. 1.2 Anatomy of the proximal femur The femur is the longest, largest, and strongest bone in the human body [15]. It spans from the hip joint, where it articulates with the pelvis, to the knee joint, where it articulates with the tibia and patella. The femur broadly consists of proximal and distal extremities joined by an almost cylindrical shaft. For hip fracture research purposes, we are mainly interested in the proximal femur, encompassing the proximal (superior) extremity and the top portion of the shaft. The most proximal portion of the femur is the femoral head (Figure 1-1), a spheroidal section of bone covered with cartilage that slides against the acetabulum in the pelvis. The femoral neck is a longitudinally concave cylinder that connects the femoral head to the shaft [16]. The greater trochanter (GT) is a large protuberance that forms the most lateral and palpable portion of the femur. It serves as the point of attachment for several muscles that connect the femur to the pelvic bone. The lesser trochanter is a smaller feature protruding from the posteromedial shaft, just inferior to the GT. It provides attachment to large muscles connected to the iliac crest of the pelvic bone and to the transverse processes of the lumbar vertebrae. The two trochanters are connected anteriorly by the intertrochanteric line and posteriorly by the intertrochanteric crest. The circle created by the two trochanters and two intertrochanteric ridges provide attachment for a number of ligaments that connect to several locations in the pelvic bone, including those that form the hip joint capsule. 1 ? Introduction 3 Internally, the proximal femur is composed mainly of trabecular bone surrounded by a thin layer of cortical bone. Cortical bone thickens inferiorly on the proximal end of the shaft, while trabecular bone becomes sparser in the same region. Gradually, the internal structure of the shaft turns into a thick cylinder of cortical bone mainly filled with marrow in the medullary cavity. The trabeculae at the head are arranged in such a way that they form a dense central wedge supported mainly against the upper and lower profiles of the neck [15]. The typical vertical loads applied when bearing the body weight are thus passed from the femoral head through the central wedge and toward the junction of the neck and shaft. On patients with low bone density, both the thickness of the cortical shell and the density of the internal trabecular structure are affected [17] (Figure 1-2). Figure 1-1: Proximal femur location and anatomy from anterior (middle) and posterior (right) views. Overview image on the left adapted from BodyParts3D, ? The Database Center for Life Science, licensed under CC Attribution-Share Alike 2.1 Japan 1 ? Introduction 4 1.3 Mechanics of hip fracture Hip fractures are fractures of the femur around the hip joint area. They are generally classified according to their relative position on the femur (Figure 1-3). When the fracture occurs in the femoral neck region, between the trochanters and the femoral head, it is referred to as a femoral neck fracture [18]. It can also be referred to as an intracapsular fracture, as this area is inside of the articular capsule of the hip joint. If the fracture is located more distally, in the area between the greater and lesser trochanters, it is referred to as an intertrochanteric fracture. A fracture located even farther distally is referred to as a subtrochanteric fracture. About 37% of hip fractures are in the femoral neck, 49% are intertrochanteric, and about 14% are subtrochanteric [19]. Femoral neck and intertrochanteric fractures are often the result of falls from standing height on the greater trochanter, particularly for elderly patients. Subtrochanteric fractures, on the other hand, are typically the result of high energy impacts such as motor vehicle accidents and falls from a height [20] Figure 1-2: CT-derived coronal section of proximal femur showing the internal structure of a bone with normal density (left) and of one classified as osteoporotic (right) 1 ? Introduction 5 The three different types of fractures have different challenges and treatments associated with them. Femoral neck fractures are likely to disrupt the blood supply to the femoral head, possibly leading to necrosis of this area. For this reason, femoral neck fractures often require hip replacement surgery, in which the proximal femur ? sometimes together with the acetabulum ? is removed and replaced with an artificial implant. Intertrochanteric fractures are less likely to disrupt blood supply, but the position of tendons and ligaments around the fractures makes the broken pieces more likely to become displaced and to have healing complications. Subtrochanteric fractures tend to heal more slowly than intertrochanteric fractures because the blood supply to the subtrochanteric region is not as good as that in the intertrochanteric region [18]. While some hip fractures occur spontaneously because of a sudden or excessive muscle contraction [21], [22], clinical studies have suggested that most hip fractures are the result of sideways falls from standing [23]. It is important to note, however, that only around 5% of falls Figure 1-3: Relative incidence and approximate location of the different types of hip fracture. Adapted from local Baltimore data collected by Michelson [19], which is consistent with estimates found by other authors [18]. 1 ? Introduction 6 in the elderly result in fracture [24]. It has been shown that bone loss in the elderly femur occurs preferentially at the superior aspect of the neck [25], and this same region has been observed to be a site of neck fracture initiation during ex-vivo experiments under sideways fall loading [26]. This fragility has been commonly explained by noting that the superior femoral neck only experiences light tensile stresses during\t ?regular\t ?standing\t ?and\t ?walking;\t ?thus,\t ?following\t ?Wolff?s\t ?Law, bone in this region tends to weaken over time if alternative loading is not regularly experienced. During a sideways fall impact to the greater trochanter, however, the same area is heavily loaded in compression (Figure 1-4). The weakened structure is often unable to support the sudden increased load, leading to higher likelihood of fracture [27]. The specific mechanism that leads to fracture initiation in the superior femoral neck under compression is a subject of debate. Mayhew [28] suggested that local buckling of the thinning cortical shell may play an important role during fractures, given that this area of the neck often becomes thinner with age. Turner [27], however, countered the argument by noting that shell buckling is unlikely given the reinforcing effects of trabecular bone and bone marrow. While fractures have been observed to initiate in areas that are subjected to compression during sideways falls and where thinning cortical bone is common [26], a finding that potentially supports the buckling hypothesis, to the best of our knowledge a buckling fracture has not been directly observed in any published study to date. Figure 1-4: Stress reversal on the femoral neck occurring during sideways falls compared to regular upright loading. Adapted from Turner [27] with permission. 1 ? Introduction 7 1.4 Mechanical testing of femurs A number of studies have performed mechanical testing of human femurs to analyse its behaviour under load configurations representative of both standing and falling [21], [26], [29]?[34]. Studies have used mechanical testing techniques to establish the importance of loading direction on fracture load during a fall [35], [36]. An increase in femoral strength and stiffness with increasing loading rate has been similarly demonstrated [30], [37]. Researchers have almost exclusively used materials testing machines to apply the loads on the femur specimens. These machines allow strict control over the displacement rate or forces throughout a test, a feature that is of great use for determining the material properties of a viscoelastic material like bone. However, materials testing machines also have limitations when experimentally simulating falls. First, many hydraulic machines ? which are the most common type ? can generally achieve maximum displacement rates around 0.1 m/s [37], far from the estimated impact speed of 3.0 m/s during a sideways fall [38], [39]. Second, when testing to fracture, the machine has access to an infinite amount of energy to transfer to the bone; this scenario is not equivalent to a sideways fall, where energy is limited by the fall height and fracture is not prescribed ? to reiterate, only about 5% of falls in the elderly result in fracture [24]. Finally, during a sideways fall, the displacement rate is not constant but is dependent ? among other things ? on the stiffness of the femur itself. Since femoral stiffness is also dependent on loading rate [30], [37], a recursive process is present, which makes it virtually impossible to perform a biofidelic simulation of a fall using a materials testing machine. Alternative methods to simulate sideways falls impacting the femur can be found in Robinovitch?s\t ?studies\t ?on\t ?hip\t ?protectors\t ?using an inverted pendulum [40] and in\t ?our\t ?lab?s\t ?own\t ?hip fracture research using a drop tower setup [41]. These experimental, inertia-driven devices overcome all the shortcomings of materials testing machines listed above at the expense of losing precise control over the loading conditions after the start of a test. The experimental data used throughout this thesis ? collected and described by Gilchrist in the above citation ? comes from experiments performed using the drop tower setup. 1.5 Explicit and implicit finite element methods The finite element method is a numerical technique used to obtain approximate solutions to differential equations. This technique, which was developed gradually by mathematicians and engineers over a period of more than a hundred years [42], is used in various fields of science and engineering to model physical problems. The process starts by dividing a large and complex 1 ? Introduction 8 system into collections of a finite, though often large, number of simple elements. The simple elements are assembled together systematically and according to the equations that govern the physical problem before boundary conditions are applied. An approximate answer to the problem is then obtained by solving ? almost always with the aid of powerful computers ? the resulting set of system equations [43]. In solid mechanics, the problem generally involves either finding the displacement of a complex structure under a given load or finding the load generated when a given displacement is applied. In either case, the structure is first divided, or meshed, into elements that can be assumed to have a simple mechanical response ? e.g. beams, struts, cables, or small solid pieces. The points of connection between multiple elements and between the structure and its surroundings ? i.e. the nodes ? are also noted at this stage. Material models that describe the constitutive relationship between internal forces and displacements on each element are assigned. The elements are then assembled, usually in an automated fashion, using equations that ensure that the physical requirements of displacement compatibility (no two portions of the structure occupy the same space) and equilibrium (internal forces add up to zero at all points) are maintained. The result here is a set of equations describing the relationships between nodal displacements and any external forces on the structure. The final step is applying boundary conditions, i.e. assigning external forces and displacement constraints. The system can then be solved to obtain the desired forces or displacements. If nonlinearities are present in the system because of complex materials, large rotations or deformations, or contact between detached elements, then the solution must be computed in steps. Since this subject is relevant for this thesis, the rest of this section describes the difference between the two main stepwise methods used to solve nonlinear systems. FE modelling of solids involves finding the solution to the equation, Ku = F, where F represents the external forces on the solid, u represents displacements, and K is the stiffness matrix of the system. On linear problems the stiffness of the system is constant, allowing the solution to be found in a single step. On nonlinear and dynamic problems, on the other hand, K and F can be functions of u and its time derivatives ??? and ???. These problems require the solution of the stiffness equation in a series of time steps. The stiffness matrix, K, must be updated on each time step to account for deformed geometries, material changes, contact changes, or any other type of nonlinearity present in the system. There are two main methods used to perform this stepwise approximation; they are referred to as implicit and explicit FE analyses [44]. 1 ? Introduction 9 Explicit methods of analysis obtain the solution for each time step by using information from the previous time step. In other words, the solution for time t + ?t is only dependent on the solution at time t. The explicit methodology allows straightforward computation of the stiffness matrix at each step, but it does not strictly enforce equilibrium, as internal forces are balanced using the displacements at the previous time step and not at the current one. It is imperative that the length of time steps is kept short to keep this force imbalance small. Explicit methods are said to be conditionally stable; that is, the error introduced on each step is only bounded as long as the length of the time step, ?t, remains below a given value. For solid dynamic analysis, this maximum value of ?t\t ?is 2/?, where ?\t ?is the natural frequency of the element with the shortest side length [45]. Implicit methods are the most commonly used type of nonlinear FE analysis. They differ from the explicit methodology in that equilibrium is strictly enforced on each time step. This is done by calculating the difference between internal and external forces on the current time step, generally called the residual, and solving the equation that reduces this residual to a specified tolerance [46]. With this procedure, however, the solution for time t + ?t becomes dependent not only on the solution at time t but also on the solution at time t + ?t (i.e. itself). For this reason, the governing equation must generally be obtained iteratively, most commonly through the use of the Newton-Raphson method. The advantage of implicit methods is that they are unconditionally stable, allowing the use of much longer time steps ? and thus fewer steps in total ? than explicit methods. The disadvantage is that the complexity of solving each time step is greatly increased, leading to longer computation time for each step and the possibility of non-convergence on large deformations. The choice between explicit and implicit methods depends on the type of simulation being performed. Implicit methods are better suited for the analysis of systems with static or slowly increasing loads, as using explicit methods on them tends to require an unreasonably large number of short time steps to maintain stability. In this case, the ability of implicit methods to use large steps easily overcomes the increased complexity of each step. However, in many cases where an FE model simulates an impact event with relatively compliant elements, explicit methods may be more appropriate. The short duration of an impact event keeps the overall number of steps small. Additionally, since natural frequency is proportional to ?? ?? ? where k is stiffness and m is mass, high compliance (i.e. low stiffness) leads to low natural frequencies and thus large time steps allowed for explicit stability. Furthermore, explicit methods more 1 ? Introduction 10 easily overcome the convergence problems associated with large sudden deformations. Explicit methods have been successfully used on car, train, and airplane collisions, as well as metal stamping and extrusion [45], and even human joint mechanics [47]. All of these applications showcase the suitability of explicit analysis on the computation of transient dynamic loads. 1.6 Finite element analysis of hip fracture The finite element method has been used on the analysis of the femur since at least the 1970s [48], focusing over time on the analysis of prosthetic implants [49]. The use of FE modelling to specifically investigate hip fracture risk of the intact proximal femur started in the early 1990s with the investigation by Lotz, Cheal, and Hayes. They used data from quantitative computed tomography (QCT) images to create models of the proximal femur and load them in both single-leg stance and sideways fall configurations using both a linear-elastic model [11] and a second model with asymmetric nonlinear material properties applied on the most highly strained regions [50]. While their results demonstrated differences of up to 56% in force and up to 550% in surface strains compared to experimental measurements, the feasibility of the method was demonstrated, providing a platform for further research. Over the past 20 years, a number of studies have created specimen-specific, QCT-based FE models in single-leg stance [12]?[14], [51]?[57] and sideways fall [12], [52], [54], [58]?[60] loading configurations. These studies have generally used heterogeneous material properties on increasingly denser femur meshes, validating the models against ex-vivo mechanical tests. All published studies we have reviewed have used static techniques and assumptions to virtually load the bones, matching the relatively low loading rates used experimentally for model validation. This is despite the fact that sideways falls are viewed and analysed as dynamic events outside of FE-based research [39], [61] and that dynamic effects have been identified as a significant limitation of existing FE models [62]. Perhaps this lack of dynamic simulations is related to the difficulties of modelling dynamic effects with traditional implicit methods. To the best of our knowledge, there are currently no published studies that use explicit FE for hip fracture research, with the exception of our own work currently in review [63]. The main variation found among existing FE studies of femoral fracture relates to the materials used in the simulations. Most studies have used linear isotropic material models [11], [12], [14], [51], [53], [55], [57]?[59], [64], and some of them have found very good correlations with empirical measurements using carefully derived properties [55], [57], [59]. Several researchers have obtained improved results against in-vitro experiments while using nonlinear properties in 1 ? Introduction 11 the form of elastoplastic models [13], [62], or with the inclusion of strain softening or damage [52], [65]. While bone is known to be anisotropic in nature, almost all FE studies have used simplified isotropic material properties. It is unclear to what degree the addition of anisotropy would affect FE results, as published investigations of the relative effect of orthotropic vs. isotropic material properties on FE models of the femur have produced conflicting results [66]?[68]. In an attempt to obtain the material properties that best work with FE models of the femur, some published studies have adjusted material property parameters to fit experimental results. This has\t ?usually\t ?been\t ?done\t ?by\t ?dividing\t ?the\t ?specimens\t ?tested\t ?mechanically\t ?into\t ??training?\t ?and\t ??testing?\t ?groups\t ?[51], [58]. While this technique is helpful to find an initial range for material parameters, it carries the risk of narrowing the applicability of developed material models, as any change in experimental conditions necessitates a new training group and new parameters. FE models must ultimately be developed using material data collected independently and systematically from bone cores to provide a proper base for validation against mechanical experiments. Furthermore, the process from CT data to FE model results should ideally be as automated as possible [69]; this is not only to provide a basis for future clinical implementation of FE models but also to ensure a fair comparison by making the specimen the only variable. 1.7 Multi-scale nested FE modelling As described in section 1.2, the internal proximal femur consists mainly of a trabecular bone lattice-like structure with marrow filling the intertrabecular pores. Structurally, the trabeculae act as beams and struts to bear load and provide bracing for the load-carrying cortical bone. The density and orientations of the trabeculae are not random but instead specifically arranged to carry the regular local loading applied to them ? a\t ?relationship\t ?commonly\t ?known\t ?as\t ?Wolff?s\t ?Law.\t ?Thus, when taken as a whole, trabecular bone has a wide range of possible mechanical properties (e.g. stiffness or yield strength) and varying degrees of anisotropy with variable principal directions. The thickness of cortical bone also varies greatly with location; compare, for example, the difference in cortical thickness between the inferior-medial neck and the superior-lateral head in Figure 1-2. These three features ? irregular trabecular density, irregular anisotropy properties, and irregular cortical thickness ? create significant obstacles to biofidelic computer modelling of bone. Researchers have used two competing modelling methodologies to negotiate these obstacles on FE models of bone. The first and most common method uses relatively large elements (with 1 ? Introduction 12 sides around 1 mm or larger) and applies heterogeneous density-based material properties to them. The second method uses much smaller elements (with sides typically smaller than 100 ?m) and simplified, homogeneous material properties. In essence, there is a trade-off between the levels of detail used to model the material and geometry of bone. In the first method, QCT data of the specimen is typically collected at clinical resolution (around 0.6 mm). Images at this resolution do not capture much of the detail in trabecular structure or thin layers of cortical bone; however, it has been shown that the grey-level voxel values in QCT are linearly related to the ash density of bone [70], [71]. Several studies have probed the mechanical properties, particularly the stiffness, of trabecular bone of varying ash densities [57], [71]?[74]; allowing the conversion from CT value to ash density and from ash density to local stiffness. Almost all the FE models described in section 1.6 were created by averaging local stiffnesses over the volume of each element. Note, however, that using this technique, material properties inside each element are typically homogeneous and isotropic. Some groups have studied the use of anisotropic properties, but the problem of finding the principal directions and degree of anisotropy appropriate for each element is far from trivial [67], [68]. Furthermore, this technique\t ?suffers\t ?from\t ?partial\t ?volume\t ?artifacts,\t ?where\t ?boundary\t ?voxels\t ??dissolve?\t ?thin layers of cortical bone in the surrounding fluid (typically water or air, depending on the technique used to scan the bone), artificially weakening the virtual structure [75]. Finally, this type of model is insensitive to small defects like holes or surface irregularities on the bone, making it possible to miss stress concentrations during FE analysis. In the second method, QCT data is collected at much higher resolutions (with voxels as small as 40 ?m or even less),\t ?obtaining\t ?a\t ?detailed\t ?map\t ?of\t ?the\t ?specimen?s\t ?trabecular\t ?structure.\t ?Voxels\t ?can\t ?then be categorized as either bone or not bone through the use of a threshold CT value. Bone voxels are directly converted to prismatic elements of the same size, and these elements are assigned equal homogeneous material properties independent of CT value. It is assumed here that, at this level of detail, real bone material is also homogeneous and isotropic. FE models that use this technique typically contain millions of elements, requiring highly specialized solvers and powerful supercomputers [76]. Some researchers have chosen to model small regions of bone in order to decrease the number of elements and the processing time required. However, the process of virtually isolating these small regions generally requires slicing the inside of a bone, and here it becomes very difficult to determine what boundary conditions are appropriate for biofidelic loading. 1 ? Introduction 13 Our lab is currently testing the idea of combining both approaches described above. The idea relies on using two types of\t ?model:\t ?an\t ??organ-level?\t ?model\t ?that\t ?follows\t ?the\t ?first\t ?method\t ?with\t ?large elements and density-based properties to model the entire proximal femur, and one or more\t ??micro-level?\t ?models\t ?that\t ?follow the second approach to model only small regions of interest in the bone. Specimen data for both models can be collected as a high-resolution CT scan that can be downsized to clinical resolution for the organ-level model. The problem of internal boundary conditions in the micro-level models is resolved by applying the solved nodal displacements of the organ-level model as boundary conditions on the micro-level models. The advantage of this technique is that one could use the organ-level model to identify regions where high strains develop, and then model these regions with micro-level models to determine whether any stress concentrations would make a fracture likely to be initiated. 1.8 Objectives The overall objective of this investigation is to explore the use of FE analysis to simulate a femur under sideways fall loading. The specific aim is to model some features of sideways fall loading ? particularly the effects of a dynamic loading rate ? that have not been considered in previous FE studies of the subject. One final objective is to probe the feasibility of multi-scale nested FE modelling, a novel technique designed to simulate the behaviour of bone at a scale that resolves its trabecular structure without incurring the great computational expense of high-resolution FE modelling of an entire proximal femur. The investigation is divided in three parts corresponding to the main chapters of this thesis. Each chapter follows a manuscript format, intended for later publication as a standalone article. x In chapter 2, I use quasi-static linear FE modelling, with the intention of determining the appropriateness of several model components and establishing the capabilities of our simulation technique in a relatively simple loading scenario similar to that used by others in the research community. x In chapter 3, I use explicit dynamic nonlinear FE modelling, a technique that is novel for hip fracture research, to develop a model of the femur under the loading rates and conditions observed in sideways falls. x In chapter 4, I present a preliminary investigation of multi-scale nested FE modelling. 1 ? Introduction 14 All three chapters include results comparisons to measurements from appropriate mechanical tests of human femurs. 1.9 Scope This thesis presents a portion of a larger investigation on hip fractures performed collaboratively by a team of researchers at the University of British Columbia (UBC) and the Swiss Federal Institute of Technology (ETH) Zurich. Some components of this investigation are necessary for the clear presentation and validation of this research, yet a detailed analysis of all aspects of the collaborative project is not within the scope of this thesis. The FE results presented here are compared against data from empirical ex-vivo fall simulations and mechanical testing of human femurs. The validity of the results for clinical hip fractures is therefore dependent on the validity of the empirical data itself. However, the validation of these data against clinical hip fractures falls outside the scope of this thesis. For further information on the experimental side\t ?of\t ?our\t ?group?s\t ?research,\t ?please\t ?refer\t ?to\t ?Gilchrist?s\t ?PhD\t ?thesis [77]. 15 2 Quasi-static linear simulation 2.1 Introduction Osteoporosis is characterized by a loss of bone density and is known to result in elevated fracture risk of the hip, spine and forearm, with hip fractures resulting in the most severe consequences. Epidemiologic evidence supports the contention that low areal bone mineral density (aBMD) is associated with increased population-based risk of fracture [7], [78]. However, there is considerable overlap of aBMD values between fracture and non-fracture patients, indicating that the method is not sensitive enough to identify many of the individuals likely to suffer a fracture [8], [9]. Subject-specific finite element (FE) models based on quantitative computed tomography (QCT) have been created in a number of studies to analyse femoral behaviour in single-leg stance [12]?[14], [51]?[57] and sideways fall [12], [52], [54], [58]?[60] loading configurations. Most of these models have been validated against ex-vivo mechanical tests and some have been found to yield better predictions than aBMD [51], [58], underscoring that these models could be ideal for identifying individuals at high risk. Despite these positive advances, no agreement on a consistent FE methodology has emerged in the hip fracture biomechanics community. Published models have used a variety of modulus-density relationships to estimate heterogeneous material properties from QCT-derived density data. Models have also used a variety of methods to define bone and element geometry, with some researchers using voxel models that directly match the QCT voxel grid [12], [51], [52], [54] and others using irregular tetrahedral or hexahedral isoparametric elements meshed over a smoothed solid geometry [13], [14], [53], [55]?[60]. Smooth-geometry models further differ in the strategies used to distribute ? or map ? material properties from a regular grid of QCT voxel values to an irregular mesh of elements. Given the lack of agreement on the appropriate materials and methodologies required for femoral FE modelling, each research group has generally developed its own technique and determined the material properties that best work with its models, validating them against experimental data from mechanical tests. A few published studies have produced FE results using multiple material-density relationships and compared them against mechanically derived data from a few femoral specimens [55], [57], [64]. Limited comparisons of material mapping methods [79] and meshing techniques [80] have also been published. All of these comparative studies lack validation against a sizeable number 2 ? Quasi-static linear simulation 16 of independent specimens. It is our objective to perform a systematic comparison of material relationships and modelling methods using validation data from a wide variety of mechanically tested specimens. We use a semi-automated and modular FE modelling system that allows the isolation and switching of specific portions of our methodology. The automation in this scheme allows us to easily build and run equivalent FE models for a range of specimens, while the modularity gives us the ability to analyse the effect of changes in isolated portions of the model. The resulting comparisons should allow us to separate the effects of modelling methodology from those of inter-specimen variability. We intend to use this strategy, validated with data from sixteen mechanically tested specimens, to compare the results of multiple material mapping strategies and modulus-density relationships. For our work, we choose to work with sideways fall loading instead of single-leg stance loading. We choose this because the majority (>90%) of hip fractures are associated with a low trauma fall [81], and falls to the side have a high fracture rate [23]. The aim of our hip fracture FE simulations is to recreate virtually the conditions of a sideways fall. While we believe that the role of dynamic behaviour on the femur response to a fall impact needs to be investigated, and it is covered in the next chapter of this thesis, we need to first isolate the mechanical response of femurs in a simplified scenario. An FE model of the femur loaded quasi-statically and validated against data from slow mechanical tests allows us to investigate the basic material properties of the femur without the additional complexity of dynamics. We choose to work with FE meshes using smoothed solid geometry, since these types of meshes have been shown to produce more accurate results than voxel-type meshes at similar densities [80]. However, smoothed geometries require mapping of heterogeneous material properties from a regular grid of QCT voxel to an irregular mesh. The most common mapping strategy is to average CT-derived values over the volume of each element. This strategy, specifically averaging voxel-based Young?s\t ?moduli,\t ?is\t ?used\t ?by\t ?the popular BoneMat V3 software [79]. An alternative strategy is to use trilinear interpolation on the CT values to find appropriate material properties at each mesh node, creating in effect a varying field of nodal properties [82]. One final mapping strategy found in published studies uses higher order (p-method) elements, allowing the application of element properties as continuous spatial functions [83]. This last strategy, however, is excluded from our current comparison, since it requires the use of specialized meshes and solvers. 2 ? Quasi-static linear simulation 17 For material models, we reviewed published modulus-density relationships developed and validated using mechanical test data and chose those that have received special attention. One relationship in particular, was developed by Morgan, Bayraktar and Keaveny [72] using a very robust experimental protocol. This relationship has been used by other research groups who have found it to produce superior surface strain predictions in a single-leg stance [55] and sideways fall [59] loading configurations. However, it has also been found to systematically overestimate femoral stiffness under sideways fall despite a good predictive correlation [58]. Another relationship, presented by Keyak et al [12], combined two previous relationships for human trabecular bone over a wide density range. Studies using this relationship in single-leg stance loading have produced contrasting results, with one reporting accurate stiffness predictions and somewhat overestimated strains [13] and another reporting very good strain agreement and overestimated stiffness [56]. Further research by Wille, Rank, and Yosibash [57] also found very good agreement for surface strains using the Keyak relationship. This last paper also developed a novel stochastic modulus-density relationship derived from pooled femoral bone data from a variety of studies. The accompanying deterministic relationship performed well on our preliminary femoral FE modelling tests, encouraging us to investigate it further. The three relationships described above: Morgan [72], Keyak [12], and Wille [57], are compared in this chapter. The investigation described below was performed as a collaborative project between the Orthopaedic and Injury Biomechanics Group at the University of British Columbia (OIBG-UBC) and the Institute for Biomechanics at the Swiss Federal Institute of Technology in Zurich (IfB-ETHZ). The contribution of Dr. Benedikt Helgason from IfB was particularly instrumental to the development of this investigation, as he not only supervised my work while writing the FE models but also helped run them and summarize the resulting data. 2.2 Methodology Sixteen proximal femurs were mechanically tested non-destructively on a materials testing machine. Specimen-specific finite element (FE) models were created and solved using an implicit solver in ANSYS Mechanical APDL (v14.0). Stiffness and surface strains on the anterior neck resulting from the FE model were compared to the values measured experimentally. 2 ? Quasi-static linear simulation 18 2.2.1 Experimental setup The methodology used to collect the experimental data for this project is described by Gilchrist [41] as validation for digital image correlation. This section provides a brief overview of Gilchrist?s\t ?experimental\t ?setup. Sixteen fresh frozen proximal femurs were evaluated (Table 2-1). The specimens were specified to have no reported history of musculoskeletal or metabolic disease. The bones, which had their soft tissues removed, were scanned using dual energy X-ray absorptiometry (DXA) (QDR 4500W, Hologic, Bedford, MA) and high resolution (41 ?m) peripheral quantitative computed tomography (XtremeCT, SCANCO Medical AG, Br?ttisellen, Switzerland). Table 2-1: General data for femoral specimens Specimen Age/Gender (years/-) Height/Weight (m/kg) Femoral Neck aBMD (g/cm2) Total aBMD (g/cm2) T-score WHO classification H1167L 50/F 1.68/54 0.658 0.841 -0.8 Normal H1168R 73/M 1.70/77 0.650 0.813 -1.1 Osteopenic H1268L 82/F NA/NA 0.536 0.698 -2.0 Osteopenic H1365R 71/F 1.52/45 0.504 0.638 -2.5 Osteoporotic H1366R 73/F NA/45 0.740 0.839 -0.8 Normal H1368R 70/F 1.65/57 0.500 0.606 -2.8 Osteoporotic H1369L 69/F 1.73/68 0.498 0.601 -2.8 Osteoporotic H1372R 79/F 1.68/113 0.638 0.793 -1.2 Osteopenic H1373R 76/F 1.52/63 0.773 0.833 -0.9 Normal H1374R 78/F 1.45/66 0.486 0.621 -2.6 Osteoporotic H1375L 83/F 1.52/90 0.633 0.782 -1.3 Osteopenic H1376L 79/F 1.63/64 0.552 0.683 -2.1 Osteopenic H1377R 80/F 1.65/64 0.503 0.627 -2.6 Osteoporotic H1380R 71/F NA/NA 0.695 0.834 -0.9 Normal H1381R 92/F 1.63/51 0.392 0.545 -3.3 Osteoporotic H1382L 96/F 1.52/61 0.441 0.546 -3.2 Osteoporotic Each bone was tested mechanically on a materials testing machine (8874, Instron, Norwood, MA). The general arrangement during mechanical testing is shown in Figure 2-1. The distal end of each femur was potted in an aluminum cylinder using poly(methyl methacrylate) (PMMA). The cylinder and its positioning fixture were designed to only allow rotation around a distal anterior-posterior axis. The positioning fixture held the bone in a position that mimics the position thought to occur during a sideways fall [38], [39]. The femoral head was constrained medially, but allowed to displace in the sagittal plane using a steel plate supported on free-sliding bearing plates. The greater trochanter was loaded in a lateral to medial direction. Forces 2 ? Quasi-static linear simulation 19 at both the head and the greater trochanter were distributed over the surfaces using PMMA pads, approximately 4 cm in diameter, to prevent local crushing. The pads were moulded to both bone and plate surfaces but not glued to them. A single strain rosette (FRA-2-11-3LT, Tokyo Sokki Kenkyujo Co., Tokyo, Japan) was attached to the anterior-superior surface of the femoral neck. This area was selected for monitoring because it was a common site of fracture initiation on earlier research performed in our lab [26]. After placing the bone in the positioning fixture, an Optotrak Certus motion capture system (Northern Digital Inc, Waterloo, ON, Canada, accuracy < 0.1 mm) was used to record the 3D locations of a number of points on the surface of the bone and several areas of the positioning fixture. These digitized data were used to accurately reproduce loading directions and constraint locations on the computer models, as described in section 2.2.4. Each bone was loaded at quasi-static loading rates (~ 0.5 mm/s) up to half of its fracture load, as predicted from the DXA-based areal bone mineral density (aBMD) [29]. Forces and displacements measured by the materials testing machine were recorded using LabView (v2010 SP1, National Instruments, Austin, TX). Surface strains were measured locally at the anterior-superior neck by the strain rosette described above. Figure 2-1: Sample femur placed in positioning fixture 2 ? Quasi-static linear simulation 20 2.2.2 FE model geometry The high-resolution peripheral QCT data was calibrated using a standard hydroxyapatite (HA) phantom and filtered using a morphological filter to isolate the bone from intertrabecular pores. This step was necessary to avoid artificially stiffening the model by using a bone modulus-density relationship on the noisy data present in marrow-filled pores. The images were then reduced from their original 41 ?m resolution to clinical resolution of 0.615 mm in Matlab (v R2012b, The Mathworks, Natick, MA, US). In order to maintain a constant amount of bone during the downsizing process, the average intensity of each voxel was kept constant. The reduced images were segmented using ITK-SNAP (v2.2.0) [84], with the resulting data converted to a NURBS-based solid in SolidWorks (2011, Dassault Systems, Concord, MA). One of the bones was selected arbitrarily to generate a series of 2nd order tetrahedral meshes from the NURBS data using Ansys Workbench (v.14.0, Altair Engineering Inc., Troy, MI, USA). These meshes were tested with the FE analysis described below to determine the\t ?model?s\t ?sensitivity to mesh size. The test, which compared resultant forces and stiffnesses from meshes with 1.5, 2.0, 3.0, and 4.0 mm mean edge length, concluded that a model with average element edge length around 2.0 mm yields converged results in a reasonable computation time. All bones were then meshed using similarly-sized elements. Summary mesh data is presented on Table 2-2. On each of the meshes, a node was placed at the location of the strain gauge according to the points digitized with the motion capture system. Table 2-2: Mesh characteristics for each bone Bone Volume (cm3) Elements Nodes H1167L 129.5 133 450 192 385 H1168R 168.5 170 777 244 504 H1268L 119.9 124 709 180 308 H1365R 105.6 108 935 157 757 H1366R 143.7 146 567 211 012 H1368R 147.1 149 390 214 701 H1369L 171.6 171 631 245 697 H1372R 121.1 212 195 296 728 H1373R 154.1 155 939 223 925 H1374R 126.4 131 452 189 801 H1375L 132.0 134 675 193 909 H1376L 119.3 122 806 177 454 H1377R 142.5 249 667 348 407 H1380R 124.8 128 428 185 447 H1381R 162.0 164 833 236 347 H1382L 147.1 149 726 214 932 2 ? Quasi-static linear simulation 21 2.2.3 Heterogeneous material stiffness mapping The voxel values from the reduced images, which represent hydroxyapatite (HA) content, were converted to ash density values by the following relationship: ???? = ????? + 0.09?/1.14 (2-1) where ?ash represents the ash density of the bone voxel in g/cm3, and ?QCT represents the corrected QCT voxel intensity in g HA/cm3. Apparent density was then calculated from the resulting voxel ash densities using the relationship ????????= 0.60 (2-2) Both of these density relationships are based on the investigation of Schileo et al [85]. Dr. Benedikt Helgason, a main collaborator in the present investigation and my co-supervisor in ETH Zurich, previously examined the suitability of eight bone density-stiffness relationships presented in published studies on FE models of the proximal femur. From this early investigation, three relationships were selected for the present study (Table 2-3 and Figure 2-2). Each relationship was used to calculate the stiffness value of each CT voxel. Table 2-3: Selected modulus-density relationships Short Name Modulus-Density Relationships (MPa) Densito-metric range Source material Morgan 49.18506 appE U Whole Range [72] Wille 45.100012 ashE U Whole Range [57] Keyak 20.290033 ashE U 4693075 \u000E ashE U 01.220010 ashE U 27.0dashU 60.027.0 \u001F\u001F ashU ashUd6.0 [12] 2 ? Quasi-static linear simulation 22 Surface voxels, which are at the boundary between bone tissue and the surrounding environment, were likely to suffer from partial volume artifacts: reduced properties due to the voxel volume being filled only partially with bone. One method that can be used to correct this is to\t ??peel? the outermost layer of bone voxels. The values of these voxels can then be replaced with those from the next inner layer if these were higher, in effect extruding the last layer of voxels that contain only bone. To test the influence of the peeling step, each bone was assigned properties with and without peeling. Having\t ?corrected\t ?Young?s\t ?moduli\t ?assigned\t ?to\t ?each\t ?CT\t ?voxel,\t ?the\t ?stiffness\t ?at\t ?each\t ?node\t ?in\t ?the\t ?mesh was then calculated by 3D linear interpolation. These nodal stiffnesses could be applied directly\t ?to\t ?the\t ?FE\t ?model\t ?using\t ?Helgason?s\t ?nodal\t ?interpolation (NI) technique [82], or they could be distributed to nearby elements to obtain more traditional element stiffnesses. Both techniques were tested as described below. Nodal stiffnesses were distributed to element stiffnesses using weighted averages with nodal weight factors that were inversely proportional to the distance from node to element centroid. The following three material mapping methods were selected for testing: Figure 2-2: Comparison of modulus-density relationships over bone density range 2 ? Quasi-static linear simulation 23 x Method 1: Element stiffnesses and no peeling correction. This case was selected for its similarity to BoneMat V3 used by other researchers [79]. x Method 2: Element stiffnesses with peeling correction. This is an intermediate case that can be easily implemented in almost any FE software. x Method 3: Nodal stiffnesses with peeling correction. This case tests the influence of Helgason?s\t ?NI\t ?method [82]. Each of the three cases was tested with each of the three modulus-density relationships. All simulations used linear elastic material properties. 2.2.4 Boundary conditions The positions, constraints, and loading directions of Gilchrist?s\t ?physical loading test (described in section 2.2.1) were closely mimicked digitally using the points digitized with the motion capture system. The digitized points included arbitrary locations over the surface of the femur and specific locations in the supporting and loading devices (Figure 2-3). The points in the bone surface together with surface nodes in the CT-derived mesh were run through an iterative closest point algorithm [86] in Matlab to register the conversion between the CT coordinates used by the mesh and the digitizer coordinates used to record boundary condition data. The rest of the digitized points were then used to calculate the locations and directions of boundary conditions to match the experimental setup. Figure 2-3: Example of digitized points used for registration (The registered CT-derived mesh and calculated geometry are shown for clarity) 2 ? Quasi-static linear simulation 24 The distal (inferior) end of the model was constrained in such a way that the bone could only rotate around a distal hinge on the frontal (coronal) plane. This was accomplished by constraining a node at the midpoint of the distal hinge. This hinge node was connected to the distal end of the solid femur model using beam elements that emulate the axial and bending stiffness of the PMMA-filled aluminum cylinder with which the real femoral shaft was potted. To simulate reality as closely as possible, three beams were used in series: The most proximal beam used dense bone properties and circular tube geometry to simulate the distal portion of the femoral shaft that was outside the CT scanner range and thus could not be modelled as a solid. The middle beam was given composite properties corresponding to a PMMA-filled aluminum cylinder, emulating the potted portion of the cylinder. The most distal beam used the properties of a hollow aluminum cylinder to simulate the portion beyond the potting. In order to simulate the experiment, the femoral head required a constraint that prevented medial movement without hindering its ability to rotate in any direction or translate in the other two axes. It was decided that the most appropriate method was to model a PMMA support pad, constrain its medial end, and apply a frictionless contact between it and the femoral head. The pad used the same 2nd order tetrahedral elements used in the femur, elastic material properties with E = 3.1 GPa, and a geometry that matched the femoral head at the contact surface. Finally, the greater trochanter was loaded medially through a PMMA loading pad model. Similarly to the head constraint, a frictionless contact was set up between PMMA and bone. The lateral end of the loading pad was loaded using a rigid plate over the PMMA pad and applying a downward displacement of 1.0 mm (Figure 2-4). The use of PMMA pads for constraints closely matched the experimental setup, which had detachable moulded PMMA pads at the same two locations. 2 ? Quasi-static linear simulation 25 2.2.5 Comparison criteria The comparison of FE to experimental results was done using two main parameters: overall bone stiffness ? measured as the displacement of the greater trochanter in the loading direction ? and principal surface strains at the strain gauge location. Overall\t ?stiffness\t ?was\t ?calculated\t ?experimentally\t ?from\t ?the\t ?linear\t ?portion\t ?of\t ?each\t ?bone?s\t ?force-displacement curve. Force was measured by the load cell in the materials testing machine and filtered in Matlab using a 2nd order Butterworth filter in forward and reverse directions with a 500Hz cut-off. Displacement was similarly measured by the machine and corrected for both the deformation due to machine compliance ? which Gilchrist measured directly on a separate test ? and the deformation of the cartilage layer around the femoral head ? based on an independent cartilage model developed by Helgason. On the numerical side, overall FE stiffness was calculated from the reaction force at the head support and the displacement applied at the greater trochanter. Surface strains were measured experimentally using the strain rosette. On the FE model, a node was placed at the digitized strain gauge location of every mesh. For the comparison to strain gauge measurements, FE strains were averaged among all surface corner nodes within 0.9 mm of this location. This radius accounted for the actual size of the strain gauge. Figure 2-4: Sample FE model and boundary conditions 2 ? Quasi-static linear simulation 26 2.3 Results All nine simulations ? for three material mapping methods and three modulus-density relationships ? were successfully run on each of the sixteen bones. The results from these simulations were compared to experimental values, resulting in correlation coefficients between 0.71 and 0.79 for stiffness and between 0.90 and 0.92 for strains. Figure 2-6 shows the validation results for stiffness, while Figure 2-7 shows the results for surface strains at the strain gauge location. When comparing overall bone stiffness (as defined in Figure 2-5),\t ?the\t ?Keyak\t ?relationship?s\t ?results produced regression slopes closest to unity and the lowest root-mean-squared error ratio among the three relationships, making this relationship the most accurate among the three for overall stiffness. All sets of stiffness results produced R2 values between 0.7 and 0.8. When comparing strains, the\t ?Morgan\t ?relationship?s\t ?results\t ?produced\t ?regression\t ?slopes\t ?closest\t ?to unity and the lowest root-mean-squared error ratio among the three relationships, making this relationship the most accurate among the three for surface strains. All sets of strain results produced R2 values greater than 0.9. This finding agrees with the results of Schileo [55], who determined that surface strains were best predicted using the Morgan relationship. Note, however,\t ?that\t ?Schileo?s\t ?study\t ?investigated\t ?femurs\t ?under\t ?single-leg stance rather than the sideways fall loading configuration of the present study. Figure 2-5:\t ?Definition\t ?of\t ?femoral\t ?stiffness.\t ?F\t ?represents\t ?the\t ?applied\t ?force\t ?and\t ??\t ?the\t ?displacement of the greater trochanter 2 ? Quasi-static linear simulation 27 MORGAN WILLE KEYAK METHOD 1 METHOD 2 METHOD 3 y = 1.4356x + 0.0846R? = 0.7915RMSE = 29.57%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmy = 1.2698x + 0.1348R? = 0.7201RMSE = 22.85%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmy = 1.0479x - 0.1885R? = 0.7253RMSE = 11.15%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmy = 1.4976x + 0.2184R? = 0.7174RMSE = 37.0%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmy = 1.2926x + 0.2105R? = 0.7136RMSE = 25.61%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmy = 1.0896x - 0.1432R? = 0.7088RMSE = 12.07%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmy = 1.4699x + 0.1275R? = 0.7239RMSE = 33.49%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmy = 1.2698x + 0.1348R? = 0.7201RMSE = 22.85%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmy = 1.0561x - 0.1841R? = 0.7117RMSE = 11.55%123456781 2 3 4 5 6 7 8FE Stiffness, kN/mmExperimental Stiffness, kN/mmTrendline 1:1 MatchFigure 2-6: Comparison of experimental to FE stiffness results for each modulus-density relationship and material mapping method 2 ? Quasi-static linear simulation 28 MORGAN WILLE KEYAK METHOD 1 METHOD 2 METHOD 3 2.4 Discussion The aim of the present study was to validate FE-models of the proximal femur in a sideways fall configuration against experimental test results for a broad range of specimens. The validation was performed in terms of overall specimen stiffness as well as measurements from strain gauges mounted at critical locations on the femoral neck. y = 1.153x + 0.0513R? = 0.9074RMSE = 17.06%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)y = 1.2191x + 0.0773R? = 0.9153RMSE = 18.58%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)y = 1.9916x + 0.0806R? = 0.9017RMSE = 56.16%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)y = 1.0475x + 0.0424R? = 0.9082RMSE = 14.20%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)y = 1.2461x + 0.0617R? = 0.9092RMSE = 20.24%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)y = 1.8411x + 0.0565R? = 0.9000RMSE = 49.00%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)y = 1.0232x + 0.0566R? = 0.9145RMSE = 13.28%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)y = 1.2191x + 0.0773R? = 0.9153RMSE = 18.58%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)y = 1.8123x + 0.0663R? = 0.9051RMSE = 47.04%-8-6-4-202-8 -6 -4 -2 0 2FE Strains, (x10-3 )Experimental Strains, (x10-3)Trendline 1:1 MatchFigure 2-7: Comparison of experimental to FE strain results for each modulus-density relationship and material mapping method 2 ? Quasi-static linear simulation 29 Two key portions of the FE model formed the focus of our comparison: 1) modulus-density relationship and 2) material mapping method. The comparison between the three modulus-density relationships ? here simply referred to as Morgan, Wille, and Keyak ? is discussed first. Over almost all of the bone density range, Morgan is the stiffest relationship and Keyak the most compliant (Figure 2-2).\t ?Since\t ?models\t ?were\t ?developed\t ?from\t ?the\t ?same\t ?CT\t ?values,\t ?Young?s\t ?moduli\t ?were generally higher in the Morgan simulations and lower in the Keyak simulations, with Wille?s\t ?values\t ?somewhere\t ?in\t ?between.\t ?As\t ?should\t ?be\t ?expected,\t ?this\t ?difference\t ?between\t ?moduli\t ?at\t ?the local level produced similar differences in overall bone stiffness, with the Keyak relationship producing the lowest overall stiffness and the other two relationships producing higher ? and overestimated ? stiffness results. The relative comparison slopes obtained for surface strains from the different material relationships can also be traced to the general differences in stiffness. The Morgan relationship produced the lowest strains, as would be expected from the stiffest material model, and thus the lowest slope. The other two relationships, being more compliant, produced higher maximum (more positive) and minimum (more negative) principal strains. It follows that the slopes of the experimental-to-FE comparison are higher for these two relationships. Our R2 values ranged from 0.71 to 0.92, which are similar to those obtained by other authors in comparable studies ? e.g. 0.91 by Grassi et al [59], 0.87 by Dragomir-Daescu et al [58]. Our study used a sizeable sample of specimens (bones from 16 individuals) but few strain measurements per specimen (two strains at a single location for each bone). This is in contrast to other investigations that have compared FE and experimental strains, like Grassi et al [59] (16 locations in 3 bones) or Trabelsi et al [56] (5 locations on 6 bone pairs). Our approach produces statistical results with generally less power than the alternative, but it ensures independence between measurements. Perhaps more telling than the correlation coefficients are the root-mean-squared errors (RMSEs) reported as a percentage of the maximum experimental measurement. In our results, these vary between 11.2 and 56.0%, showing the relative accuracy of different model types. This metric shows more clearly the ability of the model to predict strains using the Morgan material relationship and stiffness using Keyak. The Wille relationship is also shown to be a compromise for both stiffness and strain prediction. 2 ? Quasi-static linear simulation 30 Comparing the different material mapping methods on both stiffness and strain results did not produce enough variation in results to deem any of the methods superior to the rest. The difference between mapping methods is clearly overshadowed by the difference between modulus-density relationships. Any difference between the accuracy of the three methods likely requires a larger sample to be resolved. The studies by Grassi et al [59] and the pioneering work of Lotz et al [11] are the only other published research we found that include local strain validation of a FE femoral model under sideways fall loading. Grassi et al used FE models similar to the present model with the Morgan relationship and material mapping method 1. Their investigation, which explored surface strain predictions on three bones using 16 strain measurement locations and 12 loading configurations, found very good correlation and low errors between strain measurements and predictions. Lotz et al, who similarly instrumented two bones with nine strain gauges each, found a rather poor correlation between FE strains and those measured from strain gauges. On the stiffness side, Dragomir-Daescu et al [58], using data from 18 femurs, also found that the Morgan relationship leads to an overestimation of overall bone stiffness. The lower degree of overestimation found in the results here presented (regression slope around 1.5, compared to their 2.0) may be explained by the differences in FE boundary conditions, which in our model were purposely set to allow free rotation of the femoral head. Our results indicate that the Keyak material relationship leads to the best predictions of overall stiffness under sideways fall loading. To the best of our knowledge, no previous studies have found a better predictor of bone stiffness under sideways fall loading than that shown here. The ultimate goal of femoral FE research is, in general, to develop a modelling methodology that can be implemented clinically for accurate fracture risk assessment. Successful validation against cohort data is essential for further advancement toward this goal. We have found, however, that previous studies that include this validation often use voxel-based FE modelling, a technique that has been shown to be less accurate than smoothed-solid modelling [80]. Furthermore, previous studies have sometimes not given enough emphasis to boundary conditions. During our investigation, we found that, under a fall configuration, a femoral FE model can be easily overconstrained. Two boundary condition features in particular demonstrated surprising importance for our models: 1) ensuring rotational freedom for the femoral head through a sliding surface contact, and 2) proper representation of the bending and torsional stiffnesses of the distal support mechanism. Failure to account for either of these can 2 ? Quasi-static linear simulation 31 lead to displacement errors of perhaps several hundredths of a millimetre, which could significantly alter the resulting stiffness and strains in the model. We found that the Wille modulus-density relationship, while inferior to Morgan for surface strain predictions and inferior to Keyak for overall stiffness predictions, produces reasonable estimates of both stiffness and strains. For this reason, this relationship is recommended for further study and use in FE simulations of the proximal femur. It is important to note, however, that the experimental strains are likely to have better accuracy than the experimental deformations ? which required several estimated corrections ? and thus the overall stiffnesses. For this reason, the stiffness results carry somewhat less weight than the strains. The Morgan relationship should also be considered where appropriate. As far as we are aware, the present study is the first to apply the Wille relationship on a sideways fall scenario. In our results, the most compliant modulus-density relationship (Keyak) produced the most accurate overall stiffness, yet the stiffest relationship (Morgan) produced the most accurate surface strains. These results may point toward the need to model a stiffer layer of cortical bone. It can be reasoned that a stiffer cortical bone would yield lower and more accurate surface strains; while more compliant trabecular bone may reduce the overall stiffness of the bone, yielding more accurate results. The differentiation between trabecular and cortical models is supported by recent findings that suggest that cortical bone and trabecular bone perform distinct roles in the mechanical behaviour of the femur. While cortical bone typically carries the majority of the load in both standing [34] and fall [87] configurations, trabecular tissue appears to provide the structural support necessary to engage the cortex during a fall [88]. There are a number of limitations to this investigation. Given the small number of specimens, it is possible and perhaps even likely that the tested bones are not representative of the overall population. Similarly, the fact that surface strains were compared using a single strain rosette for each bone makes it likely to miss strain patterns and features present in the surface of the neck. Future studies would benefit from testing more specimens and comparing more locations for strain. This could be done with more strain gauges or, as our lab is already investigating, through the use of digital image correlation to measure strains from video data [41]. The measurement of overall stiffness has limitations as well, since it is based on deformation data that requires a number of estimations to account for the compliance of the head cartilage and non-rigid support and loading structures. While we found these corrections to be important, they could add small inaccuracies to the displacements and thus to the stiffness too. 2 ? Quasi-static linear simulation 32 Additionally, there are some limitations to the presented methodology, including the lack of specific validation data for the properties used for the distal beams. Similarly, only isotropic properties were used for the bone material ? a simplification used by most comparable studies [12], [55], [57]?[59] ? despite the well-known anisotropic behaviour of bone. In summary, the presented methodology supports the use of the Wille and Morgan relationships on femoral FE models under a sideways fall configuration. Not enough differences were noted from the three presented material mapping methods to recommend any particular one. Finally, we recommend that particular attention be paid to the boundary conditions of femoral FE models to avoid overconstraining the model. 33 3 Dynamic nonlinear simulation 3.1 Introduction It is estimated that each year over 1.6 million people suffer a hip fracture worldwide [1]. Beyond its devastating effects on mobility and quality of life [89]?[93], this type of injury has been associated with increased risk of death [94], [95]. There is ample evidence showing that the majority (>90%) of hip fractures are associated with a short, low trauma fall, and while falls in elderly adults are common [96], only 1-5% of falls result in a fracture [24], [97], [98], with falls to the side having a much higher fracture rate [23]. The biomechanics of hip fractures associated with a sideways fall have been empirically studied, most often by applying forces using materials testing machines [12], [26], [52], [58], [99]. Two aspects of this method limit its capacity to represent a sideways fall. First, loading rates of 0.1 m/s or lower have been generally used, while the impact speed associated with a sideways fall from standing is around 3.0 m/s [39]; this rate is beyond the reach of most current materials testing machines. Loading rate is important, as it has been shown that the mechanical properties of bone are strain rate dependent [100], [101]. The second limitation is that load has generally been made to increase monotonically until fracture occurs using a materials testing machine. In contrast to the limited energy available during in sideways fall, a machine like this delivers any amount of energy necessary to fracture the bone ? i.e. infinite energy is available to produce the fracture. This represents failure in an overload condition, where the load at which the bone fractures may greatly exceed the load related to a physiological fall from standing height. In other words, this method probes the load carrying capacity of a femur on its side rather than set a scenario equivalent to a sideways fall. Furthermore, given the heterogeneity and rate dependence associated with the mechanical properties of bone, loading paths could be quite different in the dynamic and quasi-static cases. In addition to utilizing these experimental methods, many investigators have studied the mechanics of the proximal femur using subject-specific finite element models based on X-ray computed tomography (CT) [11]?[13], [51]?[53], [55], [57]?[59], [65], [83], [85]. These models have been tested with a variety of material modelling methods and loading configurations. However, they have almost invariably followed a quasi-static structural approach and been validated against quasi-static in-vitro models. While this simplified approach has certainly led to improved understanding of femoral behaviour under normal and abnormal loading, there are 3 ? Dynamic nonlinear simulation 34 significant limitations to neglecting the influence of body inertia, soft-tissue damping, rate dependent material response, and other dynamic effects present in-vivo during a fall [61]. Given the parameters available during a quasi-static investigation of the proximal femur, researchers have most often chosen to regard hip fracture risk as a function of only femoral strength ? i.e. the experimentally or numerically-derived loading capacity of the bone [13], [65], [102]. In a literature review, Hayes, Piazza, and Zysset [103] demonstrated a viable method for fracture risk prediction from static strength data and also underscored the difficulties associated with it. A different approach would calculate the energy available before a mechanical test ? e.g. using a free-falling mass in a drop tower setup ? then account for energy spent until the point of fracture initiation; with this information and multiple tests, it would be possible to estimate the amount of energy a particular femur can absorb before it fractures (see Appendix for details). Such an approach could improve understanding and prediction of fracture risk. This approach would incorporate the ability of a femur to absorb and dissipate energy through elastic and plastic deformation before suffering total collapse. Only dynamic modelling can provide the necessary data to test the viability of such an approach. Unfortunately, dynamic modelling is, to the best of our knowledge, currently missing in hip fracture literature. To address the current limitations in experimental and computational simulations of sideways falls on the hip we propose to use the drop tower methodology developed by Gilchrist [41] to apply loads at speeds matching those of actual falls and to model the same scenarios using strain-rate dependent explicit FE simulations. The results are compared in terms of overall stiffness, ultimate force and energy absorbed. 3.2 Methodology Fifteen proximal femurs were destructively tested by Gilchrist in a fall simulator. These bones were the same set previously tested non-destructively (Chapter 2), excluding data from H1268L. During impact testing of this bone, there was an accidental contact that rendered that test invalid. We normally expect that only 1-5% of falls result in a fracture [24], [97], [98], however, for the purpose of methodological development, the specimens tested in the present study were subjected to input energies that would almost surely fracture them in order to verify whether the drop tower setup could induce clinically relevant fracture patterns. Specimen-specific finite element models were created and solved using the explicit LS-DYNA solver in ANSYS Mechanical APDL (v14.0). Overall stiffness, maximum force, absorbed energy, and surface strains patterns were compared to the values measured experimentally. 3 ? Dynamic nonlinear simulation 35 3.2.1 Experimental setup The methodology used to collect the experimental data for this project was described in detail by Gilchrist [41]. This section provides a brief overview of his experimental setup. Fifteen fresh frozen proximal femurs were evaluated (Table 3-1). The specimens were specified to have no reported history of musculoskeletal or metabolic disease. The bones, which had their soft tissues removed, were scanned using dual energy X-ray absorptiometry (DXA) (QDR 4500W, Hologic, Bedford, MA) and high-resolution (41 ?m) peripheral quantitative computed tomography (XtremeCT, SCANCO Medical AG, Br?ttisellen, Switzerland). Table 3-1: General data for femoral specimens Specimen Age/Gender (years/-) Height/Weight (m/kg) Femoral Neck aBMD (g/cm2) Total aBMD (g/cm2) T-score WHO classification H1167L 50/F 1.68/54 0.658 0.841 -0.8 Normal H1168R 73/M 1.70/77 0.650 0.813 -1.1 Osteopenic H1365R 71/F 1.52/45 0.504 0.638 -2.5 Osteoporotic H1366R 73/F NA/45 0.740 0.839 -0.8 Normal H1368R 70/F 1.65/57 0.500 0.606 -2.8 Osteoporotic H1369L 69/F 1.73/68 0.498 0.601 -2.8 Osteoporotic H1372R 79/F 1.68/113 0.638 0.793 -1.2 Osteopenic H1373R 76/F 1.52/63 0.773 0.833 -0.9 Normal H1374R 78/F 1.45/66 0.486 0.621 -2.6 Osteoporotic H1375L 83/F 1.52/90 0.633 0.782 -1.3 Osteopenic H1376L 79/F 1.63/64 0.552 0.683 -2.1 Osteopenic H1377R 80/F 1.65/64 0.503 0.627 -2.6 Osteoporotic H1380R 71/F NA/NA 0.695 0.834 -0.9 Normal H1381R 92/F 1.63/51 0.392 0.545 -3.3 Osteoporotic H1382L 96/F 1.52/61 0.441 0.546 -3.2 Osteoporotic Each bone was distally potted in a polyvinyl chloride (PVC) tube using PMMA. The potted bone and tube were then locked in an aluminum cylinder, part of a positioning fixture ? the same used during quasi-static testing ? that held the bone in a position that mimics the loading during a sideways fall [38], [39]. The cylinder and its positioning fixture were designed to only allow rotation around a distal quasi-anterior-posterior axis. After undergoing mechanical non-destructive testing for our quasi-static investigation (Chapter 2), each bone was destructively tested in a custom drop tower fall simulator [41]. This simulator, pictured in Figure 3-1, worked by guiding the drop of a mass on the greater trochanter (GT). The femoral head was constrained medially, but allowed to displace in the 3 ? Dynamic nonlinear simulation 36 sagittal plane using a steel plate supported on free-sliding bearing plates. The GT was loaded in a lateral to medial direction. Forces at both the head and the GT were distributed over the surfaces using PMMA pads, approximately 4 cm in diameter, to prevent local crushing. The pads were moulded to both bone and plate surfaces but not glued to them. Several pieces initially rested over the bone and represented different structures at play during a sideways fall. First, a PMMA pad lay directly over the GT, as described above. A 19 mm-thick closed-cell polyethylene foam (Evazote, Zotefoams plc, Croydon, England) rested over the PMMA and represented soft tissues. The foam supported a 51 x 51 x 19 mm aluminum impactor block and a six-axis load cell (Denton 4366J, Humanetics, Plymouth, MI) above it. Finally, a coil spring with a stiffness of 50 N/mm lay above the load cell and simulated the compliance of the pelvic bone during a fall [61]. The spring was directly impacted and compressed by the falling mass, in effect distributing the impact force over time. The 32 kg mass, selected based on the work of Robinovitch [104], was dropped from a height that, according to previously recorded calibration measurements, would produce an impact speed of 3.0 m/s. This speed was chosen from data presented by Feldman and Robinovitch [39]. Figure 3-1: Overview (left) and detail (right) of fall simulator and test setup. The white straps shown in the overview picture, meant to prevent prolonged foam compression, were removed before testing. 3 ? Dynamic nonlinear simulation 37 Two falling masses, adding up to 1.98 kg and representing the mass of a portion of the pelvic bone, were allowed to bypass the pelvis spring. These bypassing masses created an initial force spike seen in previous studies of sideways falls [105], [106]. After placing the bone in its positioning fixture, an Optotrak Certus motion capture system (Northern Digital Inc, Waterloo, ON, Canada, accuracy < 0.1 mm) was used to record the 3D locations of a number of points on the surface of the bone and several areas of the positioning fixture. These digitized data would later help to accurately reproduce loading directions and constraint locations on the computer models (see section 2.2.4). The six-axis load cell placed between the pelvis spring and the soft tissue foam measured forces during the impact event. These forces were recorded at 20 kHz using LabView (v2010 SP1, National Instruments, Austin, TX) and subsequently filtered using a Butterworth filter in forward and reverse directions with a 500 Hz cut-off. High-speed video of the impact event was recorded using four cameras at frequencies between 6 000 and 10 000 fps. Two Phantom V12 cameras (Vision Research, Wayne, NJ) observed the superior anterior femoral neck, one Phantom V9 camera observed the posterior femur, and one final Phantom V9 camera observed the impact surface at the GT from a superior perspective. This final camera observed a frame that included the PMMA-GT interface and the aluminum impactor block, allowing the subsequent digital tracking of displacement over time. The impactor block carried a high-contrast tracking marker to simplify this task. 3.2.2 Geometry and material mapping The high-resolution peripheral QCT data was calibrated using a standard hydroxyapatite (HA) phantom and filtered using a morphological filter to isolate the bone from intertrabecular pores. This step was necessary to avoid artificially stiffening the model by using a bone modulus-density relationship on the noisy data present in marrow-filled pores. The images were then reduced from their original 41 ?m resolution to clinical resolution of 0.615 mm in Matlab (v R2012b, The Mathworks, Natick, MA, US), ensuring that the average intensity of each voxel was maintained. The reduced images were segmented using ITK-SNAP (v2.2.0) [84], with the resulting data converted to a NURBS-based solid in SolidWorks (2011, Dassault Systems, Concord, MA). One of the bones was selected arbitrarily to generate a series of 2nd order tetrahedral meshes from the NURBS data using Ansys Workbench (v.14.0, Altair Engineering Inc., Troy, MI, USA). These\t ?meshes\t ?were\t ?tested\t ?with\t ?the\t ?FE\t ?analysis\t ?described\t ?below\t ?to\t ?determine\t ?the\t ?model?s\t ?3 ? Dynamic nonlinear simulation 38 sensitivity to mesh size. The test, which compared resultant forces and stiffnesses, concluded that a model with average element edge length around 3.0 mm yields similar results to those from denser meshes with a reasonable computation time. All bones were then meshed using similarly-sized elements. Summary mesh data is presented on Table 3-2. Table 3-2: Mesh characteristics for each bone 3.2.3 Material properties The voxel values from the reduced CT images were first corrected for partial volume artifacts. Surface voxels, which are at the boundary between bone tissue and the surrounding environment, were likely to suffer from partial volume artifacts: reduced properties due to the voxel volume being filled only partially with bone. To correct this, the outermost layer of bone voxels ? calculated from the mesh surface ? was replaced with the values from the next inner layer where these were higher.\t ?This\t ?technique,\t ?which\t ?our\t ?group\t ?refers\t ?to\t ?as\t ??peeling?, in effect extruded the last layer of voxels that contain only bone. The voxel values from the corrected CT images, representing hydroxyapatite (HA) content, were converted to ash density values by the following relationship: ???? = ????? + 0.09?/1.14 (3-1) where ?ash represents the ash density of the bone voxel in g/cm3, and ?QCT represents the corrected QCT voxel intensity in g HA/cm3. Apparent density was then calculated from the resulting voxel ash densities using the relationship Bone Volume (cm3) Elements Nodes H1167L 130.0 41 682 62 095 H1168R 169.2 53 359 78 626 H1365R 105.8 34 375 51 455 H1366R 144.0 45 457 67 513 H1368R 147.4 46 658 69 113 H1369L 171.9 53 373 78 676 H1372R 121.3 38 901 58 012 H1373R 154.4 48 518 71 820 H1374R 126.4 40 909 61 016 H1375L 132.6 41 746 62 037 H1376L 119.4 38 405 57 264 H1377R 142.8 74 000 105 546 H1380R 125.1 40 935 60 987 H1381R 162.8 51 548 76 148 H1382L 147.2 46 206 68 419 3 ? Dynamic nonlinear simulation 39 ????????= 0.60 (3-2) Both of these density relationships are based on the investigation of Schileo et al [85]. Stiffness was calculated for each CT voxel according to its apparent density, based on research performed by Morgan [72]: ? = 6 ?850?????.?? (3-3) where ?app is apparent density in g/cm3 and E is\t ?the\t ?material\t ?Young?s\t ?modulus\t ?in\t ?MPa. Despite, the good results we had using the Wille relationship in our quasi-static study (chapter 2), we chose to use the Morgan relationship in the dynamic models for its superior strain accuracy. The resulting\t ?Young?s\t ?modulus\t ?was used as the base parameter for the nonlinear material models that\t ?were\t ?applied\t ?on\t ?the\t ?model?s\t ?elements. Having\t ?corrected\t ?Young?s\t ?moduli\t ?assigned\t ?to\t ?each\t ?CT\t ?voxel,\t ?the\t ?stiffness\t ?at\t ?each\t ?node\t ?in\t ?the\t ?mesh was then calculated by 3D linear interpolation. These nodal stiffnesses were then distributed to nearby elements to obtain more traditional element stiffnesses. The distribution algorithm used weighted averages with nodal weight factors that were inversely proportional to the distance from node to element centroid. Since nonlinear properties were required, a piecewise linear plasticity model (#24 in LS-DYNA) was used for all elements in the femur. This model accepts a stress-strain curve defined by a series of points and it can also model strain rate-based hardening. Stress-strain\t ?curves\t ?were\t ?created\t ?based\t ?on\t ?the\t ?Young?s\t ?modulus\t ?calculated\t ?from\t ?Morgan?s\t ?modulus-density relationship. All materials had tension-compression symmetry. Two hundred materials were created between the softest and stiffest elements, assigning to each element the material with the closest modulus to its own. The initial (low-strain) response of each material was based on four points: the origin (0,0), a proportionality limit (?p,?p), a yield point (?y,?y), and an ultimate point (?ult,?ult) (Figure 3-2). ?y and ?u were set to 1.04% and 2% strains, respectively, independently from tissue density [107]. ?y was calculated using a 0.2% offset rule, as follows: ?? = ? ? (?? ? 0.2%) (3-4) where E is\t ?the\t ?Young?s\t ?modulus\t ?calculated\t ?from\t ?equation (3-3). The proportionality limit was set at 80% of ?y, while the ultimate stress was set to 110% of ?y [108]. The behaviour of bone at high compressive strains is not well understood, as it has received limited attention in previously published studies. However, preliminary testing suggested that some elements could have to withstand strains well beyond ultimate. To account for this, the 3 ? Dynamic nonlinear simulation 40 material stress-strain curve was extended to an extreme range that includes bone softening and densification. The curves presented by Hayes and Carter [109] were used as basis for the shape and strain values of the final material models, however, the stresses shown in these curves do not fit the Morgan model that was used for the low-strain part of the material, so they were scaled to Morgan levels to allow a qualitative adaptation of the behaviour, as described next. The softening and densification behaviour (Figure 3-2) was based on three points beyond the ultimate strain (?u). First, the stress-strain curve decreased linearly from the ultimate point (?u, ?u) to the first minimum (?min1, ?min) then remained at a constant stress until the second minimum (?min2,?min) and finally grew as a parabola until the maximum (?max, ?max). The position of these points was calculated using bone volume fraction (BV/TV) as a parameter, assuming a BV/TV of 1.8 g/cm3 corresponds to BV/TV = 1. The equations relating BV/TV to stresses and strains are below: ????? = c ? ?1 ?????? + ???? ????? = ??????? +(????/??)???? ? + ???? ???? = ?1 ? (1 ? BV/TV)??? ? ???? ???? = 1 ? BV/TV + ???? (3-5) Figure 3-2: Overview of nonlinear material properties 3 ? Dynamic nonlinear simulation 41 The parameters a, b, c, and d control the shape of the softening and densification curve and were set to emulate the material curves by Hayes and Carter mentioned above. The value of ?min was limited to a minimum of 0.8 ?ult. The value of ?max was set to the ultimate stress of dense bone, calculated as described above. Finally, strain rate-based hardening was applied as a factor modifying the stress values of each point calculated above. The factor was calculated using the following relationship: ? = ? ???.?????.?? (3-6) where R is the strain rate factor and ?? is the current strain rate. This equation is based on the work of Carter and Hayes [101], using the strain rate (1/s) that was used in the testing protocol of Morgan et al. [72]. This relationship was applied on the FE model via the strain rate-based hardening option in the piecewise linear plasticity material model. Ten strain rate-based material curves, dividing the range between ?? = 1?10??/? ? and ?? = 1/? using logarithmically spaced values, were fed into the model for each of the 200 materials, allowing the program to interpolate between them. 3.2.4 Boundary conditions The femur was constrained at two points: the shaft and the femoral head. It was also loaded at the greater trochanter using solid models for the PMMA, foam, and impactor block. All boundary conditions were generated using points that were digitized on the experimental setup, as described in section 2.2.4. Figure 3-3 shows an overview of the applied boundary conditions. 3 ? Dynamic nonlinear simulation 42 The distal (inferior) shaft was constrained in such a way that the bone could only rotate around a distal hinge on the frontal (coronal) plane. This was accomplished by connecting the distal surface of the shaft to a hinge node via beam elements in a similar fashion to what was done in the quasi-static model, as described in section 2.2.4. The beam elements were modelled to emulate the axial and bending stiffness as well as the mass of the PMMA-filled aluminum cylinder where the real femoral shaft was potted (Figure 3-1). To simulate reality as closely as possible, three types of beams were used in series: The most proximal beam type used dense bone properties and circular tube geometry to simulate the distal portion of the femoral shaft that was outside the CT scanner range. The middle beam type was given composite properties corresponding to a PMMA-filled aluminum cylinder, emulating the potted portion of the cylinder. The most distal beam type used the properties of a hollow aluminum cylinder to simulate the portion beyond the potting. This last beam type connected to the hinge node and continued a short distance beyond to represent the short section of aluminum cylinder that extended beyond the hinge. All dimensions accurately represented the mechanical parts used in the fall simulator. Load distribution at the femoral head and greater trochanter was achieved using solid models representing the PMMA support and loading pads. The pads were modelled as cylinders Figure 3-3: FE boundary conditions Ant Lat Inf 3 ? Dynamic nonlinear simulation 43 moulded to the bone surface and meshed with the rest of the bone model, ensuring that the elements in the bone and pads were similarly sized and had matching meshes at their contact surfaces.\t ?The\t ?material\t ?model\t ?representing\t ?PMMA\t ?was\t ?linear\t ?elastic\t ?with\t ?a\t ?Young?s\t ?modulus\t ?of\t ?3.1 GPa. The interaction between femoral head and PMMA was simulated using a sliding frictionless contact defined through the Automatic Single Surface Contact algorithm in LS-DYNA. The use of support pads and sliding contact allowed the bone to displace and rotate freely at its support and loading points, closely emulating the true experimental constraints. A solid block model was used to represent the soft-tissue foam above the loading pad at the greater trochanter. This block was modelled with regular hexahedral elements and a low-density foam material model. The foam material model used a stress-strain curve derived from force-displacement data from the foam used in the fall simulator. Force data came from the six-axis load cell and foam deformation data from the difference in tracked displacements between the impactor block and the PMMA-femur interface in the high-speed video. Using these high-speed measurements meant that the foam behaviour did not need further adjustment due to loading rate. An average foam material curve was used for all models. While the support pad at the femoral head simply had its bottom (medial) nodes constrained in the medial-lateral axis, load application at the greater trochanter required some additional complexity. A rigid plate was created above the soft-tissue foam, and this plate was assigned a velocity-time curve calculated from video tracking of the impactor block (Figure 3-4). This tracked displacement was corrected for the compliance of the drop tower according to a simple dynamic lumped-mass model that used the measured stiffness of the support structure. No corrections for the compliance of cartilage were made, assuming that the thin cartilage displacements are very small at the high loading rates of this experiment [110]. Velocity was used instead of displacement to avoid the creation of infinite accelerations in the dynamic model. Both the rigid plate and the topmost nodes of the foam were fixed in the sagittal plane, allowing only medial-lateral displacements. 3 ? Dynamic nonlinear simulation 44 3.2.5 Comparison criteria The comparison of FE to experimental results was done according to initial overall stiffness, peak force, and energy absorbed, as explained below. Additionally, a qualitative comparison of the FE areas of high strain against fracture locations was performed. Experimental forces were measured using the six-axis load cell, filtered in forward and reverse directions using a 4th-order Butterworth filter. Experimental displacements were obtained from the high-speed video (Figure 3-4) using automated tracking of the bone-PMMA interface at the greater trochanter (TEMA Automotive v3.0, Image Systems, North Hollywood,CA). The tracked GT displacement was corrected for machine compliance by subtracting the support deformation estimated using a single-degree-of-freedom model; the machine stiffness used for this model was measured independently. FE forces were measured as the contact forces between the rigid loading plate and the soft-tissue foam. FE displacements were calculated as the average displacements of all surface nodes in the greater trochanter that were in contact with the PMMA loading pad. Figure 3-4: Sample frames from displacement tracking video showing the superior surface of the neck. The arrows mark the approximate tracked position of the impactor block (red, pointing left) and PMMA-femur interface (green, pointing right) Lat Pos 3 ? Dynamic nonlinear simulation 45 Overall stiffness was calculated for both ex-vivo and FE models from the initial linear portion of each\t ?bone?s\t ?force-displacement curve. The linear range was selected by visual inspection of each curve. Force was measured by the six-axis load cell and GT displacement was measured Peak force was measured experimentally as the first large peak in the load-cell data. The first peak was used ? even when this was not the largest peak observed before collapse ? because any bone damage caused during the first peak would be present in later peaks. Given the limited ability of the FE model to simulate bone damage, it was reasoned that using the first peak to avoid damaged bone would be the most reasonable course. On the FE side, however, there was generally no peak force, since the model was not programmed to detect fracture or damage beyond some level of post-yield softening. Thus, the FE peak force was compared at two points: at the same time as the experimental force peak (time matching) and at the same GT displacement as the experimental force-peak (displacement matching). These two points were generally not the same, since it was the displacement of the impactor and not of the GT that was applied as a boundary condition; the GT displacement in the FE model was thus a resulting value that depended on the stiffness of the bone and foam. The energy absorbed by the bone was calculated as the area under the force-GT displacement curves for both experimental and FE results. Similarly to the force calculation, FE energy results were calculated up to the experimental peak using both time-matching and displacement-matching. Finally, strain comparison was achieved by comparing the plots of principal compressive and tensile strains in the neck and intertrochanteric regions from the FE model to the fracture locations seen in high-speed videos. 3 ? Dynamic nonlinear simulation 46 3.2.6 Parameters tested for sensitivity The model presented above was tested for sensitivity against several of the parameters and features it includes. The form of the model already described in the previous sections was treated\t ?for\t ?comparison\t ?purposes\t ?as\t ?a\t ??standard\t ?model.?\t ?Four bones (H1365R, H1369L, H1372R, and H1373R) that demonstrated varied behaviours in their results from the standard model (see section 3.3.1) were used in the sensitivity analysis. The use of four bones was intended to account for inter-specimen variability. Unless otherwise noted, only one aspect of the standard model described in the previous sections was changed for each test. The sensitivity of the following features was tested: Foam behaviour: The standard FE model?s foam material used the mean force-displacement behaviour of the foam pads measured during the fall-simulator tests. To test the sensitivity to this behaviour, two more force-displacement relationships were used for the foam. These relationships\t ?are\t ?described\t ?as\t ??stiff?\t ?and\t ??compliant?\t ?according\t ?to\t ?their\t ?respective\t ?behaviour\t ?at\t ?high strains. Figure 3-6 shows a comparison of the stiff, standard, and compliant foam force-displacement relationships. One further model was created eliminating the foam completely. In this last model, the video-tracked displacement of the greater trochanter ? not the impact block ? was applied to the rigid plate. The PMMA pad over the greater trochanter was still present to allow free sliding and rotation, but it was made fully rigid. Figure 3-5: Energy comparison from FE data using time and displacement matching. Experimental results are in black, FE results in blue 3 ? Dynamic nonlinear simulation 47 Nonlinear material behaviour of bone: Two aspects of the material nonlinearity were investigated. The first aspect, strain rate dependence, was tested using a rate independent model. This model used a strain rate factor R ? see equation (3-6) ? equal to 1 at all strain rates. The second aspect, post-yield behaviour, was tested by eliminating the softening and densification portions from all material curves, leaving perfectly plastic deformation beyond ?ult ? see Figure 3-2. One further model used material properties with disabled strain rate dependence and disabled post-yield behaviour. Material modulus-density relationship: As shown in chapter 2, the Morgan relationship used in the standard model is relatively stiff compared to other published relationships. To test the influence of this parameter, a model that used the more compliant Keyak-Keller relationship [12] was tested. 3.3 Results 3.3.1 Ex-vivo fall simulations The fifteen femurs were successfully tested and the measurements recorded. A thorough analysis can be found in Gilchrist?s\t ?thesis and results papers [77], [111]. The time to peak force Figure 3-6: Comparison between stiff, mean, and compliant foam behaviour 3 ? Dynamic nonlinear simulation 48 ranged from 13.7 to 19.0 ms (x?: 16.2 ms,\t ??: 1.19 ms), while the total impact duration ranged from 20.7 to 79.3 ms (x?: 38.7 ms,\t ??: 19.5 ms). These numbers demonstrate shorter impacts than those recommended by Robinovitch et al [61]. Each of the tested bones was examined by an orthopaedic surgeon who deemed their fractures clinically relevant. Force-time results for all fifteen tests showed similar features, including one or more small and short peaks ? likely due to the impact of the bypass masses ? followed by a steep rise in force up to the peak value. Each specimen roughly followed one of three behaviours (Figure 3-7): 1) Four bones (H1365R, H1366R, H1369L, and H1382L) collapsed soon after the force peak, immediately showing a rapid decrease in force. 2) For nine bones (H1167L, H1168R, H1368R, H1372R, H1374R, H1375L, H1376L, H1377R, and H1381R), there was an elongated period following the first force peak before finally collapsing with a rapid decrease in force. 3) The final two bones (H1373R and H1380R) completed the drop without total collapse, thus showing a very long period of sustained, somewhat oscillating force beyond the first peak. Both bones exhibited trochanteric crushing, with H1380R suffering a much higher degree of damage than H1373R. In fact, this last bone appeared at first glance to have withstood the drop unscathed; it was only during examination of the high-speed video that the trochanteric damaged was discovered. Figure 3-7: Samples of three behaviours observed experimentally. The dashed lines represent the time at peak force 3 ? Dynamic nonlinear simulation 49 3.3.2 Force, displacement, and time relationships This section describes the relationships between three variables: force, GT displacement, and time. Force and displacement were measured experimentally and calculated from the FE model as described in section 3.2.5. Time was used as the common frame of reference between the ex-vivo and FE models, since this ensured the comparison was done with equal impactor displacement on both models. It is important to reiterate here that the GT displacement was not an input in the standard FE model. Instead, the velocity of the impactor block measured from high-speed video (Figure 3-4) was used to load the rigid impact plate in the FE model. The GT displacement was a result value compared between experimental and FEA. It was not only controlled by the movement of the impact plate, but also by the stiffness of the bone, foam, and PMMA. The force-time plots (Figure 3-8) show similar relationships between the experimental and FE data, with the force rising at similar rates and at somewhat similar times. There is a clear time delay, however, in the FE model force rise for most bones, compared to their experimental curves. The experimental curves show a short early force peak attributed to the impact of the bypass masses; none of the FE curves displays this behaviour, however. It is also apparent that the FE curves do not all show a force peak, but instead almost all continue growing in force at a decreased rate. This monotonic behaviour could be attributed to the limited simulation of element damage in the FE model. As elements are strained, they can only develop some degree of softening; they never break or collapse, even if when subjected to extreme tensile strains. The force-GT displacement plots (Figure 3-9) are shown as an attempt to separate the response of the bone from that of the compliant structures (foam and PMMA). The force is, however, the same force described for the force-time plots above, so it does include any inertial forces that could potentially develop above the bone. Both experimental and FE results generally show a period of relative linearity followed by gradual yielding. After yielding, most experimental curves peak, showing reduced reaction force for additional deformation. Most FE curves, on the other hand, do not peak but only flatten into what appears to be a force asymptote, with constant force for increased deformation. The initial behaviour of H1167L and H1168R is noticeably different from that of the other bones (Figure 3-9). H1167L shows an early force rise under little deformation and also a very short post-yield peak. H1168R shows a decrease in deformation beyond its force peak and signs of hysteresis on the rebound. The corresponding high-speed videos show that the bones resisted 3 ? Dynamic nonlinear simulation 50 the first compressive wave largely unscathed, failing only after a slight but visible rebounding of the impactor block. This small bounce is likely the reason behind the force behaviour, as the displayed curves mainly show the effects of the first wave. The FE curves for both of these bones show similar features, supporting the idea that it is the rebound of the impact block ? which was directly used as input for the FE model ? that leads to this behaviour. Displacement-time plots (Figure 3-10) are included to give an idea of how the previous two sets of plots are related. The dashed line included shows the displacement of the impactor block, which is by design equal for experimental and FE models. Note that for most of the bones GT displacement was underestimated when compared at equal times. 3 ? Dynamic nonlinear simulation 51 Figure 3-8: Force-time comparison between experimental and FE results 3 ? Dynamic nonlinear simulation 52 Figure 3-9: Force-displacement comparison between experimental and FE results 3 ? Dynamic nonlinear simulation 53 Figure 3-10: Displacement-time comparison between experimental and FE results 3 ? Dynamic nonlinear simulation 54 3.3.3 Stiffness, peak force, and energy regressions Figure 3-11 shows the comparison between ex-vivo and FE stiffness. These values were derived from the most linear portion of the respective force-displacement curves, as shown in Figure 3-9. The regression curve shows relatively good estimation of stiffness for the lower end of the scale, but clear underestimation for the stiffest bones. The correlation coefficient, R2, and root mean square error are both moderate values. Figure 3-12 shows the comparison between experimental peak forces and the FE forces corresponding to the same times and displacements. While the time-based comparison resulted in somewhat better slope and correlation, both approaches led to similar mean errors. Generally, the current model appears to be a rather poor force predictor. H1167LH1168RH1365R H1366RH1368R H1369LH1372RH1373RH1374RH1375LH1376LH1377RH1380RH1381RH1382Ly = 0.2024x + 1.1268R? = 0.320RMSE = 22.9%p = 0.025801234560 2 4 6FE Stiffness, N/mmExperimental Stiffness, N/mmTrendline 1:1 MatchFigure 3-11: Stiffness results regression 3 ? Dynamic nonlinear simulation 55 Figure 3-13 shows the comparison between experimental and FE strain energy absorbed by the bone. As briefly described in section 3.2.5, both energy values were calculated for each bone by integrating the measured force with respect to GT displacement. The integration was done from the start of the impact event to either the time (time-matching) or GT displacement (displacement-matching) corresponding to the experimental peak force. The results from both methods show a moderate predictive ability of the model, with slopes closer to unity and slightly higher correlation coefficients than the stiffness and force regressions. It should be noted, however, that the correlation values, most notably with displacement matching, seem to be enhanced by the results from two bones ? H1375L and H1376L ? with the other results demonstrating a more scattered and less linear pattern. Mean error values are moderate, similarly to the other regressions?. H1167LH1168RH1365R H1366RH1368R H1369LH1372RH1373RH1374RH1375LH1376LH1377RH1380RH1381RH1382Ly = 0.4125x + 0.8173R? = 0.211RMSE = 23.43%p = 0.081900.511.522.533.540 1 2 3 4FE Force, kNExperimental Force, kNH1167LH1168RH1365R H1366RH1368RH1369LH1372RH1374RH1375LH1376LH1377RH1380RH1381RH1382Ly = 0.4512x + 1.1971R? = 0.1185RMSE = 23.47%p = 0.227100.511.522.533.540 1 2 3 4FE Force, kNExperimental Force, kNTrendline 1:1 MatchFigure 3-12: Peak force results regression based on time-matching (left) and displacement matching (right) 3 ? Dynamic nonlinear simulation 56 3.3.4 High strains and fracture locations Table 3-3 compares the experimental locations of initial fracture for each bone, manually identified in a qualitative fashion by Gilchrist from the high-speed videos, with the minimum (compressive) principal surface strains results from the FE model. Only fractures that developed within 5 ms of the first fracture are shown on the figure. The FE model is shown at the same GT displacement as the experimental peak force. It was determined with a simple qualitative visual check that the locations matched well for seven of the bones (H1167L, H1168R, H1369L, H1374R, H1375L, H1377R, and H1381R), matched partially for seven more (H1365R, H1366R, H1368R, H1372R, H1376L, H1380R, and H1382L), and failed to match for only one bone (H1373R). This last bone was the least damaged out of all tested bones, suffering only very minor damage that was only visible during examination of the high-speed video and not during examination of the tested specimen. Given the variability of resulting forces and deformations shown in the previous sections, it is not surprising to find that some bones have higher strain magnitudes than others at the GT displacement corresponding to peak force. Most bones in the figure are shown with a colour scale that assigns yellow to strains of 2% (corresponding to the ultimate point in the material model Figure 3-2); however, some bone models with large portions of their surface showing H1167LH1168RH1365RH1366RH1368RH1369L H1372RH1373RH1374RH1375LH1376LH1377RH1380R H1381RH1382Ly = 0.6832x - 0.2479R? = 0.3813RMSE = 27.62%p = 0.012701234567890 2 4 6 8FE Energy, JExperimental Energy, JH1167LH1168RH1365RH1366RH1368RH1369LH1372RH1373R H1374RH1375LH1376LH1377RH1380RH1381RH1382Ly = 1.354x - 1.2976R? = 0.4733RMSE = 28.20%p = 0.003801234567890 2 4 6 8FE Energy, JExperimental Energy, JTrendline 1:1 MatchFigure 3-13: Energy of deformation regression based on time-matching (left) and displacement matching (right) 3 ? Dynamic nonlinear simulation 57 high strains were coloured with a scale up to 5% strain for clarity. Those bones are marked in the figure with an asterisk. Minimum (compressive) principal strains are shown, as it has been shown that fractures typically initiate in areas loaded in compression during a sideways fall [26]. For completeness, however, a similar figure with maximum (tensile) principal strains is included as an appendix. Table 3-3: Comparison of initial experimental fracture locations ? as marked by Gilchrist ? against compressive principal FE strain patterns (continued) Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1167L H1168R H1365R 3 ? Dynamic nonlinear simulation 58 Table 3-3: Comparison of initial experimental fracture locations ? as marked by Gilchrist ? against compressive principal FE strain patterns (continued) Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1366R H1368R H1369L H1372R 3 ? Dynamic nonlinear simulation 59 Table 3-3: Comparison of initial experimental fracture locations ? as marked by Gilchrist ? against compressive principal FE strain patterns (continued) Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1373R H1374R H1375L* H1376L* 3 ? Dynamic nonlinear simulation 60 Table 3-3: Comparison of initial experimental fracture locations ? as marked by Gilchrist ? against compressive principal FE strain patterns (continued) Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1377R H1380R H1381R* 3 ? Dynamic nonlinear simulation 61 Table 3-3: Comparison of initial experimental fracture locations ? as marked by Gilchrist ? against compressive principal FE strain patterns (continued) Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1382L Standard FE colour scale: FE colour scale for bones marked with an asterisk (*): 3.3.5 Sensitivity analyses Figure 3-14 shows the force-displacement results using the different foam models. In general, the\t ?stiff\t ?foam\t ?produced\t ?very\t ?similar\t ?results\t ?to\t ?the\t ?standard\t ?model?s\t ?mean\t ?foam,\t ?while\t ?the\t ?compliant foam led to visibly stiffer results overall. The model that uses no foam ? which applies the\t ?tracked\t ?displacement\t ?of\t ?the\t ?GT\t ?instead\t ?of\t ?the\t ?impactor\t ?block?s ? shows results that are fairly similar to those from the stiff and mean foam models. This result hints at the possibility of eliminating the foam from the standard dynamic model in the future, thus eliminating the variability due to it. 3 ? Dynamic nonlinear simulation 62 Figure 3-15 compares the sensitivity of results using different nonlinear properties for the bone material. All the curves behave similarly under low levels of displacement; this is expected, as linear elastic behaviour tends to control this portion of the curves. On the opposite end, the curves diverge at high displacements in a somewhat predictable fashion, with the plastic curve ? which does not suffer from post-yield softening ? showing the highest forces and the rate independent curve ? which does not benefit from strain-rate hardening ? showing the lowest. The standard and combination (plastic & rate independent) models combine a hardening effect with a softening effect, resulting in forces in between the two extremes. It can be seen, however, that rate dependence is the more powerful factor, dividing the curves earlier and further compared to the difference created by the two post-yield behaviours. Figure 3-14: Comparison of sensitivity to foam properties 3 ? Dynamic nonlinear simulation 63 Figure 3-16 compares the force-displacement curves from two modulus-density relationships: the Morgan relationship used in the standard simulation and the Keyak-Keller relationship described in chapter 2. As expected, the Keyak relationship results in visibly decreased stiffness and lower forces overall. The modulus-density relationship appears to play an important role in determining the magnitudes of forces, stiffness, and energy Figure 3-15: Comparison of sensitivity to nonlinear material behaviour 3 ? Dynamic nonlinear simulation 64 3.4 Discussion The current investigation presents, to the best of our knowledge, the first published study using explicit dynamic FE modelling to simulate the behaviour of a femur under the high-speed loading associated with a sideways fall. This type of FE model was chosen because it allows a straightforward implementation of nonlinear properties of bone, in particular its strain-rate hardening, and because it is capable of simulating the dynamics involved in an impact. Our group is currently investigating the use of deformation energy for fracture prediction (see Appendix). It is important for this investigation that a dynamic FE model is capable of predicting the deformation energy of a bone during impact. While the current model does not yet simulate local damage or cracking reliably enough to predict energy at the final collapse, the results shown for energy at first fracture are an encouraging first step. The results suggest that deformation energy can be modelled using dynamic FE simulations. Even though peak forces derived from the model do not match the experimental values well, the energy of deformation ? which is derived in part from force values ? does show a moderate match. Figure 3-16: Comparison of results using two modulus-density relationships 3 ? Dynamic nonlinear simulation 65 Stiffness results on this dynamic model show an important pattern when compared to the quasi-static results presented in section 2.3: Even though the dynamic model used the Morgan modulus-density relationship, which was shown to overestimate stiffness in the quasi-static models (Figure 2-6), and that the model was further stiffened through the use of a strain rate factor (equation (3-6)), the dynamic results show a general underestimation of stiffness for the stiffest bones (Figure 3-11). This means that our specimens, or at least the stiffest among them, demonstrated a higher increase in stiffness going from quasi-static to dynamic loading than our FE models anticipated. While the dynamic model is complex, and there could be a number of confounding factors in it, this pattern potentially points toward a larger influence of dynamics and strain-rate hardening on bone mechanics than previously thought. It certainly points toward the need for further investigation of impact loading on the femur. Out of the results presented, perhaps the one with the closest concordance is the comparison between high principal strains and fracture locations. Despite the differences in force and stiffness magnitudes on each of the bones, the high-strain locations predicted where the initial surface fracture would develop in most of the bones. These results point to the ability of the model to detect the most vulnerable portions of bone given a particular loading scenario. A quantifiable test of this ability would compare the surface strains seen in the ex-vivo tests to those from the FE model. Our group is planning to perform this test in the near future using stereo digital image correlation (DIC) to measure strains from the images recorded by two of the high-speed cameras. The sensitivity comparisons reported here are meant to provide guidance on which aspects of the FE model deserve the most attention for future improvements. The results show, first of all, that the foam behaviour is only an important factor if the foam itself is highly compliant. It is indeed expected that having two structures with widely differing stiffnesses loaded in series, the compliant structure will largely dictate the behaviour of the system. Increasing the foam stiffness brings the two stiffnesses closer together, allowing the bone to have a larger role in the overall behaviour and diminishing the effect of further foam stiffness increases. The test that completely eliminated the foam yielded similar results to the standard model while eliminating an important source of variability; it is thus recommended that future models also eliminate the foam and rely on the GT displacement tracked from high-speed videos. It is necessary here to temper the above recommendation by pointing out that it is not always beneficial to simplify the boundary conditions on this type of model. During development, we 3 ? Dynamic nonlinear simulation 66 tested the influence of other boundary conditions and found that including more detail ? in particular, adding the sliding PMMA support pad and the distal beam models ? generally improved the match between experimental and FE results. It is, therefore, important in our opinion to model the boundary conditions in as much detail as is reasonable and then test whether a simplification is still valid, as was done here. The modulus-density relationship demonstrated a sizeable effect at all points in the simulation, while strain-rate hardening and post-yield material properties have a more prominent effect after some bone damage has accumulated. It is thus important that modulus-density relationships are thoroughly investigated to find what relationship is most appropriate for an explicit-dynamic FE model. Comparing the effects of strain-rate hardening and post-yield material properties, it appears that it is the strain-rate effect that is most prominent in the current model. There are other models, aside from that presented by Carter and Hayes [101] and used for our standard model, that incorporate the effect of strain rate on cortical and trabecular bone properties [112]?[115]. Few studies, however, have focused an experimental investigation on the post-yield behaviour of bone reaching extreme strains greater than 10%. It seems unlikely, given the sensitivity results presented, that adjustments to these two aspects of bone material behaviour will account for the observed differences between ex-vivo and FE results. Similarly to the quasi-static analysis presented in Chapter 2, it is likely that the estimate of the femur deformation introduces a level of error in the results. The compliance of surrounding structures, like PMMA or foam, makes it challenging to determine precisely what portion of the overall displacement belongs to the femur. In the standard FE model, we included models for the foam and PMMA pads to account for their deformation; while in the sensitivity model without foam, we applied the video-tracked displacement of the greater trochanter directly. In the four bones where the sensitivity models were run, the results of the two methods were quite similar, which is an encouraging sign. However, both methods are very likely to carry significant errors that would affect our results. While it is clearly very challenging to accurately measure the femur deformation during the ex-vivo impact, it is highly recommended that this area of the empirical research is analysed and improved for future investigations. There are a few other aspects that were not explored in the current investigation, but which warrant further study. The first is the role of bone marrow in the dynamic behaviour of bone. Our FE models produced overall stiffnesses spanning a visibly narrower range than the ex-vivo 3 ? Dynamic nonlinear simulation 67 models?\t ?results. It seems likely that there is a mechanism that is not currently being modelled that could account for the additional stiffness variability shown in the experimental results. One possible source of this variability is the bone marrow, which, as a fluid, has mechanical properties that are dependent on loading rate, the level of confinement provided by the cortical bone,\t ?and\t ?the\t ?marrow?s\t ?ability\t ?to\t ?flow\t ?through\t ?this\t ?confinement\t ?[116]. Another aspect that was not included in this study and would likely have a significant effect on the results is compression/tension asymmetry. The FE models in the current investigation applied the same material properties in both tension and compression, despite the use of a softening and densification model that could only be applicable to compressive behaviour. Bone is highly asymmetric, particularly in its failure modes [117], [118]. Application of asymmetry could range from the simple use of element death once ultimate tensile strain is reached to a full implementation of differing material curves and strain-rate behaviour for tension and compression. Other aspects left unexplored in our research include the role of muscles, ligaments, and synovial fluid as sources of both confinement and additional loading on the bone; and the possibility that small features like vascular canals could create stress concentrations. In summary, this investigation demonstrated how explicit dynamic FE modelling can be used to simulate femurs under sideways fall impact loading. The presented data show the current abilities of this type of model and also where improvements are necessary. This thesis takes the first steps toward a novel methodology for using FE modelling to understand and simulate the behaviour of femurs during a sideways fall. 68 4 Multi-scale nested simulation 4.1 Introduction Finite element modelling of femurs relies on the appropriate simplification of a very complex structure. The actual proximal femur is largely composed of highly porous trabecular bone. Structural support for the bone is provided by a thin layer of cortical bone together with the thousands of thin struts and beams (i.e. the trabeculae) that form the solid phase of trabecular bone. This makes the femur a highly heterogeneous environment. To model the behaviour of the entire organ using the FE method, researchers must\t ?assume\t ?that\t ?each\t ?element?s\t ?volume\t ?can\t ?be\t ?approximated with a homogeneous material. There are two ways to do this: either to use elements small enough to represent the geometry of individual trabeculae, or to use larger elements that cover both trabeculae and pores and apply a material that approximates the homogenized properties of the porous volume. Given the small size of trabeculae, a model of the whole proximal femur using the former method requires an enormous number of elements, massive computing resources, specialized solvers, and relatively simple material properties; nonetheless, the technique has been used [76]. Unsurprisingly, the majority of research on FE models of the femur have followed the latter method [11], [54], [55]. However, high-resolution FE investigations modelling the trabecular structure have been shown to produce accurate results in models of small trabecular bone samples [119], [120]. High-resolution FE models are also likely to have a role in fracture prediction, as small stress concentrations could be generated at small scales in trabeculae and around cortical bone defects. In an effort to obtain the benefits of high-resolution FE modelling without incurring the extreme computational costs associated with them, our team is developing a methodology that allows the creation of high-resolution models of small bone pieces while applying whole-bone boundary conditions. The process is done in two stages: one solid model of the entire proximal femur using relatively large, internally homogenized elements; and one high-resolution porous model of a location of particular interest. By using the deformation results of the first stage as the boundary conditions of the second, a nested multi-scale model can be formed. The relatively small size of these models allows the application of nonlinear material properties and dynamic effects at the trabecular level without incurring the extreme computing costs of whole-bone high-resolution models. 4 ? Multi-scale nested simulation 69 On this chapter, the feasibility of this methodology is explored by comparing the results of a nested model to measurements of a femoral bone tested under simulated sideways fall conditions. The data here presented is currently in review for publication as part of a larger investigation [63]. 4.2 Methodology One fresh-frozen proximal femur was tested dynamically on an early version of the fall simulator described in section 3.2.1. High-resolution QCT images previously collected were reduced to clinical resolution to create an FE model of the entire proximal femur, here referred to as the organ-level model. Data collected from high-speed video taken during the ex-vivo fall simulation and from the organ-level FE model were used to identify three volumes of interest (VOIs) that would become the micro-level models. These smaller, more detailed models were created using geometry data from the full resolution CT images and boundary conditions derived from the organ-level displacement results. The organ-level model and all three micro-level models were processed using the explicit LS-DYNA solver in ANSYS Mechanical APDL (v14.0). 4.2.1 Ex-vivo and organ-level models The experimental setup used for this chapter\t ?formed\t ?an\t ?early\t ?part\t ?of\t ?Gilchrist?s\t ?investigation\t ?and was executed by Jason Chak. The fall simulator was similar to the one described in section 3.2.1, but with several important differences: x The drop mass was only 16.5 kg x There was no hip compliance spring, bypass masses, or soft tissue foam. The drop mass, which included both the six-axis load cell and an impactor plate, fell directly on the PMMA pad placed over the greater trochanter. x Instead of a PMMA pad, a tennis ball cut in half was used under the femoral head to distribute the reaction force. The rubber in the tennis ball formed a somewhat more compliant support than PMMA. Similarly, the organ level FE model mostly followed the description in sections 3.2.2 to 3.2.4, but with several differences: 4 ? Multi-scale nested simulation 70 x The frictionless support pad under the femoral head was made very stiff. The compliance of the tennis ball was instead taken into account by estimating its displacement from an isolated mechanical test and subtracting from the input displacement. x The frictionless pad at the greater trochanter was also made very stiff. The compliance of the PMMA pad present in the ex-vivo test did not require modelling, as the input displacement applied on the GT was taken directly from tracking of the PMMA-GT interface in the high-speed video. x Since there was no soft tissue foam in the ex-vivo model, there was also no foam block in the FE model. The rigid impact plate directly contacted the stiff pad at the GT. 4.2.2 Micro geometry and material properties Three VOIs, each forming a cube with 5 mm sides, were chosen as follows (Figure 4-1): x VOI A: Location of the first crack appearing in the femoral neck region, according to high-speed video. x VOI B: Location of a crack appearing on the anterior surface of the femur about 1 ms after the crack above, according to high-speed video. x VOI C: Location of highest surface strains on the neck region, according to the organ level FE model. 4 ? Multi-scale nested simulation 71 The VOIs were located in the full 41-?m resolution QCT images. The filtered CT values in these volumes were used to calculate\t ?Young?s\t ?moduli\t ?using\t ?Morgan?s\t ?modulus-density relationship (see section 3.2.3). Every voxel with a modulus of 10 MPa or higher was treated as bone. Unconnected groups of voxels were removed using a component labelling function in Matlab (v R2012b, The Mathworks, Natick, MA, US). Then, the largest group of bone voxels was directly turned into an FE model ? referred to as the micro-level model ? with cubic elements. Each element was directly assigned the modulus of the corresponding CT voxel. Post-yield behaviour and strain-rate dependence were derived in the same way as described in section 3.2.3 and applied using a piecewise linear plasticity model. 4.2.3 Boundary conditions The results from the organ-level model were used as input for the micro-level model. First, all of the organ-level elements covering each VOI were isolated, obtaining the resultant velocities of each node in those elements. The organ nodal velocities were treated as a field of values in order to interpolate velocities to micro-level nodes. Prescribed velocities were applied only on micro-level nodes located on the exterior surface of the bounding cube ? i.e. only nodes forming part of the cut faces. Figure 4-1: Location of VOIs on organ-level model 4 ? Multi-scale nested simulation 72 4.2.4 Comparison criteria The purpose of the organ-level model in this short investigation was to provide a platform for the application of boundary conditions on the micro-level models. For this reason, only the model?s\t ?deformations\t ?were of interest. Experimentally, Gilchrist calculated surface strains fields using digital image correlation (DIC) of the high-speed video. The DIC algorithm compared portions of each video frame to calculate the motion and deformation of the bone surface. This information could be used to estimate surface strains that could be compared to those from the organ-level FE model. Before the ex-vivo test, the femoral neck was covered with a layer of white paint and a black speckle pattern intended to enhance tracking for DIC. In this early version of the methodology only one high-speed camera was used, and thus only 2D DIC could be performed. For this reason, only regions approximately parallel to the video frame could be taken into account. Comparison of the micro-level models to experimental data is rather difficult, as there is no reasonable way to measure forces or stresses during the impact. Each VOI was instead compared to post-fracture video frames to determine whether the deformation of the small volumes reflected the large deformations that resulted in fracture of the actual bone. 4.3 Results 4.3.1 Ex-vivo fall simulation The high-speed video showed that the bone first fractured at the superior end of the greater trochanter 4.88 ms after contact. The main fracture leading to total collapse was first seen at the superior femoral neck 6.99 ms after contact. Other fractures developed at later times until catastrophic failure was seen at t = 8.21 ms. For modelling purposes, fracture time was defined at 6.99 ms. Calculating strains through DIC was somewhat difficult, as the impact caused saline solution ? which had been regularly sprayed on the specimen to maintain its moisture ? to splash between the bone and the camera. The occlusions created by the flying fluid affected the strain calculation in some areas of the neck. Nevertheless, surface strains were successfully calculated in most of the\t ?neck?s\t ?anterior\t ?surface,\t ?reaching\t ?principal\t ?compressive\t ?strains\t ?up\t ?to\t ?7.1%\t ?on\t ?the\t ?superior\t ?neck. 4 ? Multi-scale nested simulation 73 4.3.2 Organ-level model Figure 4-2 compares the surface compressive principal strains from the organ-level FE model to those obtained by DIC at three points in time. The FE strains are generally of higher magnitude than the experimental measurements. Their distribution, however, resembles that obtained from DIC, with areas of high strain on either side of the intertrochanteric line. The maximum principal compressive strain developed on the organ-level FE on the anterior neck area was 24.5%. 4.3.3 Micro-level models Given that applied velocities were directly imported from the organ-level model, which models the displacements of the entire bone, the micro-level models all showed a sizeable amount of whole body displacement. Only the differences in velocities between nodes produced deformation in these models. The large displacements developed on the organ-level model resulted in convergence difficulties on the micro-level models, so none of the three micro models was able to reach the fracture time Figure 4-2: Comparison of organ-level surface strains to DIC-based strains calculated by Gilchrist 4 ? Multi-scale nested simulation 74 of 6.99 ms. Figure 4-3 shows the micro-level results at three different time points. The rightmost figure\t ?for\t ?each\t ?VOI\t ?represents\t ?the\t ?model?s\t ?state\t ?at\t ?the\t ?maximum\t ?available time. The results from VOI A show a general direction of high-strain that is consistent with the fracture of the superior GT seen in the high-speed video at 4.88 ms. Similarly, VOI B showed a pattern of high strains coming from its superior-lateral edge, consistent with the main fracture line passing through it. The model for VOI B also revealed the presence of a 1-2 mm wide vascular hole. This hole had not been discovered in the specimen due to its small size and the layer of paint covering it. Thick walls surround the hole except for its inferior side. This configuration appears to be effective leading physiological loads ? which run mostly in the inferior-superior direction ? away from the hole, but it is likely not as effective spreading the medial-lateral loads of a sideways fall. This finding points toward vascularity playing a larger role during sideways falls than anticipated. VOI C was different than the other two volumes because it was selected based on the results from the organ-level FE model instead of the ex-vivo results. The micro-level results show a clear line of high strains developing in the middle of this volume and the bone surface folding over it. Comparison of these strains to the high-speed video showed that, while the main fracture line does not directly cross VOI C, there is a secondary fracture that becomes visible after the main fracture and just superior to it. The high-strain line of the VOI C model conforms to the path of the secondary fracture line. It is possible that the secondary fracture developed under the surface, covered by the layer of paint, and only became visible at a later time point. 4 ? Multi-scale nested simulation 75 4.4 Discussion The purpose of this short investigation was to test the viability of multi-scale nested FE modelling of the femur. The results demonstrate that the technique is feasible and potentially useful for identifying the kind of bone features that lead to the creation of fractures during sideways falls. Three VOIs were analysed, with the first two being empirical points of fracture initiation and the last one being a point of strain concentration in the full-bone organ-level FE model. The FE models for all three volumes showed deformation and strain patterns consistent with those seen on the ex-vivo model. Even though the strain magnitudes reported by the FE models did not match those measured experimentally from DIC, the similarity shown by the patterns and locations of high strain demonstrates that the technique already captures the general behaviour Figure 4-3: Comparison of minimum principal strains and deformations resulting from micro-level FE models to fracture locations from high-speed video 4 ? Multi-scale nested simulation 76 of the actual bone. Further advances in the modelling of bone material on FE models, particularly of its nonlinear behaviour and tension-compression asymmetry, and further improvements in the measurement of strains through DIC will likely bring the magnitudes closer together. The results from VOI B are of particular interest, as the micro modelling of this volume brought to our attention the presence of a vascular hole where the fracture line first became visible. This region appeared to be structurally sound in the organ-level model and from visual inspection of the bone. While it is important to keep in mind that this was a single specimen, the structural weakness discovered in this bone warrants further investigation of any defects around the locations where fractures initiate during sideways fall loading. There are many limitations to the presented implementation of multi-scale FE modelling. Starting on the experimental protocol, this pilot test lacked modelling of soft tissues or hip compliance, resulting in an impact much shorter than what has been described in published studies as appropriate for sideways falls [61]. The presented DIC technique also requires refinement, as here its accuracy is limited to surfaces approximately perpendicular to the camera view. Future\t ?iterations\t ?of\t ?our\t ?group?s\t ?experimental\t ?model\t ?will\t ?use two cameras to perform stereo-DIC. Furthermore, there is a need to reduce the amount of liquid splashing in order to reduce tracking artifacts. This could be done by thoroughly wiping excess liquid from the bone and support surfaces before testing. On the FE side, there is a need for better modelling of bone material properties. In particular, it is important to implement tension-compression asymmetry; the current implementation contains post-yield softening and densification regions in tension, where they clearly do not apply. Anisotropy was similarly neglected in the current models. There is further need for validation of post-yield and rate-dependent behaviours of bone material for dynamic FE modelling. In summary, this investigation demonstrates the current state of our multi-scale nested dynamic FE modelling technique. While this early version of the technique clearly requires improvements in multiple areas of modelling, it is here shown that the method is feasible and that it has the potential to yield valuable data on the physical characteristics of bone at fracture initiation points. 77 5 Conclusions This investigation focused on the development of an innovative methodology to generate finite element models of the proximal femur under conditions present during sideways falls. Several novel features ? most significantly the use of an explicit dynamic FE solver, the modelling of materials at very high loading rates, the comparison to experimental data from dynamic impact testing of femurs, and the development of a multi-scale modelling technique ? were implemented and tested. This thesis documents the first steps on the development of a promising new approach to numerical modelling of hip fractures; an approach that complements the work done by Gilchrist for in-vitro mechanical testing of femurs. In chapter 2, a semi-automated and modular technique was developed and used to generate implicit quasi-static FE models of 16 bones. Stiffness and strain results from these models were compared against quasi-static data collected experimentally in a materials testing machine. Results were further compared using three modulus-density relationships and three techniques for mapping of heterogeneous material properties. Modulus-density relationships demonstrated to have an important role in FE results, with the stiffest relationship producing the best predictions of surface strains and the most compliant relationship producing the most accurate predictions of overall femoral stiffness. Material mapping techniques had a comparatively weak effect on FE results. Boundary conditions proved to be of great importance during model development. In particular, we found that neglecting to model a contact-based sliding support at the head or a bendable distal support at the shaft easily led to overconstrained models. Chapter 3 expanded the technique to generate explicit dynamic FE models for 15 of the bones tested in the previous chapter. The FE results for stiffness, peak force, and strain energy at peak force were compared against dynamic data collected experimentally using a drop tower device. Even though the stiffest modulus-density relationship from chapter 2 was used, and that models were further stiffened through a strain-rate factor, the dynamic results showed a general underestimation of stiffness for the stiffest specimens. FE surface strains were found to match experimental fracture patterns observed in high-speed videos. The sensitivity of the model was also tested by varying the stiffness of the foam (used experimentally to represent soft tissues), the nonlinear behaviour of the bone material, and the modulus-density relationship of the bone. The results suggest that the foam could be eliminated from the modelling methodology, while the effects of strain-rate hardening and the modulus-density relationship should be studied 5 ? Conclusions 78 further. In general, the results support the feasibility of explicit FE modelling and the need for further investigation of impact loading on the femur. Chapter 4 explored the feasibility of multi-scale nested FE modelling. A preliminary version of this technique was presented, demonstrating that the technique is not only feasible but also potentially useful to identify microstructural features of bone that might make it prone to fracture. Further research is definitely required to make compelling conclusions about the validity of this modelling technique, but the results so far certainly encourage further work. 79 References [1] O.\t ?Johnell\t ?and\t ?J.\t ?Kanis,\t ??An\t ?estimate\t ?of\t ?the\t ?worldwide\t ?prevalence\t ?and\t ?disability\t ?associated\t ?with osteoporotic\t ?fractures,?\t ?Osteoporos Int, vol. 17, no. 12, pp. 1726?1733, Sep. 2006. [2] E.\t ?A.\t ?Papadimitropoulos,\t ?P.\t ?C.\t ?Coyte,\t ?R.\t ?G.\t ?Josse,\t ?and\t ?C.\t ?E.\t ?Greenwood,\t ??Current\t ?and\t ?projected\t ?rates\t ?of\t ?hip\t ?fracture\t ?in\t ?Canada,?\t ?CMAJ, vol. 157, no. 10, pp. 1357?1363, Nov. 1997. [3] C.\t ?L.\t ?Leibson,\t ?A.\t ?N.\t ?A.\t ?Tosteson,\t ?S.\t ?E.\t ?Gabriel,\t ?J.\t ?E.\t ?Ransom,\t ?and\t ?L.\t ?J.\t ?Melton,\t ??Mortality,\t ?Disability, and Nursing Home Use for Persons with and without Hip Fracture: A Population-Based\t ?Study,?\t ?Journal of the American Geriatrics Society, vol. 50, no. 10, pp. 1644?1650, 2002. [4] M. E. Wiktorowicz, R. Goeree, A. Papaioannou, J. D. Adachi, and E. Papadimitropoulos, ?Economic\t ?Implications\t ?of\t ?Hip\t ?Fracture:\t ?Health\t ?Service\t ?Use,\t ?Institutional\t ?Care\t ?and\t ?Cost\t ?in\t ?Canada,?\t ?Osteoporosis International, vol. 12, no. 4, pp. 271?278, May 2001. [5] B.\t ?Gullberg,\t ?O.\t ?Johnell,\t ?and\t ?J.\t ?A.\t ?Kanis,\t ??World-wide\t ?Projections\t ?for\t ?Hip\t ?Fracture,?\t ?Osteoporos Int, vol. 7, no. 5, pp. 407?413, Sep. 1997. [6] J.\t ?A.\t ?Kanis,\t ??Diagnosis\t ?of\t ?osteoporosis\t ?and\t ?assessment\t ?of\t ?fracture\t ?risk,?\t ?Lancet, vol. 359, no. 9321, pp. 1929?36, 2002. [7] O. Johnell, J. A. Kanis, A. Oden, H. Johansson, C. De Laet, P. Delmas, J. A. Eisman, S. Fujiwara, H.\t ?Kroger,\t ?D.\t ?Mellstrom,\t ?P.\t ?J.\t ?Meunier,\t ?L.\t ?J.\t ?Melton\t ?III,\t ?T.\t ?O?Neill,\t ?H.\t ?Pols,\t ?J.\t ?Reeve,\t ?A.\t ?Silman,\t ?and A. Tenenhouse,\t ??Predictive\t ?Value\t ?of\t ?BMD\t ?for\t ?Hip\t ?and\t ?Other\t ?Fractures,?\t ?J Bone Min Res, vol. 20, no. 7, pp. 1185?1194, Jul. 2005. [8] R.\t ?Marks,\t ?J.\t ?P.\t ?Allegrante,\t ?C.\t ?R.\t ?MacKenzie,\t ?and\t ?J.\t ?M.\t ?Lane,\t ??Hip\t ?fractures\t ?among\t ?the\t ?elderly:\t ?causes,\t ?consequences\t ?and\t ?control,?\t ?Ageing Research Reviews, vol. 2, no. 1, pp. 57?93, Jan. 2003. [9] S.\t ?L.\t ?Greenspan,\t ?E.\t ?R.\t ?Myers,\t ?L.\t ?A.\t ?Maitland,\t ?N.\t ?M.\t ?Resnick,\t ?and\t ?W.\t ?C.\t ?Hayes,\t ??Fall\t ?Severity\t ?and\t ?Bone\t ?Mineral\t ?Density\t ?as\t ?Risk\t ?Factors\t ?for\t ?Hip\t ?Fracture\t ?in\t ?Ambulatory\t ?Elderly,?\t ?JAMA: The Journal of the American Medical Association, vol. 271, no. 2, pp. 128 ?133, Jan. 1994. [10] M.\t ?Peacock,\t ?C.\t ?H.\t ?Turner,\t ?G.\t ?Liu,\t ?A.\t ?K.\t ?Manatunga,\t ?L.\t ?Timmerman,\t ?and\t ?C.\t ?C.\t ?J.\t ?Jr,\t ??Better\t ?discrimination\t ?of\t ?hip\t ?fracture\t ?using\t ?bone\t ?density,\t ?geometry\t ?and\t ?architecture,?\t ?Osteoporosis Int, vol. 5, no. 3, pp. 167?173, May 1995. [11] J.\t ?C.\t ?Lotz,\t ?E.\t ?J.\t ?Cheal,\t ?and\t ?W.\t ?C.\t ?Hayes,\t ??Fracture\t ?prediction\t ?for\t ?the\t ?proximal\t ?femur\t ?using\t ?finite element models: Part I--Linear\t ?analysis,?\t ?J Biomech Eng, vol. 113, no. 4, pp. 353?60, 1991. [12] J. H.\t ?Keyak,\t ?S.\t ?A.\t ?Rossi,\t ?K.\t ?A.\t ?Jones,\t ?and\t ?H.\t ?B.\t ?Skinner,\t ??Prediction\t ?of\t ?femoral\t ?fracture\t ?load\t ?using\t ?automated\t ?finite\t ?element\t ?modeling,?\t ?Journal of Biomechanics, vol. 31, no. 2, pp. 125?133, May 1998. [13] M. Bessho, I. Ohnishi, J. Matsuyama, T. Matsumoto, K.\t ?Imai,\t ?and\t ?K.\t ?Nakamura,\t ??Prediction\t ?of\t ?strength and strain of the proximal femur by a CT-based\t ?finite\t ?element\t ?method,?\t ?Journal of Biomechanics, vol. 40, no. 8, pp. 1745?1753, 2007. [14] E.\t ?Schileo,\t ?F.\t ?Taddei,\t ?L.\t ?Cristofolini,\t ?and\t ?M.\t ?Viceconti,\t ??Subject-specific finite element models implementing a maximum principal strain criterion are able to estimate failure risk and\t ?fracture\t ?location\t ?on\t ?human\t ?femurs\t ?tested\t ?in\t ?vitro,?\t ?J Biomech, vol. 41, no. 2, pp. 356?367, 2008. [15] Gray?s\t ?anatomy:\t ?the\t ?anatomical\t ?basis\t ?of\t ?clinical\t ?practice,\t ?39th\t ?ed.\t ?Edinburgh,\t ?Scotland?;\t ?Toronto, ON: Elsevier Churchill Livingstone, 2005. References 80 [16] R. Drake, W. Vogl, and A. Mitchell, Gray?s\t ?Anatomy\t ?for\t ?Students, 1st ed. Churchill Livingstone, 2004. [17] R. Bartl and B. Frisch, Osteoporosis: diagnosis, prevention, therapy, 2nd ed. Berlin: Springer, 2009. [18] M. Butler, M. Forte, R. L. Kane, S. Joglekar, S. J. Duval, M. Swiontkowski, and T. Wilt, ?Treatment\t ?of\t ?Common\t ?Hip\t ?Fractures,?\t ?Agency\t ?for\t ?Healthcare\t ?Research\t ?and\t ?Quality\t ?(US),\t ?Rockville (MD), 09-E013, Aug. 2009. [19] J.\t ?D.\t ?Michelson,\t ?A.\t ?S.\t ?Myers,\t ?R.\t ?Jinnah,\t ?Q.\t ?Cox,\t ?and\t ?M.\t ?M.\t ?Van\t ?Natta,\t ??Epidemiology\t ?of\t ?Hip\t ?Fractures\t ?Among\t ?the\t ?Elderly:\t ?Risk\t ?Factors\t ?for\t ?Fracture\t ?Type,?\t ?Clinical Orthopaedics and Related Research, pp. 129?135, Feb. 1995. [20] R. B. Birrer, Field guide to fracture management. Philadelphia, PA: Lippincott Williams & Wilkins, 2005. [21] L.\t ?Cristofolini,\t ?M.\t ?Juszczyk,\t ?S.\t ?Martelli,\t ?F.\t ?Taddei,\t ?and\t ?M.\t ?Viceconti,\t ??In\t ?vitro\t ?replication\t ?of\t ?spontaneous\t ?fractures\t ?of\t ?the\t ?proximal\t ?human\t ?femur,?\t ?J Biomech, vol. 40, no. 13, pp. 2837?45, 2007. [22] M. Viceconti, F. Taddei, L. Cristofolini, S. Martelli, C.\t ?Falcinelli,\t ?and\t ?E.\t ?Schileo,\t ??Are\t ?spontaneous fractures possible? An example of clinical application for personalised, multiscale neuro-musculo-skeletal\t ?modelling,?\t ?J Biomech, vol. 45, no. 3, pp. 421?426, Feb. 2012. [23] J. Parkkari, P. Kannus, M. Palvanen, A. Natri, J. Vainio, H. Aho, I. Vuori, and M. J?rvinen, ?Majority\t ?of\t ?Hip\t ?Fractures\t ?Occur\t ?as\t ?a\t ?Result\t ?of\t ?a\t ?Fall\t ?and\t ?Impact\t ?on\t ?the\t ?Greater\t ?Trochanter\t ?of\t ?the\t ?Femur:\t ?A\t ?Prospective\t ?Controlled\t ?Hip\t ?Fracture\t ?Study\t ?with\t ?206\t ?Consecutive\t ?Patients,?\t ?Calcified Tissue International, vol. 65, no. 3, pp. 183?187, Sep. 1999. [24] N.\t ?M.\t ?Nachreiner,\t ?M.\t ?J.\t ?Findorff,\t ?J.\t ?F.\t ?Wyman,\t ?and\t ?T.\t ?C.\t ?McCarthy,\t ??Circumstances\t ?and\t ?consequences of falls in community-dwelling\t ?older\t ?women,?\t ?J Womens Health (Larchmt), vol. 16, no. 10, pp. 1437?1446, Dec. 2007. [25] T. Yoshikawa, C. h. Turner, M. Peacock, C. w. Slemenda, C. m. Weaver, D. Teegarden, P. Markwardt,\t ?and\t ?D.\t ?b.\t ?Burr,\t ??Geometric\t ?structure\t ?of\t ?the\t ?femoral\t ?neck\t ?measured\t ?using\t ?dual-energy X-ray\t ?absorptiometry,?\t ?Journal of Bone and Mineral Research, vol. 9, no. 7, pp. 1053?1064, 1994. [26] P.\t ?M.\t ?de\t ?Bakker,\t ?S.\t ?L.\t ?Manske,\t ?V.\t ?Ebacher,\t ?T.\t ?R.\t ?Oxland,\t ?P.\t ?A.\t ?Cripton,\t ?and\t ?P.\t ?Guy,\t ??During\t ?sideways falls proximal femur fractures initiate in the superolateral cortex: evidence from high-speed video of\t ?simulated\t ?fractures,?\t ?J Biomech, vol. 42, no. 12, pp. 1917?1925, Aug. 2009. [27] C.\t ?H.\t ?Turner,\t ??The\t ?biomechanics\t ?of\t ?hip\t ?fracture,?\t ?The Lancet, vol. 366, pp. 98?99, Jul. 2005. [28] P. M. Mayhew, C. D. Thomas, J. G. Clement, N. Loveridge, T. J. Beck, W. Bonfield, C. J. Burgoyne,\t ?and\t ?J.\t ?Reeve,\t ??Relation\t ?between\t ?age,\t ?femoral\t ?neck\t ?cortical\t ?stability,\t ?and\t ?hip\t ?fracture\t ?risk,?\t ?The Lancet, vol. 366, no. 9480, pp. 129?135, 2005. [29] H. Boehm, A. Horng, M. Notohamiprodjo, F. Eckstein, D. Burklein, A. Panteleon, J. Lutz, and M.\t ?Reiser,\t ??Prediction\t ?of\t ?the\t ?fracture\t ?load\t ?of\t ?whole\t ?proximal\t ?femur\t ?specimens\t ?by\t ?topological analysis of the mineral distribution in DXA-scan\t ?images,?\t ?Bone, vol. 43, no. 5, pp. 826?831, Nov. 2008. [30] A. C. Courtney, E. F. Wachtel, E. R. Myers, and\t ?W.\t ?C.\t ?Hayes,\t ??Effects\t ?of\t ?loading\t ?rate\t ?on\t ?strength\t ?of\t ?the\t ?proximal\t ?femur,?\t ?Calcif Tissue Int, vol. 55, no. 1, pp. 53?58, Jul. 1994. [31] F. Eckstein, C. Wunderer, H. Boehm, V. Kuhn, M. Priemel, T. M. Link, and E.-M. Lochm?ller, ?Reproducibility\t ?and\t ?side differences of mechanical tests for determining the structural strength\t ?of\t ?the\t ?proximal\t ?femur,?\t ?J. Bone Miner. Res, vol. 19, no. 3, pp. 379?385, Mar. 2004. References 81 [32] J.\t ?H.\t ?Keyak,\t ??Relationships\t ?between\t ?femoral\t ?fracture\t ?loads\t ?for\t ?two\t ?load\t ?configurations,?\t ?J Biomech, vol. 33, no. 4, pp. 499?502, 2000. [33] E.\t ?M.\t ?Lochm?ller,\t ?O.\t ?Groll,\t ?V.\t ?Kuhn,\t ?and\t ?F.\t ?Eckstein,\t ??Mechanical\t ?strength\t ?of\t ?the\t ?proximal\t ?femur as predicted from geometric and densitometric bone properties at the lower limb versus\t ?the\t ?distal\t ?radius,?\t ?Bone, vol. 30, no. 1, pp. 207?216, Jan. 2002. [34] G.\t ?Holzer,\t ?G.\t ?von\t ?Skrbensky,\t ?L.\t ?A.\t ?Holzer,\t ?and\t ?W.\t ?Pichl,\t ??Hip\t ?Fractures\t ?and\t ?the\t ?Contribution\t ?of\t ?Cortical\t ?Versus\t ?Trabecular\t ?Bone\t ?to\t ?Femoral\t ?Neck\t ?Strength,?\t ?Journal of Bone and Mineral Research, vol. 24, no. 3, pp. 468?474, Mar. 2009. [35] T.\t ?P.\t ?Pinilla,\t ?K.\t ?C.\t ?Boardman,\t ?M.\t ?L.\t ?Bouxsein,\t ?E.\t ?R.\t ?Myers,\t ?and\t ?W.\t ?C.\t ?Hayes,\t ??Impact\t ?direction\t ?from a fall influences the failure load of the proximal femur as much as age-related bone loss,?\t ?Calcif. Tissue Int, vol. 58, no. 4, pp. 231?235, Apr. 1996. [36] C.\t ?M.\t ?Ford,\t ?T.\t ?M.\t ?Keaveny,\t ?and\t ?W.\t ?C.\t ?Hayes,\t ??The\t ?effect\t ?of\t ?impact\t ?direction\t ?on\t ?the\t ?structural\t ?capacity\t ?of\t ?the\t ?proximal\t ?femur\t ?during\t ?falls,?\t ?Journal of Bone and Mineral Research, vol. 11, no. 3, pp. 377?383, Mar. 1996. [37] T. Weber,\t ?K.\t ?Yang,\t ?R.\t ?Woo,\t ?and\t ?R.\t ?J.\t ?Fitzgerald,\t ??Proximal\t ?femur\t ?strength:\t ?Correlation\t ?of\t ?the\t ?rate\t ?of\t ?loading\t ?and\t ?bone\t ?mineral\t ?density,?\t ?ASME Advances in Bioengineering, vol. 22, pp. 111?114, 1992. [38] A. J. van den Kroonenberg, W. C. Hayes, and T. A. McMahon, ?Hip\t ?impact\t ?velocities\t ?and\t ?body\t ?configurations\t ?for\t ?voluntary\t ?falls\t ?from\t ?standing\t ?height,?\t ?J Biomech, vol. 29, no. 6, pp. 807?811, Jun. 1996. [39] F.\t ?Feldman\t ?and\t ?S.\t ?N.\t ?Robinovitch,\t ??Reducing\t ?hip\t ?fracture\t ?risk\t ?during\t ?sideways\t ?falls:\t ?evidence in young adults of\t ?the\t ?protective\t ?effects\t ?of\t ?impact\t ?to\t ?the\t ?hands\t ?and\t ?stepping,?\t ?J Biomech, vol. 40, no. 12, pp. 2612?2618, 2007. [40] A.\t ?C.\t ?Laing\t ?and\t ?S.\t ?N.\t ?Robinovitch,\t ??The\t ?Force\t ?Attenuation\t ?Provided\t ?by\t ?Hip\t ?Protectors\t ?Depends on Impact Velocity, Pelvic Size, and Soft Tissue\t ?Stiffness,?\t ?J. Biomech. Eng., vol. 130, no. 6, pp. 061005?9, Dec. 2008. [41] S.\t ?Gilchrist,\t ?P.\t ?Guy,\t ?and\t ?P.\t ?A.\t ?Cripton,\t ??Development\t ?of\t ?an\t ?Inertia-Driven Model of Sideways Fall\t ?for\t ?Detailed\t ?Study\t ?of\t ?Femur\t ?Fracture\t ?Mechanics,?\t ?J Biomech Eng, vol. 135, no. 12, pp. 121001?121001, Oct. 2013. [42] O.\t ?C.\t ?Zienkiewicz,\t ??The\t ?birth\t ?of\t ?the\t ?finite\t ?element\t ?method\t ?and\t ?of\t ?computational\t ?mechanics,?\t ?International Journal for Numerical Methods in Engineering, vol. 60, no. 1, pp. 3?10, 2004. [43] O. Zienkiewicz, R. L. Taylor, J. Z. Zhu, and P. Nithiarasu, The Finite Element Method, 6th ed., vol. 1, 3 vols. New York: Elsevier/Butterworth-Heinemann, 2005. [44] J. N. Reddy, Introduction to Nonlinear Finite Element Analysis. Cary, NC, USA: Oxford University Press, 2004. [45] S. R. Wu and L. Gu, Introduction to the Explicit Finite Element Method for Nonlinear Transient Dynamics. Somerset, NJ, USA: Wiley, 2012. [46] A. A. Becker, An Introductory Guide to Finite Element Analysis. New York, NY, USA: ASME, 2004. [47] J.\t ?P.\t ?Halloran,\t ?A.\t ?J.\t ?Petrella,\t ?and\t ?P.\t ?J.\t ?Rullkoetter,\t ??Explicit\t ?finite\t ?element\t ?modeling\t ?of\t ?total\t ?knee\t ?replacement\t ?mechanics,?\t ?Journal of Biomechanics, vol. 38, no. 2, pp. 323?331, Feb. 2005. [48] W. A. M. Brekelmans, H. W. Poort, and T. J. J. H. Slooff,\t ??A\t ?New\t ?Method\t ?to\t ?Analyse\t ?the\t ?Mechanical\t ?Behaviour\t ?of\t ?Skeletal\t ?Parts,?\t ?ACTA Orthop. Scand, vol. 43, no. 5, pp. 301?317, 1972. References 82 [49] H.\t ?H.\t ?Vichnin\t ?and\t ?S.\t ?C.\t ?Batterman,\t ??Stress\t ?Analysis\t ?and\t ?Failure\t ?Prediction\t ?in\t ?the\t ?Proximal\t ?Femur Before and After Total\t ?Hip\t ?Replacement,?\t ?J Biomech Eng, vol. 108, no. 1, pp. 33?41, Feb. 1986. [50] J.\t ?C.\t ?Lotz,\t ?E.\t ?J.\t ?Cheal,\t ?and\t ?W.\t ?C.\t ?Hayes,\t ??Fracture\t ?prediction\t ?for\t ?the\t ?proximal\t ?femur\t ?using\t ?finite element models: Part II--Nonlinear\t ?analysis,?\t ?J Biomech Eng, vol. 113, no. 4, pp. 361?5, 1991. [51] D.\t ?D.\t ?Cody,\t ?Gross,\t ?F.\t ?J.\t ?Hou,\t ?H.\t ?J.\t ?Spencer,\t ?S.\t ?A.\t ?Goldstein,\t ?and\t ?D.\t ?P.\t ?Fyhrie,\t ??Femoral\t ?strength\t ?is\t ?better\t ?predicted\t ?by\t ?finite\t ?element\t ?models\t ?than\t ?QCT\t ?and\t ?DXA,?\t ?Journal of Biomechanics, vol. 32, no. 10, pp. 1013?1020, Oct. 1999. [52] E.\t ?Dall?Ara,\t ?B.\t ?Luisier,\t ?R.\t ?Schmidt,\t ?F.\t ?Kainberger,\t ?P.\t ?Zysset,\t ?and\t ?D.\t ?Pahr,\t ??A\t ?nonlinear\t ?QCT-based finite element model validation study for the human femur tested in two configurations\t ?in\t ?vitro,?\t ?Bone, vol. 52, no. 1, pp. 27?38, Jan. 2013. [53] L. Duchemin,\t ?D.\t ?Mitton,\t ?E.\t ?Jolivet,\t ?V.\t ?and\t ?Bousson,\t ?J.\t ?D.\t ?and\t ?Laredo,\t ?and\t ?W.\t ?and\t ?Skalli,\t ??An\t ?anatomical subject-specific FE-model\t ?for\t ?hip\t ?fracture\t ?load\t ?prediction,?\t ?Computer Methods in Biomechanics and Biomedical Engineering, vol. 11, no. 2, pp. 105?111, 2008. [54] J.\t ?H.\t ?Keyak,\t ?S.\t ?A.\t ?Rossi,\t ?K.\t ?A.\t ?Jones,\t ?C.\t ?M.\t ?Les,\t ?and\t ?H.\t ?B.\t ?Skinner,\t ??Prediction\t ?of\t ?fracture\t ?location\t ?in\t ?the\t ?proximal\t ?femur\t ?using\t ?finite\t ?element\t ?models,?\t ?Medical Engineering and Physics, vol. 23, no. 9, pp. 657?664, 2001. [55] E. Schileo, F. Taddei, A. Malandrino,\t ?L.\t ?Cristofolini,\t ?and\t ?M.\t ?Viceconti,\t ??Subject-specific finite element\t ?models\t ?can\t ?accurately\t ?predict\t ?strain\t ?levels\t ?in\t ?long\t ?bones,?\t ?Journal of Biomechanics, vol. 40, no. 13, pp. 2982?2989, 2007. [56] N. Trabelsi, Z. Yosibash, C. Wutte, P. Augat, and S.\t ?Eberle,\t ??Patient-specific finite element analysis of the human femur?A double-blinded\t ?biomechanical\t ?validation,?\t ?Journal of Biomechanics, vol. 44, no. 9, pp. 1666?1672, Jun. 2011. [57] H.\t ?Wille,\t ?E.\t ?Rank,\t ?and\t ?Z.\t ?Yosibash,\t ??Prediction\t ?of\t ?the\t ?mechanical\t ?response of the femur with uncertain\t ?elastic\t ?properties,?\t ?Journal of Biomechanics, vol. 45, no. 7, pp. 1140?1148, Apr. 2012. [58] D. Dragomir-Daescu, J. Op Den Buijs, S. McEligot, Y. Dai, R. C. Entwistle, C. Salas, L. J. Melton, K. E. Bennet, S. Khosla, and S.\t ?Amin,\t ??Robust\t ?QCT/FEA\t ?models\t ?of\t ?proximal\t ?femur\t ?stiffness\t ?and\t ?fracture\t ?load\t ?during\t ?a\t ?sideways\t ?fall\t ?on\t ?the\t ?hip,?\t ?Annals of biomedical engineering, vol. 39, no. 2, pp. 742?755, 2011. [59] L. Grassi, E. Schileo, F. Taddei, L. Zani, M. Juszczyk, L. Cristofolini, and M. Viceconti, ?Accuracy\t ?of\t ?finite\t ?element\t ?predictions\t ?in\t ?sideways\t ?load\t ?configurations\t ?for\t ?the\t ?proximal\t ?human\t ?femur,?\t ?J Biomech, vol. 45, no. 2, pp. 394?399, Jan. 2012. [60] J. E. M. Koivum?ki, J. Thevenot, P. Pulkkinen, V. Kuhn, T. M. Link, F. Eckstein, and T. J?ms?, ?CT-based finite element models can be used to estimate experimentally measured failure loads\t ?in\t ?the\t ?proximal\t ?femur,?\t ?Bone, vol. 50, no. 4, pp. 824?829, Apr. 2012. [61] S. N. Robinovitch, S. L. Evans, J. Minns, A. C. Laing, P. Kannus, P. A. Cripton, S. Derler, S. J. Birge,\t ?D.\t ?Plant,\t ?I.\t ?D.\t ?Cameron,\t ?D.\t ?P.\t ?Kiel,\t ?J.\t ?Howland,\t ?K.\t ?Khan,\t ?and\t ?J.\t ?B.\t ?Lauritzen,\t ??Hip\t ?protectors: recommendations for biomechanical testing?an international consensus statement\t ?(part\t ?I),?\t ?Osteoporos Int, vol. 20, no. 12, pp. 1977?1988, Oct. 2009. [62] M. Bessho, I. Ohnishi, T. Matsumoto, S. Ohashi, J. Matsuyama, K. Tobita, M. Kaneko, and K. Nakamura,\t ??Prediction\t ?of\t ?proximal\t ?femur\t ?strength\t ?using\t ?a\t ?CT-based nonlinear finite element method: Differences in predicted fracture load and site with changing load and boundary\t ?conditions,?\t ?Bone, vol. 45, no. 2, pp. 226?231, Aug. 2009. [63] B. Helgason, S. Gilchrist, O. Ariza, J. D. Chak, G. Zheng, R. Widmer, S. J. Ferguson, P. Guy, and P.\t ?A.\t ?Cripton,\t ??Development\t ?of\t ?a\t ?balanced\t ?experimental-computational approach to References 83 understanding\t ?the\t ?mechanics\t ?of\t ?proximal\t ?femur\t ?fractures,?\t ?Medical Engineering & Physics, no. In Press, 2014. [64] Z.\t ?Yosibash,\t ?N.\t ?Trabelsi,\t ?and\t ?C.\t ?Milgrom,\t ??Reliable\t ?simulations\t ?of\t ?the\t ?human\t ?proximal\t ?femur\t ?by high-order finite\t ?element\t ?analysis\t ?validated\t ?by\t ?experimental\t ?observations,?\t ?Journal of Biomechanics, vol. 40, no. 16, pp. 3688?3699, 2007. [65] J.\t ?H.\t ?Keyak,\t ??Improved\t ?prediction\t ?of\t ?proximal\t ?femoral\t ?fracture\t ?load\t ?using\t ?nonlinear\t ?finite\t ?element\t ?models,?\t ?Med Eng Phys, vol. 23, no. 3, pp. 165?73, 2001. [66] L.\t ?Peng,\t ?J.\t ?Bai,\t ?X.\t ?Zeng,\t ?and\t ?Y.\t ?Zhou,\t ??Comparison\t ?of\t ?isotropic\t ?and\t ?orthotropic\t ?material\t ?property\t ?assignments\t ?on\t ?femoral\t ?finite\t ?element\t ?models\t ?under\t ?two\t ?loading\t ?conditions,?\t ?Medical Engineering & Physics, vol. 28, no. 3, pp. 227?233, Apr. 2006. [67] V.\t ?Baca,\t ?Z.\t ?Horak,\t ?P.\t ?Mikulenka,\t ?and\t ?V.\t ?Dzupa,\t ??Comparison\t ?of\t ?an\t ?inhomogeneous\t ?orthotropic\t ?and\t ?isotropic\t ?material\t ?models\t ?used\t ?for\t ?FE\t ?analyses,?\t ?Medical Engineering & Physics, vol. 30, no. 7, pp. 924?930, Sep. 2008. [68] H. Yang,\t ?X.\t ?Ma,\t ?and\t ?T.\t ?Guo,\t ??Some\t ?factors\t ?that\t ?affect\t ?the\t ?comparison\t ?between\t ?isotropic\t ?and\t ?orthotropic\t ?inhomogeneous\t ?finite\t ?element\t ?material\t ?models\t ?of\t ?femur,?\t ?Medical Engineering & Physics, vol. 32, no. 6, pp. 553?560, Jul. 2010. [69] N. Trabelsi, Z. Yosibash, and\t ?C.\t ?Milgrom,\t ??Validation\t ?of\t ?subject-specific automated p-FE analysis\t ?of\t ?the\t ?proximal\t ?femur,?\t ?Journal of Biomechanics, vol. 42, no. 3, pp. 234?241, Feb. 2009. [70] M.\t ?J.\t ?Ciarelli,\t ?S.\t ?A.\t ?Goldstein,\t ?J.\t ?L.\t ?Kuhn,\t ?D.\t ?D.\t ?Cody,\t ?and\t ?M.\t ?B.\t ?Brown,\t ??Evaluation\t ?of\t ?orthogonal mechanical properties and density of human trabecular bone from the major metaphyseal\t ?regions\t ?with\t ?materials\t ?testing\t ?and\t ?computed\t ?tomography,?\t ?J Orthop Res, vol. 9, no. 5, pp. 674?82, 1991. [71] J.\t ?H.\t ?Keyak,\t ?I.\t ?Y.\t ?Lee,\t ?and\t ?H.\t ?B.\t ?Skinner,\t ??Correlations\t ?between\t ?orthogonal\t ?mechanical\t ?properties\t ?and\t ?density\t ?of\t ?trabecular\t ?bone:\t ?Use\t ?of\t ?different\t ?densitometric\t ?measures,?\t ?Journal of Biomedical Materials Research, vol. 28, no. 11, pp. 1329?1336, 1994. [72] E.\t ?F.\t ?Morgan,\t ?H.\t ?H.\t ?Bayraktar,\t ?and\t ?T.\t ?M.\t ?Keaveny,\t ??Trabecular\t ?bone\t ?modulus?density relationships\t ?depend\t ?on\t ?anatomic\t ?site,?\t ?Journal of Biomechanics, vol. 36, no. 7, pp. 897?904, Jul. 2003. [73] T.\t ?S.\t ?Keller,\t ??Predicting\t ?the\t ?compressive\t ?mechanical\t ?behavior\t ?of\t ?bone,?\t ?Journal of biomechanics, vol. 27, no. 9, pp. 1159?1168, 1994. [74] L. Duchemin, V. Bousson, C. Raossanaly, C. Bergot, J. D. Laredo, W. Skalli, and D. Mitton, ?Prediction\t ?of\t ?mechanical\t ?properties of cortical bone by quantitative computed tomography,?\t ?Medical Engineering & Physics, vol. 30, no. 3, pp. 321?328, Apr. 2008. [75] J.\t ?H.\t ?Keyak\t ?and\t ?Y.\t ?Falkinstein,\t ??Comparison\t ?of\t ?in\t ?situ\t ?and\t ?in\t ?vitro\t ?CT\t ?scan-based finite element model predictions of\t ?proximal\t ?femoral\t ?fracture\t ?load,?\t ?Med Eng Phys, vol. 25, no. 9, pp. 781?7, 2003. [76] B.\t ?van\t ?Rietbergen,\t ??Trabecular\t ?bone\t ?tissue\t ?strains\t ?in\t ?the\t ?healthy\t ?and\t ?osteoporotic\t ?human\t ?femur,?\t ?Journal of Bone Mineral Research, vol. 18, no. 10, 2003. [77] S. Gilchrist,\t ??Comparison\t ?of\t ?proximal\t ?femur\t ?deformations,\t ?failures\t ?and\t ?fractures\t ?in\t ?quasi-static and inertially-driven\t ?simulations\t ?of\t ?a\t ?sideways\t ?fall\t ?from\t ?standing,?\t ?University\t ?of\t ?British Columbia, Vancouver, Canada, 2014. [78] S. L. Greenspan, E. R. Myers, D. P. Kiel,\t ?R.\t ?A.\t ?Parker,\t ?W.\t ?C.\t ?Hayes,\t ?and\t ?N.\t ?M.\t ?Resnick,\t ??Fall\t ?Direction, Bone Mineral Density, and Function: Risk Factors for Hip Fracture in Frail Nursing\t ?Home\t ?Elderly,?\t ?The American Journal of Medicine, vol. 104, no. 6, pp. 539?545, Jun. 1998. [79] F. Taddei,\t ?E.\t ?Schileo,\t ?B.\t ?Helgason,\t ?L.\t ?Cristofolini,\t ?and\t ?M.\t ?Viceconti,\t ??The\t ?material\t ?mapping\t ?strategy influences the accuracy of CT-based finite element models of bones: An evaluation References 84 against\t ?experimental\t ?measurements,?\t ?Medical Engineering & Physics, vol. 29, no. 9, pp. 973?979, Nov. 2007. [80] M.\t ?Viceconti,\t ?L.\t ?Bellingeri,\t ?L.\t ?Cristofolini,\t ?and\t ?A.\t ?Toni,\t ??A\t ?comparative\t ?study\t ?on\t ?different\t ?methods\t ?of\t ?automatic\t ?mesh\t ?generation\t ?of\t ?human\t ?femurs,?\t ?Medical Engineering & Physics, vol. 20, no. 1, pp. 1?10, Apr. 1998. [81] K. M.\t ?Fox,\t ?J.\t ?Magaziner,\t ?J.\t ?R.\t ?Hebel,\t ?J.\t ?E.\t ?Kenzora,\t ?and\t ?T.\t ?M.\t ?Kashnei,\t ??Intertrochanteric\t ?Versus\t ?Femoral\t ?Neck\t ?Hip\t ?Fractures:\t ?Differential\t ?Characteristics,\t ?Treatment,\t ?and\t ?Sequelae,?\t ?J Gerontol A Biol Sci Med Sci, vol. 54, no. 12, pp. M635?M640, Dec. 1999. [82] B. Helgason, F. Taddei, H. P?lsson, E. Schileo, L. Cristofolini, M. Viceconti, and S. Brynj?lfsson,\t ??A\t ?modified\t ?method\t ?for\t ?assigning\t ?material\t ?properties\t ?to\t ?FE\t ?models\t ?of\t ?bones,?\t ?Medical Engineering & Physics, vol. 30, no. 4, pp. 444?453, May 2008. [83] Z. Yosibash,\t ?L.\t ?Joskowicz,\t ?C.\t ?Milgrom,\t ?and\t ?R.\t ?Padan,\t ??A\t ?CT-Based High-Order Finite Element Analysis of the Human Proximal Femur Compared to In-vitro\t ?Experiments,?\t ?J Biomech Eng, vol. 129, no. 3, pp. 297?309, Nov. 2006. [84] P. A. Yushkevich, J. Piven, H. C. Hazlett,\t ?R.\t ?G.\t ?Smith,\t ?S.\t ?Ho,\t ?J.\t ?C.\t ?Gee,\t ?and\t ?G.\t ?Gerig,\t ??User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency\t ?and\t ?reliability,?\t ?NeuroImage, vol. 31, no. 3, pp. 1116?1128, Jul. 2006. [85] E.\t ?Schileo,\t ?E.\t ?Dall?Ara, F. Taddei, A. Malandrino, T. Schotkamp, M. Baleani, and M. Viceconti, ?An\t ?accurate\t ?estimation\t ?of\t ?bone\t ?density\t ?improves\t ?the\t ?accuracy\t ?of\t ?subject-specific finite element\t ?models,?\t ?Journal of Biomechanics, vol. 41, no. 11, pp. 2483?2491, Aug. 2008. [86] P. Bergstr?m,\t ??Iterative\t ?Closest\t ?Point\t ?Method,?\t ?26-Nov-2007. [Online]. Available: http://www.mathworks.com/matlabcentral/fileexchange/12627. [87] K.\t ?K.\t ?Nishiyama,\t ?S.\t ?Gilchrist,\t ?P.\t ?Guy,\t ?P.\t ?Cripton,\t ?and\t ?S.\t ?K.\t ?Boyd,\t ??Proximal\t ?femur\t ?bone\t ?strength estimated by a computationally fast finite element analysis in a sideways fall configuration,?\t ?Journal of Biomechanics, vol. 46, no. 7, pp. 1231?1236, Apr. 2013. [88] S.\t ?Nawathe,\t ?H.\t ?Akhlaghpour,\t ?M.\t ?L.\t ?Bouxsein,\t ?and\t ?T.\t ?M.\t ?Keaveny,\t ??Microstructural\t ?failure\t ?mechanisms in the human\t ?proximal\t ?femur\t ?for\t ?sideways\t ?fall\t ?loading,?\t ?Journal of Bone and Mineral Research, 2013. [89] S.\t ?Boonen,\t ?P.\t ?Autier,\t ?M.\t ?Barette,\t ?D.\t ?Vanderschueren,\t ?P.\t ?Lips,\t ?and\t ?P.\t ?Haentjens,\t ??Functional\t ?outcome and quality of life following hip fracture in elderly women: a prospective controlled\t ?study,?\t ?Osteoporos Int, vol. 15, no. 2, pp. 87?94, 2004. [90] I.\t ?Hallberg,\t ?A.\t ?M.\t ?Rosenqvist,\t ?L.\t ?Kartous,\t ?O.\t ?L?fman,\t ?O.\t ?Wahlstr?m,\t ?and\t ?G.\t ?Toss,\t ??Health-related\t ?quality\t ?of\t ?life\t ?after\t ?osteoporotic\t ?fractures,?\t ?Osteoporos Int, vol. 15, no. 10, Mar. 2004. [91] J. Magaziner, L. Fredman, W. Hawkes, J. R. Hebel, S. Zimmerman, D. L. Orwig, and L. Wehren, ?Changes\t ?in\t ?Functional\t ?Status\t ?Attributable\t ?to\t ?Hip\t ?Fracture:\t ?A\t ?Comparison\t ?of\t ?Hip\t ?Fracture\t ?Patients to Community-dwelling\t ?Aged,?\t ?American Journal of Epidemiology, vol. 157, no. 11, pp. 1023 ?1031, Jun. 2003. [92] T.\t ?H.\t ?Nevalainen,\t ?L.\t ?A.\t ?Hiltunen,\t ?and\t ?P.\t ?Jalovaara,\t ??Functional\t ?ability\t ?after\t ?hip\t ?fracture\t ?among patients home-dwelling\t ?at\t ?the\t ?time\t ?of\t ?fracture,?\t ?Cent. Eur. J. Public Health, vol. 12, no. 4, pp. 211?216, Dec. 2004. [93] E. K. Osnes, C. M. Lofthus, H. E. Meyer, J. A. Falch, L. Nordsletten, I. Cappelen, and I. S. Kristiansen,\t ??Consequences\t ?of\t ?hip\t ?fracture\t ?on\t ?activities\t ?of\t ?daily\t ?life\t ?and\t ?residential\t ?needs,?\t ?Osteoporos Int, vol. 15, no. 7, pp. 567?574, Jul. 2004. [94] C.\t ?Cooper,\t ??The\t ?crippling\t ?consequences\t ?of\t ?fractures\t ?and\t ?their\t ?impact\t ?on\t ?quality\t ?of\t ?life,?\t ?Am. J. Med, vol. 103, no. 2A, p. 12S?17S; discussion 17S?19S, Aug. 1997. References 85 [95] S. Haleem, L. Lutchman, R. Mayahi, J. E. Grice, and M. J. Parker,\t ??Mortality\t ?following\t ?hip\t ?fracture:\t ?Trends\t ?and\t ?geographical\t ?variations\t ?over\t ?the\t ?last\t ?40\t ?years,?\t ?Injury, vol. 39, no. 10, pp. 1157?1163, Oct. 2008. [96] P.\t ?Geusens,\t ?K.\t ?Milisen,\t ?E.\t ?Dejaeger,\t ?and\t ?S.\t ?Boonen,\t ??Falls\t ?and\t ?fractures\t ?in\t ?postmenopausal\t ?women:\t ?a\t ?review,?\t ?J Br Menopause Soc, vol. 9, no. 3, pp. 101?106, Sep. 2003. [97] M.\t ?C.\t ?Nevitt,\t ?S.\t ?R.\t ?Cummings,\t ?S.\t ?Kidd,\t ?and\t ?D.\t ?Black,\t ??Risk\t ?Factors\t ?for\t ?Recurrent\t ?Nonsyncopal\t ?Falls,?\t ?JAMA: The Journal of the American Medical Association, vol. 261, no. 18, pp. 2663 ?2668, May 1989. [98] M.\t ?E.\t ?Tinetti,\t ?M.\t ?Speechley,\t ?and\t ?S.\t ?F.\t ?Ginter,\t ??Risk\t ?factors\t ?for\t ?falls\t ?among\t ?elderly\t ?persons\t ?living\t ?in\t ?the\t ?community,?\t ?N. Engl. J. Med, vol. 319, no. 26, pp. 1701?1707, Dec. 1988. [99] J.\t ?Lotz\t ?and\t ?W.\t ?Hayes,\t ??The\t ?use\t ?of\t ?quantitative computed tomography to estimate risk of fracture\t ?of\t ?the\t ?hip\t ?from\t ?falls,?\t ?J Bone Joint Surg Am, vol. 72, no. 5, pp. 689?700, Jun. 1990. [100] D.\t ?R.\t ?Carter\t ?and\t ?W.\t ?C.\t ?Hayes,\t ??Bone\t ?compressive\t ?strength:\t ?the\t ?influence\t ?of\t ?density\t ?and\t ?strain\t ?rate,?\t ?Science, vol. 194, no. 4270, pp. 1174?1176, Dec. 1976. [101] D.\t ?R.\t ?Carter\t ?and\t ?W.\t ?C.\t ?Hayes,\t ??The\t ?compressive\t ?behavior\t ?of\t ?bone\t ?as\t ?a\t ?two-phase porous structure,?\t ?J Bone Joint Surg Am, vol. 59, no. 7, pp. 954?62, 1977. [102] E. S. Orwoll, L. M. Marshall, C. M. Nielson, S. R. Cummings, J. Lapidus, J. A. Cauley, K. Ensrud,\t ?N.\t ?Lane,\t ?P.\t ?R.\t ?Hoffmann,\t ?D.\t ?L.\t ?Kopperdahl,\t ?and\t ?T.\t ?M.\t ?Keaveny,\t ??Finite\t ?Element\t ?Analysis\t ?of\t ?the\t ?Proximal\t ?Femur\t ?and\t ?Hip\t ?Fracture\t ?Risk\t ?in\t ?Older\t ?Men,?\t ?Journal of Bone and Mineral Research, vol. 24, no. 3, pp. 475?483, 2009. [103] W.\t ?C.\t ?Hayes,\t ?S.\t ?J.\t ?Piazza,\t ?and\t ?P.\t ?K.\t ?Zysset,\t ??Biomechanics\t ?of\t ?fracture\t ?risk\t ?prediction\t ?of\t ?the\t ?hip\t ?and\t ?spine\t ?by\t ?quantitative\t ?computed\t ?tomography,?\t ?Radiol Clin North Am, vol. 29, no. 1, pp. 1?18, 1991. [104] S.\t ?N.\t ?Robinovitch,\t ?W.\t ?C.\t ?Hayes,\t ?and\t ?T.\t ?A.\t ?McMahon,\t ??Distribution\t ?of\t ?contact\t ?force\t ?during\t ?impact\t ?to\t ?the\t ?hip,?\t ?Ann Biomed Eng, vol. 25, no. 3, pp. 499?508, May 1997. [105] D. P. Beason, G. J. Dakin, R. R. Lopez, J. E. Alonso, F. A. Bandak, and\t ?A.\t ?W.\t ?Eberhardt,\t ??Bone\t ?mineral\t ?density\t ?correlates\t ?with\t ?fracture\t ?load\t ?in\t ?experimental\t ?side\t ?impacts\t ?of\t ?the\t ?pelvis,?\t ?J Biomech, vol. 36, no. 2, pp. 219?27, 2003. [106] A.\t ?C.\t ?Laing\t ?and\t ?S.\t ?N.\t ?Robinovitch,\t ??Characterizing\t ?the\t ?effective\t ?stiffness\t ?of\t ?the\t ?pelvis during\t ?sideways\t ?falls\t ?on\t ?the\t ?hip,?\t ?J Biomech, vol. 43, no. 10, pp. 1898?1904, Jul. 2010. [107] H. H. Bayraktar, E. F. Morgan, G. L. Niebur, G. E. Morris, E. K. Wong, and T. M. Keaveny, ?Comparison\t ?of\t ?the\t ?elastic\t ?and\t ?yield\t ?properties\t ?of\t ?human\t ?femoral\t ?trabecular and cortical bone\t ?tissue,?\t ?Journal of Biomechanics, vol. 37, no. 1, pp. 27?35, Jan. 2004. [108] B.\t ?Helgason,\t ?M.\t ?Viceconti,\t ?T.\t ?P.\t ?R?narsson,\t ?and\t ?S.\t ?Brynj?lfsson,\t ??On\t ?the\t ?mechanical\t ?stability of porous coated press fit titanium implants: A finite element study of a pushout test,?\t ?Journal of Biomechanics, vol. 41, no. 8, pp. 1675?1681, 2008. [109] W.\t ?C.\t ?Hayes\t ?and\t ?D.\t ?R.\t ?Carter,\t ??Postyield\t ?behavior\t ?of\t ?subchondral\t ?trabecular\t ?bone,?\t ?Journal of Biomedical Materials Research, vol. 10, no. 4, pp. 537?544, Jul. 1976. [110] A.\t ?Oloyede,\t ?R.\t ?Flachsmann,\t ?and\t ?N.\t ?D.\t ?Broom,\t ??The\t ?Dramatic\t ?Influence\t ?of\t ?Loading\t ?Velocity\t ?on\t ?the\t ?Compressive\t ?Response\t ?of\t ?Articular\t ?Cartilage,?\t ?Connective Tissue Research, vol. 27, no. 4, pp. 211?224, 1992. [111] S. Gilchrist, K. K. Nishiyama, P. de\t ?Bakker,\t ?S.\t ?K.\t ?Boyd,\t ?T.\t ?Oxland,\t ?and\t ?P.\t ?A.\t ?Cripton,\t ??Impact\t ?and constant displacement rate loading changes the response of the femur in fall simulation,?\t ?Journal of Biomechanics, no. In Review, 2013. [112] J.\t ?D.\t ?Currey,\t ??The\t ?effects\t ?of\t ?strain\t ?rate,\t ?reconstruction and mineral content on some mechanical\t ?properties\t ?of\t ?bovine\t ?bone,?\t ?Journal of Biomechanics, vol. 8, no. 1, pp. 81?86, 1975. References 86 [113] U.\t ?Hansen,\t ?P.\t ?Zioupos,\t ?R.\t ?Simpson,\t ?J.\t ?D.\t ?Currey,\t ?and\t ?D.\t ?Hynd,\t ??The\t ?Effect\t ?of\t ?Strain\t ?Rate\t ?on\t ?the Mechanical Properties of\t ?Human\t ?Cortical\t ?Bone,?\t ?Journal of Biomechanical Engineering, vol. 130, no. 1, p. 011011, 2008. [114] J.\t ?H.\t ?McElhaney,\t ??Dynamic\t ?response\t ?of\t ?bone\t ?and\t ?muscle\t ?tissue.,?\t ?J Appl Physiol, vol. 21, no. 4, pp. 1231?1236, Jul. 1966. [115] F. Linde, P. N?rgaard, I.\t ?Hvid,\t ?A.\t ?Odgaard,\t ?and\t ?K.\t ?S?balle,\t ??Mechanical\t ?properties\t ?of\t ?trabecular\t ?bone.\t ?Dependency\t ?on\t ?strain\t ?rate,?\t ?Journal of Biomechanics, vol. 24, no. 9, pp. 803?809, 1991. [116] U.\t ?A.\t ?Gurkan\t ?and\t ?O.\t ?Akkus,\t ??The\t ?Mechanical\t ?Environment\t ?of\t ?Bone\t ?Marrow:\t ?A\t ?Review,?\t ?Ann Biomed Eng, vol. 36, no. 12, pp. 1978?1991, Dec. 2008. [117] W.\t ?C.\t ?W.\t ?Chang,\t ?T.\t ?M.\t ?Christensen,\t ?T.\t ?P.\t ?Pinilla,\t ?and\t ?T.\t ?M.\t ?Keaveny,\t ??Uniaxial\t ?yield\t ?strains\t ?for\t ?bovine\t ?trabecular\t ?bone\t ?are\t ?isotropic\t ?and\t ?asymmetric,?\t ?Journal of Orthopaedic Research, vol. 17, no. 4, pp. 582?585, 1999. [118] D.\t ?T.\t ?Reilly\t ?and\t ?A.\t ?H.\t ?Burstein,\t ??The\t ?elastic\t ?and\t ?ultimate\t ?properties\t ?of\t ?compact\t ?bone\t ?tissue,?\t ?Journal of Biomechanics, vol. 8, no. 6, pp. 393?405, 1975. [119] J.\t ?A.\t ?MacNeil\t ?and\t ?S.\t ?K.\t ?Boyd,\t ??Bone\t ?strength\t ?at\t ?the\t ?distal\t ?radius can be estimated from high-resolution peripheral quantitative computed tomography and the finite element method,?\t ?Bone, vol. 42, no. 6, pp. 1203?1213, Jun. 2008. [120] G.\t ?L.\t ?Niebur,\t ?M.\t ?J.\t ?Feldstein,\t ?J.\t ?C.\t ?Yuen,\t ?T.\t ?J.\t ?Chen,\t ?and\t ?T.\t ?M.\t ?Keaveny,\t ??High-resolution finite element models with tissue strength asymmetry accurately predict failure of trabecular\t ?bone,?\t ?Journal of Biomechanics, vol. 33, no. 12, pp. 1575?1583, Dec. 2000. 87 Appendix ?A ?? ?Energy ?balance ?during ?a ?drop ?tower ?test The idea of calculating the energy balance during in-vitro dynamic testing and investigating the relationship between the amount of work done on the bone and its ability to withstand fracture was developed by Dr. Benedikt Helgason at ETH. Part of this idea is presented in one of our earlier studies currently under review for publication [63]. We are motivated to look at energy because of the current difficulty identifying what constitutes a femur at risk. Most researchers studying hip fracture using numerical models have focused on determining the maximum force that can be applied on a particular femur before fracture. It is, however, very difficult to work out the amount of force that a femur has to withstand in-vivo during a fall, since this force varies greatly according to body weight, the amount of soft tissue, impact direction, etc. So, while it becomes possible with FE models to compare which bone is weaker or stronger relative to other bones, the problem remains that one cannot identify who is at risk of fracture and who is not. Looking at energies instead of forces simplifies the task of characterizing a fall from standing, as the total energy available during a fall is bounded by potential\t ?energy\t ?of\t ?a\t ?person?s\t ?body\t ?mass elevated his or her centre of gravity while standing, a figure that is easily estimated. The effects of pelvic compliance and soft-tissue damping ? which dissipate and absorb energy, hence lowering the work on the femur ? are also more easily calculated in terms of energies than with forces. In theory, working with energy would greatly simplify the calculation\t ?of\t ?a\t ?person?s\t ?risk\t ?of\t ?fracture. We hypothesize that there is a maximum amount of energy that a bone can absorb during impact, and that fracture occurs when more energy than this maximum is delivered. We can investigate this energy connection because of the feasibility of calculating the amount of work applied on the femur during an impact test in our fall simulator apparatus (described in section 3.2.1 and by Gilchrist [41]). Numerically, we can also easily calculate the deformation energy in a simulated bone using finite element analysis. During a mechanical test in our fall simulator, at any time t after the initial impactor contact but before the femur fractures or stops the falling mass, the energy in the system Etot is given by: ???? = ????????? + ???(??)?? +????? +???? (A-1) where: ff(uf): force-displacement relationship at the greater trochanter Appendix A ? Energy balance during a drop tower test 88 fs(us): force displacement-relationship at the femoral head support uf: femur deformation us: deformation of the femoral head support m: the dropped mass v: instantaneous speed of the greater trochanter ??????????? ? Wres: residual energy such as the kinetic energy gained by the femur, heat, noise, deformation of the drop tower, damping and friction. uf, us, ff, fs, v, and Wres all vary with time t. It is also possible to calculate Etot before the impact, as this is simply the kinetic and potential energies carried by the falling mass: ???? =?????? +????? + ??? (A-2) where: vi: speed of the impactor right before impact g: acceleration due to gravity (9.81 m/s2) Note that uf and us are functions of time t, as the potential energy is here calculated using a moving datum. Using this definition, Etot is thus also a function of time but always equal to the expression in equation A-1. Theoretically, one could use a simpler energy expression for Etot by calculating the potential energy of the drop mass at its maximum height, before it is dropped; however, we found that a large portion of this potential energy is lost before impact due to friction in the drop tower. Using load cells at the greater trochanter and head support, together with tracked displacement data from high-speed video, we have experimental measurements for all variables in equations A-1 and A-2, except for Wres, which we assume to be small. We can use these data to investigate whether there is a specimen-specific value for ????????? that separates fracture and non-fracture events. If this turned out to be the case, then predicting a fracture using FE models would be a straightforward work calculation. However, we must first determine the accuracy of FE-based predictions of the amount of work done on a femur during a fall; hence the emphasis placed on energy results in chapter 3 of this thesis. Appendix B ? Maximum principal strains from dynamic model 89 Appendix ?B ?? ?Maximum ?principal ?strains ?from ?dynamic ?model Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1167L H1168R H1365R Appendix B ? Maximum principal strains from dynamic model 90 Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1366R H1368R H1369L H1372R Appendix B ? Maximum principal strains from dynamic model 91 Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1373R H1374R H1375L H1376L Appendix B ? Maximum principal strains from dynamic model 92 Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1377R H1380R H1381R Appendix B ? Maximum principal strains from dynamic model 93 Bone Experimental Fractures FE Principal Strains Anterior View Posterior View Anterior View Posterior View H1382L "@en . "Thesis/Dissertation"@en . "2014-05"@en . "10.14288/1.0166898"@en . "eng"@en . "Mechanical Engineering"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "Attribution-NonCommercial-ShareAlike 2.5 Canada"@en . "http://creativecommons.org/licenses/by-nc-sa/2.5/ca/"@en . "Graduate"@en . "A novel approach to finite element analysis of hip fractures due to sideways falls"@en . "Text"@en . "http://hdl.handle.net/2429/46324"@en .