Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 49, 3, pp. 599-619, Warsaw 2011 NUMERICAL MODELLING OF IMPACT FRACTURE OF CORTICAL BONE TISSUE USING X-FEM Adel A. Abdel-Wahab Vadim V. Silberschmidt Loughborough University, Wolfson School of Mechanical and Manufacturing Engineering, Loughborough, Leicestershire, UK e-mail: a.a.mohamed@lboro.ac.uk; v.silberschmidt@lboro.ac.uk A cortical bone tissue is susceptible to fracture that can be caused by events, such as traumatic falls, sports injuries and traffic accidents. A proper treatment of bones and prevention of their fracture can be sup- portedby in-depthunderstanding of deformation and fracture behaviour of this tissue in such dynamic events. Parameters such as damage initia- tion under impact, damage progression and impact strength can help to achieve this goal. In this paper, Extended Finite-Element Method (X- FEM) implemented into the commercial finite-element software Abaqus is used to simulate the actual crack initiation and growth in a cantilever beam of cortical bone exposed to quasi-static and impact loading using the Izod loading scheme. Izod tests were performed on notched bone specimens of bovine femur tomeasure its impact strength and to valida- te simulations. The simulation results show a good agreement with the experimental data. Key words: cortical bone, impact, X-FEM, finite-element, fracture 1. Introduction Bone is the principal structural component of a skeleton: it assists the load- bearing framework of a living body. Bone fractures have significant health, economic and social consequences. Both healthy and unhealthy bones are su- sceptible to fracture as a result of low- or high-energy trauma. High-energy traumas are usually linked to car or cycling accidents, while low-energy trau- mas often occur in contact sports. In both cases they are caused by dynamic loading. Factors such as mass, material properties and geometry of bone as well as the magnitude and orientation of applied loads affect its response to 600 A.A. Abdel-Wahab, V.V. Silberschmidt such loading. Bones are fracturedwhen they are exposed to severe loads that, in their turn, generate stresses exceeding its ultimate strength. Thus, a frac- ture event occurs initially at the material level that eventually affects the load-carrying capacity of thewhole bone at its structural level (Cullinane and Einhorn, 2002). It is worth mentioning here that bone is a viscoelastic mate- rial. Therefore, this type of mechanical behaviour has to be considered when dealing with dynamic events, such as impact. By developing adequate numerical models capable of predicting and de- scribing bone deformation and fracture behaviour, a detailed study of reasons for, and ways to prevent bone fracture could be performed. To plan preven- tion therapies and treatment strategies, scientific knowledge of bone fracture mechanisms is also needed. From the experimental point of view, the impact resistance of the cortical bone tissue did not get enough attention in the lite- rature though it is essential for activities such as jumping and running as well as protecting internal organs from impact as in the case of skull and ribs (Lee et al., 2011). Still, numerous previous studies have been devoted to analysis of quasi-static mechanical properties and resistance to fracture of the cortical bone tissue. For instance, Augat and Schorlemmer (2006) demonstrated the role of structural properties of cortical bone and itsmicrostructure in its com- petence. Another study, byBonney et al. (2011), was devoted to investigation of local variations in mechanical properties of cortical bone (porcine femur). Also, various experimental studies investigating the effect of structural pro- perties of cortical bone and its mechanical properties have been conducted (Zioupos, 1998; Zioupos et al., 1999;Wachter et al., 2002; Currey, 2004; Kulin et al., 2010). They dealt with acquisition of the respective data at different levels of the bonehierarchical structure –macroscopic andmicroscopic – using differentmethodological approaches. From the fracturemechanics perspective, a review of the structure and properties of bone focusing on mechanical and deformation behaviour at different length scales was introduced by Launey et al. (2010). In another study, Nalla et al. (2005) analysed the nature of local cracking events that preceded catastrophic fracture of human cortical bone and their relation to the microstructure. Zioupos et al. (2008) studied a re- lation between the strain rate and the microcracking damage in the fracture process of humancortical bone in tensile failure.According todynamicproper- ties of this tissue, only few studies paid attention to that issue. For instance, bothdynamic and staticmaterial properties of human femurwere investigated using, respectively, a split Hopkinson bar technique and tests with a universal testing machine (Katsamanis and Raftopoulos 1990). The average dynamic Youngmodulus of 19.9GPa was found to be 23% greater than that for static loading – 16.2GPa. In terms of bone impact characteristics, only preliminary Numerical modelling of impact fracture... 601 data are available (Panagiotopoulos et al., 2005), with a Charpy impact test used to measure the energy absorbed by strips cut from proximal femur. In a related experimental work, employing an Izod impact tester, Kovan (2008) investigated the absorbed energy and the impact strength of a mandible at different positions. In our previous work, the impact properties of a cortical bone tissue were investigated using Izod tests for different cortex position of bovine femur (Abdel-Wahab et al., 2010). The obtained experimental results emphasised that bovine femur cortical bone had a nearly uniform fracture energy character with regard to the cortex position. In addition, a 2D finite- element model to simulate the test and capture its behaviour up to fracture was developed. In that model, the behaviour of three different constitutive material models – linear-elastic, elastic-plastic and viscoelastic – based on our previous experiments (Abdel-Wahab et al., 2011) were compared. Itwas found that the viscoelastic model showed a good agreement with the experimental data. In other experiments, Lee et al. (2011) tested non-mineralized and mine- ralized materials, such as cortical bone utilizing a drop-weight test to investi- gate the impact strength along with the impact damage. In a similar study, longitudinal human cortical specimens were tested in a tensile impact tester at a strain rate of 133s−1 (Saha and Hayes, 1976). A marked non-linearity was observed in the stress-strain behaviour, including plastic deformation and strain-hardening effects. Themean tensile impact strength and impact energy were 126.3±33.1MPa and 18790±7355J/m2, respectively. Often, in-vivo bone fracture is initiated and promoted by cracks; therefo- re, fracture mechanics is used as an important tool to assess its strength and fracture toughness and to improve the diagnoses and treatment of bone fractu- res (Wang and Puram, 2004). Up to now, Linear-Elastic Fracture Mechanics (LEFM) was mostly used to assess the toughness of cortical bone tissue; it yields a single-valued fracture parameter – the critical stress intensity factor or the critical strain energy release rate (Norman et al., 1995, 1996). Due to complex composition andmicrostructure of the cortical bone tissue, it has se- veral tougheningmechanisms, such as diffuse microcracking, crack deflection, and fibre bridging (Zioupos, 1998; Vashishth et al., 2000; Yeni and Norman, 2000; Nalla et al., 2004). The inadequacy of LEFM theory application to de- scribe cortical bone tissue fracturewas raised due to observed resistance curve (R-curve)behaviour (Lucksanasombool et al., 2001;Ural andVashishth, 2006). Therefore, cohesive zone models were used to analyse the initiation and pro- pagation of cortical bone cracks (Ural andVashishth, 2006; Yang et al., 2006; Cox andYang, 2007). In a recent study,Morais et al. (2010) demonstrated the 602 A.A. Abdel-Wahab, V.V. Silberschmidt adequacy of a Double-Cantilever Beam (DCB) test for determining fracture toughness under puremode-I loading of cortical bone by implementing a new data-reduction scheme based on specimen compliance. Despite this body of research, experimental and numerical studies of the dynamicbehaviour of a cortical bone tissue attracted less attention.Therefore, this study comprises two parts covering experimental and numerical aspects of such analysis. Still, analysis of the actual crack initiation and growth was hard to achieve using thementioned approaches in simulations.With the new Extended Finite-ElementMethod (X-FEM), crack initiation and propagation can be modelled more easily. Thus, the aim of this study is to develop and validate a numerical model for analysis of fracture behaviour of the cortical bone tissue under impact and quasi-static loading using X-FEM. 2. Materials and methods 2.1. Specimen preparation Longitudinal specimens for Izod tests were cut from fresh bovine femora bones (aged 1.5-2 years) along their osteons, see Fig.1a. Sixteen specimens were used to assess reproducibility of experimental results. All the specimens had the same dimensions: 50mm×8mm×4mm (length×width×thickness), see Fig.1b. A 300µm-deep notch was created perpendicular to the bone axis and along the radial axis using a razor. Specimens were stored at room tem- perature in a 0.9% saline solution until tested. Fig. 1. (a) Cortical bone with axes and direction of specimen cutting; (b) izod test specimen Numerical modelling of impact fracture... 603 2.2. Izod test An amount of energy absorbed during a suddenly applied force can be qu- antified with impact tests of materials; the required data about deformation and failure of materials during high-strain-rate loading cannot be determi- ned using quasi-static fracture tests. Impact testing is usually performedwith Charpy or Izod test machines. An Izod test system incorporates a swinging pendulum (hammer) that impacts a notched specimen fixed in a cantilever- beam position with the notch facing the hammer (Lee et al., 2011). In our study, impact tests were carried out using a CEAST Resil Impactor. In the tests, the bottomhalf of the specimenwas fixedfirmly in themachine vice and aknife-edgewedgewas used to define thenotch position.Theupperhalf of the specimenwas struck by a pendulumhammerwith a controlled level of energy. The distance between the notch and the position of hammer strike was stan- dard – 22mm. In this study, a calibrated hammer with a mass of 0.6746kg and 0.3268m long was used. The maximum nominal hammer energy of 2J corresponds to the striking position of 150◦ resulting in an impact velocity of 3.46m/s. The level of initial energy can be varied by changing the initial angle of the hammer. The energy level used in this studywas 0.5J that corre- sponded to the initial angle of 58◦. A piezoelectric force transducer was fixed rigidly to the hammer to capture the impact-force signal.When the pendulum is released from the pre-defined angle, its impact with the specimen generates a change in electrical resistance of the piezoelectric sensor that is captured by the data acquisition system –DAS 8000 – connected to the impactor. 2.3. Numerical model 2.3.1. Extended finite element method X-FEM is a numerical method that enables analysis of crack propagation without remeshing a cracked specimen in accordance to newly created crack faces. It employs a local enrichment of the approximation areas. Themethod can be useful for evolving processes with non-smooth characteristics in small parts of a computational domain, e.g. near discontinuity or singularity regions, as in the case of cracks forwhich the standard finite elementmethod is not ac- curate. TheX-FEMwas first introduced byBelytschko andBlack (1999). The enrichment is realised based on the partition-of-unity concept developed by Melenk and Babuska (1996); it allows incorporation of local enrichment func- tions intoafinite-element approximationdomain.Spatial enrichment functions with additional degrees of freedom ensure account for discontinuities. The ge- neral framework of thismethod is incorporated in the finite- element software Abaqus 6.10 (2010). 604 A.A. Abdel-Wahab, V.V. Silberschmidt 2.3.2. Geometry, mesh, and boundary conditions In this study, three finite-elementmodels (FEM)were developed:ModelA, Model B, andModel C.Model A is a 2DX-FEMmodel used to simulate the fracture of cortical bone exposed to impact loading in the Izod test, see Fig.2. Model B is a 3D formulation of Model A, whereas Model C is a 3D X-FEM model for quasi-static fracture analysis. The impact tests were simulated with the finite-element software Abaqus 6.10/Implicit using Models A and B to verify the applicability of the X-FEM to analysis of the failure behaviour of the cortical bone tissue under impact-loading conditions. In addition,ModelC was developed to elucidate fracture development in the cortical bone tissue under different loading conditions. A full description ofModel A can be found in our previous work (Abdel-Wahab et al., 2010); however, in that work only the behaviour up to failurewas studied. Figure 2 shows its geometry andmesh formulation. Fig. 2. (a) Real hammer; (b)Model A; (c) hammer-specimen interaction andmesh around the notch In Model B, the real geometry and masses of the hammer and specimen with a 300µm notch were used (see Fig.3). In Abaqus 6.10, a kinematic co- upling constraint is used to transmit rotation to the structurewhile permitting radialmotion.Hence, this featurewas employed toconstraint thehammer from radial or translational movements except around one axis only; it is z-axis in our case. To get the exact movement of the hammer, as it happens in reality, Numerical modelling of impact fracture... 605 all the nodes of the surface of the inner cylinder of the upper block of the hammerwere kinematically coupled to a reference point at themiddle of that cylinder, then the reference pointwas restrained to translate along x, y, and z and to rotate around x or y axes. The reference point was only free to rotate around the z-axis. Two sets were defined: the reference point (set1) and the inner surface (set2). The following equation was used to define the coupling between those sets ur1−ur2 =0 (2.1) where ur1 and ur2 are the rotational degrees of freedom for set1 and set2, respectively. Fig. 3. (a) Setup of Izod test; (b) Model B; (c) hammer-specimen interaction; (d) Meshing of hammer and specimen (specimen is increased) On the other hand, two surfaces were chosen to define a surface-to-surface contact between the specimen and the hammer. These surfaces are shown as S1 and S2 inFig.3d.Themaster surfacewas chosen tobe S2with S1 chosen to be the slave one. Themesh of themaster surfacewas adjusted to be finer than that of the slave surface. A finite slidingwith frictionless tangential behaviour formulation was chosen between the two surfaces. To reduce the computation time, the hammer was assembled very close to the specimen as the initial position in simulations. An angular velocity of 5.33rad/s – around the x-axis – corresponding to the initial angle of 58◦ (initial energy of 0.5J) was applied 606 A.A. Abdel-Wahab, V.V. Silberschmidt to the hammer.The support of the specimenwasmodelled as rigid; all degrees of freedomof the specimenbottompartwere constrained (seeFigs.3b). Linear tetrahedron (C3D4) elements were used for both specimen and hammer; it is currently the only element type that can be used for 3D X-FEM analysis. A total number of 53113 elements and 10967 nodes for the bone specimen and 23695 elements and 6871 nodes for the hammer were used. The force due to contact pressure between the piezoelectric force sensor and counterpart of the specimenwas recorded in the history output of the finite-element software Abaqus/Implicit.Also, the status of theX-FEMthat shows the crackpathwas used as an output alongwith the distributions of stress and strain components and their principal values. In order to compare the fracture behaviour of the cortical bone specimen under quasi-static and impact loading,Model Cwas developed, Fig.4. It con- sists of a cantilever-beam specimen of cortical bone with the same geometry, mesh andmaterial properties asModelB. InModelC, the hammerwas exclu- ded fromthe analysis and, instead, adisplacement-controlled load of 2mmwas applied at the same position of the hammer-specimen interaction, see Fig.4b. Fig. 4. (a)Meshed 3D quasi-static specimen; (b) applied displacement and boundary conditions of 3D quasi-static model 2.3.3. Material properties Elastic material properties for the hammer and cortical bone tissue used in the numerical simulations are given in Table 1. The viscous behaviour of bones was introduced into models A and B in terms of the Prony series Numerical modelling of impact fracture... 607 expansion based on the normalized creep compliance. These material con- stants are the normalized shear modulus, g =0.13256±0.01, and relaxation time, τ = 119.57± 0.5. These relaxation parameters were defined based on our experimental creepdata andusing the evaluation technique for viscoelastic materials implemented in Abaqus 6.10, as discussed elsewhere (Abdel-Wahab et al., 2011). This technique converts the creep data into relaxation data by means of convolution integrals; then it uses thenormalized shearmodulusdata to determine the Prony series parameters (Abaqus 6.10). Only a linear-elastic material model was used forModel C for quasi-static loading.Allmaterial properties for cortical bonewere obtained in our previous experiments (Abdel-Wahab et al., 2011). Table 1.Material properties for finite-element models Part Material Elastic Poisson’s Density modulus [GPa] ratio [–] [kg/m3] Hammer Carbon steel 210 0.3 7850 Specimen Bone 23.15±0.72 0.44±0.009 1860±0.9 In this study, the X-FEM-based cohesive segments methodwas used to si- mulate crack initiation andpropagationalonganarbitrary, solution-dependent path in thebalkmaterial, since thecrackpath is independentof theboundaries of the elements in the mesh (Abaqus 6.10, 2010). In the model, the enrich- ment area was chosen as the bone specimen; the crack was introduced as a planewith dimensions of 300µm×4mm to reproduce the notch depth and the specimen thickness of real experiments, respectively. Damagemodelling allows simulation of crack initiation and eventual failure of the enriched area in the solution domain. The initial response is linear, while the failure mechanism consists of a damage initiation criterion and a damage propagation law. The damage initiationwasdefinedbasedon themaximumprincipal strain of 0.25% (Pattin et al., 1996; Bayraktae et al., 2004). When the damage initiation cri- terion is met, the damage propagation law starts to take place. In this study, the damage evolution is defined in terms of fracture energy (per unit area), and a linear softeningwas chosen. Themixed-mode behaviourwas chosen and the fracture energies for those modes were introduced into X-FEM. The frac- ture toughness values were 1374N/m, 4710N/m and 4016N/m for Mode I, Mode II and Mode III, respectively (Feng et al., 2000). The following model assumptions were made: (1) homogeneous and isotropic material properties for both the specimen and the hammer; (2) frictionless contact between them. 608 A.A. Abdel-Wahab, V.V. Silberschmidt 3. Results and discussion The fracture behaviour of bone is closely coupled to its underlyinghierarchical structure; therefore, the measured fracture parameters depend on the length scale at which they are measured. Moreover, in order to correctly assess the cracking and fracture behaviour of bone, the crack direction should be clini- cally relevant, i.e. cracks propagating in the transverse direction and involve realistic flaw sizes (Koester et al., 2008). Accordingly, in the experiment part of this study, two groups of specimens were tested: the first group consists of specimens with 300µmnotch and second group consists of specimens with 600µm notch. For all the specimens the notch was generated perpendicular to the osteons, see Fig.1a. The chosen notch sizes were below 600µm accor- ding to the physiological pertained flaw size reported byKoester et al. (2008). As our quasi-static tests on cortical bone demonstrated (Abdel-Wahab et al., 2011) and as well-known from the literature, the cortical bone properties va- ry from one cortex position to another due to variation in both composition and microstructure. Thus, impact strength was measured for specimens cut from different cortex positions called anterior, posterior, medial and lateral, see Fig.1a. Figure 5 shows a comparison for data obtained for these positions for 300µm- and 600µm-notched specimens; the initial energy level used in those tests was 0.5J. The impact strength was measured as the absorbed im- pact energy divided by the un-notched cross sectional area of the specimen. The obtained results show that for the same applied energy level and cortex position, the notch size has a negative effect on impact strength of the cortical bone tissue. In general – apart from the medial cortex position with nearly the same average impact strength for both notch sizes – the specimens not- ched with 300µm required higher energy per unit area to fail compared to those with 600µm notch. Statistically, the average impact strength required to produce fracture appears to be higher at the lateral position, with diffe- rent magnitudes for different notch sizes. However, considering the spread of the mean impact energy for all cortex positions and notch depths, it is appa- rent that it is within the interval from 6kJ/m2 to 12kJ/m2. The statistical analysis for cortex positions revealed no significant difference of themean im- pact strength for both notch depths for 300µm (p = 0.862) and for 600µm (p=0.354). Also, checking the combined effect of both factors – cortex posi- tion and notch depth – on themean impact strength no significant difference was demonstrated (p = 0.642). Based on these results, bovine femur seems to have a nearly uniform impact strength. The obtained values for the im- pact strength are in a good agreement with other studies in the literature; for instance, using four-point bending setup, bovine cortical boneswere found Numerical modelling of impact fracture... 609 to have an impact strength in the range of approximately 6-24kJ/m2 (Reilly andCurrey, 2000). In a separate study, the impact strength of bovine cortical bone was found to be 10kJ/m2 in drop-weight-tower tests (Lee et al., 2010). The obtained impact strength results of our study agree well with those in the literature, even based on different techniques. Fig. 5. Impact strength of longitudinal cortical bone specimens for different cortex positions and two different notch sizes Below, results of simulationswithModelsA,B andCof the fracture beha- viour of cortical bone tissue under impact and quasi-static loading are presen- ted in comparison with experimental data. The contact-force profile obtained in our Izod tests of cortical-bone specimens was used to validate the develo- ped finite-element models A andB. A comparison of the experimental results with simulation ones (Fig.6) demonstrates thatModel A reproduces transient fracture behaviour of the cortical bone tissue. However,Model B, though sho- wing a good agreement both withModel A and the experimental results until its termination point, results in an unrealistic fracture scenario afterwards as will be discussed below. Some deviations from the experimental results can be linked to different factors that are not incorporated into the current stage of model development. Among them, there is a complex hierarchical structure, anisotropy, and heterogeneity of the cortical bone tissue contributing to its fracture behaviour. Undoubtedly, as a consequence of composition andmicro- structure, there are several experimentally observed toughening mechanisms in the fracture process of cortical bone tissue, such as diffuse microcracking, crack deflection andfibrebridging (Zioupos, 1998;Vashishth et al., 2000;Nalla et al., 2004) that can affect propagation of the crack. It was also noticed that the cortical bone specimen failed in tests catastro- phically in a brittle manner as soon as the maximum force was reached; it is 610 A.A. Abdel-Wahab, V.V. Silberschmidt Fig. 6. Comparison of evolution of contact force in impact loading (notch size 300µm) Fig. 7. Evolution of contact force and normalized crack length inModel A (notch size 300µm) presented by a nearly vertical line in Fig.6 after the peak. In the simulation results of Model A, see Fig.7, specimen fracture is represented in terms of a normalized crack length (Lcr/(W−Lnotch)) against time. It was calculated as the length of the crack (Lcr) divided by the un-notched width (W −Lnotch) measured on the specimen surface. At the beginning, the crack started to evo- lve slowly up to approximately 10% of this width at t = 0.2ms, then grew steeply up to the contact force peak (331.3N at t = 0.36ms) at which only 30% of the specimen failed. However, even when the force started to decline; it was still high enough to propagate the crack. The crack growth accelera- Numerical modelling of impact fracture... 611 ted after the maximum force position up to the specimen failure point, with the crack spending only 0.28ms to reach the opposite side of the specimen. Obviously, at the moment of complete failure of the specimen, the contact force vanished. For Model A, the STATUSXFEM output – available in Abaqus 6.10 – shows the crack evolution during the course of analysis. This parameter varies over the range between 0 and 1; when it equals to 0 there is no damage and in the case of complete failure, it equals 1. Figure 8 shows the evolution of the crack originating from the notch and propagating across the width of the specimen towards the opposite side. Fig. 8. Crack evolution at different time increments (Model A) It was noticed that the crack started to propagate immediately along an inclined plane, with the notch root indicatingmixed-mode fracture behaviour (Mode I and Mode II). That was followed by a small horizontal crack path. These changes in the crack path direction correspond to kinks in the first part of the contact force-time curve up to the contact-force peak, see Fig.7. Obvio- usly, from the point of view of the beam theory, the studied cantilever beam is exposed to two types of stresses: normal stress (bending) and transverse shear stress. In thesemodels, the beam span-to-width ratio is 2.2; so the transverse shear stresses should be significant. It was found that the shear stress level was comparable to that of normal ones (Figs.9 and 10). Also, a transverse distribution of shear stress is parabolic across the width of the cantilever and uniform from the position of specimen-hammer interaction up to the notch location, see Fig.9. This stress state causes Mode II fracture. On the other hand, normal stress S22 in the y-direction causes Mode I fracture behaviour. Therefore, immediately after the impact takes place and at the notch root, only Mode I fracture took place due to vanishing (or very small) shear stres- ses. At the neutral plane of the cantilever, where the shear stress has its peak 612 A.A. Abdel-Wahab, V.V. Silberschmidt while the bending stress vanishes, Mode II dominates fracture. The opposite side of the beam was under compression, with shear stress vanishing. Thus, between the notch side and/or the opposite side and the neutral plane, both Modes I and II took place. When the crack started to propagate, it caused stresses redistribution resulting in a complex crack path; this can be seen in Fig.11. Fig. 9. Distribution of shear stress S12 inModel A (a),Model B (b) andModel C (c) Fig. 10. Distribution of normal stress S22 in Model A (a), Model B (b) and Model C (c) While bothModels A and C showed one main crack that originated from the notch and propagated in the depth andwidth of the specimen towards the opposite face (Figs.11a and 11c), Model B demonstrated a band of random short cracks around the notched area of the specimen with no major crack percolating the specimen. In this model, due to the impact load suddenly applied to the specimen, elastic stress-waves were generated, propagating in Numerical modelling of impact fracture... 613 different directions in the specimenactivating all fracturemodes and initiating multiple cracks. As elastic waves keep moving and reflecting in a complex way inside the specimen, a band of damaged zone was formed instead of a single crack demonstrating limitations of the current 3D X-FEM routine of Abaqus 6.10 with respect to dynamic cracking problems (Fig.11b). On the otherhand, inModelC, a single crackpropagated fromthe root of thenotch to theopposite side of the specimen.That crackpathwas similar to that observed in the experiment, see Fig.11d. Still, the effect of underlying heterogeneous microstructure, such as pores andweak interfaces – cement lines – between the osteons and interstitial matrix that can deflect the crack (Nalla et al., 2004), can be responsible for some deviations from the solution obtained with the used isotropic homogeneous formulation. Fig. 11. Distributions of maximum principal strain inModel A (a), Model B (b) and Model C (c). (d) Final crack path in Izod-test specimen The damage initiation and evolution behaviour for Models A and C are shown in Fig.12. In this figure, the linear component of displacement in the x-direction at the contact position is obtained from Abaqus and used as an external parameter for the axis of abscissas. Apparently, damage started im- mediately as the beam was subjected to the load. In addition, at the same deformation level, the crack length in Model C was larger than that in Mo- del A. It is worth recalling that a linear elastic material model was assigned to the quasi-static model (Model C), whereas a viscoelastic material model was used in the impact models (Models A and B). Hence, different behaviour could be due to activated relaxation mechanisms that assisted to damp some of the applied energy of the hammer and resulted in differences between two results. However, at the beginning up to 0.15mmboth cases are close; at this stage bothmodels still behave elastically. Apparently, the current 3DX-FEM routine inAbaqus 6.10 causes some convergence problems even in quasi-static 614 A.A. Abdel-Wahab, V.V. Silberschmidt formulations: simulations with Model C terminated before the crack reached the opposite side (at the normalized length of about 0.75 (Fig.12). Model A did not demonstrate such a problem. By investigating a through-thickness va- riation in the crack length forModel C, its front demonstrated a non-uniform character (Fig.13) changing with displacement. Due to the stress state vary- ing across the thickness, the crack propagates with a different rate at various positions along its front. Fig. 12. Evolutions of crack length forModels A and C Fig. 13. Crack length variation along its front inModel C Numerical modelling of impact fracture... 615 A direct measurement of crack length as a function of time and/or defor- mation is not available in Abaqus 6.10; therefore various images at different time increments were taken and measured using Image Pro-Express softwa- re (Image-Pro Express 2005). The crack propagation rate for Model A was obtained by differentiation of the curve-fit equation for the crack length. Figu- re 14 shows both parameters as a function of time. Since the crack length has a quadratic curve-fit equation, the crack growth rate shows a linear increase with time. At t=0, the specimen was exposed to a sudden impact with the hammer with the initial velocity of 1.74m/s, causing an initial crack propa- gating rate of 2.054m/s. This rate has evolved during the fracture process of the specimen to reach 19.5m/s at the final percolation of the specimen. Fig. 14. Evolution of crack length and crack propagation rate inModel A 4. Conclusions Crack initiation and growth under quasi-static and impact loading of cortical bone tissue were studied using experimental and numerical simulations. Izod tests were performed to characterise its impact strength for different cortex positions. Three different finite-element models –Models A, B, and C – were implemented into the commercial finite-element software Abaqus 6.10 using its implicit solver. A series of simulationswas performed to study the crack in- itiation and propagation under quasi-static and impact loading. The obtained numerical results were quite close to the experimental ones, and the numeri- cal models have the capability to reproduce the failure of cortical bone tissue under both impact and quasi-static loading. The finite-element results provi- 616 A.A. Abdel-Wahab, V.V. Silberschmidt de more detailed information than the experimental tests and helped to gain a better understanding of the fracture behaviour of the cortical bone tissue. Numerical simulations showed that its fracture behaviour was reasonably well predicted using Model A in terms of the contact force profile and the crack path, while Model B exhibited unrealistic fracture scenario: formation of a damage band with multiple cracks across the specimen width and thickness around thenotch area.The stress state generatedby the applied load triggered Modes I and II in the fracture process of the bone specimen in the Izod test setup. In general, the results showed the suitability of the developed numerical approach to study the fracture of cortical bone tissue under quasi-static and impact loading. References 1. Abaqus 6.10, 2010,Analysis user’s manual, Section 10.6.1 2. Abdel-Wahab A.A., Maligno A., Silberschmidt V.V., 2010, Dynamic properties of cortical bone tissue: Izod tests and numerical study, Comput. Mater. Continua, 19, 3, 217-238 3. Abdel-WahabA.A.,AlamK., SilberschmidtV.V., 2011,Analysis of ani- sotropic viscoelastoplastic properties of cortical bone tissues, J. Mech. Behav. Biomed. Mater., DOI.org/10.1016/j.jmbbm.2010.10.001 4. Augat P., Schorlemmer S., 2006, The role of cortical bone and its micro- structure in bone strength,Age Ageing, 35-s2, ii27-ii31 5. Bayraktae H.H., Morgan E.F., Niebur G.L., Morris G.E., Wong E.K., Keaveny T.M., 2004, Comparison of the the elastic and yield pro- perties of human femoral trabecular and cortical bone tissue, J. Biomech., 37, 27-35 6. BelytschkoT., BlackT., 1999,Elastic crack growth in finite elementswith minimal remeshing, Int. J. Numer. Meth. Eng., 45, 5, 601-620 7. BonneyH.,ColstonB.J.,GoodmanA.M., 2011,Regional variation in the mechanical properties of cortical bone from the porcine femur, Med. Eng. and Phys., 33, 4, 513-520 8. Cox B.N., Yang Q., 2007, Cohesive zonemodels of localization and fracture in bone,Eng. Fract. Mech., 74, 1079-1092 9. CullinaneD.M.,EinhornT.A., 2002,Biomechanicsofbone, [In:]Principles of Bone Biology, Bilezikian J.P., Raisz L.G., Rodan A.R. (Edit.), San Diego, USA, Academic Press, 1 Numerical modelling of impact fracture... 617 10. Currey J.D., 2004, Tensile yield in compact bone is determined by strain, post yield behaviour bymineral content, J. Biomech., 37, 549-556 11. Feng Z., Rho J., Ziv I., 2000,Orientation and loading condition dependence of fracture toughness in cortical bone,Mater. Sci. Eng. C, 11, 41-46 12. KatsamanisF.,RaftopoulosD.D., 1990,Determinationofmechanicalpro- perties of human femoral cortical bone by the Hopkinson bar stress technique, J. Biomech., 23, 11, 1173-1184 13. Koester K.J., Ager III J.W., Ritchie R.O., 2008, The true toughness of human cortical bone measured with realistically short cracks, Nat. Mater., 7, 672-677 14. Kovan V., 2008, An assessment of impact strength of the mandible, J. Bio- mech., 41, 16, 3488-3491 15. Kulin R.M., Jiang F., Vecchio K.S., 2010, Effects of age and loading rate on equine cortical bone failure, J. Mech. Behav. Biomed. Mater., 4, 1, 57-75 16. Launey M.E., Buehler M.J., Ritchie R.O., 2010, On the mechanistic origins of toughness in bone,Ann. Rev. Mater. Res., 40, 25-53 17. Lee S., Novitskaya E.E., Reynante B., Vasquez J., Urbaniak R., Ta- kahashiT.,WoolleyE., TombolatoL., ChenPo-Yu,McKittrick J., 2011, Impact testing of structural biologicalmaterials,Mater. Sci. Eng. C, 31, 4, 730-739 18. Lucksanasombool P., Higgs W.A.J., Higgs R.J.E.D., Swain M.V., 2001, Fracture toughness of bovine bone: influence of orientation and stora- ge media,Biomater., 22, 3127-3132 19. Melenk J. M., Babuska I., 1996, The partition of unity finite element me- thod: basic theory and applications,Comput. Method Appl. Math., 39, 289-314 20. Morais J.J.L., de Moura M.F.S.F., Pereira F.A.M., Xavier J., Do- uradoN.,DiasM.I.R.,AzevedoJ.M.T., 2010,Thedouble cantilever beam test applied tomode I fracture characterizationof cortical bone tissue, J.Mech. Behav. Biomed. Mater., 3, 6, 446-453 21. Nalla R.K., Kruzic J.J., Ritchie R.O., 2004, On the origin of the tough- ness of mineralized tissue: microcracking or crack bridging?Bone, 34, 790-798 22. Nalla R.K., Stolken J.S., Kinney J.H., Ritchie R.O., 2005, Fracture in human cortical bone: local fracture criteria and toughening mechanisms, J. Biomech., 38, 1517-1525 23. Norman T.L., Nivargikar S.V., Burr D.B., 1996, Resistance to crack growth in human cortical bone is greater in shear than in tension, J. Biomech., 29, 1023-1031 24. Norman T.L., Vashishth D., Burr D.B., 1995, Fracture toughness of hu- man bone under tension, J. Biomech., 28, 309-320 618 A.A. Abdel-Wahab, V.V. Silberschmidt 25. Panagiotopoulos E., Kostopoulos V., Tsantzalis S., Fortis A.P., Doulalas A., 2005, Impact energy absorption by specimens from the upper end of the human femur, Injury, 36, 5, 613-617 26. PattinC.A., CaletW.E., CarterD.R., 1996,Cyclicmechanical property degradation during fatigue loading of cortical bone, J. Biomech., 29, 69-79 27. Reilly G.C., Currey J.D., 2000, The effect of damage and microcracking on the impact strength of bone, J. Biomech., 33, 337-343 28. Saha S., Hayes W.C., 1976, Tensile impact properties of human compact bone, J. Biomech., 9, 4, 243-244, IN5, 245-251 29. UralA.,VashishthD., 2006,Cohesivefinite elementmodelling of agerelated toughness loss in human cortical bone, J. Biomech., 39, 2974-2982 30. V. Image-Pro Express, V.m.c., 2005 31. Vashishth D., Tanner K.E., Bonfield W., 2000, Contribution, develop- ment andmorphology of microcracking in cortical bone during crack propaga- tion, J. Biomech., 33, 1169-1174 32. Wachter N.J., Krischak G.D., Mentzel M., Sarkar M.R., Ebinger T.,Kinzl L., Claes L.,AugatP., 2002,Correlationof bonemineral density with strength and microstructural parameters of cortical bone in vitro, Bone, 3, 90-95 33. WangX., PuramS., 2004,The toughness of cortical bone and its relationship with age,Ann. Biomed. Eng., 32, 123-135 34. Yang Q.D., Cox B.N., Nalla R.K., Ritchie R.O., 2006, Re-evaluating the toughness of human cortical bone,Bone, 38, 878-887 35. Yeni Y.N., NormanT.L., 2000, Calculation of porosity and osteonal cement line effects on the effective fracture toughness of cortical bone in longitudinal crack growth, J. Biomed. Mater. Res. B Appl. Biomater., 51, 504-509 36. Zioupos P., 1998, Recent developments in the study of failure of solid bioma- terials and bone: fracture and prefracture toughness, Mater. Sci. Eng. C, 6, 33-40 37. Zioupos P., Currey J.D., 1998, Changes in the stiffness, strength, and to- ughness of human cortical bone with age,Bone, 22, 57-66 38. ZiouposP.,CurreyJ.D.,HamerA.J., 1999,The roleof collagen in thedec- lining mechanical properties of aging human cortical bone, J. Biomed. Mater. Res., 45, 108-6 39. ZiouposP., HansenU., Currey J.D., 2008,Microcrackingdamage and the fracture process in relation to strain rate in human cortical bone tensile failure, J. Biomech., 41, 14, 2932-2939 Numerical modelling of impact fracture... 619 Modelowanie numeryczne złamania uderzeniem istoty korowej kości za pomocą rozszerzonej metody elementów skończonych Streszczenie Istota korowa kości jest tkanką podatną na złamanie wywołanem.in. takimi zda- rzeniami jak upadki, kontuzje sportowe czy wypadki samochodowe. Odpowiednie leczenie czy profilaktyka stanu kości powinny być wsparte głębokim zrozumieniem zjawisk deformacji i pękania tej tkanki zachodzących podczas takich dynamicznych obciążeń. Znajomość parametrów opisujących inicjację uszkodzenia pod wpływem obciążenia uderzeniowego, progresja tego uszkodzenia i wreszcie odporność na zła- manie może być pomocna w tej praktyce. W prezentowanej pracy przedyskutowano implementację rozszerzonej metody elementów skończonych (X-FEM) do standardo- wego i komercyjnego pakietu ABAQUS w celu symulacji faktycznych zjawisk towa- rzyszących inicjacji i propagacji pęknięcia belki wykonanej z materiału oddającego właściwości istoty korowej poddanej działaniu quasi-statycznych i uderzeniowych ob- ciążeń,wspartej próbami udarnościowymi Izoda.Testy Izoda zostały przeprowadzone naprzygotowanychpróbkachzbydlęcej kości udowej z naciętymkarbemwcelu ekspe- rymentalnej weryfikacji wytrzymałości na uderzenie. Rezultaty badań numerycznych wykazały dobrą zgodność z wynikami tych prób. Manuscript received March 10, 2011; accepted for print April 29, 2011