Fracture in porous bone analysed with a numerical phase-field dynamical model

A dynamic phase-field fracture finite element model is applied to discretized high-resolution three-dimensional computed tomography images of human trabecular bone to analyse rapid bone fracture. The model is contrasted to quasi-static experimental results and a quasi-static phase-field finite element model. The experiment revealed complex stepwise crack evolution with multiple crack fronts, and crack arrests, as the global tensile displacement load was incrementally increased. The quasi-static phase-field fracture model captures the fractures in the experiment reasonably well, and the dynamic model converges towards the quasi-static model when mechanically loaded at low rates. At higher load rates, i.e., at larger impulses, inertia effects significantly contribute to an increased initial global stiffness, higher peak forces and a larger number of cracks spread over a larger volume. Since the fracture process clearly is different at large impulses compared to small impulses, it is concluded that dynamic fracture models are necessary when simulating rapid bone fracture.


Introduction
Rapidly growing fractures in low-quality porous (osteoporotic) bone are a major cause of long-term pain and physical disability, and have an enormous impact on the health care system (Braithwaite et al., 2003).However, to study bone fracture on a tissue level is difficult; its multiscale architecture and intricate material composition makes mechanical analyses of bone a severely complicated task (Sabet et al., 2016).Damage and fracture in bone are known to be dependent on the volume density of bone as well as the quality of bone tissue (cf.Keaveny et al., 2001;Dubey and Tomar, 2018).Many studies, using a combination of experimental observations and finite element (FE) analyses, have been made to characterize the mechanical behaviour of bone, e.g.Schileo et al. (2008), Schileo et al. (2014), Grassi et al. (2016), Wu et al. (2018Wu et al. ( , 2022)), and a wide range of elastic moduli have been reported.On a global scale, FE models of bone (e.g.femur heads) using varying stiffness and toughness distributions based on density estimations within the bone from computed tomography (CT) scans have recently been employed to e.g.evaluate deformation fields (cf.Shen et al., 2019;Gustafsson et al., 2021;Kok et al., 2022).
Most bone fractures are caused by impact upon falling, and the majority of hip fractures are reported as being caused by sideways fall from standing (cf.Parkkari et al., 1999;Iolascon et al., 2013;Helgason et al., 2014).Helgason et al. (2014) performed a dynamic 3D FE simulation of femur, modelled as a homogenous solid on a global scale.Focus was on the cortical surface, and the simulations showed high strain magnitudes in similar regions as where fractures were observed in experiments.Varga et al. (2016) compared dynamic drop tower experiments on a femur head to a quasi-static non-linear 2D FE model to analyse strength, stiffness and fracture sites.More recent studies, both experimental and numerical, have addressed the issue of dynamical fracture of femurs (cf.Askarinejad et al., 2019;Jazinizadeh et al., 2020;Kok et al., 2021).However, like most previous studies, these recent analyses were conducted on a global scale, using a "smeared", or averaged, continuum representation of the bone tissue.There is, however, a need to capture the crack dynamics also on a lower scale, i.e., at the trabecular scale.In response to this need, a quasi-static 3D phase-field finite element fracture model was applied on tissue level for trabecular bone by Braesch-Andersen et al. (2022).The trabecular bone was extracted from a human femur head and the aim was to analyse trabecular failure under global compression at high resolution.The study concludes that the numerical model to a great extent could simulate stably growing fractures in the trabecular bone.
In this study, dynamic fracture in trabecular bone, cut from a human femur head, is analysed numerically, as well as in comparison with quasi-static experiments.Trabecular bone samples, with pre-fabricated notches, were loaded in tension while monitored in synchrotron-based X-ray CT to obtain 3D images of growing cracks.The crack experiments are analysed with respect to fracture evolution and compared qualitatively to numerical simulations at high-resolution using a phasefield approach to simulate brittle fracture.The simulations are performed with slow, quasi-static loads and with rapid, dynamic loads at different loading rates in order to investigate to what extent the mechanical response is altered when subjected to loading rates close to those expected in actual fracturing incidents.Only trabecular bone is considered, meaning the difference in density and elastic properties on a global scale is inherently obtained.Due to large variation and uncertainty in mechanical properties of bone, quantitative comparisons of loads and displacement are outside the scope of this study.

Experiment
Human femur heads, removed in hip arthroplasty, were collected (ethics committee approval EK-29/2007, Zürich) and stored in − 80 • C. Before experimental use the material was thawed overnight to room temperature.Two samples with rectangular cross-section 8 × 4 mm and height 20 mm were cut from the femur heads with a clinical oscillating saw.Between cutting and testing, the samples were kept in phosphate buffered saline to avoid dehydration.A blunt edge-notch of length 2.3 mm was cut into the side of each sample at its mid-height in order to localize fracture (illustrated in the schematics in Fig. 1).The CT scans were performed at level with the notch, Fig. 1.The scanned region has dimension 6 × 3 × 3 mm and contains the "tip-part" of the blunt notch with length ~1.4 mm.This can cause boundary effects in the FEsimulation as the model geometry was discretized from the CT scan.Of particular importance for the numerical study is the length of the crack which is within the scanned region.In order to observe crack propagation over a distance, part of the notch was left outside the field of view.The discretized model thus contained a notch of length in the range of the native pore size in the trabecular bone (see Fig. 1).Boundary effects are further discussed in Section 4.
Displacement-controlled quasi-static tensile tests were conducted in a synchrotron-based X-ray CT (the Tomcat beamline at the Swiss Light Source, Paul Scherrer Institute, Villigen, Switzerland).The samples were clamped at the top and bottom, such that the free height was about 10 mm.The region around the notch was scanned by CT at zero loading, then at regular intervals, until the bone sample was completely broken.To avoid interference from movement between the loading and scanning, the tensile testing machine was paused during the scanning.The load rate was 38 μm/s and holding time was about 1 min.The voxel resolution was 3.25 μm; a compromise between a high resolution and a sufficiently large field of view, i.e., 6 mm in width and 3 mm in height.Specific parameters for the scanning were: 20 keV for beam energy, mA for ring current, two aluminium filters with thickness of 200 μm and 50 μm, 0.18 • for the angular rotation step, 200 flat field images and black current images for noise reduction and, finally, 30 ms for the exposure time of each radiograph.In total, 1231 sequential projections were acquired within approximately 50 s for each scan.The projections were reconstructed into cross-sectional images using the method of Schindelin et al. (2012) and cropped to roughly 2000 × 1000 × voxels in dimensions to mainly contain the trabecular bone and reduce the amounts of debris.
A slight shift in sample surface colour was observed after the CTscans, indicating that the X-ray radiation affected the samples.Exactly how much of the samples that were affected or how it may alter their mechanical properties is still an open scientific question (Peña Fernandez et al., 2018a;Barth et al., 2011).However, the discussion is outside the scope of this work and we wish to only note down that there was an observable effect on the surface colour.

Theory
For numerical computation, phase-field models for brittle fracture implemented with a linear finite element discretization are utilized.The implementation is described for isotropic materials; however, it is possible to extend the theory to apply to materials with elastic anisotropy, cf.van Dijk et al. (2020).As derivations for the dynamic and quasi-static theories are similar, we first present the dynamic theory and then summarize key equations for the quasi-static theory.
Let a sharp crack Γ be contained in an arbitrarily shaped isotropic elastic body Ω.The body is subject to prescribed boundary displacement U on part of the boundary ∂Ω U ⊂∂Ω and traction T on ∂Ω T ⊂∂Ω (Fig. 2).The sharp crack is approximated by a crack density functional, , often referred to Ambrosio-Tortorelli 1 (Ambrosio and Tortorelli, 1990;Bourdin et al., 2000) where ∇ is the spatial differential operator, l is a regularization parameter characterizing the width of the damaged region and d is the damage field defined over the whole body Ω, d ∈ [0, 1] where d = corresponds to undamaged material and d = 1 to completely damaged material (Fig. 2).
A small deformation elasticity theory is assumed.Let the displacement field be denoted u.The Lagrangian of the fracture problem in absence of any volume forces is given by where, ψ k = 1 2 ρ u⋅ u is the kinetic energy density, ρ the mass density, u, (2) The split of the strain energy density into two parts aims to ensure that the material does not degrade under compression and stiffness is kept in the case of crack closure.However, the choice of split is not obvious, and many different splits have been suggested (cf.Hesammokri et al., 2022).Here, a hydrostatic-deviatoric split (Amor et al., 2009) is applied in which the hydrostatic and deviatoric deformation modes are the independent modes of deformation which can be split into positive and negative contributions according to where ε = 1 2 [∇u +(∇u) T ] is the linearized strain tensor, e the linearized strain deviator, the parameter α = 1 if the volume expands, i.e., tr(ε) > 0, and α = 0 otherwise, K and μ are the bulk and shear modulus, respectively.One may notice that this split assumes that all strain energy associated with shear deformation affects damage, also under compressive load.Thus, it is assumed that crack growth in bone is driven by shear and tensile bulk deformation.Further, to impose irreversibility of the crack evolution in the computer simulations, a history field H is used in place of ψ + e , such that H is the maximum positive strain energy density experienced (Miehe et al., 2010), i.e., ψ + e (u, τ). (4) The history field (4) ensures that it is the highest strain energy density experienced in the material during the loading history that determines the present damage state, i.e., the material cannot heal.Now, using the above relations, the energy functional is written as, cf.Borden et al. (2012), dΩ . (5) The first variation of ( 5) should be zero according to the principle of least action.The Euler-Lagrange equations of ( 5) with respect to the displacement and damage fields are, where the Cauchy stress tensor is given by σ The boundary and initial conditions are given by where x denote the spatial coordinates and n is the outward surface normal at material boundaries.
The dynamic displacement problem is solved using a staggered time integration scheme.At each time step the crack phase problem (i.e., d) is first solved independently, based on the displacements of the previous step.Then, the displacement is solved using a Newmark (1959) algorithm.The displacement of the new time step n + 1 is predicted by u n+1 , using the values of the previous time step n according to The acceleration at the next time step is obtained as where , M is the mass matrix, K n+1 is the stiffness matrix and F n+1 is the nodal force vector (at time step n + 1).
The updated acceleration ün+1 is then input into (8) to obtain the updated velocity and displacement fields.A fully explicit Newmark algorithm is used, in which the constants β 1 = 1/2 and β 2 = 0.A structural damping term could be included in ( 9), however, in this study, no structural damping was used which means that all energy dissipated from the dynamic system is due to fracture.
In the quasi-static case, the rate-dependency is excluded from the total energy equation ( 5).The governing equations are derived from variational minimization of ( 5) with respect to d, i.e.
The boundary conditions are given by ⎧ For quasi-static loading, i.e., when stress equilibrium prevails, it can be shown that the model for a homogeneous body gives a uniaxial (Pham et al., 2011), and a uniaxial tensile fracture strain ε 0 = [3G c /(8El)] 1/2 .For convenience, a fracture energy The quasi-static simulation is solved with a staggered Newton-Raphson scheme, such as described in, among others, Espadas-Escalante et al. (2019).At each load step n, the displacement field u and damage field d are solved and updated sequentially in an iteration process until the incremental change of the damage field Δd is sufficiently small (max{|Δd|}〈10 − 3 ) and stress equilibrium prevails.A flowchart of the staggered scheme adopted is shown in Fig. 3.
In both the dynamic and quasi-static problems, discretization is done by interpolating the displacement and phase-field variables with standard finite element shape functions.Typical discretization examples can be found in Nguyen et al. (2015), Carlsson and Isaksson (2018).All the equations were solved using well-established finite element algorithms implemented in a Matlab (2020) code.Cubic 8-node isoparametric tri-linear elements were used in all computations.

Model
CT image stacks were imported for the 3D model.In order to enable reasonable speed of the computations, the images were run through an iterative procedure to coarsen the finite element mesh and remove isolated debris.A convergence study was performed on a limited volume extracted from one of the CT scans.The volume was large enough to contain enough topological details for the convergence test, yet small enough to allow manageable computational times also on fine meshes.
With the exception of the very coarsest mesh, the global stress-strain response is reasonably consistent for all meshes considered.Fracture patterns are similar for meshes of element edge-length h smaller than ∼ 32.5 μm.In general, an element edge-length h < l/2 is necessary to obtain a sufficient resolution of the regularized crack surface (Ambati et al., 2015;Borden et al., 2012;Carlsson and Isaksson, 2018, 2020a, 2020b;Molnár and Gravouil, 2017).This is verified by keeping l fixed and varying h, which confirms that l = 2h is converged.The convergence study is described in more detail in Appendix A. With these findings in mind, an element edge-length of h = 32.5 μm was chosen, giving a final model of about half a million elements, Fig. 4. The characteristic length parameter was chosen as l = 2h.One may observe in Figs.4-5 that a typical trabecula has about 0.2 mm in thickness, i.e., about 6h or 3l.Thus, the width of a trabecula is about 3 times the length parameter l which means that crack propagation within (as opposed to across) a trabecula cannot be captured well by the model.Crack propagation across trabeculae can be captured, although not in detail due both to the lack of detail at sub-element length scales and the width of the crack.This is considered acceptable, since slender structures tend to fracture transversely.
For simplicity, the base material is assumed to be isotropic, linear elastic and subject to small strains.In reality, trabecular bone has a small degree of both anisotropy and nonlinearity.Moreover, there are slightly darker areas along the edges of the trabeculae which can be seen in e.g.Fig. 5, indicating lower density and possibly different material properties.These variances are not included in the models, but might be of interest for future studies.The material parameters were chosen as: Young's modulus E = 1 GPa, Poisson's ratio ν = 0.2, a critical energy release rate G c = 20 J/m 2 and density ρ = 1600 kg/m 3 .These parameters give a uniaxial strength σ 0 of 11 MPa and a uniaxial fracture strain ε 0 of 1% (which is in agreement with experimental studies, cf.Peña-Fernández et al., 2018b).As pointed out in the Introduction, according to the literature there are large uncertainties about the elastic material properties of bone (Wu et al., 2018(Wu et al., , 2022)).The critical energy release rate, or fracture toughness, would be even more difficult to estimate accurately.A quantitative comparison of loads and displacement is therefore outside the scope of this study, hence the somewhat arbitrary choice of material parameters.However, a parameter study, presented in Appendix A, shows that the scaled results are more or less unaffected by the exact choice of parameters.
In the quasi-static simulation, incremental displacements Δu = H/ 10 4 were given in the vertical direction (z) to the body's top and bottom boundaries at each load step, Fig. 6.All other boundaries remained stress free.In the dynamic simulations, prescribed velocities v z were given in the vertical direction (z) to the body's top and bottom boundaries at each time step while all other boundaries remained stress free.The velocity v z was linearly increased from v z = 0 to v z = v 0 during the first 100 time-steps and then kept at a constant velocity v 0 for the remainder of the simulation (Fig. 6).The time step used was Δt = 0.01 μs.Dynamic simulations with four different values of v 0 were performed: v c , v c /2, v c /4 and v c /8, where the maximum load rate v c = 4.44 m/s was chosen as the resulting velocity from a free fall from a 1 m height.This velocity provides an upper bound on the impact velocity in an actual fall, as in a real fall soft tissues surrounding the bone (such as fat tissue, muscles etc.) would absorb energy and decrease the velocity; to what extent is outside the scope of this study and will also be situationspecific.In both the static and the dynamic models, the simulations were stopped once the global reaction forces became lower than 10% of the peak forces obtained in respective simulation.

Results
The results are presented in three parts: experimental results, quasistatic numerical results contrasted to the experiment, and then the dynamic rate-dependent numerical results.Two samples were used in the study, giving very similar outcomes.For consistency, only one sample is used to present the results in the following sections.

Experiment
Resulting fractures observed in the experiment are discussed based on the 3D-image data from the CT.Analysis will be in a purely qualitative manner.Inspection of CT images in Fig. 7 indicate that on the global scale, cracks nucleate and start growing at the lower right corner of the scanned region, and stepwise grow towards the notch.It is possible that the notch length (2.3 mm) was not sufficiently large to On the scale of a single trabecula, one observes that some microcracks eventually grow and coalesce, causing complete fracture of a trabecula, Fig. 8.One also observes branching and arrest of cracks, as well as initiation of new cracks that eventually merge with other cracks.Also, cavities in the trabecula seems to localize microcracks, Fig. 8.These observed fracture mechanisms confirm earlier studies (Thurner et al., 2006).

Comparison between the experiment and the quasi-static simulation
The experiment was simulated with the quasi-static phase-field numerical model discussed above.Fig. 9 shows a comparison between experiment and simulation for a representative cross-section.The simulation captures many of the final crack paths.There are also a few regions with disagreement, e.g., close to the top boundary in Fig. 9a where one from the simulation would expect a fracture, while it is difficult to visually detect any crack in the same region in the experiment.However, the majority of fractures obtained in the numerical simulation are also observed in the experiment, cf.Fig. 9b and c.
Interestingly, while the final fracture paths are similar in the experiment and the simulation, the crack growth process itself differs slightly, as visualized in the fracture sequences in Fig. 10.In the experiment, a crack opened close to the right boundary of the sample.As load-carrying trabecula fractured, the load was re-distributed and new cracks were nucleated in a step-wise process toward the notch.The simulated fracture path is similar to the one observed in the experiment, however, the first simulated cracks appear closer to the notch and not at the remote boundary.
Fig. 11 shows close-ups of a cracked region (dashed box in Fig. 10) to visualize the coarsening effect on the numerical model, as compared to the original high-resolution image (Fig. 8a).It is noted that in the FE model, the mesh is too coarse, and the length parameter l too large, in order for the numerical model to predict the detailed crack path observed in the experiment.It is not possible for the model to capture microcracks at length-scales which are smaller than the element size (element edge-length h).Nevertheless, high-stressed region may be revealed, as indicated in Fig. 11.

Comparison between quasi-static and dynamic simulations
For the dynamic computations, increasing load rates were applied to the boundary.The aim was to investigate how the material behaves at  Fig. 6.Left: Cross-section of the mesh with schematic of boundary conditions: prescribed velocity in the dynamic simulation and stepwise displacement in the quasistatic one.The blunt notch is at the left vertical edge.Right: Profile of the prescribed velocity v z in the dynamic case.The velocity v z is linearly increasing until t = t 0 (t 0 = 1 μs) where v z = v 0 .At t ≥ t 0 , v z is held constant.higher load rates and to what extent the behaviour is changed when loaded at realistic loading rates, i.e., whether it is of interest to study fracture at loading rates close to the rates expected in actual fracturing incidents.
Fig. 12 shows predicted normalized global stress-strain curves for different loading velocities.The global stress is given by the reaction forces on the top and bottom boundaries divided by the structure's crosssection area, while the global strain is given by the displacements of the top and bottom boundaries related to the initial distance between them.One observes that the dynamic simulations converge toward the quasistatic model as the load-rate decreases, which lend confidence to the simulations.During the initial acceleration of the boundaries, at the impact, the dynamic loadings show a steep slope, i.e., relatively higher global stiffness, because of the impulse, as compared to the quasi-static loading case.The higher the impulse, the stiffer the response, because of the inertia.Once the load rate v 0 is reached, there is a relaxation window, where the mass of the porous material adjusts to the set load rate, i. e., inertia effects.During this window the global stress slightly decreases but stays relatively stable around a value just below the global stress at which v 0 was reached.The global stress increases once again when stress equilibrium starts to establish throughout in the body.This can be seen in Fig. 12 when the dynamic stress-strain curves join the quasi-static Fig. 7. Developing cracks at increasing displacement (time t 1 < t 2 < t 3 < t 4 ).Red markings indicate newly formed cracks and yellow markings indicate cracks formed at earlier load steps.Several cracks are observed that do not continue to propagate.At the last image (t 4 ), cracks that did not propagate becomes wholly or nearly invisible in the CT scan when the crack surfaces close, as the load is relieved from that region.Each cross-section has dimension 6 × 3 mm and the horizontal blunt notch can be seen at midheight at the left vertical edges.curve.Inertia effects also result in higher global peak-stress in the dynamical cases.Higher impulses result in higher peak-stresses, and at larger boundary displacements, because of inertia effects.
In the numerical models, energy can only dissipate as surface energy by fracture, i. e., at higher v 0 ).The positions of maximum elastic energy stored in the body corresponds to the global peak stress (Fig. 12).At global peak stress, the crack density in the porous body start to increase whereupon elastic energy is dissipated through creation of crack surfaces.As may be seen in both Figs. 12 and 13, lower load rates come closer to, but do not completely recover the quasi-static results.Crack growth in the quasistatic simulation is unstable, as can be seen in the vertical jumps in Figs.12-13.These vertical drops indicate incremental load steps in which significant damage growth take place in load-carrying regions, leading to substantial global stiffness loss.This is a typical feature in phase-field fracture models solved with a fully staggered scheme together with the linear Ambrosio-Tortorelli approximation of the Munford-Shah functional (Pham et al., 2011(Pham et al., , 2017;;Espadas-Escalante et al., 2019;Ambati et al., 2015).Here we interpret this phenomenon as unstable crack growth, even though the term often is associated with dynamic phenomena.In these "unstable crack growth load steps", the load in the trabecular network is redistributed while the damage field expands and stress equilibrium prevails.From one stable solution of the crack field in the minimization problem, to another stable solution quite far away from the one in the previous load step.Nevertheless, as Figs.12-13 reveals, the dynamic model asymptotically approaches the quasi-static solution as the impulse vanish, i.e., in the limit when v 0 →0.The manner in which the cracks initiate and grow, as well as their locations are somewhat different in the quasi-static and dynamic loading simulations (Fig. 14).The cracks developed in the quasi-static simulation (and most of the cracks observed in the experiment) are contained within a narrow region at approximately the same height as the notch, in a plane-like appearance despite the large porous structure.On the other hand, cracks developed at higher load rates are spread throughout the whole volume, closer to the top and bottom.At larger impulses (i.e., at higher v 0 ), stress equilibrium is not established, whereupon stresses and strains are higher closer to the top and bottom boundaries, causing   Computed global stress-strain relations at different load rates contrasted with the quasi-static simulation.The global stress is given by the reaction forces on the top and bottom boundaries divided by the structure's crosssection area, while the global strain is given by the displacements of the top and bottom boundaries related to the initial distance between them.The global stress σ is normalized with the homogeneous tensile strength σ 0 , while the global strain ε is normalized with the homogeneous fracture strain ε 0 .cracks nucleating at locations deviating from the expected ones at stress equilibrium (quasi-static model, experiment).

Discussion
The experiments showed a complex stepwise crack evolution with multiple crack fronts and crack arrests.As there were only two samples, generalizations may be limited, but similar fracture behaviour have previously been reported by e.g.Thurner et al. (2006) and Voide et al. (2009).The simulated quasi-static fractures were similar to the experimental ones, with few discrepancies.The largest normal strains were less than 4% in uncracked regions indicating that a small deformation theory is applicable.Some experimental factors that may have affected the results are important to point out.Firstly, damage from extracting, machining and mounting the samples might affect the local material properties and therefore also the load distributions and crack nucleation sites.As mentioned in Section 3.2, the location of the first trabecular fracture differed between the simulation and experiment.It is possible that the first crack observed in the experiment was always present but was closed during the initial scan (unloaded sample).The crack could have nucleated during the sample preparation or mounting in the load stage.However, this is only speculations; we don't know.Secondly, the CT scans did not include the full height, or even width, of the samples.While the final fractures occurred at the height of the notch (approximately in the samples mid-height), there is no knowing how the trabeculae were damaged outside the scanned regions.This very much relates to the question of boundary conditions.The few differences in fracture locations between the simulation and the experiment are believed to be related to boundary effects.The numerical model is loaded by continuously increasing the displacements on the top and bottom boundaries, but since the CT scan on which the model geometry is based only covers part of the sample, this approximation may deviate from the experiment.To improve agreement between experiment and simulations it would be desirable to adjust the experimental setup in order to give more reproducible boundary conditions Thirdly, the generality of the results may be limited due to the effect of radiation (Barth et al., 2010).When exposed to radiation, it has been shown that bone tissue tends to become more brittle in nature.In what way and how much this might have affected the fracture appearance and behaviour in the experiment is difficult to deduce (Peña Fernández et al., 2018a, 2018b).Thus, the application of the phase-field fracture model in this study should be considered explorative and not as a representation of the tissue mechanics.Despite the experimental limitations, the high-resolution images of the CT scans reveal interesting fracture behaviour with crack arrests, deflection and microcrack nucleation.
The quasi-static phase-field fracture method is well established and have been previously applied in studies on various heterogeneous materials (Nguyen et al., 2015;Espadas-Escalante et al., 2019;Carlsson and Isaksson, 2018) with promising results.Here too, predicted cracks are mostly found in similar positions as the experimental observed cracks.It is, however, important to note that the phase-field model in this investigation uses a relatively coarse mesh and thus cannot resolve fracture on the length scale of individual trabeculae, cf.Fig. 8 and Figs.10-11.The cavity that attracts, or localizes, a crack in Fig. 8, for instance, is not captured in the model, Fig. 11.The numerical method is not able to capture the microcracking seen in the experiment, nor is it able to capture small cavities of sub-element size in the trabecula.Coarsening the mesh was necessary for a reasonable computational cost, but as is observed in Fig. 11, morphological information may be lost.Similarly, a too coarse mesh naturally influences the width of the diffusive and regularized cracks since a crack cannot be thinner than two element edge-lengths.To improve the local agreement, one likely has to increase the resolution in the finite element mesh of the phase-field simulation, leading to a significantly increased computational cost.The cost addition may however be kept to a minimum by improvements such as local mesh refinement or multiscale strategies.An interesting possibility would be to use different mesh resolution depending on the local value of damage.It is worth mentioning that the computation time for the quasi-static case was approximately ten days (using about 18 thousand stress equilibrium iterations and roughly 2 million dofs).
A challenging problem is how to represent the material behaviour of the bone tissue (Helgason et al., 2008;Morgan et al., 2018;Du et al., 2019;Bauer et al., 2014).There are multiple studies on the building blocks of trabecular bone tissue, suggesting nonlinear and direction-dependent material properties, including composite mechanisms for fracture resistance.Bone tissue may also experience viscoelastic and viscoplastic behaviour.In this study, a simple, homogeneous brittle linear elastic isotropic material is assumed.The parameter sensitivity study described in Appendix A shows that the results are consistent over a wide range of values of material parameters; only sublime differences can be noticed.This suggests thatqualitativelythe results obtained using this simplified model should be valid also for a material with a low degree of anisotropy, non-linearity or rate-dependence.It is, however, difficult to say with certainty to what extent the constitutive assumptions have affected the final fracture paths and in what way.In summary, there are several aspects to consider to improve the reproduction of experimental results, but the basic phase-field fracture model used here shows promising results and potential.
The employed dynamic phase-field method seems to capture the dynamic events of fracture.As the load rate was decreased, fracture paths and global load-displacement curves converged toward the quasistatic case.Comparing higher load rates to lower, i.e., larger impulses to smaller impulses, one generally observes a rapid passage of events with higher amount of dissipated energy, i.e., larger amounts of damage/ crack growth, and that the damage/crack growth is spread over a larger region in the sample, a phenomenon resulting from inertia effects (the ability of the trabecular bone material to resist changes in motion varies with its mass, or density).At larger impulses, when stress waves propagate through the material, the local strain energy in a trabecular crosssection may momentarily become sufficiently high to separate material, i.e., microcrack nucleation.Such nucleated cracks may generate new stress waves, which are reflected in the trabecular bone network, and interference and resonance phenomena may take place.Such stresswave driven microcracks may continue to nucleate until steady-state is reached, i.e., when stress equilibrium prevails.A final global fracture may occur if microcracks coalesce and localize to form an unstably running crack, which in turn emits stress waves that can be reflected or scattered back to the crack-tip region.Thus, bone subject to larger impulses naturally get more energy induced and fractures nucleates and propagates at more places.
The results confirm previous studies on dynamic fracture (cf.Carlsson and Isaksson, 2018, 2020a) which also show more dispersed damage for higher loading rates.The chosen phase-field fracture theory is a convenient way to analyse bone subject to rapid loads.Higher resolution of the structural mesh will enhance morphological and fracture details, however, at a substantially higher computational cost.It is a delicate balance act to determine the mesh resolution, at the end of the day it is linked to the mechanical problem to be analysed.
Finally, here we have analysed a bone volume subject to global tension.In reality, mechanical loading of bone is intricate and highly complicated.Real trabecular bone may be subject to mixed-mode loading, i.e., complex combinations of e.g., tension, compression, torsion, or shearing.Additionally, long-time effects may be present, such as creep, healing or fatigue.These subjects are topics for future studies.

Conclusions
The quasi-static experiment revealed complex stepwise crack evolution with multiple crack fronts, and crack arrests, as the global displacement load was incrementally increased.The quasi-static phasefield brittle fracture model captures the fractures obtained in the experiment reasonably well.The dynamic model converges towards the quasi-static model when loaded at low load rates.At higher load rates inertia effects contribute to an increased initial stiffness, higher global peak stresses and larger amounts of cracks spread over a larger volume.Since the fracture process clearly is different at large impulses compared to small impulses, dynamic fracture models are necessary when simulating rapid bone fracture.
Due to a relatively coarse mesh, the model could not replicate crack arrests, deflections or branching induced by tissue microstructure (e.g., cavities) present at a smaller length scale than the mesh resolution (i.e., element edge-length).Still, the dynamic phase-field fracture model shows promising results and potential.Higher resolution of the structural mesh will enhance morphological and fracture details, however, at a substantially higher computational cost.
Figure A1 displays the global stress-strain curves and the final fracture patterns using a phase-field characteristic length l = 2h for the evaluated mesh sizes.The global stress is given by the reaction forces on the top and bottom boundaries divided by the structure's cross-section area, while the global strain is given by the displacements of the top and bottom boundaries related to the initial distance between them.Figure A2 shows the global stress-strain curves and the final fracture patterns using a constant characteristic length l = 20h 0 = 65 μm for some mesh sizes.The simulations shown in Figs.A1 and A2 were run using Young's modulus E = 1 GPa, Poisson's ratio ν = 0.2, and a critical energy release rate G c = 20 J/m 2 .Based on the observations in Figs A1-A2, a mesh resolution h = 32.5 μm and a characteristic length l = 65 μm seem to provide consistent reults and is used in the study.

Fig. A1.
Global stress-strain relations and final fracture paths when using different values of h/h 0 and l = 2h.The global stress σ is normalized with the tensile strength σ 0 , while the global strain ε is normalized with the fracture strain ε 0 .For convenience, both σ 0 and ε 0 are constant and determined using l = 2h 0 (corrsponding to the case h/h 0 = 10).

Fig. A2.
Global stress-strain relations and final fracture paths when using a constant l = 20h 0 while using different values of h/h 0 .The global stress σ is normalized with the tensile strength σ 0 , while the global strain ε is normalized with the fracture strain ε 0 .
To evaluate the models' sensitivity to chosen material parameters, several simulations were made using different parameters.The same boundary conditions as above were applied, using a mesh size h/h 0 = 10.The simulations were run using Young's modulus E = 1 or 10 GPa, a Poisson's ratio ν = 0.2 or 0.4, and a critical energy release rate G c = 20 or 200 J/m 2 .Hence, 8 different simulations were made using unique set-ups of the material parameters E, ν and G c . Figure A3 shows some global stress-strain relations and final fracture paths, revealing that the simulations give consistent reults and the precise values of the selected material parameters E, ν and G c are not critical.A Young's modulus E = 1 GPa, a Poisson's ratio ν = 0.2, and a critical energy release rate G c = 20 J/m 2 is used in the study.

Fig. 1 .
Fig. 1. 3D image of trabecular bone sample by reconstructed CT data.The scanned region has dimension 6 × 3 × 3 mm and contains a notch with length of about 1.4 mm.As indicated in the schematics, the scan does not include the full sample.

Fig. 2 .
Fig. 2. A sharp crack Γ contained in a porous domain Ω (a).Diffuse approximation of the sharp crack, γ, with prescribed displacement U and prescribed traction force T acting on the parts of the domain surface ∂Ω (b).Profile of the diffusive phase-field crack as an approximation to a sharp crack (c).Here, ξ is a spatial coordinate axis directed orthogonal to the crack plane.

Fig. 3 .
Fig. 3. Flowchart of the staggered scheme used for the quasi-static simulation.

Fig. 4 .
Fig. 4. Discretized three-dimensional body of a notched trabecular bone sample.The cubic element edge-length is h.

Fig. 5 .
Fig. 5. Cross-section of trabecular strut.A slightly darker area can be seen along the edges of the trabecula, one of them marked with an arrow.Darker colour indicates a lower density material, that seems to make up an outer layer of the trabecular bone.

Fig. 8 .
Fig. 8. a) A crack is attracted to a cavity.The crack travels through some osteocytes (indicated with arrows).b) Crack arrest and appearance of another crack.c) The two cracks merge and completely separates the top and bottom parts of the trabecula.Each cross-section has dimension 0.6 × 0.3 mm.The region is indicated with a dashed box in the upper left image in Fig. 7.

Fig. 9 .
Fig. 9. Comparison of fracture locations in three cross-sections (a-c).Top row is the experiment with cracks marked by red circles.Bottom row is numerical results with more damage the brighter the colour.Black corresponds to no damage (d = 0) and bright/white is fully damaged (d = 1).Each cross-section has dimension 6× 3 mm and the horizontal blunt notch is located at the left vertical edges.
e., Π s = W − Π k − Π e , where W is the work done by external forces, Π k = ∫ Ω ψ k dΩ is the kinetic energy and Π e = ∫ Ω ψ e dΩ is the elastic (or potential) energy.Elastic energy and dissipated energy are shown in Fig. 13.Higher values of energy are seen at larger impulses (i.

Fig. 10 .
Fig. 10.Fracture propagation in a cross-section with increasing displacement from left to right.Experimental images on top row show fracture initiation close to the boundary opposite to the notch.The dashed box in the lower right image is the area shown in Fig. 11.Each cross-section has dimension 6 × 3 mm and the horizontal notch can be seen at mid-height at the left vertical edges.

Fig. 11 .
Fig. 11.Close-up of mesh and CT image in Fig. 8a to reveal the effect of coarsening on the resolution.The region is marked with a dashed box in the lower right image in Fig. 10.

Fig
Fig. 12.Computed global stress-strain relations at different load rates contrasted with the quasi-static simulation.The global stress is given by the reaction forces on the top and bottom boundaries divided by the structure's crosssection area, while the global strain is given by the displacements of the top and bottom boundaries related to the initial distance between them.The global

Fig. A3 .
Fig. A3.Global stress-strain relations and final fracture paths when using Young's modulus E = 1 or 10 GPa, a Poisson's ratio ν = 0.2 or 0.4, and a critical energy release rate G c = 20 or 200 J/m 2 in the simulations.The global stress σ is normalized with the tensile strength σ 0 , while the global strain ε is normalized with the fracture strain ε 0 .
later occurring ü, denotes first and second derivatives of u with respect to time t, ψ e is the elastic strain energy density and ψ s =